Integrated computational and in vitro analysis reveals the toxic mechanisms of PFAS in female diminished ovarian reserve.

Information on the four PFAS compounds, including their CAS numbers and SMILES, was obtained from PubChem, as shown in Table S2 . The LD50 values of the PFAS compounds were predicted using ProTox 3.0. Specifically, the LD50 value for both PFOS and PFHxS was 154 mg/kg, whereas those for PFOA and PFNA were 518 mg/kg. As shown in Fig.  1 A and B, the results indicate that PFAS exhibited genotoxicity, eye corrosion/irritation, respiratory toxicity, skin sensitization, repeated dose toxicity, and nephrotoxicity. Toxicity predictions from admetSAR 3.0 suggest that all PFAS compounds demonstrate reproductive toxicity. Furthermore, they present a high risk of mitochondrial toxicity (48.5%-62.8%), which may impair ovarian energy metabolism and subsequently lead to diminished ovarian reserve. Although ADMETlab 3.0 does not explicitly mention the reproductive toxicity of PFAS, its prediction results indicate a high risk of genotoxicity (94.6%-98.1%), suggesting potential interference with the normal transmission of cellular genetic information. Interestingly, Vnn-ADMET predictions revealed that both PFOS and PFHxS may cross the blood-brain barrier and inhibit matrix metalloproteinase activity. These preliminary toxicity assessments established a foundation for exploring the potential toxic mechanisms of PFAS. Fig. 1 ( A ) Heatmap of the toxicity of PFAS. ( B ) Radar chart of the toxicity of PFAS. ( C ) Venn diagram of PFAS targets. ( D ) Upset plot of PFAS targets. ( E ) Venn diagram of the potential targets of PFAS and DOR. ( F ) The PPI network of potential targets. ( G–I ) GO enrichment analysis. ( G ) Circular chord diagram. Chords connect candidate genes to significantly enriched BP. ( H ) Circular chord diagram. Chords connect candidate genes to significantly enriched MF. ( I ) Circular chord diagram. Chords connect candidate genes to significantly enriched CC. ( J ) Circular chord diagram. Chords connect candidate genes to significantly enriched KEGG pathways. ( K ) Bubble plot of the GO analysis. ( L ) Sankey-bubble plot of 10 KEGG enrichment pathways and candidate genes. ( A ) Heatmap of the toxicity of PFAS. ( B ) Radar chart of the toxicity of PFAS. ( C ) Venn diagram of PFAS targets. ( D ) Upset plot of PFAS targets. ( E ) Venn diagram of the potential targets of PFAS and DOR. ( F ) The PPI network of potential targets. ( G–I ) GO enrichment analysis. ( G ) Circular chord diagram. Chords connect candidate genes to significantly enriched BP. ( H ) Circular chord diagram. Chords connect candidate genes to significantly enriched MF. ( I ) Circular chord diagram. Chords connect candidate genes to significantly enriched CC. ( J ) Circular chord diagram. Chords connect candidate genes to significantly enriched KEGG pathways. ( K ) Bubble plot of the GO analysis. ( L ) Sankey-bubble plot of 10 KEGG enrichment pathways and candidate genes. We integrated target information from the ChEMBL, STITCH, SwissTarget Prediction, and PharmMapper databases and identified 499 targets for PFOS, 468 targets for PFOA, 468 targets for PFHxS, and 510 targets for PFNA. Among these, PFOS and PFNA shared 486 common targets, whereas all four compounds shared 435 common targets, designated as PFAS toxicity targets, as shown in Fig.  1 C and D. Additionally, a total of 932 DOR-related targets were mined and consolidated from the GeneCards and OMIM databases. By comparing the previously identified 435 PFAS toxicity targets with the DOR-specific disease targets, Venn diagram analysis revealed 35 overlapping targets as potential toxicity targets for PFAS-induced DOR, as illustrated in Fig.  1 E. A PPI network was constructed for the overlapping targets of PFAS and DOR using the STRING database and visualized using Cytoscape. Analysis of the 35 potential targets associated with both DOR and PFAS via the STRING database generated a PPI network comprising 33 nodes and 140 edges, with an average node degree of eight, as shown in Fig.  1 F. Subsequently, CytoNCA and CytoHubba plugins were employed to identify the top 10 core targets based on degree, including albumin (ALB), estrogen receptor 1 (ESR1), Peroxisome Proliferator activated receptor gamma (PPARG), Erb-b2 receptor tyrosine kinase 2 (ERBB2), Heat Shock Protein 90 alpha family class A member 1 (HSP90AA1), insulin-like growth factor 1 receptor (IGF1R), glycogen synthase kinase 3 beta (GSK3B), androgen receptor (AR), Annexin A5 (ANXA5), and mitogen-activated protein kinase 8 (MAPK8), as detailed in Table S3 and Fig.  1 F. Additionally, the MCODE plugin was used to screen for core target clusters and a cluster diagram was generated, as shown in Fig.  1 F. GO and KEGG pathway analyses were performed for 35 targets of PFAS-induced DOR using Metascape. The analysis identified 482 BP, 25 CC, and 81 MF in the GO entries. Figure  1 G–I presents the top 10 entries ranked by the − lg P values for BP, CC, and MF. BP related to PFAS-induced DOR is primarily associated with the regulation of cell proliferation, apoptosis, and gland development, among others. CC analysis revealed significant enrichment of genes involved in vesicle transport and secretion as well as membrane dynamics and assembly. MF analysis indicated strong associations with small-molecule binding, kinase binding and activity, nuclear receptor activity, and related functions. In addition, KEGG pathway enrichment analysis identified 72 significant pathways. Figure  1 K and L display the top 10 entries of − lg P values for KEGG. Key targets were predominantly enriched in pathways regulating cell proliferation and differentiation (such as PI3K-Akt, MAPK, Ras, and ErbB signaling), immune and inflammatory responses (e.g., Th17 cell differentiation), and hormone signaling pathways (including prolactin, estrogen, and progesterone signaling). These results collectively implicate a cellular signaling network that is dependent on growth factor receptors and sex hormones. We obtained the GSE232306 dataset using the GEOquery package and normalized the Fragments Per Kilobase of exon model per million mapped fragment (FPKM) values using the limma package. Differential genes were screened using the thresholds of |log₂FC| > 1 and adjusted P -value ( P .adjust) < 0.05. A total of 5836 differentially expressed genes (DEGs) were identified, including 2593 upregulated and 3243 downregulated genes. The intersection between all DEGs and the 35 target genes yielded 19 overlapping genes. Differential analysis was performed using the limma package, while visualization was performed using the ggplot2 package for the volcano plot, the pheatmap package for the heatmap, and the ggvenn package for the Venn diagram (Fig.  2 A–C). Fig. 2 ( A ) Volcano plot of DEGs in the dataset (GEO accession: GSE232306 ). Purple dots represent 2593 upregulated genes, green dots indicate 3243 downregulated genes. ( B ) Venn diagram of the targets of PFAS and DEGs. ( C ) Hierarchical clustering heatmap of overlapping genes (rows: genes, columns: samples). Purple/green color intensity correlates with up/downregulation. ( D–G ) The results of Machine Learning. ( D ) Results from the RF classifier. ( E – F ) Feature selection via LASSO regression. ( G ) Results from the XGBoost algorithm. ( H ) Boxplots displays the differential expression of 3 core genes (GSTP1, PRDX5 and TGFBR1) in normal and DOR samples. ( I ) Core genes (GSTP1, PRDX5 and TGFBR1) localization in chromosomes. ( J ) ROC of 3 core genes demonstrates their ability to distinguish between DOR samples and normal samples. ( K ) External validation of the core targets in the independent GSE315266 dataset. ( L ) External validation of the core targets in the independent GSE296744 dataset. ( M ) Sample gene GSEA enrichment plots. * P  < 0.05, ** P  < 0.01 vs. control. ( A ) Volcano plot of DEGs in the dataset (GEO accession: GSE232306 ). Purple dots represent 2593 upregulated genes, green dots indicate 3243 downregulated genes. ( B ) Venn diagram of the targets of PFAS and DEGs. ( C ) Hierarchical clustering heatmap of overlapping genes (rows: genes, columns: samples). Purple/green color intensity correlates with up/downregulation. ( D–G ) The results of Machine Learning. ( D ) Results from the RF classifier. ( E – F ) Feature selection via LASSO regression. ( G ) Results from the XGBoost algorithm. ( H ) Boxplots displays the differential expression of 3 core genes (GSTP1, PRDX5 and TGFBR1) in normal and DOR samples. ( I ) Core genes (GSTP1, PRDX5 and TGFBR1) localization in chromosomes. ( J ) ROC of 3 core genes demonstrates their ability to distinguish between DOR samples and normal samples. ( K ) External validation of the core targets in the independent GSE315266 dataset. ( L ) External validation of the core targets in the independent GSE296744 dataset. ( M ) Sample gene GSEA enrichment plots. * P  < 0.05, ** P  < 0.01 vs. control. To minimize the risk of overfitting inherent in small-sample transcriptomic datasets, we employed a consensus-based machine learning pipeline integrating RF, LASSO (with L1-regularization), and XGBoost. Using the random forest package, the FPKM profiles of the intersecting genes were analyzed between the DOR and control groups. With the seed set to 2024 and ntree = 1000, 11 genes were identified based on a MeanDecreaseGini threshold greater than 0.2 (Fig.  2 D). Subsequently, the glmnet package was applied to the data matrix of these 11 RF-selected genes, with parameters set.seed(500) and family = “binomial.” The optimal model corresponding to lambda.min was selected, which retained 4 genes (Fig.  2 E and F). Finally, the xgboost package was used with a maximum tree depth of 3, a subsampling ratio of 0.9, and regularization to prevent overfitting. A 5-fold cross-validation was performed to select the optimal model. Variable importance within the model was calculated, and genes with non-zero importance were selected, yielding three key genes: GSTP1, PRDX5, and TGFBR1(Fig.  2 G) 70 . The ggpubr package was used to construct box plots to visualize the expression levels of the three key genes. Compared with the control group, the expression of GSTP1 and PRDX5 was significantly downregulated, while that of TGFBR1 was significantly upregulated in the granulosa cells of DOR patients (Fig.  2 H). Chromosomal location analysis indicated that GSTP1 and PRDX5 were located on chromosome 11 and TGFBR1 was located on chromosome 9 (Fig.  2 I). As shown in Fig.  2 J, ROC curves were generated using the pROC package and the area under the curve (AUC) values were calculated. All three key genes demonstrated a high diagnostic accuracy (AUC > 0.9). To further validate the robustness of the identified core targets, external ROC analysis was performed using two independent datasets. In GSE315266 , GSTP1, PRDX5, and TGFBR1 achieved AUC values of 0.75, 0.78, and 0.75, respectively (Fig. 2  K). In GSE296744 , the diagnostic efficacy remained consistent, with AUC values of 0.639 (GSTP1), 0.667 (PRDX5), and a notably high 0.819 (TGFBR1) (Fig. 2  L). These results demonstrate the cross-cohort stability and diagnostic potential of the selected core genes. To further investigate the potential regulatory mechanisms of the core genes in DOR, we performed a single-gene GSEA for each of the three core genes using the GSE232306 dataset. Our analysis revealed that the expression of these three core genes was closely associated with signal transduction and metabolic pathways, including the ErbB signaling pathway, insulin signaling pathway, phosphatidylinositol signaling system, and pyrimidine metabolism (Fig.  2 M). These findings further support the notion that DOR progression is a complex biological process and that the three core genes may influence DOR development by regulating distinct pathways. Molecular docking analysis was performed to investigate the interactions between PFAS and three key genes: GSTP1 (PDB ID: 5J41), PRDX5 (PDB ID: 3MNG), and TGFBR1 (PDB ID: 5E8S). Binding energy was used to evaluate the interactions between the ligands and receptors: a binding energy < 0 kcal/mol indicated spontaneous binding, while a binding energy < -5 kcal/mol suggested stable binding. The molecular docking results demonstrated that the binding energies between the three key genes and the four PFAS compounds ranged from − 9.8 to -6.5 kcal/mol, with a mean value of -7.9 kcal/mol, indicating strong binding activity between the three key genes and all four PFAS compounds (Fig.  3 A–L). Among these, the most stable interactions were observed between TGFBR1 and the four PFAS compounds, all with binding energies lower than − 8.0 kcal/mol (Fig.  3 M). These robust docking results highlight the spontaneous and strong binding capability of PFAS with these key genes, underscoring their high affinity and critical role in the molecular mechanisms affecting the ovarian reserve. To further elucidate the complex binding configurations, we performed docking visualizations to display the lowest-energy binding conformations. This visualization vividly illustrates the close interactions between PFAS and their key targets, providing valuable structural insights into their binding mechanisms and reinforcing the significance of PFAS in influencing the ovarian reserve at the molecular level. Fig. 3 ( A–L ) 3D molecular docking results of PFAS with core targets. ( A ) GSTP1-PFOS ( B ) PRDX5-PFOS ( C ) TGFBR1-PFOS ( D ) GSTP1-PFOA ( E ) PRDX5-PFOA ( F ) TGFBR1-PFOA ( G ) GSTP1-PFHxS ( H ) PRDX5-PFHxS ( I ) TGFBR1-PFHxS ( J ) GSTP1-PFNA ( K ) PRDX5-PFNA ( L ) TGFBR1-PFNA. Purple and red denotes target proteins, while green represents PFAS molecules. ( M ) Heatmap visualization of binding energy (kcal/mol) derived from molecular docking involving 3 core targets and their corresponding PFAS. ( N–R ) Molecular dynamics analysis of the core target-PFOA complex. ( N ) RMSD curve of GSTP1-PFOA complex (green line), PRDX5-PFOA complex (red line) and TGFBR1-PFOA complex (purple line). ( O ) Rg curve of the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( P ) SASA values for the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( Q ) Number of hydrogen bonds in the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( R ) RMSF curve of GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( A–L ) 3D molecular docking results of PFAS with core targets. ( A ) GSTP1-PFOS ( B ) PRDX5-PFOS ( C ) TGFBR1-PFOS ( D ) GSTP1-PFOA ( E ) PRDX5-PFOA ( F ) TGFBR1-PFOA ( G ) GSTP1-PFHxS ( H ) PRDX5-PFHxS ( I ) TGFBR1-PFHxS ( J ) GSTP1-PFNA ( K ) PRDX5-PFNA ( L ) TGFBR1-PFNA. Purple and red denotes target proteins, while green represents PFAS molecules. ( M ) Heatmap visualization of binding energy (kcal/mol) derived from molecular docking involving 3 core targets and their corresponding PFAS. ( N–R ) Molecular dynamics analysis of the core target-PFOA complex. ( N ) RMSD curve of GSTP1-PFOA complex (green line), PRDX5-PFOA complex (red line) and TGFBR1-PFOA complex (purple line). ( O ) Rg curve of the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( P ) SASA values for the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( Q ) Number of hydrogen bonds in the GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. ( R ) RMSF curve of GSTP1-PFOA complex, PRDX5-PFOA complex and TGFBR1-PFOA complex. The root mean square deviation (RMSD) serves as a reliable metric for evaluating the conformational stability of proteins and ligands, reflecting the degree of atomic positional deviation from their initial states. Lower deviations indicate better conformational stability. Therefore, the equilibrium of the simulation systems was assessed using the RMSD. As shown in Fig.  3 N, the GSTP1–PFOA and TGFBR1–PFOA complex systems reached equilibrium after 10 ns, fluctuating around 2 Å and 1.8 Å, respectively. The PRDX5–PFOA complex system reached equilibrium after 65 ns, ultimately fluctuating at approximately 2.2 Å. These results indicated that PFOA exhibits high stability when bound to GSTP1, TGFBR1, and PRDX5. Further analysis revealed that the radius of gyration (Rg) and solvent-accessible surface area (SASA) of GSTP1–PFOA, PRDX5–PFOA, and TGFBR1–PFOA complexes remained stable throughout the simulation, showing consistent fluctuations (Fig.  3 O and P). This suggests that the PFOA-bound complexes with GSTP1, TGFBR1, and PRDX5 did not undergo significant contraction or expansion during simulation. Hydrogen bonds play a crucial role in ligand-protein binding. The number of hydrogen bonds between the small molecules and the target proteins during the dynamic process is shown in Fig.  3 Q. For the GSTP1–PFOA complex, the number of hydrogen bonds ranged from 0 to 3, with approximately one hydrogen bond maintained under most conditions. The PRDX5–PFOA complex system exhibited 0–5 hydrogen bonds, with approximately two hydrogen bonds formed in most cases. In the TGFBR1–PFOA complex system, the number of hydrogen bonds varied from zero to six, and approximately three hydrogen bonds were observed under most conditions. These results indicate favorable hydrogen bond interactions in small-molecule–target protein complexes. Root mean square fluctuation (RMSF) reflects the flexibility of amino acid residues in proteins. As shown in Fig.  3 R, GSTP1–PFOA, PRDX5–PFOA, and TGFBR1–PFOA complex systems exhibited relatively low RMSF values (mostly below 3 Å), indicating lower flexibility and higher structural stability. In summary, GSTP1–PFOA, PRDX5–PFOA, and TGFBR1–PFOA complex systems exhibited stable binding and favorable hydrogen bond interactions. These findings suggest that PFOA can stably interact with GSTP1, TGFBR1, and PRDX5 target proteins at the molecular level, thus laying a computational foundation for further exploration of its potential role in the DOR. KGN cells were treated with increasing concentrations of PFOA (0, 0.2, 2, 20, and 200 µM) for 24 and 48 h, and cell viability was assessed using CCK-8 assay. PFOA significantly reduced cell viability in a dose- and time-dependent manner, with 200 µM PFOA treatment for 24 h markedly inhibiting cell proliferation (Fig.  4 A). Similarly, exposure of KGN cells to 0.2–200 µM PFOA for 48 h also suppressed proliferative activity (Fig.  4 B). To evaluate the chronic effects of PFOA at environmentally relevant levels, KGN cells were exposed to 10, 100, and 1000 nM PFOA for 72 h. The CCK-8 assay revealed a subtle but statistically significant dose-dependent suppression of cell proliferation, with cell viability remaining high at 94.2% and 92.1% for the 100 nM and 1000 nM groups, respectively (Fig.  4 C). To assess the effect of PFOA on KGN proliferation, an EdU assay was performed. As shown in Fig.  4 D and F, treatment with 20 µM and 200 µM PFOA for 24 h significantly decreased the total number of EdU-positive cells compared to that in the control group. Similar results were observed after 48-hour treatment (Fig.  4 E and G). These findings indicate that PFOA inhibits DNA synthesis and proliferation in KGN cells. Given the pronounced suppression of cell proliferation and viability by 200 µM PFOA, concentrations ranging from 0 to 20 µM were selected for subsequent mechanistic studies to avoid potential confounding effects due to the reduced cell numbers. Fig. 4 PFOA reduced cell viability, induced cytotoxicity, and promoted apoptosis in KGN cells. ( A–F ) PFOA reduced KGN cell viability in a dose- and time-dependent manner. ( A – B ) KGN cells were treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( A ) and 48 h ( B ), and cell viability was assessed using the CCK-8 assay. The control group (0 µM) was treated with an equivalent volume of DMSO (0.1%). Data are expressed as mean ± SD ( n  = 4). ( C ) KGN cells were treated with PFOA (0, 10, 100 and 1000 nM) for 72 h, and cell viability was assessed using the CCK-8 assay. ( D – E ) The EdU positive cell percent (%) were calculated using ImageJ. KGN cells were treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( D ) and 48 h ( E ). ( F – G ) EdU incorporation assay showing reduced proliferation in KGN cells treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( F ) and 48 h ( G ). Hoechst 33,342 was used for nuclear staining (blue), and EdU positive proliferating cells are labeled in red. Data are expressed as mean ± SD ( n  = 3). Scale bars = 200 μm. ( H–K ) PFOA reduced KGN cell apoptosis in a dose- and time-dependent manner. Apoptosis rates in KGN cells treated with PFOA (0, 0.2, 2 and 20 µM) for 24 h ( H ) and 48 h ( J ) were determined by Annexin V-FITC/PI, and total apoptoisis rate (%) were calculated ( I and K ). The results were evaluated statistically. ns, not significant; * P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** P  < 0.0001 vs. control. PFOA reduced cell viability, induced cytotoxicity, and promoted apoptosis in KGN cells. ( A–F ) PFOA reduced KGN cell viability in a dose- and time-dependent manner. ( A – B ) KGN cells were treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( A ) and 48 h ( B ), and cell viability was assessed using the CCK-8 assay. The control group (0 µM) was treated with an equivalent volume of DMSO (0.1%). Data are expressed as mean ± SD ( n  = 4). ( C ) KGN cells were treated with PFOA (0, 10, 100 and 1000 nM) for 72 h, and cell viability was assessed using the CCK-8 assay. ( D – E ) The EdU positive cell percent (%) were calculated using ImageJ. KGN cells were treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( D ) and 48 h ( E ). ( F – G ) EdU incorporation assay showing reduced proliferation in KGN cells treated with PFOA (0, 0.2, 2 ,20 and 200 µM) for 24 h ( F ) and 48 h ( G ). Hoechst 33,342 was used for nuclear staining (blue), and EdU positive proliferating cells are labeled in red. Data are expressed as mean ± SD ( n  = 3). Scale bars = 200 μm. ( H–K ) PFOA reduced KGN cell apoptosis in a dose- and time-dependent manner. Apoptosis rates in KGN cells treated with PFOA (0, 0.2, 2 and 20 µM) for 24 h ( H ) and 48 h ( J ) were determined by Annexin V-FITC/PI, and total apoptoisis rate (%) were calculated ( I and K ). The results were evaluated statistically. ns, not significant; * P  < 0.05, ** P  < 0.01, *** P  < 0.001, **** P  < 0.0001 vs. control. To assess the effect of PFOA exposure on apoptosis in KGN cells, flow cytometry was performed following Annexin V-FITC/PI staining. As shown in Fig.  4 H and I, compared with the control group (0 µM), exposure to 2–20 µM PFOA for 24 h significantly increased the percentage of early and late apoptotic cells ( P   0.05). Similar results were observed after 48 h of treatment (Fig.  4 J and K). These findings indicated that PFOA induces apoptosis in KGN cells. Both AMH and inhibin B (INHB) are members of the transforming growth factor-β (TGF-β) superfamily and are secreted by ovarian granulosa cells. Previous studies have identified AMH as a biomarker for DOR 77 , and INHB as an important biomarker for assessing ovarian reserve 78 , 79 . Compared to FSH, estradiol (E₂), INHB, and AFC, AMH reflects the age-related decline in ovarian reserve at an earlier stage, and its levels are not influenced by menstrual cycle, hormonal contraceptives, or pregnancy. Previous research has indicated that INHB is composed of an α-subunit and βB-subunit linked by disulfide bonds 80 . Therefore, these two subunits were validated separately in subsequent mRNA and protein expression experiments. In this experiment, we utilized KGN cells to investigate the effect of PFOA exposure on hormone synthesis and secretion in granulosa cells, thereby evaluating its impact on molecular markers associated with ovarian reserve in granulosa cells.. The results demonstrated that after 48 h of treatment, higher concentrations of PFOA significantly downregulated the mRNA expression levels of AMH, INHBA, and INHBB in KGN cells compared with the control group ( P   0.05) (Fig.  5 A). WB analysis demonstrated that PFOA reduced the expression levels of ovarian reserve-related proteins (INHBA and INHBB) in KGN cells, which was consistent with the RT-qPCR results (Fig.  5 E). ELISA analysis of AMH levels in the cell culture medium revealed that after 48 h of exposure, higher PFOA concentrations progressively reduced AMH secretion ( P  < 0.05) (Fig.  5 C). In contrast, after 24 h of exposure, although a decreasing trend in AMH levels was observed with increasing PFOA concentration, the results did not reach statistical significance (Fig.  5 B). To further evaluate the chronic impact of environmentally relevant PFOA exposure on ovarian granulosa cell function, we measured AMH protein expression levels in KGN cells via ELISA after 72 h of treatment with 0, 10, 100, and 1000 nM PFOA. PFOA exposure at 100 and 1000 nM significantly downregulated AMH expression relative to the control group ( P  < 0.05). In contrast, AMH levels remained unchanged at 10 nM PFOA (Fig.  5 F). These findings suggest that prolonged exposure to high concentrations of PFOA adversely affects the secretion of ovarian-related hormones in granulosa cells. Fig. 5 RT-qPCR analysis of mRNA expression levels of key genes involved in ovarian reserve function biomarker in KGN cells after 24 h ( A ) and 48 h ( D ) of PFOA exposure (0, 0.2, 2 and 20 µM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( B–C ) ELISA analysis of the AMH expression level in KGN cells after 24 h ( B ) and 48 h ( C ) of PFOA exposure (0, 0.2, 2 and 20 µM). Western blot analysis of ovarian reserve function-related proteins in KGN cells following PFOA treatment (0, 0.2, 2 and 20 µM) for 48 h. E Expression of ovarian reserve function-related proteins including INHBA and INHBB. GAPDH was used as a loading control. Band intensities were quantified using ImageJ. Data are shown as mean ± SD ( n  = 3). F ELISA analysis of the AMH expression level in KGN cells after 72 h of PFOA exposure (0, 10, 100 and 1000 nM). * P  < 0.05, ** P  < 0.01, *** P  < 0.01, **** P  < 0.0001 vs. control. RT-qPCR analysis of mRNA expression levels of key genes involved in ovarian reserve function biomarker in KGN cells after 24 h ( A ) and 48 h ( D ) of PFOA exposure (0, 0.2, 2 and 20 µM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( B–C ) ELISA analysis of the AMH expression level in KGN cells after 24 h ( B ) and 48 h ( C ) of PFOA exposure (0, 0.2, 2 and 20 µM). Western blot analysis of ovarian reserve function-related proteins in KGN cells following PFOA treatment (0, 0.2, 2 and 20 µM) for 48 h. E Expression of ovarian reserve function-related proteins including INHBA and INHBB. GAPDH was used as a loading control. Band intensities were quantified using ImageJ. Data are shown as mean ± SD ( n  = 3). F ELISA analysis of the AMH expression level in KGN cells after 72 h of PFOA exposure (0, 10, 100 and 1000 nM). * P  < 0.05, ** P  < 0.01, *** P  < 0.01, **** P  < 0.0001 vs. control. RT-qPCR and Western blot analyses were performed to validate the expression changes of key genes identified in the preliminary network analysis of PFOA-exposed granulosa cells. The results showed that prolonged exposure to higher concentrations of PFOA significantly downregulated the mRNA and protein expression levels of oxidative stress-related genes (GSTP1 and PRDX5) compared to the control group ( P  < 0.05), while the expression level of TGFBR1 was significantly increased ( P  < 0.05) (Fig.  6 A and G). To evaluate the chronic effects of PFOA at environmentally relevant levels, we examined the mRNA expression levels of core target genes (GSTP1, PRDX5, and TGFBR1) in KGN cells following 72 h of exposure to 0, 10, 100, and 1000 nM PFOA using RT-qPCR. At the highest concentration tested (1000 nM), PFOA exposure significantly downregulated the mRNA levels of GSTP1 and PRDX5 ( P  < 0.05) and significantly upregulated TGFBR1 mRNA expression ( P  < 0.05) compared to the control group (Fig.  6 C). These findings are consistent with network toxicology predictions regarding gene expression associated with diminished ovarian reserve. The upregulation of TGFBR1 in KGN cells indicates that PFOA triggers pro-fibrotic signaling at the cellular level, which may contribute to the fibrotic remodeling of the ovarian microenvironment. To evaluate the chronic effects of PFOA at environmentally relevant levels, we measured intracellular ROS levels in KGN cells following 72 h exposure to 0, 10, 100, and 1000 nM PFOA. ROS levels were significantly elevated starting from 100 nM compared to the control group. Notably, while our previous data showed that short-term exposure (24–48 h) to a higher concentration (0.2 µM, equivalent to 200 nM) did not induce downregulation of oxidative stress‑related markers (GSTP1 and PRDX5) at either the mRNA or protein level, prolonged exposure (72 h) to lower concentrations (100 nM and 1000 nM) resulted in a significant increase in ROS production (Fig.  6 D). These findings suggest that chronic PFOA exposure, even at environmentally relevant concentrations, can induce oxidative stress in KGN cells in a time-dependent manner. Treatment with 2–20 µM PFOA for 24 and 48 h significantly up-regulated the mRNA expression level of the pro-apoptotic gene BAX in KGN cells, while downregulating the expression of the anti-apoptotic gene Bcl-2, resulting in an increased BAX/Bcl-2 ratio ( P  < 0.05) (Fig.  6 E and F). Exposure to 0.2 µM PFOA for 24 and 48 h significantly up-regulated BAX mRNA expression ( P   0.05) (Fig.  6 E and F). Furthermore, WB analysis indicated that the expression patterns of these proteins were generally consistent with the RT-qPCR results, suggesting the activation of apoptotic signaling in granulosa cells. Specifically, PFOA treatment led to a significant increase in BAX protein levels and a notable decrease in Bcl-2 protein levels (Fig.  6 H). This provided compelling visual evidence of PFOA’s impact of PFOA on the expression of apoptosis-related proteins, thereby supporting cytotoxicity findings at the cellular level. These results suggest that the cytotoxic effects of PFOA on KGN cells may involve alterations in the pathways related to apoptosis and oxidative stress. Fig. 6 ( A–D ) RT-qPCR analysis of mRNA expression levels of core targets (GSTP1, PRDX5 and TGFBR1) and apoptosis-related genes (BAX and Bcl-2) in KGN cells after 24 h ( A and E ) and 48 h ( B and F ) of PFOA exposure (0, 0.2, 2 and 20 µM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( C ) RT-qPCR analysis of mRNA expression levels of core targets (GSTP1, PRDX5 and TGFBR1) in KGN cells after 72 h of PFOA exposure (0, 10, 100 and 1000 nM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( D ) PFOA induced oxidative stress in KGN cells. KGN cells were treated with PFOA (0, 10, 100 and 1000 nM) for 72 h, and ROS levels in cells were determined by DCFH-DA staining assay. ( G–H ) Western blot analysis of core targets and apoptosis-related proteins in KGN cells following PFOA treatment (0, 0.2, 2 and 20 µM) for 48 h. GAPDH was used as a loading control. Band intensities were quantified using ImageJ. Data are shown as mean ± SD ( n  = 3). * P  < 0.05, ** P  < 0.01, *** P  < 0.01, **** P  < 0.0001 vs. control. ( A–D ) RT-qPCR analysis of mRNA expression levels of core targets (GSTP1, PRDX5 and TGFBR1) and apoptosis-related genes (BAX and Bcl-2) in KGN cells after 24 h ( A and E ) and 48 h ( B and F ) of PFOA exposure (0, 0.2, 2 and 20 µM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( C ) RT-qPCR analysis of mRNA expression levels of core targets (GSTP1, PRDX5 and TGFBR1) in KGN cells after 72 h of PFOA exposure (0, 10, 100 and 1000 nM). Expression levels were normalized to GAPDH and presented as mean ± SD ( n  = 3). ( D ) PFOA induced oxidative stress in KGN cells. KGN cells were treated with PFOA (0, 10, 100 and 1000 nM) for 72 h, and ROS levels in cells were determined by DCFH-DA staining assay. ( G–H ) Western blot analysis of core targets and apoptosis-related proteins in KGN cells following PFOA treatment (0, 0.2, 2 and 20 µM) for 48 h. GAPDH was used as a loading control. Band intensities were quantified using ImageJ. Data are shown as mean ± SD ( n  = 3). * P  < 0.05, ** P  < 0.01, *** P  < 0.01, **** P  < 0.0001 vs. control.

Perfluorooctanoic acid PFOA (171468–5G, purity 95%) was purchased from Sigma-Aldrich (USA). Dulbecco’s modified Eagle’s medium/F-12(DMEM/F-12, CM16405) (1:1) was obtained from M&C Gene Technology (China). Fetal bovine serum (FBS, 04-001-1 A) was purchased from BIOIND (Israel). Penicillin-streptomycin (P/S, C100C5; 1%), protease and phosphatase inhibitors (P002) were obtained from NCM Biotech (China). PowerUp SYBR Green Master Mix (A25741) was purchased from Thermo Fisher Scientific (USA). The RNA-Quick Purification Kit (RN001) was obtained from ES Science (China). Cell counting kit-8 (CCK-8, K1018) was purchased from APExBIO (USA). Trypsin (CSP051, 0.25%) was obtained from ZQXZBIO (China). Bovine serum albumin (BSA, BS114) was purchased from Biosharp (China). The PrimeScript RT Master Mix (RR036A) was obtained from Takara (Japan). EdU Imaging Kit (C0071, 488), dimethyl sulfoxide (DMSO, ST2335, purity ≥ 99.9%), Triton X-100 permeabilization solution (ST1723), RIPA lysis buffer (P0013), BCA protein quantification kit (enhanced, P0009), and BeyoECL Star (ultra-sensitive ECL chemiluminescence detection reagent, P0018AS) were purchased from Beyotime (China). The annexin V-FITC/PI apoptosis kit (A211-01) was obtained from Vazyme (China). Antibodies against BAX (YP-Ab-00317), Bcl-2 (YP-Ab-00322), INHBA (YP-Ab-16029), INHBB (YP-Ab-16028), TGFBR1 (YP-Ab-13695), GSTP1 (YP-Ab-02329), PRDX5 (YP-Ab-06894), and GAPDH (YP-rAb-17721) were purchased from UpingBio (China). The human anti-Müllerian hormone (AMH) ELISA kit (SYP-H0128) was obtained from UpingBio (China). 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE, P0014B) and SDS-PAGE protein loading buffer (P0286) were purchased from Servicebio (China). Polyvinylidene difluoride (PVDF, 36126ES03) membranes were obtained from YEASEN(China). 2’,7’-dichlorofluorescin diacetate (DCFH-DA, HY-D0940) was purchased from MedChemExpress (MCE) (USA). PFOA was dissolved in DMSO to prepare a 200 mM stock solution. Prior to treatment, the stock solution was diluted 1:1000 in the cell culture medium to achieve the maximum working concentration of 200 µM, ensuring that the final DMSO concentration did not exceed 0.1% (v/v). To avoid solvent interference, the control group and all treatment groups were standardized to a final concentration of 0.1% (v/v) DMSO. The human ovarian GC tumor cell line, KGN, was purchased from Procell Life Science and Technology Co. (Wuhan, China). The cells were cultured in DMEM/F-12 (1:1), supplemented with 10% FBS and 1% P/S. Cells were routinely passaged in at 37 °C, 5% CO₂ incubator and passaged with 0.25% trypsin digestion, while the cell density was maintained at 80–90%. The KGN cell line was selected as the in vitro model for this study because it effectively retains the key physiological features of human ovarian granulosa cells, including functional follicle-stimulating hormone (FSH) receptor expression and high aromatase activity. Compared to primary granulosa cells, KGN cells offer superior phenotypic stability and experimental reproducibility, avoiding the rapid luteinization and significant inter-donor variability typically associated with primary cultures. This ensures a consistent biological background for evaluating the toxicological mechanisms of PFAS exposure. The 2D and 3D structures and Simplified Molecular-Input Line-Entry System (SMILES) of perfluorooctanesulfonic acid (PFOS), perfluorooctanoic acid (PFOA), perfluorohexanesulfonic acid (PFHxS), and perfluorononanoic acid (PFNA) were retrieved from PubChem ( https://pubchem.ncbi.nlm.nih.gov/ ). The SMILES of these four PFAS compounds were then input into ProTox 3.0 ( https://tox.charite.de/protox3/ ), ADMETlab 3.0 ( https://admetlab3.scbdd.com/ ), Vnn-ADMET ( https://vnnadmet.bhsai.org/ ), and admetSAR 3.0 ( https://lmmd.ecust.edu.cn/admetsar3/ ) for toxicity prediction and analysis 57 – 60 . The toxic characteristics of the PFASs were predicted by leveraging the structural modeling capabilities of these platforms, thereby providing a preliminary understanding of their toxicological profiles. The toxicity targets of PFAS were retrieved from the ChEMBL database ( https://www.ebi.ac.uk/chembl/ ) using “perfluorooctanoic acid (PFOS)”, “perfluorooctanoic acid (PFOA)”, “perfluorooctanoic acid (PFHXS)”, and “Perfluorononanoic acid (PFNA)” as keywords, with the species filter set to “ Homo sapiens ”. The SMILES structures of these four compounds were then uploaded to the STITCH database ( http://stitch.embl.de/ ) and the SwissTargetPrediction database ( http://www.swisstargetprediction.ch/ ), with the species specified as “ Homo sapiens ” in both platforms 61 – 63 . In the SwissTargetPrediction database, targets with a " probability > 0” were screened as potential toxicity targets. This was followed by a supplementary search of the PharmMapper database ( http://lilab.ecust.edu.cn/pharmmapper ) with the setting “Human Protein Targets Only,” where targets for the four compounds were further filtered using a “Normalized Fit Score > 0.6” as the criterion 64 . The target data from the above four databases were then integrated, consolidated with duplicate entries removed, and standardized using the UniProt database ( https://www.uniprot.org/ ) for target nomenclature normalization 65 , thus establishing a more comprehensive PFAS toxicity target database. The disease-related targets for DOR were gathered as follows: The GeneCards database ( https://www.genecards.org/ ) and the OMIM database ( https://www.omim.org/ ) were searched using “Diminished Ovarian Reserve” as the keyword to identify disease-related targets 66 . The GeneCards database selected targets with “Relevance score ≥ Median” targets and the OMIM database selected targets with *. All retrieved targets were then standardized using the UniProt database, with Homo sapiens specified for both target and gene name conversions. Finally, the targets collected from the two databases were merged and de-duplicated to obtain a consolidated set of DOR-related disease targets. To elucidate the interactions between PFAS toxicity targets and DOR disease targets, an intersection analysis was performed between the established PFAS toxicity target database and DOR disease target database, and a Venn diagram was constructed. The resulting overlapping targets were imported into the STRING database ( https://version-12-0.string-db.org/ ) for the PPI analysis. The analysis was configured with the organism set to Homo sapiens , the minimum required interaction score was set to medium confidence (0.400), and network display options were enabled for 3D bubble design and to hide disconnected nodes in the network. All other parameters were maintained at their default settings 67 . The PPI data in TSV format were then imported into Cytoscape 3.10.1 for further topological analysis. The “CytoNCA” and “CytoHubba” plugins were employed to identify the top 10 targets based on the node degree 53 . Subsequently, the MCODE plugin was used with default parameters (degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max. depth = 100) to identify the cluster with the highest clustering score as the core target cluster and the corresponding network graph was generated. To further investigate the toxicity mechanisms of PFAS-induced DOR, enrichment analysis was performed using Metascape, a well-established platform that integrates resources such as GO and KEGG 68 . GO analysis was conducted to explore the biological processes (BP), cellular components (CC), and molecular functions (MF) underlying the potential toxic mechanisms of PFAS in DOR. KEGG pathway analysis was used to identify the core toxicity pathways associated with the PFAS. The overlapping targets were uploaded to Metascape with the species set to “ Homo sapiens ” and the parameters configured as follows: P Value Cutoff: 0.01, Min Overlap: 3, Min Enrichment: 1.5. Finally, enrichment results were visualized using a bioinformatics platform ( http://www.bioinformatics.com.cn ) 69 . The keyword “Diminished Ovarian Reserve” was entered into the GEO database ( http://www.ncbi.nlm.nih.gov/geo/ ) to search for DOR-related datasets. The retrieved GSE232306 dataset, which included data from six DOR samples and six healthy control samples, was selected. After normalization of the dataset, a bioinformatics platform ( http://www.bioinformatics.com.cn ) was used to perform heatmap and volcanic map analyses of the differentially expressed genes (DEGs). The resulting DEGs were then intersected with PFAS targets and a Venn diagram was drawn to visualize the overlap. Three mainstream machine learning algorithms were employed to screen core targets: random forest (RF), least absolute shrinkage and selection operator (LASSO), and extreme gradient boosting (XGBoost). The RF algorithm was implemented using the R package “randomForest”, LASSO with the “glmnet” package, and XGBoost with the “xgboost” package 70 . We adopted an innovative consolidated analytical framework integrating RF-LASSO-XGBoost to identify the final core genes that are considered potential key targets in PFAS-induced DOR. Subsequently, the bioinformatics platform ( http://www.bioinformatics.com.cn ) was used to generate boxplots to visualize the expression patterns of the core targets, Circos plots to display their chromosomal locations, and receiver operating characteristic (ROC) curves based on the R package “pROC” to validate their diagnostic efficacy. To validate the robustness of the machine learning model, the GSE315266 and GSE296744 dataset were retrieved from the GEO database as external validation sets. Overfitting was controlled by implementing L1-penalty regularization in LASSO, 5-fold cross-validation in XGBoost, and a strict consensus-based selection across the three algorithms. To investigate the potential regulatory mechanisms of the core genes in DOR more comprehensively, we performed single-gene Gene Set Enrichment Analysis (GSEA) on the core genes identified through machine learning. GSEA software (version 4.4.0) was obtained from the GSEA website ( http://software.broadinstitute.org/gsea/index.jsp ) 71 . Based on the expression levels of the machine learning-derived core genes, samples were stratified into high-expression (≥ 50%) and low-expression (< 50%) groups. The c2.cp.kegg.v7.4.symbols.gmt subset was downloaded from the Molecular Signatures Database ( http://www.gsea-msigdb.org/gsea/downloads.jsp ) to evaluate the related pathways and molecular mechanisms. Using the gene expression profiles and phenotypic groupings, the analysis was configured with a minimum gene set size of 5, a maximum gene set size of 5000, and 1000 permutations, with P  < 0.05 considered statistically significant 72 . Docking experiments were performed between the four PFAS and the three core targets. For receptor preparation, the identified core targets were imported into the RCSB PDB database ( https://www.rcsb.org/ ) to retrieve optimal protein conformations, which were saved in the PDB format. For ligand preparation, the 3D structures of PFAS obtained from PubChem were converted from SDF to MOL2 format using OpenBabel 3.1.1. Both receptor and ligand files were submitted to CB-Dock2 ( https://cadd.labshare.cn/cb-dock2/ ) for molecular docking 73 , 74 . The minimum binding energy values from the docking results were used for heatmap analysis, and binding conformations were visually represented. In this study, we constructed molecular complexes of PFOA with three key target proteins and performed 100-ns MD simulations using Gromacs 2022 software. Charmm 36 75 was chosen as the protein force field, Gaff2 was chosen as the ligand force field, and the TIP3P water model was chosen to add solvents to the protein-ligand system and create a water box with a periodic boundary of 1.2 nm 76 . Particle grid Ewald (PME) and Verlet algorithms were used to deal with electrostatic interactions. Then 100,000 steps of isothermal isovolumetric ensemble equilibrium and isothermal isobaric ensemble equilibrium were simulated with a coupling constant of 0.1 ps and a duration of 100 ps. Both van der Waals and Coulomb interactions were calculated using 1.0 nm cutoff values. Finally, the system was simulated using Gromacs 2022 at a constant temperature (310 K) and pressure (1 bar) for a total of 100 ns. The MM/PBSA method was used to calculate the binding free energy after stabilization of the system. KGN cells were seeded in 96-well plates at a density of 2 × 10³ cells/well and incubated for 24 h. The cells were then treated with 0, 0.2, 2, 20, and 200 µM PFOA for 24 and 48 h. And the cells were then treated with 0, 10, 100, 1000nM PFOA for 72 h. After cell adherence, 10 µL of CCK-8 reagent was added to each well, followed by incubation at 37 °C for 2 h. Absorbance at 450 nm was measured using a microplate reader (TECAN, Switzerland), and cell viability was calculated relative to that of the control group. Cell proliferation was assessed using an EdU Cell Proliferation Kit. KGN cells were seeded in 6-well plates at a density of 1 × 10 6 cells per well. After treatment with PFOA (0, 0.2, 2, and 20 µM) for 24 and 48 h, the cells were incubated with 10 µM EdU solution in the growth medium for 2 h. The cells were then stained with Azide 555 (red) and Hoechst 33,342 (blue). Finally, the cells were observed under a fluorescence microscope (Olympus IX73, Tokyo, Japan) at 100× magnification. EdU-positive cells were analyzed using the ImageJ software. KGN cells were seeded in 6-well plates at a density of 1 × 10 6 cells per well. After treatment with PFOA (0, 0.2, 2, and 20 µM) for 24 and 48 h, cell apoptosis was assessed using an Annexin V-FITC/PI Apoptosis Detection Kit, according to the manufacturer’s instructions. All cells were washed with PBS and resuspended in 100 µL Annexin V-binding buffer. Next, 5 µL of Annexin V-FITC and 5 µL of propidium iodide (PI) were added and incubated at room temperature for 10 min in the dark. Apoptotic cells were detected using a flow cytometer (Beckman Coulter, USA) and the results were analyzed using FlowJo software (version 10). Intracellular ROS levels were detected using DCFH-DA fluorescent probe in strict accordance with the manufacturer’s instructions. Briefly, KGN cells were seeded into 12-well plates and exposed to PFOA at concentrations of 0, 10, 100, and 1000 nM for 72 h, respectively. After the treatment period, DCFH-DA was added to each well at a final concentration of 10 µM, and the cells were incubated in the dark at 37 °C for 30 min. Fluorescence was measured under a fluorescence microscopy. KGN cells were seeded in 6-well plates at a density of 1 × 10 6 cells per well. After treatment, total RNA was extracted from the harvested KGN cells using an RNA-Quick Purification Kit, followed by reverse transcription into cDNA and RT-qPCR. The reaction conditions were as follows: 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. Target mRNA levels were normalized to GAPDH and analyzed with the 2 −ΔΔCt method. The primer sequences used for RT-qPCR are listed in Table S1 . KGN cells were seeded in 6-well plates at a density of 1 × 10 6 cells per well. WB was performed to evaluate the expression levels of key target proteins, including INHBA, INHBB, BAX, Bcl-2, GSTP1, PRDX5, and TGFBR1 in KGN cells. After treatment, the harvested KGN cells were lysed using RIPA buffer supplemented with protease and phosphatase inhibitors on ice for 30 min. Equal amounts of protein were separated by SDS-PAGE and subsequently transferred to PVDF membranes. Following transfer, the membranes were blocked with 5% skimmed milk at room temperature for 1 h, followed by incubation with primary antibodies (1:1000 dilution) at 4 °C overnight. The membranes were washed with TBST buffer and incubated with a horseradish peroxidase (HRP)-conjugated secondary antibody (1:50000 dilution) at room temperature for 1 h. After additional TBST washes, the protein bands were visualized using an ECL reagent and scanned with a gel imaging system (Tanon). Band intensities were quantified using ImageJ software, with GAPDH serving as the internal control for normalization. KGN cells were seeded in 6-well plates at a density of 1 × 10 6 cells per well. After treatment, the culture medium from each group was collected and centrifuged at 3000 × g for 10 min at 4 °C to remove cellular debris. Subsequently, the concentration of AMH in the supernatant was determined using a human AMH ELISA kit according to the manufacturer’s instructions. The absorbance at 450 nm was measured using a microplate reader and the AMH concentration in each sample was calculated based on a standard curve. All data were obtained from at least three independent biological replicates (i.e., independent experiments performed on different cell passages and batches at different times). For each biological replicate, three technical replicates were included, and their average value was used for subsequent analysis. Data are expressed as the mean ± SD. Statistical analyses were performed using GraphPad Prism v.10. The assumptions for parametric tests were systematically verified for each dataset. Specifically, the normality of the data distribution was confirmed using the Shapiro-Wilk test ( P  > 0.05), and the homogeneity of variance was verified using the Brown-Forsythe or Bartlett’s test ( P  > 0.05). All analyzed datasets satisfied these criteria. For comparisons between two groups, Student’s t-test was employed. For comparisons among three or more groups, one-way analysis of variance (ANOVA) was performed followed by Tukey’s post-hoc test (for multiple comparisons) or Dunnett’s test (if comparing to a control). If variances were unequal (heteroscedasticity), Welch’s correction was applied. Differences were considered statistically significant at P  < 0.05. Significance levels are indicated as * P  < 0.05, ** P  < 0.01, *** P  < 0.001, and **** P  < 0.0001.

In conclusion, this study provides a systematic exploration of granulosa cell (GC) dysfunction associated with selected PFAS congeners through an integrated multi-dimensional strategy. By combining network toxicology and machine learning with predictive MD simulations, we proposed a potential mechanistic framework for molecular impairments in human GCs following PFOA exposure. In vitro validation in KGN cells demonstrated that PFOA exposure was associated with activated apoptotic signaling and inhibited cell proliferation. Crucially, PFOA exposure led to altered expression of key proteins involved in redox homeostasis (GSTP1/PRDX5) and pro-fibrotic signaling (TGFBR1) within the GC compartment, while concurrently reducing the secretion of essential GC-derived ovarian reserve markers, such as AMH and INHB. Despite these insights, several limitations warrant consideration. First, the relatively modest sample size of the discovery transcriptomic dataset ( n  = 6 per group) may influence the generalizability of our model; thus, the identified genes should be interpreted as potential mechanistic signatures rather than validated biomarkers. Second, while molecular docking and MD simulations suggest the plausibility of potential interactions, these in silico predictions do not constitute confirmed biochemical binding. Regarding chemical selection, this study prioritized PFOA as a representative long-chain congener; however, caution must be exercised when generalizing these PFOA-specific findings to the entire structurally diverse PFAS class. Additionally, our investigation was limited to the 2D KGN cell line, which cannot recapitulate the 3D complexity of ovarian tissue or tissue-level fibrotic processes. Specifically, as direct fibrotic endpoints such as collagen deposition and SMAD signaling were not assessed, the upregulation of TGFBR1 should be interpreted as an activation of pro-fibrotic molecular signatures rather than confirmed ovarian fibrosis. Finally, the absence of functional perturbation assays prevents definitive confirmation of the causal requirement of GSTP1, PRDX5, and TGFBR1. Nonetheless, their consistent dysregulation establishes these genes as high-confidence candidate signatures of GC-mediated ovarian impairment. Future functional studies and histological validation in in vivo models will be instrumental in translating these correlative molecular findings into targeted clinical interventions.

PFAS are a group of chemicals known for their multisystem toxicity and capacity to damage multiple organs. They primarily affect the endocrine system, reproductive system, immune system, and liver, while increasing the risk of cancer. In recent years, increasing research has focused on the link between PFAS exposure and reproductive system diseases. For instance, one study reported that increased exposure to PFOA and PFOS significantly elevated the risk of POI and showed a significant negative correlation with AMH levels 30 . Furthermore, PFOS and its chlorinated ether analog 6:2 Cl-PFESA can impair germ cell differentiation, reduce the number of functional oocytes, and disrupt reproductive capacity in a dose-dependent manner 81 . A study found that PFOA exposure can induce granulosa cell apoptosis by downregulating the SIRT1/FOXO1-SOD2 pathway, leading to impaired mitochondrial antioxidant capacity 82 . Granulosa cell apoptosis is a key factor in PCOS pathogenesis. Previous research has established the impact of PFAS exposure on the development and progression of female reproductive endocrine disorders such as PCOS, endometriosis, and POI 83 . However, studies investigating the mechanisms through which PFOA influences DOR development remain limited. While several recent studies have explored the impact of PFAS on ovarian function, our work provides a more integrated and experimentally validated framework. For instance, Shen et al. (2023) 32 utilized clinical metabolomics to associate PFOA with metabolic shifts in DOR patients, but did not explore the specific protein-level drivers. Gao et al. (2024) 54 investigated the PFOS-induced PI3K/AKT pathway in KGN cells through traditional toxicology, but lacked the global target discovery and structural insights provided by our integrated workflow. Most notably, while Che et al. (2025) 84 employed network toxicology and MD simulations to study PFAS-induced POI, their findings remained purely computational. The lack of wet-lab validation in their study leaves the identified mechanisms as theoretical predictions. Our research addresses this critical gap by not only employing advanced Machine Learning algorithms (RF-LASSO-XGBoost) for superior target prioritization but also providing rigorous invitro validation in KGN cells. We experimentally confirmed that the identified hub targets (GSTP1, PRDX5, and TGFBR1) are indeed modulated by PFOA, thereby bridging the gap between ‘in silico’ prediction and biological reality. This experimental confirmation, combined with our focus on DOR (an earlier clinical stage than POI), ensures that our study is a substantial advancement rather than a repetition of existing work.In this study, we preliminarily evaluated the toxicity potential of PFAS using ProTox 3.0, ADMETlab 3.0, Vnn-ADMET, and admetSAR 3.0. A network toxicology approach was subsequently employed for further analysis, involving screening of multiple databases, including ChEMBL, STITCH, SwissTargetPrediction, PharmMapper, GeneCards, and OMIM, which led to the identification of 35 candidate targets associated with PFAS-induced DOR. Subsequently, PPI network analysis was performed to construct a relational network of core targets. GO analysis revealed that the key targets were 482 BP, 25 CC, and 81 MF. KEGG enrichment analysis demonstrated that the key targets were primarily enriched in pathways regulating cell proliferation and differentiation (such as PI3K-Akt, MAPK, Ras, and ErbB signaling), immune-inflammatory responses (e.g., Th17 cell differentiation), and hormone signaling pathways (including prolactin, estrogen, and progesterone signaling). These results collectively implicate a cellular signaling network that is dependent on growth factor receptors and sex hormones. Additionally, single-gene GSEA revealed that the expression of the three core genes was associated with signal transduction and metabolic pathways, including the ErbB signaling pathway, insulin signaling pathway, phosphatidylinositol signaling system, and pyrimidine metabolism. These findings suggest that PFAS may play a significant role in the development and progression of DOR by modulating these pathways. Through RNA-seq analysis of ovarian granulosa cells from both the normal and DOR populations, we identified 5836 differentially expressed genes, comprising 2593 upregulated and 3243 downregulated genes. The intersection between these differentially expressed genes and the 35 target genes yielded 19 overlapping genes. Subsequently, progressive machine learning screening of 19 overlapping genes identified three core targets: GSTP1, PRDX5, and TGFBR1. To address the inherent risk of model overfitting often associated with small-sample transcriptomic studies, we performed a rigorous external validation of the machine-learning-identified core targets across two independent clinical cohorts. The high diagnostic performance observed in GSE315266 (AUC values of 0.75 for GSTP1, 0.78 for PRDX5, and 0.75 for TGFBR1) and the consistent results in GSE296744 underscore the remarkable cross-cohort stability and generalizability of our findings. Notably, TGFBR1 exhibited exceptional robustness, achieving a high AUC of 0.819 in the second validation dataset. This suggests that the activation of TGF-β signaling and its mediated fibrotic processes represent a fundamental and highly conserved pathogenic driver in DOR, regardless of the clinical platform or patient population. While some fluctuations in the AUC values of GSTP1 and PRDX5 were observed (0.639 and 0.667 in GSE296744 ), such variability is biologically expected due to the clinical heterogeneity of DOR patients and differences in detection sensitivities across transcriptomic platforms. Therefore, we propose that these three targets are central to the PFAS-induced DOR. Molecular docking results demonstrated that all four PFAS compounds exhibited strong binding affinity (binding energies ranging from − 9.8 to -6.5 kcal/mol) towards the three core proteins (GSTP1, PRDX5, TGFBR1), suggesting their potential involvement in PFAS-induced DOR. Furthermore, a 100-ns MD simulation of PFOA, as a representative PFAS compound, demonstrated that it maintained structural stability over time, thereby supporting the rationality of the docking predictions. Although molecular docking can preliminarily reveal potential interaction patterns between ligands and targets, its predictive reliability is limited by the absence of structurally analogous reference compounds. This study further evaluated the dynamic stability of the predicted complexes by incorporating MD simulations; however, additional comparative studies are required to assess their binding specificity and selectivity more accurately. In summary, these findings suggest that PFAS may interact with key molecular targets involved in cell proliferation and differentiation, inflammatory responses, and hormone signal transduction, thereby providing a computational foundation for understanding the potential pathogenic mechanisms of DOR. The KGN cell line selected as the research model is widely used in ovarian studies because of its physiological similarity to the ovarian granulosa cells 85 . Furthermore, we selected PFOA as the most representative PFAS compound for in vitro experimental validation. The results demonstrated that PFOA exposure significantly induced KGN cell injury in a dose- and time-dependent manner, characterized by reduced cell viability, suppressed proliferation, and increased apoptosis, indicating its cytotoxic effects. Previous studies have shown that PFAS exert cytotoxicity in Human Kidney-2 (HK-2) cells in a dose- and chain length-dependent manner 86 , while also enhancing cholesterol accumulation in macrophages, thereby triggering dysregulated inflammatory responses 87 . This demonstrating that PFOA produces similar cytotoxic effects in various human cell types. However, considering the structural diversity of the PFAS family (e.g., chain length and functional groups), these PFOA-specific effects may not be universally applicable to the entire class, necessitating future comparative studies with emerging short-chain PFAS. Li et al. established a human trophoblast organoid model with a near-physiological ratio of extravillous trophoblasts and reported that PFOA exhibited a stronger inhibition of trophoblast differentiation than PFOS 88 . Furthermore, PFOA can induce apoptosis in KGN cells by disrupting mitochondrial function via the SIRT1/FOXO1-SOD2 pathway 82 . These findings are consistent with our observations, demonstrating that PFOA produces similar cytotoxic effects in various human cell types under specific exposure conditions. To further elucidate the molecular mechanisms of PFOA-induced DOR, we integrated network toxicology predictions with experimental validation, focusing on expression changes in genes related to apoptosis and oxidative stress as well as ovarian reserve biomarkers. Apoptosis is a primary cause of cellular damage. Network toxicology analysis predicted the involvement of signaling pathways related to proliferation and apoptosis, such as PI3K-Akt, MAPK, and Ras, which was closely aligned with our experimental observations. Furthermore, our results demonstrated that PFOA suppressed KGN cell viability and proliferation in a dose- and time-dependent manner. Furthermore, PFOA significantly upregulated the expression of the apoptosis-related marker BAX while downregulating Bcl-2 expression, resulting in an increased BAX/Bcl-2 ratio, indicating that PFOA exposure activates apoptotic signaling pathways in KGN cells. Bcl-2 is a classic anti-apoptotic protein that plays a crucial role in granulosa cell survival. Under conditions of diminished ovarian function, Bcl-2 expression frequently decreases. For instance, in a DOR mouse model, suppression of Bcl-2 accelerated ovarian cell apoptosis, thereby contributing to diminished ovarian reserve 89 . Additionally, the upregulation of BAX contributes to the progression of apoptosis. In the 4-Hydroperoxy-Cyclophosphamide (4-HC)-induced DOR cell model, KGN cells exhibited elevated BAX expression levels, accompanied by reduced Bcl-2 levels, thereby inducing apoptosis 90 . These findings suggest that PFOA exposure may activate both the intrinsic and extrinsic apoptotic cascades in KGN cells, ultimately enhancing apoptosis. Glutathione S-transferase Pi 1 (GSTP1) is an enzyme widely present in humans that neutralizes harmful substances by catalyzing the conjugation of electrophilic compounds with glutathione 91 . Thus, it participates in detoxification processes including ROS clearance. An imbalance between ROS production and elimination, which leads to cellular damage, is a direct cause of oxidative stress. Previous studies have demonstrated that oxidative stress is a critical factor in ovarian aging and can accelerate the decline in ovarian function by promoting apoptosis, inflammatory responses 92 , 93 , and mitochondrial dysfunction 94 . It has been reported that GSTP1 can bind to and inhibit the activity of the pro-apoptotic protein c-Jun N-terminal kinase (JNK), thereby negatively regulating the apoptotic pathway 95 . Furthermore, studies have indicated that tryptanthrin can induce senescence in hepatocellular carcinoma cells via the GSTP1/ROS/DDR/NF-κB/SASP signaling axis 96 . On the other hand, inhibition of GSTP1 not only leads to excessive ROS accumulation but also accelerates the senescence process of regulatory T cells, thereby promoting inflammatory dysregulation and the onset of immunosenescence 97 . Therefore, the downregulation of GSTP1 may result in ROS accumulation, triggering sustained oxidative stress and immune dysfunction, which in turn accelerates ovarian aging. Our study found that the expression of GSTP1 was significantly downregulated ( P  < 0.01) in the ovarian granulosa cells of patients with DOR compared to that in normal women. Additionally, in vitro experiments revealed that prolonged exposure to high concentrations of PFOA downregulates GSTP1 expression. These findings suggest that GSTP1 may be a potential target for PFOA-induced senescence in granulosa cells. PRDXs (Peroxiredoxins) are a family of thiol peroxidases that play key roles in maintaining redox homeostasis. Peroxiredoxin V (PRDX5) is an isoform expressed in the mitochondria, peroxisomes, cytoplasm, and nucleus 98 . Mature oocytes contain approximately 100,000 mitochondria, a reserve that meets the ATP demands for oocyte maturation and provides energy support for early embryonic development post-fertilization 94 . Oocytes and their surrounding granulosa cells form a cumulus-oocyte complex (COC). Granulosa cells can deliver energy substrates, such as pyruvate and lactate, produced from glucose metabolism, to the oocyte via gap junctions, providing ATP for maturation. Furthermore, granulosa cells supply essential precursors for the synthesis of genetic material, proteins, and lipids in oocytes 99 . Dysfunction of granulosa cells inevitably leads to abnormalities in mitochondrial biogenesis and energy metabolism within the oocytes, thereby impairing oocyte maturation and development. Jin et al. demonstrated that PRDX5 overexpression eliminates intracellular ROS accumulation and regulates the classical apoptotic pathway, thereby inhibiting emodin-induced apoptosis in gastric cancer cells 100 . Conversely, studies by Wu et al. revealed that Prdx5 promotes M1 polarization and apoptosis in prostate epithelial cells via the TLR4/NF-κB axis in a ROS-dependent manner 101 , indicating its active role in oxidative stress and immune regulation during disease pathogenesis. In contrast to normal women, the present study observed significantly downregulated expression of PRDX5 ( P  < 0.01) in the ovarian granulosa cells of patients with DOR. Furthermore, our in vitro experiments demonstrated that prolonged exposure to high PFOA concentrations downregulated PRDX5 expression. Consequently, we hypothesized that PFOA may contribute to diminished ovarian reserves by suppressing PRDX5 expression in granulosa cells, thereby disrupting mitochondrial energy metabolism in oocytes. Transforming growth factor-β receptor type 1 (TGFBR1) is a key transmembrane serine/threonine kinase involved in the TGF-β superfamily signaling pathway. The heterodimer formed by TGFBR1 and TGFBR2 constitutes the TGF-β receptor. As a central mediator of both canonical Smad-dependent and non-canonical (e.g., MAPK and PI3K/Akt) signaling pathways, it regulates a series of cellular processes crucial for reproductive physiology, including cell differentiation, proliferation, apoptosis, wound healing, and tissue fibrosis, with particular importance in ovarian folliculogenesis and the maintenance of ovarian reserve 102 – 104 . In cardiac disease, TGFBR1 primarily contributes to myocardial fibrosis, and inhibition of its expression has been shown to reduce collagen deposition via the Smad signaling pathway, thereby improving cardiac function 102 , 105 . In chronic kidney disease, aberrant activation of TGF-β accelerates renal fibrosis progression, and interferon regulatory factor 5 has been found to reverse this process by suppressing TGFBR1 transcription 106 . A recent study revealed that fibrosis is not only a hallmark of aging but also serves as a primary cause of ovarian dysfunction. In the ovary, excessive deposition of extracellular matrix components (e.g., collagen) increases tissue stiffness, thereby physically impeding follicular growth and disrupting the essential biomechanical signaling required for follicle-stroma crosstalk. Moreover, pathologically elevated levels of pro-fibrotic factors, particularly TGF-β, can directly lead to impaired follicular development and ovulation 107 . Compared to normal women, this study revealed a significant upregulation of TGFBR1 expression ( P  < 0.01) in the ovarian granulosa cells of patients with DOR. Furthermore, our in vitro experiments demonstrated that prolonged exposure to high PFOA concentrations upregulated TGFBR1 expression. Thus, we hypothesized that PFOA exposure leads to aberrant activation of TGFBR1, accelerating the fibrotic process within the follicular microenvironment and consequently impairing follicular growth and development in the ovary. The Stages of Reproductive Aging Workshop + 10 (STRAW + 10) criteria recognize AMH and INHB as important biomarkers for assessing the ovarian reserve. As a member of the TGF-β superfamily, AMH primarily functions to diminish the stimulatory effects of FSH on primordial and small antral follicles, thereby inhibiting the recruitment of primordial follicles and the development of antral follicles, which helps prevent premature follicle depletion 108 . Moreover, a key advantage of AMH over other hormonal markers in predicting ovarian reserve is its independence from the menstrual cycle as AMH levels remain relatively stable throughout the menstrual cycle. INHB, also a member of the TGF-β superfamily, is a dimeric peptide (composed of α- and β-subunits) secreted by the granulosa cells of small antral follicles. As INHB levels fluctuate with the menstrual cycle, some studies have questioned its value as a potential biomarker 109 . Nevertheless, its selective inhibitory effect on FSH, along with findings from recent studies, support its utility in predicting ovarian reserve function 80 , 110 . In our study, PFOA exposure downregulated the expression of both AMH and INHB in KGN cells in a dose- and time-dependent manner, suggesting impaired hormone secretion by the granulosa cells. Although histological evidence could not be obtained in vitro, the observed changes in AMH and INHB levels, along with concurrent alterations in apoptosis and oxidative stress markers, support the notion that PFOA may contribute to the risk of diminished ovarian reserve. Further in vivo studies will help to clarify the extent of granulosa cell injury and validate these molecular findings. A key challenge in environmental toxicology is reconciling experimental dosages with real-world exposure levels. In this study, while initial mechanistic screening was performed at µM levels to capture acute apoptotic events, our subsequent low-dose validation (10–1000 nM) confirmed that functional impairment occurs at much lower concentrations. Specifically, the significant reduction in AMH secretion and the induction of ROS at 100 nM—a concentration close to those detected in clinical samples—demonstrate that PFOA directly compromises granulosa cell function before triggering cell death. This ‘low-dose functional sensitivity’ suggests that the identified GSTP1/PRDX5/TGFBR1 axis is a reliable molecular signature of PFAS-induced DOR, even under chronic, low-level environmental exposure. In summary, PFOA exposure in KGN cells disrupted multiple cellular processes, including suppression of cell proliferation, increased apoptosis, induction of oxidative stress, interference with TGFβ pathway protein expression, and impaired secretion of ovarian reserve-related hormones. Mechanistically, upregulation of BAX coupled with downregulation of Bcl-2 promotes apoptotic signaling. Reduced expression of GSTP1 and PRDX5 compromises cellular antioxidant defenses, leading to an accumulated oxidative burden. Conversely, upregulation of TGFBR1 exacerbates fibrotic processes within the follicular microenvironment, accelerating ovarian aging. The decreased expression of AMH and INHB, key biomarkers of ovarian reserve, directly reflects the diminished hormonal secretory capacity of the granulosa cells. Collectively, these molecular alterations are likely to synergistically impair granulosa cell functions. Consistent with previous reports on PFAS-related ovarian toxicity, our findings underscore the potential risk of PFOA exposure contributing to the DOR.

The rapid progression of global industrialization and urbanization has led to widespread environmental contamination, which poses a significant threat to ecosystem integrity and public health 1 . Among various environmental pollutants, endocrine-disrupting chemicals (EDCs) are of increasing concern because of their ability to perturb endocrine homeostasis. According to the Endocrine Society, EDCs are exogenous substances that interfere with hormone action, thereby disrupting homeostatic systems that allow an organism to respond to its environment 2 . Many EDCs, such as the fungicide vinclozolin 3 , pesticides, including DDT and chlorpyrifos 4 , and plasticizers, such as perfluoroalkyl substances (PFAS) and bisphenol-A 5 , 6 , are pervasive in industrial and consumer products. A substantial body of evidence indicates that these disruptors exert adverse effects on endocrine, neurological, and reproductive systems by modifying hormone receptors, disrupting signaling pathways, and altering gene expression 7 – 9 . Per- and polyfluoroalkyl substances (PFAS), a class of EDCs extensively used in industrial, agricultural, and consumer products, have attracted significant attention owing to their extreme environmental persistence. These compounds are characterized by remarkable chemical and thermal stability, rendering them highly resistant to environmental degradation 10 . Their hydrophilic nature and high environmental mobility, coupled with their persistence and strong potential for bioaccumulation, facilitate transport across various environmental media and accumulation in living organisms. These properties collectively underpin their designation as “forever chemicals” 11 . Environmental monitoring has detected PFAS in multiple pathways, including drinking water sources and the food chain, often via contamination from packaging materials 12 . Consequently, the primary human exposure routes have been identified as the ingestion of contaminated drinking water and food supplemented by cumulative intake off air and indoor dust 13 . Epidemiological studies have confirmed that PFAS exert multiple toxic effects including, but not limited to, neurotoxicity 14 , immunotoxicity 15 , potential carcinogenicity 16 , endocrine toxicity 17 , reproductive toxicity 18 , and respiratory toxicity 19 . These compounds are widely distributed in the liver, kidneys, and blood 20 , 21 , and have also been detected in human urine, amniotic fluid, follicular fluid (FF), and placental tissues 22 . Critically, PFAS cross the placental barrier and transfer through breastfeeding, posing exposure risks during critical developmental windows 23 . This dual exposure pathway may lead to sustained chemical exposure during critical developmental windows, warranting an in-depth investigation of its potential developmental toxicity. Given their ecological and human health risks, perfluorooctanoic acid (PFOA) and perfluorooctane sulfonate (PFOS) are being phased out globally and are listed as Persistent Organic Pollutants under the Stockholm Convention, a regulatory stance mirrored by China’s inclusion, and perfluorohexane sulfonic acid (PFHxS) on its Key Controlled New Pollutants List. Accumulating evidence has linked PFAS exposure to adverse female reproductive outcomes. Mechanistically, long-chain PFAS disrupt follicular maturation and hormone secretion by activating PPARγ in granulosa cells, leading to anovulation and reduced fertility 24 . At human-relevant concentrations, PFAS mixtures inhibit mitochondrial function and oxidative phosphorylation, reducing ATP production and steroidogenesis, while downregulating key enzymes such as 3βHSD and SOD1. Additionally, PFDA promotes granulosa cell necroptosis and follicular atresia via the RIPK1 pathway 25 . Epidemiologically, PFAS exposure is associated with an increased risk of ovarian disorders, such as PCOS 26 , 27 , endometriosis 28 , 29 , premature ovarian insufficiency (POI) 30 , and DOR 31 , 32 . DOR is characterized by a reduction in the quantity and/or quality of oocytes in the ovary. The diagnostic criteria typically include an AMH level < 1.1 ng/ml, an antral follicle count (AFC) < 7, or a basal follicle-stimulating hormone (FSH) level ≥ 10 IU/L. Although its exact etiology remains unclear, current evidence suggests that the development of DOR may involve multifactorial interactions, including genetic predisposition, autoimmune disorders, infections, smoking, medication use, and exposure to environmental pollutants 33 . Numerous studies have demonstrated a close association between PFAS exposure and the development and progression of DOR. A study focusing on Chinese women revealed that elevated levels of perfluorohexanoic acid (PFHxA) in follicular fluid were significantly correlated with DOR 34 . Adolescent exposure to PFAS mixtures (PFOA, PFOS, GenX/HFPO-DA, PFBS) was found to impair folliculogenesis by activating Hippo/Yap1 signaling, thereby compromising ovarian function in adulthood 35 . Furthermore, PFOA has been shown to disrupts steroid hormone synthesis and normal follicular development by interfering with the hypothalamic-pituitary-gonadal (HPG) axis 36 , 37 . Given that oocyte maturation is highly dependent on gap junctional intercellular communication (GJIC) between oocytes and granulosa cells, López et al. (2019) discovered that PFOA may induce oocyte death by blocking GJIC within the cumulus-oocyte complex (COC) 38 . Mitochondria serve as the primary energy suppliers for oocyte development and maturation, whereas reactive oxygen species (ROS) are natural byproducts of mitochondrial metabolism. Cells maintain redox homeostasis through sophisticated antioxidant systems that precisely regulate ROS concentration and protect against oxidative damage. Zhang et al. demonstrated that increased PFOA exposure compromises oocyte stability and viability by inducing mitochondrial damage, thereby reducing ATP and DNA production in oocytes 39 . In a separate study, López et al. observed that PFOA exposure significantly elevated ovarian ROS levels, and this ROS elevation triggered mitochondrial dysfunction while suppressing oocyte growth 38 . Furthermore, PFOA induces oxidative stress by inhibiting mitochondrial energy production, leading to mitochondrial dysfunction and subsequent cellular apoptosis. These processes involve multiple signaling pathways, including kisspeptin and PPAR-α signaling pathways 40 , 41 . Collectively, these studies demonstrate that PFAS significantly impair oocyte viability, thereby compromising ovarian reserve function. Studies of women undergoing in vitro fertilization-embryo transfer (IVF-ET) consistently detected PFAS in serum and FF, with PFOS and PFOA as the predominant compounds. Chinese cohorts report relatively high PFOA levels (e.g., median ~ 5.56–5.85 ng/mL in FF), whereas studies from Australia, Europe, and the United States typically show lower PFOA concentrations (e.g., mean ~ 1.80–2.40 ng/mL in FF) but higher PFOS levels in some regions (e.g., mean up to 7.5 ng/mL in Belgian FF) 42 – 47 . This geographical variation, potentially reflecting regional production and use, was corroborated by the European Food Safety Authority (EFSA) assessment, which identified PFOS (64%) and PFOA (16%) as the most prevalent PFAS in adult serum (EFSA CONTAM Panel). Computational toxicology demonstrates significant technical advantages in deciphering compound toxicity mechanisms through its capacity to integrate multisource bioinformatics data 48 . By enabling the systematic visualization of toxicological mechanisms, it transforms complex multi-target synergistic or antagonistic interactions from traditional toxicological research into network maps with well-defined topological features. This facilitates multi-level mechanistic interpretation, ranging from molecular interactions to pathway regulation, thereby aiding the discovery of biomarkers and potential therapeutic interventions 49 – 51 . To enhance the reliability of target identification and interaction prediction, molecular docking and molecular dynamics (MD) simulations have been employed as complementary validation tools 52 , 53 . In this study, we further applied MD simulations to verify the interaction patterns and structural stability of PFAS-target complexes predicted through molecular docking, thereby providing comprehensive dynamic insights into potential toxicological risks. Human ovarian granulosa cells (hGCs) are core functional units that regulate ovarian physiology. They are not only responsible for synthesizing and secreting hormones to maintain follicular development and maturation but also serve as primary targets for ovarian toxicity induced by various endocrine-disrupting chemicals, including environmental pollutants. The KGN cell line, an immortalized human ovarian granulosa cell model, has been extensively used as a representative in vitro system for investigating ovarian physiological functions, assessing ovarian toxicity, and exploring hormonal regulatory mechanisms. This widespread adoption stems from the retention of key morphological and functional characteristics inherent to the native human ovarian granulosa cells 54 – 56 . In this study, four compounds (PFOS, PFOA, PFHxS, and PFNA) were selected as representative substances for computational toxicological investigation. Among these, PFOA was specifically prioritized for in vitro experimental validation. This choice was informed by its status as a prototypical long-chain legacy PFAS with high environmental persistence and widespread biological detection. Importantly, PFOA has been identified as the most abundant PFAS congener in the follicular fluid of Chinese women, making it a highly clinically relevant model for investigating the reproductive toxicological mechanisms of this class of chemicals. It should be noted that this study primarily utilizes the human KGN granulosa cell line for experimental validation. Although KGN cells provide a robust and validated platform for exploring molecular signaling and GCs function, they do not fully capture the multicellular complexity or the follicular-stromal interactions of the intact ovary. Therefore, our findings are framed within the context of GC-mediated mechanisms of ovarian injury.

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1