Mechanisms of early- and late-feathering in Qingyuan partridge chickens

Abstract Long non-coding RNAs (lncRNAs) are important in feather development and feathering patterns, but few studies on this subject have been conducted in chickens. To understand the follicle development and feathering phenotypes, lncRNA expression profiles in chicken wing skin were determined and combined with previously determined miRNA and mRNA expression profiles of chicken wing skin. We then predicted some regulatory networks among differentially expressed mRNAs, miRNAs and lncRNAs using bioinformatics. Compared to chickens with no feathers growing out (N1 group), 778 lncRNAs were differentially expressed in early-feathering chicks with primary feathers more than 2 mm longer than primary-covert feathers (F1 group), and 443 lncRNAs were differentially expressed in late-feathering chicks with primary feathers shorter than primary-covert feathers (L2 group). Only 45 lncRNAs were differentially expressed between F1 and L2 (fold-change > 2, q < 0.01). The targets of differentially expressed lncRNAs were involved in multiple processes related to feather growth and development. Integrated analysis of lncRNAs, miRNAs and mRNAs showed that 16 pairs negatively and 17 pairs positively interacted in feather formation. XLOC_045182 might inhibit early- and late-feather formation and feather phenotype via FK1L. 107052435 might negatively regulate feather growth and development via gga-miR-31-5p. 107052611 might restrict feather development by regulating SHH expression and XLOC_235660 might have a positive effect on feather development via FGF10.


Introduction
Early-feathering (EF) and late-feathering (LF) are sex-linked phenotypes in chicken, and the sex-linked dominant K locus located on the Z-chromosome is associated with feather development and delayed emergence of primary feathers, while the k þ locus is responsible for fast emergence of feathers [1]. Previous studies have found two candidate genes associated with the K allele for chicken feathering: prolactin receptor (PRLR) and sperm flagellar 2 (SPEF2). Moreover, loss of 149C-terminal amino acids from the PRLR protein resulted in the fast emergence of feathers phenotype [2]. PRLR is a receptor of the anterior pituitary hormone, prolactin, which is involved in many processes, including hair/coat morphology [3].
SPEF2 has a role in sperm formation, normal male fertility and cell differentiation [4], and recent investigations revealed that SPEF2 may be a good candidate gene for chicken feathering [2].
Long non-coding RNAs (lncRNAs), longer than 200 nucleotides, play an important role in chromatin modification, transcription and post-transcriptional processing by regulating gene expression [5]. Several lncRNAs are highly conserved and predicted to be important for the early stages of feather formation [6]. Chen et al. [6] analyzed the expression of lncRNAs in the zebra finch and predicted the function of identified lncRNAs. However, whether feather development in birds is regulated by lncRNAs is an unanswered question. Recently, it was found that lncRNAs inhibit microRNA (miRNA) expression at the transcript level to promote gastric cancer progression [5]. However, very few studies have addressed the mechanism of how lncRNAs and miRNAs regulate the genes responsible for bird feathering.
In our previous study, we analyzed miRNA and mRNA expression profiles in chicken wing skin tissues, and discussed the molecular mechanisms of miRNAs and mRNAs in early-and late-feathering in birds [7]. The purpose of this study is to identify lncRNAs in chickens, predict their function and analyze basic lncRNA-miRNA-mRNA networks to better understand the molecular mechanism of early-and late-feathering in birds.

Chicken samples
Samples from nine one-day-old chicks were obtained from Qingyuan of Guangdong province. These samples were assigned to three groups representing different chicken feathering phenotypes: (1) L2 group (late-feathering chickens with primary feathers shorter than primary-covert feathers) (n ¼ 3); (2) N1 group (late-feathering chicks with no feathers growing out) (n ¼ 3); (3) F1 group (early-feathering chicks with primary feathers more than 2 mm longer than primarycovert feathers) (n ¼ 3). Wing skin tissues from all chickens in the three groups were collected and stored at À80 C. Three samples per group were used for RNA-seq (total ¼ nine samples).

RNA extraction and library preparation for sequencing
Total RNA was isolated from wing skin samples using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. RNA purity was checked using a NanoDrop 2000 spectrophotometer. RNA integrity was assessed using an Agilent 2100 Bioanalyzer with a RNA 6000 Nano Kit (Agilent, USA). An lncRNA library was prepared using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA) according to the manufacturer's protocols. One-hundred-base-pair paired-end reads were obtained using the HiSeq 4000 platform.
The novel transcripts were screened for novel lncRNAs using the following standards: (1) length ! 200 bp; (2) exon number ! 2; (3) Maximum RPKM ! 2; and (4) expressed in at least two samples. The transcripts meeting the above standards were retained and evaluated for coding potential using coding potential calculator [8], coding-non-coding index (CNCI) [9] and pfamscan software [9]. Only those transcripts without coding potential were recognized as lncRNAs.

Identification of differentially expressed lncRNAs and clustering analysis
Differentially expressed (DE) lncRNAs were identified using DEGseq software for each of the following comparisons: (1) N1 group vs. F1 group; (2) N1 group vs. L2 group; (3) F1 group vs. L2 group. DE lncRNAs of each comparison were identified with a threshold of q < 0.01 and fold-change > 2. Furthermore, the hierarchical clustering analysis of lncRNAs was analyzed using 'pheatmap' in the R package. For unsupervised hierarchical clustering, the expression of lncRNAs was normalized using the RPKM method [5].

Validation of differentially expressed lncRNAs by qPCR
In order to validate the results of RNA-seq data, the expression of overlapped differentially expressed lncRNAs among each group in the same skin tissues (n ¼ 6) was tested by qPCR. First, total RNA was used to produce cDNA according to the random primer method using the PrimeScript TM RT reagent Kit (Takara, Japan). b-actin was used as reference gene for lncRNA detection. All the primers are shown in Supplemental Table S1. The qPCR for lncRNA-107052435, 107054950 and 107053458 was performed on an ABI 7500 system (ABI, USA) as follows: 30 s at 95 C for initial denaturation, then 40 cycles at 95 C for 5 s, 57 C for 1 min 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 qPCR for lncRNA-101749149 and 107054020 was performed on an ABI 7500 system (ABI, USA) as follows: 30 s at 95 C for initial denaturation, then 40 cycles at 95 C for 5 s, 51 C for 1 min 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 analyzed using the 2 ÀDDct method. All data were obtained by repeating experiments three times. Student's t-test was used to compare the expression levels among each pair of groups. The threshold for significance was set at p < 0.05.
LncRNA-miRNA-mRNA interaction analysis DE mRNAs targeted by DE lncRNAs were obtained based on the expression correlation coefficient between lncRNAs and mRNAs (Pearson correlation ! 0.93). The miRNAs targeted by DE lncRNAs were predicted using miRanda software [10] with the following criteria: free energy < À20 kcal/mol and score > 140.

Functional enrichment analysis
Gene ontology (GO) enrichment analysis was performed on the DE target genes of DE lncRNAs using DAVID 6.7 (https://david.ncifcrf.gov/) [11], which is based on Fisher Exact statistical methodology. GO enrichment with p < 0.05 was regarded as statistically significant.

Identification of skin and feather developmentrelated lncRNAs
To construct lncRNA-miRNA-mRNA networks related to skin and feather development, we first selected lncRNAs involved in skin and feather development. Next, the lncRNA-mRNA interactions were analyzed, and miRNAs that target these lncRNA-interacting mRNAs were identified. The relationship among these lncRNAs, miRNAs and mRNA were reconstructed using Cytoscape version 3.6.0 [12].

Identification of lncRNAs
To identify lncRNAs, we focussed on unannotated transcripts. First, we applied the coding potential calculator [8] to assess the coding potential by considering the quality of predicted open reading frames (ORFs), and the homology with known proteins and we screened for non-coding transcripts using a cut-off score of negative value. Our second approach was to use a non-coded transcriptional prediction tool, CNCI (coding-non-coding index) [9], to estimate the coding potential of the transcripts. The third approach was to evaluate putative non-coding transcripts with similar reading frames in the Pfam HMM database (E-value < 10 À4 ). From the overlapping results of the three approaches, we identified 1820 putative lncRNAs.
Genome-wide identification of differentially expressed lncRNAs between early-feathering and late-feathering chickens The growth rate of feathers has been related to carcase weight, sexual maturity, egg production and disease resistance. To understand the molecular mechanisms controlling feather development, lncRNA expression profiles from three different groups: primary feathers longer than primary-covert feathers (F1), or no feathers (N1), or primary feathers shorter than primary-covert feathers (L2) were identified. A total of 778 lncRNAs were differentially expressed in the F1 group compared with the N1 group, including 331 upregulated and 447 down-regulated lncRNAs. In the L2 group, 443 lncRNAs were differentially expressed compared with the N1 group. Only 45 lncRNAs were differentially expressed between groups F1 and L2. A summary of the up-/down-regulation is shown in Figure 1.
Among these lncRNAs, 372 overlapped between F1 vs. N1 and L2 vs. N1, revealing that these lncRNAs were mainly responsible for feather growth and development with no influence on feather type or growth rate. To analyze the difference between early-and late-feathering groups, we identified 32 lncRNAs that overlapped between F1 vs. N1 and F1 vs. L2, hypothesizing that these lncRNAs were mainly responsible for feather growth rate and feather phenotype. Sixteen DE lncRNAs overlapped among F1 vs. N1, L2 vs. N1 and F1 vs. L2 (Table 1), indicating that these lncRNAs were mainly responsible for not only feather growth and development, but also feather types and feather phenotype. Consistent with a previous study [1], the majority of lncRNAs were significantly differentially expressed in avian feather development.

Validation of differentially expressed lncRNA
To test the RNA-seq results, the expression of five overlapped differentially expressed lncRNAs in each pair of groups was detected by qPCR in more samples (n ¼ 6 per group). The results showed that the expression trend of these lncRNAs was consistent with the RNA-seq data ( Figure 2).

LncRNA-miRNA-mRNA interaction
Many recent studies have shown that lncRNAs can act as competing endogenous RNAs to regulate target mRNAs, thereby playing an important role in diverse biological processes. miRNAs can inhibit target mRNA translation and silence target expression by binding to the 3 0 -untranslated region of the target mRNA. Moreover, lncRNAs can target miRNAs in the indirect regulation of gene expression.
To analyze the molecular function of these common DE lncRNAs (Table 2), lncRNA targets were predicted. Interestingly, our previous study showed that PRLR was found significantly decreased in the F1 and L2 groups compared to the N1 group and SPEF2 was significantly decreased in the F1 group compared to the N1 or L2 group [6]. In this study, we found that PRLR and SPEF2 are both positively regulated by lncRNA-XLOC_235618. PRLR and SPEF2 were reported to be candidate genes related to the late-feathering phenotype [12]. LncRNA-XLOC_235618 might play important role in feather development. Twelve target DE mRNAs potentially regulated by 11 common lncRNAs were identified ( Table 2). Among these, XLOC_045182 regulates the expression of feather keratin 1-like (FK1L), a member of b-keratin family involved in feather evolution, to impact early-and late-feather formation and feather phenotype. XLOC_229824 and XLOC_229822 regulate the expression of scale keratinlike, which were down-regulated in F1 compared with N1 and L2, and L2 compared with N1 groups.  XLOC_043059 is up-regulated in chicken feather development and positively regulates the expression of differentially expressed coiled-coil domain-containing protein 102B (CCDC102B) mRNA, lncRNA-107053752 regulates transcription factor COE3. All these results indicate that several lncRNAs play important roles in the control of early-and late-feather formation and feather phenotype through positively regulating the expression of several DE mRNAs. However, we found that XLOC_229820, XLOC_229997, XLOC_229891, 107054950, 107053458 and 107054020 target mRNAs were only differently expressed in F1 vs. N1 and L2 vs. N1, indicating that these lncRNAs mainly contribute to feather growth and development by regulating the differential expression of mRNAs. Next, we analyzed the regulation relationships to identify interactions among DE lncRNAs, miRNAs and mRNAs. We identified seven pairs of inter-connectivity among DE lncRNAs, miRNAs and mRNAs in F1 vs. N1 and L2 vs. N1 (Table 3), in which a lncRNA positively regulates a miRNA and, in turn, the miRNA regulates a target gene. For comparison between the F1 and N1 groups, seven pairs of lncRNA-miRNA-mRNA were identified (Table 3), including three lncRNA-miRNA-mRNA pairs with negative regulation. For comparison between the L2 and N1 groups, only three lncRNA-miRNA-mRNA pairs were identified, including one pair with negative regulation. Among these, down-regulated lncRNA-107052435 was shown to positively regulate gga-miR-31-5p. We also found that inactive phospholipase C-like 2, ADP ribosylation factor like GTPase 14 effector protein like (ARL14EPL), and latent transforming growth factor beta binding protein 2 (LTBP2) were up-regulated and cyclin dependent kinase 15 (CDK15) and TYRO3 protein tyrosine kinase (TYRO3) were down-regulated by gga-miR-31-5p in chicken skin tissue. miR-31 is a major miRNA involved in hair follicle growth, and in an miR-31 transgenic mouse model keratinocytes are impaired [8]. This leads to abnormal proliferation, apoptosis and differentiation that lead to altered hair growth, while the loss of miR-31 results in increased hair growth. These results indicate that lncRNA-107052435 might inhibit the processes of feather growth and development via gga-miR-31-5p.

GO enrichment analysis of DE lncRNAs
Here, we describe the enrichment results of DE lncRNAs targeting DE mRNAs. In early-feathering chicks with feathers growing out (F1 vs. N1), functional enrichment analysis revealed that DE lncRNAs mainly contributed to 56 key biological processes, including transcription, DNA-templated (26 genes); positive regulation of transcription from RNA polymerase-II promoter (19 genes); signal transduction (12 genes); embryonic hind-limb morphogenesis (four genes); fibroblast growth factor receptor signalling pathway (four genes) and limb bud formation (four genes), of which the latter three processes are related to skin and feather development (Table 4). These results indicate that some DE lncRNAs might be associated with follicle development.
In the late-feathering chicks with feathers growing out (L2 vs. N1), DE lncRNAs were involved in several processes, including negative regulation of canonical Wnt-signalling pathway (five genes), limb bud formation (four genes), epithelial cell differentiation (three genes), positive regulation of osteoblast differentiation (three genes), negative regulation of epithelial cell proliferation (three genes), embryonic limb morphogenesis (three genes), osteoblast development (three genes) and epithelial cell proliferation (two genes), - which are related to skin and feather development ( Table 5).

Analysis of lncRNAs related to skin and feather development
To further explore the lncRNA mechanism related to skin and feather development, the lncRNA-mRNA interactions were analyzed (Table 6 and Table 7), and lncRNA-associated genes targeting miRNAs were identified (Table 8 and Table 9). The relationships among these lncRNAs, mRNAs and miRNAs were constructed using Cytoscape version 3.4.0 ( Figure 3). We found that lncRNA-107052611, lncRNA-XLOC_ 214525 and lncRNA-101751961 are involved in six, three and three biological processes related to skin and feather development, respectively (Table 6 and  Table 7). The target genes of lncRNA-107052611 were enriched in six biological processes, including embryonic hind-limb morphogenesis, limb bud formation, negative regulation of canonical Wnt-signalling, positive regulation of osteoblast differentiation, embryonic limb morphogenesis and osteoblast development. The target genes of lncRNA-XLOC_214525 were enriched in limb bud formation, negative regulation of canonical Wnt-signalling and negative regulation of epithelial cell proliferation pathways. The target genes of lncRNA-101751961 were enriched in embryonic hindlimb morphogenesis, positive regulation of osteoblast differentiation and osteoblast development. Therefore,         Negative regulation of epithelial cell proliferation Figure 3. Network construction of lncRNA-mRNA-miRNA interaction possibly related to skin and feather development identified between groups. Red nodes represent down-regulated genes (F1 or L2 group compared to the N1), and green nodes represent up-regulated genes (F1 or L2 group compared to the N1). Large nodes indicate bigger degrees. The diamond-shaped nodes represent lncRNAs, the round rectangle-shaped nodes represent microRNAs. The Ellipse-shaped nodes represent genes. The number of edges between two nodes represents the number of their enrichment GO terms involved in skin and feather development.
these lncRNAs might play an important role in skin and feather development, but further study is needed. According to these results (Table 8 and Table 9), several mRNAs associated with lncRNAs related to skin and feather development can be regulated by miRNAs. Fibroblast growth factor 18 (FGF18) and Catenin beta 1 (CTNNB1) were identified to be regulated by four and three miRNAs, respectively. FGF18 was regulated by DE novel_mir_176, novel_mir_72, gga-miR-211 and gga-miR-204, and by lncRNA-107052525, which is enriched in embryonic hind-limb morphogenesis ( Table 6). CTNNB1 was regulated by DE novel_mir_50, novel_mir_72 and novel_mir_123, and its target lncRNA is 101748409, which is involved in the fibroblast growth factor receptor signalling pathway ( Table 6).
Overlapped GO enrichment for limb bud formation was found in F1 vs. N1 and L2 vs. N1, for 107052611, 769756, 101749678, 417447, XLOC_214525 and XLOC_235660, the latter five of which were up-regulated. Limbless chicken mutants have no limb bud outgrowth and no outgrowth of the wing and thus no flight feathers [13]. When such mutant phenotypes are restored using FGF, double lines of feathers appear at the posterior margin of the developed forelimb [13]. These results revealed a correlation between limb bud formation and feather formation. Our analysis showed that these lncRNAs potentially positively regulate four genes (SOX9, SHH, HAND2 and FGF10). The results show that 769756, 101749678, 417447, XLOC_214525 and XLOC_235660 might be responsible for feather growth and development by upregulating the expression of transcription factor SOX9, heart-and neuralcrest derivatives-expressed protein 2 (HAND2) and fibroblast growth factor 10 (FGF10), which are all involved in limb bud formation, but 107052611 might inhibit feather growth and development by upregulating the expression of Sonic hedgehog (SHH), which is enriched in limb bud formation. SHH, a determinant of the anteroposterior axis in the limb, suppresses limb development [14]. In this study, we also found that SHH gene expression was 2.8-fold and 2.5-fold lower in the F1 group and the L2 group, respectively, compared with the N1 group, but there was no significant difference between the F1 and L2 groups. lncRNA-mRNA interaction analysis showed that SHH was positively regulated by lncRNA 107052611. SHH suppression might, therefore, contribute to feather development but has no significant effects on feather type or growth rate, and lncRNA 107052611 might restrict feather development by regulating SHH expression. FGF10 is a paracrine-signalling molecule found in limb bud development [15]. Interestingly, over-expression of FGF10 promotes the expansion of proximal components in the follicle, which indicates a strong association between FGF10 and feather development [16]. LncRNA-mRNA interaction analysis identified FGF10 to be positively regulated by XLOC_235660. These results show that lncRNA XLOC_235660 might have an important effect on feather development via FGF10. We also found that HAND2 and SOX9 were negatively regulated by miRNAs (novel_mir_50, novel_mir_176, novel_mir_184 and gga-miR-449a).

Conclusions
Transcriptome sequencing of wing skin tissues revealed large differences in lncRNAs expression between F1 vs. N1 and L2 vs. N1, but few differences between F1 vs. L2. XLOC_045182 might inhibit early-and late-feather formation and feather phenotype via FK1L. 107052435 might negatively regulate feather growth and development via gga-miR-31-5p. 107052611 might restrict feather development by regulating SHH expression. XLOC_235660 might also have a positive effect on feather formation and feather phenotype via FGF10. These results provide new insights into the molecular mechanism of follicle development and the feathering phenotype.

Availability of data and materials
The raw sequencing data in this study is 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.

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 to ameliorate their suffering.