The Wanxi White Goose (WWG) is a premium indigenous Chinese goose breed renowned for its high-quality meat and superior down. However, its suboptimal reproductive performance significantly constrains the scaled development of the WWG farming industry (Yang et al., 2024). As an important reproductive organ in female poultry, the ovary directly determines egg production (Chen et al., 2024). The ovarian development and follicular morphology of poultry show significant changes during different egg production cycles. For example, at sexual maturity the ovarian volume increases, the surface gradually forms primordial follicles and the ovary enters the pre-oviposition phase (PP) (Huang 2024). With the continuous enlargement of the ovaries and the full maturation of ovarian function, primordial follicles mature to initiate ovulation (Li and Chian, 2017). When the oocyte enters the oviduct, the tubular glands progressively secrete calcium compounds to form the eggshell and this complete process of follicular development, ovulation and eggshell calcification marks the onset of the oviposition phase (OP) (Nys et al., 2001). The post-laying phase (LP) is marked by distinct ovarian atrophy and functional regression, characterized by arrested follicular development and prevalent follicular atresia, leading to cessation of egg production (Shi et al., 2025).
       
In recent years, transcriptome sequencing (RNA-seq) technology has been extensively applied to screen for genes in the ovaries associated with the egg-laying performance of poultry including chickens, ducks and geese. Mu et al., (2021) used RNA-seq to analyze ovarian tissues from chickens with divergent egg production performances (high-yield and low-yield) to investigate the expression differences of ovarian development-related genes and to elucidate the regulatory mechanisms underlying varying production traits. Bhavana et al., (2022) performed RNA-seq on ovarian tissues of high-and low-egg-producing Indian domestic ducks and identified 38 key candidate genes associated with egg production performance. Nevertheless, research on the expression dynamics of relevant genes during ovarian development across different oviposition phases remains scarce. Therefore, this study used ovarian tissues from WWG during different oviposition phases (PP, OP and LP) as research subjects, employing RNA-seq technology and protein-protein interaction (PPI) network analysis to investigate the gene expression profiles in these tissues across different oviposition phases. This work will serve to identify hub genes and signaling pathways related to the egg production perfor-mance of WWG, reveal the key molecular mechanisms affecting egg production performance in geese and provide an important theoretical basis for molecular breeding of geese.

Ethics statement
 
The experimental procedures were conducted in full compliance with the Guidelines for Laboratory Animal Administration and received official endorsement from the Anhui University of Science and Technology’s Ethics Panel on Animal Studies (Ethical Approval Code: 2024-016).
 
Experimental animals and sample preparation
 
The experimental geese were obtained from Dingyuan Junming Ecological Farm (Dingyuan, Anhui, China) and kept in the same environment with standardized feeding management. Three-year-old WWG female geese were selected as the standardized test group and nine female geese of similar weight under the same feeding conditions were randomly selected at the beginning of the PP (39 months old), the OP (41 months old) and the LP (46 months old). Each group had nine geese and three geese in each group were randomly selected for euthanasia by cervical dislocation, for the collection of organ samples. The ovaries collected during the PP, the OP and the LP were labeled as PP1-3, OP1-3 and LP1-3, rapidly frozen in liquid nitrogen and subsequently transferred to a -80^oC freezer for total RNA extraction and transcriptome sequencing.
 
Transcriptome sequencing and total RNA extraction
 
The RNAprep Pure Tissue Kit (Tiangen Biotech) was used for high-quality RNA extraction from the different ovarian tissues, with experimental workflows strictly adhering to the manufacturer-provided specifications. A Nanodrop 2000 was used to evaluate the RNA purity and concentration from OD260 / 280 readings, while RNA integrity was assessed by 1% agarose gel electrophoresis. Qualified RNA samples were subsequently sent to Gene Denovo Biotechnology Co., Ltd (Guangzhou, China) for sequencing analysis.
 
Transcriptome data analysis
 
To ensure data reliability and accuracy, raw sequencing reads were first processed to remove low-quality bases and filter out poor-quality reads (Xu et al., 2023). After quality control, sequence reads were mapped against the reference genome of WWG (NCBI accession GCF_002166845.1) with the HISAT2 software (version 2.0.4) (Kim et al., 2015) to determine transcript abundance. Transcript reconstruction was performed with StringTie (v1.3.5) (Kovaka et al., 2019), with expression quantification conducted with the FPKM metric (Fragments Per Kilobase per Million mapped reads). Significant differentially expressed genes (DEGs) were selected based on thresholds parameters of |log₂ Fold Change≥1 and P-values < 0.05 (Hu et al., 2020). The original sequencing datasets are publicly accessible through NCBI’s Sequence Read Archive under accession code PRJNA1169106.
 
GO and KEGG enrichment analysis
 
Functional enrichment analysis of the identified DEGs was performed using ClusterProfiler (v3.10.1) (Zhang, 2025) and KOBAS (v2.0) (Xie et al., 2011) for the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses (Pandian et al., 2022). The GO terms and KEGG pathways with p<0.05 were considered significantly enriched, while those with p<0.01 were classified as highly significantly enriched.
 
Protein interaction network analysis and hub gene screening
 
Protein-protein interaction ( PPI ) analysis of the DEGs was performed using the STRING database (http://string-db.org/) (Szklarczyk et al., 2010) with a high-confidence interaction score threshold of 0.700. The resulting PPI network was visualized in the Cytoscape software (v3.8.0) (Raja et al., 2023) and hub genes were identified by applying four topological algorithms (EPC, Degree, MNC and MCC) through the CytoHubba plugin (Shannon et al., 2003).
 
Gene set enrichment analysis (GSEA)
 
Gene Set Enrichment Analysis ( GSEA ) (Jas et al., 2023) was performed on the genome-wide expression profiles to systematically evaluate functional enrichment patterns, including non-significant pathways that might have been overlooked in the conventional KEGG enrichment analysis, while quantitatively assessing the degree of enrichment (Kiser et al., 2025).
 
Extraction of total protein and western blotting detection
 
Ovarian tissue samples from different egg-laying periods of WWGs were placed in pre-cooled mortels and thoroughly ground with liquid nitrogen. The Total protein was extracted from the ovarian tissues of WWGs at different laying periods using the Total protein extraction kit. The protein concentration was determined by the BCA quantitative detection kit (Vazyme, Nanjing, China). The operation was carried out in accordance with the SDS-PAGE and Western blot kit instructions (Abcam plc, Cambridge, UK). After membrane transfer, the color was developed after incubation with the primary antibody and the secondary antibody. The gray values of the protein bands were calculated using Image J 1.39U software (National Institutes of Health, NIH, USA). Data analysis was conducted using the gray value of the target protein/gray value of the internal reference protein (GAPDH) as the relative expression level of the target protein.
 
Validation by real-time quantitative PCR (RT-qPCR)
 
The primer sequences used in this study (Table 1) were designed using the Oligo 7 primer design software (v7.60) (Rychlik, 2007) and all primers were synthesized and purified by Bioengineering (Shanghai Co., Ltd.) following stringent bioinformatics criteria. We selected GAPDH as the reference gene according to previous reports (Gilsbach et al., 2006). The RT-qPCR reaction system was quantified with the SYBR Green qPCR Master Mix (EZBioscience, USA) and the reaction procedure was as follows: pre-denaturation at 95C for 5 min followed by 40 cycles of amplification (95^oC for 10 s, 60^oC for 30 s) and finally, the melting curve analysis was performed to confirm the amplification specificity. Finally, melting curve analysis was performed to confirm the amplification specificity. All samples were set up with three technical repetitions, using the 2^-ΔΔCt method (Livak and Schmittgen, 2001) for relative quantitative analysis in GraphPad Prism 8.0 software (Shannon et al., 2003) to ensure the reliability of the experimental results.

Sequencing data statistics
 
The nine target samples yielded 778,534,602 raw sequencing reads. After quality filtering and GC content distribution checks (Table 2), 775,268,402 high-quality clean reads were retained (average Q30 score: 93.85%). The GC content across samples ranged between 48.77% and 51.61% and reference genome alignment demonstrated a >75% mapping rate for all samples. These statistics suggest that the sequencing quality was good and that these reads could be used for subsequent analyses.


 
Screening and enrichment analysis of DEGs
 
Bar plots and volcano plots were generated to accurately visualize the expression trends of the DEGs in the ovarian tissues of WWG across the three reproductive phases. In the PP vs OP group, 1,701 DEGs were identified of which 667 were upregulated and 1,034 were downregulated (Fig 1A), When the OP group was compared with the LP group, 1 259 DEGs were found of which 871 were upregulated and 388 were downregulated (Fig 1B). The PP vs LP group had 652 DEGs (515 upregulated, 137 downregulated; Fig 1C). Notably, 14 DEGs were consistently differentially expressed across all three comparisons (Fig 1D).


       
The GO enrichment analysis revealed that DEGs in the PP vs OP comparison were significantly enriched in 1,315 functional terms (p<0.01) of which the major enriched categories included: Biological adhesion, Regulation of ion transport, Response to corticosteroid and Regulation of multicellular organismal process (Fig 2A). In the OP vs LP comparison, the 1,219 significantly enriched terms (p<0.01) included: Passive transmembrane transporter activity, Ion channel activity, Response to steroid hormone and Regulation of ion transport (Fig 2B). The PP vs LP comparison had 1,791 significantly enriched terms (p<0.01), dominated by: Biological adhesion, Cell adhesion, Regulation of immune system process and Regulation of cell proliferation (Fig 2C).


       
The KEGG pathway analysis revealed a significant enrichment of DEGs in 18 signaling pathways in the PP vs OP comparison  (p<0.05, Fig 3A), primarily including Neuroactive ligand-receptor interaction, ECM-receptor interaction and Cytokine-cytokine receptor interaction. There was a significant enrichment of DEGs in the OP vs LP comparison across six signaling pathways (p<0.05; Fig 3A) primarily including Neuroactive ligand-receptor, interaction Arachidonic acid metabolism and MAPK signaling pathway (p<0.05, Fig 3B). A comparative analysis of the DEGs between the PP and LP groups revealed significant enrichment in 37 signaling pathways (p<0.05), where the key enriched pathways included: Calcium signaling pathway, PI3K-Akt signaling pathway and Neuroactive ligand-receptor interaction (Fig 3C). Venn diagram analysis of significantly enriched KEGG pathways (p<0.05) revealed that the three groups of signaling pathways were all significantly enriched in Neuroactive receptor-ligand interaction, ECM-receptor interaction and cytokine-cytokine interaction (Fig 3D).


 
Protein-protein interaction networks and identification of hub genes
 
An integration of common DEGs across the three experimental groups and genes from co-enriched signaling pathways was performed, followed by the construction of a protein-protein interaction network using the STRING online database. The resultant network was visualized with Cytoscape and is depicted in Fig 4 and comprised 72 nodes and 82 edges. The CytoHubba module (version 0.1) implemented four network topology metrics for hub gene screening, namely: Maximum Neighborhood Component (MNC), Maximal Clique Centrality (MCC), Node Connectivity Degree and Edge Percolated Component (EPC), ultimately identifying the 10 most significant regulatory nodes (Table 3). Intersection analysis via the Venn diagram revealed five consensus hub genes ( ITGB3, VTN, FN1, ITGA2 and VWF ) identified by all four algorithms (Fig 4). This robust intersection demonstrates their topological centrality within the interaction network (Fig 5), suggesting their potential pivotal roles in regulating WWG egg-laying performance.






 
GSEA analysis
 
Further genome-wide transcriptional profiling was performed using GSEA on all the detectable mRNAs across the distinct egg-laying stages in the ovarian tissues, with the results visualized in Fig 6. In the PP vs OP groups, there was significant enrichment of mRNAs in 11 signaling pathways (Fig 6A), predominantly including Ribosome biogenesis, ECM-receptor interaction and Oxidative phosphorylation (Statistical thresholds: NOM p<0.05, FDR q-val<0.25, |NES|>1). In the OP vs LP comparison, there were 15 significantly enriched pathways (Fig 6B), primarily involving Oxidative phosphorylation, Steroid hormone biosynthesis and Ribosomal function (Statistical thresholds: NOM p<0.05, FDR q-val < 0.25, |NES|>1). In the PP vs LP comparison, transcriptomic profiling identified eight signaling pathways with significantly enriched mRNAs (Fig 6C), including the Ribosomal pathway, ECM-receptor interaction and Cytokine-cytokine receptor interaction (Statistical thresholds: NOM p<0.05, FDR q-val<0.25, |NES|>1).


       
Integrative analysis of the GSEA results across all three experimental groups revealed consistent and significant enrichment of the ribosome pathway ( p<0.05 ) (Fig 7 A-C). The expression of the core-enriched ribosomal genes is shown in Fig 7 (D-F). To elucidate the expression dynamics of the ribosomal genes, hierarchical clustering analysis was performed on these genes across all three experimental groups. The resulting heatmap (Fig 8) and dot plot of DEGs in the ribosomal signaling pathway revealed significantly higher expression levels during the OP compared to both the onset of the laying egg period and the LP ( fold change > 2.0, FDR < 0.01). Notably, 78.6% of the enriched genes in this pathway were ribosomal protein genes, strongly suggesting that ribosomal proteins play critical roles in goose ovarian development and egg-laying.




 
Western blotting of protein
 
The expressions of ACTB, CPEB1, CD36 and ZAR1 in ovarian tissues at different laying periods were detected by Western blotting (Fig 7, Fig 8). It was found that the expression of ACTB protein during the PP was significantly lower than that during the OP and the expression of CPEB1 protein during the PP was significantly higher than that during the OP (P<0.05). Secondly, the expression of CD36 protein during the OP was significantly higher than that during the CP and the expression of ZAR1 protein during the OP was significantly lower than that during the CP (P<0.05).
 
Real-time fluorescence quantitative PCR verification (RT-qPCR)
 
To validate the reliability of the transcriptome sequencing data, RT-qPCR analysis was performed on five randomly selected DEGs. The results demonstrated significant concordance (p<0.05) between the expression patterns in the ovarian tissues and the RNA-seq data, confirming the accuracy and reproducibility of our transcriptomic findings (Fig 9).


 
Western-blot was used to detect the expression of ACTB, CPEB1, CD36 and ZAR1 proteins in ovarian tissues at different oviposition stages
 
The expression of ACTB, CPEB1, CD36 and ZAR1 proteins in the ovary of Wanxi White Goose at different spawning stages was determined by Western-blot. Fig 10 is the result of Western-blot and the relative content of each protein in different samples is shown by bands. The results showed that the expression levels of ACTB and CPEB1 in OP were significantly higher than those in PP ( * indicated statistical difference), the expression level of CD36 in OP was significantly higher than that in CP and the expression level of ZAR1 in CP was significantly higher than that in OP, which intuitively showed the dynamic changes of each protein in different spawning periods (Fig 11).




       
The ovary is a dynamically developing organ that constitutes a critical component of the female reproductive system, which in turn plays a decisive role in regulating avian egg-laying performance (Chen et al., 2021). Empirical studies have demonstrated that during the PP in geese, the ovarian tissue enlarges significantly, with the ovarian surface exhibiting numerous developing follicles (Zhao et al., 2022). During the OP, the ovarian tissue volume of geese further increases, with the ovarian surface becoming densely populated by follicles of varying diameters (Hu et al., 2021). During the resting period, the ovarian surface of geese develops characteristic depressions with evident follicular atresia (apoptotic index > 30%), leading to the gradual cessation of oviposition (Li et al., 2024). These morphological changes confirm significant developmental differences in follicular dynamics across distinct laying phases. However, there exists a significant paucity of research elucidating the gene expression dynamics governing follicular tissue development across these distinct laying phases in geese. Therefore, this study used ovarian tissues from WWG at distinct laying phases and with RNA-sequencing technology, systematically investigated the gene expression dynamics during ovarian development across different reproductive stages. The identification of key candidate genes and signaling pathways associated with egg-laying performance in WWG will not only facilitate genetic improvement of this breed’s reproductive traits but also provide a molecular basis for precision breeding strategies. The results showed the presence of 1,701 (PP vs OP), 1,259 (OP vs LP) and 652 (PP vs LP) DEGs (Fig 2). Functional enrichment analysis using the Gene Ontology database revealed that the DEGs across OP were significantly enriched (p<0.01, FDR < 0.05) in key biological processes including Biological adhesion, Cell adhesion, Regulation of cell proliferation and Response to corticosteroid. KEGG pathway analysis of DEGs identified that three signaling pathways found to be significantly enriched coexisted in the three groups, including the neuroactive receptor-ligand interaction, ECM-receptor interaction and cytokine-cytokine interaction signaling pathways. The neuroactive ligand-receptor interaction pathway regulates the specific binding of ligands such as neurotransmitters and neuropeptides to their receptors, including G protein-coupled receptors (Gruber et al., 2010) and ion channel receptors (Su et al., 2009), thereby modulating signal transduction and physiological functions. Furthermore, the neuroactive ligand-receptor interaction pathway participates in the signal transduction of the hypothalamic-pituitary -gon -adal (HPG) axis, regulating hormone secretion and reproductive processes and thus mediating avian egg-laying performance (Yan et al., 2022). Follicular development, maturation, ovulation and corpus luteum formation are regulated by the ECM-receptor interaction pathway, which mainly facilitates specific binding and signaling between the extracellular matrix and cell surface receptors (Guo et al., 2022; Hrabia, 2021; Kulus et al., 2021). The Cytokine-cytokine receptor interaction pathway primarily regulates signal transduction between diverse cytokines, influencing ovarian tissue development and function (Yuan et al., 2025). The findings suggest a coordinated role for neuroactive ligand-receptor, ECM-receptor and cytokine-receptor interactions in regulating follicular development and egg-laying, impacting reproductive hormone levels and intercellular signaling throughout the reproductive cycle.
       
To further investigate the gene expression patterns in WWG across different oviposition phases, GSEA was performed on the three reproductive phases, with results demonstrating significant enrichment (p<0.05) of the Ribosome pathway in all three groups. The ribosome is essential for protein synthesis, translating the genetic information encoded in mRNAs into polypeptide chains through precise amino acid polymerization, followed by protein folding (Hong et al., 2024). Analysis of the core-enriched gene expression in the Ribosome pathway revealed extensive expression of ribosomal protein (RP) gene family members. Consistent with the results of Wang et al., (2023) and Jiang et al., (2024), the core-enriched genes RPL15 and RPS24 were highly expressed in ovarian tissues during the OP, with lower expression during the PP and complete transcriptional silencing in the LP. The results of GSEA further verified the reliability of KEGG pathway analysis, in which the ribosome pathway, the neuroactive ligand-receptor interaction pathway, the ECM-receptor interaction pathway and the cytokine-cytokine receptor interaction pathway were all significantly enriched, suggesting their pivotal roles in the reproductive regulation of WWG. The DEGs within these pathways exhibited dynamic expression patterns in ovarian tissues across different oviposition phases. The DEGs in the Ribosome pathway are primarily involved in protein synthesis within reproductive cells, while those in the neuroactive ligand-receptor interaction pathway mainly participate in reproductive hormone signal transduction. Genes differentially expressed in the ECM-receptor interaction pathway mainly control the follicular microenvironment via paracrine signaling. Therefore, the coordinated expression patterns of these genes underlie the molecular mechanisms regulating reproductive cycles and maintaining egg-laying performance in WWG.
       
A PPI network was constructed using the screened DEGs, with high-confidence interactions selected to generate the core network. Here, the CytoHubba plugin identified ITGB3, VTN, FN1, ITGA2 and VWF as the hub genes with crucial regulatory roles in the network topology. The ITGB3 (integrinβ3) gene is a pivotal member of the integrin family (Zhu et al., 2019) and along with its ligand VTN (vitronectin), they serve as critical receptor-ligand pairs mediating ECM-receptor interactions (Liu et al., 2015). Within the follicular microenvironment, ITGB3 activates and binds to VTN to regulate angiogenesis and cellular differentiation during ovarian development (Kulus et al., 2019). Furthermore, VTN modulates multiple reproductive developmental pathways, including ovarian morphogenesis and neurogenesis, through its interaction with the ITGB3 receptor (Chen et al., 2020). Here, VTN expression was significantly higher during the OP compared to the PP and LP. Conversely, ITGB3 exhibited its lowest expression during the OP, with intermediate expression in the PP and the highest expression during the LP. These inverse expression patterns suggest a potential regulatory mechanism whereby ITGB3 activation downregulates its expression post-VTN stimulation, while VTN-receptor binding modulates ovarian development and egg production through reproductive developmental pathways. Also, FN1 (fibronectin 1), encoding an extracellular matrix (ECM) glycoprotein ligand, plays essential roles in focal adhesion formation and ECM-receptor interactions within theca and ovarian granulosa cells, thereby regulating follicular growth, maturation and functional maintenance (Leng et al., 2024). It exhibited significantly higher expression during the OP compared to the PP and LP. This expression pattern suggests its functional role in providing nutritional support for oocyte development and maintaining follicular morphological integrity.
               
The ITGA2 gene specifically interacts with collagen and focal adhesion components within the extracellular matrix (ECM) to preserve tissue structure and granulosa cell integrity (Clark and Brugge, 1995). In this study, ITGA2 expression was significantly elevated during the OP compared to the PP and LP. These findings suggest ITGA2 may regulate granulosa cell proliferation and differentiation to promote follicular growth and development. Thus, the elevated expression of ITGB3, VTN, FN1 and ITGA2 within the ECM-receptor interaction pathway during the oviposition phase suggests that their encoded proteins act as ligands or receptors to cooperatively maintain cell-ECM homeostasis and promote ovarian development, thereby modulating reproductive behaviors in geese (Kulus et al., 2021). The VWF (von Willebrand factor) gene encodes a large multifunctional glycoprotein synthesized by endothelial cells and megakaryocytes, which serves as a biomarker for endothelial cells (Mojzisch and Brehm, 2021). Keesler et al., (2021) demonstrated that VWF directly interacts with fibronectin, a key ECM component, regulating platelet adhesion and vascular endothelial function. Significant changes in follicular morphology are preceded by an abnormal blood supply, which becomes an important factor in follicular atresia and occurrence (Isobe et al., 2001). In this study, VWF was expressed less during the PP compared to its elevated expression levels in the OP and LP. This upregulation suggests potential associations with vascular network disruption and endothelial cell impairment during follicular atresia, which may adversely affect egg-laying performance and follicular development. This study’s identification of five hub genes (ITGB3, VTN, FN1, ITGA2 and VWF) suggests they stimulate ovarian development and consequently boost goose reproductive performance and egg production by activating the ECM-receptor interaction pathway. Concurrently, GSEA revealed consistent significant enrichment of the Ribosome pathway across all the experimental groups, suggesting that ribosomal protein (RP) gene families within this pathway may drive ovarian and follicular development through enhanced protein biosynthesis, thereby regulating egg-laying performance and influencing total egg output. 

This study employed RNA-seq technology to analyze ovarian tissues of WWG across distinct oviposition phases. Five hub genes (ITGB3, VTN, FN1, ITGA2 and VWF) exhibiting strong correlations with egg-laying performance were identified. Among them, the expression levels of FN1, ITGA2 and VWF genes were the highest during the OP, the expression level of ITGB3 gene was the highest during the PP and the expression level of VTN gene was the highest during the LP. Through the genome-wide GSEA of all mRNAs, the Ribosome pathway was identified as a key signaling pathway associated with egg-laying performance. These findings provide a theoretical foundation for providing an important reference for the subsequent use of gene editing technology to breed new varieties and promote the development of WWG breeding industry.

This work was supported by grants from the talent introduction project of Anhui Science and Technology University [DKYJ202105, DKYJ202104], Anhui Province Science and Technology Major Project [17030701004], Local goose gene bank in Anhui Province, Science and Technology Project of Chuzhou City, Anhui Province [2022ZN002], Veterinary Science Peak Discipline Project of Anhui Science and Technology University [XK-XJGF002].
 
Disclaimers
 
The views and conclusions expressed in this article are solely those of the authors and do not necessarily represent the views of their affiliated institutions. The authors are responsible for the accuracy and completeness of the information provided, but do not accept any liability for any direct or indirect losses resulting from the use of this content.
 
Informed consent
 
All animal procedures for experiments were approved by the Institutional Animal Care and Use Committee of Anhui Science and Technology University, Chuzhou, China.

The authors declare that there are no conflicts of interest regarding the publication of this article.