Two highly inbred Fayoumi chicken lines (M5.1 and M15.2) obtained from Iowa State University (Ames, IA, USA) were used in the present study. The care and experimental use of these birds was approved by the Beltsville Area Institutional Animal Care and Use Committee of the US Department of Agriculture-Agricultural Research Service (USDA-ARS). The birds’ management, the experimental induction of NE, and sample collection were performed at USDA-ARS Beltsville, USA. Subsequent experiments, including high-throughput small RNA sequencing and quantitative real-time polymerase chain reaction analysis, were performed in the Republic of Korea, according to the USDA-ARS and Chung-Ang University agreement (59-1265-0-042F).

NE was experimentally induced in chickens from lines M5.1 and M15.2 by infecting them with E. maxima strain 41A (1.0×10⁴ oocysts/bird) by oral gavage on post-hatch day 14, followed by oral gavage with C. perfringens strain Del-1 (1.0×10⁹ colony-forming units/bird) on post-hatch day 18. The induction of NE was further promoted by feeding the chickens with an antibiotic-free organic starter diet containing 17% crude protein between days 1 and 18, followed by a standard grower diet containing 24% crude protein between days 18 and 20. The uninfected control chickens were maintained in separate cages and provided with standard feed and water ad libitum. To examine the establishment of NE, fecal samples were collected between post-hatch days 14 and 20 (0 to 7 days after E. maxima infection) and the oocyst numbers determined using a McMaster chamber (HK Inc., Tokyo, Japan). Gut lesion scores were determined at post-hatch day 20 (2 days after C. perfringens infection), on a scale of 0 to 4, by three independent observers in a blinded fashion. The average gut lesion scores of the NE-induced M5.1 and M15.2 chickens selected for RNA sequencing were about 1.9 and 2.3, respectively as described early [12].

The IML was collected from NE-induced and control chickens from M5.1 and M15.2 lines (five each) on post-hatch day 20. The mid-jejunum area of the small intestines in both NE-induced and control chickens was cut longitudinally and washed three times with ice-cold Hanks’ balanced salt solution supplemented with 100 U/mL penicillin and 100 mg/mL streptomycin (Sigma-Aldrich, St. Louis, MO, USA). The IML was carefully removed using a cell scraper (Nunc, Thermo Scientific Inc., Waltham, MA, USA). The total RNA was extracted from the samples using TRIzol (Invitrogen, Carlsbad, CA, USA), and purified using an RNeasy Mini Kit (Qiagen, Germantown, MD, USA). Trace DNA in the samples was removed using DNase (Promega, Madison, WI, USA) prior to high-throughput small non-coding RNA sequencing.

Approximately 5 μg of total RNA extracted from the IML of M5.1 and M15.2 NE-induced and control chickens was used by the Theragen Bio Institute (Suwon, Korea) to generate high-throughput small non-coding RNA data, using the HiSeq2000 high-throughput sequencer (Illumina, Inc., San Diego, CA, USA) according to the standard manufacturer protocols. Raw data were processed to obtain clean reads by removing low-quality reads, 5’primer and 3’primer contaminants, sequences without the insert tag, poly A tags, and reads that were shorter than 18 nucleotides. To further identify sequence reads with reliable chromosomal locations, all clean reads were mapped onto the reference chicken genome (Gallus gallus v.4) obtained from the UCSC Genome database (University of California, Santa Cruz, CA, USA). Among the clean reads, those with desired sequence lengths between 18 and 30 nucleotides were selected as potential small non-coding RNAs [10].

The annotation of small non-coding RNAs was performed following standard protocols from the Theragen Bio Institute (Korea). Specifically, miRNAs were identified by aligning the sequence reads with the precursor or mature miRNAs of chickens in miRBase (version 18.0) [15]. Additionally, the possibility of novel miRNAs, and their corresponding secondary structures, was predicted by MIREAP from unannotated reads that did not fall in any small non-coding RNA category. The miRNAs that were differentially expressed in the IML of NE-induced chickens, when compared to their respective controls, were determined based on their fold difference (≥2-fold change) and p-value (≤ 0.01). The miRNAs that were significantly upregulated or downregulated in the NE-induced chickens, when compared to their respective controls, were subjected to hierarchical clustering using the Multiexperiment Viewer (MeV) tool obtained from the TM4 microarray software suite [16]. Additionally, the miRNAs that were differentially expressed between the NE-induced chickens were determined and clustered.

The expression levels of four miRNAs that differentially expressed in the IML of NE-induced chickens when compared to their respective controls were examined by real-time qPCR. The selected miRNAs were: gga-miR-9-5p, upregulated in NE-induced M5.1 chickens; gga-miR-20b-5p, downregulated in NE-induced M5.1 chickens; gga-miR-196-5p, upregulated in NE-induced M15.2 chickens; and gga-let-7d, downregulated in NE-induced M15.2 chickens. The expression levels of the top-ranking/top-scoring target genes for each of these miRNAs (mitogen-activated protein kinase kinase kinase 3 [MAP3K3], mitogen-activated protein kinase kinase kinase 2 [MAP3K2], aquaporin 4 [AQP4], and insulin-like growth factor 2 mRNA binding protein 3 [IGF2BP3], respectively) were also examined by real-time qPCR. The predicted target genes for each miRNA were identified using the miRDB, an online database for miRNA target prediction and functional annotations [17].

Approximately 2 μg of total RNA from the IML of M5.1 and M15.2 control and NE-induced chickens was prepared using TRIzol and purified using an RNeasy Mini Kit. Traces of DNA in the samples were removed using DNase prior to cDNA synthesis for miRNA or target gene amplification. The cDNA for miRNA amplification was synthesized from the total RNA (2 μg) using the NCode miRNA First-Strand cDNA Synthesis and qRT-PCR Kit (Invitrogen, USA) according to the manufacturer’s protocols. Forward primers for gga-miR-9-5p (5′-TCT TTG GTT ATC TAG CTG TAT GA), gga-miR-20b-5p (5′-CAA AGT GCT CAT AGT GCA GGT AG), gga-miR-196-5p (5′-TAG GTA GTT TCA TGT TGT TGG), and gga-let-7d (5′-AGA GGT AGT GGG TTG CAT AGT) were designed following the guidelines of Agilent Technologies (Santa Clara, CA, USA) and synthesized by Genotech Co. Ltd. (Daejeon, Korea). Real-time qPCR was performed in triplicate with 20 μL of reaction mix containing 1 μL of cDNA (100 ng), 0.5 μL of miRNA-specific forward primer (5 pmol), 0.5 μL of NCode universal reverse primer (5 pmol), 10 μL of 2× AMPIGENE qPCR Green Mix Lo-ROX (Enzo Life Sciences Inc., Farmingdale, NY, USA), and 8 μL of nuclease free water. The thermal conditions for qPCR were as follows: pre-incubation at 50°C for 2 min and 95°C for 2 min; 45 cycles of 2-step amplification at 95°C for 15 s, and 60°C for 30 s; and a final cooling step at 37°C for 30 s. The samples were quantified using a Light Cycler 96 Real-time PCR System (Roche Diagnostics, Indianapolis, IN, USA). The expression of miRNAs in the NE-induced samples was calculated after normalizing with chicken U24 small nucleolar RNA (Z48762, Forward primer: 5′-GGG ATG TAA AAA AAT ACT TGC TAT C). The correlation coefficient of miRNAs expression by RNA sequencing and real-time qPCR was further analyzed.

The cDNA for target gene amplification was synthesized from the total RNA (2 μg) using a RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific Inc., USA) according to the manufacturer’s protocols. Forward (F) and reverse (R) primers for MAP3K3 (XM_015299375, F: 5′-TGA CCC GGA AAT ACA CTC GC; R: 5′-ACA TTC CCA GCA GAG TCA CG), MAP3K2 (XM_015290042, F: 5′-TGA GGG ACC CAC CAG AAA GA; R: 5′-GCG AGT CAC GTT CTC TGT CA), AQP4 (NM_001004765, F: 5′-GTG CCA GTA TGA ACC CTG CT; R: 5′-AAG AGC ACC AGC AAG GAC TG), and IGF2BP3 (NM_001006359, F: 5′-GCC TTG GCA GTT GGA GCT AT; R: 5′-AGC TTG GCA TCT GGT CCT TC) were designed using the Primer3 program (National Center for Biotechnology Information [NCBI] database) and synthesized by Genotech Co. Ltd. Real-time qPCR was performed in triplicate with 20 μL of reaction mix containing 1 μL of cDNA (100 ng), 0.5 μL of gene-specific forward primer (5 pmol), 0.5 μL of gene-specific reverse primer (5 pmol), 10 μL of 2× AMPIGENE qPCR Green Mix Lo-ROX, and 8 μL of nuclease free water. The thermal conditions for qPCR were as follows: pre-incubation at 95°C for 3 min; 45 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 30 s; and a final dissociation cycle at 95°C for 10 s, 65°C for 60 s, 97°C for 1 s, and 37°C for 30 s. The samples were quantified using a Light Cycler 96 Real-time PCR System (Roche Diagnostics, USA). The expression of genes in the NE-induced samples was calculated after normalizing with chicken glyceraldehyde 3-phosphate dehydrogenase (NM_204305, F: TGC TGC CCA GAA CAT CAT CC; R: ACG GCA GGT CAG GTC AAC AA).

The predicted target genes of the four miRNAs that differentially expressed in the IML of NE-induced chickens (gga-miR-9-5p, gga-miR-20b-5p, gga-miR-196-5p, and gga-let-7d) when compared to their respective controls were identified using the miRDB [17]. All target genes for each miRNA were then entered into the Kyoto encyclopedia of genes and genomes (KEGG) database [18] to search for their involvement in various pathways. Preferentially, those genes involved in signal transduction pathways and immune system-related pathways were selected. Next, the interactions of target genes that were involved in these pathways were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database program (Version 10.0) [19]. The results were retrieved based on a medium confidence view (score 0.4) and using the default prediction methods.

The expression of miRNAs and their target genes in control chickens and NE-induced chickens was analyzed using the Student’s t test in the SAS software (SAS Institute, Cary, NC, USA). The significance level of these results was analyzed using the general linear model in the SAS software. Statistical significance was indicated as follows: * p<0.05, ** p<0.01, and *** p<0.001.

Initially, high-throughput sequencing of small non-coding RNAs in the IML of control and NE-induced Fayoumi M5.1 and M15.2 chickens was performed using the HiSeq2000 high-throughput sequencer, and followed standard protocols to screen the raw sequencing data. Raw data were processed to obtain clean reads by removing low-quality reads, adaptor-ligated contaminants, sequences without the insert tag, and poly A tags. Sequencing reads shorter than 18 nucleotides were also excluded. After cleaning, there were 28,193,988 and 28,235,783 reads in the IML of control and NE-induced M5.1 chickens, respectively. There were 30,369,880 and 28,229,041 reads in the IML of control and NE-induced M15.2 chickens, respectively. These reads were mapped onto the reference chicken genome obtained from the UCSC genome database. The results revealed that 59.7% and 60.65% of the total reads were mapped for the control and NE-induced M5.1 chickens, respectively. Approximately 59.41% and 46.26% of the total reads were mapped for the control and NE-induced M15.2 chickens, respectively. Similarly, a significant proportion of unique reads in the IML of control and NE-induced Fayoumi chickens was obtained (Table 1). Among the total/unique reads, the most frequent read length was between 22 to 24 nucleotides in our test samples. The total reads returned by high-throughput small RNA sequencing include several categories of non-coding RNAs such as miRNA, transfer RNA (tRNA), ribosomal RNA (rRNA), small nucleolar RNA (snoRNA), small cytoplasmic RNA (scRNA), and small nuclear RNAs (snRNA) [10]. However, we were particularly interested in investigating miRNA distribution and expression in the IML of NE-induced Fayoumi chickens. miRNA is a strong candidate for the regulation of gene expression at the post-transcriptional level [20], and dysregulated miRNA expression has been reported in disease conditions [8]. Furthermore, the annotation of miRNAs is simple compared to the annotation of the other small non-coding RNAs that are involved in the post-transcriptional regulation of gene expression.

To identify the distribution of miRNAs across our test samples, the total reads obtained from the IML of control and NE-induced M5.1 and M15.2 chickens were aligned with precursor or mature chicken miRNAs using miRBase. According to miRBase source, 152 known miRNAs in control M5.1 chickens and 157 known miRNAs in NE-induced M5.1 chickens were found. In the M15.2 line, approximately 119 miRNAs and 144 miRNAs were found in control and NE-induced chickens, respectively (Table 2). These results indicate that there are slightly more miRNAs in the IML of M5.1 chickens, when compared to M15.2 chickens. Specifically, the number of known miRNAs was found to be higher in the NE-induced M5.1 chickens. There are clear instructions for identifying or naming miRNAs. Briefly, when two miRNAs originate from the same precursor, the predominantly expressed sequence is named as miRNA and the non-predominant sequence is named as miRNA* (from the opposite arm of the precursor). When the evidence is insufficient to determine which of the two miRNAs is predominant, they are named as miRNA from the 5’arm (miRNA-5p) and miRNA from the 3’arm (miRNA-3p) [21]. In this study, several sequencing reads from all test samples were found to be matched with known miRNAs, miRNA-5p, and miRNA-3p, but no reads were matched with miRNA*. The unannotated reads from all test samples were further examined using MIREAP to identify potentially novel miRNAs and their secondary structures. In total, 116,809 and 95,053 reads were found to be novel miRNA candidates in the control and NE-induced M5.1 chickens, respectively. Similarly, 112,573 and 151,737 reads were found to be novel miRNA candidates in the control and NE-induced M15.2 chickens, respectively.

Since several known miRNAs distributed across our test sam ples, we aimed to identify differentially expressed miRNAs in the IML of NE-induced M5.1 and M15.2 chickens, when compared to their respective controls. After setting the fold difference (≥2-fold change) and p-value (≤0.01), 50 upregulated miRNAs and 26 downregulated miRNAs in the NE-induced M5.1 chickens was found, when compared to the control (Supplementary Table S1). On the other hand, 32 upregulated miRNAs and 11 downregulated miRNAs in the NE-induced M15.2 chickens was found, when compared to the control (Supplementary Table S2). These results reveal that the total number of differentially expressed miRNAs/upregulated miRNAs is higher in the NE-induced M5.1 chickens than in the NE-induced M15.2 chickens. Figure 1 shows the hierarchical clustering of differentially expressed miRNAs in the IML of NE-induced M5.1 (Figure 1A) and M15.2 (Figure 1B) chickens. When we compared the treatment groups (NE-induced M5.1 chickens and NE-induced M15.2 chickens), 29 upregulated miRNAs in the NE-induced M15.2 chickens and 63 upregulated miRNAs in the NE-induced M5.1 chickens was found (Figure 1C, Supplementary Table S3). This result reinforce that the M5.1 chickens are more resistant to help upregulation of miRNAs that are crucial for the post-transcriptional regulation of gene expression. Upregulation of miRNAs potentially leads to downregulation of the genes that cause harmful effects to the cellular processes.

In high-throughput miRNA sequencing studies, identifying the target genes of miRNAs is helpful for their further analysis and annotation. The target genes of miRNAs can be identified by various computational methods primarily based on sequence complementarity search. In this study, the predicted target genes of miRNAs that were differentially expressed in the IML of NE-induced chickens (Supplementary Tables S1–S3) were identified using miRDB [17]. Next, four miRNAs that showed the highest fold-changes in the NE-induced chickens, compared to their respective controls, was selected for further analysis and annotation. The selected miRNAs were: gga-miR-9-5p (upregulated 37.6 fold in the IML of NE-induced M5.1 chickens), gga-miR-20b-5p (downregulated 35.27 fold in the IML of NE-induced M5.1 chickens), gga-miR-196-5p (upregulated 36.03 fold in the IML of NE-induced M15.2 chickens), and gga-let-7d (downregulated 36.8 fold in the IML of NE-induced M15.2 chickens). The miRNA gga-miR-9-5p targets 260 genes, with the top-ranking/top-scoring gene being MAP3K3 (target rank 1, score 100). The miRNA gga-miR-20b-5p targets 581 genes, with the top-ranking/top-scoring gene being MAP3K2 (target rank 1, score 100). The miRNA gga-miR-196-5p targets 92 genes, with the top-ranking/top-scoring gene being AQP4 (target rank 1, score 99). Finally, the miRNA gga-let-7d targets 207 genes, with the top-ranking/top-scoring gene being IGF2BP3 (target rank 1, score 100).

Real-time qPCR was used to further examine the expression of the four miRNAs (gga-miR-9-5p, gga-miR-20b-5p, gga-miR-196-5p, and gga-let-7d) that showed the highest fold-changes in the IML of NE-induced chickens, compared to their respective controls, as well as the expression of their top-ranking target genes (MAP3K3, MAP3K2, AQP4, and IGF2BP3, respectively) (Figure 2). The miRNAs gga-miR-9-5p and gga-miR-20b-5p and their target genes MAP3K3 and MAP3K2, respectively, were amplified with miRNA- or gene-specific cDNAs prepared from the IML of control and NE-induced M5.1 chickens. The miRNAs gga-miR-196-5p and gga-let-7d and their target genes AQP4 and IGF2BP3, respectively, were amplified with miRNA- or gene-specific cDNAs prepared from the IML of control and NE-induced M15.2 chickens. The expression of gga-miR-9-5p in the IML of NE-induced M5.1 chickens was not as expected based on the sequencing results. The expression of its target gene MAP3K3 showed significant increase in the NE-induced chickens compared to the control. The expression of gga-miR-20b-5p was significantly decreased, and the expression of its target gene MAP3K2 significantly increased, in the IML of NE-induced M5.1 chickens. The expression of gga-miR-196-5p was significantly increased, and the expression of its target gene AQP4 significantly decreased, in the IML of NE-induced M15.2 chickens. The expression of gga-let-7d was slightly decreased, and the expression of its target gene IGF2BP3 was increased, in the IML of NE-induced M15.2 chickens (Figure 2). The miRNAs expression trends determined by RNA sequencing and real-time qPCR were mostly consistent, with a correlation coefficient of R² = 0.82549 (Figure 3).

Several studies have reported that the overexpression of miRNAs downregulates their target genes, while downregulation of miRNAs upregulates their target genes, owing to their post-transcriptional gene regulatory functions [7]. This hypothesis is somewhat supported in the case of gga-miR-9-5p, gga-miR-20b-5p, gga-miR-196-5p, and gga-let-7d and the expression of their target genes MAP3K3, MAP3K2, AQP4, and IGF2BP3, respectively in the IML of NE-induced M5.1 chickens. In miRNA studies, it is important to know whether the expression and functions of selected miRNAs have been reported previously. Among the four miRNAs selected for this study, the expression and functions of miR-9-5p, miR-196-5p, and let-7d had not been previously reported in chickens. A recent study in chicken reported that miR-20b-5p involved in the proliferation and differentiation of chicken myoblast [22]. In humans and/or mice, the expression and functions of miR-9-5p, miR-20b-5p, miR-196-5p, and let-7d have been reported in various disease conditions. Among the four miRNA target genes, the expression and functions of MAP3K3, MAP3K2, and IGF2BP3 had not been previously reported in chickens. AQP4 is a member of the water channel family characterized by high water permeability, and its expression was detected in several chicken tissues including gastro-intestinal tract [23]. In humans and/or mice, MAP3K family members MAP3K3 and MAP3K2 were reported as the regulators of the NF-kB and MAPK pathways, and play an important role in enhancing the immune response in the microenvironment of the infected cells [24,25]. Since these two genes were found upregulated in the NE-induced M5.1 chickens, they might involve in enhancing the immune response in this resistant line. The immune related AQP4 autoantibodies positively involved in an inflammatory disease of the central nervous system in humans [26]. Since the AQP4 gene was found downregulated in the NE-induced M15.2 chickens, there might be a clue that this gene decreases the immune response in this susceptible line. IGF2BP3 was reported to decrease the activity of immune cells, and blocking the expression of this gene could increase the fighting ability of immune system [27]. Since IGF2BP3 was found upregulated in the NE-induced M15.2 chickens, it might decrease the immune response in this susceptible line.

The predicted target genes of the four miRNAs (gga-miR-9-5p, gga-miR-20b-5p, gga-miR-196-5p, and gga-let-7d) that showed the highest fold-changes between the IML of NE-induced chickens and their respective controls, in sequencing analyses, were first retrieved from the miRDB database. According to miRDB, gga-miR-9-5p targets 260 genes (target score: 100-50), gga-miR-20b-5p targets 581 genes (target score: 100-50), gga-miR-196-5p targets 92 genes (target score: 99-50), and gga-let-7d targets 207 genes (target score: 100-50). All target genes for each miRNA were then entered into the KEGG database [18] to search for their involvement in pathways identified in chickens. KEGG returned various chicken pathways in relation to these target genes for each miRNA. The intestinal mucosa is one of the most important tissues for the generation of the innate immune response against infectious pathogens that enter through the gastrointestinal tract [28]. Therefore, our analyses was restricted to find the miRNA target genes involved specifically in immune system-related pathways, and signal transduction pathways that essentially support the immune response [29,30]. In our results, the target genes of gga-miR-9-5p were found in eight signal transduction pathways and three immune system-related pathways (Figure 4A). In particular, maximum four genes were involved in the phosphatidylinositol signaling system when compared to the number of genes involved in other pathways. The full list of gga-miR-9-5p target genes that are involved in signal transduction and immune system-related pathways are shown in Supplementary Table S4. The target genes of gga-miR-20b-5p were found in 12 signal transduction pathways and four immune system-related pathways, with the maximum of 16 genes in MAPK signaling pathway (Figure 4B, Supplementary Table S5). The target genes of gga-miR-196-5p were found in eight signal transduction pathways and two immune system-related pathways, with the maximum of two genes in the calcium signaling pathway (Figure 4C, Supplementary Table S6). Finally, the target genes of gga-let-7d were found in nine signal transduction pathways, with the maximum of five genes in the MAPK signaling pathway (Figure 4D, Supplementary Table S7).

Next, the target genes of gga-miR-9-5p, gga-miR-20b-5p, gga-miR-196-5p, and gga-let-7d that were found to be involved in chicken signal transduction and immune system-related pathways were entered into the STRING database program [19] to find their interactions. The STRING database program enables the analysis of gene/protein-protein interactions in many species including chickens [19,30]. The gene interaction results based on a medium confidence view (score 0.4) and using the default prediction methods, which included neighborhood, gene fusion, co-occurrence, co-expression, experiments, databases, and textmining were retrieved. In total, 17 gene targets of gga-miR-9-5p were involved in signal transduction and immune system-related pathways; however, only seven genes showed interactions (Figure 5A). Among the eight target genes of gga-miR-196-5p that were involved in signal transduction and immune system-related pathways, none of the genes showed interactions with any of the others (Figure 5C). Among the 15 target genes of gga-let-7d that were involved in signal transduction and immune system-related pathways, only three genes showed interactions (Figure 5D). By contrast, of the 45 gga-miR-20b-5p target genes involved in signal transduction and immune system-related pathways, about 35 genes showed confident interactions (Figure 5B). In the STRING database program, interactions include physical and functional associations that are mainly derived from genomic context, high-throughput experiments, co-expression, and PubMed literature [19]. Overall, the miRNA gga-miR-20b-5p target genes were involved in more signal transduction and immune system-related pathways in the KEGG database and showed more confident interactions in STRING database. The principal reason underlying this finding could be that there were more than 500 target genes for gga-miR-20b-5p, compared with much fewer target genes for the other differentially expressed miRNAs in the IML of NE-induced chickens.