Identification and Profiling of MicroRNAs in the Embryonic Breast Muscle of Pekin Duck

MicroRNAs (miRNAs) regulate gene expression by fully or partially binding to complementary sequences and play important roles in skeletal muscle development. However, the roles of miRNAs in embryonic breast muscle of duck are unclear. In this study, we analyzed the miRNAs profiling in embryonic breast muscle of Pekin duck at E13 (the 13th day of hatching), E19, and E27 by high-throughput sequencing. A total of 382 miRNAs including 359 preciously identified miRNAs 23 novel miRNA candidates were obtained. The nucleotide bias analysis of identified miRNAs showed that the miRNAs in Pekin duck was high conserved. The expression of identified miRNAs were significantly different between E13 and E19 as well as between E27 and E19. Fifteen identified miRNAs validated using stem-loop qRT-PCR can be divided into three groups: those with peak expression at E19, those with minimal expression at E19, and those with continuous increase from E11 to E27. Considering that E19 is the fastest growth stage of embryonic Pekin duck breast muscle, these three groups of miRNAs might be the potential promoters, the potential inhibitors, and the potential sustainer for breast muscle growth. Among the 23 novel miRNAs, novel-miRNA-8 and novel-miRNA-14 had maximal expression at some stages. The stem-loop qRT-PCR analysis of the two novel miRNAs and their two targets (MAP2K1 and PPARα) showed that the expression of novel-mir-8 and PPARα reached the lowest points at E19, while that of novel-mir-14 and MAP2K1 peaked at E19, suggesting novel-miRNA-8 and novel-miRNA-14 may be a potential inhibitor and a potential promoter for embryonic breast muscle development of duck. In summary, these results not only provided an overall insight into the miRNAs landscape in embryonic breast muscle of duck, but also a basis for the further investigation of the miRNAs roles in duck skeletal muscle development.


Introduction
In birds, a primary and a secondary generation of fibers arise respectively during the embryonic and fetal stages of development. Following these two waves of myogenesis, the total number of fibers is fixed [1] and there are no significant changes in fiber numbers during later bird development [2,3]. Birds, such as poultry, are therefore attractive models for studying muscle development during the prenatal development and has become a research focus [4]. The family of myogenic regulatory factors (MRFs), which includes MyoD, Myf5, MRF4, and MyoG, is important for embryonic muscle development [5]. The members of MRFs coordinate the expression of genes involved in muscle growth, morphogenesis, muscle cell differentiation and contractility. Recently, it has been shown that miRNAs play important roles in skeletal muscle development [6,7].
MicroRNAs (miRNAs) are short (approximately 22 nucleotides) noncoding RNA molecules that bind to complementary mRNAs sequences, hereby promoting mRNA degradation or translational repression [8][9][10]. An essential role of miRNAs in skeletal muscle development is evidenced since the deletion of Dicer which is responsible for the maturation of miRNAs results in perinatal lethality due to skeletal muscle hypoplasia [11]. In particular, the critical roles of three muscle-specific miRNAs, miR-1, miR-133 and miR-206, in the regulation of myogenesis have been well documented [6,12] with miR-1 and miR-133 regulating different aspects of skeletal muscle development both in vitro and in vivo [13]. The miR-1 promotes muscle cell differentiation by repressing the expression of histone deacetylase 4 (HDAC4), an inhibitor of muscle differentiation. In C2C12 myoblasts, miR-133a promotes proliferation by partially repressing serum response factor (SRF). Like miR-1, miR-206 promotes differentiation of C2C12 myoblasts in vitro by repressing the expression of the DNA polymerase A subunit (Polal) [14], connexin 43 (Cx43) [15], follistatin-like 1 (Fstl1) and utrophin (Utrn) [16]. In addition, other miRNAs have also been shown to play a role in muscle development. Over expression of miR-181 during muscle cell differentiation is important in promoting myogenesis by down-regulating the homeobox protein Hox-A11, an inhibitor of myogenesis [17]. The miR-486 has been shown to induce myoblast differentiation by down-regulating Pax7 [18], while miR-27b regulates Pax3 translation and ensures myogenic differentiation [19]. Recently, studies have shown that miR-148a positively regulates myogenic differentiation via down-regulating Rho-associated coiled-coil containing protein kinase 1 (ROCK1), a known inhibitor of myogenesis and miR-214 may target the negative regulators of Myf5, MyoD and myogenin in the corresponding stages of skeletal muscle development in vivo to regulate embryonic myogenesis [5]. It has recently been evidenced that miRNAs is one of the most abundant players of gene regulatory molecules in vertebrates.
Currently, there are approximate 21264 predicted hairpin miRNAs and 25141 novel mature miRNAs from 193 species in the publicly available miRNA database miRBase (Release 19.0, August 2012) (http://www.mirbase.org). It is surprising that there is no duck miRNAs presented in the miRBase because duck not only has tremendous agricultural importance [20][21][22] but also is a natural reservoir of influenza A viruses [23,24]. A few studies have begun to explore duck miRNAs in various aspects. Zhang et al. profiled miRNAs in duck feather follicle and skin with highthroughput sequencing technology [25]. Powder et al. identified and compared the miRNAs expressed in cranial NC cells from three avian species (chicken, duck, and quail) before and after species-specific facial distinctions occur [26]. In addition, the novel and differentially expressed miRNAs in the ovaries of laying and non-laying duck have been identified by Yu [27]. However, the analysis of miRNAs in many other tissues, including skeletal muscle, is still deficient. Therefore, exploring of miRNAs in duck skeletal muscle will greatly improve the understanding of the role of miRNAs in avian skeletal development.
In a preliminary study we observed that the growth rate of embryonic breast muscle of Pekin duck reached its peak at the 19 th day of hatching (E19) and the expression level of muscle regulatory factor 4 (MRF4), coincidentally, peaked at E19 (unpublished data). These results indicated that E19 is the fastest growth stage of embryonic breast muscle of Pekin duck, yet the underlying molecular mechanism regulating this rapid growth stage is still unclear. Given the important roles of miRNAs in skeletal muscle development, identification of the differentially expressed miRNAs in different developmental stages is a critical first step to investigate the function of miRNAs in embryonic muscle development in ducks.
Here, we analyzed miRNA expression from embryonic breast muscle of Pekin duck at developing stage of E13 (the 13 th day of hatching), E19, and E27 by high throughput sequencing. With bioinformatics analysis and stem-loop qRT-PCR validation of some identified and novel miRNAs, we identified differentially expressed miRNAs in embryonic breast muscles in duck and hereby providing a basis for further investigation on the molecular mechanisms of breast muscle development in duck.

Results and Discussion
High-throughput sequencing and small RNA discovery Total RNAs from embryonic breast muscle of Pekin duck at stage E13, E19 and E27 were used to construct small RNA libraries for high throughput sequencing. A total of 14881453, 13411560 and 15775148 raw reads were obtained from the E13, E19 and E27 libraries, respectively. After quality control and adaptor removal, 14580115, 13016970 and 15549081 high-quality reads were available from E13, E19 and E27 libraries, respectively (Table 1) for further analysis. Length distribution analysis showed that most reads ranged from 21-23 nt with the percentage of the 22 nt reads of the total reads being 41.63%, 65.34% and 64.92% for the three libraries, respectively (Fig. 1). The high-quality reads were subsequently annotated to different classes of RNA categories (identified miRNAs, repeats-associated RNA, rRNA, tRNA, snRNA, snoRNA, etc) using different databases such as miRBase (V19.0), from Rfam(10.1) and Genbank ( Table 2, Fig. 2). The most abundant RNA species (based on read count) in the three libraries was classified as miRNAs, accounting for 48.41%, 81.70% and 76.34% in the three libraries, respectively. This indicates that the deep sequencing data were highly enriched for mature miRNA sequences and that the data are well suitable for expression profiling analysis of identified miRNAs and discovering of novel miRNAs. The second most abundant category was rRNAs, accounting for 26.75%, 9.48% and 12.26% in the three libraries, respectively. In addition, unknown RNAs also represent-   (Table S1).

Identification and evaluation of identified miRNAs
Currently, there is no miRNAs database of duck in the miRBase. We therefore aligned our small RNA reads against the miRNA precursor/mature miRNA of all animals in the miRBase (V19.0) and extracted the sequence and count of miRNA families (no specific species) which were presented in our samples. The criterions for identifying identified miRNAs in the breast muscle tissues of Pekin duck were as follows: Firstly, to account for the difference among species, we aligned the reads to the precursor/ mature miRNAs of all animals in miRBase allowing two mismatches and free gaps. Secondly, the highest expressed miRNA for each mature miRNA family was used to construct a miRNA database which was used to align the reads to estimate the expressions of the miRNAs by counting the number of reads aligning to each miRNA in the database. Thirdly, we predicted the precursor of the identified miRNAs and discarded the miRNAs without a predicted hairpin structure.
Among the identified miRNAs, we first investigated their base preference. It has been suggested that the bases of 1 st and 9 th position from 59 terminal which are responsible for targeting mRNAs for gene regulation, and the 39 terminal positions are enriched with a U base [29]. In the current study, we obtained similar results from E13, E19 and E27 stages. From the three stages, U accounted for 85.9%, 93.9% and 88.7%, respectively at the 1 st position, 48.2%, 38.4% and 50.0%, respectively for the 9st position and 74.0%, 86.2% and 87.6% respectively, for the 59 terminal (Fig. 3, Table S2). The bases preference of duck miRNAs are consistent with the statistical results through previous largescale genome analysis performed by Zhang [30] and indicated that miRNAs in duck, like in other species, is highly conserved.
A total of 359 mature miRNAs, belonging to 262 families, were identified from the E13, E19 and E27 libraries (Table S3). The largest miRNA family identified was Let-7 consisting 21 members.
Many other families, such as miR-1, miR-10, miR-1013, miR-124, miR-128, miR-130, miR-133, miR-206 and miR-148 contained two members, whereas most identified families had only one member. The finding that many members of identified miRNA families in all other animals were expressed in the breast muscle of Pekin duck suggested that a broad scale miRNAs are needed for muscle development of duck. Different family members also displayed drastically different expression levels. For example, the abundance of the Let-7 family varied from 3 read (let-7a-2-3p) to 2149081 reads (let-7f-5p). This was also the case for some other miRNA families, such as miR-499 (from 1to 10316 reads), miR-199 (from 1004 to 256996 reads) and miR-10 (from 17 to 6058 reads). The existences of a dominant member in a miRNA family suggest that the regulatory role of certain family is performed by the dominant member at a specific stage.
After literature mining, 24 of the identified miRNAs were involved in skeletal muscle development (Table 3). Almost all of these muscle-specific miRNAs (myomiRs) were highly expressed. The most abundant one was miR-1a-3p (Table S3), which was represented by 1906033, 5142560, and 2847317 reads in E13, E19 and E27 libraries, respectively (Table 3) which is in line with the distinct role of miR-1 in modulating skeletal muscle cell proliferation and differentiation [13]. miR-206 and miR-206-3p (4 th and 5 th abundance, Table S3) accounted for more than 760,000 reads in each library, which is consistent with their established function in skeletal muscle development [52] and reported roles in myogenesis [53,54]. Compared with the above myomiRs, different types of miR-133 were expressed at relatively low levels in the skeletal muscle libraries: there were 3715, 4912 and 7043 reads for miR-133, and 2791,4253 and 5151 reads for miR-133-3p in the three libraries, respectively (Table S3). These variations in global expression of muscle-specific miRNAs could reflect different roles of these miRNAs in myogenesis [13,38]. In addition, miR-107, miR-103a-3p, miR-128, miR-222a, miR-26a-5p, miR-221, miR-146b-5p, miR-16c, and miR-130a, which have been implicated as the muscle-related miRNAs, also showed high count number (Table 3). The fact that muscle-specific or musclerelated miRNAs are highly enriched indicated that the data extracted from the duck breast muscle are credible and that the expression patterns of miRNAs identified in this study likely reflect the expression characteristics of miRNAs in the three muscle tissues.
The main objective of the present study is to identify the miRNAs which are involved in the regulation of embryonic breast muscle development in duck. The expression of identified miRNAs in the three samples was assessed by plotting Log2-ratio and Scatter Plot (Fig. 4). The expression profiles among the three libraries are shown in Table S4. Between E19 and E13 libraries a total of 333 differentially expressed at a higher level miRNAs were observed, of which 207 miRNAs genes were expressed in E13 than in E19 and 126 miRNAs genes were expressed at a lower level in E13 than in E19. Among these differentially expressed miRNAs, 160 miRNAs were extremely significantly different and 6 were significantly different and the fold-change (log2 E13/E19) ranged between 27.26 (miR-202-5p) and 7.95 (miR-2995). For E19 and E27 libraries there were 342 differentially expressed miRNAs, of which 240 were higher expressed in E27 compared to E19 and 102 were lower expressed in E27 compared to E19. A total of 161 miRNAs were extremely significantly different and 3 were significantly different and the fold-change (log2 E27/E19) ranged  from 28.23 (miR-6087) to 8.08 (miR-2995) among the differentially expressed miRNAs. Combining the results of differentially expressed miRNAs among the three libraries and the fact that E19 is the stage with fastest growth in embryonic breast muscle development of duck, the miRNAs with significantly different expression at stage E19 compared to stages E13 and E27 may play important roles in fine-tuning the required gene expression for fast growth of muscle.

Prediction of novel miRNA
The characteristic hairpin structure of a miRNA precursor can be used to predict novel miRNA. The Mireap software [http:// sourceforge.net/projects/mireap/] was used to predict novel miRNAs by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small unassigned highquality reads. The following key conditions are used as criteria for assigning novel miRNAs: (1) the small unassigned high-quality reads must map to an intron region or an antisense exon region of a reference genome; (2) Sequence structures must satisfied the following two criteria: hairpin miRNAs can fold secondary structures and mature miRNAs are present in one arm of the hairpin precursors; (3) The mature miRNA strand and its complementary strand have 2-nucleotide 39 overhangs; (4) Hairpin precursors lack large internal loops or bulges; (5) The secondary structures of the hairpins are steady, with a free energy of hybridization lower than or equal to 218 kcal/mol. Based on the criteria described above, 23 (11, 5 and 18 from E13, E19 and E27, respectively) novel miRNA candidates were obtained (Table 4). Because genomic information of duck are not available, we predicted the novel miRNAs through mapping the unassigned small RNA sequences to the genomes of all animals in UCSC database. Therefore, the novel miRNAs identified in this study are not previously characterized in any other species under the experimental and analytical condition at that time. Li et al. reported that novel miRNAs are lowly expressed in chicken skeletal muscle [55]. In this study, we also found that the read number of each novel miRNA was much lower than the identified miRNAs. Only two novel miRNAs (novel-mir-8 with 3795 reads at E13 and 5154 at E27 and novel-mir-14 with 831 reads at E13) had read counts greater than 100 in the libraries (Table 4). These results indicate that novel miRNAs are masked by highly expressed miRNAs because of the low abundance [56]. Targets prediction showed that the novel miRNA candidates corresponded to 3867, 2680 and 4239 target genes and 8362, 3847 and 13099 genomic loci, respectively (Table 5). For predicted targets of the novel miRNAs, we performed GO (gene ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis. Among the pathways, the metabolic pathways contained the highest number of targets genes (344, 235, and 370 targets in the three libraries respectively) ( Table S5), suggesting that metabolism is a major event at different development stages. In addition, many muscle development relevant pathways, such as MAPK signaling   pathway, mTOR signaling pathway, Wnt signaling pathway and calcium signaling pathway were also discovered from the KEGG pathway analysis of the targets of novel miRNAs suggesting the novel miRNAs predicted in this study may play some roles in muscle development.
Notably, the expression of novel-mir-8 in E13 and E27 were higher than any other identified novel miRNAs and there was no expression at E19. On the contrary, novel-mir-14 was only highly expressed at E19, but was not expressed at E13 and E27. Considering the fact that the embryonic breast muscle growth of duck is fastest at E19 (Table 4), we speculated that novel-mir-8 could be a suppressor of the embryonic breast muscle growth of duck, while novel-mir-14 was a promoter.
To further explore the roles of novel-mir-8 and novel-mir-14 in the embryonic breast muscle growth of Pekin duck, we investigated the function of their targets using the DAVID database for annotation, visualization and integrated discovery (with default parameters) [57]. Because there is no genomic information available for duck, we annotated these targets against the chicken genome. Through the gene name batch viewer analysis, 547 and 878 GenBank_Accession of the targets for novel-mir-8 and novelmir-14 were linked with gene names in GenBank (Table S6). The results of GO analysis showed that there were 124, 27 and 32 terms in the three main categories of biological process, cellular component and molecular function, respectively, for the targets of novel-mir-8 and 241, 34 and 33 terms for that of novel-mir-14, respectively (Table S7). In the KEGG pathway analysis for the targets of novel-mir-8, we found six significant pathways (P,0.05). Among the six pathways, MAPK signaling pathway was the pathways with the highest number of targets (18 targets) (Table  S8). Mounting studies have shown that MAPK signaling pathway is involved in skeletal muscle development and modulation of skeletal muscle fiber type [58][59][60][61][62]. Three pathways including the adipocytokine signaling pathway, were found to be significant enriched (P,0.05) in the KEGG pathway analysis for the targets of novel-mir-14 (Table S9). Previous studies have shown that increased volume and number of adipocyte are positively correlated with leptin production [63,64] and negatively correlated with production of adiponectin (also known as ACRP30), a hormone that decreases hepatic gluconeogenesis and increases lipid oxidation in muscle [65,66]. Leptin signaling has also been implicated in the pathogenesis of arterial thrombosis [67]. Therefore, functional annotation of the novel miRNA targets suggest that novel-mir-8 and novel-mir-14 might play important roles in the embryonic breast muscle development of duck, which promoted us to validate the function of these two novel miRNAs and their targets.

Validation of biological variability between samples of a stage
To minimize the effect of biological variability, three full-sib embryos were used to collect breast muscle samples at stage E13, E19 and E27. Then, we isolated total RNA from these three embryos of a stage and pooled them together to generate a RNA pool of this stage similar to the pooling for the high-throughput sequencing. Because the genetic variability in studies is extensive and inevitable, further validation of miRNAs expression changes need be performed between individuals of a specific stage. Therefore, we randomly selected three highly expressed miRNAs (miR-1, miR-107, and miR-26a-5p) to performed stem-loop qRT-PCR analysis in each sample (Fig. 5). The results showed that there were no significant differences among samples of a stage. This indicates that the effect of biological variability is not significant in this study and the data used in this study is reliable.

Validation of identified miRNAs
Muscle-specific miRNAs are predominantly expressed in muscle-related tissues or organs and are involved in a range of processes including myogenesis (proliferation, differentiation, and fiber type specification), muscle regeneration, hypertrophy, and dystrophy [13,[68][69][70][71]. Therefore, understanding the miRNAs expression pattern can reveal the potential function of the miRNAs. To validate the identified miRNAs in embryonic breast muscle of Pekin duck, stem-loop qRT-PCR analysis of 15 identified duck miRNAs was performed in different tissues or organs (leg muscle, heart, liver, kidney, muscle stomach, small intestine, abdominal fat, skin fat) at E27 and in breast muscle at various developmental stages (E11, E13, E16, E19, E23, E27). Among the 15 miRNAs, 14 miRNAs (93.3%) were in agreement with the expression pattern found in the high-throughput sequencing data (Fig. 6), indicating the high-throughput sequenced data and analysis methods are reliable. Through comparing the 15 miRNAs expression profiles among tissues, we found that the three  muscle-specific miRNAs (miRNA-206, miRNA-1, and miRNA-133) were highly expressed exclusively in in muscle tissue or related organs (breast muscle, leg muscle, and heart), while six myogenesis-related miRNAs (miR-181a-3p, miR-103a-3p, miR-107, miR-10a-5p, miR-222a, and miR-26a-5p) and two highly expressed miRNAs (miR-152 and miR-143) could be detected in all tissues. Interestingly, the expression level of miRNA-152 was approximately equal in all tissues/organs. The remaining 4 miRNAs were not expressed in one or several tissues or organs, like let-7i which had no expression in liver, miRNA-23a were not express in liver and kidney, miRNA-24 hardly showed any expression in liver, kidney, abdominal fat and skin fat and miR-214 could not be detected in liver, kidney, and stomach. The expression of the 15 validated miRNAs were all highly expressed in muscle-related tissues (breast skeletal muscle, leg muscle, and heart) (Fig. 6) suggesting that these miRNAs might play some roles in skeletal muscles development.
To further explore the temporal expression of the 15 miRNAs validated above in the developing embryonic breast muscle of Pekin duck, we performed stem-loop qRT-PCR analysis of the miRNAs in embryonic breast muscle tissues at E11, E13, E16, E19, E23, and E27. The results showed that the analyzed miRNAs could be divided into three groups according to their expression profiling. The expression profiling of miRNAs in group 1 (miR-1, miR-181a-3p, and miR-222a) peaked at E19. Conversely, the miRNAs expression patterns in group 2 (miR-107, miR-10a-5p, miR-214, and miR-143) were lowest at E19. Given that E19 is the fastest growth stage of embryonic breast muscle of Pekin duck according to our prior research, the miRNAs in group 1 might be the potential promoters and the miRNAs in group 2 might be the potential inhibitors for embryonic breast muscle development of duck. The expression of miRNAs in group 3 (let-7i, miR-103a-3p, miR-133, miR-206, miR-26a-5p, miR-24, and miR-23a) increased in the embryonic breast muscles from E11 to E27 (Fig. 6). These results suggested that the miRNAs in group 3 may be involved in the growth and functional maintenance of embryonic breast muscles of duck.

Validation of novel-mir-8 and novel-mir-14 and their targets
Because of the high expression of novel-mir-8 and novel-mir-14 in prenatal breast muscle of duck at some specific stages (novelmir-8 is highly expressed at E13 and E27, novel-mir-14 is highly expressed at E19) and the pathways of muscle development or fat deposition involved in by their targets, we were urged to carry out the expression validation of the two novel miRNAs and some of their targets by stem-loop qRT-PCR. We carried out expression profiles analysis of novel-mir-8, novel-mir-14, mitogen-activated protein kinase kinase 1 (MAP2K1, one of the target of novel-mir-8 in MAPK signaling pathway) and peroxisome proliferatoractivated receptor alpha (PPARa, one of the target of novel-mir-14 in adipocytokine signaling pathway) in embryonic breast muscle of Pekin duck at E11, E13, E16, E19, E23 and E27. The results showed that the expression of novel-mir-8 and PPARa simultaneously reached minimal expression level at E19. In contrast, the expression of novel-mir-14 and MAP2K1 reached maximal expression at E19 (Fig. 7). Evidence showed that activated MAP2K1, also known as MEK1, is a positive regulator in the mid-stage of skeletal myogenesis [72] and that peroxisome proliferator-activated receptor alpha (PPARa), which enhances peroxisomal b-oxidation [73] and activation of lipoprotein lipase [74], can promote fat deposition. Based on the expression profiles of novel-mir-8, novel-mir-14, MAP2K1 and PPARa as well as the roles of MAP2K1 and PPARa as reported previously, we speculated that novel-mir-8 and novel-mir-14 might be the promoter and inhibitor, respectively, for the embryonic breast muscle development in Pekin duck. However, further follow up studies, such as RNAi, are needed to verify this hypothesis.

Ethics statement
All samples' collection and subsequent experiments were approved by the Ethical and Animal Welfare Committee of Beijing, China. Birds involving slaughtering in this study were humanely sacrificed as necessary to ameliorate their suffering.

Tissue collection, small RNA isolation and cDNA libraries construction
Three full-sib embryos of Pekin duck at stage E11, E13, E16, E19, E22 and E27 were collected into sterile physiological saline immediately after removal from the hatching eggs. Pekin duck tissue samples used for stem-loop real-time reverse transcription polymerase chain reaction (qRT-PCR) analysis included leg muscle, heart, liver, kidney, intestines, abdominal fat, and skin fat from stage E27 and breast muscles from stage E11, E13, E16, E19, E22 and E27. These tissues were snap-frozen in liquid nitrogen and stored at 280uC. Total RNA was isolated using TRIzol reagent (Invitrogen) from all the tissue samples. For stage E13, E19, and E27, the isolated total RNA was pooled and used for the generation of the small RNA libraries where the population of recovered small RNAs, ranging in size from 18 to 30 nucleotides, was purified using 15% polyacrylamide gel. Then, 59 adaptors (Illumina, USA) were ligated to the purified small RNAs followed by purification of ligation products on Novex 15% TBE-urea gel. The 59 ligation products were then ligated to 39 adaptors (Illumina) and products with 59 and 39 adaptors were purified using Novex10% TBE-urea gel (Invitrogen). Subsequently, reverse transcription reactions were performed using the RT primer, and PCR reactions were performed using the forward and reverse Illumina primers. The PCR product was purified via phenol/chloroform extraction and ethanol precipitation and was delivered to the Beijing Genomics Institute (BGI) (Shenzhen, China) for sequencing on an Illumina Gnome Analyzer.

Small RNA analysis
The adaptor/acceptor sequences and low quality reads were removed and the contaminated reads formed by the adaptoradaptor ligation was cleaned by a software developed by BGI, hereby the high-quality reads were obtained. Because there is no genomic information of duck available, we compared the highquality reads with GenBank noncoding RNA database (http:// www.ncbi.nlm.nih.gov/) and Rfam database (http://www.sanger. ac.uk/software/Rfam). Sequences with highest similarity to rRNA, tRNA and sncRNA were removed. Subsequently, the high-quality reads were compared with chicken miRNAs, and with other animal miRNAs in miRBase 19.0. Finally, reads number of each miRNA candidates was counted and compared between different tissues. In view of the miRNA diversities in different species, two nucleotide mismatches were allowed in these analysis. The remaining unmatched small RNA fragments were blasted against the chicken genome sequence to identify the exon, intron and repeat sequences using SOAP 2.0 program [28].

Prediction of novel miRNAs
To discover potential novel miRNAs precursor sequences, unique sequences that have more than 10 hits to the chicken genome or match to known non-coding RNAs were removed. Then the flanking sequences (150 nt upstream and downstream) of each unique sequence were extracted for secondary structure analysis with Mfold [http://www.bioinfo.rpi.edu/applications/ mfold] and then evaluated by Mireap [http://sourceforge.net/ projects/mireap/]. After prediction, the potential miRNA loci were examined carefully based on the distribution and numbers of small RNAs on the entire precursor regions. Those sequences residing in the stem region of the stem-loop structure and ranging between 20-22 nt with free energy hybridization lower than 2 18 kcal/mol were considered to be potential novel miRNA candidates [75].

Prediction of novel miRNAS and functional annotation of the target genes
RNAhybrid (http://bibiserv.techfak.uni-bielefeld.de/rnahybrid) [76] was used to predict the targets of novel miRNAs, complying with the following criteria in seed region: (1) No mismatch between 1-9 nt on the 59 end; (2) G-U was permitted, but the number can't exceeds three. The next thing is DAVID (http://david.abcc. ncifcrf.gov/) [57] being used for the functional annotation of the predicted targets. Because there is no genomic information available for duck, we annotated tgohese targets against the chicken genome using the GenBank Accession of the targets of novel miRNAs.

miRNA expression analysis
Comparison of the identified or novel miRNA expression between two samples was conducted to identify those differentially expressed miRNAs. The expression of miRNAs was shown in two samples by plotting a Log2-ratio figure and a Scatter Plot. The procedures were as follows: (1) Normalize the expression of miRNAs in two samples (E19 and E13 or E19 and E27) to get the expression of transcript per million (TPM). Normalized expression (NE) = Actual miRNA count/Total count of high-quality or novel reads. (2) Calculate fold-change and P-value from the normalized expression. Then generate the Log2-ratio plot and scatter plot. Fold-change formula: Fold change = log2 (E13/E19 or E27/E19). P-value formula: The x and y represented normalized expression levels, and the N 1 and N 2 represented total count of a given miRNA in the small RNA libraries of E19 and E13 or E19 and E27, respectively.
Some members of differentially expressed miRNAs were validated using stem-loop qRT-PCR with SYBR Green [77]. The primers for stem-loop qRT-PCR are listed in Table 6. qRT-PCR was carried out with an iCycler IQ5 Multicolor Real-Time PCR Detection System (Bio-Rad, USA) and SYBR Green PCR Master Mix (TaKaRa, Dalian, China) in a 20 ml reaction. The reaction mixtures were incubated in a 96-well plate at 95uC for  s, followed by 40 cycles of 95uC for 10 s, 60uC for 10 s, and 68uC for 20 s. All reactions were run in triplicate. The threshold cycle (Ct) was defined as the cycle number at which the fluorescence intensity passed a predetermined threshold. The quantification of each miRNA relative to U6 gene was calculated using the equation: N = 2 2DDCt .
The qRT-PCR analysis of the target genes of novel-mir-8 and novel-mir-14 The SYBR PrimeScript RT-PCR Kit (TaKaRa, Dalian, China) and a reference gene (b-actin) were used for detecting the expression of MAP2K1 (a of the target of novel-mir-8) and PPARa (a of the target of novel-mir-14). The qRT-PCR reactions were carried out with an iCycler IQ5 Multicolor Real-Time PCR Detection System (Bio-Rad, USA). The qRT-PCR reaction contained 1 mL of cDNA template, 12.5 mL of SYBR Premix ExTaq, 9.5 mL of sterile water, and 1 mL of each gene-specific primer ( Table 6). Thermal cycling parameters were 1 cycle at 95uC for 2 min, 40 cycles of 95uC for 15 s, 60uC for 34 s. Dissociation curve analysis was done after each real time reaction to ensure that there was only one product. qRT-PCR analysis of each sample was repeated for three times. The quantification of each gene relative to b-actin gene was calculated using the equation: N = 2 2DDCt .

Data deposition
Data described in this study is available in the NIH Short Read Archive (SRA) under accession number PRJNA098308.    Table S5 The pathway annotation of the targets of novel miRNAs in E13, E19 and E27.

Supporting Information
(XLSX) Table S6 The gene name batch viewer analysis of the targets of novel-mir-8 and novel-mir-14.