- Research article
- Open Access
Dynamic transcriptome profiling towards understanding the morphogenesis and development of diverse feather in domestic duck
BMC Genomics volume 19, Article number: 391 (2018)
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.
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.
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.
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 . The prosperous feather diversities have contributed significantly to the rapid and extensive radiation of birds to become the dominant vertebrate . 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 .
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 . 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 .
The morphogenesis of feather is initiated by the interplay of epithelia and subjacent mesenchayme and usually involves a series of dynamic cellular processes . 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 . 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 . 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 . Perturbing the gradient of WNT3A converts bilaterally symmetric feathers into radially symmetrical feathers . 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 . Feathers consist mainly of flexible corneous materials made of α- and β-keratin multigene families . Ng et al.  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 .
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 .
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.
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 R2 ≥ 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 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 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). Keratinization, extracellular matrix organization, cell adhesion, hair follicle development, keratinocyte differentiation and hair follicle morphogenesis were all significantly enriched in the DEGs of MP/EP, LP/MP and LP/EP. Different KEGG pathways accounting for plumulaceous feather development were determined in 3 different comparison libraries (See for full list of KEGG pathway in Additional file 6: 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) . 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 .
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 . 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 . Common polymorphisms in WNT10A have effects on the morphology of ectodermal appendages, including tooth and hair . 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) . 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 . 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 .
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 . 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 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 was constructed to ascertain the relations among the DEGs in F libraries (Fig. 3b). In this gene interaction network, we could find that few DEGs belong to BMP families were screened in F libraries, including BMP2, BMP4, BMPR2 and BMPR1A. BMP2 is expressed in the posterior region of the feather buds. BMP and SHH pathways are involved in regulating the formation and balance among the rachis and barbs. A wide variety of WNT family (WNT7A, WNT10A, WNT11, WNT16, FZD2), and its receptors of frizzled protein (FZD2, FZD3, FZD4, FZD10) and DKK protein (DKK1, DKK2), as well as DEGs belong to FGF family (FGF1, FGF2, FGFR1, FGFR3, FGF9, FGF12, FGF13, FGF14, FGF18) were determined in the F libraries. Some particular DEGs (such as HHIP and DVL1) were identified only in F libraries comparing with P libraries. HHIP interact with hedgehog family to regulate the hedgehog signaling of several cell types .
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 . FZD3 and FZD10 encode Wnt receptors, and the research has shown that FZD10 is also expressed in the feather bud . 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 . The studies on the mouse have shown that inhibiting Jak-STAT pathway makes hair of rapid growth . 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 . 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 .
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 . 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 . And NF-kappa B activation is essential for induced SHH and cyclin D1 expression and subsequent hair placode down growth . 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 . 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 , 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 . In chick embryos, early asymmetric expression of PITX2 leads to asymmetric ovarian development . 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 . 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 . 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 . 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 . 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.
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.
This study was carried out in strict accordance with the recommendations of the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004). Housing and caring of Cherry Valley ducks and collection of skin samples for use in the described experiments were conducted following the approved protocol of the Institutional Animal Care and Use Committees of college of Life Sciences, Shaanxi Normal University. All ducks in this study 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 IlluminaTruSeq™ 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  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  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  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 . Generally, Fisher’s exact test was applied to idenfity the significant GO categories and FDR was used to correct the p-values. The threshod of significance was defined by p-value < 0.05 and FDR < 0.05. Within the significant category, the enrichment was given by: Enrichment = (nf/n)/(Nf/N), where “N” was the number of BG genes (background-genes) which were achieved from the genome gtf file download from NCBI, “n” was the number of total DEGs or the genes we analyzed with GO annotation, “nf” and “Nf” represent the DEGs and BG genes with the target GO-term annotation .
Pathway analysis was used to find out the significant pathway of the differential expressed genes according to the KEGG database. A Fisher exact test was used to find the significant enrichment pathway with the threshod of significance of p-value < 0.05 and FDR < 0.05, which helps us to find those more significant pathways in this study. The calculation method of enrichment value is similar with the GO enrichment.
Pathway-act-network analysis was performed to reveal the interactive network among the pathways with enriched DEGs based on the KEGG database, while graphical representations of the pathways were generated based on the Cytoscape software . The core pathway which is related with the morphogenesis and development of feather could be identified by building the signaling pathway network.
Gene-act-net analysis was conducted to reveal the network of the DEGs based on the interactions among the genes, proteins and compounds included in the KEGG database. In the network, cycle nodes represent genes, and edges between two nodes represent interactions between genes.
Gene co-expression network analysis was performed to track the interactions among the DEGs, according to their dynamic expression changes in three developmental stages. Pearson correlation was applied to each pair of genes and the significantly correlated pairs were used to construct the network . K-core scoring of a given gene, which indicates its hub or nodal status with connection to “k” other genes in a network, was introduced to locate the core regulatory genes in the network. Accordingly, the genes with largest k-core scores and highest degrees of connection were identified as “key regulatory genes” in a network [54, 55].
Verification by quantitative real-time PCR
To assess the reliability of our sequencing and analysis by real-time quantitative PCR (RT-qPCR), we used the same RNA samples for transcriptome sequencing. For first-strand cDNA synthesis, 1 μg of total RNA was reverse-transcribed using the Transcriptor First Strand cDNA Synthesis Kit (Roche) according to the manufacturer’s protocol. In the subsequent experiments, the cDNA samples were diluted 10 times prior to qPCR. Real-time qPCR was performed using FastStart Essential DNA Green Master (Roche) on the CFX96 Real-Time System (BioRAD, USA). The reaction was performed using the following conditions: denaturation at 95 °C for 10 min, followed by 40 cycles of amplification (95 °C for 10 s, 60 °C for 15 s and 72 °C for 20 s). Relative expression was calculated by the 2-ΔΔCt method using β-actin as the reference control . All primers are listed in Additional file 10: Table S5.
Bone morphogenesis protein
Cell adhesion molecules
Differentially expressed genes
Early growth flight feathers
Early growth plumulaceous feathers
Fibroblast growth factor
Fragments per kilobase of transcript per million mapped reads
Janus kinase and the signal transducer and activator of transcription
Kyto encyclopaedia of genes and genomes
Late growth flight feathers
Late growth plumulaceous feathers
Mitogen-activated protein kinase
Middle growth flight feathers
Middle growth plumulaceous feathers
Platelet-derived growth factor C
Real-time quantitative PCR
Transforming growth factor-beta; NF-kappaB
Vascular endothelial growth factor
Chuong CM, Randall VA, Widelitz RB, Wu P, Jiang TX. Physiological regeneration of skin appendages and implications for regeneration medicine. Physiology (Bethesda). 2012;27(2):61–72.
Chen CW, Chuong CM. Avian integument provides multiple possibilities to analyse different phases of skin appendage morphogenesis. J Investig Dermatol Symp Proc. 1999;4(3):333–7.
Feng C, Gao Y, Dorshorst B, Song C, Gu X, Li J, et al. A cis-regulatory mutation of PDSS2 causes silky-feather in chickens. PLoS Genet. 2014;10(8):e1004576.
Kowata K, Nakaoka M, Nishio K, Fukao A, Satoh A, Ogoshi M, et al. Identification of a feather β-keratin gene exclusively expressed in pennaceous barbules cells of contour feathers in chicken. Gene. 2014;542(1):23–8.
Alibardi L. Cell structure of barb ridges in down feather and juvenile wing feathers of the developing chick embryo: barb ridge modification in relation to feather evolution. Ann Anat. 2006;188(4):303–18.
Chuong CM, Edelman GM. Expression of cell-adhesion molecules in embryonic induction. II. Morphogenesis of adult feathers. J Cell Biol. 1985;101(3):1027–43.
Widelitz RB, Jiang TX, Yu M, Shen T, Shen JY, Wu P, et al. Molecular biology of feather morphogenesis: a testable model for evo-devo research. J Exp Zool B Mol Dev Evol. 2003;298((1):109–22.
Jiang TX, Tuan TL, Wu P, Widelitz RB, Chuong CM. From buds to follicles: matrix metalloproteinases in developmental tissue remodeling during feather morphogenesis. Differentiation. 2011;81(5):307–14.
Lin CM, Jiang TX, Widelitz RB, Chuong CM. Molecular signaling in feather morphogenesis. Curr Opin Cell Biol. 2006;18(6):730–41.
Chuong CM, Patel N, Lin J, Jung HS, Widelitz RB. Sonic hedgehog signaling pathway in vertebrate epithelial appendage morphogenesis: perspectives in development and evolution. Cell Mol Life Sci. 2000;57(12):1672–81.
Harris MP, Williamson S, Fallon JF, Meinhardt H, Prum RO. Molecular evidence for an activator-inhibitor mechanism in development of embryonic feather branching. PNAS. 2005;102(33):11734–9.
Yue Z, Jiang TX, Widelitz RB, Chuong CM. Wnt3a gradient converts radial to bilateral feather symmetry via topological arrangement of epithelia. PNAS. 2006;103(4):951–5.
Alibardi L. Gap and tight junctions in the formation of feather branches: a descriptive ultrastructural study. Ann Anat. 2010;192(4):251–8.
Alibardi L, Toni M. Cytochemical and molecular characteristics of the process of cornification during feather morphogenesis. Prog Histochem Cytochem. 2008;43(1):1–69.
Ng CS, Wu P, Fan WL, Yan J, Chen CK, Lai YT, et al. Genomic organization, transcriptomic analysis, and functional characterization of avian α- and β-keratins in diverse feather forms. Genome Biol Evol. 2014;6(9):2258–73.
Wang Y, Gao Y, Imsland F, Gu X, Feng C, Liu R. The crest phenotype in chicken is associated with ectopic expression of HOXC8 in cranial skin. PLoS One. 2012;7(4):e34012.
Morozoca O, Marra MA. Applications of next-generation sequencing technologies in functional genomics. Genomics. 2008;92(5):255–64.
Samuelov L, Sprecher E, Ralf P. The role of P-caderin in skin biology and skin pathology: lessons from the hair follicle. Cell Tissue Res. 2015;360(3):761–71.
Rishikaysh P, Dev K, Diaz D, Qureshi WM, Filip S, Mokry J. Signaling involved in hair follicle morphogenesis and development. Int J Mol Sci. 2014;15(1):1647–70.
Reddy S, Andl T, Bagasra A, Lu MM, Epstein DJ, Morrisey EE, et al. Characterization of Wnt gene expression in developing and postnatal hair follicles and identification of Wnt5a as a target of sonic hedgehog in hair follicle morphogenesis. Mech Dev. 2001;107(1–2):69–82.
Chen X, Bai H, Li L, Zhang W, Jiang R, Geng Z. Follicle characteristics and follicle developmental related Wnt6 polymorphism in Chinese indigenous Wanxi-white goose. Mol Biol Rep. 2012;39(11):9843–8.
Kimura R, Watanabe C, Kawaguchi A, Kim YI, Park SB, Maki K, et al. Common polymorphisms in WNT10A affect tooth morphology as well as hair shape. Hum Mol Genet. 2015;24(9):2673–80.
Huang HC, Klein PS. The frizzled family: receptors for multiple signal transduction pathways. Genome Biol. 2004;5(7):234.
Chu Q, Cai L, Fu Y, Chen X, Yan Z, Lin X, et al. Dkk2/Frzb in the dermal papillae regulates feather regeneration. Dev Biol. 2014;387(2):167–78.
Noramly S, Freeman A, Morgan BA. β-catenin signaling can initiate feather bud development. Development. 1999;126(16):3509–21.
Chuong CM, Chen HM, Jiang TX, Chia J. Adhesion molecules in skin development: morphogenesis of feather and hair. Ann N Y Acad Sci. 1991;642:263–80.
Bosanac I, Maun HR, Scales SJ, Wen X, Lingel A, Bazan JF, et al. The structure of SHH in complex with HHIP reveals a recognition role for the Shh pseudo active site in signaling. Nat Struct Mol Biol. 2009;16(7):691–7.
Adolphe C, Nieuwenhuis E, Villani R, Li ZJ, Kaur P, Hui CC, et al. Patched 1 and patched 2 redundancy has a key role in regulating epidermal differentiation. J Invest Dermatol. 2014;134(7):1981–90.
Kawakami Y, Wada N, Nishimatsu S, Komaguchi C, Noji S, Nohno T. Identification of chick frizzled-10 expressed in the developing limb and central nervous system. Mech Dev. 2000;91(1–2):375–8.
Ng CS, Chen CK, Fan WL, Wu P, Wu SM, Chen JJ, et al. Transcriptomic analyses of regenerating adult feathers in chicken. BMC Genomics. 2015;16:756.
Harel S, Higgins CA, Cerise JE, Dai Z, Chen JC, Clyes R, et al. Pharmacologic inhibition of JAK-STAT signaling promotes hair growth. Sci Adv. 2015;1(9):e1500973.
Zhu B, Xu T, Yuan J, Guo X, Liu D. Transcriptome sequencing reveals differences between primary and sencondary hair follicle-derived dermal papilla cells of the cashmere goat (Capra hircus). PLoS One. 2013;8(9):e76282.
Xu T, Guo X, Wang H, Hao F, Du X, Gao X, et al. Differential gene expression analysis between anagen and telogen of Capra hirus skin based on the de novo assembled transcriptome sequence. Gene. 2013;520(1):30–8.
Owens P, Han G, Li AG, Wang XJ. The role of SMADs in skin development. J Invest Dermatol. 2008;128(4):783–90.
Morlon A, Munnich A, Smahi A. TAB2, TRAF6 and TAK1 are involved in NF-kappaB activation induced by the TNF-receptor, Edar and its adaptator Edaradd. Hum Mol Genet. 2005;14(23):3751–7.
Schmidt-Ullrich R, Tobin DJ, Lenhard D, Schneider P, Paus R, Scheidereit C. NF-kappaB transmits Eda A1/EdaR signaling to activate Shh and cyclin D1 expression, and controls post-initiation hair placode down growth. Development. 2006;133(6):1045–57.
Li A, Figueroa S, Jiang TX, Wu P, Widelitz R, Nie Q, et al. Diverse feather shape evolution enabled by coupling anisotropic signaling modules with self-organizing branching programme. Nat Commun. 2017;8:14139.
Yao Y, Nowak S, Yochelis A, Garfinkel A, Boström KI. Matrix GLA protein, an inhibitory morphogen in pulmonary vascular development. J Biol Chem. 2007;282(41):30131–42.
Torlopp A, Khan MA, Oliveira NM, Lekk I, Soto-Jiménez LM, Sosinsky A, et al. The transcription factor pitx3 positions the embryonic axis and regulates twinning. elife. 2014;3:e03743.
Ishimaru Y, Komatsu T, Kasahara M, Katoh-Fukui Y, Ogawa H, Toyama Y, et al. Mechanism of asymmetric ovarian development in chick embryos. Development. 2008;135(4):677–85.
Sakiyama J, Yokouchi Y, Kuroiwa A. HoxA and HoxB cluster genes subdivide the digestive tract into morphological domains during chick development. Mech Dev. 2001;101(1–2):233–6.
Godwin AR, Capecchi MR. Hoxc13 mutant mice lack external hair. Genes Dev. 1998;12(1):11–20.
Hrycaj SM, Dye BR, Baker NC, Larsen BM, Burke AC, Spence JR, et al. Hox5 genes regulate the Wnt2/2b/Bmp4-signaling axis during lung development. Cell Rep. 2015;12(6):903–12.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.
Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.
Dupuy D, Bertin N, Hidalgo CA, Venkatesan K, Tu D, Lee D. Genome-scale analysis of in vivo spatiotemporal promoter activity in Caenorhabditis elegans. Nat Biotechnol. 2007;25(6):663–8.
Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7:191.
Miller LD, Long PM, Wong L, Mukherjee S, McShane LM, Liu ET. Optimal gene expression analysis by microarrays. Cancer Cell. 2002;2(5):353–61.
Ramoni MF, Sebastiani P, Kohane IS. Cluster analysis of gene expression dynamics. PNAS. 2002;99(14):9121–6.
Gene Ontology Consortium. The gene ontology (GO) project in 2006. Nucleic Acids Res. 2006;34(Database issue):D322–6.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Prieto C, Risueño A, Fontanillo C, De las Rivas J. Human gene coexpression landscape: confident network derived from tissue transcriptomic profiles. PLoS One. 2008;3(12):e3911.
Barabási AL, Oltvai ZN. Network biology: understanding the cell’s functional organizaiton. Nat Rev Genet. 2004;5(2):101–13.
Ravasz E, Somera AL, Mongru DA, Oltvai ZN, Barabási AL. Hierarchical organizaiton of modularity in metabolic networks. Science. 2002;297(5586):1551–5.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25(4):402–8.
We thank Fei Ye, Miao Sha, Fei Zhou and Yang Song (Shaanxi Normal University) for assistance with the sample collections. We are also grateful for the BGI and Novelbio for assistance with the transcriptome sequencing and data analyzing.
This research was supported by the Chinese Academy of Sciences (KSZD-EW-Z-005-003), State Key Program of NSFC (31330073) and the Key Laboratory of the Zoological Systematics and Evolution of the Chinese Academy of Sciences (No. O529YX5105) to Fumin Lei.
Availability of data and materials
Our RNA high-throughput sequencing data are available at the GEO website under accession number GSE104386 [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104386].
Ethics approval and consent to participate
This study was approved by the Animal Care and Use Committees of college of Life Sciences, Shaanxi Normal University. Housing and collecting of skin samples involving duck in the described experiments were carried out in strict accordance with the recommendations of the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004). All animals used in this study were handled in strict accordance with good clinical practices and all efforts were made to minimize suffering.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Gene structure analyses of the mapped reads on different regions of the reference duck genome. (PNG 402 kb)
Figure S2. Analysis of the gene expression correlation between the two biological replicates of six groups (EF, EP, MF, MP, LF and LP). (PNG 280 kb)
Figure S3. Regression analysis of gene expression fold changes (FC) obtained from quantitative PCR and RNA-Seq. (PNG 35 kb)
Figure S4. Venn diagram of the DEGs of P, F libraries. A: Venn diagram of the DEGs of P libraries (including EP/MP, MP/LP and EP/LP); B Venn diagram of the DEGs of F libraries (including EF/MF, MF/LF and EF/LF). (PNG 329 kb)
Table S1. The detailed information regarding the identity of the expression profiles of DEGs from MPvsEP, LPvsMP and LPvsEP. (XLS 1357 kb)
Table S2. KEGG pathway functional enrichment analysis related with plumulaceous feather development of the DEGs of MP/EP, LP/MP and LP/EP. (DOCX 16 kb)
Figure S5. Pathway-act-network analysis of the plumulaceous development in Cherry Valley duck. (PNG 228 kb)
Table S3. The detailed information regarding the identity of the expression profiles of DEGs from MFvsEF, LFvsMF and LFvsEF. (XLS 1656 kb)
Table S4. KEGG pathway functional enrichment analysis related with flight feather development of the DEGs of MF/EF, LF/MF and LF/EF. (DOCX 16 kb)
Table S5. Primers used for qPCR validation. (DOCX 12 kb)