Exploration of key regulators driving primary feather follicle induction in goose skin

Journal Pre-proofs Research paper Exploration of Key Regulators Driving Primary Feather Follicle Induction in Goose Skin Xuewen Hu, Xiaokang Zhang, Zhiwei Liu, Shaomei Li, Xinting Zheng, Yangfan Nie, Yingfeng Tao, Xiaoliu Zhou, Wenqing Wu, Ge Yang, Qianqian Zhao, Yang Zhang, Qi Xu, Chunyan Mou PII: DOI: Reference:

S0378-1119(20)30007-X https://doi.org/10.1016/j.gene.2020.144338 GENE 144338

To appear in:

Gene Gene

Received Date: Revised Date: Accepted Date:

16 November 2019 13 December 2019 6 January 2020

Please cite this article as: X. Hu, X. Zhang, Z. Liu, S. Li, X. Zheng, Y. Nie, Y. Tao, X. Zhou, W. Wu, G. Yang, Q. Zhao, Y. Zhang, Q. Xu, C. Mou, Exploration of Key Regulators Driving Primary Feather Follicle Induction in Goose Skin, Gene Gene (2020), doi: https://doi.org/10.1016/j.gene.2020.144338

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier B.V.

Research paper

Exploration of Key Regulators Driving Primary Feather Follicle Induction in Goose Skin Xuewen Hu a, y, Xiaokang Zhang a, y, Zhiwei Liu a, Shaomei Li a, Xinting Zheng a, Yangfan Nie a, Yingfeng Tao a, Xiaoliu Zhou a, Wenqing Wu a, Ge Yang a, Qianqian Zhao a, Yang Zhang b, Qi Xu b and Chunyan Moua, * a

Laboratory of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, College of Animal Science and Technology, Huazhong Agricultural University, Wuhan, 430000, China b College of Animal Science and Technology, Yangzhou University, Yangzhou, 225000, China * Correspondence: [email protected] y These authors equally contributed to this work.

Highlights .Transcriptome profiling reveals Wnt candidates in primary feather follicle induction. .Six Wnt pathway genes exhibit upregulation in primary feather follicle induction. .LEF1 and DKK1 are potential regulators to specify presumptive feather tracts. .WNT16 and FZD1 are potential regulators to induce primary feather follicle primordia. .SFRP1 and FRZB are potential regulators to confine primary feather follicle primordia. Abstract The primary feather follicles are universal skin appendages widely distributed in the skin of feathered birds. The morphogenesis and development of the primary feather follicles in goose skin remain largely unknown. Here, the induction of primary feather follicles in goose embryonic skin (pre-induction vs induction) was investigated by de novo transcriptome analyses to reveal 409 differentially expressed genes (DEGs). The DEGs were characterized to potentially regulate the de novo formation of feather follicle primordia consisting of placode (4 genes) and dermal condensate (12 genes), and the thickening of epidermis (5 genes) and dermal fibroblasts (17 genes), respectively. Further analyses enriched DEGs into GO terms represented as cell adhesion and KEGG pathways including Wnt and Hedgehog signaling pathways that are highly correlated with cell communication and molecular regulation. Six selected Wnt pathway genes were detected by qPCR with up-regulation in goose skin during the induction of primary feather follicles. The localization of WNT16, SFRP1 and FRZB by in situ hybridization showed weak expression in the primary feather primordia, whereas FZD1, LEF1 and DKK1 were expressed initially in the inter-follicular skin and feather follicle primordia, then mainly restricted in the feather primordia. The spatial-temporal expression patterns indicate that Wnt pathway genes DKK1, FZD1 and LEF1 are the important regulators functioned in the induction of primary feather follicle in goose skin. The dynamic molecular changes and specific gene expression patterns revealed in this report provide the general knowledge of primary feather follicle and skin development in waterfowl, and contribute to further understand the diversity of hair and feather development beyond the mouse and chicken models. Abbreviations BMP, bone morphogenetic protein; cDNA, DNA complementary to RNA; DEGs, differentially expressed genes; E, developmental stage of the embryo; ECM, extracellular matrix; Epi, epidermis; FGF, fibroblast growth factors; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; KOG, euKaryotic Ortholog Groups; MAPK, mitogen-activated protein kinase; NCBI, National Center for Biotechnology Information; Nr, NCBI non-redundant protein sequences; Nt, NCBI non-redundant nucleotide sequences; Pfam, Protein family; qRT-PCR, quantitative reverse transcriptase polymerase chain reaction; RNA-seq, RNA sequencing; Shh, sonic hedgehog; Swiss-prot, a manually annotated and reviewed protein sequence database; Wnt, wingless-type MMTV integration site family members

Keywords skin appendage; morphogenesis; primordia; transcriptome; Wnt

1. Introduction Birds have diverse feathers with various shapes termed countour, semiplume, filoplume, bristle and down feathers that cover the body surface to assist flight and mating, and protect against the insults from the environment like cold and water (Strong, 1916; Anderson, 1975; Prum and Brush, 2002; Th and Currie, 2015; Chen et al., 2016; Ruxton et al., 2017). The feathers are the external compartments of the feather follicles which are skin appendages and represented by the primary and secondary feather follicles (Anderson, 1975; Chen et al., 2012b; Liu et al., 2018). The primary feather follicles are universally present in the skin of feathered birds, while the secondary feather follicles are only observed in specific birds like duck and goose. The morphogenesis and development of the feather follicles in chicken are extensively studied to reveal at least five developmental phases, termed macro-patterning, micro-patterning, intra-bud morphogenesis, follicle morphogenesis and regenerative feather cycling growth that are regulated by complex signaling pathways such as wingless-type MMTV integration site family members (Wnt), bone morphogenetic protein (BMP), sonic hedgehog (Shh), and fibroblast growth factors (FGF) (Lin et al., 2006). Of these processes, the de novo formation of the primary feather follicles which are composed of epidermal placode and dermal condensate is particularly interesting to unveil the principles of organogenesis and regeneration that are governed by the crosstalk among different signaling pathways. The specific location and size of the feather follicles are determined by the mechanical cell movements and the delicate molecular interaction and regulation among a series of key genes and pathways (Wessells, 1962; Novel, 1973; Dhouailly, 1975; Chuong et al., 1996; Jiang et al., 1999; Widelitz and Chuong, 1999; Ting-Berreth and Chuong, 2010; Painter et al., 2017; Ho et al., 2019). Initially, the mesenchymal cells begin to form a gradient cell density from the midline to the lateral edge across the chicken back skin. The morphogenetic stripe that is supposed to contain higher cell density compared with the adjacent regions, is emerged in the midline of the back skin indicating the feather tract field as detected by the homogeneous expression of cell adhesion molecules such as neural cell adhesion molecules (NCAM) (Chuong et al., 1991; Jiang and Chuong, 1992; Chuong et al., 1993; Jiang et al., 1999). The cells in the tract field receive the molecular activators and inhibitors to develop the micro-aggregates and later feather follicle primordia consisting of feather placode and dermal condensate (Jiang et al., 1999). The dynamic balance between the activators (CTNNB1, FGF4, SHH and Noggin) and inhibitors (BMP2, BMP4 and BMP12) to stimulate or repress the formation of feather placode and dermal condensate, is supposed to determine the size and density of the primary feather follicles in chicken back skin (Jung et al., 1998; Noramly and Morgan, 1998; Jiang et al., 1999; Widelitz et al., 2000; Ting-Berreth and Chuong, 2010; Mou et al., 2011). Interestingly, during the induction of primary feather follicles, most of the activators and inhibitors are co-localized in the feather primordium such as activators for feather placode (EDAR) and dermal condensate (SHH and FGF4) (Jung et al., 1998; Houghton et al., 2005; Ho et al., 2019), and inhibitors (BMP2, 4, 7 and 12) that aim to abolish the formation of feather primordium in feather tract (Jung et al., 1998; Noramly and Morgan, 1998; Jiang et al., 1999; Mou et al., 2011). The molecular signals display complicated crosstalk to generate the feather follicle primordia during the induction phase. For example, transforming growth factor beta 2 (TGFβ2) is a potential downstream effector of Shh signaling and mainly produced in the feather placode epithelium. TGFβ2 activates the expression of cell adhesion molecules such as NCAM and tenascin-C (TNC) and inhibits the expression of protein kinase C (PKC) to induce the emergence of dermal condensates in chicken skin (Ting-Berreth and Chuong, 1996). Fibroblast growth factor 2 (FGF2) is an early epidermal placode signal to mediate the cell migration in dermis and to induce the occurrence of dermal condensate (Song et al., 1996; Song et al., 2010), whereas fibroblast growth factor 10 (FGF10) that is enriched in mesenchymal cells stimulates the epidermal development by interacting with signaling pathways like BMP and β-catenin (Tao et al., 2002). Importantly, fibroblast growth factor 20 (FGF20) is proposed to be a core signal that is coordinated with cell movements to compress epidermis and induce cell aggregation during the feather follicle induction phase (Ho et al., 2019). EDA that is mainly localized in the inter-placode epidermis, interacts with CTNNB1 and its receptor EDAR to initiate the epidermal expression of FGF20 and attenuate the threshold of the mesenchymal cells with gradient

density to form the cell aggregation and later feather placode and bud (Houghton et al., 2005; Drew et al., 2007; Ho et al., 2019). Chicken is an animal model frequently used to investigate the development and regeneration of feather follicles owing to the existence of one type of feather follicles termed primary feather follicles in chicken back skin (Noramly et al., 1999; Prum, 1999; Yu et al., 2004; Houghton et al., 2005; Lin et al., 2006; Drew et al., 2007; Chen et al., 2015; Cheng et al., 2018; Lin and Yue, 2018). Many other birds, like goose and duck, have two types (primary and secondary) of feather follicles. The primary feather follicles develop different appearances of feathers on the body skin of chicken and goose, whereas the secondary feather follicles produce mainly the down feathers in specific regions of goose skin. Hence, goose may be used as an additional model to further understand the complexity and diversity of the feathers in nature. The morphogenesis and development of the primary feather follicles in chicken skin are profoundly studied, whereas the knowledge of feather follicles in goose skin remains largely unknown. It would be interesting and essential to study the principles governing the occurrence and arrangement of primary feather follicles in goose skin and to further understand the varieties of feathers in birds. Previously, few reports investigate the regulatory network functioned in the formation and elongation of the feather follicle buds in goose skin (Xu et al., 2007; Wu et al., 2008a; Wu et al., 2008b; Chen et al., 2012b; Liu et al., 2018; Sello et al., 2019), while the very early stage, termed the initiation or induction of the primary feather follicles in goose skin is not explored. In this study, the de novo transcriptome sequencing was performed to study the dynamic molecular changes during the primary feather follicle induction in goose skin. The enriched data were compared with those of published data in chicken, mouse and sheep back skin during the primary feather/hair/wool follicle induction stage and gained a handful of candidate genes and signaling pathways, respectively. A total of six Wnt signaling pathway genes were selected to evaluate the spatial-temporal gene expression patterns in the feather primordium formation. The outcome of this report would provide the basic knowledge regarding the feather follicle development in waterfowl, and add extra information to further understand the primary feather follicle development beyond the chicken studies. 2. Materials and methods 2.1. Ethics statement All animal experiments were performed in accordance with the “Guidelines for Experimental Animals” of China council on animal. The experimental procedure was approved by the Standing Committee of Hubei People’s Congress and the ethics committee of Huazhong Agricultural University. 2.2. Sample collection and processing The fertilized eggs of goose (Wanxi White) and chicken (White Leghorn) were collected from the goose breeding farm (Anhui) and the local farm (Wuhan) in China, respectively. The eggs were incubated in egg incubator under certain conditions of temperature (37.8C) and humidity (～65%). The embryos were then sampled from E8.5 to E11.5 (goose) and from E6.5 to E9.5 (chicken) for fixation in 4% paraformaldehyde at 4C, or for dissection the back skin to be stored in TRIzolTM Reagent (Invitrogen, USA) at -80C for further use. 2.3. RNA purification, cDNA library construction and sequencing The total RNA was extracted from TRIzolTM Reagent following the standard protocol. The quality and integrity of RNAs were evaluated with NanoPhotometer® spectrophotometer (IMPLEN, USA), and the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA), respectively. The RNAs of five independent goose back skin at E8.5 or E10 were pooled to prepare two cDNA libraries and run on the Illumina HiSeq 2000 sequencing system (San Diego, USA). Briefly, the two sequencing libraries were prepared from total RNA using the NEBNext R UltraTM Directional RNA Library Prep Kit for Illumina R (NEB, United States). The products were purified with an AMPure XP system (Beckman Coulter, USA) and gone through the quality assessment with an Agilent Bioanalyzer 2100 system in Novogene (Beijing, China). 2.4. De Novo transcriptome reconstruction and gene annotation

The clean reads were applied for de novo assembly by Trinity software (Grabherr et al., 2011). The redundancy sequences were removed by using hierarchical cluster analysis with Corset (Davidson and Oshlack, 2014). The longest transcripts of each cluster (cluster sequence) were selected as unigenes for subsequent analysis. All the unigenes were subsequently applied to search the homologous sequences and annotate the gene function by using the BLASTX program against seven public databases, including Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG (euKaryotic Ortholog Groups), Swiss-prot (A manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes), and GO (Gene Ontology). The software and parameters used in each database are as follows: NR/NT/SwissProt (NCBI blast 2.2.28+, e-value = 1e-5), PFAM (hmmscan package, HMMER3.0, e-value=0.01), KOG (NCBI blast 2.2.28+, e-value = 1e-3), KEGG (KAAS, r140224, e-value= 1e-10), and GO (blast2go, b2g4pipe_v2.5, e-value = 1e-6). 2.5. Differential expression and enrichment analyses of the unigenes The unigenes were mapped to de novo assembled transcriptome by using RSEM software (Li and Dewey, 2011) and edgeR program package through one scaling normalized factor to quantify the read counts of each sequenced library (Robinson et al., 2010) and followed with the differential expression analysis of two libraries by using the DEGseq R package (Anders and Huber, 2010; Wang et al., 2010). P-value was adjusted as qvalue (Storey and Tibshirani, 2003). The threshold of adjusted qvalue < 0.005 and |log2FoldChange| >1 was set up to enrich the differentially expressed transcripts or genes (DEGs) with significance. The DEGs were gone through GO and KEGG pathway enrichment analyses by using the GOseq R packages based on Wallenius non-central hyper-geometric distribution (Young et al., 2010) and KOBAS software (Mao et al., 2005) to adjust gene length bias and the significance of DEGs, respectively. The Gallus gallus in STRING database (http://string-db.org/) was used to analyze the protein-protein interaction (PPI) of DEGs and visualized in Cytoscape (Shannon et al., 2003). The flow chart of the bioinformatics application was shown in Fig. S1. 2.6. Validation of candidate genes with quantitative real-time PCR (qRT-PCR) A total of 6 genes related to the development of feather and hair follicles including WNT16, DKK1, SFRP1, FRZB, FZD1 and LEF1 were selected and confirmed the expression tendency by qRT-PCR. Total RNA of goose embryonic back skin was extracted and reverse transcribed by using PrimeScript™ RT reagent Kit (TaKaRa, Japan). The qRT-PCR analysis was performed by using iTaqTM Universal SYBR® Green Supermix and the Roche LightCyclerR 96 Real-Time Detection System (Bio-Rad, United States). The primers used for qRT-PCR were listed in Table 1. The amplification procedures were 95 C for 6 min initially, followed by 45 cycles of 95 C for 15 s and 60 C for 1 min. Quantification of mRNAs were performed using the 2^-ΔΔCt method with the average cycle thresholds (Ct). The qRT-PCR data were generated from seven independent samples per embryonic stage, and statistically analyzed using Student’s t-test (n = 7). The goose GAPDH was chosen as the internal control. Table 1. The list of primers for qRT-PCR validation.

Genes

Forward (5'–3')

Reverse (5'–3')

WNT16

CACATCTGAACTCCACACTA

ACCTACCGACCACATTGA

LEF1

AGTCTGTCTAGGATGTCATTC

CTGGACCTTGTGATAGTGAA

FZD1

AAGCACGATGGCACCAAGA

AATAGCAGGCAATGACGATGG

DKK1

ACGGACAGATGGGACAGA

GAGTATTCACAACACAGTAACG

SFRP1

AGCCATTGTTGAGCATCT

CTTCTTCAGGTTCTTCTTCC

FRZB

TGTAGGAATAGACTGACTTACC

TGTAGGCACAGACTGAAGA

2.7. Whole mount in situ hybridization

The linearized plasmids containing specific sequences of selected genes were used as templates for in vitro transcription to generate the desired antisense RNA probes labeled with Digoxigenin (Roche, United States). Whole mount in situ hybridization was applied to E8.5-E11.5 goose and E6.5-E9.5 chicken embryos as described previously (Mou et al., 2011). The primers for probe preparation are listed in Table 2. Table 2. The list of primers for PCR amplification.

Genes

Forward (5'–3')

Reverse (5'–3')

WNT16

GGAAGACCTGCTGTATGTG

GGTGATGACTGGTTGACTC

LEF1

CAAGAACATCCGCACAAC

CAGAGGTGGAGTGAAGAC

FZD1

GCCCATCATCTTCCTCTC

GTGAGCCTGGTGTAGAAC

DKK1

CAGCTCCAACGCCATCAAG

GCAGCACCCTCGGAATAATC

SFRP1

TCGGAAGCCATTGTTGAG

CAGCCTGAAGACACCATAA

FRZB

CGCTCCAGGTTACTCTTG

TATGTAGGCACAGACTGAAG

3. Results 3.1. Morphological characterization of primary feather follicle induction in embryonic goose skin The dynamic changes of primary feather follicle morphogenesis in goose back skin were assessed in goose embryos (E8.5-E11.5) (Fig. 1A-D), and compared with those of chicken embryos (E6.5-E9.5) (Fig. 1E-H) by detecting the expression patterns of feather follicle marker gene β-catenin (CTNNB1) with whole mount in situ hybridization technique.

Fig. 1. Morphological identification of the induction of primary feather follicles in goose and chicken back skin. The morphogenesis of primary feather follicles in goose skin is later than that of chicken skin revealed by the positive signals of primary feather placode marker gene β-catenin (CTNNB1) in goose (A-D) and chicken (E-H) back skin using whole mount in situ hybridization. (A-D) The expression of CTNNB1 mRNA in the back skin of goose embryos reveals the homogenous signals in the pre-induction of primary feather follicles in (A) at E8.5, faint signals in the primary feather primordia in (B) at E9.5, the clear feather follicle primordia in (C) at E10.5 and the short feather buds in (D) at E11.5, respectively. (E-H) The developmental patterns of primary feather follicles in chicken back skin at four different stages are shown in (E) as weak primary feather follicle primordia at E6.5, in (F) as clear feather follicle primordia at E7.5, in (G) with the proximal-distal (P-D) axis of feather buds at E8.5, and in (H) with the anterior-posterior (A-P) axis of the elongated feather buds at E9.5.

Two rows of positive hybridization signals of CTNNB1 indicating the location of the prospective primary feather follicles, were firstly detected in goose embryonic back skin at E8.5, displaying nearly 2 days delay compared with the induction of primary feather follicles in chicken embryos at E6.5. The relatively

fuzzy primary feather follicle primordia were developed in goose back skin at E9.5, followed with approximately 10 rows of follicle primordia spreading across the midline of the back skin at E10.5, and exhibiting the equivalent occurrence of feather follicles in chicken embryos at E7.5. 3.2. Transcriptome sequencing, de novo assembly and functional annotation of goose embryonic skin involved in primary feather follicle induction In order to elucidate the key molecular mechanism regulating the induction of goose primary feather follicles, the back skin of E8.5 and E10 goose embryos corresponding to the pre- (E8.5) and induction (E10) of primary feather follicle placodes was applied for cDNA library preparation and Illumina sequencing to generate 57,716,412 (E8.5) and 43,901,286 (E10) raw reads, respectively. A total of 54,858,390 and 41,604,960 clean reads were produced after quality control and applied for the de novo assembly to gain 356,854 transcripts with the length ranging mainly from 200 to 24,000 nucleotides (nt). Of those, 290,582 (81.4%) and 66,272 transcripts (18.6%) were in the length range of 200 to 1000 nucleotides (nt), and >1000 nucleotides (nt), respectively (Fig. 2A). The longest transcript of each gene was used as unigene for subsequent analyses such as annotation.

Fig. 2. De novo transcriptome reconstruction and gene annotation. (A)The length distribution of the assembled unigene transcripts. X-axis: The length of unigenes; Y-axis: The number of unigenes in each interval of length range. (B) The number and percentage of annotated unigenes compared with the listed seven databases. X-axis: The different databases; Y-axis: The number of annotated unigenes. (C) The venn diagrams of annotation compared with the listed five databases. (D) The species classification of the assembled unigenes against the Nr database. (E) The classification of annotated unigenes against Gene Ontology (GO) database. X-axis: GO term; Y-axis: The number of unigenes. (F) The classification of unigenes successfully annotated against KO database by using Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways. X-axis: The percentage of unigenes annotated to the different metabolic pathways; Y-axis: The name of the KEGG pathway. (G) The

classification of successfully annotated unigenes against euKaryotic Ortholog Groups (KOG) database. The terms are divided into 26 groups. X-axis: The 26 KOG groups; Y-axis: The percentage of unigenes.

All assembled 299,689 unigenes were aligned with 7 databases including Nr, Nt, Pfam, KOG, Swiss-prot, KO, and GO databases by using BLASTX program. Briefly, 11,378 unigenes were successfully annotated in seven databases, with the mapping of 25,880 (8.63%) in the Nr database, 106,508 (35.53%) in the Nt database, 14,762 (4.92%) in the KO database, 28,499 (9.5%) in the Swiss-Prot database, 38,358 (12.79%) in the Pfam database, 38,687 (12.9%) in the GO database, and 15,569 (5.19%) in the KOG database (Fig. 2B and 2C), respectively. In detail, the alignment of Nr database revealed that 28.8% of the unigenes were highly homologous with waterfowl Mallard ducks (Anas platyrhynchos) (Fig. 2D). The analyses of GO terms (Lev2) of the unigenes were highly enriched in several biological processes including cellular process, metabolic process and single-organism process (Fig. 2E). The annotated KEGG pathways of the unigenes were mostly involved in cellular processes, environmental information processing and organismal systems (Fig. 2F). In addition, the alignment of unigenes with KOG database was mainly categorized into signal transduction mechanisms, general function prediction only, and posttranslational modification (Fig. 2G). 3.3. Enrichment of differentially expressed genes in goose embryonic skin involved in primary feather follicle induction All the annotated genes were applied to enrich 409 differentially expressed genes (DEGs) including 207 up-regulated and 202 down-regulated genes in goose back skin between the two stages of primary feather follicle induction (E8.5 vs E10) with the standards of qvalue < 0.005 and |log2FoldChange| >1 (Fig. 3A and Table S1).

Fig. 3. De novo RNA sequencing analyses of the goose back skin involved in the induction of primary feather follicles. (A) Volcano plot of differentially expressed genes (DEGs) enriched in the induction of primary feather follicles between two stages (E8.5 vs E10). (B) The Gene Ontology (GO) classification of DEGs. (C) The

enriched top20 KEGG pathways of DEGs. (D) The interaction network of DEGs potentially involved in the induction of primary feather follicles.

The DEGs were further compared with the signature genes enriched in different compartments of E14.5 mouse back skin representing the occurrence period of hair follicle placode and associated dermal condensate (Sennett et al., 2015). A total of 53 genes were shared between mouse and goose back skin during the induction of primary hair and feather follicles (Table 3). Of those, a group of 4 and 12 genes were associated with the development of hair follicle placode and dermal condensate, respectively, whereas 5, 17, 5 and 10 genes were involved in the development of epidermis, dermal fibroblasts, melanocytes and Schwann cells in skin, respectively. Table 3. The listed candidate genes corresponding to different compartments of the goose embryonic skin compared with dataset derived from the back skin of E14.5 mice (Sennett et al., 2015). Pc: Placode; Epi: Epidermis; Mc: Melanocyte; DC: Dermal Condensate; Df: Dermal Fibroblasts; Sch: Schwann cells.

Compartments of skin

Number

Genes

Pc

4

BNC1, LRP4, NCKAP5, RASGEF1B

Epi

5

ABHD14B, CYP2W1, METRNL, PLA2G4B, S100A14

Mc

5

ATP11A, CCK, DMXL2, PMEL, SLC7A8

DC

12

BOK, BTG1, FGF14, FRZB, GLI1, HS3ST6, LGALS1, PGBD5, PLK2, PROKR2, S100A6, THY1

Df

17

C1QTNF3, COL12A1, COL16A1, COL3A1, COL5A1, FAT4, FZD1, ITGBL1, KCNE4, KDELC2, KERA, OGN, PPAP2B, RSPO1, SEC24D, SRD5A2, STRA6

Sch

10

FABP7, FOXD3, LIMCH1, MOXD1, MTSS1, NID2, OLFML2A, PPP1R18, SYNM, TIMP3

The DEGs were further compared with datasets derived from the induction stage of primary feather or wool follicles in chicken or sheep back skin to gain 25 and 19 overlapped genes, respectively (Table 4). Of the 25 genes enriched between chicken and goose comparison, there were 3 marker genes categorized for dermal condensate (S100A6), dermal fibroblast (COL6A1) and Schwann cells (MTSS1). Of the 19 genes enriched in skin between the induction of primary feather and primary wool follicles, a total of 6 marker genes were characterized for placode (LRP4, BNC1 and RASGEF1B), dermal condensate (GLI1), melanocyte (PMEL), and Schwann cells (LIMCH1). Table 4. The differentially expressed genes (DEGs) shared between the induction of primary feather/wool follicles in goose, chicken and sheep back skin (Gong et al., 2018; Nie et al., 2018).

The induction of feather/wool follicle Goose vs Chicken

Number

Genes

25

AKAP9, ALDH1L2, COL3A1, COL6A1, DACT2, EBF3, EIF4B, GPX3, HBZ, ITGAV, JARID2, KRT14, LLGL1, LUC7L3, MTSS1, PORCN, RPLP0, RPLP1, RPS19, RPS27, RPS6, S100A6, SMO, TERF1, TUBA1B

Goose vs Sheep

BNC1, CPED1, DCX, GLI1, KRT14, LEF1, LHFP, LIMCH1, LRP4, MFAP4, PLEKHA5, PMEL, RASGEF1B, RPS23, SFRP1, SULF1, TENM3, TP63, WNT16

19

3.4. GO term and KEGG pathway enrichment analyses of the differentially expressed genes (DEGs) Further details revealed that enriched GO terms of the differentially expressed genes (DEGs) were connected with cell adhesion (18 genes enriched), epithelial cell development (4 genes enriched), Wnt signaling pathway (5 genes enriched), single-organism developmental process (18 genes enriched), regulation of embryonic development (2 genes enriched) and methylation (4 genes enriched) (Fig. 3B; Table 5). Most of the enriched genes displayed upregulated expression in E10 goose back skin. The KEGG pathway analyses showed that the DEGs were enriched in different pathways relevant to skin and feather follicle development (Fig. 3C; Table 6) including cellular community (basal cell carcinoma, signaling pathways regulating pluripotency of stem cells, adherens junction and focal adhesion), signal transduction (Wnt, Hedgehog and Hippo signaling pathways), and signaling molecules and interaction (ECM-receptor interaction and MAPK signaling pathways). The enrichment of Wnt signaling pathway in both KEGG pathway and GO term analyses indicated that the Wnt pathway genes could be potential candidates involved in the primary feather follicle induction in goose back skin. Table 5. The enriched GO terms of the DEGs potentially correlated with the induction of primary feather follicles in goose skin. The setting bold fonts stick out the up-regulated DEGs.

Description

Number

Genes

Cell adhesion

18

COL6A1, EFNB1, FAT1, FAT4, FLNA, GOLGB1, ICE1, LAMB4, LY75, NCAN, NID2, PIK3R2, PTK7, SLIT1, SMC4, TRIM62, VCAN, VIM

Epithelial cell development

4

APOA1BP, CENPF, FLNA, SAMD11

Wnt signaling pathway

5

APC, DKK1, FZD1, TIMP3, WNT16

Single-organism developmental process

18

APOA1BP, ARHGEF3, CENPF, COL5A1, DBN1, DKK1, EFNB1, FLNA, ITGAV, LAMB4, PALM, PIK3R2, RNF144A, SAMD11, SREK1, TENM3, TPR, WNT16

Regulation of embryonic development

2

LAMB4, PIK3R2

Methylation

4

DMD, ICE1, PRMT5, STARD9

Table 6. The enriched KEGG pathways of the DEGs related to the induction of primary feather follicles in goose skin. The setting bold fonts stick out the up-regulated DEGs.

Signaling Pathway

Number

Genes

Basal cell carcinoma

7

APC, FZD1, GLI1, LEF1, PTCH2, SMO, WNT16

Wnt

8

APC, DKK1, FZD1, LEF1, JNK/MAPK8, PORCN, SFRP1, WNT16

Hedgehog

4

GLI1, PTCH2, SMO, WNT16

Signaling pathways regulating pluripotency of stem cells

7

APC, FGFR1, FZD1, ID4, JARID2, PIK3R2, WNT16

Hippo signaling pathway - fly

5

FAT4, LLGL1, LIX1L, JNK/MAPK8, MPDZ

ECM-receptor interaction

5

COL3A1, COL6A1, COL6A3, HMMR, ITGAV

Adherens junction

4

FGFR1, IQGAP1, LEF1, WASF1

Focal adhesion

8

COL3A1, COL6A1, COL6A3, FLNA, ITGAV, JNK/MAPK8, PDGFA, PIK3R2

MAPK

7

FGFR1, FLNA, JNK/MAPK8, PDGFA, PLA2G4A, PLA2G4B, TAOK1

3.5.The mRNA and mRNA interaction network of the DEGs involved in goose primary feather follicle induction The differentially expressed genes (DEGs) were selected for protein-protein interaction (PPI) network construction (Fig. 3D) to establish several small networks regulating cell proliferation (ASPM, MKI67, KNTC1 and CENPF) (Falini et al., 1989; Zhou et al., 2005; Williams et al., 2015; Liu et al., 2019), and hair and skin development (ITGAV, OCN and PDGFA) (Ge et al., 2004; Wang et al., 2015; Rivera-Gonzalez et al., 2016; Hu et al., 2018; Pazzaglia et al., 2019). The altered genes were also highly correlated with the development of skin and feather or hair follicles represented by the Wnt (APC, DKK1, FZD1, LEF1, SFRP1 and WNT16), Hedgehog (GLI1, PTCH2, SMO and WNT16), and Hippo (FAT4, LLGL1 and MPDZ) signaling pathways. 3.6. Validation of candidate genes functioned in the development of primary feather follicles in goose skin by quantitative real-time PCR (qRT-PCR) A total of 6 genes enriched in the DEGs were selected to validate the expression tendency by qRT-PCR during the induction and development of the primary feather follicles in goose back skin at four different embryonic stages with seven independent replicates per stage (n=7) (Fig. 4). Of those, WNT16 and DKK1 shared the similar expression tendency to have the lowest expression level at E8.5, followed with strong up-regulation at E9.5 and E10.5, and down-regulation at E11.5 (Fig. 4A and 4B). SFRP1 and FRZB both showed continuously increased expression in back skin from E8.5 to E11.5 and displayed the weakest expression at E8.5 and the sharp elevation from E10.5 to E11.5 (Fig. 4D and 4E). The expression of FZD1 and LEF1 were both up-regulated in goose back skin from E8.5 to E11.5 with slight differences. FZD1 exhibited gradually increased expression from E8.5 to E11.5, while LEF1 showed strong up-regulation from 8.5 to E9.5, followed with gradually increased expression from E9.5 to E11.5 (Fig. 4C and 4F). All detected six Wnt pathway genes showed clearly stimulated expression during the induction of primary feather follicles in goose embryonic back skin staged from E8.5 to E10.5. During the bud stage from E10.5 to E11.5, four genes (LEF1, SFRP1, FRZB and FZD1) showed continuously increased expression pattern, while other 2 genes (WNT16 and DKK1) display decreased expression. Fig. 4. The dynamic expression of selected Wnt pathway genes is identified by qPCR in goose embryonic back skin at 5 developmental stages of the primary feather follicles. The listed genes were WNT16 (A), DKK1(B), LEF1 (C), SFRP1 (D), FRZB (E) and FZD1 (F) detected in goose skin from E8.5 to E11.5 (n=7). All six genes are up-regulated during the induction of primary feather follicles from E8.5 to E10.5. During the bud stages from E10.5 to E11.5, four genes (LEF1, SFRP1, FRZB and FZD1) are up-regulated, while other 2 genes (WNT16 and DKK1) are down-regulated. X-axis: The developmental stages of goose embryos; Y-axis: Relative expression of genes of interest. The significant level of differences is indicated as * (p value≤0.05); ** and *** (p value≤0.01).

3.7. Detection of selected Wnt pathway genes in goose embryonic skin by whole mount in situ hybridization Whole mount in situ hybridization was applied to further define the dynamic expression patterns of selected six Wnt signaling pathway genes termed WNT16, FZD1, LEF1, DKK1, SFRP1 and FRZB during the induction of primary feather follicles in goose embryos from E8.5 to E11.5 (Fig. 5A-F).

Fig. 5. The expression patterns of six selected Wnt pathway genes in 4 developmental stages of primary feather follicles in goose back and neck skin by using in situ hybridization. The expression patterns of WNT16 (A);

FZD1 (B); LEF1 (C); DKK1(D) , SFRP1(E) and FRZB (F) in the back and neck skin are all localized in the primary feather follicles. FZD1, LEF1 and DKK1 show stronger expression in both back and neck skin compared with those of WNT16, SFRP1 and FRZB genes. Scale bar: Back, 200 μm; Neck, 500 μm.

WNT16, functioned as one of the Wnt ligands, was gradually increased the expression level in the back and neck skin from E8.5 to E11.5, with the faint signals across the skin at E8.5 and the restricted positive signals in the feather follicle at E10.5 and E11.5 both in back and neck skin (Fig. 5A). The expression of FZD1 was nearly undetectable in E8.5 back skin and was focally restricted in the feather follicle primordia of the midline of back skin at E9.5, and later in the feather follicles in E10.5 and E11.5 back skin. FZD1 showed the homogenous expression at E8.5 and the strongest expression at E11.5 in the feather follicles of neck skin (Fig. 5B). The positive signals of LEF1 were firstly evident with two strong rows in the middle line of E8.5 back skin, followed with clearly specified signals in the feather follicles and fuzzy signals in the inter-follicular place of the back skin at E9.5, and later with mostly focused signals in the feather follicles at E10.5 and E11.5 (Fig. 5C). The expression of LEF1 in the neck showed similar expression patterns with those of the back skin displaying strong expression in the presumptive feather tract area at E8.5 and focal expression in the feather follicles at E11.5. Though DKK1, SFRP1 and FRZB are Wnt inhibitory molecules, they showed different expression patterns during the primary feather follicle induction in goose skin (Fig. 5D-F). DKK1 were detected with strong and broad expression in both back and neck skin at E8.5, followed with gradually restricted signals in the primary feather follicle primordia and moderate expression in the skin area at E9.5, E10.5 and E11.5 (Fig. 5D). It is interesting that the signals in neck skin were much stronger than those of back skin as observed both in early and later stages. SFRP1 and FRZB both showed unspecified weak expression in the back and neck skin at E8.5 and E9.5, while the fairly stronger positive signals were detected in the feather follicles both in the back and neck skin at E10.5 and E11.5 (Fig. 5E-F). All the detected six genes were expressed in the primary feather follicles, whereas the positive signals of WNT16, SFRP1 and FRZB in both neck and back skin were much weaker compared with those of LEF1, FZD1 and DKK1 genes. 4. Discussion Goose is a domesticated waterfowl and well known for the good quality of the feathers that are developed from the feather follicles in skin. Compared with goose, chicken is widely used to investigate the morphogenesis and development of the skin and feather follicles owing to the unique presence of primary feather follicles in chicken back skin. It is not clear whether the regulatory mechanisms underlying primary feather follicle development between goose and chicken are conserved or not. In this report, the preliminary comparison between goose and chicken back skin was conducted to figure out the similarities and differences regulating the induction of primary feather follicles. The morphogenesis of the primary feather follicles in goose skin was staged between E8.5 to E11.5 by in situ hybridization of feather placode marker gene CTNNB1. Combined with the previous reports, it suggests that the primary feather follicles begin with the development of epidermal placode and associated dermal condensate, which is occurred in embryonic back skin at E10-E10.5 in goose (Xu et al., 2007; Wu et al., 2008a; Wu et al., 2008b; Chen et al., 2012b), and at E7-E7.5 in chicken (Hughes et al., 2011; Chang et al., 2015; Wu et al., 2017). During this process, the similar morphological changes are observed between chicken and goose skin (Wu et al., 2008a), including the de novo formation of the placode and dermal condensate, and the thickening of the epidermis and dermis. The molecular changes were explored to reveal the regulatory network governing the induction of the primary feather follicles in goose skin in this study, as well as the primary hair follicles in mouse skin, and the primary wool follicles in sheep skin using RNA-sequencing strategy as reported previously (Sennett et al., 2015; Nie et al., 2018). The DEGs gained in this study were compared with the signature genes enriched in the induction of primary hair follicles in mouse back skin (Sennett et al., 2015), primary feather follicles in chicken back skin (Gong et al., 2018) and primary wool follicles in sheep fetal skin (Nie et al., 2018) to obtain 53, 25 and 19 overlapped genes, respectively (Table 3 and Table 4). Of these, a total of 6 genes were shared in back skin among goose, mouse and sheep, represented by marker genes for placode (RASGEF1B, LRP4 and BNC1), for dermal condensate (GLI1), for melanocytes (PMEL), and for Schwann cells (LIMCH1). It is not surprising that the comparison of DEGs enriched in the induction of primary hair/feather follicles between goose and mouse back skin gained the most overlapped marker genes for six compartments of the skin including LRP4 for

placode, S100A14 for epidermis, FGF14, GLI1 and S100A6 for dermal condensate, and COL12A1, COL16A1, COL3A1, COL5A1, FAT4, FZD1 and ITGBL1 for dermal fibroblast. It highly suggests that from E8.5 to E10 or E10.5, the embryonic goose skin experiences great histologic changes including the thickening of epidermis and dermis, and the de novo formation of the placode and dermal condensate that are marked by the signature genes and driven by the dynamic molecular changes. It also indicates that the molecular networks regulating the initial development of the primary feather/hair/wool follicles share common features and possess species-specific characteristics to make the delicate morphological differences of the skin or follicles in different animal models. The enriched candidate genes are involved in Wnt (LRP4), FGF (FGF14) and Shh (GLI1 and SMO) signaling pathways. Of these genes, the LRP4 gene is a low-density lipoprotein (LDL) receptor-related protein, which differs from LRP5/6 co-receptors of Wnts in the intracellular domain and generally functions as negative regulator for the Wnt/β-catenin signaling pathway (Simonchazottes et al., 2006; Fliniaux et al., 2008). LRP4 is expressed in the placodes of skin appendages such as mammary glands and hair follicles in murine skin (Weatherbee et al., 2006; Fliniaux et al., 2008; Ahn et al., 2013) and is one of the regulators controlling the initiation and patterning of skin appendage placodes (Ahn et al., 2013). In LRP4 mutant skin, the primary hair follicle placodes are initiated approximately one day later at E14.5 compared with those of control mice at E13.5. It suggests that the morphogenesis of primary hair follicles is delayed rather than abolished in LRP4 mutants and LRP4 potentially interacts with Wnt signaling components to facilitate the primary hair placode initiation (Ahn et al., 2013). The enriched LRP4 in our data indicates that LRP4 is potentially involved in the initiation of primary feather placodes in goose skin. GLI1 is a transcription factor and a transcriptional target of Shh signaling which is essential during the morphogenesis of the hair follicle, particularly the progression of hair peg down-growth into the skin dermis. GLI1 is expressed in murine back skin with higher levels in dermal condensate at E14.5 and later in dermal papilla at E18.5, and with lower levels in the epidermis and dermis of inter-follicular regions (St-Jacques et al., 1998; Chiang et al., 1999). Smoothened (SMO) is another Sonic Hedgehog (Shh) signaling effector that is broadly expressed in murine skin. The targeted knock out of SMO in ventral dermis impairs hair follicle generation with reduced follicle number and arrested progressive stages from E17.5 to P0 in murine ventral skin (Woo et al., 2012; Tomoko et al., 2013). GLI1 and SMO are candidate genes involved in the formation of dermal condensate of the primary feather follicles in goose skin. A handful genes were grouped and correlated with the development of epidermis and dermal fibroblast like collagen genes. The telomere repeat-binding factor 1 (TERF1 or TRF1) conditionally mutant mice show severe skin epithelial abnormalities with the immature hair follicles, impaired epithelial barrier and skin hyperpigmentation (Martínez et al., 2009; Muñoz et al., 2009; Brossard et al., 2015), which are associated with the telomere dysfunction (Blasco et al., 1997; Lange, 2005; Okamoto et al., 2008). COL3A1 belongs to type III collagens that are present in elastic tissues like skin. The mutation of COL3A1 results in the malformed collagens with thin and translucent skin, referred to vascular Ehlers-Danlos syndrome (vEDS) (D'Hondt et al., 2018). Integrin alpha v (ITGAV) encodes one of the integrin subunits that are glycoproteins attaching cells to the extracellular matrix or to the ligands. The transgenic mice overexpressing ITGAV under the enhancer and promoter of COL1A2 exhibit thinner dermis combined with decreased protein and mRNA levels of COL1A2, COL3A1, CTGF and integrin beta 3 (ITGB3) in the skin (Wang et al., 2015). These genes are implicated to function in the formation of epidermis and dermis in goose embryonic skin. The enriched genes correlated with cell adhesion, adherens junction, focal adhesion, Wnt, Hedgehog, and Hippo signaling pathways in KEGG and GO analyses indicate that the cells in the goose embryonic back skin (E8.5 vs E10) are modulated to communicate, proliferate, and differentiate to shape the individual skin compartments from the homogeneous skin layers directed by the mechanical cell movements and molecular signals at the same time, as suggested during the induction of primary feather follicles in chicken skin or primary wool follicles in sheep skin (Gong et al., 2018; Nie et al., 2018). The Wnt/β-catenin signaling pathway genes were extensively studied in the development of feather and hair follicles in chicken and mouse models (Chen et al., 2012a; Cui et al., 2014). CTNNB1 is active in the early patterning of the feather tract to trigger the morphogenesis of feather buds (Noramly et al., 1999; Widelitz et al., 2000; Wu et al., 2008b). In this study, CTNNB1 was used as a marker gene to detect the early development of primary feather follicles in goose skin, which shows the similar expression patterns with those of chicken skin (Fig. 1A-H). Interestingly, several Wnt pathway genes were enriched in the goose back skin represented by APC, DKK1, FZD1, LEF1, SFRP1 and WNT16, and in the chicken, goose, mouse and sheep back skin during

the primary feather/hair/wool follicle development such as APC, SFRP1 and WNT16 (Gong et al., 2018; Nie et al., 2018). WNT16 acts as a Wnt ligand of the canonical Wnt/β-catenin signaling to drive the proliferation and differentiation of keratinocytes, osteoclasts and perivascular stem cells (Clements et al., 2011; Kobayashi et al., 2015; Genthe and Clements, 2017; Nalesso et al., 2017; Shen et al., 2017; Mendoza-Reinoso and Beverdam, 2018; Meyers et al., 2018). WNT16 was weakly expressed in the follicular epithelium at the placode stage in murine skin and potentially involved in primary wool follicle induction (Reddy et al., 2001; Nie et al., 2018). The role of WNT16 in the development of the primary feather follicles remains unclear. In this study, the increased expression of WNT16 was accompanied by the expansion of feather tract and the increased number of feather primordia. WNT16 may be involved in the proliferation of skin epithelial cells, thereby promoting the occurrence of the feather primordia. Interestingly, the expression of WNT16 and DKK1 in embryonic goose skin shared the similar mRNA expression trend (Fig. 4A and 4B), not the spatial expression patterns during the primary feather follicle morphogenesis (Fig. 5A and 5D). DKK1 was extensively up-regulated in goose skin before the feather follicle initiation, and then focally expressed with strong signals in feather primordium and feather bud (Fig. 4B and 5D), displaying the similar spatial patterns of DKK1 in the development of feather follicles in chicken skin (Gong et al., 2018). The role of DKK1 in the development and periodic growth of hair follicles in skin is widely studied in mouse models. The ectopic expression of DKK1 in the skin of transgenic mice results in the failure of hair follicle placode induction and blockage of hair follicle development before the bud stage during the embryonic development (Andl et al., 2002). Moreover, DKK1 promotes the transition of anagen to catagen during the postnatal cycling growth of hair follicles by partially regulating the activity of follicular keratinocytes and modulates the hair follicle size combined with WNT10b during the hair regeneration in mouse models (Lei et al., 2014; Kobayashi et al., 2015; Qiu et al., 2017). Several studies have shown that the overexpression of DKK1 mediated by RCAS-DKK1 suppresses the formation of feather tracts and dermal structure of individual feather buds in the presumptive feather region, as well as the skin thickness by affecting Wnt/β-catenin signaling in keratinocytes (Chang, 2007; Yamaguchi et al., 2008; Yamaguchi et al., 2009). DKK1 may function as a negative regulator for the primary feather follicle development in goose skin as evidenced in mouse and chicken models. The expression of FZD1 and LEF1 were gradually increased spatially from E8.5 to E11.5 (Fig. 4C and 4F), while FZD1 mRNA expression was lower than LEF1. In our whole mount in situ hybridization results, FZD1 had strong signals in the back and neck skin of goose embryos at E11.5, that is the period of emerged proximal-distal (P-D) axis and anterior-posterior (A-P) axis of the feather bud. Interestingly, FZD1 is also expressed in the distal and posterior epithelium at the early stage of feather bud in chicken embryos (Chodankar et al., 2003). In addition, FZD1 is expressed with lower levels in dermis and specifically elevated in hair follicle placodes (Reddy et al., 2004). Accumulating evidences suggest that the homolog members of frizzled gene family are identified to function in planar cell polarity (PCP) regulation in vertebrates (Daudet et al., 2002; Wang et al., 2006; Deans et al., 2007; Seifert and Marek, 2007; Devenport and Fuchs, 2008; Yu et al., 2010; Goodrich and Strutt, 2011). Together, FZD1 may be involved in the occurrence of primary feather primordium, and the formation of P-D and A-P axis of goose feather bud by affecting the planar cell polarity. Previous results indicate that LEF1 regulates the epithelial-mesenchymal interaction in the development of skin appendages, such as hair, tooth and mammary gland (van Genderen et al., 1994; Kratochwil et al., 1996; Li et al., 2002; Liu et al., 2004; Boras-Granic et al., 2006; Wang et al., 2012). These results indicate that LEF1 is a Wnt signaling molecule required for the induction of skin appendage placodes. This is consistent with the observation in this study that LEF1 showed the strong positive signals in primary feather placode and the epidermis of the feather buds, especially in the neck skin (Fig. 5C). The higher expression of LEF1 than other detected genes validated by qPCR (Fig. 4C) indicates that LEF1 is a positive regulator for the primary feather follicle development in goose skin. SFRP1 and FRZB, as inhibitory molecules of Wnt signaling pathway, have similar expression levels and expression trends during the development of primary feather follicles. Secreted frizzled-related protein 1 (SFRP1) which is expressed in the bulge cells and dermal cells of the hair follicles, induces the differentiation of bulge stem cells and inter-follicular epidermal stem cells (Lim et al., 2016; Ichijo et al., 2017), and regulates Wnt/β-catenin activity during hair follicle development (Tripurani et al.; Hawkshaw et al., 2018; Zhang et al., 2018). Studies related to FRZB gene in skin appendages are mainly focused on chickens. In the embryonic period, FRZB is first expressed in the ectodermal layer cells of the neural plate, and subsequently expressed in

the neural tube, neural crest and mesenchyme of limb bud (Duprez et al., 1999; Wada et al., 1999; Baranski et al., 2000). In postnatal period, DKK/FRZB in feather dermal papilla (DP) interacts with Wnt ligand in feather epithelium/pulp to regulate feather regeneration and feather axis formation (Chu et al., 2014). Notably, FRZB gene is identified in other studies to potentially affect fur development in rabbits and feathered-leg trait in chicken (Zhao et al., 2017; Yang et al., 2019). The weak expression of both SFRP1 and FRZB during the development of primary feather primordia in goose back skin indicates that SFRP1 and FRZB may be combined with other Wnt pathway genes to regulate the development of primary feather follicles. 5. Conclusions In summary, the transcriptome analysis revealed the complex gene network directed by signaling pathways like Wnt and Shh regulating the induction of primary feather follicles in goose skin. The comparison of the enriched data in the induction of primary feather/hair/wool follicles in skin between goose and other animal models (chicken, mouse and sheep) demonstrates that the mechanism underlying the induction of primary feather follicles in goose skin is partially conserved with other species and owns species-specific features. Our report for the first time revealed the dynamic molecular changes governing the induction of primary feather follicles in goose back skin, which not only provides the general knowledge regarding the feather follicle and skin development in waterfowl, but also contributes to further understand the diversity of skin and hair or feather development beyond the mouse and chicken models. Supplementary Materials: Figure S1. The flow chart for the bioinformatics process applied to analyze the transcriptome sequencing data in this study. Table S1. The differentially expressed genes (DEGs) between E8.5 and E10 of the induction of primary feather follicles in goose back skin. Author Contributions: Conceived and designed experiments, C.M.; performed experiments, X.Z. (Xiaokang Zhang) and Z.L.; data analysis, X.H., X.Z. (Xiaokang Zhang), X.Z. (Xinting Zheng) and Y.N.; validation, X.H., S.L., and G.Y.; visualization, X.H., X.Z. (Xiaokang Zhang), Y.T. and W.W.; resources, Q.X. and Y.Z.; data curation, X.Z. (Xiaoliu Zhou) and Q.Z; writing—original draft preparation, X.H. and X.Z. (Xiaokang Zhang); writing—revise and editing, X.H. and C.M. Funding: This work was supported by grants from National Natural Science Foundation of China (No. 31972548) and National Key R&D Program of China (2018YFD0501301). Acknowledgments: We thank all the members involved for their valuable suggestions and great help. Conflicts of Interest: The authors declare no conflict of interest.

References Ahn, Y., Sims, C., Logue, J.M., Weatherbee, S.D. and Krumlauf, R., 2013. Lrp4 and Wise interplay controls the formation and patterning of mammary and other skin appendage placodes by modulating Wnt signaling. Development 140, 583-593. Anders, S. and Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biol 11, R106. Anderson, W.D., 1975. Avian Anatomy - Integument, Parts I and II. Canadian Veterinary Journal La Revue Veterinaire Canadienne 16, 96. Andl, T., Reddy, S.T., Gaddapara, T. and Millar, S.E., 2002. WNT signals are required for the initiation of hair follicle development. Dev Cell 2, 643-53. Baranski, M., Berdougo, E., Sandler, J.S., Darnell, D.K. and Burrus, L.W., 2000. The dynamic expression pattern of frzb-1 suggests multiple roles in chick development. Dev Biol 217, 25-41. Blasco, M.A., Lee, H.W., Hande, M.P., Samper, E., Lansdorp, P.M., Depinho, R.A. and Greider, C.W., 1997. Telomere shortening and tumor formation by mouse cells lacking telomerase RNA. Cell 91, 25-34. Boras-Granic, K., Chang, H., Grosschedl, R. and Hamel, P.A., 2006. Lef1 is required for the transition of Wnt signaling from mesenchymal to epithelial cells in the mouse embryonic mammary gland. Dev Biol 295, 219-31. Brossard, M., Fang, S., Vaysse, A., Wei, Q., Chen, W.V., Mohamdi, H., Maubec, E., Lavielle, N., Galan, P. and Lathrop, M., 2015. Integrated pathway and epistasis analysis reveals interactive effect of genetic variants at TERF1 and AFAP1L2 loci on melanoma risk. International Journal of Cancer 137, 1901-1909.

Chang, H.Y., 2007. Patterning skin pigmentation via dickkopf. J Invest Dermatol 127, 994-5. Chang, K.W., Huang, N.A., Liu, I.H., Wang, Y.H., Wu, P., Tseng, Y.T., Hughes, M.W., Jiang, T.X., Tsai, M.H. and Chen, C.Y., 2015. Emergence of differentially regulated pathways associated with the development of regional specificity in chicken skin. Bmc Genomics 16, 22. Chen, C.F., Foley, J., Tang, P.C., Li, A., Jiang, T.X., Wu, P., Widelitz, R.B. and Chuong, C.M., 2015. Development, regeneration, and evolution of feathers. Annu Rev Anim Biosci 3, 169-95. Chen, C.K., Chen, S.N., Wu, S.M., Chen, J.J., Cheng, P.L., Wu, P., Lu, M.Y.J., Chen, D.R., Chuong, C.M. and Cheng, H.C., 2016. Regulatory Differences in Natal Down Development between Altricial Zebra Finch and Precocial Chicken. Molecular Biology & Evolution 33, 2030. Chen, D., Jarrell, A., Guo, C., Lang, R. and Atit, R., 2012a. Dermal beta-catenin activity in response to epidermal Wnt ligands is required for fibroblast proliferation and hair follicle initiation. Development 139, 1522-33. Chen, X., Bai, H., Li, L., Zhang, W., Jiang, R. and Geng, Z., 2012b. Follicle characteristics and follicle developmental related Wnt6 polymorphism in Chinese indigenous Wanxi-white goose. Molecular Biology Reports 39, 9843-9848. Cheng, D., Yan, X., Qiu, G., Zhang, J., Wang, H., Feng, T., Tian, Y., Xu, H., Wang, M. and He, W., 2018. Contraction of basal filopodia controls periodic feather branching via Notch and FGF signaling.

9, 1345.

Chiang, C., Swan, R.Z., Grachtchouk, M., Bolinger, M., Litingtung, Y., Robertson, E.K., Cooper, M.K., Gaffield, W., Westphal, H., Beachy, P.A. and Dlugosz, A.A., 1999. Essential role for Sonic hedgehog during hair follicle morphogenesis. Dev Biol 205, 1-9. Chodankar, R., Chang, C.H., Yue, Z., Jiang, T.X., Suksaweang, S., Burrus, L.W., Chuong, C.M. and Widelitz, R.B., 2003. Shift of Localized Growth Zones Contributes to Skin Appendage Morphogenesis: Role of the Wnt/β-catenin Pathway. Journal of Investigative Dermatology 120, 20-26. Chu, Q., Cai, L., Fu, Y., Chen, X., Yan, Z., Lin, X., Zhou, G., Han, H., Widelitz, R.B. and Chuong, C.M., 2014. Dkk2/Frzb in the dermal papillae regulates feather regeneration. Developmental Biology 387, 167-178. Chuong, C.M., Chen, H.M., Jiang, T.X. and Chia, J., 1991. Adhesion molecules in skin development: morphogenesis of feather and hair. Annals of the New York Academy of Sciences 642, 263-80. Chuong, C.M., Widelitz, R.B. and Jiang, T.X., 1993. Adhesion molecules and homeoproteins in the phenotypic determination of skin appendages. The Journal of investigative dermatology 101, 10S-15S. Chuong, C.M., Widelitz, R.B., Tingberreth, S. and Jiang, T.X., 1996. Early events during avian skin appendage regeneration: dependence on epithelial-mesenchymal interaction and order of molecular reappearance. Journal of Investigative Dermatology 107, 639-646. Clements, W.K., Kim, A.D., Ong, K.G., Moore, J.C., Lawson, N.D. and David, T., 2011. A somitic Wnt16/Notch pathway specifies haematopoietic stem cells. Nature 474, 220-4. Cui, C.Y., Yin, M., Sima, J., Childress, V., Michel, M., Piao, Y. and Schlessinger, D., 2014. Involvement of Wnt, Eda and Shh at defined stages of sweat gland development. Development 141, 3752-60. D'Hondt, S., Guillemyn, B., Syx, D., Symoens, S., De Rycke, R., Vanhoutte, L., Toussaint, W., Lambrecht, B.N., De Paepe, A., Keene, D.R., Ishikawa, Y., Bachinger, H.P., Janssens, S., Bertrand, M.J.M. and Malfait, F., 2018. Type III collagen affects dermal and vascular collagen fibrillogenesis and tissue integrity in a mutant Col3a1 transgenic mouse model. Matrix Biol 70, 72-83. Daudet, N., Ripoll, C., Molès, J.P. and Rebillard, G., 2002. Expression of members of Wnt and Frizzled gene families in the postnatal rat cochlea. Brain Res Mol Brain Res 105, 98-107. Davidson, N.M. and Oshlack, A., 2014. Corset: enabling differential gene expression analysis for de novo assembled transcriptomes. Genome Biol 15, 410. Deans, M.R., Dragana, A., Kaye, S., Scott, M.P., Axelrod, J.D. and Goodrich, L.V., 2007. Asymmetric distribution of prickle-like 2 reveals an early underlying polarization of vestibular sensory epithelia in the inner ear. Journal of

Neuroscience 27, 3139-3147. Devenport, D. and Fuchs, E., 2008. Planar polarization in embryonic epidermis orchestrates global asymmetric morphogenesis of hair follicles. Nat Cell Biol 10, 1257-68. Dhouailly, D., 1975. Formation of cutaneous appendages in dermo-epidermal recombinations between reptiles, birds and mammals. Wilhelm Rouxs Archives of Developmental Biology 177, 323-340. Drew, C.F., Lin, C.M., Jiang, T.X., Blunt, G., Mou, C., Chuong, C.M. and Headon, D.J., 2007. The Edar subfamily in feather placode formation. Dev Biol 305, 232-45. Duprez, D., Leyns, L., Bonnin, M.A., Lapointe, F., Etchevers, H., De Robertis, E.M. and Le, D.N., 1999. Expression of Frzb-1 during chick development. Mechanisms of Development 89, 179-183. Falini, B., Flenghi, L., Fagioli, M., Stein, H., Schwarting, R., Riccardi, C., Manocchio, I., Pileri, S., Pelicci, P.G. and Lanfrancone,

L.,

1989.

Evolutionary

conservation

in

various

mammalian

species

of

the

human

proliferation-associated epitope recognized by the Ki-67 monoclonal antibody. J Histochem Cytochem 37, 1471-8. Fliniaux, I., Mikkola, M.L., Lefebvre, S. and Thesleff, I., 2008. Identification of dkk4 as a target of Eda-A1/Edar pathway reveals an unexpected role of ectodysplasin as inhibitor of Wnt signalling in ectodermal placodes. Developmental Biology 320, 60-71. Ge, G., Seo, N.S., Liang, X., Hopkins, D.R., Hook, M. and Greenspan, D.S., 2004. Bone morphogenetic protein-1/tolloid-related metalloproteinases process osteoglycin and enhance its ability to regulate collagen fibrillogenesis. J Biol Chem 279, 41626-33. Genthe, J.R. and Clements, W.K., 2017. R-spondin 1 is required for specification of hematopoietic stem cells through Wnt16 and Vegfa signaling pathways. Development 144, 590-600. Gong, H., Wang, H., Wang, Y., Bai, X., Liu, B., He, J., Wu, J., Qi, W. and Zhang, W., 2018. Skin transcriptome reveals the dynamic changes in the Wnt pathway during integument morphogenesis of chick embryos. Plos One 13, e0190933. Goodrich, L.V. and Strutt, D., 2011. Principles of planar polarity in animal development. Development 138, 1877-1892. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N. and Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29, 644-52. Hawkshaw, N.J., Hardman, J.A., Haslam, I.S., Shahmalak, A., Gilhar, A., Lim, X. and Paus, R., 2018. Identifying novel strategies for treating human hair loss disorders: Cyclosporine A suppresses the Wnt inhibitor, SFRP1, in the dermal papilla of human scalp hair follicles. Plos Biology 16, e2003705. Ho, W.K.W., Freem, L., Zhao, D., Painter, K.J., Woolley, T.E., Gaffney, E.A., McGrew, M.J., Tzika, A., Milinkovitch, M.C., Schneider, P., Drusko, A., Matthäus, F., Glover, J.D., Wells, K.L., Johansson, J.A., Davey, M.G., Sang, H.M., Clinton, M. and Headon, D.J., 2019. Feather arrays are patterned by interacting signalling and cell density waves. PLoS biology 17, e3000132. Houghton, L., Lindon, C. and Morgan, B.A., 2005. The ectodysplasin pathway in feather tract development. Development 132, 863-72. Hu, X., Li, Y.Q., Li, Q.G., Ma, Y.L., Peng, J.J. and Cai, S.J., 2018. Osteoglycin (OGN) reverses epithelial to mesenchymal transition and invasiveness in colorectal cancer via EGFR/Akt pathway. J Exp Clin Cancer Res 37, 41. Hughes, M.W., Wu, P., Jiang, T.X., Lin, S.J., Dong, C.Y., Li, A., Hsieh, F.J., Widelitz, R.B. and Chuong, C.M., 2011. In search of the Golden Fleece: unraveling principles of morphogenesis by studying the integrative biology of skin appendages. Integrative Biology Quantitative Biosciences from Nano to Macro 3, 388-407. Ichijo, R., Kobayashi, H., Yoneda, S., Iizuka, Y., Kubo, H., Matsumura, S., Kitano, S., Miyachi, H., Honda, T. and

Toyoshima, F., 2017. Tbx3-dependent amplifying stem cell progeny drives interfollicular epidermal expansion during pregnancy and regeneration. Nature Communications 8, 508. Jiang, T.X. and Chuong, C.M., 1992. Mechanism of skin morphogenesis. I. Analyses with antibodies to adhesion molecules tenascin, N-CAM, and integrin. Developmental biology 150, 82-98. Jiang, T.X., Jung, H.S., Widelitz, R.B. and Chuong, C.M., 1999. Self-organization of periodic patterns by dissociated feather mesenchymal cells and the regulation of size, number and spacing of primordia. Development (Cambridge, England) 126, 4997-5009. Jung, H.S., Francis-West, P.H., Widelitz, R.B., Jiang, T.X., Ting-Berreth, S., Tickle, C., Wolpert, L. and Chuong, C.M., 1998. Local Inhibitory Action of BMPs and Their Relationships with Activators in Feather Formation: Implications for Periodic Patterning ☆. Developmental Biology 196, 11-23. Kobayashi, Y., Thirukonda, G.J., Nakamura, Y., Koide, M., Yamashita, T., Uehara, S., Kato, H., Udagawa, N. and Takahashi, N., 2015. Wnt16 regulates osteoclast differentiation in conjunction with Wnt5a. Kratochwil, K., Dull, M., Farinas, I., Galceran, J. and Grosschedl, R., 1996. Lef1 expression is activated by BMP-4 and regulates inductive tissue interactions in tooth and hair development. Genes & Development 10, 1382-94. Lange, T.D., 2005. Shelterin: the protein complex that shapes and safeguards human telomeres. Genes & Development 19, 2100-10. Lei, M., Guo, H., Qiu, W., Lai, X., Yang, T., Widelitz, R.B., Chuong, C.M., Lian, X. and Yang, L., 2014. Modulating hair follicle size with Wnt10b/DKK1 during hair regeneration. Exp Dermatol 23, 407-13. Li, B. and Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323. Li, B., Mackay, D.R., Dai, Q., Li, T.W., Nair, M., Fallahi, M., Schonbaum, C.P., Fantes, J., Mahowald, A.P., Waterman, M.L., Fuchs, E. and Dai, X., 2002. The LEF1/beta -catenin complex activates movo1, a mouse homolog of Drosophila ovo required for epidermal appendage differentiation. Proceedings of the National Academy of Sciences of the United States of America 99, 6064-9. Lim, X., Tan, S.H., Yu, K.L., Lim, S.B. and Nusse, R., 2016. Axin2 marks quiescent hair follicle bulge stem cells that are maintained by autocrine Wnt/β-catenin signaling. Proceedings of the National Academy of Sciences of the United States of America 113, 201601599. Lin, C.M., Jiang, T.X., Widelitz, R.B. and Chuong, C.M., 2006. Molecular signaling in feather morphogenesis. Curr Opin Cell Biol 18, 730-41. Lin, J. and Yue, Z., 2018. Coupling of apical-basal polarity and planar cell polarity to interpret the Wnt signaling gradient in feather development. Development 145, dev162792. Liu, C., Sello, C.T., Sun, Y.F., Zhou, Y.X., Lu, H.T., Sui, Y.J., Hu, J.T., Xu, C.G., Sun, Y., Liu, J., Li, S.Y., Zhang, Y.M. and Zhang, K.Y., 2018. De Novo Transcriptome Sequencing Analysis of Goose (Anser anser) Embryonic Skin and the Identification of Genes Related to Feather Follicle Morphogenesis at Three Stages of Development. International Journal Of Molecular Sciences 19. Liu, C.T., Min, L., Wang, Y.J., Li, P., Wu, Y.D. and Zhang, S.T., 2019. shRNAmediated knockdown of KNTC1 suppresses cell viability and induces apoptosis in esophageal squamous cell carcinoma. Int J Oncol 54, 1053-1060. Liu, X., Driskell, R.R., Luo, M., Abbott, D., Filali, M., Cheng, N., Sigmund, C.D. and Engelhardt, J.F., 2004. Characterization of Lef-1 Promoter Segments that Facilitate Inductive Developmental Expression in Skin. Journal of Investigative Dermatology 123, 264-274. Mao, X., Cai, T., Olyarchuk, J.G. and Wei, L., 2005. Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787-93. Martínez, P., Thanasoula, M., Muñoz, P., Liao, C., Tejera, A., McNees, C., Flores, J.M., Fernández-Capetillo, O., Tarsounas, M. and Blasco, M.A., 2009. Increased telomere fragility and fusions resulting from TRF1 deficiency lead

to degenerative pathologies and increased cancer in mice. Genes & Development 23, 2060. Mendoza-Reinoso, V. and Beverdam, A., 2018. Epidermal YAP activity drives canonical WNT16/β-catenin signaling to promote keratinocyte proliferation in vitro and in the murine skin. Stem Cell Research 29, 15-23. Meyers, C.A., Shen, J., Lu, A. and James, A.W., 2018. WNT16 induces proliferation and osteogenic differentiation of human perivascular stem cells. Journal of Orthopaedics 15, 854-857. Mou, C., Pitel, F., Gourichon, D., Vignoles, F., Tzika, A., Tato, P., Yu, L., Burt, D.W., Bed'hom, B., Tixier-Boichard, M., Painter, K.J. and Headon, D.J., 2011. Cryptic patterning of avian skin confers a developmental facility for loss of neck feathering. PLoS Biol 9, e1001028. Muñoz, P., Blanco, R., Carcer, G.D., Schoeftner, S., Benetti, R., Flores, J.M., Malumbres, M. and Blasco, M.A., 2009. TRF1 Controls Telomere Length and Mitotic Fidelity in Epithelial Homeostasis. Molecular & Cellular Biology 29, 1608-25. Nalesso, G., Thomas, B.L., Sherwood, J.C., Yu, J., Addimanda, O., Eldridge, S.E., Thorup, A.S., Dale, L., Schett, G., Zwerina, J., Eltawil, N., Pitzalis, C. and Dell'Accio, F., 2017. WNT16 antagonises excessive canonical WNT activation and protects cartilage in osteoarthritis. Ann Rheum Dis 76, 218-226. Nie, Y., Li, S., Zheng, X., Chen, W., Li, X., Liu, Z., Hu, Y., Qiao, H., Qi, Q., Pei, Q., Cai, D., Yu, M. and Mou, C., 2018. Transcriptome Reveals Long Non-coding RNAs and mRNAs Involved in Primary Wool Follicle Induction in Carpet Sheep Fetal Skin. Front Physiol 9, 446. Noramly, S., Freeman, A. and Morgan, B.A., 1999. beta-catenin signaling can initiate feather bud development. Development 126, 3509-21. Noramly, S. and Morgan, B., 1998. BMPs mediate lateral inhibition at successive stages in feather tract development. Development 125, 3775-3787. Novel, G., 1973. Feather pattern stability and reorganization in cultured skin. Journal of Embryology & Experimental Morphology 30, 605. Okamoto, K., Iwano, T., Tachibana, M. and Shinkai, Y., 2008. Distinct Roles of TRF1 in the Regulation of Telomere Structure and Lengthening. Journal of Biological Chemistry 283, 23981-8. Painter, K.J., Ho, W. and Headon, D.J., 2017. A chemotaxis model of feather primordia pattern formation during avian development. Journal of Theoretical Biology 437, 225-238. Pazzaglia, I., Mercati, F., Antonini, M., Capomaccio, S., Cappelli, K., Dall’Aglio, C., La Terza, A., Mozzicafreddo, M., Nocelli, C. and Pallotti, S., 2019. PDGFA in Cashmere Goat: A Motivation for the Hair Follicle Stem Cells to Activate. Animals 9, 38. Prum, R.O., 1999. Development and evolutionary origin of feathers. J Exp Zool 285, 291-306. Prum, R.O. and Brush, A.H., 2002. The evolutionary origin and diversification of feathers. Quarterly Review of Biology 77, 261-295. Qiu, W., Lei, M., Zhou, L., Bai, X., Lai, X., Yu, Y., Yang, T. and Lian, X., 2017. Hair follicle stem cell proliferation, Akt and Wnt signaling activation in TPA-induced hair regeneration. Histochem Cell Biol 147, 749-758. Reddy, S., Andl, T., Bagasra, A., Min, M.L., Epstein, D.J., Morrisey, E.E. and Millar, S.E., 2001. 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. Mechanisms of Development 107, 69-82. Reddy, S.T., Thomas, A., Min-Min, L., Morrisey, E.E. and Millar, S.E., 2004. Expression of Frizzled genes in developing and postnatal hair follicles. Journal of Investigative Dermatology 123, 275-282. Rivera-Gonzalez, G.C., Shook, B.A., Andrae, J., Holtrup, B., Bollag, K., Betsholtz, C., Rodeheffer, M.S. and Horsley, V., 2016. Skin Adipocyte Stem Cell Self-Renewal Is Regulated by a PDGFA/AKT-Signaling Axis. Cell Stem Cell 19, 738-751. Robinson, M.D., McCarthy, D.J. and Smyth, G.K., 2010. edgeR: a Bioconductor package for differential expression

analysis of digital gene expression data. Bioinformatics 26, 139-40. Ruxton, G.D., Persons, W.S. and Currie, P.J., 2017. A continued role for signalling functions in the early evolution of feathers. Evolution; international journal of organic evolution 71. Seifert, J.R.K. and Marek, M., 2007. Frizzled/PCP signalling: a conserved mechanism regulating cell polarity and directed motility. Nature Reviews Genetics 8, 126-138. Sello, C.T., Liu, C., Sun, Y., Msuthwana, P., Hu, J., Sui, Y., Chen, S., Zhou, Y., Lu, H., Xu, C., Sun, Y., Liu, J., Li, S. and Yang, W., 2019. De Novo Assembly and Comparative Transcriptome Profiling of Anser anser and Anser cygnoides Geese Species' Embryonic Skin Feather Follicles. Genes (Basel) 10, undefined. Sennett, R., Wang, Z., Rezza, A., Grisanti, L., Roitershtein, N., Sicchio, C., Mok, K.W., Heitman, N.J., Clavel, C., Ma'ayan, A. and Rendl, M., 2015. An Integrated Transcriptome Atlas of Embryonic Hair Follicle Progenitors, Their Niche, and the Developing Skin. Dev Cell 34, 577-91. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B. and Ideker, T., 2003. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 13, 2498-504. Shen, J., Chen, X., Jia, H., Meyers, C.A., Shrestha, S., Asatrian, G., Ding, C., Tsuei, R., Zhang, X. and Peault, B., 2017. Effects of WNT3A and WNT16 on the Osteogenic and Adipogenic Differentiation of Perivascular Stem/Stromal Cells. Tissue Engineering Part A 24. Simonchazottes, D., Tutois, S., Kuehn, M., Evans, M., Bourgade, F., Cook, S., Davisson, M.T. and Guénet, J.L., 2006. Mutations in the gene encoding the low-density lipoprotein receptor LRP4 cause abnormal limb development in the mouse. Genomics 87, 673-677. Song, H., Wang, Y. and Goetinck, P.F., 1996. Fibroblast growth factor 2 can replace ectodermal signaling for feather development. Proceedings of the National Academy of Sciences of the United States of America 93, 10246-10249. Song, H.K., Lee, S.H. and Goetinck, P.F., 2010. FGF-2 signaling is sufficient to induce dermal condensations during feather development. Developmental Dynamics 231, 741-749. St-Jacques, B., Dassule, H.R., Karavanova, I., Botchkarev, V.A., Li, J., Danielian, P.S., McMahon, J.A., Lewis, P.M., Paus, R. and McMahon, A.P., 1998. Sonic hedgehog signaling is essential for hair development. Curr Biol 8, 1058-68. Storey, J.D. and Tibshirani, R., 2003. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A 100, 9440-5. Strong, R.M., 1916. A Study of the Structure of Feathers, with Reference to Their Taxonomic Significance by Asa C. Chandler. Dental Materials 29, e27. Tao, H., Yoshimoto, Y., Yoshioka, H., Nohno, T., Noji, S. and Ohuchi, H., 2002. FGF10 is a mesenchymally derived stimulator for epidermal development in the chick embryonic skin. Mech Dev 116, 39-49. Th, P.W. and Currie, P.J., 2015. Bristles before down: a new perspective on the functional origin of feathers. Evolution 69, 857-862. Ting-Berreth, S.A. and Chuong, C.M., 1996. Local Delivery of TGF β2 Can Substitute for Placode Epithelium to Induce Mesenchymal Condensation during Skin Appendage Morphogenesis. Developmental Biology 179, 347. Ting-Berreth, S.A. and Chuong, C.M., 2010. Sonic Hedgehog in feather morphogenesis: induction of mesenchymal condensation and association with cell death. Dev Dyn 207, 157-170. Tomoko, T., Itaru, I., Ichiro, T., Kiyomi, H., Hiroko, T., Juichi, I. and Ryoichiro, K., 2013. Hedgehog signaling regulates prosensory cell properties during the basal-to-apical wave of hair cell differentiation in the mammalian cochlea. Development 140, 3848-3857. Tripurani, S.K., Yan, W., Ying-Xin, F., Massod, R., Lily, W., Min-Hyung, L., Starost, M.F., Rubin, J.S., Johnson, G.R. and Ginsberg, M.H. Suppression of Wnt/β-catenin signaling by EGF receptor is required for hair follicle

development. Molecular Biology of the Cell. van Genderen, C., Okamura, R.M., Farinas, I., Quo, R.G., Parslow, T.G., Bruhn, L. and Grosschedl, R., 1994. Development of several organs that require inductive epithelial-mesenchymal interactions is impaired in LEF-1-deficient mice. Genes Dev 8, 2691-703. Wada, N., Kawakami, Y., Ladher, R., Francis-West, P.H. and Nohno, T., 1999. Involvement of Frzb-1 in mesenchymal condensation and cartilage differentiation in the chick limb bud. International Journal of Developmental Biology 43, 495-500. Wang, H.D., Yang, L., Yu, X.J., He, J.P., Fan, L.H., Dong, Y.J., Dong, C.S. and Liu, T.F., 2012. Immunolocalization of β-catenin and Lef-1 during postnatal hair follicle development in mice. Wang, L., Feng, Z., Wang, X., Wang, X. and Zhang, X., 2010. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136-8. Wang, Y., Guo, N. and Nathans, J., 2006. The role of Frizzled3 and Frizzled6 in neural tube closure and in the planar polarity of inner-ear sensory hair cells. J Neurosci 26, 2147-56. Wang, Z., Jinnin, M., Kobayashi, Y., Kudo, H., Inoue, K., Nakayama, W., Honda, N., Makino, K., Kajihara, I., Makino, T., Fukushima, S., Inagaki, Y. and Ihn, H., 2015. Mice overexpressing integrin alphav in fibroblasts exhibit dermal thinning of the skin. J Dermatol Sci 79, 268-78. Weatherbee, S.D., Anderson, K.V. and Niswander, L.A., 2006. LDL-receptor-related protein 4 is crucial for formation of the neuromuscular junction. Development 133, 4993-5000. Wessells, N.K., 1962. Tissue interactions during skin histodifferentiation. Developmental biology 4, 87-107. Widelitz, R.B. and Chuong, C.M., 1999. Early Events in Skin Appendage Formation: Induction of Epithelial Placodes and Condensation of Dermal Mesenchyme. Journal of Investigative Dermatology Symposium Proceedings 4, 302-306. Widelitz, R.B., Jiang, T.X., Lu, J. and Chuong, C.M., 2000. beta-catenin in epithelial morphogenesis: conversion of part of avian foot scales into feather buds with a mutated beta-catenin. Dev Biol 219, 98-114. Williams, S.E., Garcia, I., Crowther, A.J., Li, S., Stewart, A., Liu, H., Lough, K.J., O'Neill, S., Veleta, K., Oyarzabal, E.A., Merrill, J.R., Shih, Y.Y. and Gershon, T.R., 2015. Aspm sustains postnatal cerebellar neurogenesis and medulloblastoma growth in mice. Development 142, 3921-32. Woo, W.M., Zhen, H.H. and Oro, A.E., 2012. Shh maintains dermal papilla identity and hair morphogenesis via a Noggin-Shh regulatory loop. Genes Dev 26, 1235-1246. Wu, P., Yan, J., Lai, Y.C., Ng, C.S., Li, A., Jiang, X., Elsey, R.M., Widelitz, R., Bajpai, R. and Li, W.H., 2017. Multiple Regulatory Modules Are Required for Scale-to-Feather Conversion. Molecular Biology & Evolution 35, 417-430. Wu, W., Xu, R.F., Li, C.H. and Wu, C.X., 2008a. Characterization of Embryonic Feather Follicle Development in the Chinese Indigenous Jilin White Goose. Asian Australasian Journal of Animal Sciences 21, 346-352. Wu, W., Xu, R.F., Xiao, L., Xu, H. and Gao, G., 2008b. Expression of the beta-catenin gene in the skin of embryonic geese during feather bud development. Poultry Science 87, 204. Xu, R.F., Wu, W. and Xu, H., 2007. Investigation of feather follicle development in embryonic geese. Poult Sci 86, 2000-7. Yamaguchi, Y., Morita, A., Maeda, A. and Hearing, V.J., 2009. Regulation of skin pigmentation and thickness by Dickkopf 1 (DKK1). J Investig Dermatol Symp Proc 14, 73-5. Yamaguchi, Y., Passeron, T., Hoashi, T., Watabe, H., Rouzaud, F., Yasumoto, K., Hara, T., Tohyama, C., Katayama, I., Miki, T. and Hearing, V.J., 2008. Dickkopf 1 (DKK1) regulates skin pigmentation and thickness by affecting Wnt/beta-catenin signaling in keratinocytes. Faseb j 22, 1009-20. Yang, S., Shi, Z., Xiaoqian, O.U. and Liu, G., 2019. Whole-genome resequencing reveals genetic indels of feathered-leg traits in domestic chickens. Journal of Genetics. Young, M.D., Wakefield, M.J., Smyth, G.K. and Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for

selection bias. Genome Biol 11, R14. Yu, H., Smallwood, P.M., Wang, Y., Vidaltamayo, R., Reed, R. and Nathans, J., 2010. Frizzled 1 and frizzled 2 genes function in palate, ventricular septum and neural tube closure: general implications for tissue fusion processes. Development 137, 3707-17. Yu, M., Yue, Z., Wu, P., Wu, D.Y., Mayer, J.A., Medina, M., Widelitz, R.B., Jiang, T.X. and Chuong, C.M., 2004. The biology of feather follicles. Int J Dev Biol 48, 181-91. Zhang, H., Nan, W., Wang, S., Si, H. and Li, G., 2018. Balance between fibroblast growth factor 10 and secreted frizzled-relate protein-1 controls the development of hair follicle by competitively regulating β-catenin signaling. Biomedicine & Pharmacotherapy 103, 1531-1537. Zhao, B., Chen, Y., Yan, X., Hao, Y., Zhu, J., Weng, Q. and Wu, X., 2017. Gene Expression Profiling Analysis Reveals Fur Development in Rex Rabbits (Oryctolagus cuniculus). Genome 60, 1. Zhou, X., Wang, R., Fan, L., Li, Y., Ma, L., Yang, Z., Yu, W., Jing, N. and Zhu, X., 2005. Mitosin/CENP-F as a negative regulator of activating transcription factor-4. J Biol Chem 280, 13973-7.