The fertilized eggs of mallards were collected from the waterfowl breeding farm of Sichuan Agricultural University, Ya’an. After incubation, the hatched mallards were maintained for 14 weeks in a comfortable environment. These ducks could freely feed and drink throughout their growth stages. The head and back skin tissue of drakes and the head skin tissue of females were sampled at the ages of 7, 11, and 14 weeks. We considered the difficulty of sampling feather follicles, the feather follicles to stay in the skin, and tested the skin tissue directly. An equal amount of skin tissue in each sample and feather follicles remained in each skin sample. In addition, feather samples were collected from 6- to 7-week-old and 14-week-old mallards. After sample collection was completed, some of the skin tissue samples were stored at −80°C, and the remaining skin tissue samples were stored in 4% paraformaldehyde, while the feather samples were stored in 3% glutaraldehyde.

The skin samples were dehydrated with EthOH in successive concentrations of 70, 80, 90, 95, 99, and 100%. The samples were further infiltrated three times in xylene. Each dehydration step was performed in a dehydration system (Leica TM 1020, Germany) for 22 h, and the tissues were then paraffin-embedded. We investigated skin anatomy in cross-sections and vertical sections. The skin tissue samples were cut using microtome blades on a Leica RM 2135 microtome. The sections were stained with L-Dopa (Solarbio, China). Images of the tissue sections were obtained on an Olympus SZX16 stereomicroscope (Japan). Finally, the melanin content of the barb ridges was analyzed using IPP 6.0 (Image-Pro Plus 6.0). The melanin content was expressed by calculating the area occupied by the melanin, which was stained dopa in the cross-section of the feather follicles.

Barbule anatomy was investigated by transmission electron microscopy (TEM) using standard methods (Thorpe et al., 2008). Briefly, a feather piece was subjected to a series of dehydration in acetone at increasing concentrations of 30, 50, 70, 80, 90, 95, and 100% and then infiltrated with epoxy resin at 3:1, 1:1, and 1:3 proportions. Next, we cured the resin blocks by heating them and cut them into approximately 50 nm sections using a Leica EM UC7 ultramicrotome (Leica Microsystems GmbH, Germany). The sections were stained using uranyl acetate and lead citrate. Next, we checked the sections on a JEM-1400PLUS TEM (Japan). Finally, the number of melanosomes was analyzed by using IPP 6.0.

Total RNA was extracted from duck skin tissues using TRIzol (RNAiso Plus, Takara, Japan). Illumina platform sequencing of all samples was conducted at Novogene Bioinformatics Technology Co., Ltd., Beijing. In total, 15 sample libraries were constructed, including 12 libraries from head skin samples (three repeats from each stage of drake sampling in the 7^th, 11th^, and 14th weeks and three repeats from female mallards sampled in the 14th week) and three from dorsum skin samples collected from drakes in the 14th week. Low-quality reads (i.e., tags containing only adaptors and ambiguous bases) were removed from all samples.

HISAT 2.0 software (Kim et al., 2015) was used to map the sequencing data with the assembly of the mallard reference genome (IASCAAS_PekingDuck_PBH1.5, GCF_003850225.1) to generate read count results. Second, we used STEM (Short time-series expression miner) software to analyze the genome-wide expression trends in the partial FPKM results, including three growth stages (7th weeks, 11th weeks, and 14th weeks) of drakes. Third, differential gene expression was tested using DESeq 2.0 from the R package (Love et al., 2014). The log₂ (fold-change) values of the genes were tested using DESeq to evaluate the expression patterns based on all sequencing data. Fourth, the gene ontology (GO) pathways of differentially expressed genes were analyzed using the R package’s cluster profile (Yu et al., 2012). Finally, we used GraphPad 8.0, BioVenn, and MeV 4.9 software to generate the volcano plots, Venn diagram, and heat map of gene expression, respectively.

The online transcription factor (TF) binding site prediction tool (http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/tfbs_predict) from the Animal TFBD database was used to predict the TF binding sites in the promoter regions of the TYR (Chr1- 8,264,483∼82,694,830 bp) and TYRP1 genes (ChrZ- 307,517,598∼30,752,598 bp).

In the synteny analysis, we analyzed six avian assemblies from NCBI (National Center for Biotechnology Information), including assemblies for mallard (IASCAAS_PekingDuck_PBH1.5), Anna’s hummingbird (bCalAnn1_v1.p, GCF_003957555.1), chicken (GRCg6a, GCF_000002315.6), Japanese quail (Coturnix japonica 2.1, GCF_001577835.2), turkey (Turkey 5.1, GCF_000146,605.3), and zebra finch (bTaeGut1_v1.p, GCF_003957565.1). These six species are representative avian species, and their genomes have been well annotated.

We obtained the coding regions (CDs) of five genes (TYRP1, PTPRD, LURAP1L, MPDZ, and NFIB) of 6 avian species from the NCBI database (Supplementary Table S1). We used Clustal Omega of MEGA 7.0 (Sudhir et al., 2016) to perform multisequence alignment of these sequences. Revised sequence alignments were submitted to MEGA 7.0 to construct a bootstrap (1000 replicate) tree (Felsenstein 1985) for each gene. The maximum likelihood (ML) method was used for phylogenetic analysis. We identified the best DNA models for each gene (TYRP1: K2 + G, PTPRD: T92 + G, LURAP1L: TN93 + G, MPDZ: K2 + I, and NFIB: K2). The ML search started with the tree generated by BioNJ (Make initial tree automatically) (Gascuel 1997), and the optimal tree was determined via a heuristic search using the NNI (nearest-neighbor interchange) algorithm.

We performed all statistical analyses using R packages (R Core Team, 2014). A p-value < 0.01 was considered highly significant, while a p-value < 0.05 was considered statistically significant.

The villi on the heads of the mallard ducklings began to transform into contour feathers at the age of 5 or 6 weeks. Initially, new feathers grew in areas near the eyes of drakes (Figure 1B). Then, at the 6th to 14th weeks, new feathers gradually developed in the head and neck (Figure 1B). At the 14th week, the feathers finally turned green and covered the drakes’ heads entirely (Figure 1B). Throughout the growth period of female mallards, the feathers in most areas were light in color (Figure 1B) and did not change further.

We performed histological observations to obtain a better understanding of green feather development. In short, feathers develop from the feather follicles. The feather follicles consist of the dermal papilla and epidermal collar, and the proliferation of the epithelial cells is active at the bottom of the feather follicles. During the development of the feathers, epithelial cells differentiate to the barb ridges at the top of the feather follicles firstly. And the barb ridges further proliferates and differentiates into the marginal plates, barbule plates, and axial plates, and finally form the mature feathers. During this process, mature melanocytes located in keratinocytes release the melanosomes. And the melanosomes transfer from the bottom of the feather follicles to the small barbules (Yu et al., 2004). In the barb ridges of mallard drakes (Figure 1D), melanin began to deposit at the age of 7 weeks. The melanin content increased significantly at the ages of 11th and 14th weeks (Figures 1C,D). In female ducks, melanin deposition was not observed at 14 weeks (Figures 1C,D).

According to the observations of the vertical sections of the follicle pulp of males (Figure 1D), melanin mainly transferred from the dermal papilla to the barb ridges. At the age of 7 weeks, intense melanin deposition was observed in the follicle pulp. At the age of 14 weeks, less melanin was observed in follicle pulp areas. However, there was still less melanin in females’ follicle pulp areas (Figure 1D).

In the feathers, only the barbules at the top showed black coloration at the duck age of 5 weeks (Figure 1F), which gradually changed from black to green at the ages of 6–14 weeks. We further investigated cross-sections of the small barbules between the black parts of feathers in the 6th∼7th weeks and the green parts of feathers in the 14th week (Figures 1E,F). There were few melanosomes in the barbules of the black regions, and they were irregularly arranged under the keratin layer. In contrast, in the barbules of the green parts, there were more melanosomes located under and enwrapped by the keratin layer, and they were arranged in an orderly pattern (4 layers) (Figures 1E,F). Taken together, we suggested that the development of green feathers may be caused by the abundant melanin transferred to the small barbules. We further assumed that the excessive deposition of melanin/melanosomes caused the orderly arrangement of melanosomes and eventually led to green feather development through coherent emitted light scattering.

A summary of the sequencing reads and matched genes is shown in Supplementary Table S2. After the low-quality reads were removed, an average of 5,000,000 clean reads remained for each sample. In addition, an average of 93.2% reads was mapped to the reference genome of mallard (IASCAAS_PekingDuck_PBH1.5, GCF_003850225.1).

To understand the genes involved in developing green feather color on the drake head, we performed gene expression profile cluster analysis based on the transcriptome data obtained at the 7th week, 11th week, and 14th week during green head feather development. The results showed that all genes could be divided into 16 categories (0–15), among which five significant types were found, including 14, 10, 11, 15, and 9. Furthermore, among the five categories, the expression patterns of categories 14 and 10 and those of categories 11 and 15 were similar (Figure 2A).

In these five significant categories, the expressions of all genes were upregulated in the 7th–11th weeks and downregulated in the 11th–14th weeks. In category 10, there were 16 genes whose functions were related to melanocyte development (EDNRB, KIT, LYST, SOX10), melanin production (POMC, OCA2, TYR, TYRP1, TRPM1, SLC45A2, SLC24A4, PAH), and the melanin distribution (MYO5A, RAB27B, MLANA) (Jackson 1994; Bellone 2010; Cieslak et al., 2011; Sturm and Duffy, 2012). Two other genes in category nine were related to melanin production (MC1R, RAF1) (Jackson 1994; Bellone 2010; Cieslak et al., 2011; Sturm and Duffy, 2012), while three genes in category 14 were related to melanin production (CLCN7, HPGDS, STX17) (Jackson 1994; Bellone 2010; Cieslak et al., 2011; Sturm and Duffy, 2012). GO analysis further showed that (Figure 2B; Supplementary Table S3) signaling pathways associated with BMP and Wnt were significantly enriched. These signaling pathways were reported to play roles in melanin biosynthesis (Huang et al., 2011) and feather development (Rouzankina et al., 2004). Hence, the transcriptomic analysis revealed that most pigmentation genes were highly expressed in feather follicles during green head feathers, contributing to melanin deposition. Meanwhile, the results also showed that 7–11 weeks might be a more rapid period for melanin biosynthesis than 11–14 weeks in the males’ heads.

To screen causative genes, we screened the differentially expressed genes (DEGs) among three comparisons: male head vs. female head, male head vs. male back, and male head in the 7th week vs. male head in the 11th week. The results showed that 1147, 297, and 177 DEGs were identified (Figure 3A). GO analysis of the DEGs further revealed that the TYR and TYRP1 genes were enriched in the pigment biosynthesis process, pigmentation, pigment metabolic, melanosome, and pigment granule pathways (Supplementary Tables S4–S6).

The candidate genes responsible for the male green head feathers should explain the differences between all green and nongreen feather samples; thus, we checked the consensus DEGs among the three comparison pairs. Venn analysis identified 18 genes as potential causative genes determining green male head feathers (Figure 3B,C). Among these genes, only TYR and TYRP1 were involved in melanin biosynthesis, which are critical downstream genes of the three melanin synthesis pathways (cAMP, Wnt, and MAPK pathways) (Figure 3C). In contrast, others were related to receptor proteins, molecular binding, and cellular structure. In the three comparison groups, the expression levels of TYR, SLC38A11, REELD1, SYNPR, LOC106019585, and TYRP1 were all up-regulated in green samples, while only the LOC110352708 was down-regulated (Figure 3D). In comparing males vs. females, the TYRP1 and TYR gene expression were 256 and 32 times higher in the head feather follicles of males than in females’ head feather follicles and the males’ back feathers, respectively. Among the 18 candidate genes, TYRP1 was located on the Z-chromosome, and the other genes (including TYR) were located on euchromosomes.

Finally, as males in birds have two Z-chromosomes while females only have one, we compared the gene expression patterns of the head feather follicles of male and female mallards. First, we found that the expression levels of most Z-chromosome-linked genes (Accounted for about 63% of the expressed genes) were higher in males than in females (Figure 4A; Table S11). Then, we analyzed the average log₂ fold-change values of all genes on each chromosome, and we found that the Z-chromosome genes exhibited higher log₂ fold-changes than those on all euchromosomes (Figure 4B). Therefore, we believed it is one reason that the dosage effect of the Z-chromosome in males may lead to higher expression of the TYRP1 gene in male head feather follicles than in females.

The transcription factors (TFs) binding sites of the duck TYR and TYRP1 genes were predicted using online prediction software with the promoters as input. All predicted TFs are provided in Supplementary Tables S7, S8. Binding sites of MAFA (MAF bZIP transcription factor A), GSC (goosecoid homeobox), OTX1 (orthodenticle homeobox 1), and FOXF2 (forkhead box F2) TFs were predicted in the TYR promoter region. In contrast, binding sites of MAFA, ARX (aristaless related homeobox), and FOXF2 were predicted in the TYRP1 promoter region (Figure 5A). Among these TFs coding genes, MAFA and ARX are members of 18 potential causative genes in the three comparison groups At the same time, GSC, OTX1, and FOXF2 were differentially expressed only between the head follicles and back follicles of drakes (Supplementary Table S9). Based on gene expression cluster analysis (Figure 5B) of the drake samples, we observed that ARX and MAFA expression patterns were most similar to those of the TYRP1 gene. Thus, the differential expression of TYRP1 between different body parts and time points in males was possibly due to differences in the cis-regulation of potential transcription factors.

Sexual dimorphism in plumage color is common in avian species. Therefore, we further checked the pigmentation genes on the Z-chromosomes of six avian species showing sexual dimorphism of plumage color characteristics. The genomes of these six avian species were well assembled and annotated. Three pigmentation genes, MLANA, SLC45A2, and TYRP1, are located on the Z-chromosome (Supplementary Table S10). Synteny analysis (Figure 6) showed that the TYRP1 gene and its neighboring genes were positioned according to their relative locations on the Z-chromosome. Four genes adjacent to TYRP1, including PTPRD, LURAP1L, MPDZ, and NFIB, were conserved in the six avian species and shared a similar transcription direction. Phylogenetic analysis (Figure 6) further suggested that the above five genes presented similar phylogenetic relationships among the six avian species. These findings implied that the chromosomal segment carrying the TYRP1 gene and its neighboring genes were conserved in avian species.