Temporal transcriptomic changes in long non-coding RNAs and messenger RNAs involved in the host immune and metabolic response during Toxoplasma gondii lytic cycle

Background Long non-coding RNAs (lncRNAs) are important regulators of various biological and pathological processes, in particular the inflammatory response by modulating the transcriptional control of inflammatory genes. However, the role of lncRNAs in regulating the immune and inflammatory responses during infection with the protozoan parasite Toxoplasma gondii remains largely unknown. Methods We performed a longitudinal RNA sequencing analysis of human foreskin fibroblast (HFF) cells infected by T. gondii to identify differentially expressed long non-coding RNAs (lncRNAs) and messenger RNAs (mRNAs), and dysregulated pathways over the course of T. gondii lytic cycle. The transcriptome data were validated by qRT-PCR. Results RNA sequencing revealed significant transcriptional changes in the infected HFFs. A total of 697, 1234, 1499, 873, 1466, 561, 676 and 716 differentially expressed lncRNAs (DElncRNAs), and 636, 1266, 1843, 2303, 3022, 1757, 3088 and 2531 differentially expressed mRNAs (DEmRNAs) were identified at 1.5, 3, 6, 9, 12, 24, 36 and 48 h post-infection, respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DElncRNAs and DEmRNAs revealed that T. gondii infection altered the expression of genes involved in the regulation of host immune response (e.g., cytokine–cytokine receptor interaction), receptor signaling (e.g., NOD-like receptor signaling pathway), disease (e.g., Alzheimer's disease), and metabolism (e.g., fatty acid degradation). Conclusions These results provide novel information for further research on the role of lncRNAs in immune regulation of T. gondii infection. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-021-05140-3.

there is no effective human vaccine to prevent infection, and drugs used to treat infected individuals have limitations [1,2].
Toxoplasma gondii divides by endodyogeny inside an intracellular parasitophorous vacuole (PV) [5]. When the number of tachyzoites within the PV reaches a certain level, the tachyzoites egress from infected cells and infect new host cells in a process known as the lytic cycle [6]. Following host cell invasion, the proliferating tachyzoites of T. gondii establish an acute infection. In response to the pressure conferred by the host immune response, tachyzoites differentiate to slowly replicate bradyzoites and establish a lifelong latent infection, which can be reactivated in immunodeficient hosts [7]. The interaction of T. gondii with the host cells occurs at the transcriptional and post-transcriptional levels, which is fundamental for establishing a lytic or latent infection, but the molecular mechanisms that mediate these interactions remain poorly understood.
Several approaches, including transcriptomics [8][9][10], proteomics [11][12][13], and metabolomics [14], have been used to reveal the mechanisms that mediate the interaction between T. gondii and the host cell. Long non-coding RNAs (lncRNAs) have recently received increased attention because they play vital roles in cell development [15], chromatin modification [16], and immune regulation [17]. Changes in lncRNA expression have been associated with developmental defects [18], tumorigenesis [19], and autoimmune diseases [20]. There are a few studies on the expression of lncR-NAs and their role in infectious diseases, such as those caused by viruses [21,22]. Interestingly, deregulation of the expression of 996 lncRNAs and 109 messenger RNAs (mRNAs) has been demonstrated in human foreskin fibroblast (HFF) cells infected by T. gondii by using microarray analysis [23].
In the present study, we used RNA sequencing to obtain greater insight into the temporal changes in the expression of lncRNAs and mRNAs during the lytic cycle of T. gondii infection in HFF cells.

Parasites and cell culture
The HFF cells, purchased from the American Type Culture Collection (ATCC), were cultured in Dulbecco's modified Eagle medium (DMEM) supplemented with 15% (vol/vol) fetal bovine serum (FBS), 100 U/ml penicillin, and 100 µg/ml streptomycin. The cell culture was maintained at 37 °C and 5% CO 2 . The RH (Type I) strain used in the present study was maintained in HFF cells as described previously [24].

Sample collection, RNA extraction, and RNA sequencing (RNA-seq)
The HFF cells were infected with 10 × 10 6 freshly egressed tachyzoites using a multiplicity of infection (MOI) of 1 (1 tachyzoite to 1 HFF cell). After 1.5 h, the culture medium was removed and replaced with fresh DMEM. The cultured cells were further incubated at 37 °C with 5% CO 2 . Then, infected cell samples were collected at 1.5, 3, 6, 9, 12, 24, 36 and 48 h post-infection (hpi) with T. gondii. Mock-infected cells were collected at 0 h. The experiment was carried out three times. The collected cell samples were stored at −80 °C until use for RNA extraction. Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions.
The quantity and quality of the extracted RNA were determined by using the NanoDrop spectrophotometer and Agilent 2100 Bioanalyzer (Thermo Fisher Scientific, MA, USA), respectively. Samples with an RNA integrity number (RIN) ≥ 8.0 were used for RNA sequencing. The construction of sequencing libraries and sequencing were performed at Beijing Genomics Institute in Shenzhen using the BGISE500 platform.

Data processing and differential expression analysis
The raw sequencing data were analyzed using SOAPnuke (v1.5.2) [25] to obtain high-quality reads. HISAT2 was then used to map the clean reads to the reference genome [26], followed by the application of Bowtie 2 (V2.2.5) [27] to align the clean reads to the gene set. A novel, coding and non-coding transcript database built by BGI (Beijing Genomic Institute in Shenzhen). RSEM (v1.2.12) was used to calculate gene expression levels [28], and differential expression analysis was performed using DESeq2 (v1.4.5) with Q value ≤ 0.05 [29].

Function prediction of the differentially expressed lnRNAs
Gene Ontology (GO) (http:// www. geneo ntolo gy. org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http:// www. genome. jp/ kegg/) enrichment analysis of the annotated host genes corresponding to differentially expressed lnRNAs (DElncRNAs) and DEm-RNAs were performed by Phyper (https:// en. wikip edia. org/ wike/ Hyper geome tric_ distr ibuti on) based on the hypergeometric test. The significance of the enriched GO terms and KEGG pathways was corrected by Q value with a threshold (Q value ≤ 0.05) using the Bonferroni method [30].

Quantitative real-time PCR analysis
Validation of the differential expression of lnRNAs and mRNAs identified by RNA-seq was performed by using quantitative real-time PCR (qRT-PCR) analysis of eight randomly selected RNA samples. qRT-PCR was performed using the SYBR assay according to the manufacturer's instructions in a total of 20 µl reaction volume, containing 1 µl of each primer, 2 µl cDNA, 6 µl nucleasefree H 2 O, and 10 µl master mix (Takara, Japan). Three replicates were performed for each gene, and the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used as an internal control. The level of gene expression was performed by using the 2 −∆∆CT method. The primers used in the study (Table 1) were designed according to the sequences available in the NCBI database and synthesized at Sangon Biotech (Shanghai, China).

Statistical analysis
All statistical analyses were performed using GraphPad Prism 7 software (GraphPad Software, Inc., San Diego, CA, USA) and a value of P < 0.05 was considered statistically significant.

The gene targets of DElncRNAs
To identify the biological functions of DElncRNAs, GO enrichment analysis was performed using the gene target of DElncRNAs. The top 20 most-enriched GO terms in the biological processes are shown in Fig. 3. The biological processes of DElncRNAs at different groups were mainly involved in DNA recombination (GO:0006310), reverse transcription involved in RNA-mediated transposition (GO:0032199), transcription, DNA-templated (GO:0006351), and nucleic acid phosphodiester bond hydrolysis (GO:0090305) (Fig. 3). The top 20 mostenriched KEGG pathways of DElncRNAs are shown in Table 1 The qRT-PCR primers used in the study a The GAPDH gene was used for mRNA and lncRNA normalization  Fig. S4 and Additional file 9: Fig. S5).

Functional analysis of the DEmRNAs
The top 20 most-enriched biological processes are shown in Fig. 5. The biological processes of DEmRNAs at different time points were mainly involved in immune-or inflammation-related GO terms, such as the immune response (GO:0002376), inflammatory response (GO:0006954), defense response to virus (GO:0051607) and cytokine-mediated signaling pathway (GO:0019221). The top 20 most-enriched pathways of DEmRNAs are shown in Fig. 6. Most of the DEmRNAs were involved in immune-related pathways, such as cytokine-cytokine   CCL2 and HOXA107). The x-axis shows the names of the selected RNAs and the time points after infection. The y-axis shows the relative expression levels. The GADPH gene was used for mRNA and lncRNA normalization. qRT-PCR data are mean of three replicates ± standard deviation. Results obtained by RNA-seq analysis and qRT-PCR produced similar gene expression patterns receptor interaction (e.g., IL6, CXCL10 and CSF2), IL-17 signaling pathway (e.g., CCL2, NFKB1 and FOSL1), TNF signaling pathway (e.g., MAP3K8, IRF1 and ICAM1), and NOD-like receptor signaling pathway (e.g., GBP1, GBP2 and IFNAR2). The upregulated mRNAs were mainly involved in immune-related pathways (Additional file 10: Fig. S6), and the downregulated mRNAs were mainly involved in metabolism-related pathways and neurodegenerative diseases-related pathways, such as Alzheimer's disease and Huntington's disease (Additional file 11: Fig.  S7).

Functional analysis of the gene targets of the co-expressed DElncRNAs
GO analysis and KEGG pathway enrichment were performed at different time points to identify the function of host genes corresponding to the DElncRNAs. The top 20 most-enriched biological processes are shown in Fig. 7a. They were mainly involved in reverse transcription linked to RNA-mediated transposition (e.g., EGOT, LYRM4-AS1 and SLCO4A1-AS1), nucleic acid phosphodiester bond hydrolysis (e.g., LYRM4-AS1, MSC-AS1 and GNG12-AS1), and DNA recombination (e.g., PLCE1-AS1, ZNF667-AS1 and LYRM4-AS1). The top 20 mostenriched pathways of DElncRNAs are shown in Fig. 7b. Most of DElncRNAs were involved in autophagy and mitophagy (e.g., LINC00173, LUCAT1 and ITGB1-DT).

Functional analysis of the co-expressed DEmRNAs
The top 20 most-enriched biological processes are shown in Fig. 8a. The biological processes of DEm-RNAs were mainly involved in immune-related GO terms, such as the immune system process (e.g., IFI30, TRIM13, and DHX36), defense response to virus (e.g., TRIM22, IFIT5, and IFI6), and inflammatory response (e.g., NLRP3, IFI16, and IL-6). The top 20 mostenriched pathways of DEmRNAs are shown in Fig. 8b. Most DEmRNAs were involved in immune-related pathways such as cytokine-cytokine receptor interaction (e.g., CSF1, CXCL1, and IL1A), TNF signaling

Competing endogenous RNA (ceRNA) networks of lncRNAs-mRNAs
To reveal the correlations between DElncRNAs and DEmRNAs, the ceRNA network was constructed using Cytoscape v3.6.0 (2015), based on the potential target relationship between DElncRNAs and DEmRNAs. One lncRNA was associated with one or more mRNAs (Fig. 9).
As shown in Fig. 9, lncRNAs can interact with innate immune-related mRNAs and apoptotic-related mRNAs; for example, HOXA10-AS can interact with MAVS, STAT1, IL10 and BCL2L1. However, the underlying regulatory mechanisms of these DEmRNAs and DElncRNAs need to be elucidated in further studies.

Discussion
The progressive proliferation of T. gondii tachyzoites in the host tissue causes acute infection, disrupting physiological homeostasis and resulting in clinical manifestations such as fever and encephalitis [31]. Understanding how the expression of cellular lncRNAs and mRNAs is altered during T. gondii infection can provide new insight into the mechanisms underlying the early stages of the interaction between T. gondii and its host. In this study, a total of 7722 DElncRNAs and 16,446 DEmRNAs were identified, and DElncRNAs were mainly involved in immuno-inflammation and metabolism. These results suggest that lnRNAs play a significant role in the pathogenesis of T. gondii infection.
Recent studies have revealed that pathogen invasion can alter the host lnRNA expression, which seems to help pathogens evade the host immune response [35,36] or enhance the host innate immune response [37]. In the present study, GO and KEGG analysis revealed that many of the DElncRNAs play a role in immune regulation (Fig. 4). Toxoplasma gondii can inactivate human plasmacytoid dendritic cells by functional mimicry of IL-10 [38], and lncRNA TMC3-AS1 can negatively regulate the expression of IL-10 [39]. In the present study, lncRNA TMC3-AS1 was downregulated. The deletion of lncRNA LUCAT1 in myeloid cells increases the expression of type I interferon-stimulated genes in response to LPS, and overexpression of LUCAT1 reduces inducible ISG response via their interaction with STAT1 in the nucleus [40]. In the present study, the expression of LUCAT1 was upregulated. Previous studies showed that MAVS plays a role in inducing an innate immune response and many viruses employ various strategies to suppress its activity [41][42][43]. In the present study, T. gondii downregulated the expression of MAVS in infected cells. These findings indicate that T. gondii infection can manipulate the host's immune response by altering the expression of host lncRNAs (e.g., LUCAT1 and MAVS); however, the exact     Latent T. gondii infection has been associated with various psychiatric disorders in humans [44,45]. However, the mechanism underlying this connection remains a subject of further investigation [12]. In the present study, T. gondii infection increased the expression of genes associated with Alzheimer's disease, Parkinson's disease, and Toxoplasmosis, and genes that play a role in the neurotransmitter release cycle and transmission across chemical synapses. Alterations in solute carrier family genes, such as SLC2A6, SLC1A5, and SLC22A4, have been implicated in the impairment of the synaptic transmission, which seems relevant to the pathology of schizophrenia [46]. Other studies have revealed that altered expression of DRD1, which was downregulated in the present study, may contribute to the pathophysiology of schizophrenia and affective disorders [47].
Recent studies have enriched our understanding of the metabolism of both T. gondii and its host cells [14,48,49]. In the present study, we have shown that T. gondii infection changes the host cell's metabolism by altering the expression of mRNAs and lncRNAs involved in carbon metabolism, metabolic pathways, fatty acid degradation, biosynthesis of secondary metabolites, and pyruvate metabolism (Figs. 4, 6). The tricarboxylic acid (TCA) cycle is essential for T. gondii growth [48]. In the absence of glucose, T. gondii can utilize the gluconeogenic enzyme fructose bisphosphatase 2 to provide carbon for gluconeogenesis [49,50]. In the present study, the downregulated lnRNAs were mainly involved in metabolismrelated pathways (Additional file 3: Fig. S3), indicating that lncRNAs play a role in regulating host metabolism during T. gondii infection.

Conclusions
Our data revealed significant changes in lncRNA and mRNA expression in HFF cells following T. gondii infection. The identified DElncRNAs and DEmRNAs were mainly involved in metabolism, signal transduction, and immune responses. Further investigation of these DEl-ncRNAs and DEmRNAs is warranted to reveal the exact role of these transcripts in the pathophysiology of T. gondii infection.