Impact of the Genetic Background on the Composition of the Chicken Plasma MiRNome in Response to a Stress

Circulating extra-cellular microRNAs (miRNAs) have emerged as promising minimally invasive markers in human medicine. We evaluated miRNAs isolated from total plasma as biomarker candidates of a response to an abiotic stress (feed deprivation) in a livestock species. Two chicken lines selected for high (R+) and low (R−) residual feed intake were chosen as an experimental model because of their extreme divergence in feed intake and energy metabolism. Adult R+ and R− cocks were sampled after 16 hours of feed deprivation and again four hours after re-feeding. More than 292 million sequence reads were generated by small RNA-seq of total plasma RNA. A total of 649 mature miRNAs were identified; after quality filtering, 148 miRNAs were retained for further analyses. We identified 23 and 19 differentially abundant miRNAs between feeding conditions and between lines respectively, with only two miRNAs identified in both comparisons. We validated a panel of six differentially abundant miRNAs by RT-qPCR on a larger number of plasma samples and checked their response to feed deprivation in liver. Finally, we evaluated the conservation and tissue distribution of differentially abundant miRNAs in plasma across a variety of red jungle fowl tissues. We show that the chicken plasma miRNome reacts promptly to the alteration of the animal physiological condition driven by a feed deprivation stress. The plasma content of stress-responsive miRNAs is strongly influenced by the genetic background, with differences reflecting the phenotypic divergence acquired through long-term selection, as evidenced by the profiles of conserved miRNAs with a regulatory role in energy metabolism (gga-miR-204, gga-miR-let-7f-5p and gga-miR-122-5p). These results reinforce the emerging view in human medicine that even small genetic differences can have a considerable impact on the resolution of biomarker studies, and provide support for the emerging interest in miRNAs as potential novel and minimally invasive biomarkers for livestock species.


Introduction
MicroRNAs (miRNAs) are small endogenous RNAs that pair to sites in mRNAs to direct post-transcriptional repression [1]. Recent work indicates that cells release miRNAs in the extra-cellular environment, predominantly in association with either vesicles or protein complexes that protect them from RNases [2]. These miRNAs can be passively released as a result of tissue damage or actively released from healthy cells, from which they may subsequently reach the bloodstream and constitute what it is now referred to as the ''blood-circulating extra-cellular miRNome''. Because extra-cellular miRNAs can be easily extracted from any body fluid and profiled through microarrays, real time quantitative PCR or sequencing, blood-circulating miRNAs are currently regarded as being among the most promising clinical biomarkers for the diagnosis, prognosis, and therapeutic treatment of a variety of pathological conditions including cancer, cardiovascular diseases, diabetes, liver pathologies, and sepsis [3][4][5].
Minimally invasive biomarkers which can be profiled by tiny amounts of body fluids are important for animal breeding applications. Livestock species are often subjected to a variety of stress conditions, and extra-cellular miRNAs could be used in tandem with other phenotypic measurements to monitor the responses of individual animals or populations [6]. For example, in human medicine an increase in the abundance of 'tissue specific' or 'organ specific' miRNAs in blood plasma (or other body fluids) could serve as an indication of toxicity or injury in a particular tissue/organ. Moreover, extra-cellular miRNAs could serve as specific markers for the diagnosis of diseases caused by viruses able to encode miRNAs from their genome (like several herpesviruses), as viral miRNA should be preserved in the extracellular space after the infected host cells die [7].
Intense artificial selection for phenotypic traits of economic importance has produced a large variety of livestock breeds and populations worldwide, and several experimental populations have been raised for research purposes. Among these, two chicken lines have been divergently selected since 1975 for high (R+) or low (R2) residual feed intake (RFI) at constant egg production and body weight, starting from a common unselected population of Rhode Island Red layers produced from six sires and fifty dams [8]. Since then, the two lines have been maintained as closed populations by within-line mating (nine sires and 45 dams per line, with one year generation interval), and the between-line difference in RFI is currently equivalent to five phenotypic standard deviations [9].
The R+ and R2 lines differ in only a small proportion of their genome. This has been earlier reported by DNA fingerprinting analyses and interpreted as the combined effect of divergent selection and genetic drift [10]. This appears to be confirmed by preliminary whole-genome sequencing analyses (from pools of seven individuals per line), which identified roughly 850,000 SNPs segregating between the two lines, of which 15,000 are differentially fixed (Lagarrigue et al., unpublished data). Conversely, strong differences between the two lines are found at the phenotypic level. R+ chickens are characterized by higher feed intake and lower adiposity compared to R2, as well as by increased thermogenesis and reduced liver lipogenesis [11][12][13]. An additional striking difference is the excessive appetite in the R+ and the reduced appetite in the R2, which corresponds to an 89% increase in feed intake in the R+ compared to the R2.
Here we evaluated extra-cellular miRNAs circulating in plasma as potential biomarker candidates of a response to a feed deprivation stress in a livestock species. We chose the R+ and R2 chicken lines challenged for feed deprivation as an experimental model because of their extreme energy metabolism. Our results indicate that the plasma miRNome of the R+ and R2 lines reacts promptly to a feed deprivation stress. This response reflects the phenotypic divergence that these chicken lines acquired through long term artificial selection, implying that even low levels of genetic variation can affect the miRNome considerably. This study is, to our knowledge, the first attempt to characterize the plasma miRNome of a bird species. Samples obtained for the validation of the sampling procedure: blood samples were obtained from two adult cocks (one R+ and one R2) to extract total RNA from each different blood component (plasma, white cells and red cells) and to check the levels of a subset of miRNAs in each blood component by RT-qPCR (see below).

Animals and sampling procedures
Samples obtained for deep sequencing and RT-qPCR of miRNAs: nine cocks from the R+ line and nine from the R2 line (from five sire families for each line), were blood sampled at 38 weeks of age after 16 hours of feed deprivation (denoted as FD). Five cocks per line were sacrificed before re-feeding to obtain liver samples in the FD condition. The rest of the cocks were blood sampled again four hours after re-feeding (denoted as RF, which represents the control condition) and sacrificed to obtain liver samples in the RF condition. At each blood sampling 2 ml of blood were taken by syringe from the wing vein and gently transferred into EDTA blood tubes. These samples were used for RT-qPCR validations (see below).
In order to obtain sufficient yields of total RNA for the preparation of small RNA libraries, five supplementary cocks per line (belonging to the same sire families as the first sampled group) were blood sampled in the FD and RF conditions. For this group RNAs were extracted from 2 ml of plasma. Six of these samples (three FD and three RF, issued from three cocks per line) were chosen for deep sequencing; the rest of the samples were used for RT-qPCR.
RNA extractions: the blood samples were centrifuged after sampling at 4˚C for 10 minutes, at 500 x g; the supernatants were transferred to fresh tubes, centrifuged at 2,000 g x for 30 minutes at 4˚C and re-transferred into a new tube with an equal volume of PBS before filtration with a 0.22 mm filter. All RNAs from plasma were extracted using TRIzol LS Reagent according to the manufacturer's instructions (Life Technologies). Frozen liver samples (about 100 mg each) were homogenized in TRIzol Reagent (Life Technologies) using an Ultra-Turrax (IKA-Werk) and extracted according to the manufacturer's instructions. RNA yield was assessed using an Agilent 2100 Bioanalyzer and RNA 6000 pico kits. The analyses of small RNA profiles were carried out using the Agilent 2100 expert software (Rev. B.02.08.SI648). RNA samples were stored at 280˚C until required for analysis.

Small RNA library preparation and Sequencing protocol
Small RNA libraries were constructed starting from 25-50 ng of total plasma RNA using the TruSeq SBS Kit v5-GA, while two additional libraries were constructed using the TruSeq SBS v3-HS kit, both from Illumina, according to the manufacturer's instructions. The quantities of RNA used for small RNA libraries preparation were further checked with a Qubit Fluorometer (Life Technologies). Libraries were sequenced (single read) on a GA-IIx or on a HiSeq1000 Illumina sequencer. The raw reads have been deposited at the European Nucleotide Archive (ENA) with accession number PRJEB6619.

Post-sequencing analysis
First, a set of unique chicken miRNA precursors was built from the miRBase (version 19) [14] and Ensembl (version 72) databases. While miRBase is composed of experimentally identified miRNAs, Ensembl also includes precursor predictions based on stem-loop structure and sequence homology. The secondary structures of these precursors were then computed using the RNAfold tool from the Vienna RNA package [15]. To assign read count to mature miRNAs, all miRBase mature miRNAs were mapped on the set of precursors without allowing any mismatches. When miRNAs from multiple species, including chicken, mapped at the same position on a precursor, the chicken miRNA annotation was retained. Using this approach we annotated 985 mature miRNAs (including 791 known chicken miRNAs from miRBase and 194 putative orthologs) based upon 1080 non-redundant precursors. For the final set of 41 differentially abundant miRNAs (Table 1), the mature miRNA names were updated using the current miRBase release (version 20).
Next, cutadapt v1.2.1 [16] was used to perform 39 adaptor sequence removal and trim low-quality ends of reads. Reads between 19 nt and 24 nt in length were then collapsed to remove redundancy using an in-house python script. Bowtie v0.12.8 [17] was used to map collapsed reads to the set of chicken precursors, allowing at maximum one mismatch. To discard hairpins with a read distribution inconsistent with Drosha and Dicer processing sites (i.e., reads tilled across the precursor), we expected at least a 3:1 ratio between reads that matched on any of the stem-loop arms and reads located in the loop. For the remaining hairpins, reads that mapped inside the loop were not considered for quantification. Putative new mature miRNAs were named based upon the name of the hairpin on which they were located when no known mature miRNA matched the same precursor, or from the name of the miRNA mapping on the opposite strand of the precursor. The suffixes ''-5p'' and ''-3p'' were added if the new miRNA mapped on the 59 arm or the 39 arm, respectively.
Sequences that did not map to any precursor were then successively re-aligned with the same tool and criteria against the Gallus gallus genome release 4.0, the RFAM database depleted from the miRNAs, and chicken cDNAs retrieved from the Ensembl database.

Differential expression (abundance) analyses
The term ''differential abundance'' was used when referring to miRNAs circulating in plasma. Differential analysis of the small RNA-seq expression data was done using the R/Bioconductor packages edgeR (version 3.2.3) [18] and HTSFilter (version 1.0.0) [19]. Briefly, a negative binomial generalized linear model (GLM) was fit for each miRNA, including a sample-specific normalization factor to adjust for systematic differences among samples due to variable sequencing depth and RNA composition, a fixed line effect (R2 and R+), and a fixed feeding effect (FD and RF). As recommended in [19], after estimating pergene parameters for the full dataset, miRNAs with weak abundance levels across all samples were filtered using HTSFilter, leaving 148 for further analysis. Differential abundance was assessed between lines and feeding groups using a GLM likelihood ratio test, where P-values were adjusted for multiple testing using the method by [20] to control the false discovery rate (FDR) at 5%. Log-fold changes of miRNAs identified as differentially abundant in at least one comparison were subsequently visualized using heatmaps. Hierarchical clustering using the Euclidean distance and complete-linkage were used to identify groups of miRNAs based on their log-fold changes from the comparisons previously described in the differential analyses. We note that the interaction contrast between the line and feeding groups was found to be statistically insignificant for all miRNAs. The analysis of the qPCR data was performed using the R/Bioconductor package limma (version 3.16.5) [21]. For plasma samples, the qPCR abundance data were obtained for each animal from the two feeding conditions (FD, RF), while qPCR data for the liver samples were obtained from different animals in each feeding condition because of the necessity to euthanize the animals for tissue sampling. Using the log 2 -transformed abundance relative to the threshold cycle (CT) value (2 2DC T ), a linear model was fit for each miRNA. For the plasma data the linear model included fixed effects for line (R2 and R+), feeding group (FD and RF), and batch (as plasma samples were collected in two batches), as well as a random effect for each animal to account for the correlation between measurements made on each individual before and after re-feeding. For the liver data the animal random effect and sampling batch effect were omitted because not necessary. Differential abundance was assessed between lines and feeding groups using a moderated t-test, and P-values were adjusted to control the false discovery rate at 5% [20]. As before, the interaction contrast between the line and feeding groups was found to be statistically insignificant for all miRNAs in both the plasma and liver.

cDNA synthesis and qPCR quantification
Reverse transcription was done using miRNA-specific stem loop RT primers [22] and the TaqMan microRNA Reverse Transcription Kit according to the manufacturer's instructions (Applied Biosystem). In each reaction 10 ng of total RNA from liver or 600 pg of RNA extracted from plasma were mixed with 50 nM specific stem-loop RT primer. RT reactions were carried out at 16˚C for 30 min, 42˚C for 30 min and 85˚C for 5 min. Real-time quantitative PCR was done using standard TaqMan PCR protocols on an Applied Biosystems 7900HT System. PCR primers and probes were optimized to improved stability and mismatch discrimination using locked nucleic acid nucleotides [23] (Table S1).

Co-abundance analysis of miRNAs
As for the differential analysis (see above), we adopted the term of ''coabundance'' when referring to plasma miRNAs. Groups of co-abundant miRNAs were identified from the miRNA-seq data using the clustering approach implemented in the R package HTSCluster (version 1.2) [24]. Briefly, we assumed that the population of miRNAs arises from several distinct subpopulations or clusters, each of which can be modeled separately. The filtered population of miRNAs used for the differential analysis were thus modeled using a Poisson mixture model, where parameter estimation was performed using an Expectation-Maximization algorithm and the appropriate number of clusters present in the data (K55) was identified using the Integrated Completed Likelihood (ICL) model selection criterion. The cluster-specific parameters may be interpreted as the proportion of total reads attributed to each line (R2 and R+) and feeding (FD and RF) combination, which may be useful in identifying global trends in the data.

Comparison to miRNAs identified in red jungle fowl tissues
We made use of the miRNAs identified by the Chickspress database (available at http://geneatlas.arl.arizona.edu; SRA Bioproject Accession PRJNA204941) to evaluate the relative tissue abundance of differentially expressed miRNAs across a variety of red jungle fowl tissues. The Chickspress miRNA database comprises 659 million single-end, 50 bp reads mapped to the galGal4 version of the chicken genome. miRNA data were collected from two, 20-month-old red jungle fowl (one male, one female) across a variety of tissues: male adrenal gland, male and female adipose tissue, cerebellum, cerebrum, gonad, hypothalamus, heart, kidney, liver, lung, breast-muscle, sciatic nerve, proventriculus and spleen. The database contains mappings to known miRNAs (miRBase v20; [14]) as well as novel miRNAs identified using the miRTrap software [25] for each tissue. A FPKM (fragments per kilobase per million mapped reads) quantitative value is provided for each miRNA feature in each tissue.

Target gene predictions and functional analysis
An in silico search for putative binding sites of differentially abundant miRNAs was performed using TargetScan 6.0 [26]; note that non-conserved miRNAs (ENSGALT00000042483-3p, ENSGALT00000043002-3p, ENSGALT00000043002-5p, gga-miR-1736-3p) and those with borderline differential abundance (adjusted P-values of 0.0498; ccr-miR-133a-5p, mmu-miR-144-5p, gga-miR-20a, aca-miR-499-3p) were not included in the functional analysis. Perl scripts and UTR sequences from 23-species alignments were downloaded from the TargetScan website (http://www.targetscan.org). Predictions were ranked according to the context+ score, which models the contribution of different context parameters on site efficiency including site-type, 39-pairing, local AU, position, target site abundance, and seed-pairing stability [27]. Only the putative miRNA binding sites belonging to the upper quartile of ranked predictions and present in at least five species including Gallus gallus were retained.
The potential target genes for the differentially abundant miRNAs by comparison were then classified according to their biological function determined by the Gene Ontology (GO) System (http://www.geneontology.org). The enrichment analysis to identify over-represented GO categories was performed using the DAVID Functional Annotation Tool (version 6.7; http://david.abcc. ncifcrf.gov). Enrichment was assessed using a modified Fisher's Exact Test (EASE Score) and P-values were adjusted for multiple testing to control the FDR at 5% [28]. The DAVID clustering feature was used to identify relationships among enriched terms and to cluster them into functional groups.

Validation of the sampling procedure
In humans, blood cells are a major contributor to circulating miRNA, and factors such as hemolysis can alter plasma miRNA biomarker levels by up to 50-fold [29]. These effects might be exacerbated in birds due to the presence of nucleated erythrocytes. The profiling of small RNAs from chicken blood cell components indicated that small RNAs between 45-55 nt are particularly abundant in red cells followed by a similar peak in PBMCs, and in both these compartments RNAs of larger sizes are detected ( Figure 1). In plasma the only RNAs detected are 20-30 nt in length (corresponding to circulating miRNAs), thus suggesting there is no evidence of larger RNA species contaminating our plasma preparations ( Figure 1).

Identification of miRNAs circulating in chicken plasma
Individual libraries of small RNAs were produced from 3 R+ and 3 R2 animals, in both the food-deprivation (FD) and re-fed (RF) conditions (R+ RF; R+ FD; R2 RF and R2 FD). More than 292 million sequence reads were generated by small RNAseq. We focused on reads between 19 nt and 24 nt in size which cover the range of sequence lengths for miRNAs. Approximately 45% of the reads were discarded from subsequent analyses, mostly because of their short size (94% of the discarded reads were shorter than 19 nt). The remaining reads were successively aligned to chicken miRNA precursors, the chicken genome, other non-coding RNAs and chicken cDNAs. As expected, most sequences (79%) mapped to miRNA precursors ( Figure 2). We did not find any difference in the percentage of small reads corresponding to miRNAs between the different conditions. Finally, 649 mature miRNAs were identified (Table S2). Among these miRNAs, 410 have already been described in chicken, 98 showed a perfect sequence homology with a mature miRNA from another species, 94 were putative new mature miRNAs mapping on the opposite strand of a known chicken miRNA, and 47 miRNAs were Ensembl annotated novel miRNA precursors. The most abundant miRNA was miR-2188-5p, which represented 19% of total reads, followed by gga-miR-10b (gga-miR-10b -5p in miRBase v20), gga-miR-10a-5p, and gga-miR-146c-5p, which together represented about 36% of the total reads (13.3%, 13.9% and 9.2% respectively).

Effects of genetic divergence and feed deprivation on the plasma miRNAome
Differential abundance analysis The differential abundance analysis of miRNAs circulating in plasma was conducted by comparing the two chicken lines (R+ versus R2) to observe  differences between lines (''Line''), and by comparing feed deprived and re-fed birds (FD versus RF) to assess the effect of feed deprivation vs. the re-fed control (''Condition'').
We identified 23 and 19 miRNAs with significant differential abundance in the Condition and Line comparisons, respectively, with only two miRNAs (gga-let-7f and gga-miR-203) found in common between comparisons. The miRNAs representing the Condition contrast were on average more abundantly expressed than those representing the Line contrast ( Figure S1).
In addition, we examined the effect of feed deprivation vs. re-fed conditions within each chicken line (''Condition within R+'' and ''Condition within R2''), which resulted in zero and 19 differentially abundant miRNAs, respectively; the differences between lines within each feeding condition (''R+ versus R2 within RF'' and ''R+ versus R2 within FD'') resulted in two and eleven differentially expressed miRNAs, respectively.
The list of differentially abundant miRNAs is reported in Table 1.

Clustering analysis of miRNAs
First, a model-based analysis was carried out on the 148 miRNAs retained after filtering [19] to identify co-abundant groups. Based on the Integrated Completed Likelihood criterion the model with five clusters was retained ( Figure 3). The largest cluster (Cluster 3, containing 87 miRNAs) is characterized by underabundance in R2 animals in the feed deprived condition with respect to the other groups. Cluster 1 (containing six miRNAs) and Cluster 5 (containing 12 miRNAs) are largely characterized by over-abundance in the feed deprived condition, particularly for the R2 line. Finally, Clusters 2 and 4 (containing 18 and 25 miRNAs, respectively) represent clusters with more balanced abundance among groups. The miRNAs identified as differentially abundant between lines were assigned to Cluster 3 (14 of 19 differentially abundant miRNAs), while those exhibiting differential abundance between feeding conditions were divided among Clusters 1, 3, and 5 (6, 11, and 4 of 23 differentially abundant miRNAs, respectively). Cluster 1 was entirely composed of miRNAs (6 of 6) found to be differentially abundant between feeding conditions (Table 1). In addition, we examined a heatmap of the estimated log-fold changes for each contrast for all miRNAs identified as differentially abundant in at least one comparison, where hierarchical clustering was applied to both the miRNAs and the comparisons. This analysis highlighted groups of miRNAs which are either largely over-abundant or under-abundant in the R+ line as compared to the R2 line (''R+ vs. R2'' column) and groups of miRNAs which are either moderately over-abundant or under-abundant in the RF versus FD groups (''RF vs. FD'' column) (Figure 4). The ''R+ vs. R2 (RF)'' and ''R+ vs. R2 (FD)'' comparisons showed that the miRNA plasma content in the re-fed condition is rather similar between the two chicken lines, while larger differences are observed when animals are feed deprived. The comparison ''RF vs. FD (R+)'' and ''RF vs. FD (R2)'' indicated that the response to feed deprivation is stronger for the R2 line (no miRNAs reached significance in the response to feed deprivation of the R+ line; Table 1).

Plasma
We validated the sequencing results on a larger number of individual birds by RT-qPCR of six miRNAs found to be differentially abundant in plasma. We used plasma samples obtained from eight R+ and seven R2 birds, each one sampled in both the FD and RF condition. The miRNAs were selected among those found differentially abundant in the Condition comparison (gga-miR-204, gga-miR-2188-5p and gga-miR-365-3p), in the Line comparison (gga-miR-2188-3p, and gga-miR-122-5p) or in both (gga-let-7f-5p). These miRNAs were found to be present in plasma at very different levels of abundance, from an average of 397 normalized counts (gga-miR-122-5p) to an average of 2.1 million normalized counts (gga-miR-2188-5p).
The RT-qPCR confirmed the relative levels of abundance identified by deep sequencing for all six miRNAs. However, only two (gga-miR-204 and gga-miR-2188-5p) of the four miRNAs that were previously identified as differentially abundant between feeding conditions (Condition) were found to be significant ( Figure 5). Conversely, high significance values were found for all miRNAs (gga-miR-122-5p, gga-miR-2188-3p and gga-let-7f-5p) previously identified as differentially abundant between R+ and R-animals (Line). Remarkably, gga-miR-122-5p was confirmed to be significantly more abundant in the R+ animals than in the R2 ones despite the high variability observed in the two feeding conditions. Furthermore, gga-miR-204, which was expected to be significant only between feeding conditions, also exhibited significant differential abundance for the Line contrast ( Figure 5).

Liver
Since the R+ and R2 lines differ significantly for metabolic traits such as body fat content and liver lipid metabolism [11,12], we examined miRNA expression in liver from feed deprived and re-fed animals.
The analysis showed an opposite pattern compared to that observed in plasma. No significant differences were observed in the Line comparison, while all miRNAs were found to be highly significant in the Condition comparison ( Figure 6). With the exception of gga-miR-122-5p, which was significantly downregulated, all miRNAs were strongly up-regulated in response to food Challenging R+ and R2 Chicken Lines for Feed Deprivation deprivation. This pattern suggests that liver miRNAs contribute to the extensive gene downregulation pattern that has previously been observed in the liver of chickens in response to feed deprivation [30].
Interestingly, a significant difference between R+ and R2 lines (p50.03) was observed for gga-miR-122-5p in the feed deprived condition alone and not in the re-fed one (p50.82), with gga-miR-122-5p more expressed in the R2 line.

Conservation and patterns of tissue distribution of differentially abundant miRNAs in plasma
Thirty-five of the 41 mature miRNAs identified as differentially abundant in one or more of the Condition and Line comparisons (Table 1)
Three of the six remaining mature miRNAs (mmu-miR-143-3p, mmu-miR-145a-5p and dre-miR-210-5p) have no sequence homology with Gallus gallus annotated miRNAs, and three (ENSGALT00000043002-5p and -3p, and ENSGALT00000042483-3p) are encoded by two miRNA genes predicted in the chicken genome (http://www.ensembl.org). ENSGALT00000043002 maps to a non coding region downstream from 5.8 S rRNA, while ENSGALT00000042483 is contained within an intron of a lipoxygenase homology domains 1 gene (LOXHD1), and is transcribed in the same orientation. ENSGALT0000004248323p is the most highly significant miRNA in most differential expression analyses because it was undetected in R+ chickens ( Figure 4).

Comparison to miRNAs identified in red jungle fowl tissues
The expression of conserved miRNAs across tissues has not been extensively verified in chicken, with only partial data reported so far [31][32][33][34][35]. We examined the 35 miRNAs annotated in Gallus gallus for their expression levels in 14 tissues of jungle fowl ( Table 2). Most of these miRNAs are highly conserved, with conserved patterns of tissue distribution or tissue-specific enrichment. For example, the three members of the let-7 family (let-7a, let-7f, let-7k) are broadly expressed across tissues [36] and tissue enrichment has been found for miR-499-5p and 23p in heart [37], miR-122-5p in liver [38], miR-202-5p in testis [39] and gga-miR-107-3p in brain tissues [40] ( Table 2).
The mir-2188 gene is absent in the mammalian lineage. Since its discovery in fish [41] it has been detected in reptiles, amphibians and birds. Both gga-miR-2188-5p and miR-2188-3p are detected in all adult tissues, with the exception of the breast muscle. They are particularly abundant in the lung, and gga-miR-2188-5p is largely predominant across all tissues ( Table 2). Gga-mir-2954 and gga-mir-1736 are also absent in the mammalian lineage. The gga-mir-2954 family is avian specific [42]. Gga-miR-2954 has 100% homology with tgu-miR-2954-3p (zebra   finch) and is highly expressed in brain (cerebrum, cerebellum), liver, heart and breast muscle, while it is undetected in other tissues. Gga-mir-1736 has been described only in the chicken [35,43]. It is ubiquitous in jungle fowl (Table 2), and its low levels of expression across tissues are in agreement with previous findings [35,43]. This pattern is consistent with the expectation that less evolutionarily conserved miRNA genes are generally expressed at lower levels compared with broadly conserved miRNAs (reviewed by [44]). Interestingly, gga-miR-1736-3p is found in considerable abundance in the plasma of R2 animals (approximately 3000 sequence reads), and is six fold less abundant in R+ chickens, independently of the feed deprived or re-fed condition. This miRNA is contained in intron 3 of CARS2 (a mitochondrial aminoacyl-tRNA synthetase encoded by the nuclear genome) and is transcribed in the same orientation.

miRNA expression in blood cells
The data of the Atlas support the current view that most organs and tissues release miRNAs in the bloodstream; however, blood cells can be major contributors of the extracellular miRNA content in plasma [2]. Indeed, 19 of the miRNAs found to be differentially abundant in chicken plasma have been reported to be expressed during the different phases of hematopoiesis [45] and T cell development [46] in mammals ( Table 2).
As part of our validation of the procedure to obtain chicken plasma miRNAs (  Figure 1), we profiled four miRNAs (gga-let-7f-5p, gga-miR-365-3p, gga-miR-2188-5p and gga-miR-2188-3p) by RT-qPCR from whole chicken blood, PBMCs, plasma and red cells. In mouse, the mature miRNAs encoded by mir-let-7f are predominantly expressed in hematopoietic cells, while the homolog of gga-miR-365-3p is ubiquitous and has been detected at low levels in most blood cells [45]. Interestingly, in chicken all four miRNAs are predominantly expressed by red cells and gga-miR-365-3p is undetected in PBMCs ( Figure S2). These data are in agreement with the average number of sequence reads found in the plasma of R+ and R2 animals for gga-miR-2188-5p (average 2.1 million reads), gga-let-7f -5p (average 143000 reads), gga-miR-365-3p (average 4600 reads) and gga-miR-2188-3p (average 1700 reads). Given that red cells in avian blood are on average 150 times more abundant than PBMCs [47], it can be expected that red cells contribute significantly more than PBMCs to the circulating miRNome of plasma in birds, as our data suggest.

Target gene predictions and functional analysis
We performed TargetScan predictions to identify the potential gene targets of the conserved differentially abundant miRNAs in the Condition and Line comparisons. As expected, a large number of potential target genes was found (2261 genes, of which 979 were targeted by multiple miRNAs). A great variability in the number of target genes per miRNA was observed, ranging from six for gga-miR-184-3p to 276 for hsa-let-7a-3p.
In the Condition and Line comparisons 145 and 52 significantly overrepresented GO terms were identified, respectively, of which 32 were in common between the two comparisons. Among the most enriched biological functions in both comparisons we found the regulation of DNA and RNA metabolic processes (Table S3, Table S4). To identify functional clusters of overrepresented GO terms, we used the clustering algorithm in DAVID to group similar, redundant and heterogeneous annotation terms. As expected, the top ranked annotation groups for both comparisons included GO terms related to transcriptional regulation and RNA metabolic processes. Moreover, in the Condition comparison, functional clusters with high enrichment scores grouped GO terms related to gene expression and macromolecule biosynthesis and cell morphogenesis (Table S5). In the Line comparison, the top ranked functional clusters included biological processes such as cell motion, vasculature development and epidermis development (Table S6).

Discussion
The pool of miRNAs circulating in chicken plasma reacts to 16 hours of feed deprivation, with effects that are quickly recovered (4 hours) after re-feeding. Feed deprivation leads to varied changes in abundance of the 148 plasma miRNAs retained after filtering, which grouped into five clusters of co-abundance ( Figure 3). However, only the six miRNAs (miR-2954, gga-miR-2188-5p, gga-miR-365-3p, gga-miR-193b-3p, gga-miR-204 and mmu-miR-145a-5p) which compose Cluster 1 are found to be almost three-fold more abundant when concurrently considering R2 and R+ animals (Table S2).
Our further results of differential abundance analysis indicate that i) the miRNA response to feed deprivation is indeed strongly influenced by the different genetic backgrounds and that ii) the miRNA divergence between lines is stronger under feed deprived conditions. The net result is that only a limited overlap is found between the sets of differentially abundant miRNAs in the two main comparisons (Table 1). Similarly, a limited overlap was observed in the enrichment analysis of potential target genes, where a broad variety of functional categories was present for both comparisons. This large heterogeneity in the enrichment analyses is expected for miRNAs circulating in plasma. As confirmed by our data (Table 2), plasma miRNAs can originate from several body tissues. This implies that putative targets will necessarily represent several broad GO classes, whose terms cannot be matched a priori to the physiology of a specific tissue. An interesting future line of research would be to cross data on the plasma miRNome with a transcriptome analysis from a large variety of tissues from the same animals/conditions to obtain more precise information on the functional impact of circulating miRNAs.
The response to feed deprivation is mostly driven by the response of the R2 animals, which appear to be highly reactive to this stress (Table 1, Figure 4). This is reflected by the relatively low significance obtained in RT-qPCR validations in plasma when jointly considering R+ and R2 animals ( Figure 5). Conversely, the miRNAs that differentiate lines are generally less abundant in the R2 animals and this signature is validated even when considering both feeding conditions (  Figure 4, Figure 5). However, the difference between the R+ and R2 chickens is mostly driven by feed deprivation. Indeed, when animals are fed, only two chicken-specific miRNAs (gga-miR-1736-3p and ENSGALT00000042483-3p) differentiate the two lines, while larger differences are observed when animals are feed deprived, with most miRNAs decreasing in abundance in the R2 line (  Table 1, Figure 4). From a physiological point of view, similar contrasted responses under fed and feed deprived conditions have been described for these two lines for several plasmatic metabolites [13]. Significant differences were observed only in the feed deprived state for plasmatic concentrations of glucose, non-esterified fatty acids, uric acid, T4 thyroxine hormone and the T3:T4 ratio, implying differences in the control of energy expenditure and the endocrine balance between the two lines, exacerbated in the feed deprived condition [13,48].
These results can be interpreted in light of the documented physiological divergence that these two lines acquired upon long term divergent selection. The higher feed intake and lower adiposity in the R+ chickens compared to the R2 has been justified as a possible alteration of the glucose-insulin axis [13]. This hypothesis is supported by the different plasmatic contents found in the R+ and R2 lines of three highly conserved miRNAs with a key regulatory role in energetic metabolism (gga-miR-204, gga-miR-let-7f-5p and gga-miR-122-5p) ( Table 2). These three miRNAs are more abundant in the plasma of R+ animals compared to R2 ( Figure 5) and all respond to feed deprivation in the liver ( Figure 6).
Gga-miR-204 has been shown in mammals to be transcribed in the pancreatic beta-cells to block insulin production by down regulating MAFA, an insulin transcription factor [49]. In chicken, the plasma levels of gga-miR-204 increase significantly under the feed deprived condition in both chicken lines. This miRNA is detected at low levels in liver ( Table 2) where it is strongly upregulated in response to feed deprivation ( Figure 6). Together with other members of the let-7 family, miR-let-7f regulates the glucose metabolism in multiple organs [50] and has an important role in the control of fasting glucose concentration [51]. The levels of gga-miR-let-7f-5p in plasma are five-fold decreased in R-animals after feed deprivation, while in the R+ no significant differences in abundance levels were observed (Table S2, Table 1). However, at least in liver, this miRNA is strongly upregulated by feed deprivation (Figure 6). It is remarkable that the three other mature miRNAs encoded by the let-7 family (gga-miR-let-7a-5p and -3p, and gga-miR-let-7k-5p) were found to be more abundant in the R+ than in the R2 line ( Figure 4, Table S1). Given that in mice the overexpression of the let-7 family leads to decreased fat mass and body weight [50], our data suggest a fundamental role of the let-7 family in response to intense selection for metabolic traits in these lines. Finally, an interesting example is provided by gga-miR-122-5p. This miRNA is known to be implicated in cholesterol biosynthetic pathway and fatty acid metabolism and makes up 70% of all the liver miRNAs [52]. When miR-122-5p is downregulated the hepatic synthesis rate of fatty acids and cholesterol decrease [52,53], which is in agreement with the finding that this miRNA is downregulated in response to feed deprivation in the liver ( Figure 6). However, in plasma, the level of gga-miR-122-5p decreases after feed deprivation only in the R2 line. This apparent paradox may be explained by a higher constitutive expression in the R+ line of gga-miR-122-5p by other tissues, with subsequent release in plasma. In particular, in mammals the expression of this miRNA in the adipose tissue has been found to be about 200-fold less than in liver [54] and this is fully confirmed by the jungle fowl Atlas (Table 2). However, this interpretation would contradict the very low amounts of abdominal fat in the R+ animals [11]. Further investigation on other tissues is required to explain these observations, also considering that this miRNA has recently been found to be down-regulated in adipose tissue during diet-induced development of obesity in mice [55].
Our results indicate that even small differences in the genetic background can have a considerable impact on the resolution of biomarker studies. To date most studies of circulating miRNAs have been focused on their clinical relevance as non invasive biomarkers for early diagnosis of disease and monitoring of treatment. To this end, panels of biomarkers need to be ''robust'' towards sources of variation due to genetic background, sample processing, age of patients, or other causes [5]. Interestingly, recent studies to define extra-cellular miRNA panels for breast cancer detection revealed differences in miRNA expression between different ethnic groups, with little overlap between Caucasian and African women [56]. Genetic background has been indicated as one of the possible reasons for the ostensible lack of reproducibility in published data regarding circulating miRNAs as markers for breast cancer [57,58]. Furthermore, it has been reported that a choline-and folate-deficient diet causing nonalcoholic fatty liver disease (NAFLD) determined a different extent of modulation of some miRNAs, including miR-122, in both liver and plasma of divergent strains of mice [59]. These changes in circulating miRNAs based upon genetic variation and diet corroborates the observations in our chicken model.
Overall, this study gives support to the emerging interest on miRNAs as potential novel and non invasive biomarkers for livestock species [60], particularly for the study of phenotypes that reflect the robustness of animals towards environmental challenges in addition to disease. The dynamics of the plasma miRNome upon feed deprivation and recovery reflect the fast kinetics and controlled reversibility of the multi-level transcriptional regulation of stress response, adding evidence to the key role of miRNAs in these processes [61,62]. It has recently been reported that an Argonaute 2 switch regulates circulating miR-210 to coordinate hypoxic adaptation across anatomically distinct cells [63], adding to the theory that some miRNAs may be specifically secreted in the blood stream and function in intercellular communication between distant tissues [2]. Finally, the observation that the genetic background is an important factor in shaping the extra-cellular miRNome has important consequences for potential field applications. Because the R+/R2 lines are an example of extreme phenotypic divergence, it could be expected that the sets of stress-responsive miRNA in plasma will show more overlap between breeds or populations exposed to a stress. However, identifying common panel biomarkers will be challenging in outbred and admixed populations, due to the heterogeneity of genetic backgrounds and to the presence of several sources of environmental variation to which miRNAs are highly responsive concomitantly with other genome regulators [64].