Transition from natal downs to juvenile feathers: conserved regulatory switches in Neoaves

The transition from natal downs for heat conservation to juvenile feathers for simple flight is a remarkable environmental adaptation process in avian evolution. However, the underlying epigenetic mechanism for this primary feather transition is mostly unknown. Here we conducted time-ordered gene co-expression network construction, epigenetic analysis, and functional perturbations in developing feather follicles to elucidate four downy-juvenile feather transition events. We discovered that LEF1 works as a key hub of Wnt signaling to build rachis and converts radial downy to bilateral symmetry. Extracellular matrix reorganization leads to peripheral pulp formation, which mediates epithelial -mesenchymal interactions for branching morphogenesis. ACTA2 compartments dermal papilla stem cells for feather cycling. Novel usage of scale keratins strengthens feather sheath with SOX14 as the epigenetic regulator. We found this primary feather transition largely conserved in chicken (precocious) and zebra finch (altricial) and discussed the possibility that this evolutionary adaptation process started in feathered dinosaurs.


Introduction
Evolutionary innovations have enabled birds to occupy different ecological niches.In particular, feathers show a very high degree of diversity, providing an excellent model for studying how animals adapt to different environments 1,2,3,4 .Different developmental stages of a bird exhibit different types of feathers.
Most hatchling plumages are either naked or have natal downs, which confer heat conservation for the hatchlings.When juvenile birds are ready to leave the nest, most radially symmetric downs are replaced by bilaterally symmetric juvenile feathers in the same follicle to form basic plumage.After several rounds of molting, the adult feathers that can respond to environmental changes (hormones, seasons, etc.) are eventually formed (Fig. 1a) 5,6,7,8 .These remarkable morphological transitions are based on the stem cells and dermal niches from the same feather follicles, and can give plumages different morphologies and functions at different times of a bird's life for optimized function, sometimes giving birds totally different appearances; this process is known as "organ level metamorphosis" 9,10 .
The natal down to juvenile (rachidial) feather conversion is here referred to as the primary feather transition process, while the formations of diverse adult feathers after several rounds of molting is referred to as the secondary feather transition process.These two processes set up the foundation to generate regional and timing differences in feather cycling, and to achieve optimal functions required at the different stages of adult bird's life 8,10 .The topological changes of feather follicle architectures in natal and adult feathers have been well-characterized 11,12 .The roles of morphogens in forming rachidial and barb branch in adult feathers have also been elucidated 13,14,15,16,17,18,19 .The mechanism of secondary transition of feather follicles have been explored in several birds 20,21,22 .Yet, the evolution and underlying mechanism of the primary transition, which occurs only once in life in most birds, has rarely been studied.
Based on previous studies and our ndings, rachis, strong feather sheath, regeneration ability, and feather vane are major distinctions of juvenile feathers from natal downs (Fig. 1a) 5,11,12,23 .The rachis serves as the backbone of the feather, providing support and structural integrity 1,16,17,19,24 .Wnt signaling pathways are the major regulators of both hooklet and rachis development 16,17,19 .Feather sheath is a structure that maintains feather's cylindrical shape until it starts to disintegrate near the tip, allowing the mature part of the feather to unfurl 1 .In adult chickens, feathers need to undergo molting to maintain their normal functions and the molting process is known to be controlled by the dermal papilla and regulated by Wnt inhibitors 15,25,26 .Feather hooklets are found on the barbules of bird feathers, which interlock onto the proximal barbules of the immediately adjacent barb and transform the barb branches into a planar vane 19,27 .Here we will focus on the epigenetic controls of the primary feather transition.
To identify signaling pathways and transcription factor (TF) genes involved in the primary feather transition, we obtained embryonic and juvenile time series transcriptomes of chicken posterior dorsal skins and then used two approaches to analyze the data.First, we applied the time-ordered gene coexpression network (TO-GCN), a transcription factor prediction method, 28 to infer the TO-GCN of the TF genes expressed in the transcriptomes.Second, in birds, several β-keratin gene subfamilies have been classi ed by sequence similarities and tissue speci c expression, but little has been done on the regulation of their expression.Whether diversi ed β-keratin genes are involved in the primary feather transition and how they are regulated are interesting questions.Here, we annotated β-keratin genes in the newly released chicken genome (GRCg6a) to analyze the transcriptomes and applied ATAC-seq to elucidate the epigenetic regulation of the speci c β-keratin gene subfamilies during the primary feather transition 29,30 .We then functionally validated the predicted molecules in chicken ight feathers.
To test the molecular regulatory conservation in birds, we compared the primary feather transitions in chicken (white leghorn) and zebra nch (Taeniopygia guttata).Chicken, a precocial bird, belongs to Galliformes, which is a basal lineage of Neoaves.Zebra nch belongs to Passeriformes, which is the most derived Neoaves and the largest avian order; all species in this order are altricial birds 31,32,33 .We chose chicken and zebra nch as two model birds for our study because they represent two distant species in Neoaves, so that a trait conserved in both species is likely conserved in most or all Neoaves.

Results
The primary feather transition in chicken and zebra nch: from natal downs to juvenile feathers To study the morphological changes during the primary feather transition in the posterior dorsal skin, we compared the histology of chicken feather follicles between embryos (E14) and hatchlings (posthatch day 3 to day 7, i.e., D3 to D7), and then compared the structures with that of the posterior dorsal feather follicles in D7 zebra nches.In this region, juvenile feathers grew out from the skin around D6 in both species (Fig. 1a and Supplementary Figs. 1 and 2), suggesting that the time point of the transition initiation has been conserved in both species.
To understand the morphogenesis of feather follicles during the primary transition, both cross and longitudinal para n sections with hematoxylin and eosin (H&E) staining were applied to the whole skins or feather follicles of both species (Fig. 1b-m); muscle bundles were used as the anchor for the comparisons among different stages and species.In cross sections, we found that the structures of the feather germ were similar while the sizes increased slowly from E14 to D4, suggesting a slow growing phase (Fig. 1b-d, 1b'-d', and g).Interestingly, muscle development was relatively obvious in this period, and the shapes of muscle bundles between follicles changed from square in embryonic skins to diamond in juvenile skins (Fig. 1b'-d').The juvenile feather follicles started to enlarge quickly and formed rachis after D5, showing a fast growth phase (Fig. 1d'-e', 1k-l, and Supplementary Fig. 1).In D7, the juvenile feathers in chicken and zebra nch showed similar maturities: thickened feather sheathes, rachis and barb ridges were clearly visible (Fig. 1e'-f' and l-m).The clear bard ridges suggest the formation of organized barbs and barbules (Fig. 1g).
To con rm the above ndings, longitudinal sections were conducted in the central (spinal) feather follicles and the results supported the conclusion that the maturity of D7 feather follicles is similar between chicken and zebra nch: curved collars, enlarged pulps, and biconcave dermal papilla are clearly visible (Fig. 1h-m).The peripheral pulp regions show condensed cells, suggesting the formation of peripheral pulps (Fig. 1h).Based on these data and the literature 1,3,11,12,14,15,16,25 , we propose ve major evolutionary morphogenesis events during the primary feather transition in Neoaves: biconcave dermal papilla formation, peripheral pulp formation, rachis formation, vane formation, and feather sheath thickening.

Coordination of multiple signaling pathways during the primary feather transition
As we wanted to know the overall skin tissue changes during the primary feather transition, we sampled the whole posterior dorsal skins from day 3 to 7 hatchlings in which complex gene interactions from different cell populations and different timings are involved.We applied the time-ordered gene coexpression network (TO-GCN) method 28 , which was designed to decipher the molecular regulations from time (or developmental) course transcriptomes of complex tissues, to analyze the transcriptomes from embryonic and juvenile skins.We focused mainly on transcription factor (TF) genes, because they are the major players of gene regulation.
RNA-seq analysis of the transcriptomes from embryonic and juvenile samples was conducted to assess the quality of sequencing libraries and generate normalized read counts as the input for the TO-GCN analysis (Supplementary Fig. 3 and Table 1).A TO-GCN of 11 levels was constructed in both feather types based on the expression pro les of the TF genes.The lower (earlier) levels included the TF genes expressed at the earlier feather developmental stages, while the higher (later) levels included those expressed at the later developmental stages (Fig. 2a and b; Supplementary Table 2).The expression pro les of TF genes in each of the 11 levels were similar between the two feather types (as shown in the bar chart, Fig. 2a and b; Supplementary Table 2).
In early embryonic skin development (levels 1 to 5), Rho GTPases and cell cycle related pathways are enriched (Supplementary Table 3).Rho GTPases are best known for their roles in regulating cytoskeletal rearrangements, cell motility, cell polarity, axon guidance, vesicle tra cking and the cell cycle 34,35 , suggesting that Rho GTPase-mediated cell proliferation and polarization are the major event at this stage.Unexpectedly, polycomb repressive complex 2 (PRC2) related pathways were also enriched.PRC2 is a multiunit epigenetic protein complex that silences gene expression by catalyzing trimethylation of histone H3 at lysine 27 (H3K27me3).How methylation is involved in natal down growth regulation is an interesting question.In late embryonic skin samples (levels 6 to 11), corni ed envelope and keratin formation pathways were enriched (Supplementary Table 3).
In the juvenile TO-GCN, there are two distinct level groups each with a large gene member, suggesting two morphogenesis groups in post-hatched feathers.The rst group is at levels 3 (110 TF genes) and 4 (102 TF genes) in which the major enriched pathways are extracellular matrix (ECM) organization and metabolism (Fig. 2b and e, Supplementary Tables 3 and 4).The second group is at levels 9 and 10 in which 72 and 138 TF genes were assigned, respectively.Many molecular pathways were enriched in this period but WNT/β-catenin related pathways were dominant (Fig. 2b and e, Supplementary Tables 3 and   4).ECM is vital for determining and controlling the most fundamental behaviors and characteristics of cells such as proliferation, adhesion, migration, polarity, differentiation, and apoptosis 36,37 .The WNT/βcatenin related pathway is known to control rachis and hooklet formation 16,17,19 .Both of them are highly associated with the speci c structures of juvenile feathers and are therefore the focus in this study.
In addition to ECM and the WNT/β-catenin related pathways, many other known molecular mechanisms were also enriched.In both early embryonic and juvenile skin development (levels 1 to 5), muscle formation genes were enriched, consistent with the previous nding that muscle development is important for feather positioning (Supplementary Tables 3 and 4) 38 .In late juvenile skin development (levels 6 to 11), planar cell polarity (PCP) signaling pathway genes were enriched.Coupling of apicalbasal polarity and PCP was identi ed to interpret the WNT signaling gradient, which controls the bilateral symmetric feather formation (Supplementary Table 4) 18 .
Next, we wanted to know speci c TF genes that are involved in the primary feather transition.In this study, we used the following criterion to compare the embryonic and juvenile TO-GCNs: for a TF, if d (= embryonic level -juvenile level) ≥ 3, we say it is an embryo speci c TF; if d ≤-3, we say it is a juvenile speci c TF; if d ≤ |2|, we say it is a common TF (Fig. 2c).Among the 638 expressed TF genes, 443 showed d ≤ |2|, suggesting that most TF genes are commonly used in both feather types (Fig. 2c; Supplementary Table 2, 3, and 4).The embryonic or juvenile speci c TF genes at levels 8 to 10 in both feather types were selected as the embryonic or juvenile keratinization regulators, respectively.We matched these TF genes with those in two previous studies in which keratin regulators were well characterized (Supplementary Table 5) 39,40 .Interestingly, compared to our TF genes identi ed from embryonic TO-GCN, our TF genes identi ed from juvenile TO-GCN overlapped more with TF genes controlling scale or claw keratins (embryonic TFs: 0/9 = 0% vs juvenile TFs: 3/8 = 37.5%, Supplementary Table 5), suggesting that scale/claw keratins and their regulators may contribute more to juvenile feathers than to embryonic feathers.Wnt signaling is the major regulator of rachis formation and LEF1 is a key molecular hub converting radial downy to bilaterally symmetric juvenile feathers The molecular gradient in feather follicle from anterior to posterior end is crucial for regulating the angles of barb ridges for rachis formation 14 .During rachis formation, the anterior to posterior Wnt3a gradient is known to convert radial to bilateral feather symmetry via convergence of barb ridges toward the rachis region 16 .Moreover, multiple Wnt genes showed gradient expressions during ight feather regeneration 26 , suggesting that the Wnt-based regulation could be redundant or regional speci city.
In the RNA-seq analysis, the expression levels of Wnt2b, Wnt5a, Wnt5b, Wnt7a, Wnt9a, Wnt9b, and Wnt16 were increased with the juvenile dorsal feather development, corresponding to the time point of the rachis formation (Fig. 1b'-e' and Fig. 3a).Wnt3a showed constant expression and might not be a rachis regulator of the juvenile dorsal feathers (Fig. 3a).In our TO-GCN analysis, the lymphoid enhancer binding factor 1 (LEF1) was identi ed as a level 10 TF in juvenile feather development and is known to be a key mediator of Wnt signaling in diverse biological processes (Fig. 3a and Supplementary Table 2).The section in situ hybridization (SISH) of LEF1 and CTNNB1 in D7 chicken posterior dorsal skins revealed that LEF1 was expressed in an anterior to posterior gradient in the feather epidermis but CTNNB1 was expressed evenly (Fig. 3b), suggesting that LEF1 is the key factor responding to Wnt signaling.In zebra nch, although LEF1 was also expressed in feather epidermis, the feather pigments prevented us from visualizing the gradient (Fig. 3b, lower panel).Interestingly, both CTNNB1 and LEF1 showed higher expression in the chicken central feather follicles than in the anking follicles, suggesting that the embryonic feather lateral propagation factor β-catenin could also be the activator for the secondary feather lateral propagation and LEF1 is also involved (Fig. 3b, upper panel).
To validate the effect of LEF1 on rachis formation, a dominant negative form of LEF1 (dnLEF1) was cloned into RCAS virus and injected into cavities of the plucked ight feathers.Previous studies showed that in ight feathers, the misexpression of Wnt inhibitor DKK1 could slightly disturb the rachis formation, and the overexpression of Wnt3a could disrupt the rachis formation and also cause abnormal barbs 16 .Here, all the ight feathers injected with RCAS-dnLEF1 showed defects in rachis (Supplementary Fig. 4a).Half of the feathers lost part of the rachis while the others lost the entire rachis without in uencing the surrounding barbs (Fig. 3c, Supplementary Fig. 4b), suggesting that the function of LEF1 is speci c in feather follicles.SHH is a key morphogen for barb formation via their expression in bard ridges and absence in rachis forming regions 14,16,41 .In rachidial feather, barb ridges insert into the rachidial ridge with the helical insertion angle (Fig. 3d, indicated by θ).However, in both natal down and RCAS-dnLEF1 misexpressed ight feather, all the barb ridges were formed in parallel (θ = 0, Fig. 3d), suggesting that barb ridges were not able to insert into rachidial ridge.
ECM re-organization generates peripheral pulp for feather branching morphogenesis ECM organization was speci cally enriched in three levels of the juvenile feather TO-GCN (L2 to L4, Fig. 2e).Components of ECM link together to form a structurally stable composite, contributing to the mechanical properties of tissues.ECM is also a reservoir of growth factors and bioactive molecules 36, 37 .In chicken embryonic feathers, ECM regulates mesenchymal mechanics which can spontaneously break skin symmetry 42 .In chicken adult feathers, ECM reorganization enables peripheral pulp formation 25 .Therefore, we asked what is the function of ECM in primary feather transition and which molecules control and maintain the ECM reorganization process?
In the juvenile feather TO-GCN, levels 2 to 4 basically correspond to the D3 to D5 stages in feather development (Fig. 2b).A novel morphogenesis in this stage is the generation of peripheral pulp, which has about ve layers of mesenchymal cells closely attached to the feather lament and basement membrane.The peripheral pulp expands the epithelial-mesenchymal interactive interface for barb patterning 25,43 .The gain of ECM mediated pulp differentiation is therefore our hypothesis for the primary feather transition.Tenascin C (TNC) is frequently used as ECM and differentiated pulp markers 11, 12, 15, 25, 26 .The immunochemistry (IHC) signals of TNC reveal that the peripheral pulp was gradually differentiated from embryonic pulp along with the growth of embryonic to juvenile feathers in both chicken and zebra nch (Fig. 4a'-e' and 4a"-e").We further picked up two ECM reorganization related transcription factors, TWIST2 and ZEB2, from levels 2 to 4 of TO-GCN to examine their expressions.The SISH of TWIST2 showed the initial expression in the collar of embryonic feathers, the induction in whole pulp in early juvenile feathers, and the restricted expression in apical dermal papilla and peripheral pulp in late juvenile feathers in both chicken and zebra nch (Fig. a1-e1), suggesting its role in peripheral pulp formation.Interestingly, the IHC of ZEB2 showed an almost opposite pattern in that the signals were enriched in the central pup in both species (Fig. a2-e2), suggesting that it may maintain the mature tissues.
ACTA2 shapes adult dermal papilla to compartment dermal papilla stem cells for cyclic renewal ACTA2 (encoding α-SMA) is a major feather dermal papilla marker and our TO-GCN analysis showed that ACTA2 co-expressed with levels 2 to 4 TF genes (Supplementary Table 6) 15,25 .We therefore paid additional attention to the role of α-SMA in primary feather transition.Dermal papilla of downy feather is long and slender, while dermal papilla in the juvenile and adult feathers are biconvex-shaped (Fig. 4 left side schemes) 11,12,15,23,25 .We found that, in the IHC on skin longitudinal sections, the expression of α-SMA in dermal papilla increased with the juvenile follicle development and eventually lled up the biconvex-shaped dermal papilla in D7 juvenile feathers in both chicken and zebra nch (Fig. 4a-e).Most interesting, the study of feather cycling showed that in resting phase, dermal papilla stem cells are located in the apical part of the bi-concaved dermal papilla 21 .In the growth phase, the activation of apical part dermal papilla generates pulp progenitors, giving rise to both central pulp for nutrition purpose and peripheral pulp for continuous epithelial -mesenchymal interactions required for feather branching morphogenesis.
We then tested the importance of ACTA2 in dermal papilla formation.Since the protein sequence of ACTA2 is identical between human and chicken, we cloned a published human dominant negative ACTA2 form into RCAS virus (RCAS-dnACTA2) and injected it into the cavities of the plucked ight feather in chicken 44 .The knock-down virus caused such severe effects that most (7/8) of the injected follicles could only regenerate tiny or no visible feathers two weeks after injection (Supplementary Fig. 4d).All the eight injected follicles could not be regenerated after the second plucking.Histologically, the dermal papilla and part of the collar structures were also disrupted (Fig. 4f and f').To further understand the molecular changes in those abnormal follicles, immunostaining of several known ECM factors as well as ACTA2 were applied to feather sections (Fig. 4g-j') 25 .Dermal papilla markers α-SMA and vimentin (VIM) were enriched in whole dermal papilla in the control feather follicles but attenuated and dispersed in the virus injected follicles (Fig. 4g and g', i and i').Unexpectedly, TNC and Neural Cell Adhesion Molecule (NCAM), which were found to be located in papilla ectoderm in the control feather follicles, were attenuated and dispersed in the virus injected follicles (Fig. 4h and h', j and j'), suggesting the simultaneous disruption of dermal papilla and pulp structures.It could be that the dermal papilla of juvenile feathers contains myo broblast cells and knock-down of α-SMA may then affect the integrity of these cells.These results suggest that α-SMA not only is a structural protein but also establishes an environment for the primary feather transition in Neoaves.
Many scale keratins are speci cally upregulated in juvenile feather sheath by SOX14 β-keratin genes are mainly classi ed into feather, scale, claw keratins, and keratinocytes 45,46,47 .A previous study revealed that chicken dorsal natal downs mainly express feather-β-keratin genes on Chr1, Chr10, Chr25 and also some members on Chr27 during keratinization, whereas feather-β-keratin genes on Chr2 and Chr6 are exclusively enriched in adult wing feathers 45 .However, subsequent studies revealed that β-keratin genes are differentially regulated in different skin regions 39,40,47 .Since our TO-GCN analysis along with the previous studies suggested the importance of scale keratins in the juvenile feather formation 39,40 , we employed several modi cations in the transcriptomic analysis to decipher the β-keratin gene regulation: (1) Our embryonic and juvenile feather tissues were only from chicken posterior dorsal skins.(2) In addition to β-keratin genes, we analyzed the whole epidermal differentiation complex (EDC), which was known to participate in feather functional evolution 48, 49, 50, 51 .(3) We used the newly published chicken genome (GRCg6a), which has a much-improved assembly in micro-chromosomes, where EDC gene clusters are located.(4) We manually annotated the EDC genes, especially the β-keratin genes, which are still poorly annotated in GRCg6a.(5) We conducted an ATAC-seq analysis of embryonic and adult feather tissues to look for epigenetic regulators during the primary feather transition.
We analyzed the co-expression pro le of EDC genes in both feather types during their development.Our data were partially inconsistent with the previous nding, showing that, regardless of chromosomal locations, most of the β-keratin genes were expressed during both embryonic and juvenile feather keratinization (Clusters 7 and 11, Fig. 5a).Genes in Cluster 11were highly expressed at the keratinizing stage in both feather types, which include most feather β-keratin genes located on Chr1, 2, 7, 10, 25, and 27 (Fig. 5a; Supplementary Table 7).Genes in Cluster 7 are highly expressed at both the keratinizing stage and a stage ahead of it in both feather types, which include other β-keratin genes, such as scale, claw, feather-like β-keratin, and feather-β-keratin genes on Chr10 (Fig. 5a; Supplementary Table 7).In addition to β-keratin genes, most other EDC genes also were in cluster 7 or 11 (Fig. 5a; Supplementary Table 7).
Only small co-expression clusters and individual genes showed distinct pro les between the two feather types.Two β-keratin genes in Clusters 8 and 9 (GG6AChr25Ktn11 and GG6AChr25FK5) were highly expressed in keratinizing juvenile feathers but not in natal downs (Fig. 5a; Supplementary Table 7).In Cluster 7, CRNN was highly expressed in embryonic but not in juvenile feathers (Fig. 5a; Supplementary Table 7).Interestingly, when we analyzed the β-keratin genes based on their subfamilies (β-keratin related proteins, feather keratins, feather keratin like proteins, scale keratins, claw keratins, and keratinocytes), βkeratin gene expression differences between the two feather types could be found.Chr25 feather and feather-like keratin genes expressed higher in natal downs than in juvenile feathers, while Chr25 scale keratin genes expressed higher in juvenile feathers (Fig. 5a and b, Supplementary Fig. 6, Supplementary Table 7), suggesting that: (1) the ratio differences of Chr25 β-keratin gene expression in different chromosomes may contribute to feather type differences, and (2) some scale keratin genes could be important speci cally for juvenile feather formation.To understand the expression differences of the scale keratin genes in the two feather types, we conducted SISH with GG6AChr25Scale2 and GG6AChr25Scale10 in cross sections of embryonic and juvenile feather follicles (Fig. 5c-f and Supplementary Table 7).Although the two keratin genes were faintly expressed in the barbs of both feather types, they were highly enriched in the feather sheath in juvenile feathers but not in natal downs.
What TFs control the scale keratins and whether this regulation is conserved in Neoaves were the next questions.All the scale keratins are on chromosome 25 for both chicken and zebra nch.We rst compared the synteny of chromosome 25 between the two species and found that they are basically conserved along the whole chromosome.Next, we conducted ATAC-seq to compare the differential genomic accessible regions between embryonic feather laments and the epidermis of adult ight feathers.Figure 5g shows the ATAC signals of the two tissues surrounding scale keratin gene cluster.The differential accessible regions overlapped with the conserved chromosome 25 sequences were extracted for footprint analysis to detect the enriched binding motifs and TFs.The results revealed that TFs SOX14, ESRRB, ESRRG, PRDM4, SREBF2, and SMAD5 were enriched in juvenile feathers; the corresponding binding motifs are shown in Fig. 5g'.Related studies are limited but SOX14 is essential for the initiation of neuronal differentiation in the chick spinal cord 52 .
Since SOX14 is the most signi cant TF and more highly expressed in juvenile feather than in embryonic feather transcriptomes (Fig. 5g and Supplementary Table 1), we used SISH to detect its expression in the natal and juvenile feather follicles of chickens and zebra nches.The results showed that, the same as the scale keratins, SOX14 was detected in barbs of both feather types in both species, but the expression was only enriched in the feather sheath of juvenile follicles in both species (Fig. 5h-k).The overlapped expression of scale keratins and SOX14 suggests that SOX14 is an upstream transcription factor of scale keratins.To functionally validate the prediction, we overexpressed SOX14 in embryonic feather using the RCAS-virus system and 1/3 of the embryos (N = 12) generated abnormal feather follicles with either round or thickened and shortened phenotypes (Fig. 5l).Since the skin regions included both abnormal and normal feather follicles, SISH was applied to detect the expression of SOX14 and scale keratins in the same skin sections.Follicles with higher SOX14 expressions showed abnormal phenotypes and enhanced expressions of GG6AChr25Scale2 and GG6AChr25Scale10 (Fig. 5l), suggesting that SOX14 is the regulator of the scale keratins in feather sheath.

Discussion
Feather diversity appears in different body regions and at different developmental stages of a bird 5,23 .Previous studies revealed regulatory differences among different body feathers or among different parts of a chicken feather 3,47,53,54,55 .Feather transition represents a novel type of timing control evolved for better adaptation, depending on different needs at different life times, but the underlying molecular controls have not been well characterized.Feather transition is at least achieved by ve tissue reorganizations and our study revealed four of them: 1. Wnt is the common signaling for rachis formation and LEF1 is the downstream key hub for diverse Wnt proteins.2. ECM reorganization is essential for peripheral pulp formation for further branching morphogenesis.3. ACTA2 is a key factor to compartment dermal papilla formation for feather cycling.4. Scale keratin is recruited from scale to strengthen the sheath of juvenile feathers and SOX14 appears to be a major activator for this co-option process.Feather vane formation achieved by hooklets is the primary feather transition property that was not included in this study.Although our transcriptomes included the initial time points for hooklet development (D7), we did not identify regulators of the known morphogen WNT2B 19 .Since hooklets only contribute to a feather fraction, the ner analysis may be worth pursuing in the future.The evolutionary morphogenesis events and the underlying molecular regulators during the primary feather transition are summarized in Fig. 6.All the chicken primary feather transition factors showed the same expression patterns in zebra nch feather follicles, suggesting that the primary feather transition has been conserved in Neoaves.
Multiple Wnt genes seem to participate in the primary feather transition.Whether it is simply functional redundancy or regional speci city is an interesting question.Hox genes represent the best example that elicit distinct developmental programs along the head-to-tail axis of animals and are the upstream signals of many Wnt proteins 56,57,58,59,60,61 .Whether distinct HOXs activate different WNTs and generate regional speci c feather rachis is another interesting question.In feather regeneration, our previous study has uncovered a cell movement from apical dermal papilla to peripheral pulp along the collar bulge 25 .This remarkable feature allows prolonged interactions between the dermal niche and epidermal stem cells, providing a tunable interface that is essential for barb patterning 25 .The path of the cell movement corresponds to the position of ECM molecule expression and so ECM reorganization could be a key feather regeneration factor.ACTA2 is a dermal papilla marker but is not known to regulate ECM molecules like TNC and NCAM.While the knock-down of ACTA2 disrupted their expressions, we proposed that ACTA2 could be used to establish a competent environment for molecular interaction during peripheral pulp and apical dermal papilla formation.Interestingly, we identi ed TWIST2, whose expression during feather transition corresponds to the position of apical dermal papilla and peripheral pulp, suggesting its possible role in feather regeneration.TWIST2 is known to regulate ECM in some cases 62,63,64 .
Previous studies have indicated that SOX genes are important for scale keratin regulation and overexpression of SOX18 enhanced the expression of a scale keratin in embryonic feather 39,40,65 .However, different SOX genes were identi ed as β-keratin regulators in different skin appendages (SOX10 in scale, SOX18 in embryonic feather sheath, and SOX14 in juvenile feather sheath in this study) 39,40,65 .It implies that SOX genes perform regional or temporal speci c regulations.Moreover, based on the differential SOX14 regulation and ATAC-seq analysis, we hypothesize that the embryonic feather sheath can be suppressed epigenetically (Fig. 5).The primary feather transition activates the juvenile feather sheath cells to recruit SOX14 to trigger the expression of β-keratins to enhance keratinization (Fig. 6).
Notably, a juvenile bird shows various feather types 23 .Although we focus on the major one, i.e., the contour feathers with vane structure, whether other feather types, like bristles and loplumes, are achieved by the same molecular mechanisms is another interesting question.
Precocial chickens are covered with natal downs, while altricial zebra nches have only limited natal downs in hatchlings 31,32,33 (Supplementary Figs. 1 and 2).Despite these differences between chicken and zebra nch hatchlings, their juvenile feather developmental pro les are quite similar 31 (Fig. 1).After several rounds of molting, diverse feathers are derived from the juvenile feather follicles in the secondary transition process, as seen in the sail-feathers of the mandarin duck (Aix galericulata) or colorful contour feather of the pheasants 8, 20, 66, 67 .The conserved primary feather transition suggests its functional importance for survival while the following diverse feather transitions are mainly for mating choices.Interestingly, although egg incubation in chickens is around 21 days while that in zebra nches is around only 14 days, their juvenile feathers are visible at similar post-hatch stages (Fig. 1e' and f'), suggesting that the primary feather transition is sensitive to the post hatch stimulus, but are regardless of nesting habitats and parental cares.
The conserved primary feather transition among Neoaves implies conserved regulatory sequences in their genomes.How to identify the regulatory regions and how the regulations are controlled are important questions.Although many chicken mutants have been created, the phenotype without the primary feather transition has never been found 68 , suggesting that the loss of feather transition is highly detrimental.The avian conserved non-exonic elements (CNEEs) database could be used to screen the regulatory sequences and the combination of epigenetic marks can be used to understand the epigenetic changes of the regulatory sequences 69,70 .On the other hand, when the mechanism evolved is a fascinating question.The evolution of feather structures in Dinosauria was proposed in the previous studies in which the basal family Compsognathidae (such as Sinosaurpteryx) was covered with unbranched mono laments (down feather without barbules) while the derived family Caudipteridae (such as Caudipteryx) showed feather types similar to current birds 71,72,73,74 .Whether Dinosauria encompassed feather transitions during their development is not clear, but both the natal and juvenile feather forms had been discovered, suggesting that the primary feather transition in current birds could either be inherited or modi ed from their dinosaur ancestry.
The success of the Aves to venture into new eco-spaces relied much on the exibility to convert their feather phenotypes during the lifetime of an individual, which is based on the module-based complexity formation 9 .These modules are feather follicles consisting of stem cells and their niches.In this study, we tried to decipher the strategies used to evolve these regulatory switches in these modules.We found no new molecules needed to be created.Instead, the transcription regulation redeployed existing molecular module and cellular components, which include the topological positioning of Wnt signaling center, the tissue remodeling to extend the interface for epidermal-dermal interactions, the re-shaping dermal papilla to create dermal papilla stem cell compartment and the co-optional use of scale keratins in feather sheath.Conversions of integument phenotypes from the young to the adult are also seen in some mammals and other animals.These transitions may have been driven by environmental adaptation 75 .Some strategies we learn here may help understand this fundamental mechanism in Evo-Devo.

Ethics statement
All the animals used in this study were processed following the approved protocols of the Institutional Animal Care and Use Committees of the University of Southern California (USC; Los Angeles, CA, USA) and Institutional Animal Care and Use Committees of National Chung Hsing University (NCHU, Taichung, Taiwan).

Sample collection
Fertilized pathogen-free White Leghorn chicken eggs were purchased from Charles River Laboratories.
Eggs were incubated at100℉ and 65% humidity until embryos reached the desired developmental stages.Chicken embryos at14 days (E14), posthatched 3 days (D3), posthatched 4 days (D4), posthatched 5 days (D5), posthatched 6 days (D6), and posthatched 7 days (D7) were used.Chicken hatchlings are mostly covered by natal downs.To visualize the juvenile feather growth, we plucked posterior dorsal natal downs of the hatched chicken and observed the growth patterns of juvenile feathers (Fig. 1).Although a previous mammalian study found that hair plucking can stimulate early regeneration of underneath hairs 76 , this phenomenon is not observed in the primary feather transition because the plucked and un-plucked region showed similar juvenile feather growth in both species (Fig. 1 and supplementary Figs. 1 and 2).The zebra nch eggs were hatched and the hatchlings were raised in the bird breeding room at NCHU.E12 and D7 Zebra nches were used.

Quantitative PCR
To quantify the candidate gene expressions, the cDNAs were synthesized from the total RNA by QuaniTect Reverse Transcription kit (Qiagen).Each cDNA sample containing SYBR green (KAPA SYBR FAST qPCR kit) was run on LightCycler 480 (Roche) under the appropriate conditions.Quanti cation of the TATA box binding protein (TBP) RNA was used to normalize target gene expression levels.SHH forward primer: CTGGTGAAGGACCTGAGCCCT; reverse primer: GCCCAACTGTGCTCCTCGAT, TBP forward primer: CACAGCAAGCGACACAGGGA; reverse primer: AGGTGTGGTTCCCGGCAAAG.

Collection and construction of the transcriptomic libraries
To choose the appropriate time points representing juvenile feather development, we screened several feather morphogens and found that SHH is a good marker gene because its expression corresponded to the phenotypic changes of the juvenile feather follicles (Fig. 1 and Supplementary Fig. 3a).Based on the expression pro le of SHH, we selected posthatch day 3 to 7 posterior dorsal skins of the chickens to represent the primary feather transition.All the chicken transcriptomic libraries from the posterior dorsal area used in this study and their sources are listed below: For posthatch samples we collected the whole skin because the feather follicles cannot be quickly dissected to ensure the RNA quality at these stages.The dissected tissues were preserved in RNALater solution (Invitrogen), incubated at 4°C overnight, and then transferred to -20°C until processing for isolation of total RNA.Total RNAs were isolated using the RNeasy Fibrous Tissue Mini Kit (Qiagen).The 30 min DNaseI treatment was carried out at room temperature.The RNA quantities and qualities of each individual were analyzed by NanoDrop 1000 (Thermo Scienti c) and BioAnalyzer II (Agilent Technologies).The RNA samples from the same litter that passed the quality control (RNA integrity number [RIN] > 8.0, A260/280 and A260/230 > 1.9) were used for sequencing library constructions.Pairedend 2 × 101 nt sequencing for juvenile samples and paired-end 2 × 150 nt sequencing for embryonic samples were conducted by the High Throughput Genomics Core Facility, Biodiversity Research Center, Academia Sinica, Taiwan, and Novogen, United States, respectively.

Manual annotation of keratin genes
The annotations of the latest chicken genome assembly (GG6a) were downloaded from Ensembl 77 and RefSeq 78 .The Hidden Markov Models (HMMs) that represent alpha-keratins and beta-keratins were downloaded from Pfam 33.1 79 .The putative alpha-keratins and beta-keratins in Ensembl and RefSeq annotations was predicted using the hmmsearch function in HMMER 3.1 (https://hmmer.org) 80.For alpha-keratins, we further eliminated the genes that were not located in the typical type I and type II alphakeratin gene clusters.Finally, the predicted alpha-keratin genes and beta-keratin genes in Ensembl and RefSeq annotations were further compared to the annotation in Ng, et al. ( 2014) 46 and the nal annotation of alpha-keratins and beta-keratins was obtained by manual judgement for the most appropriate gene models.

RNA-seq analysis
For each transcriptome, low-quality reads were trimmed and adapters were removed using Trimmomatic 81 .The processed paired-end reads were mapped to the reference genomes of chicken (GRCg7w for TO-GCN and manually-annotated GRCg6a for keratin analysis) using Hisat2 with default settings 82 .Two genomes were used because the manually annotated GRCg6a was made before GRCg7w was released.
From the aligned reads, the transcripts were assembled using StringTie with annotation-based settings 83 .The gene expression level was calculated in terms of Transcripts Per Kilobase Million (TPM).For TPM calculations, uniquely mapped reads were rst assigned to the genes.Multiple-hit reads were then redistributed to genes based on their relative abundances of uniquely mapped reads.A gene is considered expressed if its TPM is ≥ 1 in at least one of the transcriptomes.To compare expression pro les of the selected genes across the conditions, we applied the upper quartile normalization procedure 84 .Differential expression analysis between different samples was calculated using the R package DESeq2 85 .
Construction of TO-GCNs.
Our comparative transcriptomics method 28 was designed to analyze time-course or developmental-stage transcriptomes that have ve or more different conditions between the two feather types: E7 epidermis, E7 dermis, E9 epidermis, E9 dermis, and E14 feather laments in posterior dorsal regions for embryonic samples; D3, D4, D5, D6, and D7 posterior dorsal skins for juvenile samples.The method consists of three steps: determining co-expression cutoffs, constructing GCNs, and determining the time order of TF gene expression to transform a GCN into a TO-GCN.First, the Pearson correlation coe cient (PCC) values of all TF-TF gene pairs were calculated and used to determine the cutoffs of co-expression (for p < 0.05, PCC > 0.92 and 0.9 for natal down and juvenile feather, respectively).Second, using the co-expression relationships between TF genes, we determined the GCN for connected TF genes.Third, the time order of TF genes in the GCN was assigned by the breadth-rst search (BFS) algorithm 86 initiated from the selected node that should be the rst upregulated TF in the GCN.BFS is an algorithm for searching a network graph.It starts with an initial seed and searches all its neighbors (nodes with connecting edges) to form a set of nodes (level 1).Then, the process proceeds from all nodes in level 1 and searches their neighbors (excluding level 1 nodes) to form the second set of nodes (level 2) and so on, until all nodes in the network are assigned.The computer programs for the method are available at https://github.com/petitmingchang/TO-GCN(31).SP5 and MYOD1 TF genes were selected as the initial nodes in E7 embryonic samples and D3 juvenile samples, respectively.
Co-expressed gene sets and overrepresented functions at each TO-GCN level For the TF genes at each level of a TO-GCN, a corresponding set of co-expressed genes (usually non-TF genes) can be identi ed with the same co-expression relationship for adding the genes to the TO-GCN.
Since a gene may be co-expressed with TF genes in multiple levels, two neighboring gene sets may have some overlapping genes.For each set of genes corresponding to a level in a TO-GCN, a functional enrichment analysis was conducted with the background set of all expressed genes in this study.Fisher's exact test with the false discovery rate (FDR) < 0.05 was applied with functional annotations from Reactome (https://reactome.org/).

ATAC-seq libraries construction and analysis
The epidermis tissues of adult ight feather follicles were dissected using 2X CMF solution, following a previous method 40 .The dissected tissues were dissociated and the ATAC-seq libraries were constructed following the same study.Cutadapt was used to remove adapter sequences with the following parameters: -q 20 -cut 1 -length 75 -minimum-length 36 89 .Hisat2 with the following parameters was applied to reads mapping: --no-temp-splicesite -no-spliced-alignment -X 2000 82 .Samtools was used to remove low quality and non-unique hits with the parameters: view -b -f 3 -F 4 -F 8 -F 256 -F 2048 -q 30 90 .

Construction and misexpression of dominant negative form of the candidate genes
To perform tissue speci c gene knock down, the misexpression of sequence modi ed genes were used to achieve the dominant negative (dn) effect.The RCAS-dnLEF1 has been applied in chicken studies and was obtained from Addgene (Plasmid #14022) 93 .The dnACTA2 function was published previously 44 , and the sequence was cloned into RCAS-cherry plasmid in this study.Virus was made according to Jiang 87 and concentrated by ultracentrifugation.The ight feathers of SPF chickens were plucked and around 50ul of concentrated RCAS virus was injected into each follicle cavity.All the experiments were conducted in both sexes of chickens to eliminate the probable sexual speci c effects.The ight feathers that were plucked at the same time without virus injection were used as the control.Feather follicles including dermal papilla were dissected after two weeks of the virus injection and xed immediately in 4% PFA. The Construction of TO-GCNs and the enriched pathway in the primary feather transition.a and b The TO-GCNs of constructed from normalized gene expression levels during embryonic (a) and juvenile feather (b) developments.In a TO-GCN, each blue node represents a TF gene and each grey line represents a co-expression relationship between two TF genes.The number in each level represents the number of TF genes at that level.The mean z-score of TF gene expression at each level was used to generate the bar charts and heatmap.c The proportion of TF genes in TO-GCN level differences.Most TF genes show less than 3 level differences between embryonic and juvinile feathers (blue bars).dGene coexpression analysis and the heatmap.All the expressed protein coding genes are used for co-expression and cluster analyses after gene expression level normalization, and only the top expressed 6,000 genes were used to constructed the heatmap.e The pathways enriched in more than one level were selected from Supplementary Table 4. Epi: epidermis; Der: dermis; FF: feather lament; Skin: whole skin; L: level of TO-GCN; GO: gene ontology; FDR: false discovery rate.
The immunochemistry (IHC) and section in situ hybridization (SISH) of embryonic and juvenile feather in chickens and zebra nches, and the functional peturbation of ACTA2 in adult chicken ight feather follicles.The scheme of embryonic (upper) and juvenile (lower) feather follicles are shown in the left side.a-e" the immunostaing of α-SMA (ACTA2), TNC, and DAPI in embryonic chicken (a-a"), juvenile chicken (b-d"), and juvenile zebra nch (e-e") feather follicles.a1-e1 the SISH with TWIST2in the same tissues with a-e".a2-e2 the IHC with ZEB2 in the same tissues with a-e".f-j H&E staining (f) and immunostaing of α-SMA (g), TNC (h), VIM (i), and NCAM (j) in control follicles.f'-j' H&E (f') staining and immunostaing of α-SMA (g'), TNC (h'), VIM (i'), and NCAM (j') in RCAS-dnACTA2 injected follicles.CK: chicken; ZF: zebra nch; PP: pulp; pPP: peripheral pulp; cPP: centrral pulp; DP: dermal papilla; aDP: apical dermal papilla; bDP: basal dermal papilla; CL: collar epidermis.Scale bar: 100 um.Summary of the phenotypic and molecular changes during the primary feather transition in Neoaves.Five evolutionary morphogenesis events are shown their indicated by TO-GCN The identi ed genes for the primary feather transition are listed after the events and shown in either black font color, which are identi ed in this study, or gray color, which were identi ed in previous studies 15,16,17,19,25 .Peripheral pulp and apical dermal papilla formations for feather regeneration are the early major morphology and histology of posterior dorsal feathers at different developmental of chicken and zebra nch. a Scheme of the primary feather transitions from natal down (only one barb is shown with barbules) to adult feather in chickens and zebra nches was modi ed from previous studies 3, 5, 15, 19, 23, 25 .b-f H&E stainings of cross sections of the entire posterior dorsal skins.Dorsal views of developing chicken and zebra nch were embeded.A chestnut-anked white zebra nch mutant was used to visualize the juvenile feathers (Supplementary Figure 1).Yellow scale bar: 1 cm.b'-f' a close view of an individual follicle from b-f (indicated by red arrows), respectively.g A cross-section view of feather maturation.h, A longitudinal section view of the feather maturation.i-k H&E stainings of longitudinal sections of feather follicles on central posterior dorsal skins.FS: feather sheath; JF: juvenile follicle; PP: pulp; cPP: central pulp; pPP: peripheral pulp; FF: feather lament; DP: dermal papilla; BR: barb ridge; CL: collar; E: embryonic incubation days.D: posthatch days.A: anterior (head); P: posterior (tail).Black scale bar: 100 um.

Figure 3 The
Figure 3

Figure 5 Expression
Figure 5