All of the animal experiments were conducted in accordance with the Guidelines for Experimental Animals established by the Ministry of Science and Technology (Beijing, China). Animal experiments were approved by the Animal Management Committee (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS, Beijing, China). Ethical approval on animal survival was given by the animal ethics committee of IAS-CAAS (approval number: IASCAAS-AE20140615).

Specific-pathogen-free White Leghorn chickens were supplied by the Beijing Laboratory Animal Research Center (BLARC, Beijing, China) and were treated as described in previous studies (Li et al., 2010; Gou et al., 2012). In brief, the SPF chickens were raised in climate-controlled, fully enclosed isolation facilities at the experimental center of China Agriculture University (Beijing, China) under identical management conditions. At 3 d of age, a total of 150 SPF chickens were orally challenged with 1 ml PBS containing 10⁸ cfu of S. enteritidis (50041) and another 75 birds received 1 ml PBS as controls. Blood samples from each of 30 challenged and 15 control chickens were taken at 24 h post infection and birds were killed and the spleens were dissected, snap frozen and held at −80°C. Bacterial burden (expressed as cfu/10 μL blood) was determined indirectly by serovar-specific quantitative real-time PCR (qPCR), and along with clinical severity, was used to evaluate the resistance/susceptibility to SE challenge, as described in previous studies (Deng et al., 2008; Gou et al., 2012). 10 μL EDTA-anticoagulated blood was used for DNA extraction using MiniBEST Whole Blood Genomic DNA Extraction Kit (Takara, Code No. 9781) according to the manufacturer's instructions. Amplification was carried out in a total 25 μL reaction mixture, containing 0.6 μL of each primer (10 μM), 0.75 μL of dNTPs (10 mM), 1.25 U of ExTaq DNA Polymerase (Takara), 5.5 μL of 5 × PCR buffer (Mg²⁺), 0.8 μL of TaqMan probe (5 μM), and 2 μL of template, with deionized water to 25 μL. Each PCR consisted of a 5 min hot start at 95°C followed by 40 cycles of 30 s at 94°C, 30 s at 55°C, and a fluorescence read step. The probe (5′-FAM-TGCAGCGAGCATGTTCTGGAAAGC-TAMRA-3′) and primers set (forward primer, 5′-TCCCTGAATCTGAGAAAGAAAAACTC-3′; reverse primer, 5′-TTGATGTGGTTGGTTCGTCACT-3′) were designed from the SdfI gene (Gen-Bank Accession No. AF370707.1), as described in Gou et al. (2012). The qPCR assay was calibrated by relating threshold cycle (Ct) values to cfu, as determined by enumeration after plating serial dilutions of S. enteritidis and standard culture.

In this study, the clinical symptoms (diarrhea, drooping wings, and dying) and bacterial load of SE at 24 h after challenge were used together to discriminate susceptible (S, SE-challenged, slight clinical symptoms and >10⁷ cfu/10 μL blood) from resistant (R, SE-challenged, severe clinical symptoms and <10⁵ cfu/10 μL blood) birds. No SE was detected in the Controls (C). Total splenic RNA was extracted from three birds in each of the three groups, S, R, and C, using miRNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. RNA was quantified using the NanoDrop ND-2000 spectrophotometer (NanoDrop Products, Wilmington, DE) and purity was assessed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA) with RNA Integrity Number (RIN) number >7.0. Total RNA was stored at −80°C until used.

Total RNA of each sample (~1 μg) was used to prepare the miRNA sequencing library, which included the following steps: (1) 3′-adapter ligation with T4 RNA ligase 2 (truncated); (2) 5′-adapter ligation with T4 RNA ligase; (3) cDNA synthesis with an RT primer; (4) PCR amplification; and (5) extraction and purification of 120–140 bp PCR amplified fragments (corresponding to ~15–25 nt small RNAs) from polyacrylamide gels. An Agilent 2100 Bioanalyzer quantified the libraries, after which the samples were diluted to a final concentration of 8 pM and cluster generation was performed on the Illumina using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA), following the manufacturer's instructions. The 9 miRNA libraries were constructed and single-end sequenced (36 bp) on an Illumina Hiseq 2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol. The raw data of each sample was not <10 M reads. The raw reads were subjected to the Illumina Pipeline filter (Solexa v0.3), and then the dataset was further processed with ACGT101-miRv4.2 (LC Sciences, Houston, TX) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. Subsequently, the 18–25 nt length unique sequences were BLASTed to chicken precursors in miRBase 20.0 (Kozomara and Griffiths-Jones, 2014) (http://www.mirbase.org/) to detect known miRNAs and novel 3p- and 5p- derived miRNAs. One mismatch inside the sequence and length variation at both 3′ and 5′ ends were allowed in the alignments. The unique sequences were mapped to chicken mature miRNAs in hairpin arms recognized as known miRNAs, and mapped to the other arm of known chicken precursor hairpins opposite the annotated mature miRNA-containing arm considered to be novel 5p- or 3p-derived miRNAs. The remaining sequences were mapped to other selected species in miRBase 20.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the chicken genomes to identify their genomic positions. The aforementioned miRNAs were considered to be known miRNAs. To identify the novel predicted miRNAs, the unmapped sequences were BLASTed against the chicken genome database, and the hairpin RNA structures comprising sequences were identified using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). Modified reads per million (RPM) reads was used to quantify the normalized reads, the formula was: Normalized Expression (NE) = Actual miRNA count/Total count of clean reads. MicroRNAs were regarded as being differentially expressed (DE) based on normalized deep-sequencing levels (with the exclusion of 3 RPM) in S, R and C groups, respectively. The DE miRNAs based on normalized counts were analyzed using Student t-tests (Huang et al., 2015; Li et al., 2016) according to the experimental design and the significance threshold was set as P < 0.05. The normalized read counts of some miRNAs were set to be 0.01 for further calculation if they had no reads in the library.

Nine cDNA libraries were also constructed from splenic RNA (1 μg) of these same birds and sequenced by LC-BIO (Hangzhou, China) on an Illumina HiSeq 2500 platform and 125 bp paired end reads were generated. The raw reads were first processed through FastQC to obtain the clean data, by removing the reads that contain sequencing adapter contaminations or poly-N and the low quality reads, Q-values for which were <20. Some potential residual ribosome RNA data were also removed from the remaining data by alignment. Clean reads were then mapped to the Gallus gallus database using TopHat (Trapnell et al., 2009), and the mapped reads were assembled de novo using Cufflinks (Trapnell et al., 2010). Expression levels of mRNAs were quantified as fragments per kilobase of exon per million mapped reads (FPKM) using the Cufflinks package (Trapnell et al., 2010). Analysis of DE genes between the three groups of chickens was performed using the Cuffdiff with a P < 0.05 and |log₂fold change| > 0.58.

Only target DE genes that were predicted by both TargetScan 6.2 and miRanda 3.3 for all of the DE miRNAs were considered further. Gene Ontology (GO) and KEGG pathway enrichment of target DE genes were analyzed by DAVID 6.8 (http://david.abcc.ncifcrf.gov/), which is based upon a Fisher Exact statistic methodology similar to that previously described (Huang et al., 2009). GO and KEGG results were filtered using P < 0.05.

In order to build the miRNA-mRNA interaction network, the following method was used, as described in previous studies (Ye et al., 2016): A target gene was identified by the direction of change in a pairwise comparison, for example S to C, being the reverse of changes in the miRNAs. The miRNA-mRNAs interaction network was constructed using Cytoscape v2.8.3 software (http://www.cytoscape.org/).

The pmiR-RB-Report™ (RiboBio, Guangzhou, China) including double luciferase reporter genes was used to test and validate the target sites for gga-miR-155 and gga-miR-101-3p. The 3′ UTR of IRF4 and LRRC59 containing gga-miR-101-3p and gga-miR-155 binding sites were amplified from chicken genomic DNA. The primers for PCR are provided, as follows: IRF4: GGCGGCTCGAGGATCCTCAGAATAAGTGTT (forward) and AATGCGGCCGCGTTAGAAG-TCCCTAGAAAA (reverse); and LRRC59: GGCGGCTCGAGATGCTACAGCAGAACTCGC (forward) and AATGCGGCCGCCAGACAAATTGATGCGAAA (reverse). All PCR products were cloned into the pmiR-Repor Vector using Xhol and NotI restriction enzymes. Luciferase reporter experiments were performed in 293T (human embryonic kidney) cells, obtained from ATCC. Cells were seeded in 96-well plates at a density of 5 × 10⁴ cells/well and cultured under routine conditions with 10% fetal bovine serum. When the cells reached 70 to 80% confluence, pmiR-3′ UTR (100 ng) was co-transfected with 50 nM of a negative control or a gga-miR-101-3p mimic (GenePharma, Shanghai, China) using 0.30 μL of FugeneHD (Promega, Madison, WI) according to the manufacturer's instructions. The relative luciferase activity was measured 48 h after transfection by the Dual-Glo Luciferase Assay System (Promega).

To further validate the biological function of gga-miR-155 and gga-miR-101-3p in a chicken macrophage-like line HD11, 100 μM mimic (gga-miR-155), inhibitor (gga-miR-101-3p) and control oligos (gga-miR-NC) were transfected into HD11 cells using 12-well plates and TransIT®-2020 (Mirus Bio, Madison, WI) per the manufacturer's instructions. HD11 cells were grown at 37°C with 5% CO₂ in RPMI-1640 medium that contained 10 mM HEPES, 1 mM sodium pyruvate, 1% glutamine, 1% MEM NEAA, 10% fetal bovine serum, and 5% chicken serum (all reagents from Gibco). After 36-h transfection, the cells were harvested using MiniBEST Universal RNA (Takara, Code No. 9767) to extract the total RNA. For the LPS exposure, macrophages were challenged with 1 mg/ml LPS and harvested at different times for RNA extraction. Cells with no stimulation were collected as the control, and each experiment had three biological replicates.

To validate and characterize the DE miRNA and DE transcripts identified via high-throughput sequencing, qPCR analyses were performed in an ABI 7500 Detection System (Applied Biosystems, Foster, CA). The miScript SYBR Green PCR kit (Qiagen, Valencia, CA) and PCR Master Mix (SYBR Green) Kit (Toyobo, Osaka, Japan) were used in qPCR to determine the abundance of mRNAs and miRNAs, using β-actin and U6 genes as reference genes, respectively. The relative mRNA and miRNA expression level was calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). All primers are described in Supplementary Tables S8, S9. Three independent replications were used for each assay and data are presented as means ± SD.

An average of 6,348,747 high quality clean reads per miRNA sample, which represented 359,232 unique reads in the range of 18–26 nt in the nine libraries were obtained from splenic samples via next generation sequencing (Supplementary Table S1). These high-quality reads were mapped to chicken precursors in miRBase to identify known and novel miRNAs for further analysis. Low levels of large fragments, such as mRNA and rRNA, were also found, which indicated the high-quality and minimal degradation of RNA samples in the present study. For all nine samples, the distribution of the small RNA sequence length was mainly concentrated at 22 nt, followed by 23 and 21 nt (Figure 1A), which is consistent with the typical size range for Dicer-derived products and in agreement with most of the previous reports from other animal species.

A total of 2238 miRNAs, classified into five categories (Supplementary Table S2), were detected via BLAST in miRBase. After removing the less expressed miRNAs, i.e., the expression levels were <3 after the normalization of dataset (in at most 3 samples), 744 miRNAs were identified including 439 known chicken miRNAs, and 62 potentially novel miRNAs (defined as PC-3p or PC-5p) in chicken spleen after oral challenge with SE (Supplementary Table S2).

A total of 32 miRNAs exhibited significantly different expression (DE) among the C, R and S groups. The results showed that, for S vs. C 16 DE miRNAs (7 up- and 9 down-regulated); for R vs. C 13 DE miRNAs were found (4 up- and 9 down-regulated) and 13 were found in the R vs. S comparison (10 up- and 3 down-regulated; Figure 1B, Table 1 and Supplementary Table S3). To validate the expression profiles from sequencing, 6 miRNAs were also examined by q-PCR (Figure 1C). Except for gga-miR-92-5p with a slight difference in the R group, the expression patterns of gga-miR-101-3p, gga-miR-126-3p, gga-miR-155, gga-miR-103-5p, and gga-miR-455 were comparable by both methods. The expression profiles from the deep sequencing were therefore considered as being reliable and appropriate for further analysis.

The differences in splenic expression between the controls, resistant and susceptible birds were examined. Four miRNAs were significantly differently expressed in both S vs. C and R vs. C, and 5 in both S vs. C and S vs. R, as well as 3 in both R vs. C and R vs. S (Figure 1B, Table 1). Only 1 miRNA (gga-mir-1677) was significantly differently expressed in all three groups of birds (Figure 1B). Several miRNAs previously reported to be involved in immune responses such as miR-155, miR-9, miR-30, miR-126, and miR-29 families were identified. Also identified here were several new candidate miRNAs associated with SE infection, such as gga-miR-29c-5p (up-regulated, P = 0.01) and gga-miR-137-3p (down-regulated, P = 0.009).

In order to validate the roles of DE miRNAs in affecting expression of their target genes, mRNA in the same samples was also profiled by sequencing. Based on both TargetScan and miRanda systems, a total of 273 DE genes can be targeted by the 32 DE miRNAs in the three groups (Supplementary Table S4). As shown in Figure 2A, 148 (S vs. C), 40 (R vs. C), and 85 (R vs. S) DE miRNA-targeted genes were identified with fold change (FC) > 1.50 or < 0.67 and P < 0.05 (Supplementary Table S5). The heat map and hierarchical clustering demonstrated distinct profiles of the unique miRNA-targeted genes in the three comparisons (Figures 2B–D). Several immune-related genes were found to be significantly DE in spleen after challenge with SE. For example, the expression of IL8, CXCR4, and IRF4 were significantly up-regulated following SE challenge (FC 5.21, 3.69, and 2.02, respectively). To validate the expression profiles from sequencing, transcript abundances of eight genes were measured by qPCR (Figure 2E); overall, there was good concordance between the two methods.

The ultimate function of miRNAs is at the level of the activity of target genes. In this study, functional annotation and pathway enrichment analysis of 273 target DE genes were performed using Gene Ontology (GO) and KEGG. Potential function analysis of these genes showed that 2 immune-related KEGG pathways and 1 biological process were significantly enriched (P < 0.05), including Apoptosis, NOD-like receptor signaling pathway, and adaptive immune response (GO:0002250) (Table 2). The present results suggest that the changed miRNAs may regulate these immune-related targets in chicken spleen during SE infection.

Potential functional analyses for host immune responses to SE infection between R and S chickens were further characterized, based on the target genes of significant DE miRNAs between these two groups and the controls. For S vs. C, 4 pathways were enriched (P < 0.05), viz. Apoptosis, Spliceosome, mTOR signaling pathway, Insulin signaling and Jak-STAT signaling pathway; 2 biological processes were significantly enriched (P < 0.05); regulation of inflammatory response and heart looping. In the R vs. C comparisons, NOD-like receptor signaling pathway was significantly enriched (P < 0.05); defense response to bacterium (GO:0042742), immune-related biological processes, was enriched but not significantly (P = 0.06). Except for regulation of proteasomal protein catabolic process (GO:0061136) being enriched (P < 0.05), no pathways were found when R was compared with S (Table 3). These results are consistent with the susceptible birds being more likely to exhibit apoptosis due to an inflammatory response, while the resistant birds showed more of an innate immune response to SE infection.

Most descriptions of miRNA function have focused on their roles as post-transcriptional regulators of target mRNAs. Based on the putative miRNA-mRNA regulatory pairs, it was found that 91 SE-related genes can be targeted by 29 of the 32 DE miRNAs (Supplementary Table S6). The potentially important interaction networks for immune-related miRNA-mRNA pairs are shown in Figure 3. The relative expression of innate/inflammatory marker genes such as PIK3CD was significantly up-regulated following SE infection. Some mRNAs are highly connected and regulated by multiple miRNAs. For example, CXCR4 is involved in cytokine-cytokine receptor interaction and was identified as a potential target of gga-miR-155 and gga-miR-9-3p. IFR4 was predicted to be regulated by gga-miR-30d and gga-miR-101-3p. LRRC59 was predicted as a potential target of gga-miR-103-5p and gga-miR-155. One Salmonella-regulated miRNA of particular interest identified through the present study is gga-miR-101-3p. Although, the expression levels of miR-101-3p were relatively moderate, it is highly connected (>8 SE-related target genes) within the miRNA-mRNA network. These have not been previously reported to be associated with Salmonella infection, and are predicted here to regulate several immune-related genes.

The luciferase reporter gene system was used to validate the above-stated predicted interactions. The 3′ UTRs of IRF4 and LRRC59 were cloned into luciferase reporter plasmids to test gga-miR-101-3p and gga-miR-155 functions in vitro. Transfection with a gga-miR-101-3p mimic resulted in significant (P < 0.01) reduction in relative luciferase activity for IRF4 plasmids (Figure 4), compared with negative control miRNA (random miRNA sequence) and a no-insert control. Similarly, transfections with mimics resulted in significant (P < 0.05) reduction in relative luciferase activity for LRRC59 (Figure 5) compared with the negative miRNA and no-insert controls. These results indicate that similar responses are likely to be happening in the host during SE infection, that is, the down-regulation of gga-miR-101-3p may result in increased expression of IRF4 during Salmonella infection, and up-regulation of gga-miR-155 may inhibit expression of LRRC59.

After 36-h treatment with mimic, elevating gga-miR-155 significantly repressed the mRNA expression levels of LRRC59 compared to the miR-NC and negative controls (P < 0.05); In contrast, after 36-h treatment with gga-miR-101-3p inhibitor, the mRNA expression levels of IRF4 were significantly increased (P < 0.05) compared to the controls (Figure 6). In order to address the effect of miR-155 and miR-101 on the induction of pro-inflammatory cytokines in response to LPS, the expression levels of IL-6 and TNF-α were measured in a macrophage inflammatory response model. The results showed that miR-155 overexpression markedly decreased the expression of IL-6 and TNF-α compared with control miRNA or miR-155 inhibitor (Figure 7A; P < 0.01), while miR-101 knockdown significantly decreased the expression of IL-6 and TNF-α compared with control miRNA inhibitor (Figure 7B; P < 0.05). These data demonstrate that gga-miR-155 and gga-miR-101 could regulate the production of pro-inflammatory cytokines, IL-6 and TNF-a, which may play a negative role in response to LPS stimulation in chickens.

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.