Hundreds of putative lncRNAs are differentially regulated during infection with T. gondii

We infected mouse bone marrow-derived macrophages (BMDM) with high (Type I, RH) and low (Type II, PTG) virulence T. gondii at a 4:1 ratio of parasites to cells. This MOI resulted in approximately 80% infection for both strains (Supplementary Fig. S1). Samples of RNA were prepared 6 hr after infection, a time prior to the first parasite mitotic division when on average there was a single parasite per infected cell (Supplementary Fig. S1). Using a commercial microarray (Arraystar), we examined the differential gene expression of 35,923 putative lncRNAs after infection. The array was constructed using public transcriptome databases (Refseq, UCSC knowngenes, Ensembl, etc.), as well as relevant publications23,24. Differentially expressed lncRNAs and mRNAs were identified through p-value (<0.05) and fold change filtering (>2 fold up or down). The complete dataset is provided in Supplementary Data S1. Hierarchical clustering of the differentially regulated lncRNAs in three replicate experiments indicated that both RH and PTG infections regulated lncRNA expression and that the lncRNA landscape differed between the two infections (Fig. 1a). We also used volcano plot analysis to visualize the differentially expressed lncRNAs. This technique reveals significant changes in large sets of replicate data by plotting fold changes versus p-values. We identified substantial populations of lncRNAs that were up- or down-regulated in a statistically significant manner by RH and PTG infection (Fig. 1b). As shown in Fig. 1c, 770 lncRNAs were upregulated during RH infection, and a similar number were down-regulated. Infection with PTG had less extensive effects on lncRNA expression, with a total of 277 species upregulated and 251 down-regulated. 133 lncRNAs were upregulated by PTG relative to RH, while 174 lncRNAs were down-regulated. A total of 282 lncRNAs were either up- or down-regulated during infection with each strain type as depicated by Venn diagrams (Fig. 1d and Supplementary Data S2). However, a substantial number of lncRNAs were regulated in a parasite-strain-specific manner. Taking together all three comparisons (RH vs. uninfected, PTG vs. uninfected, and PTG vs. RH), 1,902 unique lncRNAs were identified in the microarray as significantly differentially regulated. We examined how these lncRNAs were distributed according to protein coding genes (Fig. 1e). The majority of the lncRNAs were either intergenic (>1 kb from the promoter of a protein coding gene) or exon-sense overlapping (overlapping an exon of a protein-coding gene in the sense direction).

Figure 1 Hundreds of putative lncRNAs are differentially regulated during infection with Toxoplasma, as determined using a microarray for mouse lncRNAs. Mouse BMDM were infected with either the highly-virulent RH or the less-virulent PTG strain, and 6 hr later, RNA was isolated for microarray analysis. (a) Hierarchical clustering of differentially expressed lncRNAs for uninfected vs. RH-infected, uninfected vs. PTG-infected, and PTG vs. RH-infected. Values in the color scale are normalized intensities. Red bands indicate high relative expression, and green bands indicate low relative expression. (b) Volcano plot filtering to visualize fold regulation and statistical significance in lncRNA populations. Statistically significant (p-value < 0.05) up (red) or down (green) regulated expression changes (>2-fold) are shown in pairwise comparisons of uninfected, RH-infected, and PTG-infected samples. (c) Total number of lncRNA species up- or down-regulated by infection. (d) Venn diagrams of up- and down-regulated lncRNA populations reveal shared and unique expression patterns. (e) Classification of lncRNA differentially regulated by Toxoplasma infection. Experiments were performed in triplicate with BMDM from three separate mice. Full size image

Not only were many lncRNAs differentially regulated, but also the magnitude of the observed fold changes was large (Table 1). For all three strain comparisons, the top 20 hits had fold changes larger than ±7, and 31 lncRNAs had fold changes larger than ±20. Notably, infection with high-virulence RH had greater impact on the number and magnitude of lncRNA responses than low-virulence PTG. For example, the average magnitiude of upregulated responses shown in Table 1 was 54.8 and 18.7 for RH and PTG, respectively. These effects could not be attributed to differences in replication rate (and therefore antigen load per cell) because the analyses were carried out 6 hr after infection, prior to the first parasite cell division. Also, the difference in the number of lncRNAs regulated between RH and PTG was not due to invasion efficiency, as percent infection was approximately 80% for each strain.

Table 1 lncRNAs with the 20 largest fold changes for each comparison: RH vs. uninfected, PTG vs. uninfected, and PTG vs. RH. Full size table

We also note that 17 putative lncRNAs that were regulated by infection mapped to ultraconserved regions of the genome25, suggesting that they may function across species (Table 2). This result is potentially significant because lncRNA are generally regarded as poorly conserved through evolution. Using qRT-PCR, we validated expression for one ultraconserved region uc.70 (Supplementary Fig. S2). We note that these ultraconserved regions were identified through genomic sequence comparisons and are generally not based on transcriptomic data. Therefore, the transcripts for these ultraconserved regions would need to be mapped out before further study.

Table 2 Differentially regulated lncRNAs that map to ultraconserved genomic regions. Full size table

Hundreds of mRNAs are differentially regulated during infection with Toxoplasma

In addition to lncRNAs, the microarray contains probes for 24,881 mRNAs, enabling the measurement of changes in expression of protein-encoding transcripts in response to infection. Consistent with results from Fig. 1, hierarchical clustering (Fig. 2a) and volcano plot filtering (Fig. 2b) revealed parasite-strain-specific control of mRNA expression. 1,583 mRNAs were differentially regulated during infection with the high-virulence RH strain compared to an uninfected control (Fig. 2c; the complete data set is shown in Supplementary Data S3). 582 mRNAs were differentially expressed during infection with the less-virulent PTG strain compared to an uninfected control. 358 mRNAs were differentially regulated between RH and PTG. Venn diagrams (Fig. 2d) show that there is overlap (342 total) in the up- or down-regulated mRNAs by RH and PTG infection. A compiled list of mRNAs regulated during both RH and PTG infection is included in Supplementary Data S4. In addition, and particularly for RH, there was substantial parasite-strain-specific regulation of mRNA expression. Altogether, 1,927 unique mRNAs were identified in the microarray as significantly differentially regulated across all three comparisons (RH vs. uninfected, PTG vs. uninfected, and PTG vs. RH). Table 3 shows the top 20 differentially regulated mRNAs for pairwise comparisons of RH vs. uninfected, PTG vs. uninfected, and PTG vs. RH.

Figure 2 Differential expression patterns of BMDM mRNA during infection with Toxoplasma. Changes in mRNA expression were assessed 6 hr after infection of BMDM. (a) Hierarchical clustering of pairwise combinations of uninfected, RH-infected, and PTG-infected samples. Values in the color scale are normalized intensities. Red bands indicate high relative expression, and green bands indicate low relative expression. (b) Volcano plot filtering to reveal statistically significantly (p-value < 0.05) up (red) or down (green) regulated (>2-fold) expression changes in pairwise comparisons of uninfected, RH-infected, and PTG-infected samples. (c) Total number of mRNAs up- or down-regulated by infection. (d) Venn diagrams of up- and down-regulated mRNA expression reveal unique and shared mRNA species regulated by RH and PTG infection. Experiments were performed in triplicate with BMDM from three separate mice. Full size image

Table 3 Magnitude of the 20 largest mRNA fold changes during RH and PTG infection in each indicated comparison: RH vs. uninfected, PTG vs. uninfected, and PTG vs. RH. Full size table

Gene ontology (GO) analysis identified many immune-related biological processes as differentially regulated in the microarray, including regulation of T-helper cell differentiation, interleukin 2 production, macrophage differentiation, JAK-STAT cascade, and cytokine secretion (Supplementary Fig. S3). Several of the biological pathways identified were regulated in a parasite strain-specific manner. For example, infection with RH was associated with regulation of macrophage differentiation and negative regulation of the JAK-STAT cascade (Supplementary Fig. S3e and f). KEGG Pathway analysis identified many immune-related and infection pathways as differentially regulated in the microarray, including cytokine-cytokine receptor interaction, TGF-beta signaling pathway, Toll-like receptor signaling pathway, RIG-I-like receptor signaling pathway, JAK-STAT signaling pathway, chemokine signaling pathway, Wnt signaling pathway, RapI signaling pathway, measles, malaria, HTLV-1 infection, Hepatitis B, Hepatitis C, and tuberculosis (Supplementary Fig. S4). KEGG analysis also revealed parasite strain-specific pathway enrichment. For example, RH infection was associated with selective enrichment of pathways involved in cytokine-receptor interaction, the JAK-STAT pathway and the p53 signaling pathway (Supplementary Fig. S4e and f).

Four co-regulated lncRNAs and mRNAs were validated by qRT-PCR

Many known lncRNAs regulate adjacent or overlapping protein-coding genes in a cis-regulatory manner26,27. For this reason, we examined the 282 lncRNAs that were co-regulated with a nearby protein-coding gene (Supplementary Data S5). Of these 282, we identified ~60 lncRNAs that were co-regulated with infection- or immune-related protein-coding genes (Supplementary Data S5), suggesting a potential regulatory role for these lncRNAs in the immune response. From this list, we chose 4 lncRNAs to validate by qRT-PCR: Socs2-lnc, Csf1-lnc, Il1rn-lnc, and Ifi44-lnc. Socs2-lnc (uc007gwf.2) is associated with the protein-coding gene Socs2, a suppressor of cytokine synthesis possibly involved in the anti-inflammatory response during T. gondii infection28. Csf1-lnc (ENSMUST00000155557) is associated with the protein-coding gene csf1, a cytokine that drives stem cell differentiation into the macrophage lineage29. Ifi44-lnc (ENSMUST00000133888) is associated with Ifi44, an interferon-alfa inducible protein connected with viral infection30. Il1rn-lnc (ENSMUST00000143423) is associated with Il1rn, an antagonist cytokine that inhibits activities of interleukin 1 and modulates the interleukin 1 immune and inflammatory responses31.

We designed primers to bind a unique region of each lncRNA and co-regulated mRNA. qRT-PCR largely confirmed that lncRNA transcripts were up- or down-regulated (≥or ≤2-fold) during infection with either the RH or PTG in a manner similar to the microarray (Fig. 3a and b). The microarray was considered valid if the directionality (up, down, or not regulated) of transcription changes compared to uninfected samples was consistent between the qRT-PCR data and the microarray data. Likewise, mRNA transcripts were similarly regulated in both the microarray and qRT-PCR (Fig. 3c and d). Nevertheless, during infection with the PTG strain, we identified two loci, Ifi44-lnc and Socs2, that were not expressed in the microarray but did appear to be regulated in the qRT-PCR.

Figure 3 Validation of lncRNA and mRNA microarray data by qRT-PCR. RNA from mouse BMDM were collected 6 hr after infection with either RH or PTG strains of T. gondii, and qRT-PCR was subsequently performed. Fold changes represent the comparison of infected samples to uninfected samples. Regulation of representative lncRNAs as determined by qPCR (a) and microarray (b) analysis. Regulation of immune-related mRNA is shown by qPCR (c) and microarray (d) analysis. Microarray fold change values differ slightly from those listed in Supplementary Spreadsheets 1 and 2 (which are the geometric means) because arithmetic means were calculated from the microarray data to directly compare to qRT-PCR data. Experiments were completed a minimum of three times with BMDM from three separate mice and were obtained independently of experiments used for microarray analysis. Full size image

Toxoplasma rhoptry kinase ROP16 manipulates expression of host lncRNA

During intracellular infection, T. gondii secretes a subset of effector molecules whose activities modify host responses to infection32. One example is the rhoptry protein ROP16, a kinase that directly activates host cell STAT3 and STAT617,21,33,34. While Type I ROP16 possesses strong STAT kinase activity, Type II ROP16 is catalytically inactive. Previous studies in BMDM indicate that 538 protein-coding genes are regulated by Type I (RH) ROP1618.

We hypothesized that the activity of ROP16 extends to control of lncRNAs, particularly those strongly up-regulated by the Type I RH strain relative to the Type II PTG. By using engineered ROP16 deletion and complementation mutants created on the RH background17, we examined expression of three lncRNAs that showed evidence of Type I strain-specific up-regulation: Socs2-lnc, Csf1-lnc, and Il1rn-lnc. Expression of Csf1-lnc and Socs2-lnc was significantly reduced during infection with the ROP16 deletion mutant (RH∆16) compared to RH (Fig. 4a and b). Complementation of RH∆16 with a functional copy of ROP16 (RH∆16:1) restored the up-regulation phenotype. In contrast, Il1rn-lnc was not controlled by ROP16, and the expression of Il1rn-lnc RH∆16 in Fig. 4c is similar to that of RH and RH∆16:1. As a positive control for a gene controlled by ROP16, we examined the expression of an mRNA (Il4i1) that we suspected was also controlled by ROP16 based upon the strain-specific expression pattern. Il4i1 protein-coding mRNA expression was also significantly lower in the RH∆16 strain (Fig. 4d). Our data demonstrate that the secreted parasite effector molecule ROP16 regulates several putative host lncRNAs.

Figure 4 lncRNAs are targeted for manipulation by T. gondii rhoptry kinase ROP16. BMDM were infected with RH, a ROP16 deletion mutant on the RH background (RHΔ16), a ROP16 complementation strain (RHΔ16:1), and PTG. Then, RNA was isolated 6 hr later for qRT-PCR analysis. Expression patterns of Socs2-lnc (a), Csf1-lnc (b), Il1rn-lnc (c) and Il4i1 mRNA (d) were determined. Fold changes represent the comparison of infected samples to uninfected samples. p < 0.05 was considered significant when comparing RH and RH∆16. The experiments were completed a minimum of three independent times. Full size image

Additionally, to isolate effects of live infection from other potential sources of gene regulation (such as the presence of dead parasites or soluble parasite protein), we examined expression of Socs2-lnc after addition of live, dead, and soluble tachyzoite antigen (STAg). Addition of dead parasites and STAg did not result in increased expression of Socs2-lnc, suggesting that Socs2-lnc is induced by live T. gondii infection (Supplementary Fig. S5).