- Research article
- Open Access
Transcriptomics analysis of Toxoplasma gondii-infected mouse macrophages reveals coding and noncoding signatures in the presence and absence of MyD88
BMC Genomics volume 22, Article number: 130 (2021)
Toxoplasma gondii is a globally distributed protozoan parasite that establishes life-long asymptomatic infection in humans, often emerging as a life-threatening opportunistic pathogen during immunodeficiency. As an intracellular microbe, Toxoplasma establishes an intimate relationship with its host cell from the outset of infection. Macrophages are targets of infection and they are important in early innate immunity and possibly parasite dissemination throughout the host. Here, we employ an RNA-sequencing approach to identify host and parasite transcriptional responses during infection of mouse bone marrow-derived macrophages (BMDM). We incorporated into our analysis infection with the high virulence Type I RH strain and the low virulence Type II strain PTG. Because the well-known TLR-MyD88 signaling axis is likely of less importance in humans, we examined transcriptional responses in both MyD88+/+ and MyD88−/− BMDM. Long noncoding (lnc) RNA molecules are emerging as key regulators in infection and immunity, and were, therefore, included in our analysis.
We found significantly more host genes were differentially expressed in response to the highly virulent RH strain rather than with the less virulent PTG strain (335 versus 74 protein coding genes for RH and PTG, respectively). Enriched in these protein coding genes were subsets associated with the immune response as well as cell adhesion and migration. We identified 249 and 83 non-coding RNAs as differentially expressed during infection with RH and PTG strains, respectively. Although the majority of these are of unknown function, one conserved lncRNA termed mir17hg encodes the mir17 microRNA gene cluster that has been implicated in down-regulating host cell apoptosis during T. gondii infection. Only a minimal number of transcripts were differentially expressed between MyD88 knockout and wild type cells. However, several immune genes were among the differences. While transcripts for parasite secretory proteins were amongst the most highly expressed T. gondii genes during infection, no differentially expressed parasite genes were identified when comparing infection in MyD88 knockout and wild type host BMDM.
The large dataset presented here lays the groundwork for continued studies on both the MyD88-independent immune response and the function of lncRNAs during Toxoplasma gondii infection.
The intracellular apicomplexan Toxoplasma gondii is a globally distributed parasitic microorganism infecting both humans and animals. In humans alone, Toxoplasma is conservatively estimated to be present in over a billion individuals . After ingestion of tissue cysts or oocysts, an acute phase commences characterized by parasite dissemination throughout the host as rapidly dividing tachyzoites. This is followed by establishment of latent infection, in which tachyzoites differentiate into slowly replicating bradyzoites that form cysts in tissues of the central nervous system and skeletal muscle [2, 3]. Latent, or chronic, infection is asymptomatic in most cases, but the parasite may reactivate in immunocompromised populations leading to life-threatening disease . Primary infection during pregnancy can lead to major birth defects and sequelae of infection later in life .
Toxoplasma is well known for its ability to stimulate strong Th1 immunity that has as its origin early production of IL-12 by dendritic cells [6, 7]. The IFN-γ produced during infection confers resistance to the parasite, and indeed this cytokine is central in the ability to survive acute Toxoplasma infection . While protective, IFN-γ production can result in host pathology if not appropriately regulated by counter-inflammatory cytokines such as IL-10 [9, 10]. A major function of IFN-γ is to elicit inflammatory macrophages that are major anti-microbial effectors during in vivo infection [11,12,13]. Paradoxically, macrophages along with dendritic cells also serve as early cells targeted for infection, and it has been suggested that they act as Trojan horses to enable establishment of T. gondii in the host [14,15,16,17,18]. For these reasons, macrophages are an especially important cell type to study both the host immune response and T. gondii behavior during intracellular infection.
Substantial work in mouse models has revealed an important role for Toll-like receptor (TLR) and the adaptor molecule MyD88 in innate immune recognition of T. gondii [19, 20]. The invasion-associated parasite protein profilin functions as a ligand for TLR11 and TLR12, initiating MyD88-dependent immunity [21,22,23,24]. Given the central role of the MyD88 protein in the early innate immune response in mice to T. gondii infection, it is important to understand how deletion of MyD88 impacts transcription of downstream immune genes in infected cells. In humans, the basis of immune recognition is less clear because TLR11 is present as a pseudogene and TLR12 is absent from the genome . Furthermore, a study of a pediatric population with an autosomal recessive MyD88 deficiency revealed that these individuals retain resistance to all but a minimal number of pyogenic bacterial infections [26, 27]. Thus, determining MyD88-independent responses to infection with Toxoplasma and other microbial pathogens is an important avenue of investigation in both humans and mice.
We therefore employed RNA sequencing (RNA-seq) to determine the transcriptome of MyD88+/+ and MyD88−/− bone marrow-derived macrophages (BMDM) following infection with T. gondii. In addition to yielding information on protein coding responses, RNA-seq provides insight into responses of long noncoding RNA (lncRNA), defined as transcripts greater than 200 nucleotides with no protein coding potential. lncRNAs are widely involved in gene regulation, and their study is an emerging area of interest in infection and immunology [28,29,30,31].
Our approach involved infection with high virulence (Type I strain RH) and low virulence (Type II strain PTG) isolates of Toxoplasma. Amongst Type II strains, some differences in the intensity of cytokine responses have been noted with different isolates but we employed a strain that has been extensively used in previous studies . In mice, Type I strains induce a hyperinflammatory cytokine response rapidly culminating in host death. The immune response is more restrained during Type II infection, enabling host survival and parasite establishment of latent infection [33, 34]. In vitro studies have revealed that infection with these strains activates partially nonoverlapping host signaling pathways leading to distinct responses. For example, infection with Type I parasites triggers strong and sustained activation of STAT3 and STAT6 resulting in the generation of macrophages with an M2 phenotype [35, 36]. Type II infection triggers NFκB activation and robust IL-12 production . The present study provides important information on global transcriptional changes in macrophages infected with these two Toxoplasma strains in the presence and absence of MyD88.
In addition to examining the transcriptional changes in macrophages, use of RNA-seq technology enabled us to simultaneously harvest data on the transcriptomes of high and low virulence Toxoplasma during initial stages of intracellular infection. This allowed us to compare gene expression differences between T. gondii strains, as well as examine differences in parasite gene expression when infecting MyD88+/+ versus MyD88−/− macrophages. Together, this dataset provides a host and parasite genomic framework for understanding the interactome that emerges during intracellular infection with Toxoplasma.
Dual RNA sequencing of Toxoplasma-infected macrophages identifies host and parasite transcripts
We infected wild type and MyD88 knockout (KO) bone marrow-derived macrophages with both Type I (RH) and Type II (PTG) Toxoplasma tachyzoites, then collected samples 6 h later for high throughput RNA-seq (Fig. 1). We selected 6 h for our analysis because this time point occurs prior to the first parasite mitotic division, controlling for differences in replication rate between Type I and Type II strains. This time also enabled us to determine the earliest macrophage responses to infection. We used a multiplicity of infection (MOI) of 4:1 or 5:1. The percent infection, measured via fluorescence microscopy at 6 h over multiple biological replicates, ranged from 70.4–95.7%. Sequencing was performed on 4 biological replicates for uninfected samples and 3 biological replicates for infected samples. We mapped mouse reads to GENCODE version M21, a database containing sequences for 58,899 protein-coding transcripts and approximately 30,462 lncRNA transcripts. We mapped Toxoplasma reads to ME49 strain sequences in the ToxoDB database. As expected, mouse sequences comprised the vast majority of reads in infected samples (Fig. 2). The percentage of reads mapping to the T. gondii genome ranged from 2.49 to 18.20% between samples. The variability in percentage of mapped reads between replicates correlated only weakly, at best, with percent infection measured via fluorescence microscopy (RH R2 = 0.18; PTG R2 = 0.01) but did correlate more strongly with the number of parasites per infected cell (RH R2 = 0.15; PTG R2 = 0.68). Factors in addition to the number of parasites present likely contribute to the variability in the number of parasite transcripts produced. The background level of T. gondii reads in uninfected samples was determined to not affect the parasite gene expression results presented herein and likely represent mis-mapping of mouse genes, since housekeeping genes were among the top Toxoplasma genes in uninfected samples. Principal component analysis (PCA) plots of mapped mouse reads demonstrate that the treatment (in this case parasite strain) accounted for data variability, but also that the biological replicate contributed substantially to the variability observed (Additional files 1, 2 and 3). PCA plots of mapped Toxoplasma reads demonstrate that the parasite strain largely accounted for data variability (Additional files 4, 5).
Type I strain parasites trigger stronger protein-coding gene expression effects compared to type II parasites
We defined differentially expressed (DE) transcripts using a p-value of equal or less than 0.05 and a fold change of 2 or greater. Among the wild type mouse samples, DE transcripts were primarily protein-coding (51%) and non-coding (43%), with pseudogenes and TEC (To be Experimentally Confirmed) comprising 4 and 2% of the hits respectively (Fig. 3a). A complete listing of DE transcripts for the three wild type comparisons (RH vs uninfected, PTG vs uninfected, and RH vs PTG) is shown in Additional file 6.
Among the protein-coding sequences, substantially more DE transcripts were differentially expressed with the Type I RH infection versus with the less-virulent Type II PTG infection (335 and 74, respectively). This indicates that RH has a stronger impact on the host macrophages relative to PTG. 57 transcripts were differentially expressed between RH and PTG, including previously known immune genes ccl24, csf1, socs2 and ccl17 (Fig. 3b and Additional file 6). Venn diagrams reveal that 46 DE transcripts (62%) are shared between RH and PTG infection (Fig. 3c).
Heat maps demonstrate that many genes related to the immune response, cell cycle, DNA replication, DNA recombination, DNA repair, growth, cell adhesion, and cell migration were differentially expressed during T. gondii infection (Fig. 3d). Numerous immune genes were of higher or lower abundance during RH infection, confirming that activation and suppression of immunity during Toxoplasma infection extends to the cellular level (Fig. 3d). In confirmation of previous studies [35, 36, 38], immune-related genes Arg1 and ccl17 were more abundant in RH versus uninfected cells. Many cell adhesion and migration genes were more abundant in both the RH versus uninfected and RH versus PTG comparisons. Many genes related to the cell cycle were differentially expressed in both RH versus uninfected and PTG versus uninfected. The DE genes for cell cycle include several genes relating to microtubule organizing center and DNA replication, recombination, and repair. This is of interest because Toxoplasma is thought to co-opt microtubules for its own survival [39, 40]. Many cell growth genes were of higher abundance, particularly for RH versus uninfected. Gene ontology analysis and KEGG pathway analysis results support the data shown in the heat maps (Additional file 7). In addition to the functional categories displayed in the heat maps, gene ontology analysis indicates that metabolic processes are also strongly differentially expressed in the RH strain (Additional file 7). Additional heat maps with each replicate displayed individually reveal that sample wtRH_4 had stronger effects than other replicates but the same overall trends (Additional Files 8 and 9).
Many noncoding transcripts are differentially expressed during infection of wild type BMDM
Among the mouse reads mapped to GENCODE/Ensembl, we analyzed the noncoding transcripts separately from the protein-coding transcripts. The majority (63%) of the differentially expressed noncoding transcripts identified in wild type BMDM are classified in GENCODE/Ensembl as retained_intron noncoding RNAs, defined as alternatively spliced transcripts believed to contain intronic sequences relative to other coding transcripts of the same gene (Fig. 4a). Many of these intronic transcripts may be pieces of pre-mRNAs or excised introns that are targeted for degradation. Since we cannot rule out a function for them as regulatory RNAs, they were included in this analysis. 17% of the noncoding RNAs are classified as processed_transcript, a general term for a gene/transcript that lacks an open reading frame. Nine percent of the noncoding RNAs fell into the category of lincRNA, defined as long intergenic noncoding RNA. Four percent are classified as nonsense_mediated_decay, transcripts that contain sequences tagging them for destruction. While not specifically defined as noncoding in Ensembl, nonsense_mediated_decay transcripts could possibly have functions as noncoding RNAs, so were included in the analysis. Five percent are classified as antisense or “transcripts that overlap the genomic span of a protein-coding locus on the opposite strand”. Bidirectional_promoter lncRNAs, sense_intronic, and snoRNA comprised 1% or less of the noncoding transcripts. Small RNAs, such as snoRNAs, were included in the analysis but constitute an exceedingly small portion of the overall noncoding DE transcripts identified. The reason for their underrepresentation is likely because small RNAs were not selected for in the initial RNA isolation process or in the polyA tail selection step of the library preparation process. Therefore, almost the entirety of the noncoding transcripts identified are lncRNAs, but this was by study design.
During RH infection, we identified 155 noncoding transcripts that were of higher abundance, and 94 transcripts that were of lower abundance (Fig. 4b). In comparison, 70 noncoding transcripts were of higher abundance and 13 were of lower abundance during PTG infection. When comparing RH to PTG infection, 22 noncoding transcripts were of higher abundance and 11 were of lower abundance. These 33 lncRNAs (Fig. 4d) are prime candidates to determine the role of lncRNA in parasite strain specific responses during infection. 31 noncoding transcripts were shared between RH and PTG infection (Fig. 4c and e). These 31 transcripts are the most likely candidates to be important for infection since they are similarly regulated in a strain-independent manner. A full list of noncoding DE transcripts for the three comparisons in MyD88+/+ BMDM (RH vs uninfected, PTG vs uninfected, and RH vs PTG) can be found in Additional file 10. Using qRT-PCR, differential expression of 4 lncRNAs (mir17hg, D43Rik, Loc105, and Gm19434) strongly validated the RNA-seq results (Additional file 11).
While at least three lncRNAs (Ftx, Snhg5 and Snhg15) have known functions, most of the DE long noncoding transcripts we identified have unknown function. However, many lncRNAs are associated with immune-related protein coding genes. Ftx, which is more abundant during RH infection, is a well-studied lncRNA with roles in cancer and X-chromosome inactivation [41, 42]. Two lncRNAs more highly abundant during RH infection, Snhg15 and Snhg5, are host genes for snoRNA. With roles in cancer, they appear to function as molecular sponges for microRNAs [43,44,45,46]. The conserved mir17hg lncRNA is a host gene for the mir17 microRNA cluster and is more abundant during both RH and PTG infection. Mir17 microRNAs are known to have a role in regulating apoptosis during T. gondii infection [29, 47, 48]. Interestingly, two Siva1 intronic lncRNAs (Siva1–203 and Siva1–205) were more abundant during RH infection, but the Siva1 protein-coding gene, an apoptosis-inducing factor, was not a DE. Similarly, three Nfkb1 intronic lncRNAs (Nfkb1–208, Nfkb1–210, and Nfkb1–202) were more abundant during RH infection but not the Nfkb1 protein-coding gene itself. While possibly an artifact, the differential abundance of these intronic sequences from their protein-coding counterpoints possibly suggests a function of these RNA species separate from the protein-coding gene. STAT1–210 (nonsense mediated decay) and STAT5a-207 (processed transcript) are differentially expressed in RH versus PTG and RH versus uninfected but the protein-coding transcript for these immune genes were not regulated. Several immune-related protein-coding genes (Socs2, ifi44, ifi203, ifi213) were co-regulated with overlapping lncRNAs (Socs2–209, Ifi44–203, ifi203–207, and ifi213–203).
Sequencing of MyD88 KO mouse Transcriptome reveals few gene expression differences between MyD88+/+ and MyD88−/− BMDM during T. gondi infection
The breakdown of all DE transcripts by type was remarkably similar for infection of wild type and MyD88 KO mice (Fig. 3a and Fig. 5a). This was also true for the noncoding transcripts more specifically (Figs. 4a and Fig. 5d). In general, the overall number of protein-coding and noncoding DE transcripts for RH and PTG infection was similar for wild type and MyD88 KO BMDM (Figs. 3b, 4b, 5b, and e). Likewise, the number of DE transcripts shared between RH and PTG was similar for wild type and MyD88 KO macrophages (Figs. 3c, 4c, 5c, and f). A full list of the protein-coding and noncoding DE transcripts for MyD88 KO macrophages is provided in Additional file 6 and Additional file 10.
We next focused on direct comparisons between MyD88+/+ and MyD88−/− BMDM responses regarding protein-coding RNA (Fig. 6) and noncoding (Fig. 7) profiles. Only a few protein-coding and noncoding transcripts were differentially expressed between wild type and MyD88 KO BMDM (Figs. 6a and 7a). 10 protein-coding and 5 noncoding transcripts were differentially expressed between uninfected wild type and MyD88 KO macrophages. Similarly, 9 protein-coding and 7 noncoding transcripts were differentially expressed between RH infected wild type and MyD88 KO macrophages. Only 10 protein-coding and 12 noncoding transcripts were differentially expressed between PTG infected MyD88+/+ and MyD88−/− BMDM. The degree of overlap between these three comparison groups is shown in Fig. 6b (coding transcripts) and Fig. 7b (noncoding transcripts).
In general, there was substantial overlap between the DE transcripts for wild type and MyD88 KO macrophages ranging from ~ 20–50% depending on the infection treatment (coding transcripts, Fig. 6c; noncoding transcripts, Fig. 7c). Functional analysis of the similarities between RH-infected wild type and MyD88 KO macrophages revealed enrichment of genes for cell cycle, DNA recombination, DNA replication, DNA repair, and cell migration (Additional file 7). Interestingly, immune genes were not functionally enriched in both MyD88+/+ and MyD88−/− macrophages, despite differential expression of several immune genes such as ccl17, arg1, socs2, and ccl24 during RH infection of wild type and MyD88 KO cells. Regarding lncRNAs, greater abundance of mir17hg was common to both wild type and MyD88 KO infection by RH and PTG. Siva1–205 and Nfkb1–210 lncRNAs were of higher abundance during RH infection of wild type and Myd88 KO macrophages.
Several immune-related protein-coding genes were among the DE genes when comparing MyD88+/+ and MyD88−/− macrophages, and a subset of these are known to be relevant in the response to Toxoplasma. Transcripts encoding the macrophage receptor Marco and the cytokine Csf3 were of lower abundance in the MyD88 KO BMDM relative to wild type cells for both uninfected and PTG treatment (Fig. 6d). Transcripts for Ifng (encoding IFN-γ) were less abundant in MyD88−/− cells during PTG infection. IFN-γ is well known to be essential for survival of mice to acute T. gondii infection . Transcripts for Arg1 (encoding arginase that catalyzes arginine hydrolysis) were less abundant in uninfected MyD88 KO macrophages. Toxoplasma is an arginine auxotroph and upregulation of Arg1 in wild type cells is thought to be a host defense mechanism to reduce levels of arginine [35, 49]. Ifi208 (interferon activated gene 208) is more abundant in MyD88 KO BMDM compared to wild type during RH infection.
Less is known regarding the differentially expressed noncoding RNAs identified. However, some results stand out. Two lncRNAs, Gm42047–201 and Gm47601–201, were differentially expressed among all three conditions (MyD88 KO uninfected vs. wild type uninfected, MyD88 KO RH versus wild type RH, and MyD88 KO PTG versus wild type PTG) suggesting these two lncRNAs may be important in the difference between the two macrophage strains (Fig. 7d). 5 of the 18 noncoding transcripts (28%) are long intergenic noncoding RNAs meaning they are not near or overlapping any protein-coding genes. This is of practical interest because intergenic lncRNAs are more amenable to genetic knock down, making them easier targets to study. One lncRNA, Sltm-207, is intronic to the protein-coding gene SLTM and is less abundant in MyD88 KO compared to wild type during PTG infection. When overexpressed, SLTM acts a general inhibitor of transcription that eventually leads to apoptosis .
Sequencing of the Toxoplasma Transcriptome revealed no differences in gene expression between infection of wild type and MyD88 KO macrophages
Toxoplasma reads were mapped to the ToxoDB ME49 strain and differentially expressed genes were identified (p-value equal or less than 0.05 and fold change of 2 or greater). Because both Type I and Type II strains were mapped to ME49 (Type II), genetic differences between the strains could result in issues mapping to specific genes. While the overall mapping rate to RH was higher than for PTG, discrepancies for specific genes could exist. Additionally, the number of T. gondii reads obtained for many of the samples was below our projected target of 10 M reads per sample (actual range was between 1.2 M and 10 M reads per sample). This limited our ability to detect lowly expressed genes and could result in detecting fewer DE genes overall. Because we also did not analyze a non-infection control for T. gondii sequences (which is difficult to do), we conducted fewer comparisons than for the mouse transcriptome data. We directly compared RH and PTG infection for both wild type and MyD88 KO macrophages. We also compared wild type versus MyD88 KO, for RH and PTG, respectively.
For both wild type and MyD88 KO, approximately 120 parasite genes were more abundant in RH infection compared to PTG infection (Fig. 8a). Unexpectedly, a larger number of genes (approximately 440) were less abundant in RH compared to PTG parasites (Fig. 8a). The opposite was true when comparing mouse differential gene expression between RH and PTG (Fig. 3b). More than half (approximately 55%) of the parasite DE transcripts are hypothetical proteins with no known function or homology to other proteins with known functions (Additional file 12). Surprisingly, no differentially expressed Toxoplasma genes were identified when comparing MyD88 KO RH versus wild type RH and MyD88 KO PTG versus wild type PTG host strains (Fig. 8a). This indicates that Toxoplasma behaves the same regardless of the host macrophage MyD88 genotype. Most (approximately 80%) of the differentially expressed genes between RH and PTG were shared between wild type and MyD88 KO strains (Fig. 8b and Additional File 12). This high degree of similarity is consistent with the absence of significantly differentially expressed genes between T. gondii infection of wild type and MyD88 KO strains.
Pathway analysis identified metabolism of purines and several amino acids as different between RH and PTG infection (Fig. 8c-f). Overall, these transcripts were largely less abundant during RH infection compared to PTG. Gene ontology analysis identified microtubule-based movement as less abundant during RH infection compared to PTG infection (Fig. 8e and f). Several motor proteins (myosin, dynein, intraflagellar transport protein, and kinesin) were less abundant in RH versus PTG in both wild type and MyD88 KO (Additional file 12). However, two dynein proteins were more abundant in RH versus PTG for both wild type and MyD88 KO. T. gondii is known to use motor proteins for invasion into the host and for replication [51, 52]. Several SAG-related sequences (which encode a major family of tachyzoite surface proteins) were more abundant in RH versus PTG but a subset was also less abundant (Additional file 12). Many secretory proteins (ROP38, GRA11, ROP28, ROP2A, ROP23, ROP46, MIC17B, and ROP4) were less abundant in RH versus PTG in both wild type and MyD88 KO. GRA15 was more abundant in RH versus PTG but interestingly only during infection of MyD88−/− BMDM. The rhoptry protein ROP8 was less abundant in RH versus PTG only in MyD88 KO host cells. Oocyt wall protein (OWP1) was less abundant in RH versus PTG in both wild type and MyD88−/− macrophages. Transcripts for the bradyzoite antigen BAG1 were less abundant during RH infection compared to PTG infection in both wild type and MyD88 KO. This is consistent with the RH strain being unable to form cysts during infection. We validated the lower abundance of Bag1 in RH versus PTG using qRT-PCR (Additional File 13).
Independently of the DE analysis, we examined which Toxoplasma genes were most highly expressed during infection. Using normalized expression values, we compared the top 100 expressed genes between all four samples (RH infection of wild type and MyD88−/− cells; PTG infection of wild type and MyD88−/− BMDM). 75 out of 100 most expressed genes were shared between all four samples and these are listed in Fig. 9. The most highly expressed gene was that encoding rhoptry protein ROP1, a secretory molecule whose function is unclear. Many other T. gondii secreted proteins were represented, including dense granule proteins (8 in total), rhoptry proteins (2 in total), and microneme proteins (7 in total). Together, they comprise 22.7% of the most highly expressed genes. 23 out of 75 genes (30.7%) were ribosomal proteins. Only six out of the 75 genes (8%) were hypothetical, substantially less than the overall percentage of hypothetical proteins.
In this study, we undertook a global analysis of the coding and noncoding mouse transcriptome during Type I and Type II T. gondii infection of wild type and MyD88 KO BMDMs. Numerous previous studies have examined the host protein-coding response to infection using RNA-sequencing in various parasite strains, cells/organs, and host species, including bone-marrow derived macrophages and bone-marrow derived dendritic cells [49, 53,54,55,56,57,58,59]. While not the primary focus of our study, we identified many similarities to previous studies in the wild type mouse protein-coding response, including dysregulation of immune and metabolic genes. One of our more interesting findings is that significantly more host genes were differentially expressed during Type I rather than Type II infection. Our studies were conducted on cells infected for 6 h, prior to the first parasite mitotic event. Thus, strain-specific differences in replication do not account for the effect. This indicates that Type II infection is overall more silent than Type I infection, a finding that is consistent with in vivo evidence in mice that Type II parasites overall cause less severe disease than Type I strains [60, 61]. Of interest is to consider the transcriptomic changes during RH and PTG infection in light of the polymorphic effector molecules expressed by these strains. Type I parasite strains such as RH express an active form of ROP16. This rhoptry protein is a host-directed kinase that activates STAT3, 5 and 6 signaling pathways that down-regulate several immune response genes [35, 62]. In our results, we found that many immune genes were both more and less abundant during T. gondii infection, clearly indicating that other factors, in addition to ROP16, are relevant to macrophage infection. Conversely, PTG, a type II strain, expresses GRA15, a secreted dense granule protein that activates several proinflammatory genes through an NFκB pathway . This is perhaps in contradiction to our findings, where most host immune genes were only weakly affected by PTG infection. The influence of these parasite effector molecules on the transcriptional responses reported here await further investigation.
Another interesting finding was the higher abundance of cell migration genes in both the RH versus uninfected and RH versus PTG comparisons. This is consistent with data indicating that T. gondii enhances the motility of dendritic cells, a phenomenon hypothesized to facilitate dissemination throughout the host [18, 64]. Our results suggest that this phenomenon extends beyond dendritic cells to macrophages as well.
In this study, the parasites were centrifuged and washed prior to infecting BMDM, but parasite preparations likely contain small amounts of fibroblast debris. We cannot rule out that the presence of some fibroblast debris could possibly have their own effects on gene expression. If so, this would most likely apply to mouse genes differentially expressed during both RH and PTG infection and less likely to affect only one strain, as the two strains were treated identically with regard to centrifugation and washing. The possibility that fibroblast factors rather than parasites per se mediate some of the effects seen is minimized by the fact that the fibroblasts are human in origin whereas the macrophages used for infection are mouse derived.
While many studies have assessed the protein-coding transcriptomic response to T. gondii infection, very few have looked at the noncoding RNA response to infection. Two previous studies examined expression of lncRNAs in human cells during T. gondii infection. One, using microarrays, found that at least 1000 human lncRNAs were differentially expressed during Type II infection of HFF cells . Another used qRT-PCR to assess a panel of lncRNAs with known immune regulatory activity and found 31 lncRNAs were differentially expressed during T. gondii infection of retinal müller glial cells . Despite the low degree of sequence conservation between human and mouse lncRNAs, two lncRNAs, Mir17hg and Snhg15, differentially expressed in the Rochet et al. study were also differentially expressed in this study.
We previously looked at the expression of lncRNAs in mouse cells during T. gondii infection using microarrays . In that study, we found 1522 lncRNAs differentially expressed during infection with the RH strain and 528 with the PTG strain. Here, we found 249 noncoding transcripts as differentially expressed during Type I infection and 83 noncoding transcripts as differentially expressed during Type II infection of wild type mouse macrophages. While we identified a smaller number of lncRNAs with RNA-sequencing than with microarrays, we found significant overlap including Mir17hg, Ftx, and Socs2–209. While difficult to directly compare, we identified approximately 55 lncRNAs shared between the two datasets based on gene names. We identified a larger number (155) of shared DE transcripts for mRNAs, possibly due to either more consistent naming or higher expression of protein-coding genes.
In this study, we identified 33 lncRNAs that were differentially expressed between RH and PTG. These would be important to pursue for their role in strain specific differences between Type I and Type II T. gondii. Additionally, 31 noncoding transcripts were differentially expressed in both RH and PTG infection. These 31 transcripts are perhaps most likely to be essential for infection as they are regulated in both strains. This study, combined with previous studies, provides a sense of the vast landscape of lncRNAs regulated during T. gondii infection whose functions we do not yet understand. Given the large number of lncRNAs identified in this study, as well as the emerging role of lncRNAs in regulating the immune response [30, 68], we believe this is an important avenue to investigate further in order to fully understand the immune response to T. gondii and other intracellular pathogens impacting human health. Functional studies employing approaches such as knockout or overexpression of lncRNAs is a necessary next step.
The present study is the first to examine the transcriptional response of MyD88 deficient mouse cells during T. gondii infection. The MyD88 protein plays an important role in the early innate immune response in mice to T. gondii infection . The effect of MyD88 deletion on expression of downstream immune genes in infected cells was previously unknown. Additionally, this is an important avenue of research because the human immune response to Toxoplasma is not thought to depend on the TLR/MyD88 axis, and there is much interest in determining other avenues of immune recognition in humans and mice. Perhaps surprisingly, only a few genes were differentially expressed between wild type and MyD88 KO mice. This suggests that MyD88 has only minor effects on cellular transcriptional activity, at least within the time frame of our study. Several immune genes and two meiotic genes were among the differences. Two lncRNAs (Gm42047–201 and Gm47601–201) were differentially expressed among all three conditions (RH, PTG, and noninfected), suggesting that these two lncRNAs may be especially important in the MyD88-dependent immune response. Further investigation of these protein-coding genes and non-coding transcripts in the context of the MyD88 KO immune response is important.
While not the primary focus of our study, we also examined parasite gene expression for RH and PTG infection of wild type and MyD88 KO macrophages. More than 550 genes were differentially expressed between Type I and Type II T. gondii. In both the RNA-sequencing dataset and qRT-PCR validation, the Bag1 gene was less abundant in the RH strain versus the PTG strain. This result is unsurprising as Bag1 expression is specific to the bradyzoite stage of the parasite, a life cycle stage which the RH strain is incapable of achieving. Several microtubule and motor proteins (myosin, dynein, intraflagellar transport protein, and kinesin) were less abundant in RH versus PTG in both wild type and MyD88 KO. This result is surprising given the much faster growth rate of the RH strain, however, two dynein proteins were more abundant during RH infection. It is interesting that the most highly expressed Toxoplasma gene across samples was the rhoptry protein ROP1, a secretory molecule whose function is unknown. Determining the function of this protein could be interesting and important. Previous studies have used RNA-sequencing to examine T. gondii gene expression patterns [54,55,56, 58, 70]. In 2014, Pittman et al. compared acute and chronic Type II T. gondii infection of mouse brains. While difficult to directly compare the two studies, we also found that metabolic, biosynthetic, and translation genes were among the most highly expressed T. gondii genes regardless of strain type. In 2016, Zhou et al. examined RH gene expression during infection of porcine PK-15 cells. They found that metabolic genes were largely less abundant in intracellular RH tachyzoites compared to egressed parasites. We also found that many metabolic genes were less abundant in RH when compared to PTG infection. One of our more interesting findings is that no T. gondii genes were differentially expressed between wild type and MyD88 KO macrophages, suggesting that T. gondii behaves the same regardless of which host strain it infects. Nevertheless, it is possible that at later infection time points differences in the parasite transcriptome in wild type and KO cells would emerge.
While this study covered numerous disparate aspects of early intracellular infection, we were largely focused on two facets: (1) identifying long noncoding RNAs differentially expressed during T. gondii infection and (2) examining the transcriptional response of MyD88 KO cells. Host lncRNAs represent a largely unexplored layer in the immune response to T. gondii. The MyD88-independent transcriptional response may also be important in unraveling other mechanisms of immune regulation in the response to T. gondii in both humans and mice.
Mouse strains and generation of BMDM
C57BL/6 mice (The Jackson Laboratory) were used as the wild type strain in this study. The MyD88 KO mice contain a deletion of exon 3 (The Jackson Laboratory, B6.129P2(SJL)-Myd88tm1.1Defr/J). Animals were housed (3–5 per cage) and cared for according to the Guide for the Care and Use of Laboratory Animals (8th edition) under an IACUC approved protocol (19–200,854-MC). Mice were housed in 3 single-sided ventilated cage racks (RAIR HD Super Mouse 1800 Interchangeable Micro Isolator Unit 9 tier). The mice were housed in a specific pathogen-free facility under a 12:12 light/dark cycle with controlled temperature (22–24 °C) and humidity (50–60%). Mice were given free access to autoclaved standard mouse chow and water and were checked every 24 h for general health. Animals were rested for 1 week prior to generating BMDM. We used a total of 8 wild type mice and 8 MyD88 KO mice in this study. All mice were females between the ages of 6 and 8 weeks (17–20 g). The mice were euthanized using CO2 asphyxiation, as recommended by the American Veterinary Medical Association. Mice were euthanized in the morning in a separate room from caged mice. Bone marrow from a single mouse was used for each biological replicate for each mouse strain, with four biological replicates from four separate mice shown (4 replicates for uninfected mice samples and 3 replicates for infected samples). Femur and tibia were used as a source of bone marrow cells. Macrophages were generated from single cell suspensions of bone marrow cells by 5- or 6-day culture in L929-containing supernatants, as previously described . One day prior to infection, BMDM were harvested, counted, and plated in 12-well tissue culture plates in human fibroblast medium (HFM) at a concentration of 1 × 106 cells per well. HFM consists of DMEM (Life Technologies) supplemented with 1% heat-inactivated bovine growth serum (Thermo Fisher Scientific), 100 U/mL penicillin (Life Technologies), and 0.1 mg/mL streptomycin (Life Technologies).
Parasite strains and infections
Wild type Toxoplasma strains RH (Type I) and PTG (Type II) were used in this study. Tachyzoites of both strains were maintained by approximately twice-weekly passage on human foreskin fibroblast monolayers in HFM. Flasks of fully lysed fibroblasts were utilized for infections. If fibroblasts were only partially lysed, they were passed twice through a 27 Ga needle to release the parasites. All parasites were centrifuged and washed prior to infection to remove soluble material. Infections were accomplished by the addition of tachyzoites to mouse BMDM (4:1 or 5:1 ratio of parasites to cells) on 12-well plates (Falcon, non-tissue culture treated). Plates were briefly centrifuged (3 min, 200×g) to initiate contact between tachyzoites and macrophages. Cultures were incubated 6 h (37 °C, 5% CO2), then the cells were harvested for RNA extraction. Percent infection was calculated using an Olympus BX51 immunofluorescence microscope. DAPI was used to stain the nucleus. Toxoplasma-specific polyclonal antibody conjugated to FITC (ThermoFisher Scientific) was used to stain the parasites. Percent infection of at least 70% was deemed acceptable.
RNA-seq sample preparation and Illumina NextSeq 500 sequencing
Total RNA was prepared from BMDM by RNeasy Mini Kit purification (Qiagen). The samples were subjected to Turbo DNase treatment (Life Technologies). Total RNA from each sample was quantified using a Qubit 3.0 Fluorometer. At least 1.5 μg per sample was submitted to the CETI Molecular Biology Core at the University of New Mexico, and RNA integrity was assessed using an Agilent 2100 Bioanalyzer. All samples had RIN values > 9.6, indicating the RNA was high quality. Libraries were prepared using the KAPA mRNA HyperPrep Kit with PolyA selection. In brief, mRNA was captured using magnetic oligo-dT beads. RNA was fragmented using heat and magnesium, and first strand was synthesized by random priming. Second strand synthesis and A-tailing were combined into one step, which involves converting the cDNA:RNA hybrid to double-stranded cDNA while incorporating dUTP into the second cDNA stand and adding dAMP to the 3′ ends of the double-stranded cDNA fragments. Next, dsDNA adapters with 3′-dTMP were ligated to the A-tailed library fragments, and the library fragments with adapters at both ends were amplified using PCR. The strand marked with dUTP was not amplified, allowing strand-specific sequencing. Library quality was assessed using a Qubit and a Bioanalyzer. qRT-PCR of adapters was performed to determine the amount of library in each sample. Libraries were sequenced using an Illumina NextSeq 500 Sequencer generating paired-end, 75-bp read length sequences. 1 mid-throughput run (150 M theoretical maximum clusters) and 3 high-throughput runs (400 M theoretical maximum clusters) were performed. Figure 1 indicates the total number of reads obtained per sample.
RNA-seq data analysis
Raw reads off sequencer were trimmed and filtered using Trimmomatic v0.36 (Bolger et al. 2014) with slide window of 4 nt, average score above 20 and minimum length of 36 nt. Filtered high quality reads were mapped to the annotated mouse genome release M21 from GenCode  and T. gondii ME49 (GenBank Assembly ID GCA_000006565.2) using RNA-Seq reads mapping tool Spliced Transcripts Alignment to a Reference (STAR) 2.5.3a . Reads mapping were performed using STAR options: --genomeLoad LoadAndKeep --limitBAMsortRAM 50,000,000,000 --runMode alignReads --runThreadN 17 --outBAMsortingThreadN 7. No major changes were made to default settings other than to tune CPUs and memory. Gene expression levels were estimated using software featureCounts with default settings . Differential gene expression (DE) analysis was performed using EBSeq v1.22.0 default settings . A posterior probability of differential expression (PPDE) 0.95 or higher for EBSeq was set as cutoff for DE analysis. Workflow of read trimming and mapping was built using Unix shell commands with application GNU-Parallel  to perform jobs in parallel. Table summary and figure generation were performed with the R statistical computing environment  and Bioconductor , including the following packages: ggplot2 , reshape2 , and pheatmap v1.0.12 . Heat maps were constructed using complete-linkage clustering using log2 fold change values or the z-score transformed FPKM values. PCA plots were constructed using SARTools v1.6 . For mouse genes, functional enrichment analysis, including KEGG pathway analysis & Gene Ontology analysis, was performed using g:Profiler (g:Profiler version e96_eg43_p13_563554d) with g:SCS multiple testing correction method applying significance threshold of 0.05 . For T. gondii genes, functional enrichment analysis was performed using g:Profiler (for gene ontology data) and ToxoDB Identify Metabolic Pathway tool (for KEGG pathway data) . Venn diagrams, pie charts, and other data analysis was performed using Microsoft Excel software.
Total RNA was prepared using a RNeasy Mini Kit (Qiagen), and the samples were subjected to Turbo DNase treatment (Life Technologies). RNA was converted to cDNA using the SuperScript IV VILO Master Mix (ThermoFisher Scientific). Quantitative PCR was performed on target genes and normalized to the expression of the housekeeping gene Ppia (mouse gene) or β-tubulin (T. gondii gene) using the SYBR green method (SsoAdvanced Universal SYBR Green Supermix, Bio-Rad) and the Bio-Rad CFX96 RT-PCR machine. Expression relative to uninfected control samples was calculated using the ∆∆Ct method. A control with no added reverse transcriptase was included for each sample. Primer sequences used in this study are Ppia (F-GCATGTGGTCTTTGGGAAGGTG and R-GGGTAAAATGCCCGCAAGTCAA), Mir17hg (F-GGTGGCCACTCTGTTAATGTGC and R-TAACTGCAGCTTCTCCAGACCC), D43Rik (F-TTCTCAATACAACGCCCCAGGT and R-AAAAGGGGGAAGGATTGGGGAG), LOC105 (F-GGCACACCCATGATAGGCTGTA and R-CCTCCATACCCAGCACTGTCAA), Gm19434 (F-TGACAGGGAAAAACAGAGCCCA and R-GTCACCCAGGGGAAAGCAATTG), β-tubulin (F-CGCCACGGCCGCTACCTGACT and R-TACGCGCCTTCCTCTGCACCC), and Bag1 (F-GACGTGGAGTTCGACAGCAAA and R-ATGGCTCCGTTGTCGACTTCT). Some of our uninfected samples were occasionally at the limit of detection of our qRT-PCR assay. In those cases, we used 35 cycles as the Cq (quantification cycle).
Availability of data and materials
The datasets generated and/or analysed during the current study are available in the NCBI Sequence Read Archive (SRA) repository, PRJNA695200 at https://www.ncbi.nlm.nih.gov/bioproject/695200. Reads were mapped to the annotated mouse genome release M21 from GenCode (GenBank assembly ID GCA_000001635.2) and T. gondii ME49 (GenBank Assembly ID GCA_000006565.2).
Bone marrow-derived macrophages
Encyclopedia of DNA elements
Human fibroblast medium
Kyoto encyclopedia of genes and genomes
long intergenic noncoding RNA
long noncoding RNA
Principal component analysis
Posterior probability of differential expression
quantitative reverse transcriptase polymerase chain reaction
small nucleolar RNA
To be experimentally confirmed
Robert-Gangneux F, Darde ML. Epidemiology of and diagnostic strategies for toxoplasmosis. Clin Microbiol Rev. 2012;25(2):264–96.
Dubey JP. The history and life-cycle of Toxoplasma gondii. In: Weiss LM, Kim K, editors. Toxoplasma gondii the model apicomplexan: perspective and methods. 2nd ed. San Diego: Academic Press; 2013. p. 1–17.
Weiss LM, Kim K. The development and biology of bradyzoites of toxoplasma gondii. Front Biosci. 2000;5(1):D391–405.
McLeod RM, van Tubbergen C, Montoya JG, Petersen E. Human toxoplasma infection. In: Weiss LM, Kim K, editors. Toxoplasma gondii the model Apicomplexan: perspectives and methods. 2nd ed. San Diego: Academic Press; 2013. p. 100–59.
Pfaff AW, Liesenfeld O, Candolfi E. Congenital toxoplasmosis. In: Ajioka JW, Soldati D, editors. Toxoplasma molecular and cellular biology. Norfolk: Horizon Bioscience; 2007. p. 93–110.
Dupont CD, Christian DA, Hunter CA. Immune response and immunopathology during toxoplasmosis. Semin Immunopathol. 2012;34(6):793–813.
Sasai M, Pradipta A, Yamamoto M. Host immune responses to toxoplasma gondii. Int Immunol. 2018;30(3):113–9.
Scharton-Kersten TM, Wynn TA, Denkers EY, Bala S, Grunvald E, Hieny S, et al. In the absence of endogenous IFN-gamma, mice develop unimpaired IL-12 responses to toxoplasma gondii while failing to control acute infection. J Immunol. 1996;157(9):4045–54.
Gazzinelli RT, Wysocka M, Hieny S, Scharton-Kersten T, Cheever A, Kuhn R, et al. In the absence of endogenous IL-10, mice acutely infected with Toxoplasma gondii succumb to a lethal immune response dependent upon CD4+ T cells and accompanied by overproduction of IL-12, IFN-g, and TNF-a. J Immunol. 1996;157(2):798–805.
Suzuki Y, Sher A, Yap G, Park D, Ellis Neyer L, Liesenfeld O, et al. IL-10 is required for prevention of necrosis in the small intestine and mortality in both genetically resistant BALB/c and susceptible C57BL/6 mice following peroral infection with Toxoplasma gondii. J Immunol. 2000;164:5375–82.
Cohen SB, Maurer KJ, Egan CE, Oghumu S, Satoskar AR, Denkers EY. CXCR3-dependent CD4(+) T cells are required to activate inflammatory monocytes for defense against intestinal infection. PLoS Pathog. 2013;9(10):e1003706.
Dunay IR, Damatta RA, Fux B, Presti R, Greco S, Colonna M, et al. Gr1(+) inflammatory monocytes are required for mucosal resistance to the pathogen toxoplasma gondii. Immunity. 2008;29(2):306–17.
Dunay IR, Sibley LD. Monocytes mediate mucosal immunity to toxoplasma gondii. Curr Opin Immunol. 2010;22(4):461–6.
Bierly AL, Shufesky WJ, Sukhumavasi W, Morelli A, Denkers EY. Dendritic cells expressing plasmacytoid marker PDCA-1 are Trojan horses during toxoplasma gondii infection. J Immunol. 2008;181(12):8445–91.
Caamano J, Alexander J, Craig L, Bravo R, Hunter CA. The NF-kB family member RelB is required for innate and adaptive immunity to Toxoplasma gondii. J Immunol. 1999;163(8):4453–61.
Cohen SB, Denkers EY. Impact of toxoplasma gondii on dendritic cell subset function in the intestinal mucosa. J Immunol. 2015;195(6):2754–62.
Courret N, Darche S, Sonigo P, Milon G, Buzoni-Gatel D, Tardieux I. CD11c and CD11b expressing mouse leukocytes transport single Toxoplasma gondii tachyzoites to the brain. Blood. 2006;107(1):309–16.
Sangare LO, Olafsson EB, Wang Y, Yang N, Julien L, Camejo A, et al. In vivo CRISPR screen identifies TgWIP as a toxoplasma modulator of dendritic cell migration. Cell Host Microbe. 2019;26(4):478–92 e8.
Egan CE, Sukhumavasi W, Butcher BA, Denkers EY. Functional aspects of toll-like receptor/MyD88 signalling during protozoan infection: focus on toxoplasma gondii. Clin Exp Immunol. 2009.
Pifer R, Yarovinsky F. Innate responses to toxoplasma gondii in mice and humans. Trends Parasitol. 2011;27:388–93.
Andrade WA, Souza MD, Ramos-Martinez E, Nagpal K, Dutra MS, Melo MB, et al. Combined action of nucleic acid-sensing toll-like receptors and TLR11/TLR12 heterodimers imparts resistance to toxoplasma gondii in mice. Cell Host Microbe. 2013;13(1):42–53.
Raetz M, Kibardin A, Sturge CR, Pifer R, Li H, Burstein E, et al. Cooperation of TLR12 and TLR11 in the IRF8-dependent IL-12 response to toxoplasma gondii profilin. J Immunol. 2013;191(9):4818–27.
Yarovinsky F, Zhang D, Anderson JF, Bannenberg GL, Serhan CN, Hayden MS, et al. TLR11 activation of dendritic cells by a protozoan profilin-like protein. Science. 2005;308(5728):1626–9.
Plattner F, Yarovinsky F, Romero S, Didry D, Carlier MF, Sher A, et al. Toxoplasma profilin is essential for host cell invasion and TLR11-dependent induction of an IL-12 response. Cell Host Microbe. 2008;3:1477–87.
Gazzinelli RT, Mendonca-Neto R, Lilue J, Howard J, Sher A. Innate resistance against toxoplasma gondii: an evolutionary tale of mice, cats, and men. Cell Host Microbe. 2014;15(2):132–8.
von Bernuth H, Picard C, Jin Z, Pankla R, Xiao H, Ku CL, et al. Pyogenic bacterial infections in humans with MyD88 deficiency. Science. 2008;321(5889):691–6.
von Bernuth H, Picard C, Puel A, Casanova JL. Experimental and natural infections in MyD88- and IRAK-4-deficient mice and humans. Eur J Immunol. 2012;42(12):3126–35.
Agliano F, Rathinam VA, Medvedev AE, Vanaja SK, Vella AT. Long noncoding RNAs in host-pathogen interactions. Trends Immunol. 2019;40(6):492–510.
Menard KL, Haskins BE, Denkers EY. Impact of toxoplasma gondii infection on host non-coding RNA responses. Front Cell Infect Microbiol. 2019;9(1):132.
Atianand MK, Caffrey DR, Fitzgerald KA. Immunobiology of long noncoding RNAs. Annu Rev Immunol. 2017;35(1):177–98.
Chen YG, Satpathy AT, Chang HY. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol. 2017;18(9):962–72.
Araujo FG, Slifer T. Different strains of Toxoplasma gondii induce different cytokine responses in CBA/Ca mice. Infect Immun. 2003;71:4171–4.
Gavrilescu LC, Denkers EY. IFN-gamma overproduction and high level apoptosis are associated with high but not low virulence toxoplasma gondii infection. J Immunol. 2001;167(2):902–9.
Mordue DG, Monroy F, La Regina M, Dinarello CA, Sibley LD. Acute toxoplasmosis leads to lethal overproduction of Th1 cytokines. J Immunol. 2001;167(8):4574–84.
Butcher BA, Fox BA, Rommereim LM, Kim SG, Maurer KJ, Yarovinsky F, et al. Toxoplasma gondii rhoptry kinase ROP16 activates STAT3 and STAT6 resulting in cytokine inhibition and arginase-1-dependent growth control. PLoS Pathog. 2011;7(9):e1002236.
Jensen KD, Wang Y, Wojno ED, Shastri AJ, Hu K, Cornel L, et al. Toxoplasma polymorphic effectors determine macrophage polarization and intestinal inflammation. Cell Host Microbe. 2011;9(6):472–83.
Robben PM, Mordue DG, Truscott SM, Takeda K, Akira S, Sibley LD. Production of IL-12 by macrophages infected with toxoplasma gondii depends on the parasite genotype. J Immunol. 2004;172(6):3686–94.
Patil V, Zhao Y, Shah S, Fox BA, Rommereim LM, Bzik DJ, et al. Co-existence of classical and alternative activation programs in macrophages responding to toxoplasma gondii. Int J Parasitol. 2014;44(2):161–4.
Melo EJ, Carvalho TM, De Souza W. Behaviour of microtubules in cells infected with toxoplasma gondii. Biocell. 2001;25(1):53–9.
Walker ME, Hjort EE, Smith SS, Tripathi A, Hornick JE, Hinchcliffe EH, et al. Toxoplasma gondii actively remodels the microtubule network in host cells. Microbes Infect. 2008;10(14–15):1440–9.
Hosoi Y, Soma M, Shiura H, Sado T, Hasuwa H, Abe K, et al. Female mice lacking Ftx lncRNA exhibit impaired X-chromosome inactivation and a microphthalmia-like phenotype. Nat Commun. 2018;9(1):3829.
Yang Y, Zhang J, Chen X, Xu X, Cao G, Li H, et al. LncRNA FTX sponges miR-215 and inhibits phosphorylation of vimentin for promoting colorectal cancer progression. Gene Ther. 2018;25(5):321–30.
Li Y, Guo D, Zhao Y, Ren M, Lu G, Wang Y, et al. Long non-coding RNA SNHG5 promotes human hepatocellular carcinoma progression by regulating miR-26a-5p/GSK3beta signal pathway. Cell Death Dis. 2018;9(9):888.
Tong J, Ma X, Yu H, Yang J. SNHG15: a promising cancer-related long noncoding RNA. Cancer Manag Res. 2019;11(1):5961–9.
Zanovello P, Vallerani E, Biasi G, Landolfo S, Colavo D. Monoclonal antibody against IFN-g inhibits Maloney murine sarcoma virus-specific cytotoxic T cell differentiation. J Immunol. 1988;140(4):1341–4.
Zhang Y, Zhang D, Lv J, Wang S, Zhang Q. LncRNA SNHG15 acts as an oncogene in prostate cancer by regulating miR-338-3p/FKBP1A axis. Gene. 2019;705(1):44–50.
Cai Y, Chen H, Jin L, You Y, Shen J. STAT3-dependent transactivation of miRNA genes following toxoplasma gondii infection in macrophage. Parasit Vectors. 2013;6(1):356.
Rezaei F, Daryani A, Sharifi M, Sarvi S, Jafari N, Pagheh AS, et al. miR-20a inhibition using locked nucleic acid (LNA) technology and its effects on apoptosis of human macrophages infected by toxoplasma gondii RH strain. Microb Pathog. 2018;121(1):269–76.
Hargrave KE, Woods S, Millington O, Chalmers S, Westrop GD, Roberts CW. Multi-Omics studies demonstrate toxoplasma gondii-induced metabolic reprogramming of murine dendritic cells. Front Cell Infect Microbiol. 2019;9(1):309.
Chan CW, Lee YB, Uney J, Flynn A, Tobias JH, Norman M. A novel member of the SAF (scaffold attachment factor)-box protein family inhibits gene expression and induces apoptosis. Biochem J. 2007;407(3):355–62.
Carruthers V, Boothroyd JC. Pulling together: an integrated model of toxoplasma cell invasion. Curr Opin Microbiol. 2006;10(1):83–9.
Del Rosario M, Periz J, Pavlou G, Lyth O, Latorre-Barragan F, Das S, et al. Apicomplexan F-actin is required for efficient nuclear entry during host cell invasion. EMBO Rep. 2019;20(12):e48896.
He JJ, Ma J, Wang JL, Zhang FK, Li JX, Zhai BT, et al. Global Transcriptome profiling of multiple porcine organs reveals toxoplasma gondii-induced transcriptional landscapes. Front Immunol. 2019;10(1):1531.
Melo MB, Nguyen QP, Cordeiro C, Hassan MA, Yang N, McKell R, et al. Transcriptional analysis of murine macrophages infected with different toxoplasma strains identifies novel regulation of host signaling pathways. PLoS Pathog. 2013;9(12):e1003779.
Olson WJ, Martorelli Di Genova B, Gallego-Lopez G, Dawson AR, Stevenson D, Amador-Noguez D, et al. Dual metabolomic profiling uncovers toxoplasma manipulation of the host metabolome and the discovery of a novel parasite metabolic capability. PLoS Pathog. 2020;16(4):e1008432.
Pittman KJ, Aliota MT, Knoll LJ. Dual transcriptional profiling of mice and toxoplasma gondii during acute and chronic infection. BMC Genomics. 2014;15(1):806.
Tanaka S, Nishimura M, Ihara F, Yamagishi J, Suzuki Y, Nishikawa Y. Transcriptome analysis of mouse brain infected with toxoplasma gondii. Infect Immun. 2013;81(10):3609–19.
Zhou CX, Elsheikha HM, Zhou DH, Liu Q, Zhu XQ, Suo X. Dual identification and analysis of differentially expressed transcripts of porcine PK-15 cells and toxoplasma gondii during in vitro infection. Front Microbiol. 2016;7(1):721.
Zhou CX, Zhou DH, Liu GX, Suo X, Zhu XQ. Transcriptomic analysis of porcine PBMCs infected with toxoplasma gondii RH strain. Acta Trop. 2016;154:82–8.
Gavrilescu LC, Denkers EY. IFN-g overproduction and high level apoptosis are associated with high but not low virulence Toxoplasma gondii infection. J Immunol. 2001;167:902–9.
Mordue DG, Monroy F, La Regina M, Dinarello CA, Sibley LD. Acute toxoplasmosis leads to lethal overproduction of Th1 cytokines. J Immunol. 2001;167:4574–84.
Saeij JP, Coller S, Boyle JP, Jerome ME, White MW, Boothroyd JC. Toxoplasma co-opts host gene expression by injection of a polymorphic kinase homologue. Nature. 2007;445(7125):324–7.
Rosowski EE, Lu D, Julien L, Rodda L, Gaiser RA, Jensen KD, et al. Strain-specific activation of the NF-kappaB pathway by GRA15, a novel toxoplasma gondii dense granule protein. J Exp Med. 2011;208(1):195–212.
Fuks JM, Arrighi RB, Weidner JM, Kumar Mendu S, Jin Z, Wallin RP, et al. GABAergic signaling is linked to a Hypermigratory phenotype in dendritic cells infected by toxoplasma gondii. PLoS Pathog. 2012;8(12):e1003051.
Liu W, Huang L, Wei Q, Zhang Y, Zhang S, Zhang W, et al. Microarray analysis of long non-coding RNA expression profiles uncovers a toxoplasma-induced negative regulation of host immune signaling. Parasit Vectors. 2018;11(1):174.
Rochet E, Appukuttan B, Ma Y, Ashander LM, Smith JR. Expression of long non-coding RNAs by human retinal Muller glial cells infected with clonal and exotic virulent toxoplasma gondii. Non-coding RNA. 2019;5(4):48.
Menard KL, Haskins BE, Colombo AP, Denkers EY. Toxoplasma gondii manipulates expression of host long noncoding RNA during intracellular infection. Sci Rep. 2018;8(1):15017.
Fitzgerald KA, Caffrey DR. Long noncoding RNAs in innate and adaptive immunity. Curr Opin Immunol. 2014;26(1):140–6.
Scanga CA, Aliberti J, Jankovic D, Tilloy F, Bennouna S, Denkers EY, et al. Cutting edge: MyD88 is required for resistance to Toxoplasma gondii infection and regulates parasite-induced IL-12 production by dendritic cells. J Immunol. 2002;168:5997–6001.
Croken MM, Ma Y, Markillie LM, Taylor RC, Orr G, Weiss LM, et al. Distinct strains of toxoplasma gondii feature divergent transcriptomes regardless of developmental stage. PLoS One. 2014;9(11):e111297.
Kim L, Butcher BA, Denkers EY. Toxoplasma gondii interferes with lipopolysaccharide-induced mitogen-activated protein kinase activation by mechanisms distinct from endotoxin tolerance. J Immunol. 2004;172(5):3003–10.
Frankish A, Diekhans M, Ferreira AM, Johnson R, Jungreis I, Loveland J, et al. GENCODE reference annotation for the human and mouse genomes. Nucleic Acids Res. 2019;47(D1):D766–D73.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
Leng N, Dawson JA, Thomson JA, Ruotti V, Rissman AI, Smits BM, et al. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 2013;29(8):1035–43.
Tange O. GNU parallel: the command-line power. USENIX Mag. 2011;36(1):42–7.
Team RC. R: a language and environment for statistical computing; 2014.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5(10):R80.
Wickham H. ggplot: Elegant Graphics for Data Analysis Springer; 2016.
Wickham H. reshape2: Flexibly reshape data: a reboot of the reshape package. R Packag. Version. ed; 2012.
Kolde R. Pheatmap: pretty heatmaps. R Packag. Version 61, 915 ed; 2012.
Varet H, Brillet-Gueguen L, Coppee JY, Dillies MA. SARTools: a DESeq2- and EdgeR-based R pipeline for comprehensive differential analysis of RNA-Seq data. PLoS One. 2016;11(6):e0157022.
Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, et al. G:profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019;47(W1):W191–W8.
Gajria B, Bahl A, Brestelli J, Dommer J, Fischer S, Gao X, et al. ToxoDB: an integrated toxoplasma gondii database resource. Nucleic Acids Res. 2008;36(Database issue):D553–6.
Susan Carpenter encouraged us to perform sequencing and provided guidance on best practices. Breanne Haskins performed qRT-PCR for lncRNAs. Anthony Colombo assisted with quantitative analysis and was helpful for discussions. George Rosenburg and Melissa Sanchez, managers in the UNM CETI Molecular Biology Core Facility, assisted with performing the sequencing.
This work was supported by grants from the National Institute of Allergy and Infectious Disease (AI139628) and the National Institute of General Medical Sciences (GM110907). The funders had no role in design, execution or interpretation of the study.
Ethics approval and consent to participate
All experimental protocols were approved by the University of New Mexico Institutional Animal Care and Use Committee (IACUC). The methods were carried out in accordance with the Panel on Euthanasia of the American Veterinary Medical Association.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Principal component analysis (PCA) plots of mouse RNA-sequencing data reveal associations between samples: PCA plots of mouse transcripts for wild type (wt) BMDM comparisons.
Principal component analysis (PCA) plots of mouse RNA-sequencing data reveal associations between samples: PCA plots of mouse transcripts for MyD88 KO (k88) BMDM comparisons.
Principal component analysis (PCA) plots of mouse RNA-sequencing data reveal associations between samples: PCA plots of mouse transcripts for MyD88 KO BMDM versus wild type BMDM comparisons.
Principal component analysis (PCA) plots of T. gondii RNA-sequencing data reveal associations between samples: PCA plots of T. gondii transcripts for MyD88 KO (k88) BMDM versus wild type (wt) BMDM comparisons.
Principal component analysis (PCA) plots of T. gondii RNA-sequencing data reveal associations between samples: PCA plots of T. gondii transcripts for RH versus PTG comparisons.
Complete list of differentially regulated mouse protein-coding genes in MyD88+/+ and MyD88−/− BMDM. Each tab in the spreadsheet corresponds to a specific pairwise comparison (e.g., RH vs. uninfected). Transcript name, Gene IDs, raw values, normalized values, PPDE values, fold changes, and other information are included.
Functional analysis of differentially expressed protein-coding genes. Each tab in the spreadsheet corresponds to a specific pairwise comparison (e.g., RH vs. uninfected). Results were obtained using g:Profiler. Gene Ontology: Biological Process and KEGG Pathway analysis results are included.
Heatmaps displaying individual replicates of differentially expressed immune and cell cycle genes in wildtype mice. Heatmaps displaying z-score transformed FPKM values of individual replicates among functionally related genes.
Heatmaps displaying individual replicates of differentially expressed growth, cell adhesion, and migration genes in wildtype mice. Heatmaps displaying z-score transformed FPKM values of individual replicates among functionally related genes.
Complete list of differentially expressed mouse non-coding genes. Each tab in the spreadsheet corresponds to a specific pairwise comparison (e.g., RH vs. uninfected). Transcript name, Gene IDs, raw values, normalized values, PPDE values, fold changes, and other information are included.
RNA-seq lncRNA data validation by qRT-PCR. RNA from mouse BMDM were collected 6 h after infection with either RH or PTG strains of T. gondii, and qRT-PCR was performed. Fold changes represent the comparison of infected samples to uninfected samples. (a) qRT-PCR fold change values for 4 lncRNAs. (b) RNA-seq fold change values for the same 4 lncRNAs for comparison purposes. Experiments were completed at least three times with BMDM from three separate mice.
Complete list of differentially expressed Toxoplasma gondii genes. The first tab corresponds to wild type RH versus wild type PTG. The second tab corresponds to Myd88KO RH versus Myd88KO PTG. The third tab lists the differentially regulated T. gondii genes common to both comparisons. Transcript ID, raw values, normalized values, PPDE values, fold changes, and other information are included.
RNA-seq Toxoplasma gondii data validation by qRT-PCR. RNA from mouse BMDM were collected 6 h after infection with either RH or PTG strains of T. gondii, and qRT-PCR was performed. Fold changes represent the comparison of RH infected samples to PTG infected samples. RNA-sequencing and qRT-PCR data are displayed in one graph. Experiments were completed at least three times with BMDM from three separate mice.
About this article
Cite this article
Menard, K.L., Bu, L. & Denkers, E.Y. Transcriptomics analysis of Toxoplasma gondii-infected mouse macrophages reveals coding and noncoding signatures in the presence and absence of MyD88. BMC Genomics 22, 130 (2021). https://doi.org/10.1186/s12864-021-07437-0
- Toxoplasma gondii
- Noncoding RNA