Expression Profiles of mRNA and lncRNA in HCT-8 Cells Infected With Cryptosporidium parvum IId Subtype.

Cryptosporidium parvum is one of the most important enteric protozoan pathogens, responsible for severe diarrhea in immunocompromised human and livestock. However, few effective agents were available for controlling this parasite. Accumulating evidences suggest that long non-coding RNA (lncRNA) played key roles in many diseases through regulating the gene expression. Here, the expression profiles of lncRNAs and mRNAs were analyzed in HCT-8 cells infected with C. parvum IId subtype using microarray assay. A total of 821 lncRNAs and 1,349 mRNAs were differentially expressed in infected cells at 24 h post infection (pi). Of them, all five types of lncRNAs were identified, including 22 sense, 280 antisense, 312 intergenic, 44 divergent, 33 intronic lncRNAs, and 130 lncRNAs that were not found the relationship with mRNAs' location. Additionally, real-time polymerase chain reactions of 10 lncRNAs and 10 mRNAs randomly selected were successfully confirmed the microarray results. The co-expression and target prediction analysis indicated that 27 mRNAs were cis-regulated by 29 lncRNAs and 109 were trans-regulated by 114 lncRNAs. These predicted targets were enriched in several pathways involved in the interaction between host and C. parvum, e.g., hedgehog signaling pathway, Wnt signaling pathway, and tight junction, suggesting that these differentially expressed lncRNAs would play important regulating roles during the infection of C. parvum IId subtype.

Introduction

Cryptosporidium, one of the most important enteric parasites to cause diarrhea in human and animals (Holland, 1990; Nguyen et al., 2013; Koinari et al., 2014; Delafosse et al., 2015; Deshpande et al., 2015; Nakamura and Meireles, 2015), has been recognized as the leading cause of chronic diarrhea in HIV patients and was the second contributor of moderate-to-severe diarrhea in children during the first 2 years of life (Checkley et al., 2015). After ingested, the sporozoites within Cryptosporidium oocysts were released and colonized into epithelial cells of the gastrointestinal tract, damaged the intestinal barriers to affect the host nutrition absorption, impaired immune response, and persistently retarded growth (Guerrant et al., 1999; Gookin et al., 2002; Mondal et al., 2009; Squire and Ryan, 2017).

Currently, 33 valid Cryptosporidium species have been confirmed and scattered around the world (Xiao, 2010; Fletcher et al., 2012; Hu et al., 2014; Ryan et al., 2014, 2015; Li et al., 2015; Holubová et al., 2016; Jezkova et al., 2016; Kváč et al., 2016; Feng and Xiao, 2017; Zahedi et al., 2017). Of them, Cryptosporidium parvum has been identified as the most common zoonotic species infecting humans and most animals (Wang et al., 2014). Intra-species genetic diversity of C. parvum was observed based on the nucleotide sequence of 60 kDa glycoprotein (GP60) gene, and at least 15 subtypes (IIa–IIi, IIk–IIp) were identified (Xiao, 2010; Insulander et al., 2013; Wang et al., 2014). In China, C. parvum has been detected in humans and 18 animal species from 19 provinces and one region (Qinling Mountain; Huang et al., 2014; Qi et al., 2015; Zhang X.X. et al., 2015; Cai et al., 2017). Of them, both subtypes IIa and IId were identified in China (Mi et al., 2013; Cui et al., 2014; Zhao et al., 2015). The subtype IId was widely detected in sheep, goats, and calves from China (Feng et al., 2013, 2017). Besides, the whole genome of two subtypes was different in sizes and contents. For example, the total assembly length of C. parvum IIaA15G2R1 was longer than IIdA19G1, and the gene gains/losses (e.g., SKSRa, MEDLE family of secreted proteins, Insulinase-like proteases) and single nucleotide variants (SNVs) were also found between the genomes of two subtypes (Widmer et al., 2012; Feng et al., 2017).

In the process of C. parvum invasion and parasitism, the host–parasite interaction occurred. The expression profiles of both mRNA and non-coding RNA (ncRNA) during the infection of C. parvum Iowa isolate (IIa subtype) have been investigated (Abrahamsen et al., 1996; Deng et al., 2004; Zhou et al., 2014). The differentially expressed mRNAs were predicted to be associated with the promoter enrichment of suppressive epigenetic marker, while ncRNAs may modulate epithelial immune responses and epithelial anti-microbial defense against Cryptosporidium infection (Abrahamsen et al., 1996; Zhou et al., 2009, 2014; Wang et al., 2017b). However, the gene expression in hosts infected with C. parvum subtype IId was not available. Additionally, a new RNA molecule, long non-coding RNA (lncRNA), was discovered and studied recently (Okazaki et al., 2002). Increasing evidences have certificated that lncRNAs could interact with mRNAs through cis- and trans-regulation to play key roles in tumorigenesis, tumor development (He et al., 2014; Lee et al., 2016), and infections of viruses (e.g., hepatitis C virus, influenza A virus, and severe acute respiratory syndrome coronavirus) and parasites (e.g., Leishmania, Plasmodium falciparum; Josset et al., 2014; Broadbent et al., 2015; Zhang H. et al., 2015; Pawar et al., 2017). To deeply understand the interaction between C. parvum and host, herein, we systematically investigated the expression profiles of mRNA and lncRNA in human cells infected with the C. parvum IIdA19G1 subtype.

Materials and Methods

C. parvum Isolates

The oocysts used in the present study were obtained from a pre-weaned dairy calve with diarrhea in China and molecularly identified as C. parvum IIdA19G1 subtype based on the sequence of the gp60 gene locus. This isolate was passaged by pre-weaned dairy calves in the laboratory with the specific pathogen-free condition. C. parvum oocysts were purified by using the Sheather’s sugar flotation technique and cesium chloride density gradient centrifugation, and stored in PBS with the penicillin-streptomycin (100 U/ml penicillin and 0.1 mg/ml streptomycin) and amphotericin B solutions (0.25 μg/ml).

In Vitro Infection Model of C. parvum

The human adenocarcinoma (HCT-8) cell lines were purchased from JENNIO Biological Technology (Guangzhou, China). 2 × 105 HCT-8 cells were seeded in each well of a fresh 24-well plate and cultured for 24 h (or with 80% confluence) at RMPI 1640 medium and supplemented with 10% fetal bovine serum (FBS) under 5% CO2 at 37°C. C. parvum oocysts were treated with 2% bleach for 20 min on ice and incubated into HCT-8 cells with the ratio of oocysts: cells = 2–10:1. The infection burden was measured using quantitative real-time polymerase chain reaction (qRT-PCR) targeting the small subunit ribosomal RNA (SSU rRNA) previously described (Zhao et al., 2018).

Sample Collection, RNA Extraction, and Microarray Analysis

The cell samples were collected from both experimental (C. parvum infection, O) and control (without parasites, C) groups at 24 h post infection (pi) of C. parvum oocysts. Three biological repeats were included in each group. Total RNA was extracted using TRIzol reagent (500 μl) and the chloroform-isopropyl alcohol method in accordance with the manufacturer’s instructions and stored at −80°C. The concentration and purity of total RNA samples were measured by the Smart Spec Plus spectrophotometer. The complementary DNA (cDNA) was generated using the PrimeScript TM RT reagent Kit with the gDNA Eraser (TaKaRa Shuzo Co., Ltd., Liaoning, China) following the manufacturer’s instructions for reverse transcription of the total RNAs (1 μg). Then, the samples were taken and sent to the company (CapitalBio Technology Corporation, Beijing, China) for microarray analysis. The expression profiles of mRNAs and lncRNAs were detected by using LncRNA+mRNA Human Gene Expression Microarray V4.0 (4 × 180 K). The Agilent Feature Extraction v10.7 was used to analyze and extract data, and these data were then normalized and analyzed by using Agilent GeneSpring GX software (13.1 revisions, 2015[1]).

GO and KEGG Enrichment Analysis

To investigate the biological functions of mRNAs and lncRNAs, the Gene Ontology (GO) analysis was executed. Three parts were involved in the GO terms, including biological process (BP), cellular component (CC), and molecular function (MF). The enrichment analysis of GO terms with P < 0.05 was considered significantly. Besides, Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to test the statistical enrichment of differentially expressed mRNAs and lncRNAs to predict the possible pathways involved. Both GO and KEGG were carried out by KOBAS 2.0 software (2016[2]).

Validation of Microarray Data by qRT-PCR

In order to validate the microarray results, five up- and five down-regulated genes of both mRNA and lncRNA were, respectively, selected for qRT-PCR validation, with GAPDH as an internal control, and the sequences of primers listed in Supplementary Table . The qRT-PCR was carried out using the SYBR assay in a 10 μl reaction volume, containing 0.2 μl Forward Primer, 0.2 μl Reverse Primer, 1 μl cDNA, 3.6 μl nuclease-free H2O, and 5 μl Master Mix (CWBIO, Beijing, China). The reaction program was initiated at 95°C for 10 min, then at 95°C for 10 s, 60°C for 60 s for a total 40 cycles. Three replicates were conducted for each gene, and the data were expressed as 2−ΔΔCt to value the expression of mRNA and lncRNA.

Statistical Analysis

All statistical analysis in this study was performed using the software GraphPad Prism 7.0 (2016[3]), with P < 0.05 considered as statistically significant.

Results

Identification of mRNA and LncRNA Differentially Expressed in C. parvum Infection Model

To explore the impact of Chinese prevalent C. parvum subtype (IIdA19G1 subtype) on human cells, HCT-8 cells were exposed to C. parvum IIdA19G1 subtype for 24 h and collected for LncRNA+mRNA Human Gene Expression Microarray analysis. The result revealed the differential expression profiles of mRNA and lncRNA in HCT-8 cells after C. parvum infection. Because of mild influence of C. parvum to host cell gene transcription (Deng et al., 2004; Zhou et al., 2009; Ming et al., 2017), the fold change ≥ 1.2 and P < 0.05 were used as the standard to map differentially expressed genes. All microarray data obtained from our study were deposited into GEO database with the number GSE111565. The expression patterns of mRNA (Supplementary Figure ) and lncRNA (Supplementary Figure ) between experimental and control groups were found to be significantly different by hierarchical clustering plot. A total of 1349 mRNAs (including 535 up- and 814 down-regulated; Figure and Supplementary Table ) and 821 lncRNAs (including 557 up- and 264 down-regulated; Figure and Supplementary Table ) were found to be differentially expressed by volcano plot and scatter plot filtering (Supplementary Figure ). These lncRNAs were grouped into five types reported previously, including 22 sense, 280 antisense, 312 intergenic, 44 divergent, and 33 intronic lncRNAs. Additionally, 130 lncRNAs that were not found the relationship with mRNAs’ location were also found to be differentially expressed.

Bioinformatics analysis of differentially expressed mRNAs and lncRNAs in HCT-8 cells infected with Cryptosporidium parvum IId subtype. (A) The volcano plot shows the distributions of mRNAs. (B) The volcano plot shows the distributions of lncRNAs. The significantly up- and down-regulated RNAs are presented as red or green dots, respectively, and the expression of RNAs not significantly differentially expressed is presented as black dots (fold change ≥ 1.2 and P < 0.05).

Validation of Differentially Expressed mRNAs and LncRNAs

To validate the microarray data, 10 differentially expressed mRNAs (CALCA1, CALCA2, CDR1, EHF, Cxorf56, PRKCA, INHBB, KRT13, MYADM, and PER2) and 10 lncRNAs (CDR1AS, XLOC_001265, ENSG00000238005.2, LOC636429, XLOC_005104, CRIP2, ENSG00000257497.1, BC062328, ENSG00000232220.2, and AB488780) were randomly selected for qRT-PCR. The consistent results were found between qRT-PCR and microarray data for 10 mRNAs and nine lncRNAs (Supplementary Figure ). Although the difference in the expression of the lncRNA AB488780 was not statistically significant by qRT-PCR, a consistent expression trend of this gene was observed between qRT-PCR and microarray (Figure ).

Validation for the expression of 10 significantly differential expressions of mRNAs and 10 lncRNAs by qRT-PCR. Three biological repeats were included in each gene. ∗P < 0.05.

GO and KEGG Pathway Analysis of Differentially Expressed mRNAs

To explore the potential biological functions of the differentially expressed mRNAs, the GO and KEGG pathway enrichments were carried out. The GO analysis indicated that these mRNAs were significantly enriched in chromosome organization, chromatin organization, organelle organization, and chromatin modification in BP; intracellular organelle, organelle, and membrane-bounded organelle in CC; chromatin binding, binding, transition metal ion binding, and protein binding in MF (Figure and Supplementary Table ). The pathways and molecular interactions associated with significantly differentially expressed mRNAs were then predicted by KEGG pathway enrichment analysis. The top 20 pathways were depicted in Figure , with a great number of mRNAs enriched into hedgehog signaling pathway, tight junction, and Wnt signaling pathway (Supplementary Table ).

Gene Ontology (GO) analysis of the differentially expressed mRNAs. Go annotation of differentially expressed mRNAs with top 10 enrichment scores covering domains of biological processes, cellular components, and molecular functions. The GO terms with corrected P-value less than 0.05 were considered significantly enriched by differentially expressed genes.

KEGG pathway analysis of the differentially expressed mRNAs with top 20 enrichment scores.

Co-expression of LncRNA and mRNA

To reveal the correlation between differentially expressed lncRNAs and mRNAs and figure out the possible mechanisms of lncRNAs during C. parvum infection, the co-expression network was constructed based on the mathematical relevance (Correlation > 0.99, Correlation < −0.99, and P-value < 0.05) to search similar expression profiles of lncRNAs and mRNAs. The co-expression network was constructed by using Cytoscape v3.6.0 (2015[4]). In this network, one mRNA could be correlated with one or more lncRNAs, and one lncRNA could also be associated with one or more mRNAs (Supplementary Table ). For example, the mRNA MYADM was related to the two lncRNAs (CDR1AS and ENSG00000238005.2), and the lncRNA LOC636429 corresponded to more than 300 mRNAs (Figure and Supplementary Table ).

Co-expression network of four significantly differentially expressed lncRNAs with their associate mRNAs. The network was based on the mathematical relevance (Correlation > 0.99, Correlation < –0.99, and P-value < 0.05) to search similar expression profiles of lncRNAs and mRNAs using cytoscape software (v3.6.0). The yellow triangle represents lncRNAs while the gradual color of red to green circular represents mRNAs. The black solid line indicates the correction of lncRNAs and mRNAs.

Functional Prediction of LncRNAs During C. parvum Infection

Long non-coding RNAs have been identified to function as regulators by cis- and trans-patterns (Huang et al., 2016). To predict the target genes of differentially expressed lncRNAs, the co-expressed neighboring coding genes located within 10 kb of these lncRNAs were selected for analysis. A total of 27 coding genes corresponding to 29 lncRNAs were predicted (cis-regulation). Additionally, 114 lncRNAs were also identified to indirectly regulate the expression of 109 distant genes through binding miRNAs (trans-regulation). All these cis- and trans-targets were predicted and listed in Supplementary Table . The GO enrichment analyses indicated that these target genes were significantly associated with BP in membrane raft organization, membrane organization, membrane assembly and regulation of signal transduction, CC in intracellular organelle, dihydrolipoyl dehydrogenase complex, tricarboxylic acid cycle enzyme complex andintracellular part, MF in binding, cytoskeletal protein binding, and histone threonine kinase activity (Figure and Supplementary Table ). The KEGG pathway enrichment analyses showed that these targets were also significantly enriched in hedgehog signaling pathway, tight junction, and Wnt signaling pathway (Figure and Supplementary Table ).

Gene Ontology (GO) analysis of the differentially expressed lncRNAs. Go annotation of differentially expressed lncRNAs with top 10 enrichment scores covering domains of biological processes, cellular components, and molecular functions. The GO terms with corrected P-value less than 0.05 were considered significantly enriched by differentially expressed genes.

KEGG pathway analysis of the differentially expressed lncRNAs with top 20 enrichment scores.

Discussion

Although some advances were achieved in biology, pathogenicity, and genetic characterization of Cryptosporidium (McDonald et al., 2001; Lee et al., 2017), no effective measures were developed to control cryptosporidiosis. The key challenge is that the interaction mechanism of host-Cryptosporidium has not been fully understood and appreciated. Unveiling the nature of host non-coding RNA world (e.g., lncRNA and miRNA) in last decades provided novel targets and strategies for preventing and treating infections of Theiler’s virus and Salmonella inflammatory bowel disease (IBD), diabetes, and multiple sclerosis in animals and humans (Atianand and Fitzgerald, 2014). In the present study, we systemically investigated the expression profiles of lncRNAs and mRNAs in HCT-8 cells infected with C. parvum IId subtype using microarray.

A total of 1,349 mRNAs were differentially expressed after the infection of C. parvum IId subtype. Among them, several inflammatory factors, e.g., IL-8, PTGS2, TCL-4, and CCL5 (RANTES), were up-regulated, and some genes associated with the cell proliferation and apoptosis were also significantly differentially expressed, e.g., up-regulated genes of thymidylate kinase, Cyclin A2, TM4SF1, IL1RN, Bcl2, and DUSP4, and down-regulated genes of Cyclin D1, Cyclin G2, BTG1, LAMB1, and LGALS1. These findings were consistent with previous studies using HCT-8 cells infected with C. parvum Iowa stain (IIaA15G2R1; Deng et al., 2004; Mele et al., 2004; Liu et al., 2009; Yang et al., 2015). Previous studies indicated that these genes would be involved in processes of inflammation, anti-apoptosis effect, and initiation and regulation of mucosal response during C. parvum, suggesting the important role of these genes within the interaction of host-Cryptosporidium. However, divergent expression patterns of mRNAs were also observed in HCT-8 cells infected with two subtypes of C. parvum. For example, the cell proliferation-related gene stratifin was down-regulated in our study, while it was up-regulated after the infection of C. parvum Iowa stain (Deng et al., 2004). Furthermore, the opposite trend of mRNA expression was detected for some negatively regulated genes (e.g., MT1B, MT1X, and MT1G) of apoptosis during infections of two different C. parvum subtypes (Liu et al., 2009), suggesting different pathogenic mechanisms of two subtypes.

Additionally, a total of 821 lncRNAs were found to be differentially expressed after the infection of C. parvum IIdA19G1. Co-expression and target prediction revealed 27 coding genes cis-regulated by 29 lncRNAs and 109 mRNAs trans-regulated by 114 lncRNAs. Function prediction of these differentially expressed transcripts was mainly involved in various pathways related to the infection and pathogenicity of Cryptosporidium, e.g., hedgehog signaling pathway, Wnt signaling pathway, and tight junction. In previous studies, the integrity of tight junction Zonula-Occludens-1 (ZO-1) was disrupted by Cryptosporidium infection (Buret et al., 2003). Wnt signaling plays a crucial role in the process of maintenance of intestinal epithelium. A particular hypothesis was undergoing investigated that Wnt signaling pathway was attenuated in intestinal epithelium infected with C. parvum (Zhang et al., 2016). The hedgehog signaling pathway showed negative effect within the Wnt signaling pathway and inhibition role in intestine proliferation (Katoh and Katoh, 2006). These findings suggested the possible regulating roles of host lncRNAs in these pathways during Cryptosporidium infection. Additionally, the qRT-PCR validation of 10 deregulated lncRNAs was consistent with the microarray data. Of them, the lncRNA XLOC_001265 was predicted to target RNF125. Previous studies have proved that mi-15b could regulate the Japanese Encephalitis Virus (JEV)-induced inflammatory cytokine (TNF-α, IL-1β, IL-6, CCL5, and IL-12p70) expression by targeting RNF125 (Zhu et al., 2015) in the JEV mouse model. Therefore, XLOC_001265 may be involved in the process of proinflammation caused by C. parvum IId subtype by regulating the expression of RNF125 because C. parvum had been proved to induce the expression of several inflammatory factors (IL-12, IL-17, IL-18, TNF-α, and TNF-γ; Ehigiator et al., 2005; Petry et al., 2010; O’Hara and Chen, 2011).

Conclusion

The expression profiles of mRNA and lncRNA were investigated in the present study, and a total of 1,349 mRNAs and 821 lncRNAs were significantly differentially expressed in the HCT-8 cells infected with C. parvum IId subtype. Co-expression analysis revealed that these differentially expressed lncRNAs would potentially cis- and trans-regulate the expression of mRNAs during C. parvum infection. Findings in the present study would provide novel insights for exploring the control measures for diagnosis and control of cryptosporidiosis in humans and animals.

Author Contributions

G-HZ conceived and designed the experiments. T-LL and X-CF performed the experiments and drafted the manuscript. T-LL, X-CF, Y-HL, Y-JY, and Y-LY analyzed the data. X-TW and L-XZ contributed to reagents and materials. All authors read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.