High-Throughput Sequencing Reveals miRNAs Affecting Follicle Development in Chicken

As the derivative of chicken skin, hair follicle is capable of self-renew. Its proliferation and differentiation result in hair formation. MicroRNAs (miRNAs) can effectively regulate gene expression at the post-transcriptional level and play a critical role in tissue growth, development. In this study, we used next generation sequencing technology sequenced miRNAs of the hair follicle derived from the 13 day-old chicken (Gallus gallus) embryos in which from Kirin chicken and Huaixiang chicken that feathers having morphogenesis with significantly different curling. A population of conserved miRNAs was identified. These conserved miRNAs were derived from 638 homologous hairpin precursors across 5 animal species. We identified a total of 645 miRNAs in the chicken embryos. Among them, 11 differentially expressed miRNAs were identified (>±2 Fold, p value <0.05) by comparing Kirin chicken and Huaixiang chicken. Several gene ontology (GO) biology processes and the WNT, BMP and TGF-β signaling pathways were found to be differentially expressed miRNAs as part of hair follicle development process. The miR-1623 has an effect on WNT4 and involved in hair follicle cell development. This study has identified miRNAs that associated with the chick embryonic hair follicle development and identified some target miRNAs for further research into their role played in feather growth.


Introduction
MicroRNAs (miRNAs) are small non-coding RNA molecules that suppress gene expression post-transcriptionally, and function important roles in diverse biological processes [1]. Hundreds of miRNA genes have been found in diverse animals, and many of these are phylogenetically conserved [2]. In addition to endogenous presence in cells, miRNAs can also be actively released into extracellular fluids through exosomes or microvesicles [3,4]. Consequently, miRNA research has become a hot spot in the field of biological for explaining molecular formation mechanisms [5] and important traits of animals [6][7][8][9]. The skin plays an important protection role in animal existence and it evolves with the animal bifurcation.
The feather is one of the most complex integumentary appendages due to the extensive diversity in shape, size, arrangement and pigmentation, and is therefore an excellent model for evolutionary and developmental biology as variations can occur at each step of development and differentiation [10][11][12][13]. The Kirin chicken can adapt to high-temperature environment because of unique frizzled feather branching structure characteristic, rachis stout and outwardly curved, barbs short sparse, feather hook can't connect with the back edge of the adjacent twig lead to pinna can not closed. Recent advances that frizzled feather is caused by KRT75 mutation reside in autosomal, belong to incomplete dominant inheritance [14]. Feathers develop from the hair follicle, therefore the hair follicles numbers, diameter with feathers growth have a direct relationship [15,16]. And miRNAs connected with hair follicles developmental processes [5,17] and regulate hair follicle development and hair growth [18,19] In this study, we investigated the expression profile of miRNAs in the follicles of 13-day chicken embryos from the Kirin chicken (KRC) and Huaixiang chicken (HXC). The results demonstrate that chicken embryonic follicle contains large amounts of miRNAs.

Ethics Statement
All chicken embryos experiments were approved and reviewed by the local ethical committee and the procedures in this study followed the guidelines of the Guangdong Ocean University Animal Care and Use Committee. To minimize the suffering of animals, sodium pentobarbital anesthesia was used before the collection of chicken skin hair follicles samples.

Collection of Chicken Embryonic Follicle Samples
Fertile eggs were collected from 45-wk-old KRC and HXC. Fertile eggs were incubated at 37.8°C with 65 to 75% humidity and intermittent rotations, which to provided 2-3cm skin tissue were obtained from the back of the body of 13-day chicken embryos, separately. Samples of hair follicle stored at −80°C until used.

Small RNA Library Preparation and Sequencing
Collected feather follicle from six chicken embryos at the age of 13-day were used to construct two small RNA libraries in this study. These samples included three KRC ones and three HXC, respectively.
Total RNA was extracted from the follicle using TruSeq Small RNA Sample Pre Kits (Illumine, San Diego, USA) according to the manufacturer's instructions. Total RNA quality was checked with a Bioanalyzer 2100 (Agilent Technologies, USA). The RIN was > 8.0 and A260/A280 was > 2.1 for all samples. The equal concentration total RNAs of six samples were constructed small RNA libraries by TruSeq Small RNA Sample prep Kit of illumina. The overall flow of the sequencing procedure is as follows: small RNAs ranging from 18 to 35nt in length was purified from 15% polyacrylamide gels, then ligated to 5' and 3' adapters. Reverse transcription was performed, and followed by PCR amplification. The purified PCR products (~140bp) were used directly for cluster generation and sequencing analysis using the Illumina's 2000 Sequencer according to the manufacturer's instructions (Personalbio, ShangHai, China).

Sequence Data Analysis
Sequence data analysis was done using AGGT101-miR tool. After deleting poor quality reads, adaptor pollution reads and reads less than 18nt, the clean reads were obtained.
The clean reads of small RNAs were aligned to the reference chicken (G. gallus) genome to identify known miRNAs. The sequences that matched perfectly to known miRNAs (miRBase V21.0) were determined as conserved miRNAs. Other small RNAs (rRNA, tRNA, snRNA and snoRNA) were annotated by blasting against the Rfram, Repbase and ncRNA databases.
To identify differentially expressed miRNAs, the number of conserved miRNAs was normalized to the total number of reads in each sample that matched the chicken (Gallus Gallus) genome. P-values for differentially expressed miRNAs (KRC/HXC) were calculated by Fisher's exact-test and Chi square (2×2) test.

Mirna Target Prediction and Functional Analysis
Target genes of differentially expressed miRNAs were predicted by Target Scan and miRanada. To acquire higher prediction accuracy, only common target genes were considered. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were retrieved using DAVID (http://david.abcc.ncifcrf.gov/).

Quantitative RT-PCR
Total RNAs of sampled follicle were reverse-transcribed by PrimeScript® RT reagent Kit (TAKARA, DRR037A). The primers (Table 1) were designed by Primer 5.0 (ABI). 5ul RT reaction system included: denatured RNA and RT primer (2 uM) 3.0ul, 5×PrimeScript®Buffer 1.0ul, RNase Free dH 2 O 0.6ul, PrimeScript® RT Enzyme Mix I 0.4ul. The RT reactions were performed as follows: 42°C for 15 minutes, 85°C for 5 seconds and hold at 4°C. 20ul real-time PCR reaction system included: 2×SYBR Green Mix with ROX 10.0ul, ddH 2 O 8.2ul, Primer mix (10 uM) 0.8ul, cDNA 1ul. The PCR reactions were performed as follows: 50°C for 2 minutes, 95°C for 2 minutes, then 40 cycles with 94°C for 15 seconds and 60°C for 30 seconds. All experiments were performed on ABI 7900 HT sequence detection system. Each reaction was carried out with 3 replicates. snRNA U6 was used as the control for RT-qPCR. The relative expression level of each miRNA to U6 snRNA was normalized as ∆Cp = Cq miRNA -Cq U6RNA [20]. Comparison of relative expression level in different stages was determined using the 2-∆∆Cp method [21]. Statistical significance analysis of the expression change was performed by one-way ANOVA in SPSS 20.0.

Small RNA Library Construction and Sequencing
To investigate the miRNA expression profile in chicken follicle, High-throughput HiSeq 2000 sequencing yielded 9,714,611 (HXC) and 10,875,477 (KRC) raw reads on average for each group small RNA libraries. After filtered low quality sequences, 9,134,499 (HXC) and 7,919,788 (KRC) clean reads for each group were obtained respectively ( Table  2). The histograms of the reads length distribution showed majority were 20nt ~ 24nt ( Figure 1). Of these, 35,705 (HXC) and 52,201 (KRC) unique small RNAs were identified.

Identification of Conserved Mirnas
To identify conserved miRNAs in chicken follicle, the small RNAs were aligned to current miRBase (Release V21.0). Sequences with perfect matching to known chicken (Gallus Gallus) miRNAs were considered as conserved miRNAs. In total, 645 conserved sequences were annotated as chicken miRNAs. To obtain higher reliable results, only the miRNAs with raw reads >10 at least were considered. Then 290 were sorted as common miRNAs, with only one KRC-specific miRNA (gga-miR-1682). All 291 conserved miRNAs detected by sequencing were listed in Table 3. Table 3 The conserved miRNAs expressed in the chicken follicle miRNA_name miRNA_seq gga-miR-6586-5p TGCTGCCAGATAGAAGTTCACCT gga-miR-1467-5p TCTCAGCTACATCGGTGTAAATC gga-miR-101-1-5p CGGTTATCATGGTACCGGTGCTGT gga-miR-1329-3p CCTCGTAGCTTGATCACGATAT Table 3 gga

Differential Expression Profiles of Conserved Mirnas Between HXC and KRC
To compare the differential expression of miRNAs in the follicle of HXC versus KRC chickens, the numbers of miRNAs in each group samples were normalized to the total number of reads. The expression of one KRC-specific miRNAs was not significant in KRC chickens compared to HXC counterparts. So only the differentially expressed miRNAs in follicle were showed in table 4. In total, 11 miRNAs were considered to be differentially expressed (P< 0.05), with 5 up-regulated and 6 down-regulated. Nine miRNAs had more than two fold expression changes (|log2(fold-change)|> = 1.0) from KRC to HXC (Table 4). We also confirmed the expression patterns by RT-qPCR which of four differentially expressed miRNAs in follicle (Figure 2), the results of which correlated to the RPKM values estimated by RNA sequencing (r = 0.81).

Target Prediction and Functional Analysis of Differential Expression Mirnas
To further explore the roles of differentially expressed miRNAs, putative target genes of the most differentially expressed 9 miRNAs (|Log2 (fold-change)|> = 1.0) were predicted by integrating TargetScan and miRanda. In total, 8250 common target genes were found (data is not shown) which include FZD4, WNT4, BMP and EGF of which related hair follicle development (Table 5). GO annotation showed the putative target genes were significantly enriched (counts > 30, P< 0.05) in biological processes (BP) ( Table 6), Cellular component (CC) (Table 7) and Molecular function (MF) ( Table 8). The KEGG analysis (Table  9) suggested that Focal adhesion, Phosphatidylinositol signaling system, ECM-receptor interaction, Inositol phosphate metabolism, Oocyte meiosis and Ubiquitin mediated proteolysis were the most enriched pathways (counts > 50, P< 0.01). Table 6. GO analysis of the putative target genes in biological processes.

Discussion
In this study, we detected 11 differential expressed miRNAs that were enriched in the KRC and HXC libraries and obtained several predicted target genes that may play different roles in hair follicles formation and development. We also identified several pathways associated with hair follicle cell development, including Focal adhesion, Phosphatidylinositol signaling system, ECM-receptor interaction.
Our analysis of the most abundant mature miRNAs with raw reads>=10 in the KRCs and HXCs identified 290 miRNAs that were common in the two groups. Our analysis of the differential expressed miRNAs revealed 11 miRNAs (table 2) in two groups. These differential expressed miRNAs including miR-1623, miR-184-3p, miR-199b. The miR-184 inhibition argonaute 2 protein expression [22] and the growth of hair follicles [23]. And miRNA-199b has an important role in skin and hair follicle development [24,25]. The miR-1623 target genes that played important roles in Wnt/β-catenin pathway [19]. Through qRT-PCR, we confirmed miR-1623 target gene the WNT4 in chicken and WNT4 showed the higher expression in HXCs ( Figure 3). Thus, the miR-1623 has an effect on WNT4 and involved in hair follicle cell development. We also identified several target genes, such as those encoding FZD4 and EGF, which are related to hair follicle development in chicken. These miRNAs target genes involved several signaling pathways, including WNT, TGF-β, EGF, FGF, BMP, Hox signaling pathway [26][27][28][29][30]. These signaling pathways regulated and transformed hair follicle development in different stages. The WNT and BMP signaling pathways related to differentiation of keratin cell and regulation of the hair shaft formation [31]. And the TGF-β pathway control growth of hair follicles [27]. Many complex factors support hair follicle cell development. To further investigate the functions of miRNAs and their target genes more experiments need to be performed and miRNA knockdown to identify target genes expression levels.

Conclusion
In conclusion, we identified several miRNAs such as miR-1623, miR-184-3p, miR-199b, and their target genes WNT4, FZD4 and EGF which may be involved in hair follicle development in chicken. These data provide a strong foundation for the study of hair follicle development in chicken at the molecular levels.