Phenotype, Histological, and Bone Structure Observations.
Morphological, histological, and anatomical structures were observed to determine the histological and bone structures of goose knobs. Morphological analysis showed that the knob accentuated the facial contours of the Yangzhou geese (Fig. 1A, 1B). The average knob length, width, and height were approximately 31.5 mm, 36.5 mm, 32.4 mm at 380 days of age, respectively, and these parameters varied significantly (P < 0.05). The maximum length, width, and height were approximately 36.4 mm, 40.3 mm, and 37.2 mm, and the corresponding minimum values were close to 26.6 mm, 32.8 mm, and 27.6 mm, respectively (Fig. 1D). Specific measurements standards are shown in Fig. 1C, and these standards were applied when measuring the bony crest.
Anatomically, the knob was found to include the skin derivative and bony crest. Although the knob was located at the base of the beak, it was joined to a typically frontal crest, rather than on the upper beak (Fig. 2A–2C). Unlike the whole knob, the bony crest varied greatly in terms of the length and height (P < 0.05), but not the width (Fig. 2D). Furthermore, increased crest sizes did not correlate with the hardness of the peripheral bones. Instead, the outer layer of the bone became thinner as the hump became larger, in a manner reminiscent of blowing a balloon. In addition, bony pores were often observed in large bony crests.
Based on the observations and measurements above, we analyzed the correlation between the knob size and bony crest size. The knob length, width, and height each significantly correlated with the same parameters in the bony crest (Table 1).
Histological measurements were performed to determine the thickness of five layers in the skin, including the stratum corneum (Sc), stratum spinosum (Ss), stratum epidermis (El), stratum reticular (Sr), and stratum corium (Sc). The stratum corneum is composed of flat, dead cells (corneocytes) embedded in a matrix of lipids3, which were observed as highly eosinophilic cells with few nuclei following HE staining. Denser structures and significant enrichment for eosinophilic substances were observed in skin tissues from large knobs versus those from small knobs. Greater thicknesses of stratum corneum, stratum spinosum, and stratum reticular were observed in geese with large knobs than in those with small knobs (P < 0.05) (Fig. 3A–3C).
RNA Library Construction and Sequencing
In this study, an average of 50,136,474 raw reads was obtained from the L and S samples, and the average number of clean reads was 49,500,375. All downstream analyses were based on high-quality, clean sequence data. The error rates were all less than 0.025%. Approximately 92.01%–96.37% of the clean reads in the libraries were mapped to the Anser cygnoides reference genome and the percentage of phred quality scores of >20 (Q20) was >98% in all samples. The sequence read statistics are summarized in Table S3.
DEG Screening in Geese With Different Knob Sizes
DEGs were screened in large and small goose knobs. Scatter and volcano plots showed variations in mRNA-expression levels between geese with small (S group) and large (L group) knobs (Fig. 4A, 4B). We identified 415 differentially expressed mRNAs (Fig. 4C and Table S4; FC ≥ 2), including 357 upregulated genes and 58 downregulated genes. Most upregulated genes were related to hormone secretion, skeletal system development, and tissue composition, including adrenoceptor alpha 1A (ADRA1A); AGT; bromodomain, testis-specific (BRDT); prostaglandin F2 receptor inhibitor (PTGFRN); relaxin/insulin-like family peptide receptor 3 (RXFP3); SRD5α2; TSHR; IGF1; insulin-like growth factor 2 (IGF2), insulin-like growth factor binding protein 4 (IGFBP4); insulin-like growth factor binding protein 5 (IGFBP5); bone-morphogenetic protein/retinoic acid inducible neural-specific 3 (BRINP3); calcium channel, voltage-dependent subunit alpha1 H (CACNA1H); calneuron 1 (CALN1); collagen-related genes (COL12A1, COL15A1, COL1A2, COL20A1, COL21A1, COL26A1, COL3A1, COL5A2, COL6A1, and COL6A2); DIO3; fibulin 1 (FBLN1); fibulin 2 (FBLN2); fibulin 7 (FBLN7); FBN1; fibrillin 2 (FBN2); fibroblast growth factor 7 (FGF7); fibroblast growth factor 9 (FGF9); NPPC; natriuretic peptide receptor 1 (NPR1); and OGN. The downregulated DEGs accounted for a small proportion of the total. These genes are known to be mainly related to enzymes, such as aldehyde dehydrogenase 3 family (ALDH3B1); mitogen-activated protein kinase 15 (MAP3K15); monooxygenase, DBH-like 1 (MOXD1); and protein kinase, X-linked (PRKX).
GO Enrichment and KEGG Pathway Analyses for DEGs
To further elucidate the functional roles of DEGs in determining the knob size, we performed GO and KEGG pathway-enrichment analysis for the DEGs using the online programs, Goatools and KOBAS. The DEGs were categorized into three main GO categories, namely biological process, cellular component, and molecular function. We identified 455 significantly enriched GO terms (Table S4) in 94 molecular functions, 52 cellular components, and 309 biological processes (P-value ≤ 0.05). The cellular component category mainly relates to extracellular components and chemical composition. Among the hormone-related GO terms, 13 were associated with insulin, and four were thyroid hormone-related. ACTA2, ACTA1, and ACTG2 are enriched in pathways related to tissue migration and mesenchymal migration. Some DEGs were associated with calcium ion binding, such as SULF1, SPARC, FBN1, and FBN2.
In addition, the 10 genes with the greatest significant differences (|log2 FC|) for each GO term were plotted (Fig. 5). GO-term analysis showed that significantly upregulated DEGs in tissues from large knobs were highly related to certain categories, including extra-cellular structure organization, extra-cellular matrix organization, tissue migration, collagen fibril organization, mesenchyme migration, and skeletal system development. Furthermore, we mapped the DEGs in the KEGG pathway database and classified all pathways into five categories.
We mapped the 415 DEGs to 210 KEGG pathways, and the 20 most enriched pathways are shown in the bubble chart in Fig. 6 The associated KEGG pathway terms include protein digestion and absorption, PI3K-Akt signaling pathway, and ECM-receptor-interaction. We also identified many signaling hormone-related pathway terms, such as relaxin signaling pathway; thyroid hormone synthesis (including CREB3L1, GPX7, ADCY3, and TSHR); endocrine resistance (including CDKN2C,IGF1, and ADCY3); cortisol synthesis and secretion (including CREB3L1, ADCY3, and CACNA1H); insulin secretion (including ADCY3 and CREB3L1); pancreatic secretion (including ATP2B2 and ADCY3); insulin resistance (including GFPT2 and CREB3L1); steroid hormone biosynthesis (including SRD5α2); dopaminergic synapse (including GRIA3 and CREB3L1); parathyroid hormone synthesis, secretion and action (including CREB3L1 and ADCY3); renin secretion (including NPR1 and AGT); and TGF-β signaling pathway (including DCN, FMOD, LTBP1, TGFB2, and BMP5).
Protein–Protein Interaction (PPI) Analysis
A PPI network was generated using the STRING program in Cytoscape 4 software, based on the DEGs identified in this study (Fig. 7). FBN1 and DCN were associated with the most genes; thus, they were in the center of the PPI network. The SPARC, IGF1, LUM, and COL1A2 genes were the four genes with high degree.
Validation of DEGs in Geese With Different Knob Sizes
To validate the RNA-seq results, 13 relatively important DEGs (ADCY3, AGT, BMP5, DCN, DIO3, FBN1, IGF1, LUM, NPPC, OGN, SPARC, SRD5α2, and TSHR1) were selected for RT-qPCR analysis. All of these DEGs showed concordant expression patterns between the RNA-seq and RT-qPCR results (Fig. 8).