The expression profiling dataset GSE87201 for Homo sapiens, pertinent to DOR, was retrieved from the GEO database (Barrett et al., 2007) utilizing the R package GEOquery (Davis and Meltzer, 2007). The dataset comprised thirty-five samples, categorized into human oocytes from individuals aged 20–22 years and 33–35 years. Specifically, the 20–22 age group included seven samples with low ovarian reserve (LOR) and 9 with high ovarian reserve (HOR), while the 33–35 age group comprised nine LOR and nine HOR samples. The dataset probes were annotated using the microarray platform files GPL17586 to incorporate expression profiling data from the LOR and HOR samples aged 33–35 years for subsequent analyses. Hypoxia-related genes were identified using “hypoxia” as a search term in the GeneCards database (https://www.GeneCards.org/) and corroborated by articles published in PubMed. Supplementary Tables S1 and S2 provide a comprehensive list of genes implicated in DOR and hypoxia. This study utilized exclusively publicly available de-identified genomic data from the GEO database (accession GSE87201). No primary human or animal samples were collected by the authors. The original data generation by Barragán et al. obtained appropriate ethical oversight as stated in their publication (Barragan et al., 2017), including Institutional Review Board approval and participant informed consent. To account for potential batch effects across sample processing batches, the ComBat function from the R package was applied to the normalized expression matrix. Batch correction was performed prior to downstream differential expression and PCA analysis. The effectiveness of batch effect removal was evaluated using principal component analysis (PCA), comparing sample clustering before and after correction. The experimental workflow is depicted in Figure 1.

The GSE87201 dataset was standardized utilizing the R package limma (Ritchie et al., 2015), applying a threshold of |log fold change (FC)| > 1 and P < 0.05. Subsequently, the dataset underwent standardization and differential expression analysis to identify differentially expressed genes (DEGs) between LOR and HOR conditions. Volcano plots were generated to visualize upregulated genes (logFC >0) and downregulated genes (logFC <0). The standardized dataset was further analyzed for differential expression, identifying Co-DEGs between LOR and HOR, which were then intersected with hypoxia-related genes to derive HRDEGs. Venn diagrams were employed to illustrate the HRDEGs. The outcomes of the differential expression analysis were visualized through volcano plots and heatmaps, created using the ggplot2 and pheatmap packages in R.

Genes within the GSE87201 dataset were categorized into two groups based on the sign of their logFC values. GSEA was conducted on the DEGs using the R package clusterProfiler, with significance criteria set at p. adj <0.05 and a false discovery rate (FDR) q-value <0.25 (Subramanian et al., 2005; Wu T. et al., 2022; Liu et al., 2023). All genes within the LOR and HOR regions were subjected to analysis utilizing GSVA enrichment (Kowalski et al., 2016). The screening criteria employed were a p-value of less than 0.05 and an absolute log fold change (|logFC|) greater than 0.5, aimed at investigating biological processes, cellular composition, and molecular functions (MF). GSEA and GSVA were conducted on the entire gene expression dataset rather than DEGs alone, to capture pathway-level changes that may arise from coordinated but modest shifts in gene expression, in line with established best practices.

ROC curves for the HRDEGs were generated using the pROC package in R software through univariate logistic regression models (Mandrekar, 2010). For each HRDEGs, gene expression levels served as the predictor variable, with ovarian reserve status (LOR/HOR) as the binary outcome. Predicted probabilities derived from these models were used to construct ROC curves. The area under the curve (AUC) was computed to evaluate the diagnostic efficacy of individual HRDEGs for DOR, with 95% confidence intervals calculated via 2,000 bootstrap resamples. AUC values range from 0.5 to 1, where values approaching one indicate superior diagnostic performance: AUC 0.5–0.7 denotes low accuracy, 0.7–0.9 moderate accuracy, and >0.9 high accuracy.

PPI analysis of HRDEGs was performed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (Szklarczyk et al., 2019), with the confidence level set to medium (0.400). Cytoscape (Shannon et al., 2003) (version 3.9.1) was employed to visualize the PPI network model, designating HRDEGs as hub genes.

The miRDB database (Chen and Wang, 2020) was utilized to predict micro-ribonucleic acids (miRNAs) that interact with these hub genes. The interaction network between messenger ribonucleic acid (mRNA) and miRNA was constructed by extracting segments of the database with a target score exceeding 80. The prediction of RNA binding protein (RBP) interactions with hub genes was performed using the ENCORI database (https://starbase.sysu.edu.cn/) (Li et al., 2014). The mRNA–RBP interaction network was then delineated with parameters set at ClusterNum > 6 and clipExpNum > 6. Transcription factors (TFs) interacting with hub genes were identified through the CHIPBase database (version 3.0) (https://rna.sysu.edu.cn/chipbase/) (Zhou et al., 2017) and the hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget) (Zhang et al., 2020). Furthermore, potential drugs or small molecules interacting with the hub genes were identified via the Comparative Toxicogenomics Database (CTD) (http://ctdbase.org/) (Davis et al., 2021), and the mRNA-TF and mRNA–drug interaction networks were visualized using Cytoscape software.

The GO (Gaudet et al., 2021) and KEGG (Kanehisa et al., 2017) analyses are extensively utilized databases for the storage of genomic, biological pathway, disease, and drug information. GO annotation analysis of hub genes was conducted using the R package clusterProfiler (Yu et al., 2012). The criteria for entry screening were set at a p-adjusted value of less than 0.1 and an FDR q-value of less than 0.05.

The matrix data from the DOR datasets, specifically LOR and HOR, were analyzed for immune cell enrichment scores greater than zero utilizing the CIBERSORT package (Newman et al., 2015) in conjunction with the LM22 feature gene matrix. This approach facilitated the acquisition and presentation of specific results pertaining to the immune cell infiltration abundance matrix. Correlations between various immune cells within the LOR and HOR datasets were visualized using the Spearman algorithm and the R package ggplot2. Additionally, correlations between immune cells and key genes were determined by integrating the gene expression matrix from the DOR dataset, which enabled the creation of hotspot maps via the R package ggplot2.

KGN cells were seeded in 96-well plates at a density of 5 × 10⁵ cells per well and incubated overnight to allow for wall attachment. Following passaging and cell attachment, 4-HC was introduced into the target culture medium at final concentrations of 0.2, 2.0, 6.0, and 10.0 μmol/L, and the cells were subsequently mixed and incubated for 48 h. Cell viability was assessed using the CCK-8 assay, while apoptosis was evaluated using Annexin-V-FITC/PI-PE staining. This process facilitated the identification of optimal concentrations of 4-HC for the model of 4-HC-induced granulosa cell response.

In this study, distinct groups of cells in the logarithmic growth phase were seeded at a density of 2,000 cells per well in 96-well plates. The cells were incubated at 37°C in a humidified atmosphere containing 5% CO₂ for 6 hours to allow for adhesion to the well surfaces. Following attachment, various treatments were administered according to the designated experimental groups, with six replicate wells established for each condition. Two hours prior to the conclusion of the incubation period, 10 μL of CCK-8 solution was added to each well. At the end of the incubation, optical density at 450 nm (OD450) was measured using a microplate reader.

Subsequently, the cell suspension was transferred to a centrifuge tube containing an appropriate volume of culture medium, centrifuged at 1,000 rpm for 5 min, and the supernatant was discarded. The cell pellet was washed once with phosphate-buffered saline (PBS). Annexin V Binding Solution was then added to achieve a final cell concentration of 1 × 10⁶ cells/mL. The cell suspension was mixed with Annexin V-FITC conjugate and propidium iodide (PI) solution, followed by the addition of 1X Annexin V Binding Solution. The samples were analyzed using flow cytometry within 1 hour.

For further analysis, cells were cultured in 6-well plates until they reached 60%–70% confluence. The cells were washed twice with PBS, collected, fixed, and subsequently pelleted. Approximately 1 mL of ice-cold PBS was added to resuspend the cells. Centrifuge the cells once more and carefully remove the supernatant, ensuring approximately 50 μL of PBS remains to prevent cell aspiration. Gently tap the bottom of the centrifuge tube to adequately disperse the cells and prevent clumping. For optimal staining results, it is recommended to complete the flow cytometry assay within 24 h. Red fluorescence should be detected at an excitation wavelength of 488 nm using flow cytometry, concurrently measuring light scattering.

Total RNA was extracted from tissue samples utilizing TRIzol reagent (Invitrogen) and quantified via NanoDrop spectrophotometry. High-quality RNA was subsequently reverse-transcribed into complementary DNA (cDNA) employing the iScript™ cDNA Synthesis Kit (Bio-Rad) in accordance with the manufacturer’s instructions. The qRT-PCR was conducted using SYBR^® Green Master Mix (Applied Biosystems) on a StepOnePlus™ Real-Time PCR System. The reaction mixture, with a total volume of 20 μL, comprised 1 µL of cDNA, 10 µL of SYBR Green Master Mix, and 0.5 µM of each primer. The cycling conditions were as follows: initial denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Specificity was verified through melting curve analysis. Primer sequences were designed utilizing Primer-BLAST software and subsequently validated for efficiency. The expression levels of target genes were normalized to GAPDH employing the 2^−ΔΔCt method. Each sample was analyzed in triplicate.

Proteins were extracted from tissues or cells using RIPA buffer supplemented with protease inhibitors. Protein concentrations were quantified using the BCA Protein Assay Kit. Equal amounts of protein were resolved by SDS-PAGE on 10% polyacrylamide gels and transferred onto PVDF membranes. Membranes were blocked with 5% non-fat milk in TBST for 1 hour at room temperature, followed by overnight incubation at 4°C with primary antibodies against the target proteins at a 1:1,000 dilution. Following washing, membranes were incubated with HRP-conjugated secondary antibodies at a 1:5,000 dilution for 1 hour at room temperature. Detection of signals was performed using an ECL chemiluminescence kit, and imaging was conducted with a ChemiDoc imaging system. β-actin served as an internal control to normalize protein loading. Each experiment was conducted at least three times to ensure reproducibility.

The data analysis was rigorously performed utilizing R software (4.3.3), a robust platform for statistical computing. Continuous variables were reported as means with their corresponding standard deviations, offering a comprehensive depiction of the data’s variability. The Wilcoxon rank-sum test was applied to compare different groups, thereby accommodating non-parametric data distributions. Additionally, all results not explicitly specified were calculated using the Spearman test, which is particularly advantageous for evaluating correlations among ranked variables. A p-value of less than 0.05 was considered statistically significant, suggesting that any observed differences or associations were unlikely to be due to random chance, thereby enhancing the credibility of the findings.

The DOR GSE87201 dataset underwent normalization utilizing the R package ‘limma.’ Following the removal of intersample batch effects, nine HOR and nine LOR samples were retained (Figures 2A,B). Validation of the expression matrix through PCA indicated that the batch effect was substantially mitigated after ComBat correction (Figures 2C,D). Specifically, Figure 2C shows batch-driven clustering before correction, while Figure 2D shows the dispersion of samples no longer follows batch grouping, reflecting effective batch adjustment. A total of 28,136 DEGs were identified, with 1,095 genes exhibiting a |logFC| greater than 0 and a p-value less than 0.05. Among these, 618 genes were upregulated, showing higher expression in LOR (and lower expression in HOR with LogFC > 0), while 477 genes were downregulated, exhibiting lower expression in LOR (and higher expression in HOR with LogFC < 0), as illustrated in the volcano plot (Figure 2E).

Based on the intersection of HRGs with Co-DEGs, we identified the following genes: AMPD3, FANCI, FDPS, FOXD1, GNPDA1, GPC3, HMOX1, JMJD6, KAT2A, MRTO4, SLC22A5, SMAD5-AS1, TACC3, TPX2, TTC13, TTF2, VHL, and WSB1. In total, 18 HRDEGs were depicted in a Venn diagram (Figure 2F).

From the heatmap (Figure 2G), it can be observed that there are differences in the expression patterns of various genes between the LOR and HOR groups, which may be related to ovarian function and egg quality. For instance, some genes such as GNPDA1, JMJD6, SLC22A5, FDPS, HMOX1, and FOXD1 exhibit higher expression in the HOR group, which may be associated with better ovarian reserve and egg quality. Conversely, other genes such as WSB1, VHL, FANCI, TACC3, TTC13, AMPD3, GPC3, MRTO4, KAT2A, TPX2, and TTF2 show higher expression in the LOR group, potentially correlating with reduced ovarian reserve and decreased egg quality (Figure 2G).

The GSEA of all genes from the GSE87201 dataset, utilizing the FDR p-value (q-value <0.25 and adjusted p-value < 0.05) as the criterion for enrichment screening, identified significant enrichment of DEGs in five biological processes: response to arsenite uptake (Pan et al., 2020; Chen et al., 2022; Eslami et al., 2024; Adeogun et al., 2025), zinc homeostasis, oxidative phosphorylation (Verma et al., 2018; Nasiadek et al., 2020; Camp et al., 2023; Liao et al., 2023; Liu et al., 2024), respiratory electron transport coupled with adenosine triphosphate synthesis (Anderson et al., 2018; Bao et al., 2024; Kobayashi et al., 2024), and hypoxia metagene (Lim et al., 2021; Huang et al., 2023; Gutzeit et al., 2024) (Figure 3). The GSVA identified 11 gene set pathways exhibiting distinct expression differences between the low and high outcome risk (LOR and HOR) groups, which were visualized using a heat map generated with the R package ‘pheatmap’ (Figure 4A). A comparative analysis of the degree of difference across these 11 pathways, conducted using the Mann–Whitney U test, indicated that eight pathways were statistically significant (p < 0.05) (Figure 4B). Detailed information on the five biological processes and the 11 gene set pathways is provided in Supplementary Tables S3 and S4.

Furthermore, the analysis of expression differences in 18 HRDEGs between the LOR and HOR groups, performed using the Wilcoxon signed-rank test, revealed that 12 HRDEGs were statistically significant (Figure 5A). ROC curves and AUC values for the 18 HRDEGs within the GSE87201 dataset were plotted individually (Figure 5B). The findings indicated a significant correlation between 17 HRDEGs and DOR.

To further elucidate the underlying mechanisms of hypoxia-related markers in DOR, we performed a protein-protein interaction (PPI) analysis on 18 HRDEGs within the DOR dataset, utilizing the STRING database with a medium confidence threshold of 0.400. The resultant data were imported into Cytoscape software to generate an interrelationship network diagram, which identified six hub genes: FANCI, KAT2A, TACC3, TPX2, VHL, and WSB1 (Figure 6A).

The miRNAs interacting with these hub genes were predicted using the miRDB database, leading to the construction of an mRNA-miRNA regulatory network comprising 67 miRNA molecules and 70 mRNA-miRNA interaction pairs, as detailed in Supplementary Table S5 and depicted in Figure 6B. Additionally, the interactions between hub genes and RNA-binding proteins (RBPs) were predicted using the ENCORI database, revealing an mRNA–RBP network consisting of 23 RBP molecules and 77 mRNA-RBP interaction pairs (Figure 6C), with TPX2 interacting with 17 RBP molecules, as shown in Supplementary Table S6. Transcription factors binding to the hub genes were identified through the CHIPBase (version 2.0) and hTFtarget databases. The intersection facilitated the construction of a visualization network comprising 28 transcription factors and hub genes (Figure 6D), with KAT2A exhibiting 12 interaction pairs with transcription factors (Supplementary Table S7). Subsequently, the CTD database was utilized to predict 111 potential drugs or molecular compounds targeting the hub genes (Figure 6E), among which 72 were identified as targeting TPX2 (Supplementary Table S8).

Collectively, these multi-omics analyses delineate a comprehensive regulatory landscape for the six hub genes, encompassing post-transcriptional control (miRNAs), RNA processing (RBPs), transcriptional regulation, and druggable targets. Notably, TPX2 emerges as the most interconnected node with both RBPs (17 partners) and therapeutic compounds (72 agents), highlighting its central role in DOR pathogenesis and potential as a priority target for intervention.

GO analysis, employing a screening threshold of p < 0.1 (Lu et al., 2004; Katsel et al., 2005; Alam-Faruque et al., 2014) and an FDR <0.05, revealed the enrichment of hub genes in DOR (Figure 7A). This includes genes associated with the regulation of microtubule-based processes, mitotic spindle organization, and microtubule cytoskeleton organization during mitosis, as well as other biological processes (Figure 7B). Additionally, the analysis identified cellular components (CCs) such as the mitotic spindle, spindle pole, transcription factor target coactivator complex, Ada2a-containing complex, and spindle cellular components (Figure 7C). Furthermore, MFs such as H4 histone acetyltransferase activity, DNA polymerase binding, ubiquitin-protein transferase activity, ubiquitin-like protein transferase activity, and DNA-binding transcription factor binding were enriched (Figure 7D).

The KEGG analysis revealed significant enrichment of hub genes within five specific pathways: the Fanconi anemia pathway, Notch signaling pathway, viral life cycle-HIV-1, renal cell carcinoma, and HIF-1 signaling pathway. These findings are illustrated through bar graphs (Figure 7E) and network graphs (Figure 7F). Detailed results of the GO and KEGG analyses can be found in Supplementary Tables S9 and S10.

Additionally, we employed the CIBERSORT package alongside the Pearson algorithm to assess the correlation of 22 immune cell types with the DOR GSE87201 dataset (Figures 8A–C) and to examine the relationship between hub genes and 20 immune cell types in LOR and HOR. In LOR, the strongest positive correlation was observed between M2 macrophages and gamma delta-T cells, while the strongest negative correlation was found between plasma cells and CD8 T cells. Notably, FANCI and TACC3 exhibited significant positive associations with eosinophils, naïve CD4 T cells, and follicular helper T cells. Conversely, KAT2A showed a significant negative correlation with resting dendritic cells. Furthermore, CD4 T memory cells were significantly negatively correlated with resting dendritic cells and activated CD4 memory T cells, while demonstrating a positive correlation with follicular helper T cells. TPX2 exhibited a positive association with eosinophils, naïve CD4 T cells, and follicular helper T cells, while demonstrating a negative correlation with activated CD4 memory T cells (Figure 8D). Within the HOR group, WSBI showed a positive correlation with resting mast cells and regulatory T cells. Conversely, VHL was negatively associated with eosinophils and resting CD4 T memory cells, yet positively correlated with regulatory T cells. KAT2A expression was positively correlated with naïve B cells and negatively correlated with M1 macrophages. Additionally, TPX2 expression was negatively correlated with M1 macrophages (Figure 8E). In summary, FANCI, TACC3, and TPX2 were significantly associated with various immune cells. These robust correlations delineate distinct immune microenvironment patterns associated with key hub genes in both LOR and HOR contexts. Notably, TPX2 consistently demonstrated significant associations across both groups, particularly with T cell subsets and eosinophils in LOR and M1 macrophages in HOR, suggesting its potentially broader role.

To validate the results of the bioinformatics analysis, a 4-HC-induced granulosa cell injury model was employed. Among the six identified hub genes, FANCI and KAT2A were randomly selected for preliminary validation to provide a representative functional assessment. Flow cytometry analysis revealed that the apoptosis rate of KGN cells increased with higher drug concentrations, while cell viability decreased, as determined by the CCK8 activity assay (Figures 9A,B). Compared to the Control group, the mRNA expression levels of FANCI and KAT2A were significantly elevated in the DOR group. In contrast, the mRNA expression level of FANCI was significantly reduced in the DOR + Si-FANCI group, and the mRNA expression level of KAT2A was significantly reduced in the DOR + Si-KAT2A group, relative to the DOR + Si-NC group (Figures 9C,D). These findings were corroborated at the protein level (Figure 9F). Furthermore, cell viability in the DOR group was significantly lower compared to the Control group, whereas cell viability in the DOR + Si-FANCI and DOR + Si-KAT2A groups was significantly higher compared to the DOR + Si-NC group (Figure 9E).

Furthermore, the percentage of apoptosis was significantly elevated in the DOR group compared to the Control group, and significantly reduced in the DOR + Si-FANCI and DOR + Si-KAT2A groups compared to the DOR + Si-NC group (Figure 10A). In the DOR group, there was a significant increase in the number of cells in the G0/G1 phase and a significant decrease in the number of cells in the S/G2 phase compared to the Control group. Conversely, in the DOR + Si-NC group, the number of cells in the G0/G1 phase significantly decreased, while the number of cells in the S/G2 phase significantly increased compared to the DOR + Si-FANCI and DOR + Si-KAT2A groups (Figure 10B).

Collectively, these rigorous in vitro functional assays robustly validate our bioinformatic predictions, demonstrating that both FANCI and KAT2A are functionally implicated in granulosa cell injury associated with DOR. Silencing either gene effectively mitigated the detrimental effects of 4-HC, as evidenced by significantly reduced apoptosis, restored cell viability, and normalized cell cycle progression (reduced G0/G1 arrest and increased S/G2 phase cells) compared to the non-targeting siRNA control in the injury model. These findings solidify FANCI and KAT2A as critical drivers of granulosa cell dysfunction in DOR and underscore their potential as therapeutic targets for mitigating ovarian damage.