Skip to main content
Download PDF

Genome-wide analysis of differentially expressed profiles of mRNAs, lncRNAs and circRNAs during Cryptosporidium baileyi infection

BMC Genomics volume 19, Article number: 356 (2018) Cite this article

Abstract

Background

Cryptosporidium baileyi is the most common Cryptosporidium species in birds. However, effective prevention measures and treatment for C. baileyi infection were still not available. Long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) play important roles in regulating occurrence and progression of many diseases and are identified as effective biomarkers for diagnosis and prognosis of several diseases. In the present study, the expression profiles of host mRNAs, lncRNAs and circRNAs associated with C. baileyi infection were investigated for the first time.

Results

The tracheal tissues of experimental (C. baileyi infection) and control chickens were collected for deep RNA sequencing, and 545,479,934 clean reads were obtained. Of them, 1376 novel lncRNAs were identified, including 1161 long intergenic non-coding RNAs (lincRNAs) and 215 anti-sense lncRNAs. A total of 124 lncRNAs were found to be significantly differentially expressed between the experimental and control groups. Additionally, 14,698 mRNAs and 9085 circRNAs were identified, and significantly different expressions were observed for 1317 mRNAs and 104 circRNAs between two groups. Bioinformatic analyses of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway for their targets and source genes suggested that these dysregulated genes may be involved in the interaction between the host and C. baileyi.

Conclusions

The present study revealed the expression profiles of mRNAs, lncRNAs and circRNAs during C. baileyi infection for the first time, and sheds lights on the roles of lncRNAs and circRNAs underlying the pathogenesis of Cryptosporidium infection.

Background

Cryptosporidiosis is one of the most prevalent parasitic diseases in domestic and wild birds, and it has been reported in more than 30 avian species worldwide [1, 2]. Currently, four species of avian Cryptosporidium, namely Cryptosporidium meleagridis [3], Cryptosporidium baileyi [4], Cryptosporidium galli [5] and Cryptosporidium avium (previous avian genotype V) [6], and 12 genotypes (black duck genotype, Eurasian woodcock genotype, avian genotypes I–IV, VI and goose genotypes I–V) have been recognized [7, 8]. Among them, C. baileyi, the most common Cryptosporidium species in birds, has been reported in a wide range of avian hosts worldwide, especially in chickens [9, 10]. Infection of C. baileyi causes respiratory diseases, presenting symptoms of dyspnoea, coughing, rales and sneezing [4, 10, 11], and sometimes is associated with high morbidity and mortality, especially in broiler chickens [10], leading to considerable economic losses to the poultry industry [12,13,14]. In addition, C. baileyi has been observed in the stool of an immunodeficient man, suggesting the zoonotic potential of this species [15]. However, no effective prevention measures or treatments against C. baileyi infection have been developed [16].

In the process of attachment, internalization, and formation of parasitophorous vacuoles [17,18,19,20], Cryptosporidium molecules (eg. rhoptry protein, microneme protein, dense granules) were released into host cells [21, 22] and the transcriptome of host was dysregulated [23]. Seventy-three genes were upregulated and 74 were downregulated in human HCT-8 cells infected with C. parvum [24], and these differentially expressed genes are associated with cell communication, signal transduction and amino acid metabolic processes. Previous studies have indicated that non-coding RNA (ncRNA) molecules were also involved in host-Cryptosporidium interaction process, and aberrant expressions of microRNAs (miRNAs) were detected in C. parvum infecting epithelial cells [25, 26].

Recently, two new classes of ncRNA, namely long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), were discovered [27, 28], and they were found to play an important role in regulating occurrence and progression of many diseases [29,30,31]. LncRNAs are a group of non-coding RNAs longer than 200 nt [32] and dysregulated lncRNAs play important roles in progression of cells proliferation, invasiveness and metastasis of breast cancer [33], lung cancer [34] and colorectal cancer [35] through cis- and trans-regulation of gene expressions [36, 37] and miRNA sponges [38, 39]. In the abnormal activity of biological processes, such as cell proliferation, cell motility, immune, or inflammation response from diseases like cardiovascular disorders, inflammatory and autoimmune disease, lncRNAs also participate in the disease processes and contribute to controlling the gene regulatory network of host-pathogen interactions [40,41,42]. On the other hand, the function of circRNAs remains largely unknown. Previous studies have showed that it could harbor specific miRNAs as miRNA sponges [43,44,45] to suppresses miRNA activity, resulting in increased levels of miRNA targets, influence the expression of cytokines [46], promote cell cycles and inhibit cell apoptosis to act as a candidate oncogene [47, 48]. Additionally, both lncRNAs and circRNAs have been identified as effective biomarkers for diagnosis and prognosis of diseases [49, 50]. However, the genome-wide expression and functional roles of lncRNAs and circRNAs in Cryptosporidium infection are unclear. Here, we investigated the expression profiles of mRNAs, lncRNAs and circRNAs associated with C. baileyi infection for the first time. Our findings would provide baseline data for developing novel diagnostic and therapeutic targets for avian cryptosporidiosis.

Results

Identification of lncRNAs in chicken tracheal tissue

A total of 559,574,746 raw reads were produced from the Illumina HiSeq 4000 platform. After abandoned adaptor sequences and low-quality sequences, 545,479,934 clean reads (accounting for 81.82 Gb) were obtained, and the percentage of clean reads among raw data in each library ranged from 96.30 to 98.20%. Subsequently, the clean reads were mapped to the latest Gallus gallus reference genome. Firstly, considering the characteristics of lncRNA sequences (≥ 200 nt, exon count ≥2) and their differences from other classes of RNAs (eg. mRNAs, rRNAs, tRNAs, snRNAs, snoRNAs, pre-miRNAs, and pseudogenes), the transcripts were classified into different subtypes using both Scripture beta2 and Cufflinks (v2.1.1). 93.98% of the identified transcripts (50,857) were known as the reference transcripts, whereas 6.02% (3061) were the presumed lncRNAs. To further confirm these 3061 lncRNAs, 2823 transcripts were obtained after filter of the low expression transcripts with FPKM < 0.5. Finally, coding potential analysis was performed using the software CNCI, CPC, Pfam-scan, and PhyloCSF. After being screened by rigorous criteria and four analytic tools, a total of 1376 lncRNAs from the tracheal tissue of chickens were identified and subjected to further analysis (Fig. 1a). The 1376 lncRNAs were composed of 1161 (84.4%) long intergenic non-coding RNAs (lincRNAs) and 215 (15.6%) anti-sense lncRNAs (Additional file 1), nevertheless, other types of lncRNAs such as intronic lncRNAs, sense lncRNAs and bidirectional lncRNAs were not detected in our study (Fig. 1b).

Fig. 1
figure1

Screening and classification of the candidate lncRNAs in chicken trachea transcriptome. a Venn diagrams of coding potential analysis according to strict criteria. Four tools (CPC, CNCI, PFAM, and PhyloCSF) were used to analyze the coding potential of lncRNAs. Those simultaneously shared by four analytical tools were designated as candidate lncRNAs and used in subsequent analyses. b Classification of the two subtypes of lncRNAs

Full size image

Differentially expressed profiles of mRNAs, lncRNAs and circRNAs by RNA sequencing

The number of overlap mRNAs, lncRNAs and circRNAs in the experimental group compared with the control group is displayed in the Venn diagram, respectively (Fig. 2a–c). A total of 1317 mRNAs and 124 lncRNAs (Additional file 2) were found to be differentially expressed with the q-value < 0.05 and FDR < 0.05. Among them, 862 mRNAs were upregulated and 455 were downregulated in three infectious tracheal tissues compared with the controls (Fig. 3a). Meanwhile, 58 and 66 lncRNAs were upregulated and downregulated, respectively (Fig. 3b). According to the criteria of fold change > 2.0, 656 mRNAs were up-regulated and 292 were down-regulated. Similarly, 53 and 57 lncRNAs were up-regulated and down-regulated. For circRNAs, 104 remarkably differentially expressed genes (fold change > 2.0, p-value < 0.05 and FDR < 0.05) between two groups were identified (Additional file 2). Among them, 65 circRNAs were up-regulated and 39 were down-regulated (Fig. 3c).

Fig. 2
figure2

Venn diagram displaying the number of overlap mRNAs (a), lncRNAs (b) and circRNAs (c) in the experimental group (S group) compared with the control group (N group)

Full size image
Fig. 3
figure3

Differentially expressed mRNAs (a), lncRNAs (b) and circRNAs (c) in the two groups. Volcano plot of the p-values as a function of fold-change for mRNAs, lncRNAs and circRNAs indicate the dysregulated genes in the three normal and three infected tissues. Blue dots represent RNAs not significantly differentially expressed (p-value > 0.05) and the other dots represent RNAs differentially expressed (p-value < 0.05). Up-regulated RNAs were presented as red dots and down-regulated RNAs were presented as green dots

Full size image

In order to further understand the similarities among tracheas of chickens at the transcriptomic level, the data of all the differentially expressed genes were used for cluster analysis. The heat map clearly showed that the experimental group and the control group were separately clustered together (Fig. 4). The expression patterns of lncRNAs and circRNAs between the experimental group and the control group were distinguishable. These data suggested that the expressions of lncRNAs and circRNAs in chicken tissues infected with C. baileyi were significantly different from those in non-infected chicken tissues.

Fig. 4
figure4

Heatmap of expression profiles for the lncRNAs (a) and circRNAs (b). Red through blue color indicates high to low expression level. Each column represents one sample, and each row indicates a transcript

Full size image

Comparison of mRNAs and lncRNAs features

In this study, a total of 14,698 mRNAs, 1616 annotated lncRNAs and 1376 novel lncRNAs were identified from the tracheal tissues of chickens after infection with C. baileyi. To examine the differences among these three kinds of transcripts, comparative analysis of gene structure and sequence conservation was conducted. Our results showed that most lncRNAs contained two or three exons, which was significantly less than that of mRNAs (Fig. 5a). However, a slightly discrepancy was observed in the distribution of transcript length between mRNAs and lncRNAs (Fig. 5b). Compared with mRNAs, a relatively shorter open reading frame (ORF) was one of the main features of most lncRNAs (Fig. 5c). In addition, the majority of lncRNAs were less conserved which was different with mRNAs, although the difference is not statistically significant (Fig. 5d).

Fig. 5
figure5

Comparison of genomic architecture between mRNAs and lncRNAs. a The number of exons in the mRNAs and lncRNAs. In present study, single-exon lncRNAs were filtered out from the chicken genome due to the limitations of the algorithm. b Distribution of transcript lengths in the mRNAs and lncRNAs. The horizontal axis of indicates the length of transcripts, and the vertical axis represents density. c The number of open reading frames (ORFs) in the mRNAs and lncRNAs. d Conservation of the sequence in mRNAs and lncRNAs were evaluated using phastCons

Full size image

Validation of dysregulated mRNAs, lncRNAs and circRNAs

Ten mRNAs, ten lncRNAs and six circRNAs were randomly selected from the dysregulated genes to be verified by qRT-PCR. It was demonstrated that the expression levels from qRT-PCR data were consistent with the results of RNA sequencing (RNA-seq) (Fig. 6) thus verified the facticity of RNA-seq analysis. Consequently, efficient evidence was provided by this finding that these mRNAs, lncRNAs and circRNAs could be used to investigate the pathogenesis of C. baileyi infection in the future.

Fig. 6
figure6

Comparison of results of RNA sequencing and qRT-PCR. The blank column represents the fold change values obtained from RNA sequencing and the black column displays the fold change values of qRT-PCR. The figures represent the comparison results of mRNAs, lncRNAs and circRNAs from left to right. GAPDH was used as the internal control gene and three biological replicates were used. The expression trends of the results of RNA-seq were consistent with qRT-PCR

Full size image

Co-expression analysis and target prediction

Generally, the functions of both lncRNAs and circRNAs were performed by inter-acting with their targets [27, 28]. In our study, the potential cis- and trans-targets were predicted (Additional file 3). The mRNAs 100 kb upstream and downstream of the lncRNAs were searched for cis analysis, with 2308 lncRNAs that corresponded to 7679 protein-coding genes. The trans analysis of lncRNAs was performed by constructing co-expression networks (Fig. 7) of dysregulated mRNAs and lncRNAs based on expression correlation coefficient (Pearson correlation > 0.95 or < − 0.95) (Additional file 3). A total of 14,418 mRNAs and 2570 lncRNAs containing 634,780 relationships were detected. More than one mRNAs were predicted to be regulated by one lncRNA, and one mRNA corresponded to several lncRNAs (Fig. 7).

Fig. 7
figure7

Co-expression network of the representative lncRNAs and their partial target mRNAs. The relationships of mRNA-lncRNA were re-constructed based on the expression correlation coefficients (Pearson correlation > 0.95 or < − 0.95) using the cytoscape software (v3.6.0). Different colors were used to show different genes, with yellow for lncRNAs and red for mRNAs

Full size image

Previous studies have indicated that miRNA sponge was the common function of circRNAs [43, 44]. However, we have not examined the expression profile of miRNA in chicken tracheas during C. baileyi infection. Therefore, we just predicted the potential miRNA targets based on sequence complementarity in the present study, and the possible interaction relationships between 17,607 circRNAs and the 978 miRNAs in the chicken genome were observed (Additional file 4).

To further predict the function of mRNAs, lncRNAs and circRNAs from tracheal tissue of chickens, the GO (http://www.geneontology.org/) (Additional file 5) and KEGG (http://www.genome.jp/kegg/) pathway analyses (Additional file 6) were performed to analyze the dysregulated mRNAs and the target genes of differentially expressed lncRNAs and the source genes of circRNAs between the two groups.

In the present study, the top 10 results of GO enrichment analysis were selected as the master node of directed acyclic graph (DAG). The DAG of biological process (BP), cellular component (CC) and molecular function (MF) of differentially expressed mRNAs and lncRNAs are shown in Figs. 8a–c and 9a–c, respectively. Based on the GO analysis of differentially expressed mRNAs, the most significantly enriched BP were immune response, immune system process and positive regulation of immune system process, and external side of plasma membrane, cell surface and plasma membrane part were the most enriched CC. Furthermore, microtubule motor activity, motor activity and cytokine receptor activity were identified to be the most significantly enriched MF (Fig. 8d). On the basis of GO analysis of the target genes of differentially expressed lncRNAs, the most significantly enriched BP were immune system process, immune response and leukocyte activation, and the non-membrane-bounded organelle, intracellular non-membrane-bounded organelle and external side of plasma membrane were the most remarkably enriched CC. In addition, protein binding, purine nucleotide binding and adenyl nucleotide binding were the most significantly enriched MF (Fig. 9d). The most prominent category of gene function will be the focus of future research. Since the numbers of up-regulated and down-regulated circRNAs were small, there were no significant GO terms enriched in the two groups.

Fig. 8
figure8

GO analysis of mRNAs in trachea tissue of chickens. The significant enriched biological process, cellular component and molecular function with changed mRNAs in trachea tissue of chicken were displayed. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes (a–c). The number of genes in GO term were showed in histograph (d)

Full size image
Fig. 9
figure9

GO analysis of lncRNAs in trachea tissue of chickens. The significant enriched biological process, cellular component and molecular function with changed lncRNAs in trachea tissue of chicken were shown. Directed Acyclic Graph (DAG) is the graphical display of GO enrichment results with candidate targeted genes. The branch represents the relationship of inclusion, which defines the scope from increasingly small from top to bottom. The color depth represents the enrichment degree (a–c). The number of genes in GO term were showed in histograph (d)

Full size image

As a major public database of pathways, KEGG has been used to determine the significantly enrichment pathways for candidate target genes compared with the entire genome background [51, 52]. The top 20 pathways associated with mRNAs, lncRNAs and circRNAs were showed in an enriched scatter diagram, respectively. For mRNAs, the cytokine-cytokine receptor interaction, intestinal immune network for IgA production and cell adhesion molecules (CAMs) were identified as the top enriched KEGG pathways (Fig. 10a). KEGG pathway enrichment analysis for significantly differentially expressed lncRNAs revealed pathways for intestinal immune network for IgA production, cytokine-cytokine receptor interaction and cell cycle (Fig. 10b). Whereas, the significantly enriched pathways were amino sugar and nucleotide sugar metabolism, tight junction and glycerolipid metabolism for circRNAs (Fig. 10c).

Fig. 10
figure10

mRNAs (a), lncRNAs (b) and circRNAs (c) enriched KEGG pathway scatterplot showing statistics of pathway enrichment in the tracheal tissue of chickens between the two groups. The vertical axis represents the pathway name and the horizontal axis represents the rich factor. The size of the dot represents the number of differentially expressed genes in the pathway and the color of the point corresponds to the different q-value range

Full size image

Discussion

C. baileyi has been considered to be the most common avian Cryptosporidium species worldwide with the broadest host range [7, 53]. It can reside in the epithelial cells of respiratory tract, causing clinical respiratory disorders [4, 10, 11, 54] in birds (eg. chickens, turkeys and ducks). Considering the distinct morphological and biological features and large oocyst production, C. baileyi has been suggested as a model for the characterization of cryptosporidia [55, 56]. Moreover, C. baileyi is able to establish its infection in the mucosal epithelium of a wide variety of organs (eg. respiratory tract, bursa of Fabricius and cloaca) [4, 53], and the trachea is the most common predilection site to present the inflammation and clinical signs. Therefore, the trachea from the experimental and control groups were selected to study the differentially expressed profiles of mRNAs, lncRNAs and circRNAs. The results of our study indicated that mRNAs and ncRNAs may play a significant role in the infection process of C. baileyi, representing a potential therapeutic target for C. baileyi infection as revealed by the sequencing analysis.

The chicken genome has been sequenced and assembled using the whole-genome mapping technology in 2004 [57, 58], with high-quality genome assembly and annotation. The novel mRNAs, lncRNAs and circRNAs obtained in our study greatly improved the annotation of the chicken genome. A total of 14,698 mRNAs were obtained from the tracheal tissue of chickens infected with C. baileyi and 1317 were significantly dysregulated. The pathways of cytokine-cytokine receptor interaction, CAMs and toll-like receptor signaling were significantly enriched by GO and KEGG enrichment analyses of these differentially expressed mRNAs, and previous studies also indicated that the infection of zoonotic C. parvum was closely associated with these pathways in the epithelial cells [59,60,61], indicating that these dysregulated genes may be involved in the interactions between the host and C. baileyi.

The potential functions of lncRNAs were commonly predicted by their target genes. In our study, the relationships between lncRNAs and their target genes co-located (within 100 kb) and co-expressed (the Pearson’s correlation coefficients > 0.95 or < − 0.95) were analyzed (Additional file 3, Fig. 7). Analyzing potential target genes of dysregulated lncRNAs revealed their important roles in regulation of the interaction between chickens and C. baileyi. For example, one lncRNA, ALDBGALG0000001306, was significantly upregulated between the infected and non-infected chickens, and its potential targets IL-6, IL-12 and IL-17 have been identified as important cytokines against Cryptosporidium infection [53, 62,63,64], suggesting that this lncRNA would be a regulator of immune response against C. baileyi infection in chickens.

Additionally, circRNAs, a newly discovered class of ncRNAs, may act as miRNA sponges to regulate the activity of target genes and participates in the regulation of gene transcription and protein production [65, 66] to affect prognosis of diseases, especially tumors [47, 49]. Therefore, circRNAs were tested in our study using RNA sequencing and only 104 circRNAs genes (65 were upregulated and 39 were downregulated) were demonstrated to be differentially expressed in response to the infection of C. baileyi. Compared with mRNAs and lncRNAs, the expression levels of circRNAs were generally lower (Fig. 6) and no significant GO terms were enriched in our study. The reasons for these differences may be due to the drawbacks of RNA sequencing protocol in our study or the expression of circRNA might be closely related to immune-specific organs and tissues (eg. brain, stem cells, testis and discs) predicted in previous studies [67,68,69]. The source genes of these circRNAs were mainly enriched in the metabolic pathway of amino sugar and nucleotide sugar metabolism, tight junction and glycerolipid metabolism. Some miRNAs (eg. let-7d, let-7f and let-7i) have been proven to regulate the expression of toll-like receptor four signal in respond to C. parvum infection [25, 26, 70, 71], and their interaction circRNAs (Additional file 4) were also found to be differentially expressed, suggesting that these circRNAs may implicate in the pathogenesis and progression of C. baileyi infection.

Conclusions

Differentially expressed mRNAs, lncRNAs and circRNAs were screened from chickens after C. baileyi infection. A total of 1317 mRNAs, 124 lncRNAs and 104 circRNAs were differentially expressed, and these RNAs could regulate expression of their related genes and play a key role in the pathogenesis of C. baileyi infection. The relevant signaling pathway of these predicted mRNAs and ncRNAs would be the focus of future study to fully reveal the pathogenesis of C. baileyi infection. The current work also provides new insight into the pathways and mechanisms that mediate the host immune response to C. baileyi (albeit in birds only).

Methods

C. baileyi infection model

A total of six newly hatched chickens were purchased from the Giant Long Company (Shaanxi, China), and reared in a pathogen-free laboratory and given constant light and adequate food and water during the entire experimental period. Chickens were divided into two groups randomly with three chickens per group and each chicken was kept in a separate cage. After normal feeding for 3 days, each chicken in the experimental group (S group) was orally infected with 1 × 106 C. baileyi oocysts, while the control group (N group) was given the same volume of PBS. Faecal samples were collected daily post infection (pi). The hemacytometer was used to record the numbers of oocyst per gram (OPG) of faeces. The successful infection was confirmed by the oocyst excretion (Additional file 7) and parasites in histological observations by haemotoxylin & eosinstaining and by electron microscopy (Additional file 8).

Sample collection and preparation

Since C. baileyi mainly parasitizes in the respiratory tract, causing a series of respiratory diseases, at the first peak of oocyst excretion (10 dpi), the tracheas of three chickens in both experimental and control groups were collected and sent to the company (Novogene, Beijing) for sequencing. Chickens were firstly kept into inverted beakers with infiltrated ether ball for several minutes. When they presented symptoms of slow breathing, absent corneal reflexes and skin sensation, these animals were then anesthetized. The trachea from each chicken was rapidly sampled and washed with nuclease-free PBS. All experiments were permitted by ethics committee of Northwest A&F University and conducted in the fume hood. All the instruments and reagents were treated with DEPC in advance and without RNase.

RNA isolation, library preparation and sequencing

Total RNA was extracted from each trachea sample using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The concentration and purity of RNA were examined using the Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA) and the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), respectively. The RNA integrity was measured using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 3 μg RNA per sample was then used for RNA sequencing. Firstly, the ribosomal RNA (rRNA) was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA). Subsequently, the RNAs were subjected to generate sequencing libraries by NEBNext® Ultra™ Directional RNA Library Prep Kit (NEB, USA). After that, the first and second strand of the complementary DNA (cDNA) were synthesized and the AMPure XP system (Beckman Coulter, Beverly, USA) was used to purify the library fragments for firstly selecting cDNA fragments with a length of 150–200 bp. PCR was then performed and the Agilent Bioanalyzer 2100 system was used to evaluate the quality of the library. HiSeq PE Cluster Kit v4 cBot (Illumina) was used to perform the clustering of index-coded samples according to the manufacturer’s instructions. The libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were produced after cluster generation. For circRNA sequencing, the RNase R (Epicentre, USA) was additionally used to deal with the rRNA-depleted RNAs and remove the linear RNAs before sequencing library generation.

Quality control

The raw reads were processed firstly by the in-house perl scripts. To obtain clean reads, low quality reads, containing adapter and ploy-N reads were removed from raw data in this step. Meanwhile, the GC content, Q20 and Q30 of the clean data were also calculated. The clean data with high quality was used for further analysis.

Mapping to the reference genome

The reference genome (ftp://ftp.ensembl.org/pub/release-83/fasta/gallus_gallus/dna/) and gene model annotation files (ftp://ftp.ensembl.org/pub/release-83/gtf/gallus_gallus/) were downloaded from website directly. Bowtie (v2.0.6) [72] was used to build index of the reference genome, and TopHat (v2.0.9) [73] was used to align paired-end clean reads to the reference genome. In addition, the software find_circ [74] was used to extend the anchor sequences and the back-spliced reads containing at least two supporting reads were considered to be circRNAs.

Transcriptome assembly of lncRNAs

The softwares Scripture beta2 [75] and Cufflinks (v2.1.1) [76, 77] were used to assemble mapped reads of each sample. These two methods determined exons connectivity by using spliced reads in different ways. With setting other parameters as default, Scripture and Cufflinks were run with ‘min-frags-per-transfrag = 0’ and ‘–library-type’, respectively.

LncRNAs coding potential and conservation analysis

Four analytic tools, namely coding-non-coding-index (CNCI) [78], coding potential calculator (CPC) [79], phylogenetic codon substitution frequency (PhyloCSF) [80] and Pfam-scan [81], were used to effectively distinguish protein-coding and non-coding sequences. Ultimately, the transcripts predicted with coding potential by any of the four tools were filtered out, and those without coding potential would to be our candidate set of lncRNAs.

The software phast (v1.3) [82] was usually used for phylogenetic analysis to evaluate the sequence conservation of transcripts. Then, we used the program phyloFit to calculate phylogenetic models of the conserved and non-conserved regions among species. PhastCons was used to compute the conservation scores of coding genes and lncRNAs.

Target gene prediction and differential expression analysis

Both cis and trans roles of target genes for lncRNA were predicted. Here, we searched coding genes within 100 kb upstream and downstream of each lncRNA and then analyzed their function. Fragments per kilobase for a million reads (FPKMs) of both lncRNAs and coding genes in each sample [77] were calculated by the Cuffdiff (v2.1.1). Unlike the former, psRobot [83] was used to predict the miRNA binding sites of circRNAs. Besides, the transcript per million (TPM) was used to normalize the expression level of circRNAs according to the criteria described by Zhou et al. [84]. Between two groups of chickens, transcripts or genes with p-value adjusted by the Benjamini & Hochberg method (q-value) < 0.05 were considered significant.

GO and KEGG enrichment analysis

To explore the function of mRNAs, lncRNAs and circRNAs, the Gene Ontology (GO) and Kyoko Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were conducted. The GOseq (v1.18.0) was used to perform the GO enrichment analysis of differentially expressed genes or target genes of lncRNAs or source genes of differentially expressed circRNAs, in which gene length bias was corrected. All three Go categories, namely cellular component, biological process, and molecular function, were included, and GO terms with the q-value < 0.05 was considered significantly enriched. KOBAS [51] software was used to examine the statistical enrichment analysis of differential expression genes or lncRNA target genes or source genes of differentially expressed circRNAs in KEGG pathways. The gene sets were firstly mapped to database genes, and then compared with reference chicken genome (ftp://ftp.ensembl.org/pub/release-83/fasta/gallus_gallus/dna/) to scan enriched pathways, diseases and functions. The enriched information was then evaluated by the statistical test and correction. The EASE score was calculated to test the relevance, and p-value < 0.05 was considered significantly enriched by differentially expressed genes.

Validation by quantitative real-time polymerase chain reaction (qRT-PCR)

Three trachea RNA samples from both infected and non-infected chickens were analyzed by qRT-PCR. Total cDNA of each sample was synthesized using two-step reverse transcriptase Kit (Vazyme Biotech Co., Ltd., Nanjing, China) according to the manufacturer’s instructions. qRT-PCR were performed using LightCycler 480II Real-Time PCR System (Roche, Indianapolis, Indiana) and UltraSYBR Mixture (Qiagen, Shanghai, China). Each reaction (in 10 μL) contained 5 μL 2 × QuantiFast® SYBR® Green PCR Master Mix, 0.4 μL primers (5 μM) (Additional file 9), and 1 μL cDNA. The protocol included an initial single cycle denaturing at 95 °C for 10 min, followed by 40 cycles of denaturing at 95 °C for 10 s and annealing at 60 °C or 55 °C (Additional file 9) for 30 s. All amplifications were followed by dissociation curve analysis of the amplified products. Specific primers were designed using the Primer Premier (v5.0), and specificities were confirmed with BLAST. The gene expression levels were normalized with the reference gene GAPDH by using 2−ΔΔCt value methods. The data were presented as the means ± SEM.

Statistical analysis

The Pearson’s correlation coefficients were used to calculate expression correlation between lncRNAs and mRNAs (r > 0.95 or r < − 0.95). The statistical difference in gene expression of qRT-PCR results was analyzed by Student’s t-test and false discovery rate (FDR) was also calculated to correct the p-value. It was considered to be statistically significant when p-value < 0.05 or q-value < 0.05.

Abbreviations

BP:

Biological process

C. baileyi :

Cryptosporidium baileyi

C. parvum :

Cryptosporidium parvum

CAMs:

Cell adhesion molecules

CC:

Cellular component

cDNA:

Complementary DNA

circRNAs:

Circular RNAs

CNCI:

Coding-non-coding-index

CPC:

Coding potential calculator

DAG:

Directed acyclic graph

FDR:

False discovery rate

FPKMs:

Fragments per kilobase for a million reads

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

lincRNAs:

Long intergenic non-coding RNAs

lncRNAs:

Long non-coding RNAs

MF:

Molecular function

miRNA:

microRNA

ncRNA:

Non-coding RNA

OPG:

Oocyst per gram

ORF:

Open reading frame

PhyloCSF:

Phylogenetic codon substitution frequency

pi:

Post infection

qRT-PCR:

Quantitative real-time polymerase chain reaction

RNA-seq:

RNA sequencing

rRNA:

Ribosomal RNA

SEM:

Standard error of mean

TPM:

Transcript per million

References

  1. 1.

    Amer S, Wang C, He H. First detection of Cryptosporidium baileyi in Ruddy Shelduck (Tadorna ferruginea) in China. J Vet Med Sci. 2010;72(7):935–8.

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Qi M, Wang R, Ning C, Li X, Zhang L, Jian F, Sun Y, Xiao L. Cryptosporidium spp. in pet birds: genetic diversity and potential public health significance. Exp Parasitol. 2011;128(4):336–40.

    Article  PubMed  Google Scholar 

  3. 3.

    Slavin D. Cryptosporidium meleagridis (sp. nov.). J Comp Pathol. 1955;65(3):262–6.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Current WL, Upton SJ, Haynes TB. The life cycle of Cryptosporidium baileyi n. Sp. (Apicomplexa, Cryptosporidiidae) infecting chickens. J Protozool. 1986;33(2):289–96.

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Ryan UM, Xiao L, Read C, Sulaiman IM, Monis P, Lal AA, Fayer R, Pavlasek I. A redescription of Cryptosporidium galli Pavlásek, 1999 (Apicomplexa: Cryptosporidiidae) from birds. J Parasitol. 2003;89(4):809–13.

    CAS  Article  PubMed  Google Scholar 

  6. 6.

    Holubová N, Sak B, Horčičková M, Hlásková L, Květoňová D, Menchaca S, McEvoy J, Kváč M. Cryptosporidium avium n. sp. (Apicomplexa: Cryptosporidiidae) in birds. Parasitol Res. 2016;115(6):2243–51.

    Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Nakamura AA, Meireles MV. Cryptosporidium infections in birds – a review. Rev Bras Parasitol Vet. 2015;24(3):253–67.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Chelladurai JJ, Clark ME, Kváč M, Holubová N, Khan E, Stenger BL, Giddings CW, McEvoy J. Cryptosporidium galli and novel Cryptosporidium avian genotype VI in North American red-winged blackbirds (Agelaius phoeniceus). Parasitol Res. 2016;115(5):1901–6.

    Article  PubMed  Google Scholar 

  9. 9.

    Huber F, da Silva S, Bomfim TC, Teixeira KR, Bello AR. Genotypic characterization and phylogenetic analysis of Cryptosporidium sp. from domestic animals in Brazil. Vet Parasitol. 2007;150(1–2):65–74.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Ryan U. Cryptosporidium in birds, fish and amphibians. Exp Parasitol. 2010;124(1):113–20.

    Article  PubMed  Google Scholar 

  11. 11.

    Molina-Lopez RA, Ramis A, Martin-Vazquez S, Gomez-Couso H, Ares-Mazas E, Caccio SM, Leiva M, Darwich L. Cryptosporidium baileyi infection associated with an outbreak of ocular and respiratory disease in otus owls (Otus scops) in a rehabilitation Centre. Avian Pathol. 2010;39(3):171–6.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Abbassi H, Coudert F, Chérel Y, Dambrine G, Brugère-Picoux J, Naciri M. Renal cryptosporidiosis (Cryptosporidium baileyi) in specific-pathogen-free chickens experimentally coinfected with Marek’s disease virus. Avian Dis. 1999;43(4):738–44.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Blagburn BL, Lindsay DS, Hoerr FJ, Davis JF, Giambrone JJ. Pathobiology of cryptosporidiosis (C. baileyi) in broiler chickens. J Protozool. 1991;38(6):25S–8S.

    CAS  PubMed  Google Scholar 

  14. 14.

    Goodwin MA, Brown J. Effect of Cryptosporidium baileyi on broilers infected at 26 days of age. Avian Dis. 1990;34(2):458–62.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Ditrich O, Palkovic L, Stĕrba J, Prokopic J, Loudová J, Giboda M. The first finding of Cryptosporidium baileyi in man. Parasitol Res. 1991;77(1):44–7.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Cabada MM, White AC Jr. Treatment of cryptosporidiosis: do we know what we think we know? Curr Opin Infect Dis. 2010;23(5):494–9.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Bonnin A, Lapillonne A, Petrella T, Lopez J, Chaponnier C, Gabbiani G, Robine S, Dubremetz JF. Immunodetection of the microvillous cytoskeleton molecules villin and ezrin in the parasitophorous vacuole wall of Cryptosporidium parvum (Protozoa: Apicomplexa). Eur J Cell Biol. 1999;78(11):794–801.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    De Souza W. Microscopy and cytochemistry of the biogenesis of the parasitophorous vacuole. Histochem Cell Biol. 2005;123(1):1–18.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Huang BQ, Chen XM, LaRusso NF. Cryptosporidium parvum attachment to and internalization by human biliary epithelia in vitro: a morphologic study. J Parasitol. 2004;90(2):212–21.

    Article  PubMed  Google Scholar 

  20. 20.

    Zhang H, Guo F, Zhu G. Cryptosporidium lactate dehydrogenase is associated with the Parasitophorous vacuole membrane and is a potential target for developing therapeutics. PLoS Pathog. 2015;11(11):e1005250.

    Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Sam-Yellowe TY. Rhoptry organelles of the apicomplexa: their role in host cell invasion and intracellular survival. Parasitol Today. 1996;12(8):308–16.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Valentini E, Cherchi S, Possenti A, Dubremetz JF, Pozio E, Spano F. Molecular characterisation of a Cryptosporidium parvum rhoptry protein candidate related to the rhoptry neck proteins TgRON1 of Toxoplasma gondii and PfASP of Plasmodium falciparum. Mol Biochem Parasitol. 2012;183(1):94–9.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    Mauzy MJ, Enomoto S, Lancto CA, Abrahamsen MS, Rutherford MS. The Cryptosporidium parvum transcriptome during in vitro development. PLoS One. 2012;7(3):e31715.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Yang YL, Serrano MG, Sheoran AS, Manque PA, Buck GA, Widmer G. Over-expression and localization of a host protein on the membrane of Cryptosporidium parvum infected epithelial cells. Mol Biochem Parasitol. 2009;168(1):95–101.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Zhou R, Hu G, Liu J, Gong AY, Drescher KM, Chen XM. NF-kappaB p65-dependent transactivation of miRNA genes following Cryptosporidium parvum infection stimulates epithelial cell immune responses. PLoS Pathog. 2009;5(12):e1000681.

    Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Zhou R, Feng Y, Chen XM. Non-coding RNAs in epithelial immunity to Cryptosporidium infection. Parasitology. 2014;141(10):1233–43.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Bhan A, Mandal SS. Long noncoding RNAs: emerging stars in gene regulation, epigenetics and human disease. ChemMedChem. 2014;9(9):1932–56.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365(2):141–8.

    CAS  Article  PubMed  Google Scholar 

  29. 29.

    Hirata H, Hinoda Y, Shahryari V, Deng G, Nakajima K, Tabatabai ZL, Ishii N, Dahiya R. Long noncoding RNA MALAT1 promotes aggressive renal cell carcinoma through Ezh2 and interacts with miR-205. Cancer Res. 2015;75(7):1322–31.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Liu B, Sun L, Liu Q, Gong C, Yao Y, Lv X, Lin L, Yao H, Su F, Li D, Zeng M, Song E. A cytoplasmic NF-κB interacting long noncoding RNA blocks IκB phosphorylation and suppresses breast cancer metastasis. Cancer Cell. 2015;27(3):370–81.

    CAS  Article  PubMed  Google Scholar 

  31. 31.

    Yuan JH, Yang F, Wang F, Ma JZ, Guo YJ, Tao QF, Liu F, Pan W, Wang TT, Zhou CC, Wang SB, Wang YZ, Yang Y, Yang N, Zhou WP, Yang GS, Sun SH. A long noncoding RNA activated by TGF-β promotes the invasion-metastasis cascade in hepatocellular carcinoma. Cancer Cell. 2014;25(5):666–81.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Chen Z, Luo Y, Yang W, Ding L, Wang J, Tu J, Geng B, Cui Q, Yang J. Comparison analysis of dysregulated LncRNA profile in mouse plasma and liver after hepatic ischemia/reperfusion injury. PLoS One. 2015;10(7):e0133462.

    Article  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Gupta RA, Shah N, Wang KC, Kim J, Horlings HM, Wong DJ, Tsai MC, Hung T, Argani P, Rinn JL, Wang Y, Brzoska P, Kong B, Li R, West RB, van de Vijver MJ, Sukumar S, Chang HY. Long non-coding RNA HOTAIR reprograms chromatin state to promote cancer metastasis. Nature. 2010;464(7291):1071–6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Ma MZ, Chu BF, Zhang Y, Weng MZ, Qin YY, Gong W, Quan ZW. Long non-coding RNA CCAT1 promotes gallbladder cancer development via negative modulation of miRNA-218-5p. Cell Death Dis. 2015;6:e1583.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Shi Y, Liu Y, Wang J, Jie D, Yun T, Li W, Yan L, Wang K, Feng J. Downregulated long noncoding RNA BANCR promotes the proliferation of colorectal Cancer cells via downregulation of p21 expression. PLoS One. 2015;10(4):e0122679.

    Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    He Y, Meng XM, Huang C, Wu BM, Zhang L, Lv XW, Li J. Long noncoding RNAs: novel insights into hepatocellular carcinoma. Cancer Lett. 2014;344(1):20–7.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Lee S, Kopp F, Chang TC, Sataluri A, Chen B, Sivakumar S, Yu H, Xie Y, Mendell JT. Noncoding RNA NORAD regulates genomic stability by sequestering PUMILIO proteins. Cell. 2016;164(1–2):69–80.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Liu W, Ma C, Yang B, Yin C, Zhang B, Xiao Y. LncRNA Gm15290 sponges miR-27b to promote PPARγ-induced fat deposition and contribute to body weight gain in mice. Biochem Biophys Res Commun. 2017;493(3):1168–75.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Yu SY, Dong B, Zhou SH, Tang L. LncRNA UCA1 modulates cardiomyocyte apoptosis by targeting miR-143 in myocardial ischemia-reperfusion injury. Int J Cardiol. 2017;247:31.

    Article  PubMed  Google Scholar 

  40. 40.

    Li Z, Rana TM. Decoding the noncoding: prospective of lncRNA-mediated innate immune regulation. RNA Biol. 2014;11(8):979–85.

    Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Mirza AH, Kaur S, Brorsson CA, Pociot F. Effects of GWAS-associated genetic variants on lncRNAs within IBD and T1D candidate loci. PLoS One. 2014;9(8):e105723.

    Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Yan B, Yao J, Liu JY, Li XM, Wang XQ, Li YJ, Tao ZF, Song YC, Chen Q, Jiang Q. lncRNA-MIAT regulates microvascular dysfunction by functioning as a competing endogenous RNA. Circ Res. 2015;116(7):1143–56.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Li H, Yang J, Wei X, Song C, Dong D, Huang Y, Lan X, Plath M, Lei C, Ma Y, Qi X, Bai Y, Chen H. CircFUT10 reduces proliferation and facilitates differentiation of myoblasts by sponging miR-133a. J Cell Physiol. 2018;233(6):4643-51.

  45. 45.

    Wilusz JE, Sharp PA. Molecular biology. A circuitous route to noncoding RNA. Science. 2013;340(6131):440–1.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Deng T, Yang L, Zheng Z, Li Y, Ren W, Wu C, Guo L. Calcitonin gene-related peptide induces IL-6 expression in RAW264.7 macrophages mediated by mmu_circRNA_007893. Mol Med Rep. 2017;16(6):9367–74.

    Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Kong Z, Wan X, Zhang Y, Zhang P, Zhang Y, Zhang X, Qi X, Wu H, Huang J, Li Y. Androgen-responsive circular RNA circSMARCA5 is up-regulated and promotes cell proliferation in prostate cancer. Biochem Biophys Res Commun. 2017;493(3):1217–23.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Lai Z, Yang Y, Yan Y, Li T, Li Y, Wang Z, Shen Z, Ye Y, Jiang K, Wang S. Analysis of co-expression networks for circular RNAs and mRNAs reveals that circular RNAs hsa_circ_0047905, hsa_circ_0138960 and has-circRNA7690-15 are candidate oncogenes in gastric cancer. Cell Cycle. 2017;16(23):2301-11.

  49. 49.

    Cui X, Niu W, Kong L, He M, Jiang K, Chen S, Zhong A, Li W, Lu J, Zhang L. hsa_circRNA_103636: potential novel diagnostic and therapeutic biomarker in major depressive disorder. Biomark Med. 2016;10(9):943–52.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Zhou M, Zhao H, Wang Z, Cheng L, Yang L, Shi H, Yang H, Sun J. Identification and validation of potential prognostic lncRNA biomarkers for predicting survival in patients with multiple myeloma. J Exp Clin Cancer Res. 2015;34:102.

    Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Mao X, Cai T, Olyarchuk JG, Wei L. Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics. 2005;21(19):3787–93.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36(Database issue):D480–4.

    CAS  PubMed  Google Scholar 

  53. 53.

    Zhao GH, Cheng WY, Wang W, Jia YQ, Fang YQ, Du SZ, Yu SK. The expression dynamics of IL-17 and Th17 response relative cytokines in the trachea and spleen of chickens after infection with Cryptosporidium baileyi. Parasit Vectors. 2014;7:212.

    Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Wang R, Jian F, Sun Y, Hu Q, Zhu J, Wang F, Ning C, Zhang L, Xiao L. Large-scale survey of Cryptosporidium spp. in chickens and Pekin ducks (Anas platyrhynchos) in Henan, China: prevalence and molecular characterization. Avian Pathol. 2010;39(6):447–51.

    Article  PubMed  Google Scholar 

  55. 55.

    Egyed Z, Sréter T, Széll Z, Beszteri B, Dobos-Kovács M, Márialigeti K, Cornelissen AW, Varga I. Polyphasic typing of Cryptosporidium baileyi: a suggested model for characterization of cryptosporidia. J Parasitol. 2002;88(2):237–43.

    CAS  Article  PubMed  Google Scholar 

  56. 56.

    Huang L, Zhu H, Zhang S, Wang R, Liu L, Jian F, Ning C, Zhang L. An in vitro model of infection of chicken embryos by Cryptosporidium baileyi. Exp Parasitol. 2014;147:41–7.

    Article  PubMed  Google Scholar 

  57. 57.

    International Chicken Genome Sequencing Consortium. Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature. 2004;432(7018):695–716.

    Article  Google Scholar 

  58. 58.

    Wong GK, Liu B, Wang J, Zhang Y, Yang X, Zhang Z, Meng Q, Zhou J, Li D, Zhang J, Ni P, Li S, Ran L, Li H, Zhang J, Li R, Li S, Zheng H, Lin W, Li G, Wang X, Zhao W, Li J, Ye C, Dai M, Ruan J, Zhou Y, Li Y, He X, Zhang Y, Wang J, Huang X, Tong W, Chen J, Ye J, Chen C, Wei N, Li G, Dong L, Lan F, Sun Y, Zhang Z, Yang Z, Yu Y, Huang Y, He D, Xi Y, Wei D, Qi Q, Li W, Shi J, Wang M, Xie F, Wang J, Zhang X, Wang P, Zhao Y, Li N, Yang N, Dong W, Hu S, Zeng C, Zheng W, Hao B, Hillier LW, Yang SP, Warren WC, Wilson RK, Brandström M, Ellegren H, Crooijmans RP, van der Poel JJ, Bovenhuis H, Groenen MA, Ovcharenko I, Gordon L, Stubbs L, Lucas S, Glavina T, Aerts A, Kaiser P, Rothwell L, Young JR, Rogers S, Walker BA, van Hateren A, Kaufman J, Bumstead N, Lamont SJ, Zhou H, Hocking PM, Morrice D, de Koning DJ, Law A, Bartley N, Burt DW, Hunt H, Cheng HH, Gunnarsson U, Wahlberg P, Andersson L, Kindlund E, Tammi MT, Andersson B, Webber C, Ponting CP, Overton IM, Boardman PE, Tang H, Hubbard SJ, Wilson SA, Yu J, Wang J, Yang H, International Chicken Polymorphism Map Consortium. A genetic variation map for chicken with 2.8 million single-nucleotide polymorphisms. Nature. 2004;432(7018):717–22.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Lean IS, McDonald V, Pollok RC. The role of cytokines in the pathogenesis of Cryptosporidium infection. Curr Opin Infect Dis. 2002;15(3):229–34.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    McDonald V, Korbel DS, Barakat FM, Choudhry N, Petry F. Innate immune responses against Cryptosporidium parvum infection. Parasite Immunol. 2013;35(2):55–64.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Yang Z, Fu Y, Gong P, Zheng J, Liu L, Yu Y, Li J, Li H, Yang J, Zhang X. Bovine TLR2 and TLR4 mediate Cryptosporidium parvum recognition in bovine intestinal epithelial cells. Microb Pathog. 2015;85:29–34.

    CAS  Article  PubMed  Google Scholar 

  62. 62.

    Ehigiator HN, Romagnoli P, Borgelt K, Fernandez M, McNair N, Secor WE, Mead JR. Mucosal cytokine and antigen-specific responses to Cryptosporidium parvum in IL-12p40 KO mice. Parasite Immunol. 2005;27(1–2):17–28.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Pagani M, Rossetti G, Panzeri I, de Candia P, Bonnal RJ, Rossi RL, Geginat J, Abrignani S. Role of microRNAs and long-non-coding RNAs in CD4(+) T-cell differentiation. Immunol Rev. 2013;253(1):82–96.

    Article  PubMed  Google Scholar 

  64. 64.

    Petry F, Jakobi V, Tessema TS. Host immune response to Cryptosporidium parvum infection. Exp Parasitol. 2010;126(3):304–9.

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    Chen LL, Yang L. Regulation of circRNA biogenesis. RNA Biol. 2015;12(4):381–8.

    Article  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Li J, Yang J, Zhou P, Le Y, Zhou C, Wang S, Xu D, Lin HK, Gong Z. Circular RNAs in cancer: novel insights into origins, properties, functions and implications. Am J Cancer Res. 2015;5(2):472–80.

    CAS  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Ma CJ, Liu X, Che L, Liu ZH, Samartzis D, Wang HQ. Stem cell therapies for intervertebral disc degeneration: immune privilege reinforcement by Fas/FasL regulating machinery. Curr Stem Cell Res Ther. 2015;10(4):285–95.

    CAS  Article  PubMed  Google Scholar 

  68. 68.

    Liu ZH, Sun Z, Wang HQ, Ge J, Jiang TS, Chen YF, Ma Y, Wang C, Hu S, Samartzis D, Luo ZJ. FasL expression on human nucleus pulposus cells contributes to the immune privilege of intervertebral disc by interacting with immunocytes. Int J Med Sci. 2013;10(8):1053–60.

    Article  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Lan PH, Liu ZH, Pei YJ, Wu ZG, Yu Y, Yang YF, Liu X, Che L, Ma CJ, Xie YK, Hu QJ, Wan ZY, Wang HQ. Landscape of RNAs in human lumbar disc degeneration. Oncotarget. 2016;7(39):63166–76.

    Article  PubMed  PubMed Central  Google Scholar 

  70. 70.

    Chen XM, Splinter PL, O'Hara SP, LaRusso NF. A cellular micro-RNA, let-7i, regulates toll-like receptor 4 expression and contributes to cholangiocyte immune responses against Cryptosporidium parvum infection. J Biol Chem. 2007;282(39):28929–38.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  71. 71.

    Hu G, Gong AY, Roth AL, Huang BQ, Ward HD, Zhu G, Larusso NF, Hanson ND, Chen XM. Release of luminal exosomes contributes to TLR4-mediated epithelial antimicrobial defense. PLoS Pathog. 2013;9(4):e1003261.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  Google Scholar 

  74. 74.

    Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, Kocks C, le Noble F, Rajewsky N. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    CAS  Article  PubMed  Google Scholar 

  75. 75.

    Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Adiconis X, Fan L, Koziol MJ, Gnirke A, Nusbaum C, Rinn JL, Lander ES, Regev A. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs. Nat Biotechnol. 2010;28(5):503–10.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn JL. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011;25(18):1915–27.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  77. 77.

    Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.

    Article  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Lin MF, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics. 2011;27(13):i275–82.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, Pang N, Forslund K, Ceric G, Clements J, Heger A, Holm L, Sonnhammer EL, Eddy SR, Bateman A, Finn RD. The Pfam protein families database. Nucleic Acids Res. 2012;40(Database issue):D290–301.

    CAS  Article  PubMed  Google Scholar 

  82. 82.

    Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, Weinstock GM, Wilson RK, Gibbs RA, Kent WJ, Miller W, Haussler D. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005;15(8):1034–50.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  83. 83.

    Wu HJ, Ma YK, Chen T, Wang M, Wang XJ. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40(Web Server issue):W22–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  84. 84.

    Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, Zhao X, Liang C, Wang Y, Sun L, Shi M, Xu X, Shen F, Chen M, Han Z, Peng Z, Zhai Q, Chen J, Zhang Z, Yang R, Ye J, Guan Z, Yang H, Gui Y, Wang J, Cai Z, Zhang X. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27.3 associate with clear cell renal cell carcinoma. PLoS One. 2010;5(12):e15224.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Prof. Xian-Ming Chen from Creighton University and Prof. Xing-Quan Zhu from Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences for their assistance in reviewing and copy-editing the language of this article.

Funding

This work was supported by National Key Research and Development Program of China (Grant number 2017YFD0501305), National Natural Science Foundation of China (Grant number 31572509), Shaanxi Science and Technology Project (Grant number 2016NY-113), China Postdoctoral Science Foundation funded project (Grant number 2016 M592848), and Shaanxi Postdoctoral Science Foundation.

Availability of data and materials

We promise that all the data supporting the set of our conclusions are included in this article and its additional files. The complete RNA sequencing data generated and analyzed in this study have been submitted to the Genome Expression Omnibus (GEO) database, and the accession number is GSE108147 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108147).

Author information

Affiliations

  1. Department of Parasitology, College of Veterinary Medicine, Northwest A&F University, Yangling, 712100, China

    Guan-Jing Ren, Xian-Cheng Fan, Ting-Li Liu, Sha-Sha Wang & Guang-Hui Zhao

Authors
  1. Guan-Jing Ren
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xian-Cheng Fan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ting-Li Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Sha-Sha Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Guang-Hui Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

The experiments were conceived and designed by GHZ. GJR and XCF established the C. baileyi infection model. The samples were collected and prepared by GJR, XCF, TLL and SSW. GJR conducted the data analysis and qRT-PCR validation, and wrote the manuscript under the help and guidance of GHZ. All authors have read and approved the final manuscript.

Corresponding author

Correspondence to Guang-Hui Zhao.

Ethics declarations

Ethics approval

This study was carried out strictly in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the Ministry of Health, China. Our protocol with all animal experiments was approved by the Research Ethics Committee of Northwest A&F University and they agreed to anesthetise animals with the ether that was used for rats in previous studies. All efforts were made to minimize the pain of animals used in this study.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

The information of all novel lncRNAs identified in this study. (XLSX 106 kb)

Additional file 2:

All the differentially expressed mRNAs, lncRNAs and circRNAs in this study. (XLSX 150 kb)

Additional file 3:

LncRNA-protein coding gene pairs with co-location and co-expression relationships. (XLSX 13272 kb)

Additional file 4:

The results of circRNA-binding miRNAs analysis. (XLSX 33940 kb)

Additional file 5:

GO enrichment analysis of mRNAs and lncRNAs. (XLS 18518 kb)

Additional file 6:

KEGG pathway analysis of mRNAs, lncRNAs and circRNAs. (XLS 255 kb)

Additional file 7:

The oocyst excretion of chickens infected with C. baileyi. The horizontal axis represents the day post infection (dpi) with C. baileyi oocysts of chickens and the vertical axis shows the number of OPG. The oocysts were detected in faeces of infected chickens from 5 dpi, and the peaks of oocyst excretion were observed at 10 and 16 dpi. (TIF 6095 kb)

Additional file 8:

The parasites in histological observations by haemotoxylin & eosin staining (A and B) and by electron microscopy (C and D). No parasite was detected in the tracheas of the control chickens with haemotoxylin & eosinstaining (A, × 400) and electron microscopy (C, × 2000), but there were a lot of parasites in the tracheas of the experimental chickens (B and D). (TIF 4763 kb)

Additional file 9:

Primers used in qRT-PCR analysis. (XLSX 10 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ren, G., Fan, X., Liu, T. et al. Genome-wide analysis of differentially expressed profiles of mRNAs, lncRNAs and circRNAs during Cryptosporidium baileyi infection. BMC Genomics 19, 356 (2018). https://doi.org/10.1186/s12864-018-4754-2

Download citation

Keywords

  • Differentially expressed profiles
  • mRNAs
  • lncRNAs
  • circRNAs
  • Cryptosporidium baileyi

BMC Genomics

ISSN: 1471-2164

Contact us