This experiment was approved by the Animal Health Care Committee of Animal Science and Technology College of Jilin Agricultural University (Approval No. GR (J) 18–003). Breeder eggs of Zi goose (Anser cygnoides) were obtained from the geese breeding research center, Jilin Agricultural University (Changchun, China), and incubated in an incubator in laboratory according to the routine procedure.

Candling procedure was conducted on day 6 of incubation to identify and remove the unfertilized eggs and the dead embryos. In the pre-experiment, 60 fertilized eggs were used to observe the feather follicle morphogenesis (from E9 to E25). Then, 400 fertilized eggs were divided into five groups of 80 eggs each in ovo injection experiments. The different doses of CHIR-99021 solutions were prepared and preheated to 37.9°C before injecting (E9). Then, the CHIR-99021 solutions were injected into the egg yolk sac. The group names were abbreviated according to the dose they received, consisting of 1) non-injection group (BLANK group); 2) injection of 100 μl of 0.9% normal saline (NS)/egg (CK group); 3) injection of 100 μl of 0.9% NS/egg containing 1,000 ng CHIR-99021/egg (1,000 ng group); 4) injection of 100 μl of 0.9% NS/egg containing 2,000 ng CHIR-99021/egg (2000 ng group); 5) injection of 100 μl of 0.9% NS/egg containing 5,000 ng CHIR-99021/egg (5,000 ng group).

CHIR-99021 (CHIR-99021 monohydrochloride, MedChemExpress, Shanghai, China) was dissolved in 0.9% NS. The center of the large end of each egg was located under the egg candler and an alcohol swab was used to disinfect the target injection site, then a small hole was drilled at the injection site. As shown in Figure 1A, the volume of 100 μl of the CHIR-99021 solution was rapidly injected into the egg yolk sac using a 0.7 × 30 mm needle. After paraffin sealing, the eggs were returned to the incubator immediately. The non-injection group eggs also stayed outside the incubator environment for 15 min for standardization.

The dorsal skin tissue samples (about 1.5 cm², the cross region of the dorsal midline) were obtained during 0–1 h at E13, E17, E21, and E25. Survivability was the percentage of surviving embryos in each group of the total number of fertilized eggs. A part of the skin tissue was frozen in liquid nitrogen and stored at −80°C for western blotting and RNA sequencing (RNA-seq) analyses. The samples were mixed from three embryos to ensure sufficient material for molecular biology assays. Other part of the skin tissue was fixed in 4% paraformaldehyde for further histological observations.

The feathers were obtained at E17, E21, and E25. The sampling region was the cross region of the dorsal midline. The collected feathers were cleaned with PBS. A part of the feathers was used for feather length measurement using a vernier caliper (precision 0.01 mm). Other part of feathers was stuck on slides for feather width measurement. The width of feather center was measured using CU-6 fiber fineness meter (Olympus, Tokyo, Japan).

Feather follicle morphological observations were conducted using hematoxylin and eosin (H&E) staining. The fixed tissues were dehydrated in different concentrations of ethanol (50,75, 85, 90, 95, and 100%) and were soaked twice in xylene for clearing. The tissues were embedded in paraffin according to conventional methods.

Paraffin-embedded 5 μm sections of dorsal skin tissues from E9 to E25 were prepared. The sections were stained using H&E according to standard methods and were photographed using a microscope (Olympus, Tokyo, Japan).

The distribution of feather follicles was observed and photographed from the longitudinal skin sections to measure feather follicle diameter using microscope software (cellSens v.3.1, Olympus, Tokyo, Japan). The distance of the outermost layer is the maximum in feather follicle and was measured to assess the feather follicle diameter (red arrows in Figure 3A). A total of 10 regions of each transverse skin section were photographed to calculate the number of feather follicles, which was used to calculate the feather follicle density.

The total protein samples for electrophoresis were extracted using RIPA buffer (Beyotime Biotechnology, Shanghai, China), and the total protein concentration was determined using a BCA protein concentration detection kit (Beyotime Biotechnology). For each sample, the volume containing 10 μg total proteins was pipetted onto a 10% SDS-PAGE. The proteins were transferred to the surface of the PVDF membrane. The membranes were blocked with 5% skim milk and then incubated with primary antibodies at 4°C overnight. The antibodies identifiers and dilutions were as follows: rabbit anti-β-catenin antibody (Cat# bs-23663R, bioss, Beijing, China), rabbit anti-GSK-3 β antibody (Cat# bs-0023M, bioss), 1:2,000,1:2,000; rabbit anti-c-Myc antibody (Cat# bs-4963R, bioss), 1:2,000; rabbit anti-TCF4 antibody (Cat# bs-1280R, bioss), 1:2000; rabbit anti-LEF1 antibody (Cat# bs-1843R, bioss), 1:1,500; rabbit anti-FZD4 antibody (Cat# bs-13217R, bioss), 1:2,000; ACTB rabbit monoclonal antibody (Cat# AC026, Abclonal, Wuhan, China), 1:50,000. At the end, the membranes were incubated with a goat anti-rabbit IgG (H + L)/HRP antibody (Biosynthesis Biotechnology, Beijing, China, 1:50,000 dilution) for 1 h at RT. Goose is not a common experimental animal, which causes a limitation of commercial antibodies. The antibodies reacting with chicken were used in the study. A prior experiment verified that these antibodies can specifically react with goose. The molecular weights between some proteins are nearby, stripping, and thereby was used to avoid errors cross gels. A stripping buffer (Solarbio, Beijing, China) was used for 12 min to strip the antibodies bound on the membranes, then blocking, and antibody incubating were conducted as described above. Finally, the membranes were visualized with an ECL Test Kit (Millipore, Darmstadt, Germany) under a Bio-Rad imaging system (Bio-Rad, Hercules, CA, United States). The chemiluminescence of each protein blotting was quantified using the ImageJ software (NIH, Bethesda, MD, United States). The protein levels were normalized by β-actin as an internal control.

The slides with 5-μm thick paraffin-embedded dorsal skin tissues sections were baked at 65°C for 3 h, dewaxed twice in xylene, and hydrated in descending grades of alcohol. Then, the slides succumbed to antigen retrieval by boiling in sodium citrate buffer (pH 6.0) for 20 min after deparaffinization and hydration and were rinsed with PBS for 5 min of 3 times at the end of each step. Following the steps of deparaffinization, rehydration, and antigen retrieval, the sections were blocked with 3% hydrogen peroxide at RT for 20 min and with 3% normal goat serum at RT for 20 min, and were incubated overnight at 4°C with 50 μl of rabbit anti-β-catenin antibody (Cat# bs-23663R, bioss, Beijing, China) at 1:100 diluent. Then, the sections were incubated, respectively, with a secondary antibody linked with biotin and the HRP-linked streptomycin according to UltraSensitive™ SP Rabbit IHC Kit (MXB, Fuzhou, Fujian, China) instructions. DAB (MXB, Fuzhou, Fujian, China) and hematoxylin were, respectively, used for chromogen/substrate and nuclear staining. Finally, the slides were observed and photographed using a microscope system (Olympus, Tokyo, Japan) after coverslipping.

The dorsal skin tissue samples from CK and the 5,000 ng group at E13 and E25 were used for RNA-seq. Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol. The harvested RNA was quantified using a nanodrop spectrophotometer (BIO-DL, Shanghai, China). The RNA quality was assessed on an Agilent 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and agarose gel electrophoresis. The RNA integrity numbers (≥7.3) and concentrations (≥211 ng/μl) met the requirements of library construction (Supplementary Table S1). After total RNA extraction, eukaryotic mRNA was enriched by oligo (dT) beads, while prokaryotic mRNA was enriched by removing rRNA using Ribo-Zero™ Magnetic Kit (Epicentre Biotechnologies, Madison, WI, United States). Then, the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP, and buffer. Then, the cDNA fragments were purified with a QiaQuick PCR extraction kit (QIAGEN, Dusseldorf, Germany), and were end-repaired, poly(A) added and ligated to Illumina sequencing adapters. The size of ligation products was selected using agarose gel electrophoresis. The purified amplicons were sequenced using Illumina HiSeq™ 4,000 platform by Gene Denovo Biotechnology Co. (Guangzhou, China).

The qualities of the raw reads were checked using fastp software (version 0.18.0) (Chen, et al., 2018). The reads containing adapters, reads containing more than 10% unknown nucleotides (N), and low-quality reads containing more than 50% of low-quality (Q-value ≤ 10) bases were removed. Then, a total of 538502706 raw reads and 536003454 clean reads were obtained (clean reads rates ≥99.49%). The bases Q20 (≥98.28%), Q30 (≥95.13%), and GC content (53.42–56.98%) in clean reads were within reasonable ranges (Supplementary Table S2). Next, the clean reads were mapped to RNA to identify and remove residual rRNA reads. Then, the rRNA removed high-quality clean reads and was further used for assembly and gene abundance calculations. The high-quality clean reads were mapped to the reference transcriptome using Bowtie2 software by default parameters (version 2.2.8) (Langmead and Salzberg, 2012). The clean reads assembled 91,236 unigenes from 93525236 bases by using Trinity rna-seq software (version 2.13.2) (Grabherr, et al., 2011). To get the protein function annotation information of the unigenes with the highest sequence similarity (Ashburner, et al., 2000; Conesa, et al., 2005), the unigenes were annotated to Nr, SwissProt, KEGG, and COK/KOG databases by using Blastx (NCBI, Bethesda, MD, United States). The number of homologous sequences between species, and these unigenes’ function annotations are shown in Supplementary Figure S1.

The gene abundance was calculated and normalized by RPKM (reads per kb per million reads) using RESM software (version: 1.3.3) (Mortazavi, et al., 2008). To identify differentially expressed genes between groups, the edgeR package was used. We identified genes with a fold change ≥2 and a p < 0.05 in comparison as differentially expressed genes (DEGs). The DEGs were then subjected to enrichment analysis of GO terms and KEGG pathways (Kanehisa, et al., 2008). The GO terms and KEGG pathways with p < 0.05 were defined as significantly enriched terms and pathways.

The differently expressed transcription factors were classified and calculated according to their families. The shared and unique differently expressed transcription factors between E13 and E25 were ranked by -log₁₀ (p-value) and presented as a heatmap. The PPI network of the differently expressed transcription factors was constructed using STRING (https://www.string-db.org) with confidence >0.9.

The number of 16 DEGs screened by de novo transcriptome was randomly selected to conduct q-PCR validation. The primers are listed in Supplementary Table S3. The cDNA was synthesized for q-PCR using a PrimeScript RT Reagent Kit with gDNA Eraser (Morey Biosciences, Shanghai, China). q-PCR was performed with SYBR Green Reagent (Morey Biosciences, Shanghai, China) in a CFX Connect Real-Time PCR Detection System (Bio-Rad, Hercules, CA, United States). The total volume of q-PCR reactions was 20 μl comprising 2 μl of cDNA, 0.6 μl of both forward, and reverse primers, 6.8 μl of RNase-free H₂O, and 10 μl SYBR Green Reagent. The amplification conditions were as follows: pre-denaturation at 95°C for 5 min, then 40 cycles of amplification (95°C for 15 s and 59°C for 60 s). Expression was quantified using the 2^–ΔΔCt method. The β-actin gene was chosen as the control gene to normalize the mRNA expression levels. The q-PCR results supported the reliability of de novo transcriptome in this study (Supplementary Figure S2).

The statistical tests were performed with SPSS 23.0 software (IBM, Armonk, NY, United States). The data were visualized using the GraphPad Prism 8 software (GraphPad, San Diego, CA, United States). The statistical significance was determined using one-way ANOVA by Duncan’s multiple range test. The significant difference in the data was considered as p < 0.05. The results were expressed as mean ± SEM in the figures.

To determine the time nodes for in ovo injection and sampling, feather follicle morphology was observed during E9 to E25. As shown in Figure 1A, the dermal cell condensation appeared from E9 and that was considered the outset of feather follicle development. The short bud structure started to format at E13. By E17, the stage from the cavity to follicle started and the long bud structure started sinking into the epidermis. Feather follicles have sunk into the epidermis and developed basic structures by E21. At E25, the feather follicles were completely developed, and at this period, the muscles, nerves, glands, and blood vessels connected with mature feather follicles gradually become abundant in the dermis. Therefore, E9 was chosen as the time node for in ovo injection; E13, E17, E21, and E25 were chosen as the time nodes for sampling. The survival rates in the embryos that received CHIR-99021 injection were above 95% (Figure 1B), which suggested that CHIR-99021 and the doses were in a safe range for embryos in this study.

The feather length results (Figure 2A) and feather width results (Figures 2B,C) suggested that in ovo injection of CHIR-99021 had the ability to promote feather growth. Compared with the BLANK and CK groups, in ovo injection of 5,000 ng CHIR-99021 significantly (p < 0.05) promoted the feather length at E17, E21, and E25, while in ovo injection of 2,000 ng CHIR-99021 also had significantly (p < 0.05) promoted effect on feather length at E17. Moreover, feather width in the 5,000 ng CHIR-99021 group at E21 and E25 and in the 2,000 ng CHIR-99021 group at E25 were significantly (p < 0.05) increased.

The feather follicle diameter (Figure 3) and feather follicle density (Figure 4) results suggested that CHIR-99021 also promoted feather follicle development. Morphological observations were conducted to investigate feather follicle development during morphogenesis between the groups (Figures 3A, 4A). The results (Figure 3B) showed that feather follicle diameter in the 5,000 ng CHIR-99021 group was significantly (p < 0.05) wider than that in the BLANK and CK groups at E13, E17, E21, and E25. The feather follicle diameter in the 2,000 ng group was significantly (p < 0.05) promoted at E17, E21, and E25, and that in the 1,000 ng group was significantly (p < 0.05) increased at E25. The results of primary feather follicles (Pfs) density (Figure 4B) showed that injection of 5,000 ng CHIR-99021 significantly (p < 0.05) increased Pfs density at E13 and E17, and injection of 2,000 ng CHIR-99021 significantly (p < 0.05) increased Pfs density at E13. However, the Pfs densities in the 1,000 ng group at E21 and in the 5,000 ng group at E25 were significantly (p < 0.05) decreased and the Pfs densities in other injected groups were at a common level (p > 0.05) with the BLANK and CK groups. These were possibly associated with their wider feather follicle diameter leading to a decreased number of feather follicles in the same field of view.

The expression of hub molecules in the Wnt/β-catenin signaling pathway in the dorsal tissues was investigated using western blotting between groups (Figure 5). At E13 (Figure 5A), injection of CHIR-99021 significantly (p < 0.05) inhibited the expression of GSK-3β in the embryonic dorsal tissues. A dose of 5,000 ng had more significantly inhibited the effect compared with the doses of 1,000 ng and 2,000 ng. In addition, CHIR-99021 had the ability to upregulate the expression of FZD4, β-catenin, the transcription factors LEF1 and TCF4, and the target c-Myc in the Wnt signaling pathway. Moreover, similar regulations of these Wnt molecules were found at E17 and E21 (Figures 5B,C). However, by E25 (Figure 5D), the expression of GSK-3β in the 1,000 ng and 2,000 ng groups was at a common level (p > 0.05) compared with the BLANK and CK groups and the level of GSK-3β in the 5,000 ng group were significantly (p < 0.05) increased. This might associate with complete feather follicle morphogenesis and injection timing. The levels of FZD4, β-catenin, LEF1, TCF4, and c-Myc maintained an upregulated trend in the groups that received CHIR-99021 injection. These results suggested that, at the level of protein, in ovo CHIR-99021 injection promoted feather growth and feather follicle development by inhibiting GSK-3β and upregulating FZD4, β-catenin, LEF1, TCF4, and c-Myc. The raw western blotting images are attached in Supplementary Material S1.

At the same time, IHC-P was used to investigate the expression and distribution of β-catenin which is considered the core molecule in the Wnt signaling pathway (Figure 6). The positive expression area of β-catenin in treated embryos was visually higher compared to the BLANK and CK groups (Figure 6A). β-Catenin was mainly expressed in the epidermis and placode at E13, in the inner root sheath and outer root sheath at E17, and in the epidermal collar and epidermis at E21 and E25 (Figure 6B). The IHC-P results supported that the promoting effects of CHIR-99021 on feather growth and feather follicle development were associated with the regulation of the Wnt signaling pathway by CHIR-99021.

RNA-seq was conducted to comprehensively assess the effects of CHIR-99021 injection on transcription. A total of 3,806 and 2,637 DEGs (fold change ≥2 and p < 0.05) were screened between the CK and 5,000 ng groups, respectively, at E13 and E25. There were 2,172 upregulated DEGs and 1,634 downregulated DEGs at E13 and 1,540 upregulated DEGs and 1,097 downregulated DEGs at E25 (Figure 7A).

Next, GO enrichment (Figure 7B) and KEGG (Figure 7C) enrichment analyses were performed to predict the biological functions conducted by these DEGs. The top 10 enriched GO terms in the BP category at E13 were as follows: translation, peptide biosynthetic process, peptide metabolic process, amide biosynthetic process, organonitrogen compound biosynthetic process, cellular amide metabolic process, protein targeting to ER, SRP-dependent co-translational protein targeting to membrane, co-translational protein targeting to membrane, and establishment of protein localization to the endoplasmic reticulum. In addition, the top 10 enriched GO terms in the BP category at E25 were annotated in the organonitrogen compound biosynthetic process, translation, peptide biosynthetic process, cytoplasmic translation, amide biosynthetic process, SRP-dependent co-translational protein targeting the membrane, co-translational protein targeting the membrane, peptide metabolic process, protein targeting the ER, nuclear-transcribed mRNA catabolic process, and nonsense-mediated decay.

At E13, the top 10 enriched KEGG pathways were as follows: ribosome, glycolysis/gluconeogenesis, antigen processing and presentation, carbon metabolism, biosynthesis of amino acids, fructose and mannose metabolism, Staphylococcus aureus infection, complement and coagulation cascades, systemic lupus erythematosus, and legionellosis. The top 10 pathways at E25 with the highest representation of the DEGs were annotated in the ribosome, glycolysis/gluconeogenesis, biosynthesis of amino acids, carbon metabolism, HIF-1 signaling pathway, fructose and mannose metabolism, pentose phosphate pathway, metabolic pathways, starch and sucrose metabolism, and galactose metabolism. The detailed results of DEGs, GO enrichment, and KEGG enrichment are, respectively, listed in Supplementary Materials S2–S4.

The enrichment analyses illustrated that the DEGs were associated with transcription and translation. Therefore, further analyses were performed to screen transcription factors. All of the transcription factors (TFs) detected in this study were annotated in 1,055 families, and mainly in zf-C2H2, Hompbox, and bHLH families (Figure 8A). The shared and unique differentially expressed TF genes between groups at E13 and E25 were ranked by -log₁₀ (p-value) and shown in Figure 8B. Then, the protein-protein interaction (PPI) network of these differentially expressed TF genes was constructed and predicted that CEBPB, CTNNB1, and SOX2 were the hub nodes (Figure 8C).