Confirmation of Em infection in BALB/c mice and identification of DEMs and DELs
Em infection was confirmed in all challenged mice by observing gross pathological lesions. The infected livers were enlarged, and the surface of the liver lobes was characterized by translucent or whitish alveolar vesicles.
To explore the mRNA and lncRNA expression profiles of mouse liver in response to Em during the whole infection period, we selected eight time points to estimate the transcript expression levels between CA and Con groups by microarray analyses. At each time point, CA and Con groups were both tested using three biological replicates. The workflow for the transcript analysis is shown in Fig. 1a. The mRNAs and lncRNAs with an expression change of more than two-fold (FC ≥ 2) and p < 0.05 were selected as differentially expressed candidates. As shown in Table 1, the numbers of DEMs and DELs were detected for the corresponding Con at the eight time points. Notably, for DEMs, the proportion of downregulated genes was from 51.3% to 68.2% at seven time points except for 30 dpi (downregulated genes were 25.1%). For DELs, there were more downregulated genes than upregulated genes at all the eight time points, with a proportion of 58.4% to 83.7%. In general, the changes of DELs and DEMs displayed a similar trend. Furthermore, the volcano plots (See Additional Fig. S1a, b) and the hierarchical clustering plots (See Additional Fig. S1c, d) showed that the expression patterns of mRNAs and lncRNAs between the CA and Con groups at the eight time points were significantly different. These results suggested that Em is a strong trigger that can alter the expression of numerous mRNAs and lncRNAs at different stages during the infection course. More DEMs and DELs identified at 4, 15, 60, 90, 150 dpi indicates that these time points are representative of transcriptome change times.
To investigate the biotype of identified DELs, according to the position of lncRNA in the genome relative to protein-encoding gene, DELs at each autopsy time point were divided into four categories (http://www.noncode.org/). As shown in Fig. 1b, more than 78% lncRNA were lincRNAs, antisense lncRNA and sense no exonic lncRNA accounted for the nearly equal proportion (7–9%), while exonic lncRNAs were fewest (~3%). The data demonstrated that biotypes of the DELs and their corresponding ratios at different time points during the Em infection course were almost the same.
Validation of candidate DEMs and DELs using qRT-PCR
To validate the DELs/DEMs identified in mouse livers at each autopsy time point by microarray, twenty DEMs and DELs (10 each) randomly selected were verified by qRT-PCR. As shown in Fig. 2, the qRT-PCR results of 10 mRNAs (six upregulated and four downregulated) and 10 lncRNAs (six upregulated and four downregulated) confirmed the microarray data, showing similar trends in the up- or downregulated mRNAs (Fig. 2a) and lncRNAs (Fig. 2b).
Functional annotation of thedifferentially expressed mRNAs and overall gene expression changes in Em-infected mouse liver
To better understand the biological roles of differentially expressed mRNAs at each autopsy time point, the up- and downregulated DEMs were estimated by KEGG pathway analysis. In Fig. 3a, it was found that two innate immunity pathways, Toll-like (TLRs) and RIG-I-like receptor (RLRs) signaling pathway, were included in the significantly enriched KEGG terms (top 30) of the DEMs at 30 dpi. The results showed that IRAK4, IKKβ, MKK, JNK, and IFN-α were upregulated in TLRs signaling pathway (Additional Fig. S2a), and NLRX1, IKKβ, JNK, and IFN-α were upregulated in RLRs signaling pathway (Additional Fig. S2b). Additionally, the KEGG analysis of downregulated DEMs showed that the ubiquitin-mediated proteolysis pathway was included in the significantly enriched KEGG terms at seven time points (2, 4, 8, 15, 30, 60, 150 dpi). Notably, there were more dysregulated genes in this pathway at 15 dpi (Additional Fig. S2c). This result indicated that Em infection had a broad and long-lasting impact on this pathway of the host.
Additionally, to search the genes with continued aberrant expression during the infection course, 68 DELs (downregulated) and 31 DEMs (21 downregulated and 10 upregulated) were found overlapped at all time points by Venn analysis (Additional Table S1a, b). Then, GO analysis was performed to understand the potential functions of these shared DEMs during the whole Em infection period. In Fig. 3b, the results showed that these DEMs were significantly enriched in GO terms related to BP, such as “antigen processing and presentation of peptide, polysaccharide or exogenous peptide antigen via MHC class II,” “response to interferon-gamma” and “negative regulation of T cell proliferation”; most shared DEMs were associated with GO terms related to CC, including “MHC class II protein complex,” “brush border” and “lysosomal membrane”; these shared DEMs were enriched in MF associated with transportation, including “peptide antigen binding,” “symporter activity” and “transporter activity.”
Next, we conducted 31 shared DEMs by the KEGG database to determine the signal transduction and/or metabolic pathways of these dysregulated mRNAs to understand the host–parasite interplay further. As shown in Fig. 3c, these overlapped DEMs throughout the infection course were mainly involved in “antigen processing and presentation,” “intestinal immune network for IgA production,” “Th1 and Th2 cell differentiation,” and “Th17 cell differentiation,” suggesting that Em infection had an enduring influence on these signaling pathways. They were continuously active throughout the infectious course.
Construction of co-expression network between overlapped DELs and DEMs
To demonstrate the interactive relationship between the overlapped DELs and DEMs, 31 shared DEMs and 68 shared DELs containing 2,108 relationships were estimated by calculating PCC values. According to the cutoff criteria (p < 0.05 and |PCC| > 0.7), Fig. 4 showed the lncRNA-mRNA co-expression network of the top five lncRNAs that had the most co-expression relationships with mRNAs. From the network, it was observed that NONMMUT099399.1 had the most relationships with DEM, showing that more than one mRNA was predicted to be regulated by one lncRNA. In turn, one mRNA could correspond to several lncRNAs, for example, mRNA Cers6 may be regulated by three lncRNAs (NONMMUT037651.2, NONMMUT099399.1, and NONMMUT125804.1). Furthermore, these five lncRNAs were co-expressed with 26 DEMs, which accounted for 84% of the total shared DEMs, indicating that these lncRNAs are the candidates with the highest potential for being involved in the signaling pathways enriched by the 31 shared DEMs as displayed by KEGG analysis.
Cis and trans regulating target prediction of overlapped lncRNAs and construction oflncRNA-TF-mRNA network
Generally, lncRNAs perform functions by interacting with their targets, so we predicted the cis- and trans- target gene potential of these lncRNAs. We searched for protein-coding genes (among the 31 shared DEMs) 10-kb up- and downstream of the 68 shared DELs and found no lncRNAs that were transcribed near (< 10 kb) the 31 DEMs neighbors. Then, the trans regulatory functions of lncRNAs were predicted by evaluating the relationship between the TFs and lncRNAs based on the cumulative hypergeometric test. We found that among the 68 shared lncRNAs, 41 lncRNAs corresponded to 16 TFs (EGR1, EGR2, EGR3, EGR4, STAT1, STAT5A, SMAD3, MYC, HINFP, ONECUT2, P2RX5, NR5A2, TRP53, MAX, PITX2, HIC1), which could regulate the expression of lncRNAs. In Fig. 5a, among the 16 TFs, EGR2 had the most trans relationships with DELs. EGR1, EGR3, and EGR4 were also related to the shared DELs, indicating that the EGR family may play an important role on the host DELs during Em infection. Then, since each lncRNA-TF pair results from several gene enrichment, we selected the top 10 co-expressed mRNAs by p-value ranking to further build a lncRNA-TF-mRNA ternary network (Fig. 5b).
Validation of six up-regulated DEMs involved in Toll-like and RIG-I-like receptor signaling pathways by qRT-PCR and WB
To further investigate the influence of Em infection on host innate immunity, six upregulated DEMs (IRAK4, IKKβ, MKK, JNK, IFN-α, and NLRX1) enriched in TLRs and RLRs signaling pathway in KEGG analysis were validated by qRT-PCR. In Fig. 6a, the results showed that expression levels of all these genes increased in the liver of CA mice compared to the control at 15 and 30 dpi, consistent with the microarray data (Additional Fig. S3). Then, the protein expression of IKKβ, JNK, MKK, IRAK4, and NLRX1 was determined by WB. Compared with the Con groups, the protein expression of IKKβ, MKK, IRAK4, and NLRX1 were unaltered at 15 and 30 dpi except for c-Jun amino-terminal kinase (JNK), which decreased at 30 dpi. And importantly, the protein expression of phospho-nuclear factor Kappa B (p-NF-κB) increased compared to the control at 15 and 30 dpi (Fig. 6b). Taken together, these results demonstrated that Em infection induced mRNA expression of some key components in TLRs and LRLs signaling pathway, but their protein expression levels were not altered, suggesting that the changes in the expression level of these genes during transcription and translation were inconsistent and may be due to the parasite effect on post-transcriptional regulation. In any case, downstream these two innate immune pathways, p-NF-κB, a master regulator of key inflammatory gene expression, was induced at 15 and 30 dpi.
Determination of cytokines and chemokines downstream TLRs and RLRs signaling pathway
To assess the impact of Em infection on the cytokines and chemokines downstream of TLRs and RLRs signaling pathway shown by KEGG database analyses (See Additional Fig. S2a, b), the serum levels of IL-1β, IL-6, CCL5 (RANTES), TNF-a, IL-12, CCL3 (MIP1-a), IFN-a, IFN-β, CXCL10 (IP10) were determined at 2, 8, 15 and 30 dpi. The results showed that IL-1β, IL-6, TNF-α, and CXCL10 increased significantly at 30 dpi (p < 0.05). IFN-α, IFN-β, IL-12, CCL5, and CCL3 demonstrated an increasing trend but the differences were not statistically significant compared to the control (Fig. 7). Although not all the factors were significantly elevated, there was a clear underlying trend. Furthermore, among these, as the most important proinflammatory cytokines, levels of IL-1β and IL-6 in serum increased markedly at 30 dpi, which may be due toincreased levels of p-NF-κB at 15 and 30 dpi as evident in the WB results.