Based on a comprehensive literature review using databases such as PubMed and Web of Science, we selected structural model of Bisphenol A (BPA) to gain precise insights into their reproductive toxicity. Prior literature searches indicate that BPA demonstrate strongly associated with hormonal axis dysregulation and exhibit marked nonlinear correlations with PCOS (Xu et al., 2024). The Comparative Toxicogenomics Database (CTD; http://ctdbase.org/) was utilized to identify potential hub genes of BPA and screen for shared genes between BPA and PCOS (CBP) (Supplementary Table S1). The research steps are illustrated in Supplementary Figure S1.

CBP genes identified from the CTD database were submitted to the STRING platform (https://string-db.org) to retrieve high - confidence protein interactions (confidence score >0.7). The resultant Protein - Protein Interaction (PPI) network was visualized using Cytoscape software (v3.9.1), enabling graphical analysis of molecular associations.

To characterize the biological functions and pathways linking BPA to PCOS pathogenesis, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on the CBP gene set using the ClusterProfiler R package (v.4.10.1). Terms with Bonferroni - adjusted p < 0.05 were considered statistically significant, delineating key molecular mechanisms and signaling cascades. Details are provided in Supplementary Table S2.

PCOS - related datasets from granulosa cells were retrieved from GEO. Data normalization and standardization were done using the limma R package. Three datasets (GSE34526, GSE10946, GSE102293) formed the primary analysis dataset, while two others (GSE106724, GSE98595) were consolidated and normalized as validation datasets. The limma package (v.3.58.1) was used for differential expression analysis between control and PCOS groups; genes with |logFC| ≥ 0.5 and p < 0.05 were selected as DEGs for further study. Oxidative stress- and apoptosis-related genes were downloaded from GeneCards. Venn diagrams visualized overlaps among CBP, DEGs, and the two phenotype-specific gene sets. Differential analysis results were shown via a ggplot2-generated volcano plot and pheatmap-created heatmap.

A nomogram for predicting PCOS risk was developed using hub genes with the rms (v.6.7.1) package. Its predictive accuracy was assessed through receiver operating characteristic (ROC) curve analysis, which was generated by the pROC (v.1.18.5) and ggplot2 (v.3.5.2) packages. Calibration curve analysis (rms and nomogramFormula packages, v.1.2.0.0) assessed model accuracy, while decision curve analysis evaluated clinical utility and net benefit across threshold probabilities.

Single-gene GSEA was performed using GSEA software to explore the biological functions of BPA-exposure-related genes in PCOS progression. Samples were divided into high- and low-expression groups based on median gene expression levels. Molecular pathways associated with gene expression profiles and phenotypic groupings were analyzed using the “c2.cp.kegg.v7.4.symbols.gmt” subset from the Molecular Signatures Database (MSigDB), with results visualized via ggplot2.

Immune cell subtype abundances were quantified using single-sample gene set enrichment analysis (ssGSEA) with predefined gene signatures. Spearman correlation analysis was performed to assess relationships between key gene expression levels and ssGSEA-derived immune cell scores, with results visualized using the ggplot2 R package. Additionally, the CIBERSORT algorithm was applied to deconvolute immune cell compositions from bulk gene expression data, comparing PCOS patients to controls. Statistical significance was set at p < 0.05, and intercellular correlations were visualized using heatmaps generated with ggplot2.

Molecular docking was utilized to evaluate the binding affinity between BPA and five proteins. Three-dimensional structures of BPA were acquired from PubChem (https://pubchem.ncbi.nlm.nih.gov/), whereas crystal structures of five proteins were downloaded from the RCSB PDB (https://www.rcsb.org/). After removing water molecules and small ligands, the CB-DOCK platform (https://cadd.labshare.cn/cb-dock2/) was employed to predict binding sites and affinity for protein-ligand complexes. A protein-ligand blind docking strategy was implemented, integrating cavity detection, docking simulations, and homologous template fitting. Two-dimensional receptor-ligand interaction maps were generated using Discovery Studio software to improve visualization of docking results. Outcomes were visualized and interpreted with PyMOL to clarify structure-activity relationships at the molecular level.

The scRNA-seq dataset GSE268919 was obtained from a published mouse ovarian PCOS-like model (DHEA-induced) including PCOS-like and control mice (Luo et al., 2024). After standard quality control and preprocessing in Seurat (v5.3.0), cells were filtered according to the published pipeline (e.g., thresholds for UMI/gene counts, mitochondrial content, and dissociation-related gene signatures), followed by normalization and dimensionality reduction. Cell clusters were annotated based on canonical marker genes to define major ovarian cell types, and hub-gene expression patterns were then examined in a cell-type-resolved manner. Notably, gene symbols in the scRNA-seq section follow mouse nomenclature (Cyp19a1, Dmd). Differential expression analysis was performed using FindAllMarkers (min.pct = 0.25, logfc.threshold = 0.25) to identify cluster-specific markers.

Six PCOS patients (diagnosed per Rotterdam criteria) and six healthy controls (undergoing oocyte retrieval for tubal infertility) were enrolled. Exclusion criteria included other ovarian disorders, non-PCOS hyperandrogenism, metabolic/gynecological diseases, and organ dysfunction. Baseline demographic and clinical characteristics of the clinical validation cohort are summarized in Supplementary Table S3. The two groups were generally comparable, with no statistically significant differences in baseline variables except AMH. The study was approved by the Ethics Committee of Fujian Maternity and Child Health Hospital (No. 2025KY222) with written informed consent from all participants, complying with the Declaration of Helsinki. Follicular fluid was collected during oocyte retrieval, centrifuged at 300 × g for 10 min. The supernatant was mixed with 50% Percoll (Solarbio, Beijing) and re-centrifuged at 400 × g for 20 min. Isolated granulosa cells were stored at −80 °C for subsequent experiments.

The human ovarian granulosa cell line KGN (Shuochengbio, Shanghai) was cultured in DMEM/F-12 (Gibco) supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin (Invitrogen) at 37 °C in a 5% CO₂ humidified incubator. Bisphenol A (BPA; CAS Number: 80–05-7; purity ≥99%) was acquired from Sigma-Aldrich (St. Louis, MO, United States), and dimethyl sulfoxide (DMSO) was sourced from Sorabi (Beijing, China). For experimental use, BPA was first dissolved in DMSO to prepare a stock solution, which was then further diluted with cell culture medium to achieve the target working concentrations.

Logarithmic-phase KGN cells were seeded into 96-well plates, then treated with BPA at gradient concentrations (0, 10, 50, 100, 150, 200, 250, 300, or 350 μM) for exposure durations of 0, 24 h, 48 h, 72 h, and 96 h, respectively. Post-treatment, 10 μL of CCK-8 reagent (G-Clone, Beijing) was added to each well, followed by a 4-h incubation. The absorbance at 450 nm was measured using a microplate reader, and this reading was used to assess cell viability. We acknowledge that the micromolar concentrations used in vitro exceed typical biomonitoring-derived internal exposure levels. In the present study, this concentration gradient was applied as a toxicological/mechanistic in vitro dosing range to establish a dose-time response and estimate an IC50 within a short experimental window, which is a common approach in granulosa cell/KGN models when evaluating BPA-related cellular stress, viability and apoptosis phenotypes (Celar Sturm et al., 2025; Tang et al., 2024; Huang et al., 2024; Li et al., 2024). Importantly, this design is not intended to equate acute micromolar exposure with chronic low-dose endocrine disruption; low-dose and long-term exposure paradigms may elicit distinct biological responses and warrant dedicated experimental models (Liu et al., 2021; Rajkumar et al., 2021).

KGN cells (1 × 10⁵/well) were seeded in 6-well plates, treated with BPA(0, 10, 50, or 100 μM, for 24 h), rinsed with PBS, and digested with 0.25% trypsin (Gibco). The concentrations were selected to cover a lower-to-higher toxicological range based on the CCK-8 dose–response and prior literature. Cells were centrifuged at 1000 rpm for 5 min, resuspended in 500 µL 1× Binding Buffer, and stained with 5 µL Annexin V-APC and 10 µL PI (Elabsciences, Wuhan) for 5 min at room temperature in the dark. Apoptotic cells were quantified using a BD Accuri™ C6 flow cytometer.

Total RNA was extracted from KGN cells and primary granulosa cells using TRIzol (Invitrogen), with purity and concentration determined by the A260/A280 ratio. cDNA was synthesized via a Takara Reverse Transcription Kit. qPCR was performed with SYBR GREEN Master Mix (Takara) on an ABI 7500 system. Relative gene expression was calculated using the 2^-△△Ct method, normalized to GAPDH. The primer sequences are shown in Table 1.

KGN cells or primary granulosa cells were lysed in RIPA buffer on ice for 10 min, and supernatants were collected after centrifugation. Protein concentration was measured via BCA assay. Equal amounts of protein were separated by 10% SDS-PAGE, transferred to PVDF membranes (Millipore), blocked with 5% non-fat milk for 1 h at room temperature, incubated with primary antibodies overnight at 4 °C, and then with secondary antibodies for 1 h at room temperature. Bands were visualized using ECL reagent.

Statistical analyses were performed using R (v4.3.3), GraphPad Prism (v9.0), and IBM SPSS Statistics (v26.0). Normality and homogeneity of variance were assessed using Shapiro–Wilk and Levene’s tests, respectively. For two-group comparisons, an unpaired two-tailed Student’s t-test was applied to normally distributed data; otherwise, the Mann–Whitney U test was used. For comparisons involving ≥3 groups, one-way ANOVA followed by Dunnett’s multiple-comparisons test versus the control group was used; when distributional assumptions were not met, the Kruskal–Wallis test with Dunn’s post hoc correction was applied. Time-course proliferation data were analyzed by two-way ANOVA (time × treatment) with Dunnett’s multiple-comparisons test. Dose–response curves and IC50 values were estimated by nonlinear regression using a four-parameter logistic model. Correlations were evaluated using Spearman’s rank correlation. For bulk transcriptome-derived ssGSEA scores and hub-gene expression comparisons, two-sided Wilcoxon rank-sum tests were used, with Benjamini–Hochberg correction applied where multiple immune cell types were tested simultaneously. ROC AUCs and 95% confidence intervals were computed using the DeLong method. All tests were two-sided, and p < 0.05 was considered statistically significant.

We screened 139 shared genetic markers between BPA and PCOS (Figure 1A). PPI network visualization of molecular interactions among these genes revealed an intricate regulatory network reflecting key biological pathways in BPA - induced PCOS pathogenesis. Nodes were color-coded by connectivity degree, with hub genes exhibiting high centrality, suggesting their pivotal roles in network architecture. Functional enrichment analysis (Figures 1B–E) uncovered notable enrichment in hormonal and reproductive pathways. Significantly enriched biological processes included steroid hormone biosynthesis, metabolic processes, secretion, translocation, and regulatory mechanisms. Predominant regulatory pathways such as lipid metabolism, cell growth and death, endocrine metabolic disease, signaling molecules interaction, and endocrine system regulation were also prominent, underscoring their relevance to PCOS pathophysiology. Remarkably, the cell growth and apoptosis signaling pathway emerged as one of the most significantly enriched, highlighting that apoptosis modulation serves as a pivotal mechanism through which BPA drives PCOS pathogenesis. These findings validate that BPA exposure potentially disrupts endocrine homeostasis and reproductive signaling, aligning with the hypothesis that BPA contribute to PCOS progression by interfering with hormonal regulatory networks.

To identify PCOS-related genetic markers, we first integrated three GEO datasets into a unified training set and evaluated batch-effect correction by comparing gene expression distributions and principal component analysis before and after integration (Figures 2A–D). Differential expression analysis of the merged dataset identified 448 differentially expressed genes (DEGs) between PCOS and control samples (Figure 2E). Given that pathway analyses highlighted apoptosis as a key process associated with BPA exposure, we focused on apoptosis-related mechanisms for subsequent analyses. We then intersected the 448 DEGs with 139 BPA-PCOS hub genes retrieved from the CTD and with two phenotype-related gene sets (oxidative stress and apoptosis), yielding five overlapping hub genes: PTAFR, RACGAP1, CYP19A1, FSHR, and DMD (Figure 2F). These genes may represent key molecular links between BPA-associated perturbations and PCOS pathophysiology. A heatmap further illustrated the expression patterns of these five hub genes in PCOS and control samples within the combined dataset (Figure 2G). Boxplot analysis revealed significant between-group differences in hub gene expression (Figure 2H): PTAFR was significantly upregulated in PCOS samples, whereas RACGAP1, CYP19A1, FSHR, and DMD were downregulated. These expression patterns were further supported in an independent validation dataset (Supplementary Figure S2D). Collectively, these findings identify five hub genes as candidate mediators linking BPA-associated molecular perturbations to PCOS-related granulosa cell dysfunction.

Cox regression analyses explored associations between five hub genes (PTAFR, RACGAP1, CYP19A1, FSHR, DMD) and PCOS risk. Results (Figure 3A) showed PTAFR linked to increased risk, while RACGAP1, FSHR, and DMD were associated with reduced risk. Multivariate Cox regression (Figure 3B) confirmed PTAFR and FSHR as independent predictors in risk stratification. A nomogram integrating the five hub genes was developed to predict PCOS risk in BPA-exposed individuals (Figure 3C). ROC curve analysis (Figure 3D) demonstrated strong discriminative performance with an AUC of 0.878 in the primary set and 0.802 in the validation set (Supplementary Figure S2C), confirming its ability to differentiate cases. The calibration curve (Figure 3E) showed good alignment between predicted and observed outcomes, and decision curve analysis (Figure 3F) validated clinical utility via consistent net benefit across threshold probabilities. These findings support the nomogram as a potential screening tool for BPA-mediated PCOS risk, aiding early intervention. ROC curves for the training set (Figure 3G) assessed the five hub genes’ risk-discriminating ability, and a chromosomal localization map (Figure 3H) illustrated their genomic positions, contextualizing their roles in PCOS pathogenesis.

To elucidate the potential pathobiological roles of five hub genes in PCOS development, GSEA was conducted (Supplementary Figure S3). PTAFR, enriched in immune cell activation, cytokine signaling, and inflammatory pathways (processes tightly intertwined with apoptosis), may mediate follicular inflammation and hormonal disruption in PCOS by promoting granulosa cell apoptosis via immune dysregulation-driven inflammatory cascades (Supplementary Figure S3A). RACGAP1’s enrichment in cell cycle/proliferation pathways highlights its critical role in regulating cell cycle progression, whose dysregulation disrupts granulosa cell proliferation-apoptosis balance (Supplementary Figure S3B). CYP19A1, a core regulator of steroid hormone biosynthesis and ovarian function linked to ovarian endocrine balance pathways, may indirectly mediate granulosa cell apoptosis via endocrine disruption (Supplementary Figure S3C). FSHR, essential for follicle development and gonadotropin signaling, enriched in ovarian function pathways, may contribute to PCOS anovulation and hormonal dysfunction via impaired FSH-FSHR interactions that disrupt folliculogenesis and trigger granulosa cell apoptosis (Supplementary Figure S3D). DMD—enriched in cancer-related pathways (cell cycle dysregulation, aberrant growth)—may perturb PCOS ovarian cell dynamics by dysregulating granulosa cell proliferation/survival and shifting cell fate toward apoptosis (Supplementary Figure S3E). Collectively, five hub genes converge on granulosa cell apoptosis as a central PCOS pathogenic node, mediating apoptosis via interconnected immune regulation, cell cycle control, and ovarian endocrine signaling (Supplementary Figure S3F) and providing a basis for functional studies on their roles in BPA-induced granulosa cell apoptosis in PCOS.

Immune cell infiltration patterns and their associations with five hub genes were analyzed to explore immune-related alterations in PCOS (Figure 4). A boxplot (Figure 4A) showed distinct immune profiles between PCOS and control cohorts, with increased estimated infiltration of activated effector memory CD8 T cells, central memory CD4 T cells, immature dendritic cells, B cells, natural killer T cells, neutrophils, and regulatory T cells in PCOS, indicating altered immune signatures in PCOS. Correlation analysis (Figure 4B) revealed significant associations between hub gene expression and predicted immune cell infiltration. Specifically, PTAFR expression was positively correlated with effector memory CD8 T cells, activated CD4 T cells, and regulatory T cells, suggesting a potential link with T-cell activation and immune regulation. RACGAP1 expression showed negative associations with central memory activated B cells and central memory CD4/CD8 T cells, indicating a possible relationship with B-cell activation and T-cell memory status. CYP19A1 expression was negatively correlated with eosinophils and natural killer T cells, suggesting an association with inflammatory and cytotoxic immune components. DMD expression was negatively correlated with activated dendritic cells, implying a potential relationship with antigen-presenting cell activity. In addition, FSHR expression was positively associated with T follicular helper cells, which are involved in follicular immune regulation (Supplementary Figure S4C3). Targeted correlation plots (Figures 4C–F) further illustrated key gene–immune cell associations, highlighting coordinated patterns between hub genes and immune cell signatures in PCOS. Because immune cell infiltration was inferred computationally from bulk transcriptomic data, these results represent association-based, hypothesis-generating evidence rather than experimentally validated immune-mediated mechanisms.

To characterize the crosstalk between five hub genes (PTAFR, RACGAP1, CYP19A1, FSHR, DMD) and BPA compounds, molecular docking analysis was performed. All five proteins demonstrated robust binding to BPA, with docking scores below −6 kcal mol^-1 (Table 2). Notably, PTAFR exhibited the highest binding affinity to BPA, followed by FSHR and CYP19A1. These findings indicate that BPA may exert the strongest regulatory impact on these hub genes, thereby substantially influencing the pathobiological disruptions characteristic of PCOS. The interaction interfaces between BPA and the five hub proteins were visualized using PyMOL (Figure 5), yielding structural insights at the molecular level into BPA-driven regulation of PCOS-associated gene functions. The results emphasize the potential for BPA exposure to directly influence the activity of five hub genes through specific binding interactions, shedding light on their roles in PCOS pathophysiology.

To systematically characterize the ovarian cellular landscape in a mouse PCOS-like model and to localize hub genes dysregulation at single-cell resolution, we analyzed the scRNA-seq dataset GSE268919 (mouse; DHEA-induced PCOS-like model; detailed quality control in Supplementary Figures S5A–D). UMAP clustering identified 22 cell clusters annotated into 9 major types (Supplementary Figures S5E, F), with heterogeneous cell composition between PCOS and control groups (Supplementary Figure S5G, H). Feature plots (Figures 6A–D) revealed cell-type-specific expression: Cyp19a1 was selectively enriched in granulosa and theca cells, while Dmd was specific to germ cells (with reduced expression in PCOS). Dot plots (Figures 6E,F) further quantified this dysregulation: Cyp19a1 expression was significantly elevated in PCOS granulosa cells, whereas Dmd was downregulated in PCOS germ cells. Phenotype analysis underscored the critical connection between these hub genes and granulosa cell apoptosis. Supplementary Figures S7A, B showed distinct clustering of high-stress/apoptosis cells in PCOS, and violin plots (Supplementary Figure S7C, D) confirmed differential apoptotic signature levels between groups. Notably, the high oxidative stress/apoptosis cell populations strongly overlapped with Cyp19a1-expressing granulosa cells—a finding that directly links Cyp19a1 dysregulation to redox imbalance-driven granulosa cell apoptosis. Additionally, elevated AUcell scores for stress/apoptosis phenotypes in PCOS immune and interstitial cells coincided with disrupted Cyp19a1 expression (Supplementary Figure S7A–F), suggesting a synergistic effect of hub gene dysregulation and microenvironmental stress on promoting granulosa cell apoptosis. Collectively, the scRNA-seq analysis demonstrates that hub genes (Cyp19a1 and Dmd) exert cell-type-specific regulatory effects in PCOS, with Cyp19a1 dysregulation in granulosa cells serving as a key molecular link to oxidative stress and apoptotic pathways—a core pathogenic mechanism driving granulosa cell dysfunction and PCOS progression. Dmd downregulation in germ cells, while indirectly related, further supports the broader role of hub gene networks in mediating cellular apoptosis and reproductive impairment in PCOS. Because these scRNA-seq data are derived from a mouse model, we interpret the findings as cell-type-resolved, hypothesis-generating evidence; translational relevance is primarily supported by independent human bulk transcriptomic datasets and validation in human primary granulosa cells.

To validate the biological effects of BPA on ovarian granulosa cells, we first compared hub genes expression in primary granulosa cells from PCOS patients and healthy controls. Results showed significant discrepancies in mRNA and protein levels of PTAFR, RACGAP1, CYP19A1, FSHR, and DMD between the two groups (Figures 7A,B). Functional assays on KGN cells revealed dose- and time-dependent effects of BPA. Time-course proliferation assays demonstrated suppressed cell growth with increasing BPA concentrations (Figure 7C), while CCK-8 viability analysis yielded an IC50 value of 129.3 μM (Figure 7D), indicating reduced viability under toxicological in vitro BPA exposure. Flow cytometry showed a significant dose-dependent increase in the apoptosis rate of BPA-treated KGN cells (Figure 7E), with marked elevation at 50 and 100 μM. Molecular analyses of BPA-exposed KGN cells showed dysregulated expression of hub genes: mRNA and protein levels of hub genes (PTAFR, RACGAP1, CYP19A1, FSHR, DMD) were significantly altered in a dose-dependent manner (Figures 7F–H). Additionally, BPA treatment modulated the expression of apoptosis-related markers: BAX and Caspase3 mRNA and protein levels were upregulated, while BCL-2 was downregulated, consistent with the increased apoptotic rate (Figures 7I,J). Collectively, these data demonstrate that toxicological in vitro BPA exposure reduces viability, inhibits proliferation, and promotes apoptosis in ovarian granulosa cells, accompanied by dysregulated expression of PCOS-related hub genes and apoptosis signaling pathway components.