Dynamic transcriptome profiling towards understanding the morphogenesis and development of diverse feather in domestic duck

Background Feathers with complex and fine structure are hallmark avian integument appendages, which have contributed significantly to the survival and breeding for birds. Here, we aimed to explore the differentiation, morphogenesis and development of diverse feathers in the domestic duck. Results Transcriptome profiles of skin owing feather follicle from two body parts at three physiological stages were constructed to understand the molecular network and excavate the candidate genes associated with the development of plumulaceous and flight feather structures. The venn analysis of differentially expressed genes (DEGs) between abdomen and wing skin tissues at three developmental stages showed that 38 genes owing identical differentially expression pattern. Together, our data suggest that feather morphological and structural diversity can be possibly related to the homeobox proteins. The key series-clusters, many candidate biological processes and genes were identified for the morphogenesis, growth and development of two feather types. Through comparing the results of developmental transcriptomes from plumulaceous and flight feather, we found that DEGs belonging to the family of WNT, FGF and BMP have certain differences; even the consistent DEGs of skin and feather follicle transcriptomes from abdomen and wing have the different expression patterns. Conclusions Overall, this study detected many functional genes and showed differences in the molecular mechanisms of diverse feather developments. The findings in WNT, FGF and BMP, which were consistent with biological experiments, showed more possible complex modulations. A correlative role of HOX genes was also suggested but future biological verification experiments are required. This work provided valuable information for subsequent research on the morphogenesis of feathers. Electronic supplementary material The online version of this article (10.1186/s12864-018-4778-7) contains supplementary material, which is available to authorized users.


Background
Avian feathers are branched integumentary appendages that play important roles in flight, protection, heat retention, communication, and mate attraction, due to they could form different structures in various body parts or at different phases during their life stages [1]. The prosperous feather diversities have contributed significantly to the rapid and extensive radiation of birds to become the dominant vertebrate [2]. As the most intricate integumentary appendages with diversiform shapes, arrangements and pigmentations, feathers are regarded as an excellent model for evolutionary and developmental biology research [3].
Feathers are divided into three main types of contour feathers (pennaceous), downy feathers (plumulacuous) and filoplumes. Contour feathers not only include the ordinary body contour feathers, but also flight feathers (remiges) and rectrices. Feather branching at three levels, the rachis, barbs and barbules offer more opportunities for diversity. Modulating their size, angle and symmetry generates complex feather forms. The downy feathers are radial symmetry, and their barbules do not possess the hooklets, hence a fluffy structure could be formed to keep body warm [4]. After processing, the downy feather is well suited for the production of down jacket. While strong flight feathers are left-right asymmetry, which have a larger pennacuous vane (barbules interlock with each other through hooklets), and a longer and thicker rachis for insertion deeper into the follicle. The structure of the flight feathers could anchor more securely to hold its aerodynamic function [5].
The morphogenesis of feather is initiated by the interplay of epithelia and subjacent mesenchayme and usually involves a series of dynamic cellular processes [2]. There are several key molecules that controlled the fundamental aspects of these processes, including growth factors and their receptors, cell adhesion molecules and their ligands, signal transduction molecules and transcription factors [6][7][8][9][10]. Recent studies of feather morphogenesis have concentrated on the molecular mechanisms underlying placode induction and feather bud formation. Little attention has been paid to the regulatory networks that pattern and define the morphogenesis of the elaborate feather structure. Only a handful of signal transduction molecules, cell junctions, feather keratins and other few genes have been reported to possibly involved in the formation of diverse feathers. Previous studies revealed that BMP and SHH signaling are involved in regulating the formation and balance among the rachis and barbs [7]. A feather morphogenesis model suggests that plumulaceous feather structure evolved by the establishment of activator-inhibitor interactions between SHH and BMP2 signaling in the basal epithelium of the feather germ [11]. In addition, the mis-expressed of BMP4 would enhance the rachis formation and barbs fuse. When noggin (a BMP antagonist) is mis-expressed, the rachis is split and increased barb branching ensues [7]. Perturbing the gradient of WNT3A converts bilaterally symmetric feathers into radially symmetrical feathers [12]. The complex branching pattern of feathers may derive from the establishment of specific cell junctions among barb/ barbules cells. Gap junctions serve in cell communication while tight junctions stabilize the complex branching of keratinized feather cells [13]. Feathers consist mainly of flexible corneous materials made of αand β-keratin multigene families [14]. Ng et al. [15] suggested that feather keratins on chromosome 2 of Gallus gallus may have significant effects on the formation of stiff feather structures, and feather keratin on chromosome 25 may be required for softer textures. Furthermore, the crest phenotype is caused by a cis-acting regulatory mutation underlying the ectopic expression of HOXC8 [16].
Domestic duck feathers provide an excellent system for studying the formation and development of morphological complexity because it has diverse forms in different body parts or at different phases during its whole life. Comprehensive cataloguing of gene expression changes of skin and follicle tissues of domestic duck at different physiological stages or body parts would be beneficial to further understand the complex structure of feathers. The availability of transcriptomic analysis provides an excellent opportunity to study gene expression patterns that potentially account for the diverse feather forms [17].
Our objective was to analyze the biological information regarding the transcriptional profiles of skin and follicles to further identify differentially expressed genes and key gene regulatory networks affecting the differentiation and development of diverse feathers through the domestic duck model. We characterized and quantified mRNAs that are expressed in the skin and follicles during feather development in Cherry Valley ducks at two parts of body (abdomen and inner side of the wing) during three stages of feather development. Two body parts of skin and follicles were selected to represent different feather types. We made three analyses: 1) the developmental transcriptomes of abdomen skin and follicles for capturing the functional genes and molecular events accounting for the morphogenesis and development of plumulaceous feather; 2) the developmental transcriptomes of inner side of the wing skin and follicles for obtaining the functional genes and molecular events associated with the morphogenesis and development of flight feather. 3) the comparisons of abdomen and wing skin transcriptomes at the same physiological stages for understanding the differentiation of plumulaceous feathers and flight feathers; Our work has firstly presented a large-scale genomics resource for understanding the evolution and morphogenesis of various feather types.

Results
Overview of RNA-Seq data A total of 12 libraries were sequenced from skin owing feather follicle tissues of six groups (n = 2 for each), including skin owing feather follicle tissues of the abdomen and inner side of the wing at three different physiological stages (3-day-old, 27-day-old and 6-month-old), representing the early growth plumulaceous feathers (EP), middle growth plumulaceous feathers (MP), late growth plumulaceous feathers (LP), early growth flight feathers (EF), middle growth flight feathers (MF) and late growth flight feathers (LF), respectively. About 4.00-4.36 Gb raw reads were produced for each library. After discarding low-quality reads, RNA-seq  yielded from 42,914,954 to 47,459,658 clean reads with  average about 97% Q20 bases for each sample, which were used for all further expression analysis. Among the total number of clean reads from 12 samples, 69.45 to 80.22% were successfully mapped against the reference Peking duck genome ( Table 1). The percentage of the unique mapping reads is approximate 71.93% in each sample. Two issues that may be considered as possible reasons for the lower mapped rates are the draft assembling and annotation of Peking duck genome and evolutionary divergence between Cherry Valley duck and Peking duck. Gene structure analysis was performed for each group (Additional file 1: Figure S1), most of the clean reads (43.63-54.50%) were aligned to the exon regions of the reference genome, followed by the number of the clean reads matched to the 3'UTR and intergenic regions. Correlation of transcript expression level is a crucial indicator for the reliability of the experimental results and the rationality of sampling. The Pearson correlation coefficient between two biological replicates of six groups in this study had very high repeatability (i.e., all R 2 ≥ 0.8804; Additional file 2: Figure S2).

Pairwise differential expression analyses and validation of Solexa sequencing data
In this study, we used DEGseq to screen the differentially expressed genes (DEGs) to determine the differences in different feather types from different body parts or developmental stages. In the plumulaceous feather development, a total of 4756 genes were found to be DEGs; of these, 1333, 2614, 3231 genes were obtained in the comparison of MP versus EP, LP versus MP; LP versus EP, respectively. While in the process of the flight feather development, transcriptomic analyses identified a total of 5823 DEGs; of these, 2978, 2088, 4036 genes were screened in the comparison of MF versus EF, LF versus MF, LF versus EF, respectively. Transcriptomic analyses identified a total of 727, 96, 375 DEGs over 2-fold change that are possibly related to the differentiation of downy feather and flight feather in the comparisons of EF versus EP, MF versus MP, LF versus LP, respectively (Fig. 1). Ten DEGs were randomly selected for qRT-PCR quantification from the same RNA sample of MF and MP for the purpose of validating whether our sequencing and analysis were reliable. In general, high linear correlations and Pearson correlation coefficient of fold-change (FC) values from the two methods indicated that our transcriptome sequencing was reliable and we can make reasonable deductions from the gene expression values generated from RNAseq (Additional file 3: Figure S3).

Functional genes and molecular events accounting for the development and morphogenesis of plumulaceous feather
The downy feathers of ducks are widely used as the filling materials of the insulation clothing. It is of great significance to study the morphogenesis and development of the plumulaceous feathers and identify the candidate genes participated in the process of plumulaceous production. In this part, we did some further analyses about three abdomen skin and follicle transcriptomes at different physiological stages, including series-cluster analysis, functional annotation, gene-act-network and gene co-expression analysis.

Series-cluster analysis and functional annotation of the clusters
The expression patterns not only indicate the diverse and complex interactions among genes, but also suggest that genes with similar expression patterns may have Unique mapped reads, reads that matched the reference genome in only one position similar functions in the feather molting and growth. In the plumulaceous feather libraries (P libraries), a total of 4756 genes were found to be DE (Additional file 4: Figure S4). Then 8 series-clusters were obtained based on the 4756 DEGs ( Fig. 2a; Additional file 5: Table S1). Each gene cluster exhibited a distinctive expression pattern. Largest group of P libraries is cluster 4 with 1245 (26.2%) genes, which maintained a relative stable expression from stage 1 to stage 2, and then rose at stage 3. GO enrichment analyses about 8 clusters were conducted for further understanding the gene expression pattern associated with the feather growth and development. Cluster 4 Fig. 1 The number of differentially expressed genes between the comparison libraries. Total DEGs (grey), up-regulated genes (red), and down-regulated genes were presented by histogram contained the largest number of DEGs that participate in the biological process of feather growth and development; these genes were associated with enriched GO terms including cell adhesion, epidermal growth factor receptor signaling pathway, axon guidance, extracellular matrix organization and BMP signaling pathway. Seven GO terms were enriched in cluster 3, including hair follicle development, epidermis development, WNT signaling pathway, keratinocyte differentiation, keratinization, skin development and epidermis morphogenesis. Meanwhile, no GO term that is associated with feather growth and development was enriched in cluster 6 ( Fig. 2c).

GO and pathway analysis of DEGs in P libraries
The significant enriched GO terms (p-value < 0.05) involved in the plumulaceous feather development were determined from the DEGs of MPvsEP, LPvsMP and LPvsEP (See for full list of GO terms in Table 2 Table S2). These pathways are suggested to be important during the morphogenesis and growth of the plumulaceous feather in duck, including Adherens junction, MAPK signaling pathway, Hippo signaling pathway, Hedgehog signaling pathway, cell adhesion molecules (CAMs), ECM-receptor interaction, TNF signaling pathway, Jak-STAT signaling pathway, and Wnt signaling pathway. A global pathway net was constructed based on the significant pathways to illustrate the key pathways in the process of plumulaceous feather development (Additional file 7: Figure S5). Adherens junction and Wnt signaling pathway were considered to be the most important node in the net because the component exchanges with other pathways were strongly dependent on its existence. Adherens junction, one of the cell junctions, exists in the different types of epithelial cell (including hair follicle epithelium) [18]. In the process of hair follicle morphogenesis and development, Wnt pathway plays an essential role during hair follicle induction and is considered to be master regulator during hair follicle morphogenesis [19].

Gene-act-network
After functional analysis, it is important to explore the relationships among the DEGs involved in the aforementioned biological processes. In the gene interaction network, a wide variety of Wnt family (WNT5A, WNT5B, WNT6, WNT10A, WNT11, WNT16), and its receptors of frizzled protein (FZD1, FZD2, FZD3, FZD5, FZD7, FZD8, FZD10) and DKK protein (DKK1, DKK2), CTNNB1, AXIN2, BAMBI, WIF1, TCF7L1, LEF1 were involved in the key pathways that previously mentioned including Wnt signaling pathway and adherens junction (Fig. 3a). WNT5A is confirmed to express almost in the mesenchymal and epidermal cell before the feather bud formation, and is a target of SHH in hair follicle morphogenesis [20]. The quantity of downy feather is mainly affected by the follicles development in birds, WNT6 plays a key role in follicular development as an intercellular signaling molecule, and its polymorphism is also related to the follicle development in Chinese indigenous Wanxi-white goose [21]. Common polymorphisms in WNT10A have effects on the morphology of ectodermal appendages, including tooth and hair [22]. Frizzled proteins, a G protein-coupled receptor family, mainly have functions in three distinct signaling pathways (canonical Wnt/β-catenin pathway, Wnt/calcium pathway and planar cell polarity (PCP) pathway) [23]. In this part, seven of ten Frizzled proteins have expressed differently in the three stages of down feather development. Furthermore, we found that FZD1, FZD3, FZD7 and FZD8 have complex interactive relationships with other genes through gene-act-network analysis. Whether Frizzled proteins are related to the growth and development of downy feather is an interesting question that will be explored in future study. Dkk2/Frzb, an inhibitor for Wnt signaling, could regulate the feather regeneration [24]. In Fig. 3a, CTNNB1, which encode β-catenin, is located in the center of gene-act-network. In the morphogenesis of feather, β-catenin plays a significant role in the formation of feather bud [25].

Gene-co-expression
Alternatively, we performed the gene co-expression net-work analysis of DEGs in the P libraries (Fig. 4a, Table 3). There are 35 genes with highest K-core values (K-core = 34) in the gene co-expression analysis, which possibly are the key genes involved in the plumulaceous feather; of these, 30 of 35 are keratins or keratin-associated proteins (including 4 α-keratins, 2 KAPs, 1 β-keratin and 23 feather keratins), thus indicating that feather keratins play significant roles in the growth and development of plumulaceous feather. In addition, five other genes also have highest K-cores in this part, including BMP4, GJA1, PDGFC, WNT16 and BAMBI. PDGFC (platelet derived growth factor C) could enhance the feather growth and its receptor expressed in feather collar [26]. Gene co-expression analysis could find out the key genes that may regulate the growth and development of downy feather, which enriched the comprehensive understanding of plumulaceous growth from the aspect of omics.

Functional genes and molecular events accounting for the development and morphogenesis of flight feather
The flight feather is an unparalleled structure, which has two contradictory characters; flight feather consists of light material and it could withstand strong air resistance when birds fly. Flight feather is the best research material for evolutionary biology and developmental biology. In this part, same methods were used to clarify the functional genes and molecular events associated with the development and morphogenesis of flight feather.

Series-cluster analysis and functional annotation of the clusters
To better understand the significant expression pattern related to the growth and development of flight feather, same analysis was conducted based on the 5823 DEGs (F libraries) and 8 clusters were also achieved (Additional file 8: Table S3). Largest group of F libraries is cluster 1 with 1317 (22.6%) genes; the expression levels in this cluster declined gradually from stage 1 to stage 2, and maintained a relative stable expression from stage 2 to stage 3 (Fig. 2b). However, there is no GO term concerned with feather morphogenesis that was enriched in cluster 1 and 6. While cluster 5 contained the largest number of DEGs that participate in the biological process of feather growth and development, as well as the most significant GO terms; these GO terms included hair follicle development, cell adhesion, epidermis development, WNT signaling pathway, keratinocyte Fig. 3 Gene-act-network analyses of candidate genes involved in the feather development. a Gene-act-network of candidate genes from abdomen skin and follicle transcriptomes; b Gene-act-network of candidate genes from wing skin and follicle transcriptomes; Different colours of genes indicate they belong to different clusters Fig. 4 Co-expression analyses of candidate genes involved in the feather development. a Co-expression analysis of candidate genes from P libraries; b Co-expression analysis of candidate genes from F libraries. Different colours of genes indicate they belong to different clusters differentiation, BMP signaling pathway, hair follicle morphogenesis, keratinization and skin development (Fig. 2d).

GO and pathway analysis of DEGs in F libraries
The significant enriched GO terms were identified through different comparisons (See for full list of GO terms in Table 4). Intracellular signal transduction was significantly enriched in the DEGs of MF/EF, LF/MF, LF/ EF. Likewise, KEGG pathway related to flight feather development were also recognized (See for full list of KEGG pathway in Additional file 9: Table S4), including MAPK signaling pathway, VEGF signaling pathway, Jak-STAT signaling pathway, Focal adhesion, ECM-receptor interaction, TNF signaling pathway, cell adhesion molecules (CAMs), adherens junction and NF-kappa B signaling pathway.

Gene-act-network
Gene-act-network was constructed to ascertain the relations among the DEGs in F libraries (Fig. 3b). In this  [27].

Gene-co-expression
Gene co-expression net-work analysis of DEGs of F libraries was constructed (Fig. 4b, Table 5). Core genes with the highest degrees connect with most adjacent genes in the network and are frequently identified as key indicators. There are 29 genes with highest k-core values in the gene co-expression analysis; of these, 24 are keratins or keratin-associated proteins, thus indicating that feather keratins also take important roles in the growth and development of flight feather, as well as in the plumulaceous feather. In addition, five other genes also have highest K-cores in this part, including GJA1, FZD3, FZD10, PICH2 and EFNA3. Repression of Hh (Hedgehog) signaling through a dynamic PTCH1 and PTCH2 regulatory network is a crucial event in lineage fate determination in the skin [28]. FZD3 and FZD10 encode Wnt receptors, and the research has shown that FZD10 is also expressed in the feather bud [29]. In addition, FZD3 and FZD10 also have relationships with several the other genes in the gene-act-network; hence we speculate that FZD3 and FZD10 may have certain effect on the growth and development of the flight feather.

Developmental transcriptomes comparisons between plumulaceous feather and flight feather
The developmental transcriptomes of two different body parts were compared to do further research of the growth and development differences between plumulaceous feather and flight feather. Many biological processes were both enriched in the DEGs from P and F libraries. Cell adhesion molecules may regulate feather morphogenesis by straining cell motion and forming borders [30]. The studies on the mouse have shown that inhibiting Jak-STAT pathway makes hair of rapid growth [31]. ECM-receptor interaction plays a fundamental role in the morphogenesis of tissues and organs, and their main functions are the maintenance the structure of organs and functional homeostasis, and in the control of the gene expression [32]. The studies on Cashmere goat showed that high expression of ECM and cell surface proteins was essential for the rapid growth of hair follicles during the anagen phase [33]. However, in the P developmental transcriptomes, 289 GO terms are enriched in the DEGs of LPvsEP, such as SMAD protein signal transduction (GO: 0060395), which is not an enriched GO term in F libraries. SMADs are a group of signaling mediators and antagonists of the transforming growth factor-beta (TGF-beta) superfamily, as well as responding to Activin and BMPs, which play important roles in skin development. SMAD4 affects hair follicle differentiation through regulating BMP signaling. SMAD7 significantly has impact on the hair follicle development and differentiation by blocking the TGFbeta/Activin/BMP pathway and by inhibiting WNT/ beta-catenin signaling [34]. NF-kappa B signaling pathway is significantly enriched only in the F libraries. Activation of the NF-kappaB pathway by the Edar and Edaradd is required for the development of hair follilces [35]. And NF-kappa B activation is essential for induced SHH and cyclin D1 expression and subsequent hair placode down growth [36]. In addition, there are lots of DEGs between P and F libraries, especially the DEGs from BMP, WNT and FGF families. The differences of DEGs are not only the numbers and varieties, even if the same DEGs possibly have different expression trends in P and F libraries. For example, CRABP1 has been detected belonging to the profile 3 and 1 in P library and F library, respectively. The differential levels of CRABP1 and some other genes over different body regions and time would shape the anisotropic RA landscape, in order to introduce a new dimension of vane shape variations [37]. Therefore, we speculated that the existing genes construct the different gene expression networks in the skin and follicles of duck different body parts and developmental stages, so as to further differentiate into diverse feathers.

Functional genes involved in the differentiation of plumulaceous feather and flight feather
After we thoroughly investigated the DEGs of abdomen and wing skin transcriptomes at the same physiological stages, two key genes (MGP and PITX2) were screened in the two of the three comparisons of EFvsEP, MFvsMP and LFvsLP. The interaction of matrix GLA protein (MGP, an inhibitory morphogen) and BMP4 (an activating morphogen) are suggested to be important for vascular branching [38], but has not been reported to play any roles in feather morphogenesis. MGP is likely to facilitate rachis and barb branching in chicken feathers. The transcription factor PITX2, is a key factor in left-right asymmetry during the process of vertebrate development [39]. In chick embryos, early asymmetric expression of PITX2 leads to asymmetric ovarian development [40]. Another intriguing finding of this work is the identification of a total of 39 DEGs were shared across all three comparisons (Fig. 5); of these, 38 genes have consistent expression pattern in abdomen or wing skin transcriptome (Table 6). It is very interesting to find that HOXD10, HOXD11 and FMOD are all up-regulated in the EF, MF and LF. While the other 35 genes are up-regulated in the EP, MP and LP, including eight genes belonging to HOX cluster (HOXB2, HOXB3, HOXB4, HOXB5, HOXB6, HOXB8, HOXC9 and HOXC10), one transcription factor TBX4, two feather keratin (LOC101804574 and LOC106016676), as well as other genes.
In the growth and development of hair follicle, HOX genes have been confirmed to express in the mice hair follicle among different developmental periods, and multiple members of HOX genes have been considered to play vital roles in the hair follicle. The spatial expression patterns of HOXA and HOXB cluster were correlated with morphological subdivision of the digestive tract along the anteroposterior axis [41]. HOX13 is the first homeobox gene shown to have overt phenotypic effects on hair development. Abnormal expression of HOXC13 could cause the defect of hair growth to generate hairless mice [42]. HOXC8 was the only candidate gene examined with a strikingly altered expression pattern, which may directly influence the development of feathers, especially in terms of the morphology of the cranial feathers and thus cause the Crest phenotype in chicken [16]. HOXB5 regulates mesodermal-epithelial crosstalk during development, and loss of its function leads to the downregulation of previously identified downstream targets of Wnt2/2b signaling, including LEF1, AXIN2, and BMP4 [43]. In general, different HOX clusters have differentially expression levels between plumulaceous feather and flight feather at all three physiological stages, and they possibly participated in the formation of diverse feathers.

Conclusions
In conclusion, we provide a new insight into transcriptional profiles of two feather types development at three stages using a genome-wide deep sequencing method. The morphogenesis and development of avian feathers are possibly regulated by a complicated process, including a series signal transduction molecules, growth factors and transcription factors, several types of cell connections, abundant members from HOX and feather keratin families. Furthermore, this study of feather transcriptional profiles is helpful in understanding the differences between the genetic mechanisms in plumulaceous and flight feathers, and for future work exploring genetic basis of feather divergence and development.

Ethics statement
This study was carried out in strict accordance with the recommendations of the Regulations for the  were given continuous access to a standard commercial feed ration and water. All surgery was performed under combination anesthesia, and all efforts were made to minimize suffering of animals.
Sample rearing, sampling of and RNA extraction from skin tissues The samples of Cherry Valley ducks were selected from Sichuan Xinmianying Farming Corporation. And we maintained all of these individuals in the same raising environment (fed with pellet feed and green vegetables) at the Institute of Zoology, Shaanxi Normal University, Xi'an, Shaanxi Province, China, in 2014. Male ducks at three different physiological stages (3-day-old, 27-day-old and 6-month-old) were selected and divided into early growth phase (E), middle growth phase (M) and late growth phase (L) groups. First, the feathers born in the abdomen and inner side of the wing at each stage were sheared and further shaved (we have not pulled out the whole feathers). Then the ducks were locally anaesthetized with pentobarbital sodium (3%, 1 ml/kg) through intraperitoneal injection to minimize the animal suffering. A piece of skin (1 cm in diameter) owing feather follicle from the abdominal and wing side were collected with scissors and forceps and immediately washed with 10 mL PBS (PH7.2) and 0.5 mM EDTA. Each sample was then placed in RNAlater (Ambion) and stored at 4°C for 24 h first and then at − 80°C until RNA extraction. RNAlysis was performed using TRIzol Reagent (Invitrogen, USA) followed by isolation and purification with the RNeasy kit (Qiagen, Germany). RNA degradation and contamination was assessed on 1% agarose gel electrophoresis. RNA purity and integrity were evaluated using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), respectively.

Library preparation and sequencing
Twelve samples with RNA integrity number (RIN) values above 7 were used for libraries construction. Sequencing libraries were generated using the IlluminaTru-Seq™ RNA Sample Preparation Kit (Illumina, San Diego, USA) following the manufacturer's recommendations. Briefly, Oligo(dT) magnetic beads were used to isolate Poly(A) mRNA from total RNA. Then purified mRNA was interrupted into fragments with divalent cations. Double-stranded cDNA was synthesized using random hexamer-primers, reverse transcriptase, DNA polymerase I and RNaseH by taking short fragments as templates. Subsequently, double-stranded cDNA was further subjected to end-repair and ligation with adapters. These modified products were purified and enriched with PCR to construct the final cDNA library. After test the quality of libraries, they were loaded onto the flow cell channels of an Illumina HiSeq™ 2000 platform and 90 bp pair-end reads were generated at BGI, Shenzhen, China.

Sequence reads mapping and assembly
Be assembly, raw reads of fastq format were firstly filtered by removing reads containing adaptors, reads containing poly-N and low quality reads to obtain high-quality clean reads under following criteria: a. Reads contain 20% base quality lower than Q20; b. Reads length ≥ 50 bp; c. Trim N-end. At the same time, quality parameters of clean data including Q20 and GC-content were obtained. All the succeeding analyses were carried out using high quality clean reads. Clean reads were mapped to the Peking duck genome assembly BGI_-duck_1.0 (GCA_000355885.1) by using TopHat2 [44] with following parameter (tophat -r 272 -a 10 -m 0 -i 31 -I 500000 -read-mismatches 3 -max-insertion-length 6 -max-deletion-length 6 -read-gap-length 7 -read-edit-dist 12 -b2-sensitive -solexa1.3-quals -p 8 -min-coverage-intron 31 -min-segment-intron 31).

Quantification and identification of differentially expressed genes (DEGs)
The parameter FPKM (Fragments per kilobase of transcript per million mapped reads) was applied to quantify the gene expression levels by using BGI_duck_1.0 genome annotation gtf file. Gene expression calculated through FPKM could eliminate the influence of gene length and sequencing data on gene expression. HTseq [45] was used for count calculation and FPKM was calculated with NCBI gtf file through gene length annotation by domestic code. We used the calculated gene expression to compare the differences in gene expression between the skin transcriptome of abdomen and wing during three developmental stages. The DEGseq package was applied to filter the differentially expressed genes with a fold change > 2 or fold change < 0.5, and false discovery rate (FDR) < 0.05 [46,47].

Classification, annotation and co-expression of DEGs
Series cluster analysis was performed using STEM [48] to classify the differentially expressed genes in eight clusters based on the FPKM change tendency of genes in three developmental stages. Fisher's exact test and the multiple comparison tests were used to calculate the significant levels of profiles [49,50].
Gene ontology (GO) analysis was performed to facilitate elucidating the biological implications of unique genes in the significant or representative profiles, which helps us to find those GOs with more concrete function description in this study [51]. Generally, Fisher's exact test was applied to idenfity the significant GO categories