Skip to main content
Download PDF

Characterization of microRNA and mRNA expression profiles in skin tissue between early-feathering and late-feathering chickens

BMC Genomics volume 19, Article number: 399 (2018) Cite this article

Abstract

Background

Early feathering and late feathering in chickens are sex-linked phenotypes, which have commercial application in the poultry industry for sexing chicks at hatch and have important impacts on performance traits. However, the genetic mechanism controlling feather development and feathering patterns is unclear. Here, miRNA and mRNA expression profiles in chicken wing skin tissues were analysed through high-throughput transcriptomic sequencing, aiming to understand the biological process of follicle development and the formation of different feathering phenotypes.

Results

Compared to the N1 group with no primary feathers extending out, 2893 genes and 31 miRNAs displayed significantly different expression in the F1 group with primary feathers longer than primary-covert feathers, and 1802 genes and 11 miRNAs in the L2 group displayed primary feathers shorter than primary-covert feathers. Only 201 altered genes and 3 altered miRNAs were identified between the N1 and L2 groups (fold change > 2, q value < 0.01). Both sequencing and qPCR tests revealed that PRLR was significantly decreased in the F1 and L2 groups compared to the N1 group, whereas SPEF2 was significantly decreased in the F1 group compared to the N1 or L2 group. Functional analysis revealed that the altered genes or targets of altered miRNAs were involved in multiple biological processes and pathways related to feather growth and development, such as the Wnt signalling pathway, the TGF-beta signalling pathway, the MAPK signalling pathway, epithelial cell differentiation, and limb development. Integrated analysis of miRNA and mRNA showed that 14 pairs of miRNA-mRNA negatively interacted in the process of feather formation.

Conclusions

Transcriptomic sequencing of wing skin tissues revealed large changes in F1 vs. N1 and L2 vs. N1, but few changes in F1 vs. L2 for both miRNA and mRNA expression. PRLR might only contribute to follicle development, while SPEF2 was highly related to the growth rate of primary feathers or primary-covert feathers and could be responsible for early and late feather formation. Interactions between miR-1574-5p/NR2F, miR-365-5p/JAK3 and miR-365-5p/CDK6 played important roles in hair or feather formation. In all, our results provide novel evidence to understand the molecular regulation of follicle development and feathering phenotype.

Background

Early feathering (EF) and late feathering (LF) are sex-linked phenotypes that are usually determined in chickens at one day of age [1]. Based on the length of primary feathers and primary-convert feathers at hatching, the birds with the early-feathering phenotype are classified into two subtypes: the length of the primary feathers is more than 4.5 mm longer than that of the primary-convert feathers, or the length difference is in the range of 2 mm to 4.5 mm. The birds with the late-feathering phenotype are classified into four subtypes: the length of the primary feathers is longer than primary-coverts within 2 mm, the length of the primary feathers is shorter than that of the primary-covert feathers, the length of the primary feather and primary-covert feathers are the same, or the feathers are missing. In poultry production, the feather phenotype is utilized for sexing at one-day old.

However, the genetic mechanism of late-feathering is unclear. Previous studies have found that the K locus located on the Z chromosome was responsible for the feather phenotype, which was closely linked to the integration of endogenous retrovirus 21 (ev21) [2,3,4,5]. The latest investigations revealed that a tandem duplication of 176,324 bp linked to the K locus results in a partial duplication of Prolactin Receptor (PRLR) and Sperm Flagellar Protein 2 (SPEF2) [6]. Therefore, the function of PRLR and SPEF2 genes was destroyed in late-feathering birds, both of which were regarded as the major candidates impacting feather performance. PRLR is a ligand of PRL, which was involved in more than 300 separate biological activities, including reproduction, metabolism, water and electrolyte balance, growth and development, neurotransmission and behaviour, and immunoregulation [7]. PRLR was also reported to impact hair replacement in mice [8]. SPEF2, an unknown gene, plays a critical role in spermatogenesis and ciliary dyskinesia [9]. Subsequently, researchers focused on examining PRLR and SPEF2 gene expression between early and late feathering birds, and the results revealed that there were significantly different expression profiles in the skin tissue in some breeds [10, 11]. At the same time, many studies began to analyse the association between feather performance and economic production. The results showed that the existence of ev21 caused a reduction in egg production, an increase in infection by lymphoid leucosis virus and an increase in the mortality rate [12]. These negative effects indicated that it was necessary to exploit the genetic mechanism of various feather phenotypes and to clarify how the late feathers are formatted.

MicroRNA (miRNA) is a class of non-coding small RNAs with a length of approximately 22 nt, which plays its post-transcription regulation roles mainly by degrading mRNA or inhibiting the translation process [13]. MiRNA is involved in almost all biological processes [14, 15]. Recent research has revealed that miRNA also participated in follicle development and feather formation [16, 17]. These results indicated that it would be helpful to perform genome-wide transcriptome analysis on miRNA and mRNA in chicken skin tissues, aiming to identify the major genes controlling feather formation and feather phenotypes. Our results would highlight novel genes or pathways to better understand the molecular mechanism of early and late feathering in birds.

Results

Overview of the sequencing data

The small and long RNA sequencing provided by Illumina technology generated an average of 16–19 million single-end high-quality reads and 65–69 million paired-end high-quality reads from each sample. For small RNA sequencing data, 21–24 nt length reads were the most abundant reads, and genome mapping results showed that approximately 44.0% high-quality reads were mapped to exon and intron regions, approximately 30.7% high-quality reads were mapped to annotated miRNAs, and the remaining 25.3% were mapped to other small RNAs (Fig. 1a). Characterization of these 21–24 nt long small RNAs reveals an obvious bias for uracil (U) at their 5′ ends and 3′ ends (Fig. 1b). For long RNA sequencing data, more than 81.9% reads were mapped to the chicken genome in the paired-end model, yielding approximately 56 million reads for gene expression analysis.

Fig. 1
figure 1

Characterization of small RNA sequencing data. a The distribution of raw reads mapped to the chicken genome. b The nucleotide bias of small RNA reads

Full size image

Differentially expressed miRNAs in skin tissue from early-feathering and late-feathering birds

Four hundred and fifty-eight miRNAs were detected in one-day-old chicken wing skin tissues, which accounted for approximately half of the miRNAs reported from chickens. Thirty-three miRNAs were significantly differently expressed between each group (Table 1). Compared to late-feathering birds with no primary feathers extending out (N1 group), the expression of 31 miRNAs was significantly altered in the early-feathering birds, including 11 upregulated miRNAs and 20 downregulated ones. Compared to another type of late feathering birds with primary feathers shorter than primary-covert feathers (L2 group), only 3 miRNAs were different in the early-feathering birds (F1 group). However, there were still 11 differentially expressed miRNAs between the two subtypes of late-feathering birds. Interestingly, all 11 altered miRNAs overlapped between F1 vs. N1 and L2 vs. N1. Gga-miR-31-5p was the only miRNA that overlapped among all the possible comparisons, which was downregulated in F1 vs. N1 and F1 vs. L2, while upregulated in N1 vs. L2.

Table 1 Differentially expressed miRNAs among early and late feathering groups
Full size table

Differentially expressed mRNAs in skin tissue between early-feathering and late-feathering birds

The analysis of differentially expressed genes revealed a significant difference in skin tissues between early-feathering and late-feathering birds. There were 1802 significantly differently expressed mRNAs between the late-feathering L2 group and the N1 group, including 932 upregulated and 870 downregulated genes. In addition, 2893 differently expressed mRNAs were identified between the late-feathering N1 group and the early-feathering F1 group, including 1682 upregulated and 1211 downregulated genes in the F1 group. However, there were only 201 differentially expressed mRNAs between the early-feathering F1 group and the late-feathering L2 group, including 172 upregulated and 29 downregulated genes (Fig. 2a). Among these differentially expressed genes, 114 genes overlapped between the early-feathering F1 group and the late-feathering N1 or L2 group. However, 1637 genes overlapped between F1 vs. N1 and L2 vs. N1, and 55 genes overlapped between F1 vs. L2 and L2 vs. N1. In all, there were 51 genes overlapping among each group (Fig. 2b).

Fig. 2
figure 2

The differentially expressed genes in each group. a The numbers of upregulated and downregulated genes in three pairwise comparisons. b The numbers of differentially expressed genes overlapped in three pairwise comparisons. DEGs of each comparison were selected with a setting of q value < 0.01 and fold change > 2

Full size image

Based on functional annotation, the expression patterns of most feather-follicle-development related genes, except EGF and SPEF2, were similar in the F1 group with the primary feathers longer than primary-covert feathers and the L2 group with primary feathers shorter than primary-covert feathers, compared to the N1 group with no feathers grown out. The results indicated that these genes might generally contribute to feather growth and development, without special roles in primary or primary-covert feathers. Interestingly, PRLR and SPEF2 were reported to be candidate genes related to the late-feathering phenotype [6]. In this study, the expression level of PRLR was significantly lower in the F1 group and the L2 group than in the N1 group, but there was no significant change between the F1 and L2 groups (Table 2). However, SPEF2 had a significantly lower expression level in the F1 group compared to the N1 or L2 group, suggesting that the expression pattern of SPEF2 was strongly associated with the early- and late-feathering phenotypes. In addition, EGF exhibited a completely similar expression pattern to SPEF2, which indicated that EGF might be a novel candidate gene associated with early- and late-feathering phenotypes.

Table 2 The expression profiles of genes related to feather follicle development
Full size table

Validation of differentially expressed miRNAs and mRNAs

To validate the sequencing data, qPCR was employed to test the relative expression profiles of 5 differentially expressed miRNAs and mRNAs among each group in the same skin tissues. The results showed that the alteration of these miRNAs and mRNAs from the RNA-seq data were confirmed by qPCR, and their altered expression patterns among different groups were well-matched with the RNA-seq data, which guaranteed the accuracy of subsequent functional analysis (Fig. 3).

Fig. 3
figure 3

Validation of miRNA and mRNA expression profiles by qPCR. a The expression pattern of five genes were tested between L2 and N1 group. b The expression pattern of five genes were tested between F1 and N1 group. c The expression pattern of five miRNAs were tested between L2 and N1 group. d The expression pattern of five miRNAs were tested between F1 and N1 group. U6 and β-actin were used as reference genes for miRNA and mRNA testing, respectively. Fold change values of qPCR were calculated using the comparative 2–∆∆CT (∆∆Ct = ∆Ct (target gene) – ∆Ct (reference gene)) from at least three independent experiments. Fold change values of sequencing data were calculated using the DESeq2 system

Full size image

Functional analysis of differentially expressed miRNAs and mRNAs

To analyse the molecular function of these differentially expressed miRNAs, both of RNA hybrid and miRanda systems were utilized to improve the prediction of miRNA targets, resulting in 2665 targets potentially regulated by total 33 miRNAs. Functional enrichment analysis revealed that 22 differentially expressed miRNAs mainly contributed to several key molecular processes, including cell differentiation (24 genes), cell fate commitment (9 genes), the Wnt signalling pathway (11 genes), epithelial cell differentiation (8 genes), peptidyl-tyrosine autophosphorylation (7 genes), limb development (7 genes), embryonic limb morphogenesis (6 genes) and negative regulation of osteoblast differentiation (6 genes), of which the latter 6 processes are highly related to skin and feather development (Table 3). In addition, gga-miR-34a-5p and gga-miR-365-2-5p nearly participated in nearly all these processes by targeting 12 and 11 genes, respectively. All these results suggested that several altered miRNAs played important roles in the control of follicle development.

Table 3 Functional enrichment of altered miRNAs
Full size table

The 2893 differentially expressed genes between the F1 and N1 groups were significantly enriched in 17 signalling pathways, including the TGF-beta signalling pathway (24 genes), Tyrosine tyrosine metabolism (11 genes), the MAPK signalling pathway (47 genes) and the Wnt signalling pathway (28 genes; Table 4). However, 1802 differentially expressed genes between the L2 and N1 groups were significantly enriched in 12 signalling pathways, which completely overlapped with the results of the comparison of F1 and the N1 group, except for the PPAR signalling pathway and nicotinate and nicotinamide metabolism (Table 5). Also, even in the same enriched pathway, there were obvious difference between the comparisons of F1 vs. N1 and L2 vs. N1. For example, there were 24 altered genes significantly enriched in TGF-beta signalling pathway from the comparison of F1 vs. N1, and whereas only 13 altered genes were enriched in the comparison of L2 vs. N1.

Table 4 Pathway enrichment analysis on differentially expressed genes between F1 and N1 group
Full size table
Table 5 Pathway enrichment analysis on differentially expressed genes between L2 and N1 group
Full size table

miRNA-mRNA interaction and signalling pathways related to feather development

MiRNAs have roles in biological regulation mainly through interactions between miRNA and its target mRNA, which commonly causes downregulated mRNA expression. Here, integrated analysis of both miRNA and mRNA expression profiles in the same skin tissues was performed to characterize miRNA-mRNA pairs with inverse expression levels and analyse their biological processes. For the comparison between the F1 and N1 groups, 14 pairs of miRNA-mRNA with negative regulation were identified (Table 6). Among these, gga-miR-216a, gga-miR-365-2-5p and gga-miR-130b-5p participate in epithelial cell differentiation through negatively regulating the expression of KRT6A, UPK1B and TCF21, respectively. Gga-miR-1574-5p and gga-miR-365-2-5p are involved in limb development by negatively regulating NR2F2 and PDLIM4. Gga-miR-193a-5p impacts embryonic limb morphogenesis by negatively regulating PRRX1. For the comparison between the L2 and N1 groups, only 5 pairs of miRNA-mRNA were identified, which mainly participate in cell differentiation, limb development and negative regulation of osteoblast differentiation (Table 7). In all, these miRNA-mRNA pairs cooperatively controlled multiple embryo development processes related to wing and follicle development.

Table 6 MiRNA-mRNA pairs with inverse expression relationship between the F1 and N1 groups
Full size table
Table 7 MiRNA-mRNA pairs with inverse expression relationship between the L2 and N1 groups
Full size table

Discussion

The feathering phenotype is a sex-linked trait and can be utilized to distinguish the sexes at hatching based on differences in the rate of feather growth between primary and primary-covert feathers. To understand the molecular mechanism controlling feather development, chicken wing skin tissues at hatching were used for transcriptome analysis from three different groups: primary feathers longer than primary-covert feathers, or no feathers, or primary feathers shorter than primary-covert feathers. In this study, compared to the N1 group, the F1 and L2 groups exhibited faster growth in feather development, and the main differences between F1 and L2 were the growth rate in primary and primary-covert feathers. In addition, the N1 and L2 groups exhibited the late-feathering phenotype, and the F1 group exhibited the early-feathering phenotype. Thus, the candidate genes or miRNAs related to feather development would be identified through comparing the expression profiles among the F1, N1 and L2 groups.

Transcriptome and expression profiling analysis revealed large changes in F1 vs. N1 (2893 DEGs) and L2 vs. N1 (1802 DEGs) but few changes in F1 vs. L2 (201 DEGs). Similar patterns were also observed for the miRNA expression profiles. Among these genes, 1637 genes overlapped between F1 vs. N1 and L2 vs. N1, indicating that these genes were mainly responsible for feather growth and development with no impact on feather types and growth rate. To analyse the difference between early and late feathering groups, 114 genes overlapped between F1 vs. N1 and F1 vs. L2, suggesting that these genes were mainly responsible for feather rate and feather phenotype. These results also revealed a larger number of genes controlling feather development from follicles but a smaller number of genes controlling feather growth rate and feather phenotype. This result is consistent with current genetic mapping of early and late feathering phenotypes. The late feathering phenotype was shown to be controlled by the K allele located on the Z chromosome [2]. Most studies observed that provirus ev21 was closely associated with the late feathering locus [3, 4]. Recent reports revealed a tandem duplication close to the K allele, resulting in the partial duplication of the PRLR and SPEF2 genes [6]. Subsequently, PRLR and SPEF2 were regarded as candidate genes impacting the early feathering and late feathering phenotypes [10, 11, 18]. PRLR is a member of the cytokine receptor family and is closely related to the growth hormone receptor [19]. PRLR knockout investigations in mice revealed a change in the timing of hair follicle cycling events, which indicated a high association between PRL signalling and follicle development [8]. In this study, we also focused on the expression profiles among the F1, L2 and N1 groups. The PRLR gene was 2.1-fold and 1.6-fold lower in the F1 group and the L2 group compared to the N1 group, but no significant change was seen between the F1 and L2 groups. This result suggested that PRLR suppression might contribute to follicle development but has no huge impacts on the growth rate of primary feathers or primary-covert feathers. Thus, the transcriptome data from skin tissues indicated that the feather phenotype was not determined by the PRLR gene. Previous studies showed that the expression level of PRLR in the wing skin of late feathering chickens was 1.7-fold higher than that in the early feathering type of Wenchang chickens at hatching, while a similar investigation in Suqin green-egg-shelled chickens showed no change between early and late feathering types [10, 11]. However, the expression profiles of different types of early feathering birds or late feathering birds were not compared, which should be undertaken in the future. Differently to PRLR, SPEF2 had 7.0-fold and 4.0-fold lower expression in the F1 group compared to the N1 and L2 groups, both of which were of late feathering birds, with no significant difference between each other. The same expression pattern was verified in other studies [10]. It is reasonable to predict that SPEF2 is strongly related to early and late feather formation.

The formation of feathers begins in the follicles. The feather filament grows out of the follicle via epithelial cell proliferation and differentiation. Complex feather development is controlled by multiple genes involved in gene regulation and signalling transition. Through comparing F1 or L2 to the N1 group, many differentially expressed genes and miRNAs were identified in our research that might contribute to the biological processes of feather development. Functional annotation revealed genes enriched in cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction, which is related to PRLR signalling. They were also enriched in growth-related pathways, including the Wnt signalling pathway, the TGF-beta signalling pathway and the MAPK signalling pathway. Additionally, functional analyses of targets of miRNA uncovered 22 differentially expressed miRNAs involved in epithelial cell differentiation and cell fate commitment, as well as embryonic limb development and morphogenesis, which have been interpreted in the linkage of limb development, skin cell differentiation, follicle development and feather formation [20].

MiRNAs, as post-transcriptional factors and key regulators of transcriptional gene networks, play an important role in hair cell generation [21, 22]. Several miRNAs, including the miR-183 family, miR-96, miR-15, miR-99, miR-100, miR-125, and miR-133, all might contribute to hair cell development and maintenance [23,24,25,26]. In domestic animals such as duck, goats and sheep, miRNAs play important roles in the development of follicles [27,28,29]. Here, miRNA-mRNA interaction pairs were analysed to identify the key factors impacting feather formation. Fourteen pairs between the F1 and N1 groups and 5 pairs between the L2 and N1 groups were negatively regulated in the gene regulatory network. We found that NR2F was inhibited by gga-miR-1574-5p in chicken skin tissue to impact limb development. NR2F was reported to be a key factor controlling hair cell reprogramming [30]. Another study found that miR-302 could increase the reprogramming efficiency of hair follicle cells via repression of NR2F2 [31]. In this study, miR-302 exhibited no significant change, indicating that different miRNA-mRNA regulatory interactions controlled feather formation. Gga-miR-365-2-5p was higher in the N1 group compared to the F1 and L2 groups, suggesting that it might inhibit feather growth and development by downregulating the expression of JAK3, PDLIM4, LIMD1, CDK6 and UPK1, which are involved in epithelial cell differentiation and limb development. Interestingly, JAK inhibition was shown to be essential to hair regrowth [32], and CDK6 plays important roles in hair cell differentiation through controlling the cell cycle [33, 34]. All these results revealed that miR-365-5p might have an important effect on feather or hair generation via JAK3 and CDK6 signals.

Conclusions

Transcriptome sequencing of wing skin tissues uncovered different expression profiles among chickens with various feather phenotypes at hatching, resulting in large changes in both miRNA and mRNA expression in F1 vs. N1 and L2 vs. N1 but few changes in F1 vs. L2. PRLR might only contribute to follicle development, while SPEF2 was highly related to the growth rate of primary feathers or primary-covert feathers and could be responsible for early and late feather formation. miR-1574-5p/NR2F, miR-365-5p/JAK3 and miR-365-5p/CDK6 interactions played important roles in hair or feather formation. In conclusion, our results provide new insights to understand the molecular regulation of follicle development and the feathering phenotype.

Methods

Animals

The wing skin tissues were collected in RNAlater (Ambion, Thermo Fisher), solution from 18 Qingyuan partridge chickens at one-day old, provided by the Guangdong Tinoo’s Foods Co., Ltd., including three groups with 6 individuals in each: later-feathering L2 group (L2) with primary feathers shorter than primary-covert feathers, later-feathering N1 group (N1) with no primary feathers extending out, and early-feathering F1 group (F1) with primary feathers longer than primary-covert feathers more than 2 mm.

RNA libraries preparation and sequencing

Nine small RNA libraries and 9 long RNA libraries were prepared with three replicates in each group. Briefly, total RNA was isolated from each skin tissue (approximately 80 mg) using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol. The RNA purity and quantity were tested using a NanoDrop 2000 spectrophotometer based on the OD230/260/280 value. Then, DNA treatment was performed using a DNA-free kit (Ambio, USA). The RNA integrity was checked by an Agilent 2100 Bioanalyser with RNA 6000 Nano Kits (Agilent, USA) with an RNA integrity number (RIN) more than 8. Small RNA libraries were prepared from 1 μg total RNA for each sample using the NEB Next, Multiplex Small RNA Library Prep Set kit for Illumina (NEB, USA) according to standard protocols. Briefly, the 3′ adapters and 5′ adapters were ligated to the 18–30 nt small RNAs selected by 15% PAGE gels from the total RNA. Several procedures were performed, including reverse transcription, and PCR amplification and purification, for constructing the small RNA libraries. Total RNA was purified to remove ribosomal RNA (rRNA) using the Ribo-zero rRNA Removal Kit (NEB, USA). The concentration and quality of the cDNA libraries was assessed by a Qubit 2.0 Fluorometer with the Qubit dsDNA HS Assay Kit (Life Technologies, USA) and an Agilent 2100 Bioanalyser with the High Sensitive DNA kit (Agilent, USA). Nine small RNA libraries and 9 long RNA libraries were generated for high-throughput sequencing using a HiSeq 4000 platform, yielding 50-bp single-end reads and 100-bp paired-end reads, respectively.

Long RNA sequencing data analysis

FastQC (Version 0.11.5; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Trimmomatic (Version 0.36; default setting except MINLEN:50) programs were used in the paired-end mode for obtaining high-quality reads through trimming the low-quality reads (Q20) and contaminating 5′ and 3′ adapters reads [35]. The resulting high-quality reads were aligned using TopHat2 to the latest referenced Ensembl chicken genome (Gallus gallus 5.0) with paired-end mode and fr-firststrand options [36]. Then, the paired-end mapped reads were proceeded by the featureCounts (Version 1.5.1) program for counting the number of fragments per gene per library [37]. In this study, only known coding genes were selected for further analysis such as differentially expressed gene analysis and functional annotation.

Small RNA sequencing data analysis

The raw data of miRNA sequencing was first proceeded by FastQC (Version 0.11.5; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and FASTX-Toolkit (Version 0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/index.html) programs were used with the single-end mode for obtaining high-quality reads through filtering low-quality reads (Q20), contaminating adapters and longer (> 23 nt) or shorter (< 18 nt) reads. Then, the high-quality reads were proceeded by three steps to identify miRNAs. Firstly, Blast software (version 2.2.31) was used to annotate the ncRNAs based on Rfam database (version 12.0) and removed rRNA, scRNA, snoRNA, snRNA, tRNA and other ncRNAs with E-value at 0.01. Then the output reads were aligned to the latest referenced Ensembl chicken genome (Gallus gallus 5.0) and available miRNA dataset (miRBase 21; http://www.mirbase.org/) by using BWA software (version 0.6) subsequently with one mismatch across the entire reads [38]. SAMtools (Version 1.4) was subsequently utilized to count the amounts of read numbers of all the miRNAs detected in each sample [39].

Differentially expressed miRNAs and mRNAs

To identify differentially expressed miRNA and mRNA (DEGs), DESeq2 [40], which is based on the negative binomial distribution, was employed to normalize and evaluate the significant changes in the comparisons of F1 vs. N1, L2 vs. N1, and F1 vs. L2. DEGs of each comparison were selected with a setting of q value < 0.01 and fold change > 2 for further analysis.

Validation of differentially expressed miRNAs and mRNAs by qPCR

To validate the results of high-throughput sequencing data, qPCR was performed to detect the expression patterns of altered miRNAs and mRNAs among each group. Among 33 altered miRNAs, there were 11 miRNAs overlapped between the comparisons of F1 vs. N1 and L2 vs. N1. About half of these overlapped miRNAs were randomly chosen for validation by qPCR. For the confirmation experiment of mRNA expression data, three genes were randomly chosen from overlapped mRNAs with differently expression patterns between the comparison of F1 vs. N1 and L2 vs. N1. Besides, another two reported candidate genes, PRLR and SPEF2, were also used for validation by qPCR. First, total RNA was used to generate the first strand cDNA library according to the stem-loop primer method for miRNA and the random primer method for mRNA using the PrimeScript™ RT reagent Kit (Takara, Japan). U6 and β-actin were used as reference genes for miRNA and mRNA testing, respectively. All the primers are listed in Additional file 1. The 20-μl PCR reaction solution was mixed according to the standard protocol for the SYBR Premix Ex Taq™ II kit (Takara, Japan), including 10 μl SYBR Premix Ex Taq™ II (2×), 0.8 μl forward primer, 0.8 μl reverse primer, 0.4 μl ROX Reference Dye II (50×), 2.0 μl cDNA and 6.0 μl H2O. The PCR reaction was performed on an ABI 7500 system (ABI, USA) as follows: 5 min at 95 °C for initial denaturation, then 35 cycles at 95 °C for 5 s, 57 °C for 30 s and 72 °C for 34 s, finally followed by the dissociation stage (95 °C for 15 s, 60 °C for 1 min and 95 °C for 15 s, then 60 °C for 15 s). The results were calculated and analysed by the 2-ΔΔct method, in which ∆∆Ct = ∆Ct (target gene) – ∆Ct (reference gene) [41]. All data were obtained by repeating experiments three times. Student’s t-test was used to compare expression levels among different groups. The threshold for significance was set at p value < 0.05 and for high significance was set at p value < 0.01.

miRNA target genes and bioinformatic analysis

MiRNA regulation mostly occurs through the interaction of miRNA and mRNA 3’UTR regions. We first downloaded all the chicken UTR sequences from BioMart (http://www.ensembl.org/biomart/martview/). Then, both RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/) and miRanda (http://www.microrna.org/microrna/home.do) software systems were used for predicting the targets of differentially expressed miRNAs [42, 43]. The seed sequences were of complete complementary base pairing without the G:U base pair, and the minimum free energy of the RNA secondary structure was less than − 20 kcal/mol. Only the overlapped targets identified by RNAhybrid and miRanda were selected for functional analysis.

In this study, thousands of altered genes and dozens of altered miRNAs were identified in the comparisons of F1 vs. N1 and L2 vs. N1. Comparing to N1 group with no feathering throw out, F1 and L2 indicates different feathering growth status with longer or shorter primary feathering. To understand the functional mechanism of these altered genes and miRNAs on feathering development, the differently expressed genes or putative target genes of differently expressed miRNAs generated from each paired comparison among F1, L2 and N1 groups were used for Gene Ontology and KEGG pathway enrichment by DAVID 6.7 (http://david.abcc.ncifcrf.gov/), respectively. KEGG results were filtered using p value < 0.05 based upon a Fisher Exact statistic methodology that previously described [44].

Since miRNA commonly acts as a negative regulator, up-regulated miRNA resulted in down-regulated target mRNA, and vice versa. Pearson’s correlation analysis was applied to identify the miRNA-mRNA pairs with inversed expression relationship from two independent data sets of miRNA and mRNA expression profiles in the F1, L2 and N1 groups.

Abbreviations

EF:

Early-feathering

ev21:

Avian endogenous virus 21

F1:

Primary feathers are longer than primary-covert feathers up to 4.5 mm (mm)

L2:

Primary feathers are shorter than primary-covert feathers

LF:

Late-feathering

N1:

Primary feathers do not extend out

PRLR:

Prolactin receptor

qPCR:

Quantitative real-time PCR

RPKM:

Reads per kilobases per million reads

SPEF2:

Sperm flagellar protein

TPM:

Tag per million

References

  1. Radi MH, C. D. Studies on the physiology and inheritance of feathering in the growing chick. Avian Dis. 1992;56:679–705.

    Google Scholar 

  2. Levin I, Smith EJ. Molecular analysis of endogenous virus ev21-slow feathering complex of chickens. 1. Cloning of proviral-cell junction fragment and unoccupied integration site. Poult Sci. 1990;69(11):2017–26.

    Article  PubMed  CAS  Google Scholar 

  3. Levin I, Smith EJ. Association of a chicken repetitive element with the endogenous virus-21 slow-feathering locus. Poult Sci. 1991;70(9):1948–56.

    Article  PubMed  CAS  Google Scholar 

  4. Boulliou A, Le Pennec JP, Hubert G, Donal R, Smiley M. The endogenous retroviral ev21 locus in commercial chicken lines and its relationship with the slow-feathering phenotype (K). Poult Sci. 1992;71(1):38–46.

    Article  PubMed  CAS  Google Scholar 

  5. Bacon LD, Smith E, Crittenden LB, Havenstein GB. Association of the slow feathering (K) and an endogenous viral (ev21) gene on the Z chromosome of chickens. Poult Sci. 1988;67(2):191–7.

    Article  PubMed  CAS  Google Scholar 

  6. Elferink MG, Vallee AA, Jungerius AP, Crooijmans RP, Groenen MA. Partial duplication of the PRLR and SPEF2 genes at the late feathering locus in chicken. BMC Genomics. 2008;9:391.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Goffin V, Bernichtein S, Touraine P, Kelly PA. Development and potential clinical uses of human prolactin receptor antagonists. Endocr Rev. 2005;26(3):400–22.

    Article  PubMed  CAS  Google Scholar 

  8. Craven AJ, Ormandy CJ, Robertson FG, Wilkins RJ, Kelly PA, Nixon AJ, Pearson AJ. Prolactin signalling influences the timing mechanism of the hair follicle: analysis of hair growth cycles in prolactin receptor knockout mice. Endocrinology. 2001;142(6):2533–9.

    Article  PubMed  CAS  Google Scholar 

  9. Sironen A, Kotaja N, Mulhern H, Wyatt TA, Sisson JH, Pavlik JA, Miiluniemi M, Fleming MD, Lee L. Loss of SPEF2 function in mice results in spermatogenesis defects and primary ciliary dyskinesia. Biol Reprod. 2011;85(4):690–701.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Zhao J, Yao J, Li F, Yang Z, Sun Z, Qu L, Wang K, Su Y, Zhang A, Montgomery SA, et al. Identification of candidate genes for chicken early- and late-feathering. Poult Sci. 2016;95(7):1498–503.

    Article  PubMed  CAS  Google Scholar 

  11. Luo C, Shen X, Rao Y, Xu H, Tang J, Sun L, Nie Q, Zhang X. Differences of Z chromosome and genomic expression between early- and late-feathering chickens. Mol Biol Rep. 2012;39(5):6283–8.

    Article  PubMed  CAS  Google Scholar 

  12. Harris DL, Garwood VA, Lowe PC, Hester PY, Crittenden LB, Fadly AM. Influence of sex-linked feathering phenotypes of parents and progeny upon lymphoid leukosis virus infection status and egg production. Poult Sci. 1984;63(3):401–13.

    Article  PubMed  CAS  Google Scholar 

  13. Ambros V. The functions of animal microRNAs. Nature. 2004;431(7006):350–5.

    Article  PubMed  CAS  Google Scholar 

  14. Conrad R, Barrier M, Ford LP. Role of miRNA and miRNA processing factors in development and disease. Birth Defects Res Part C, Embryo Today: Rev. 2006;78(2):107–17.

    Article  CAS  Google Scholar 

  15. Buchan JR, Parker R. Molecular biology. The two faces of miRNA. Science. 2007;318(5858):1877–8.

    Article  PubMed  CAS  Google Scholar 

  16. Andl T, Murchison EP, Liu F, Zhang Y, Yunta-Gonzalez M, Tobias JW, Andl CD, Seykora JT, Hannon GJ, Millar SE. The miRNA-processing enzyme dicer is essential for the morphogenesis and maintenance of hair follicles. Curr Biol. 2006;16(10):1041–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Bostjancic E, Glavac D. Importance of microRNAs in skin morphogenesis and diseases. Acta Dermatovenerol Alp Panonica Adriat. 2008;17(3):95–102.

    Google Scholar 

  18. Bu G, Huang G, Fu H, Li J, Huang S, Wang Y. Characterization of the novel duplicated PRLR gene at the late-feathering K locus in Lohmann chickens. J Mol Endocrinol. 2013;51(2):261–76.

    Article  PubMed  CAS  Google Scholar 

  19. Bole-Feysot C, Goffin V, Edery M, Binart N, Kelly PA. Prolactin (PRL) and its receptor: actions, signal transduction pathways and phenotypes observed in PRL receptor knockout mice. Endocr Rev. 1998;19(3):225–68.

    Article  PubMed  CAS  Google Scholar 

  20. Byrne C, Hardman M, Nield K. Covering the limb--formation of the integument. J Anat. 2003;202(1):113–23.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Kopecky B, Fritzsch B. Regeneration of hair cells: making sense of all the noise. Pharmaceuticals (Basel Switzerland). 2011;4(6):848–79.

    Article  CAS  Google Scholar 

  22. Marson A, Levine SS, Cole MF, Frampton GM, Brambrink T, Johnstone S, Guenther MG, Johnston WK, Wernig M, Newman J, et al. Connecting microRNA genes to the core transcriptional regulatory circuitry of embryonic stem cells. Cell. 2008;134(3):521–33.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Soukup GA, Fritzsch B, Pierce ML, Weston MD, Jahan I, McManus MT, Harfe BD. Residual microRNA expression dictates the extent of inner ear development in conditional dicer knockout mice. Dev Biol. 2009;328(2):328–41.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Mencia A, Modamio-Hoybjor S, Redshaw N, Morin M, Mayo-Merino F, Olavarrieta L, Aguirre LA, del Castillo I, Steel KP, Dalmay T, et al. Mutations in the seed region of human miR-96 are responsible for nonsyndromic progressive hearing loss. Nat Genet. 2009;41(5):609–13.

    Article  PubMed  CAS  Google Scholar 

  25. Friedman LM, Dror AA, Mor E, Tenne T, Toren G, Satoh T, Biesemeier DJ, Shomron N, Fekete DM, Hornstein E, et al. MicroRNAs are essential for development and function of inner ear hair cells in vertebrates. Proc Natl Acad Sci U S A. 2009;106(19):7915–20.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Weston MD, Pierce ML, Rocha-Sanchez S, Beisel KW, Soukup GA: MicroRNA gene expression in the mouse inner ear. Brain Res 2006, 1111(1):95–104.

  27. Zhang L, Nie Q, Su Y, Xie X, Luo W, Jia X, Zhang X. MicroRNA profile analysis on duck feather follicle and skin with high-throughput sequencing technology. Gene. 2013;519(1):77–81.

    Article  PubMed  CAS  Google Scholar 

  28. Ahmed MI, Alam M, Emelianov VU, Poterlowicz K, Patel A, Sharov AA, Mardaryev AN, Botchkareva NV. MicroRNA-214 controls skin and hair follicle development by modulating the activity of the Wnt pathway. J Cell Biol. 2014;207(4):549–67.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Bai WL, Dang YL, Yin RH, Jiang WQ, Wang ZY, Zhu YB, Wang SQ, Zhao YY, Deng L, Luo GB, et al. Differential expression of microRNAs and their regulatory networks in skin tissue of Liaoning cashmere goat during hair follicle cycles. Anim Biotechnol. 2016;27(2):104–12.

    Article  PubMed  CAS  Google Scholar 

  30. Tornari C, Towers ER, Gale JE, Dawson SJ. Regulation of the orphan nuclear receptor Nr2f2 by the DFNA15 deafness gene Pou4f3. PLoS One. 2014;9(11):e112247.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Hu S, Wilson KD, Ghosh Z, Han L, Wang Y, Lan F, Ransohoff KJ, Burridge P, Wu JC. MicroRNA-302 increases reprogramming efficiency via repression of NR2F2. Stem Cells. 2013;31(2):259–68.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Xing L, Dai Z, Jabbari A, Cerise JE, Higgins CA, Gong W, de Jong A, Harel S, DeStefano GM, Rothman L, et al. Alopecia areata is driven by cytotoxic T lymphocytes and is reversed by JAK inhibition. Nat Med. 2014;20(9):1043–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. He Y, Tang D, Cai C, Chai R, Li H. LSD1 is required for hair cell regeneration in zebrafish. Mol Neurobiol. 2016;53(4):2421–34.

    Article  PubMed  CAS  Google Scholar 

  34. Grossel MJ, Hinds PW. From cell cycle to differentiation: an expanding role for cdk6. Cell Cycle. 2006;5(3):266–70.

    Article  PubMed  CAS  Google Scholar 

  35. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Liao Y, Smyth GK, Shi W. The subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41(10):e108.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Li H, Durbin R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. Genome project data processing S: the sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3(6):1101–8.

    Article  PubMed  CAS  Google Scholar 

  42. Kruger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006;34(Web Server):W451–4.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in drosophila. Genome Biol. 2003;5(1):R1.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Huang da W, Sherman BT, Zheng X, Yang J, Imamichi T, Stephens R, Lempicki RA. Extracting biological meaning from large gene lists with DAVID. Curr Protoc Bioinformatics/Editor Board Andreas D Baxevanis [et al]. 2009;Chapter 13:Unit 13–1.

    Google Scholar 

Download references

Acknowledgements

We would like to thank the Hengchuang Gene Company for high throughput sequencing and technical support.

Funding

This work was supported by grants from the National Key Training Project (2014GKXM057), the Key Project of Guangdong Province (2016B020233007), and the Engineering Centre of Animal Breeding in Guangdong province (2017.1649). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

All data supporting the conclusions of this article are included in this article and its Additional file (Additional file 1: Table S1). The sequencing raw data in this study was deposited into the Sequence Read Archive repository (https://www.ncbi.nlm.nih.gov/bioproject/) with accession number PRJNA376472. The N1 group included N11, N12 and N13 individuals; the L1 group included L21, L22 and L23 individuals; and the F1 group included F11, F12 and F13 individuals.

Author information

Author notes

  1. Guijun Fang and Xinzheng Jia contributed equally to this work.

Authors and Affiliations

  1. School of Life Science and Engineering, Foshan University, Foshan, 528231, Guangdong, China

    Guijun Fang, Xinzheng Jia, Hua Li, Shuwen Tan, Hui Yu & Ying Yang

  2. College of Animal Science, South China Agricultural University, Guangzhou, 510642, Guangdong, China

    Guijun Fang, Xinzheng Jia & Qinghua Nie

  3. Guangdong Tinoo’s Foods Limited Company, Qingyuan, 511827, Guangdong, China

    Hua Li, Shuwen Tan & Hui Yu

Authors
  1. Guijun Fang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xinzheng Jia
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Hua Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Shuwen Tan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Qinghua Nie
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Hui Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Ying Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HL and GF designed the experiments, HL, GF, ST and HY performed sample collection and phenotypic experiments. GF prepared RNA samples for sequencing and performed qPCR testing. GF and XJ performed the bioinformatic analysis for transcriptomic data and wrote the manuscript. XJ interpreted the biological data and revised the manuscript. QN and YY provided many good suggestions. HL supervised this project. All authors reviewed, read and approved the final manuscript.

Corresponding author

Correspondence to Hua Li.

Ethics declarations

Ethics approval

The study was approved by the Animal Care Committee of Foshan University (Foshan, People’s Republic of China; Approval ID: FOSU#001, 3 April, 2015). Animals involved in this study were humanely euthanized as necessary to ameliorate their suffering.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Primers for qPCR. (DOCX 21 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fang, G., Jia, X., Li, H. et al. Characterization of microRNA and mRNA expression profiles in skin tissue between early-feathering and late-feathering chickens. BMC Genomics 19, 399 (2018). https://doi.org/10.1186/s12864-018-4773-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-4773-z

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us