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.

Keywords:

• LncRNAs
• skin tissues
• chicken

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.

Materials and methods

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 primary-covert 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.

Data processing, reads mapping and identification of novel transcripts

Quality was checked by FastQC (Version 0.11.5; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Trimmomatic (Version 0.36; default settings except MINLEN: 50) and clean paired-end reads were mapped to the latest referenced Ensembl chicken genome (Gallus gallus 5.0) using HISAT software [7]. Restructured transcriptional transcripts were generated by StringTie software [8].

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™ RT reagent Kit (Takara, Japan). β-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^−ΔΔct 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.

To predict mRNA targets of miRNAs, we searched for targets in the 3′-UTR regions of mRNAs. Overlapping target prediction was performed by RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/) and miRanda (http://www.microrna.org/microrna/home.do) [10] software with the following criteria: free energy < −20 kcal/mol.

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 development-related 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].

Results and discussion

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⁻⁴). 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 up-regulated 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.

Figure 1. Number of differentially expressed lncRNAs identified from each comparison.

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.

Table 1. Overlapped lncRNAs among F1 vs. N1, F1 vs. L2, and L2 vs. N1.

lncRNAs	F1 vs. N1	L2 vs. N1	F1 vs. L2
107054950	Up	Up	Up
107053458	Up	Up	Up
101749149	Up	Up	Up
107053752	Up	Up	Up
107054020	Up	Up	Up
107056782	Down	Down	Down
107052435	Down	Down	Down
XLOC_045182	Down	Down	Down
XLOC_229820	Down	Down	Down
XLOC_229997	Down	Down	Down
XLOC_229891	Down	Down	Down
XLOC_229824	Down	Down	Down
XLOC_229822	Down	Down	Down
XLOC_144890	Up	Up	Up
XLOC_145466	Up	Up	Up
XLOC_043059	Up	Up	Up

Note: XLOC represents novel lncRNA, without XLOC represents known lncRNA.

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).

Figure 2. Validation of lncRNAs expression profiles by qPCR. The expression pattern of five lncRNA were tested between F1 and N1 group (a); between L2 and N1 group (b) and between F1 and L2 group (c).

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′-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 β-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 keratin-like, 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.

Table 2. Overlapped lncRNAs that target mRNAs among F1 vs. N1, F1 vs. L2, and N1 vs. L2.

lncRNAs	Target mRNA in F1 vs. N1	Target mRNA in L2 vs. N1	Target mRNA in F1 vs. L2
107054950	COE2	COE2	–
107053458	PAMR1	PAMR1	–
101749149	–	–	–
107053752	COE3	COE3	COE3
107054020	FOXD3	FOXD3	–
107056782	–	–	–
107052435	–	–	–
XLOC_045182	ARHGAP39, FK1L	ARHGAP39, FK1L	FK1L
XLOC_229820	scale keratin-like, feather beta-keratin-like, beta-keratin-related protein	feather beta-keratin-like, Beta-keratin-related protein	–
XLOC_229997	scale keratin-like, beta-keratin-related protein	scale keratin-like, beta-keratin-related protein	–
XLOC_229891	epidermal differentiation protein	epidermal differentiation protein	–
XLOC_229824	scale keratin-like	scale keratin-like	scale keratin-like
XLOC_229822	scale keratin-like	scale keratin-like	scale keratin-like
XLOC_144890	–	–	–
XLOC_145466	–	–	–
XLOC_043059	TMX3, CCDC102B	CCDC102B	CCDC102B

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.

Table 3. Interconnectivity among miRNAs, lncRNAs and mRNAs identified in F1 vs. N1 and L2 vs. N1.

Groups	lncRNA	Expression	miRNA	Expression	Gene	Expression
F1 vs. N1	107052435	Down	gga-miR-31-5p	Down	Inactive phospholipase C-like 2	Up
	107052435	Down	gga-miR-31-5p	Down	CDK15	Down
	107052435	Down	gga-miR-31-5p	Down	ARL14EPL	Up
	107052435	Down	gga-miR-31-5p	Down	LTBP2	Up
	107052435	Down	gga-miR-31-5p	Down	TYRO3	Down
	107054929	Down	gga-miR-200b-5p	Down	SLC22A23	Down
	XLOC_077901	Up	gga-miR-199b	Up	GLPT	Up
L2 vs. N1	107052435	Down	gga-miR-31-5p	Down	TYRO3	Down
	107052435	Down	gga-miR-31-5p	Down	CDK15	Down
	107052435	Down	gga-miR-31-5p	Down	LTBP2	Up

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.

Table 4. Functional enrichment of altered lncRNA highly related to skin and feather development between F1 and N1 group.

Term	Count	P value	LncRNA
GO:0035116∼Embryonic hind-limb morphogenesis	4	0.017	107052525(down);107055065(up);107055066(down); 101751961(down);101751994(down); 107054518(down);107054519(down);XLOC_ 199390(down); 107052611(down)
GO:0008543∼Fibroblast growth factor receptor signalling pathway	4	0.036	107052990(up);107052991(up);107054189(up); XLOC_235660(up); 101748409(up)
GO:0060174∼Limb bud formation	4	0.002	769756(up);101749678(up);417447(up); XLOC_214525(up);107052611(down); XLOC_235660(up)

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).

Table 5. Functional enrichment of altered lncRNA highly related to skin and feather development between L2 and N1 group.

Term	Count	P value	LncRNA
GO:0090090∼Negative regulation of canonical Wnt-signalling pathway	5	0.008	XLOC_214525(up);107051937(up); 107052526(up);107051762 (down); 107052611(down);
GO:0060174∼Limb bud formation	4	0.000	769756(up);XLOC_214525(up); XLOC_235660(up);107052611 (down);
GO:0030855∼Epithelial cell differentiation	3	0.072	XLOC_249634(down);107053920 (down);101752007(down)
GO:0045669∼Positive regulation of osteoblast differentiation	3	0.066	101751961(down);107052611(down); 107053071(down);107053072(down); 107053073(down)
GO:0050680∼Negative regulation of epithelial cell proliferation	3	0.045	XLOC_214525(up);101749492 (up);XLOC_249634(down)
GO:0030326∼Embryonic limb morphogenesis	3	0.037	107052611(down); 107051762(down);101752007(down)
GO:0002076∼Osteoblast development	3	0.007	107052526(up);101751961(down); 107052611(down)
GO:0050673∼Epithelial cell proliferation	2	0.023	XLOC_235660(up);107052611(down)

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).

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.

Table 6. lncRNA-mRNA pairs possibly related to skin and feather development identified between F1 and N1 groups.

Gene	Log2FC	lncRNA	log2FC	Biological process
CTNNB1	−1.3	107052525	−1.5	Embryonic hind-limb morphogenesis
ZBTB16	1.8	107055065	2.0
ZBTB16	1.8	107055066	2.1
MSX-2	−2.5	101751961	−2.5
MSX-2	−2.5	101751994	−2.4
MSX-2	−2.5	107054518	−2.0
MSX-2	−2.5	107054519	−1.7
MSX-2	−2.5	XLOC_199390	−3.8
SHH	−2.8	107052611	−3.4
CTGF	1.8	107052990	2.5	Fibroblast growth factor receptor signalling pathway
CTGF	1.8	107052991	2.6
FGF7	2.5	107054189	1.7
FGF10	2.8	XLOC_235660	3.0
FGF18	1.8	101748409	1.8
HAND2	1.7	769756	1.3	Limb bud formation
HAND2	1.7	101749678	2.1
SOX9	1.4	417447	1.5
SOX9	−1.4	XLOC_214525	−1.1
SHH	2.8	107052611	3.4
FGF10	−2.8	XLOC_235660	−3.0

Table 7. lncRNA-mRNA pairs possibly related to skin and feather development identified between L2 and N1 groups.

Gene	Log2FC	lncRNA	log2FC	Biological process
SOX-9	1.2	XLOC_214525	1.2	Negative regulation of canonical Wnt-signalling pathway
SHH	−2.8	107052611	−2.9
LIMD1	1.4	107052526	1.0
LEF1	−2.0	107051762	−1.6
ISL-1	2.3	107051937	2.4
HAND2	1.2	769756	1.3	Limb bud formation
SOX9	1.2	XLOC_214525	1.2
SHH	−2.5	107052611	−2.9
FGF10	1.8	XLOC_235660	2.5
KRT6A	−1.3	XLOC_249634	−1.8	Epithelial cell differentiation
NR5A2	1.1	107053920	−1.5
DLX6	−2.4	101752007	−2.5
MSX-2	−2.4	101751961	−2.4	Positive regulation of osteoblast differentiation
SHH	−2.5	107052611	−2.9
GJA1	−3.6	107053071	−3.2
GJA1	−3.6	107053072	−2.7
GJA1	−3.6	107053073	−2.3
SOX-9	1.2	XLOC_214525	1.2	Negative regulation of epithelial cell proliferation
AR	1.7	101749492	1.6
KRT6A	−1.3	XLOC_249634	−1.4
SHH	−2.5	107052611	−2.9	Embryonic limb morphogenesis
LEF1	−1.7	107051762	−1.6
DLX6	−2.4	101752007	−2.5
MSX-2	−2.4	101751961	−2.5	Osteoblast development
SHH	−2.5	107052611	−2.9
LIMD1	1.1	107052526	1.0
COL8A1	−1.8	107054903	−1.1	Epithelial cell proliferation
FOXN1	−1.9	107054799	−3.3

Table 8. miRNA-mRNA pairs correlated to lncRNA possibly related to skin and feather development identified between F1 and N1 groups.

Gene	Log2FC	miRNA	log2FC	Biological process
CTNNB1	−1.3	novel_mir_50	−4.9	Embryonic hind-limb morphogenesis
CTNNB1	−1.3	novel_mir_72	−3.8
CTNNB1	−1.3	novel_mir_123	−1.8
FGF18	1.8	novel_mir_176	−3.9	Fibroblast growth factor receptor signalling pathway
FGF18	1.8	novel_mir_72	−3.8
FGF18	1.8	gga-miR-211	−3.0
FGF18	1.8	gga-miR-204	−3.0
HAND2	1.7	novel_mir_50	−4.9	Limb bud formation
HAND2	1.7	novel_mir_176	−3.9
SOX9	1.4	novel_mir_184	−2.9
SOX9	1.4	gga-miR-449a	−2.0

Table 9. miRNA-mRNA pairs correlated to lncRNA possibly related to skin and feather development identified between L2 and N1 groups.

Gene	Log2FC	miRNA	log2FC	Biological process
LIMD1	1.4	gga-miR-365-2-5p	−2.0	Negative regulation of canonical Wnt-signalling pathway
DLX6	−2.4	gga-miR-1552-5p	−1.0	Epithelial cell differentiation
DLX6	−2.4	novel_mir_321	−4.4
MSX-2	−2.4	novel_mir_121	−3.2	Positive regulation of osteoblast differentiation
DLX6	−2.4	gga-miR-1552-5p	−1.0	Embryonic limb morphogenesis
DLX6	−2.4	novel_mir_321	−4.4
LIMD1	1.1	gga-miR-365-2-5p	−2.0	Osteoblast development
MSX-2	−2.4	novel_mir_121	−3.2
AR	1.7	novel_mir_321	−4.4	Negative regulation of epithelial cell proliferation

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 hind-limb morphogenesis, positive regulation of osteoblast differentiation and osteoblast development. Therefore, 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 neural-crest 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.

Authors’ contributions

HL designed the experiments, HL, ST and HY performed sample collection and phenotypic experiments. GF prepared RNA samples for sequencing. LS performed the bioinformatic analysis of transcriptome data and wrote the manuscript. LS interpreted the data. YY provided many good suggestions. HL supervised this project. All authors reviewed, read and approved the final manuscript.

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.

Abbreviations
EF	=	early-feathering
LF	=	late-feathering
F1	=	primary feathers are longer than primary-covert feathers up to 4.5 mm
L2	=	primary feathers are shorter than primary-covert feathers
N1	=	primary feathers do not extend out
FK1L	=	feather keratin 1-like
SHH	=	Sonic hedgehog
FGF10	=	fibroblast growth factor 10

Acknowledgements

We thank the Hengchuang Gene Company for high throughput sequencing and technical support. We thank Jeremy Allen, PhD, from Liwen Bianji, Edanz Group China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

Disclosure statement

The authors declare that they have no competing interests.

Additional information

Funding

This work was supported by the Key Project of Guangdong Province (2016B020233007), the Engineering Centre of Animal Breeding in Guangdong province (2017.1649) and Guangdong Provincial key Laboratory of Animal Molecular Design and Precise Breeding (2019B030301010). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.