Causality and potential drug targets for hyperandrogenism and anovulation-related infertility: two-sample Mendelian randomization and colocalization analyses.

Filter out SNPs not corresponding to exposure at the genome-wide level ( P  < 5.0 × 10 − 8 ), and exclude SNPS with linkage imbalances, there were 299, 262, 4, and 439 SNPs used as instruments for TT, BIOT, DHEAS, and SHBG, respectively, conforming to the correlation hypothesis and independence hypothesis. The F statistics for each instrument were found to be greater than 10, and except for TT (Leinonen) and DHEAS, the power of other instruments was greater than 99%, demonstrating the significant robustness of the genetic instruments utilized, as depicted in Table S3. The IVW analysis results showed that TT (OR1 = 0.918; 95% CI: 0.707 to 1.192; P1  = 0.521; OR2 = 0.982; 95% CI: 0.747 to 1.293; P2  = 0.899), BIOT (OR1 = 1.331; 95% CI: 0.928 to 1.909; P1  = 0.121; OR2 = 1.279; 95% CI: 0.972 to 1.684; P2  = 0.079), DHEAS (OR = 1.488; 95% CI: 0.980 to 2.260; P  = 0.062) were not associated with the anovulation-related female infertility; SHBG had negative causality with the incidence of anovulation-related female infertility (OR1 = 0.561; 95% CI: 0.363 to 0.866; P1  = 0.009; OR2 = 0.878; 95% CI: 0.726 to 1.062; P2  = 0.180)(Fig.  2 ). Sensitivity analysis confirmed a strong relationship between SHBG and infertility, with no significant pleiotropy or heterogeneity (Table  1 ). Leave-one-out analysis yielded similar results (Fig. S1-S4) and the funnel plots were shown as Fig S5-S8. Other methods (MR-Egger, WM, simple mode, weighted mode) were also used. Fig. 2 Forest plots of the association between TT, BIOT, DHEAS, SHBG and anovulation-related infertility. The results of the UK-Biobank database showed the corrected ORs and 95% CIs. P- value<0.05 indicates statistical significance. OR: odds ratio; CI: confidence interval; TT: total testosterone; BIOT: bioavailable testosterone; DHEAS: dehydroepiandrosterone sulfate; SHBG: sex hormone-binding globulin; WM: weighted median; IVW: inverse variance weighted Forest plots of the association between TT, BIOT, DHEAS, SHBG and anovulation-related infertility. The results of the UK-Biobank database showed the corrected ORs and 95% CIs. P- value<0.05 indicates statistical significance. OR: odds ratio; CI: confidence interval; TT: total testosterone; BIOT: bioavailable testosterone; DHEAS: dehydroepiandrosterone sulfate; SHBG: sex hormone-binding globulin; WM: weighted median; IVW: inverse variance weighted Table 1 Sensitivity analysis Outcome Study Exposure P for heterogeneity test P for pleiotropy test Female infertility, associated with anovulation Ruth et al. (2020) TT 0.09 0.305 Leinonen et al. (2023) TT 0.387 0.017* Ruth et al. (2020) BIOT 0.241 0.316 Leinonen et al. (2023) BIOT 0.131 0.632 Prins et al. (2017) DHEAS 0.702 0.721 Ruth et al. (2020) SHBG 0.339 0.583 Leinonen et al. (2023) SHBG 0.848 0.183 TT  total testosterone, BIOT  bioavailable testosterone, DHEAS  dehydroepiandrosterone sulfate, SHBG  sex hormone-binding globulin Sensitivity analysis TT  total testosterone, BIOT  bioavailable testosterone, DHEAS  dehydroepiandrosterone sulfate, SHBG  sex hormone-binding globulin Androgen indices with two or more reliable MR results were integrated by meta-analysis. The meta-analysis results showed a potential positive association between BIOT (OR = 1.298; 95% CI: 1.043 to 1.615; P  = 0.019) and female infertility associated with anovulation, as well as a potential negative association between SHBG (OR = 0.817; 95% CI: 0.686 to 0.972; P  = 0.023) and such infertility (Fig.  3 ). Fig. 3 Forest plots of meta-analysis about the causality between the androgen indices and anovulation-related infertility. The results of the UK-Biobank database showed the corrected ORs and 95% CIs. P- value<0.05 indicates statistical significance. OR: odds ratio; CI: confidence interval; TT: total testosterone; BIOT: bioavailable testosterone; SHBG: sex hormone-binding globulin; WM: weighted median; IVW: inverse variance weighted; P h : P for heterogeneity Forest plots of meta-analysis about the causality between the androgen indices and anovulation-related infertility. The results of the UK-Biobank database showed the corrected ORs and 95% CIs. P- value<0.05 indicates statistical significance. OR: odds ratio; CI: confidence interval; TT: total testosterone; BIOT: bioavailable testosterone; SHBG: sex hormone-binding globulin; WM: weighted median; IVW: inverse variance weighted; P h : P for heterogeneity To support associations and understand the biological basis, we performed Hetnet connectivity search. Using the “Androgen Receptor Signaling Pathway” as the source and the “Ovarian Infertility Gene” pathway as the target, we retrieved 200 meta-paths. We selected meta-paths with statistically significant adjusted P values and paths with a path score ≥ 10 for visualization (Fig.  4 ). Fig. 4 Using the connectivity search web to explore the relationship between androgen receptor signaling and anovulation-related infertility. This figure shows an example user workflow for https://het.io/search/ . A Node search ( B ) The web returns meta-paths between the androgen receptor signaling pathway and the pathway of ovarian infertility gene. Only the adjusted P -value for the first source path is statistically significant. P- value<0.05 indicates statistical significance. C Paths with path score ≥ 10 for the first meta-paths. D  A subgraph displays the previously selected paths. PW: pathway; G: gene; p: participates; i: interacts; r: regulates; c: covaries Using the connectivity search web to explore the relationship between androgen receptor signaling and anovulation-related infertility. This figure shows an example user workflow for https://het.io/search/ . A Node search ( B ) The web returns meta-paths between the androgen receptor signaling pathway and the pathway of ovarian infertility gene. Only the adjusted P -value for the first source path is statistically significant. P- value<0.05 indicates statistical significance. C Paths with path score ≥ 10 for the first meta-paths. D  A subgraph displays the previously selected paths. PW: pathway; G: gene; p: participates; i: interacts; r: regulates; c: covaries A comprehensive Bayesian colocalization analysis in 117 gene regions (Table S4) revealed that TT, BIOT, DHEAS, SHBG are associated with anovulation-related infertility at the gene loci of rs45446698, rs727428, rs148982377, and rs545206972, respectively. These genes share a common causal variant within the ± 200 kb region (Fig.  5 A-D). Fig. 5 Locus comparing plots for the shared causal variant for the associations of TT, BIOT, DHEAS, SHBG and anovulation-related infertility. A Colocalization analysis results for the association between TT and female infertility (associated with anovulation) in the gene region (Chr7: 99332948), located within ± 200 kb from rs45446698. B Colocalization analysis results for the association between BIOT and female infertility (associated with anovulation) in the gene region (Chr17:7537792), located within ± 200 kb from rs727428. C Colocalization analysis results for the association between DHEAS and female infertility (associated with anovulation) in the gene region (Chr7:99075038), located within ± 200 kb from rs148982377. D Colocalization analysis results for the association between SHBG and female infertility (associated with anovulation) in the gene region (Chr17:7491331), located within ± 200 kb from rs545206972 Locus comparing plots for the shared causal variant for the associations of TT, BIOT, DHEAS, SHBG and anovulation-related infertility. A Colocalization analysis results for the association between TT and female infertility (associated with anovulation) in the gene region (Chr7: 99332948), located within ± 200 kb from rs45446698. B Colocalization analysis results for the association between BIOT and female infertility (associated with anovulation) in the gene region (Chr17:7537792), located within ± 200 kb from rs727428. C Colocalization analysis results for the association between DHEAS and female infertility (associated with anovulation) in the gene region (Chr7:99075038), located within ± 200 kb from rs148982377. D Colocalization analysis results for the association between SHBG and female infertility (associated with anovulation) in the gene region (Chr17:7491331), located within ± 200 kb from rs545206972 The generated and merged Seurat object comprised 28,806 TCs and 27,692 GCs. We filtered cells based on UMI count (nFeature > 200) and mitochondrial gene read percentage (TCs<10%, GCs<20%). After excluding non-GCs and non-TCs, 22,712 TCs and 18,465 GCs remained for further analysis (Fig.  6 A-C). Dividing the dataset into normal and PCOS groups, GSEA identified 18 pathways with opposing enrichment directions in GCs and TCs (Table S5). The Rho/ROCK signaling pathway was significantly enriched in PCOS GCs, with the opposite enrichment in TCs (Fig.  6 E). Fig. 6 Single-cell type expression in ovary tissues for the protein-coding genes of the identified candidate targets. A Uniform manifold approximation and projection (UMAP) plot of 22,712 granulosa cells (GCs), 18,465 theca cells (TCs), colored by cell type. B UMAP plot colored by sample, indicating effective batch effect removal. C Dot plot showing average expression and cell fraction of selected genes in GCs and TCs. Colors indicate average expression levels (green for high, red for low), and dot sizes represent the fraction of cells. D The protein-protein interaction (PPI) network of potential drug target proteins. E GSEA enrichment for the REACTOME_RHO_GTPASE_CYCLE pathway in GCs and TCs. ( F , G ) Comparison of different pathways in the normal group and the polycystic ovary syndrome (PCOS) group in GCs F and TCs G . P- value<0.05 and the asterisks indicate statistical significance of the differences between groups, with *** representing P  < 0.001, and ns representing no significance. UMAP: uniform manifold approximation and projection; GC: granulosa cell; TC: theca cell; PPI: protein-protein interaction; GSEA: gene set enrichment analysis; PCOS: polycystic ovary syndrome Single-cell type expression in ovary tissues for the protein-coding genes of the identified candidate targets. A Uniform manifold approximation and projection (UMAP) plot of 22,712 granulosa cells (GCs), 18,465 theca cells (TCs), colored by cell type. B UMAP plot colored by sample, indicating effective batch effect removal. C Dot plot showing average expression and cell fraction of selected genes in GCs and TCs. Colors indicate average expression levels (green for high, red for low), and dot sizes represent the fraction of cells. D The protein-protein interaction (PPI) network of potential drug target proteins. E GSEA enrichment for the REACTOME_RHO_GTPASE_CYCLE pathway in GCs and TCs. ( F , G ) Comparison of different pathways in the normal group and the polycystic ovary syndrome (PCOS) group in GCs F and TCs G . P- value<0.05 and the asterisks indicate statistical significance of the differences between groups, with *** representing P  < 0.001, and ns representing no significance. UMAP: uniform manifold approximation and projection; GC: granulosa cell; TC: theca cell; PPI: protein-protein interaction; GSEA: gene set enrichment analysis; PCOS: polycystic ovary syndrome Comparing pathways in normal and PCOS groups of GCs and TCs, Ucell scores revealed that in GCs, the PCOS group demonstrated increased circulating androgen levels, as well as unregulated androgen synthesis pathways, PI3K/AKT pathway, MAPK1/ERK2 pathway, TGF-β signaling, AMH pathway, and heightened activity of the AR (Fig.  6 F). In TCs, the PCOS group exhibited elevated circulating androgen levels, unregulated androgen synthesis pathways, but down-regulated AR signaling and related pathways (Fig.  6 G). PPI network analysis identified 29 pairs of interactions among 21 proteins with a confidence score threshold set at 0.4. Significant interactions were observed between CD74 and HLA-DRA (score = 0.999), CD74 and HLA-DRB1 (score = 0.995), ATP5F1B and ATP5MG (score = 0.999), CYP11A1 and STAR (score = 0.979), CYP11A1 and FDXR (score = 0.997), RPS15 and RPS9 (score = 0.999) (Fig.  6 D). In the druggability assessment, three of 27 PCOS-associated plasma proteins, including CD74, CYP11A1, and HLA-DRB1 were identified as potential candidates for pharmaceutical research and development. Repotrectinib targeting CD74, mitotane targeting CYP11A1, and apolizumab targeting HLA-DRB1 have been used in treating non-small cell lung carcinoma, adrenal cortex neoplasm, and lymphoma, respectively (Table  2 ). Table 2 Genetic information and findings for key proteins associated with PCOS Protein Name Gene Name Co-regulated Theca cells Granulosa cells Drug Development Target Development Level avg_log2FC p_val_adj avg_log2FC p_val_adj Drug Name Outcomes Trial Phase Drug-Tclin a Ligands-Tchem b ATP5F1B ATP5F1B Up 0.297886946 < 0.01 0.393176332 < 0.01 - - - 0 0 MSMO1 MSMO1 Up 0.290422522 < 0.01 0.331270587 < 0.01 - - - 0 0 ATP5MG ATP5MG Up 0.408448077 < 0.01 0.396877014 < 0.01 - - - 0 0 CD74 CD74 Up 0.615847607 < 0.01 0.526016586 < 0.01 Repotrectinib Non-small Cell Lung Carcinoma III 1 48 CYP11A1 CYP11A1 Up 0.479704339 < 0.01 0.392701027 < 0.01 Mitotane Adrenal Cortex Neoplasm IV 2 0s FDXR FDXR Up 0.283557659 < 0.01 0.25353556 < 0.01 - - - 0 0 HLA-DRA HLA-DRA Up 0.496442517 < 0.01 0.490027772 < 0.01 - - - 0 0 HLA-DRB1 HLA-DRB1 Up 0.587321294 < 0.01 0.885017617 < 0.01 Apolizumab Lymphoma II 0 20 LDLR LDLR MT1M MZT2B PAPPA PARK7 RNH1 RPS15 RPS9 SARAF STAR UBL5CHRDL1 Up 0.375920075 < 0.01 0.648042398 < 0.01 - - - 0 18 MT1M MT1M Up 0.29120842 < 0.01 0.335626478 < 0.01 - - - 0 0 MZT2B MZT2B Up 0.321545084 < 0.01 0.37077364 < 0.01 - - - 0 0 PAPPA PAPPA Up 0.785420042 < 0.01 0.358438843 < 0.01 - - - 0 0 PARK7 PARK7 Up 0.271376168 < 0.01 0.297989066 < 0.01 - - - 0 3 RNH1 RNH1 Up 0.259488668 < 0.01 0.273877319 < 0.01 - - - 0 0 RPS15 RPS15 Up 0.251464022 < 0.01 0.392164157 < 0.01 - - - 0 0 RPS9 RPS9 Up 0.309495896 < 0.01 0.405983232 < 0.01 - - - 0 0 SARAF SARAF Up 0.321075605 < 0.01 0.252252028 < 0.01 - - - 0 0 STAR STAR Up 1.183475152 < 0.01 0.259296741 < 0.01 - - - 0 0 UBL5 UBL5 Up 0.446544953 < 0.01 0.25479225 < 0.01 - - - 0 0 AHNAK AHNAK Down −0.416022103 < 0.01 −0.284696259 < 0.01 - - - 0 0 DSTN DSTN Down −0.267746826 < 0.01 −0.275020057 < 0.01 - - - 0 0 NIBAN1 NIBAN1 Down −0.363688859 < 0.01 −0.277190973 < 0.01 - - - 0 0 NREP NREP Down −0.305892285 < 0.01 −0.345200999 < 0.01 - - - 0 0 PXDN PXDN Down −0.414480896 < 0.01 −0.311726853 < 0.01 - - - 0 0 RBP1 RBP1 Down −0.337117044 < 0.01 −0.327552578 < 0.01 - - - 0 0 SPARC SPARC Down −0.370558256 < 0.01 −0.282401856 < 0.01 - - - 0 0 TXNIP TXNIP Down −0.660691169 < 0.01 −0.472651412 < 0.01 - - - 0 0 ATP5F1B  ATP Synthase F1 Subunit Beta, MSMO1  Methylsterol Monooxygenase 1, ATP5MG  ATP Synthase Membrane Subunit G, CD74  CD74 Molecule, CYP11A1  Cytochrome P450 Family 11 Subfamily A Member 1, FDXR  Ferredoxin Reductase, HLA-DRA  Major Histocompatibility Complex, Class II, DR Alpha, HLA-DRB1  Major Histocompatibility Complex, Class II, DR Beta 1, LDLR  Low Density Lipoprotein Receptor, MT1M  Metallothionein 1 M, MZT2B  Mitotic Spindle Organizing Protein 2B, PAPPA  Pappalysin 1, PARK7  Parkinsonism Associated Deglycase, RNH1  Ribonuclease/Angiogenin Inhibitor 1, RPS15  Ribosomal Protein S15, RPS9   Ribosomal Protein S9, SARAF  Store-Operated Calcium Entry Associated Regulatory Factor, STAR  Steroidogenic Acute Regulatory Protein, UBL5  Ubiquitin Like 5, AHNAK  AHNAK Nucleoprotein, DSTN  Destrin, Actin Depolymerizing Factor, NIBAN1  Niban Apoptosis Regulator 1, NREP  Neuronal Regeneration Related Protein, PXDN  Peroxidasin, RBP1  Retinol Binding Protein 1, SPARC  Secreted Protein Acidic And Cysteine Rich, TXNIP  Thioredoxin Interacting Protein a Target has at least 1 approved drug b Target has at least 1 ChEMBL compound with an activity cutoff of 80% Genetic information and findings for key proteins associated with PCOS LDLR MT1M MZT2B PAPPA PARK7 RNH1 RPS15 RPS9 SARAF STAR UBL5CHRDL1 ATP5F1B  ATP Synthase F1 Subunit Beta, MSMO1  Methylsterol Monooxygenase 1, ATP5MG  ATP Synthase Membrane Subunit G, CD74  CD74 Molecule, CYP11A1  Cytochrome P450 Family 11 Subfamily A Member 1, FDXR  Ferredoxin Reductase, HLA-DRA  Major Histocompatibility Complex, Class II, DR Alpha, HLA-DRB1  Major Histocompatibility Complex, Class II, DR Beta 1, LDLR  Low Density Lipoprotein Receptor, MT1M  Metallothionein 1 M, MZT2B  Mitotic Spindle Organizing Protein 2B, PAPPA  Pappalysin 1, PARK7  Parkinsonism Associated Deglycase, RNH1  Ribonuclease/Angiogenin Inhibitor 1, RPS15  Ribosomal Protein S15, RPS9   Ribosomal Protein S9, SARAF  Store-Operated Calcium Entry Associated Regulatory Factor, STAR  Steroidogenic Acute Regulatory Protein, UBL5  Ubiquitin Like 5, AHNAK  AHNAK Nucleoprotein, DSTN  Destrin, Actin Depolymerizing Factor, NIBAN1  Niban Apoptosis Regulator 1, NREP  Neuronal Regeneration Related Protein, PXDN  Peroxidasin, RBP1  Retinol Binding Protein 1, SPARC  Secreted Protein Acidic And Cysteine Rich, TXNIP  Thioredoxin Interacting Protein a Target has at least 1 approved drug b Target has at least 1 ChEMBL compound with an activity cutoff of 80%

The study design is shown in Fig.  1 . In brief, androgen-related hormone indices including TT, BIOT, DHEAS, and SHBG as exposure factors, female infertility associated with anovulation as the outcome factor, GWAS data from three publicly available datasets were utilized for a two-sample MR analysis in the initial study phase. A meta-analysis was performed to evaluate the combined causal effect with multiple MR results. We conducted Cochran’s Q test, Egger intercept tests, leave-one-out analysis, and funnel plots to assess the reliability of the MR estimates. Subsequently, Hetnet connectivity search, Bayesian colocalization analysis, single-cell type expression analysis, GSEA enrichment, single-cell gene signature (UCell) scoring, protein-protein interaction (PPI) network, and mediation analysis were utilized to characterize the function and speculate on the potential pathogenic pathways of these candidates associated with HA and PCOS. And then, the druggability of these therapeutic targets was evaluated. The P -values were adjusted using the Benjamini-Hochberg false discovery rate for multiple tests in our study. Fig. 1 Flowchart of the Mendelian randomization study and single-cell RNA sequencing analysis design. MR: Mendelian randomization; TT: total testosterone; BIOT: bioavailable testosterone; DHEAS: dehydroepiandrosterone sulfate; SHBG: sex hormone-binding globulin; scRNA-seq: single-cell RNA sequencing; GSEA: gene set enrichment analysis Flowchart of the Mendelian randomization study and single-cell RNA sequencing analysis design. MR: Mendelian randomization; TT: total testosterone; BIOT: bioavailable testosterone; DHEAS: dehydroepiandrosterone sulfate; SHBG: sex hormone-binding globulin; scRNA-seq: single-cell RNA sequencing; GSEA: gene set enrichment analysis The GWAS data of androgen-related indices were obtained from European participants in the UK Biobank by utilizing single-nucleotide polymorphism (SNP) information. Based on the data published by RUTH et al. [ 15 ] in 2020 and Leinonen et al. [ 16 ] in 2023, statistical clinical androgen indicators included SNP sites significantly related to TT, BIOT, DHEAS, and SHBG, involving 405,304, 346,596, 9722, and 377,351 samples respectively. Summary-level data for 117,098 cases of female infertility associated with anovulation from the FinnGen Consortium ( https://www.finngen.fi/ ). Exposure GWAS and outcome GWAS data were used from European populations to reduce the problem of population stratification. The detailed GWAS information of the exposure is displayed in Table S1-S2. Instrumental variables (IVs) need to satisfy three conditional assumptions: (1) Correlation hypothesis: IVs are highly associated with exposure factors; (2) Independence hypothesis: IVs cannot be correlated with confounding factors; (3) Exclusivity hypothesis: exposure factors are the only way that IVs affect outcome [ 17 ]. To satisfy the correlation hypothesis, the SNPs linked to exposure were identified at a genome-wide scale ( P < 5.0 × 10 − 8 ) [ 18 ]. Subsequently, the interference of linkage imbalance was excluded with a clustering threshold of R²10 were considered to be strong IVs [ 19 ]. Finally, outcome data were extracted according to expose-related SNPs in outcome GWAS, the exposure data and outcome data were integrated, and the palindromic SNP sequences were eliminated. We assessed the statistical power accounting for instrument strength (mean F-statistic >10), sample size, and the proportion of variance explained (R²) by the genetic instruments. Based on these parameters, our study had >80% power to detect an odds ratio (OR) of [1.20] at conventional significance thresholds (α = 0.05). Statistical analysis was conducted using R Studio (version 4.3.1, www.rstudio.com ) in conjunction with TwoSampleMR [ 20 , 21 ]. This research employed a variety of MR analysis techniques, such as inverse variance weighted (IVW) and MR-Egger regression, weighted median (WM), simple, and weighted mode methods. IVW was utilized to assess the primary impact of two-sample MR analyses, a weighted regression of SNP-outcome on SNP-exposure associations combined. We conducted sensitivity analyses using MR-Egger regression, presenting the results as odds ratios (ORs) and 95% confidence intervals (CIs). To assess the reliability of our findings, we evaluated the heterogeneity of individual SNPs with Cochran’s Q value. P values >0.05 were considered to indicate no significant heterogeneity or pleiotropy among the SNPs. In cases where there was notable heterogeneity in the MR results, we utilized the IVW random-effects model for correction [ 22 ]. Additionally, a leave-one-out test was performed to examine the impact of outlying and pleiotropic SNPs on causal estimates [ 23 ]. We conducted a meta-analysis to evaluate the combined relationship between androgen levels and infertility related to anovulation using MR results from both discovery and replication stages [ 24 , 25 ]. The choice of effect model was determined by the variability of results. When there was minimal significant heterogeneity (I2 ≤ 50%), we used the fixed-effects model for combining the results. In cases of substantial heterogeneity (I2 >50%), we utilized the random-effects model. The final causality was based on the meta-analysis results [ 26 ], unless only one reliable MR result was available, in which case it served as the basis for final causality determination. All statistical analyses were carried out using R language’s “meta” packages (version 4.3.1) [ 27 ]. Hetnet connectivity search integrates data from 29 public databases into a single, accessible network with a common format. This network consists of 47,031 nodes of various types and 2,250,197 edges of different types. Hetnet connectivity search can determine the relationship between two nodes in an unsupervised manner, capturing the semantic richness of edge type and providing results in both metapaths and paths [ 28 ]. The Coloc R package was employed for the analysis of colocalization and to ascertain whether the relationships between HA and anovulation-related infertility were impacted by linkage disequilibrium. The Bayesian method was utilized to assess five distinct hypotheses for each locus, yielding posterior probabilities for hypothesis testing (H0, H1, H2, H3, and H4). Colocalization between two traits in a specific region was deemed robust if the posterior probability of H4 (PPH4) was ≥ 0.8, with moderate support evidence defined as PPH4 less than 0.8 but greater than 0.5. Default parameters were employed for these analyses [ 29 ]. We downloaded the single-cell RNA sequencing (scRNA-seq) data of theca cells (TCs) in 5 PCOS and 5 non-PCOS samples from the website ( https://zenodo.org/records/7942968 ) and the scRNA-seq data of granulosa cells (GCs) in 3 PCOS and 3 non-PCOS samples from the Gene Expression Omnibus (GEO) database (Registration number: GSE240688 ). Seurat R package was used for preprocessing the downloaded matrix. Batch correction was performed using the Integration_SCP function from the SCP R package, with “Harmony” as the correction algorithm. Cell-type annotation was conducted manually. Differential expression analysis was carried out using the FindMarkers function (min. pct = 0, thresh. use = 0.01), and the resulting differentially expressed genes were utilized for identifying target genes and conducting gene set enrichment analysis (GSEA) [ 30 ]. GSEA was performed based on scRNA-seq data to analyze pathway enrichment using the Reactome database. Statistical significance for the GSEA analysis was determined with logistic regression adjusted P -value < 0.05 [ 31 ]. The “UCell” R package was employed to calculate gene feature enrichment scores in the single-cell dataset, with the corresponding gene sets sourced from the MsigDB database [ 32 ]. In order to evaluate the effectiveness of the proteins we identified, we conducted a search in the DrugBank database ( https://go.drugbank.com ) for target proteins associated with HA and PCOS [ 33 ]. We recorded details of the drugs and drug-gene interactions. The PPI network was established using the Search Tool for the Retrieval of Interacting Genes (STRING) database ( https://string-db.org ), with a minimum interaction confidence score set at 0.4 [ 34 ]. The findings from the PPI analysis were visualized using Cytoscape (v3.9.1).

To our knowledge, this is the first MR study to systematically evaluate the causal relationships between androgen-related biomarkers and female infertility associated with anovulation. Our findings suggest that SHBG appears to exert a protective effect on the incidence of anovulatory infertility. These results advance our understanding of the underlying mechanisms of these conditions and highlight potential therapeutic targets, which may facilitate early diagnosis and targeted interventions. Elevated levels of androgens in follicles have been demonstrated to stimulate the recruitment of small follicles, resulting in excessive recruitment that hinders the selection process of dominant follicles [ 35 , 36 ]. In luteinized GCs of PCOS patients, PGK1 and AR are highly expressed, and it is likely that an excessive PGK1-AR-MAP2K6 p38 axis regulates metabolic flux in GCs. This further stimulates glycolysis in GCs, disrupting the balance of the GC microenvironment and leading to dysfunctional ovulation [ 37 ]. Another study showed that excess exposure to androgens significantly increased expression levels of reactive oxygen species, Foxo1, Caspase3, and γ-H2AX while decreasing mitochondrial membrane potential in oocytes from PCOS mice. This resulted in more DNA damage and apoptosis [ 38 ]. A previous study found that the expression level of GLUT4 glucose transporter was significantly lower in endometrial epithelial cells of PCOS subjects with HA compared to the control group [ 39 ], suggesting that GLUT4 dysfunction may contribute to insulin resistance in PCOS [ 40 , 41 ]. Hyperinsulinemia, HA, and changes in paracrine signaling within follicles may disrupt activation, growth, and selection processes of follicles leading to impaired development and ovulation [ 42 ]. Furthermore, the presence of HA may impact the typical maturation of ovarian follicles through immunology-mediated pathways. Chen et al. discovered that the hyperandrogenic environment in the ovaries of PCOS patients can stimulate ovarian macrophages to M1 activation, potentially linked to changes in macrophage glucose metabolism [ 43 ]. Increased androgen levels can also raise chemokine levels within ovarian follicles, attracting chemokine receptor 1-expressing macrophages from inflammatory monocytes [ 44 ]. The heightened release of pro-inflammatory cytokines by macrophages in the hyperandrogenic environment could disrupt normal GC function, leading to premature luteinization and apoptosis, interfering with proper ovarian follicle growth and development, ultimately resulting in premature follicle activation and potential anovulation [ 43 , 44 ], aligning with our findings suggesting a potential causal relationship between elevated DHEAS and BIOT levels and anovulatory infertility. SHBG is a homodimeric polypeptide composed of amino acid residues, featuring a single binding site for sex hormones [ 9 ]. The binding of steroid to the SHBG-R SHBG complex leads to the activation of R SHBG , which in turn stimulates the production of cAMP. This subsequently initiates downstream androgen and estrogen signaling pathways [ 45 ]. SHBG can specifically bind to 80% of testosterone to form complexes to prevent it from adsorption and biochemical degradation by blood vessels, thus only free testosterone exhibits biological activity in the body [ 9 ]. Decreased serum SHBG level is one of the most significant features of PCOS patients [ 10 ]. Hyperinsulinemia and HA in PCOS patients inhibit hepatic synthesis of SHBG, resulting in decreased SHBG concentration and increased free testosterone, thereby impacting insulin metabolism and inducing insulin resistance. Consequently, this leads to compensatory overproduction of insulin and glucose metabolism disorder, ultimately forming a vicious cycle [ 10 ]. Some researches revealed that neurokinin 3 receptor antagonist fezolinetant [ 46 ] and SGLT1-2 inhibitor licogliflozin [ 47 ] could improve biochemical and clinical hyperandrogenism, as well as ovulatory dysfunction and fertility. Liu et al. found artemisinins directly targeted lon peptidase 1 (LONP1), enhanced LONP1-CYP11A1 interaction, and facilitated LONP1-catalyzed CYP11A1 degradation to block androgen overproduction [ 48 ]. It is a pity that repotrectinib, mitotane, and apolizumab are antineoplastic agents with substantial toxicity profiles, whose applicability in treating ovulation disorders is likely limited. However, our findings at least approve the druggability of potental androgen-related targets, and call for an in-depth investigation of safer, less toxic agents for facilitating effective therapeutic interventions of HA and related infertility in the future. Our findings highlight the role of androgens in ovulatory dysfunction, but infertility is a multifactorial condition influenced by diverse biological pathways. Emerging MR studies have provided valuable causal insights into additional mechanisms, including inflammation [ 49 ], obesity [ 50 ], lipid metabolism [ 51 , 52 ], reproductive endocrine disorders (including PCOS and endometriosis) [ 53 ], and genetic risk factors [ 54 ]. For instance, inflammatory proteins have been causally linked to infertility [ 49 ], suggesting potential crosstalk between androgen signaling and immune pathways in ovarian dysfunction. Similarly, dysregulated lipid metabolism—implicated in both female and male—may interact with androgen activity, as lipids are precursors for steroid hormone synthesis [ 51 , 52 ]. Notably, MEGF9 and MLLT11, identified as potential therapeutic targets for male infertility [ 54 ], underscore the importance of genetic validation in reproductive disorders. While our study focuses on androgens, these complementary findings encourage future work to explore integrative mechanisms (e.g., androgen-lipid-inflammation axes) in anovulatory infertility. Some researches revealed that neurokinin 3 receptor antagonist fezolinetant [ 46 ] and SGLT1-2 inhibitor licogliflozin [ 47 ] could improve biochemical and clinical hyperandrogenism, as well as ovulatory dysfunction and fertility. Liu et al. found artemisinins directly targeted lon peptidase 1 (LONP1), enhanced LONP1-CYP11A1 interaction, and facilitated LONP1-catalyzed CYP11A1 degradation to block androgen overproduction [ 48 ]. It is a pity that repotrectinib, mitotane, and apolizumab are antineoplastic agents with substantial toxicity profiles, whose applicability in treating ovulation disorders is likely limited. However, our findings at least approve the druggability of potental androgen-related targets, and call for an in-depth investigation of safer, less toxic agents for facilitating effective therapeutic interventions of HA and related infertility in the future. Our findings highlight the role of androgens in ovulatory dysfunction, but infertility is a multifactorial condition influenced by diverse biological pathways. Emerging MR studies have provided valuable causal insights into additional mechanisms, including inflammation [ 49 ], obesity [ 50 ], lipid metabolism [ 51 , 52 ], reproductive endocrine disorders (including PCOS and endometriosis) [ 53 ], and genetic risk factors [ 54 ]. For instance, inflammatory proteins have been causally linked to infertility [ 49 ], suggesting potential crosstalk between androgen signaling and immune pathways in ovarian dysfunction. Similarly, dysregulated lipid metabolism—implicated in both female and male—may interact with androgen activity, as lipids are precursors for steroid hormone synthesis [ 51 , 52 ]. Notably, MEGF9 and MLLT11, identified as potential therapeutic targets for male infertility [ 54 ], underscore the importance of genetic validation in reproductive disorders. While our study focuses on androgens, these complementary findings encourage future work to explore integrative mechanisms (e.g., androgen-lipid-inflammation axes) in anovulatory infertility. This research has several notable strengths. As far as we know, it is the first MR study investigating the causal relationship between androgen indices and anovulation-related infertility. We systematically evaluated multiple androgen indices using GWAS data from two independent sources, enabling a comprehensive risk assessment. Our integrative approach—combining MR with Hetnet connectivity, colocalization, enrichment, PPI analyses, and UCell scoring—enhanced the reliability of identified biomarkers and therapeutic targets for HA and PCOS. Notably, the discovered targets demonstrated druggable potential, laying the foundation for future research on safer, less toxic agents. However, there were still a few of limitations. Despite the MR approach being superior in causality inference, it was seldom possible to completely eliminate confounding bias or residual pleiotropy or heterogeneity. Our MR analysis was well-powered, but smaller power of DHEAS may require larger samples. Furthermore, as the participants in this study were of European descent, our findings limit generalizability and may not be applicable to other ethnic populations. As a result, large-scale epidemiological and clinical trials containing diverse populations are needed to support the above results, and we are currently planning to gradually conduct cell experiments and animal experiments to prove the druggability of other potential targets and search less or non-toxic promising drugs.

This comprehensive MR analysis indicated that SHBG may provide protection against anovulation-related female infertility. No exact causality between TT and ovulatory disorders suggested that TT may not be an optimal indicator for reflecting HA-related infertility. The druggabilities of CD74, CYP11A1, and HLA-DRB1 indicate potential targets for early detection of HA and shed light on future research on safer, less toxic agents for the treatment of anovulation-induced infertility associated with HA.

Infertility currently impacts around 15%−18% of couples in the reproductive age group globally [ 1 , 2 ], and its prevalence continues to rise steadily. It has been recognized as a public health priority, influenced by lifestyle and environmental factors, such as obesity and smoking, and even mental disorders [ 3 , 4 ]. Ovulatory disorders represent a major cause of infertility, accounting for about 25%~35% of female infertility cases [ 4 , 5 ]. There are numerous potential factors and contributors to ovulatory dysfunction that have complex effects on the quality of life for patients, resulting in disruptions in menstrual function, subfertility, and a range of hormone-related symptoms. It presents challenges for clinicians, trainees, educators, as well as those involved in basic, translational, clinical, and epidemiological research. Androgens are involved in folliculogenesis, but excessive androgen levels known as hyperandrogenism (HA) can lead to premature activation of follicles and an increase of granulosa cell apoptosis, and then cause anovulation and an increased risk of anovulation infertility [ 6 , 7 ]. About 70% of women with anovulation have polycystic ovary syndrome (PCOS) - characterized by ovulatory dysfunction, HA, and polycystic ovarian changes [ 4 , 8 ]. In clinical practice, serum total testosterone (TT) is commonly used to diagnose HA, however, the measurement of TT alone is not enough to accurately assess the presence of HA in women with PCOS, as only 20–30% of women with this condition have been found to have biochemical HA based solely on TT measurement [ 9 ]. While free testosterone is generally considered the biologically active form, sex hormone-binding globulin (SHBG) can specifically bind with testosterone to regulate its activity, avoid adsorption and biochemical degradation by vessels, and thus stabilize the concentration of testosterone [ 10 ]. Therefore, compared to TT, bioavailable testosterone (BIOT), dehydroepiandrosterone sulfate (DHEAS), and SHBG may accurately reflect androgen levels in the body. HA seems to be closely associated with ovulatory disorders, however, the causal relationship between androgens and anovulation has not been comprehensively assessed yet. Mendelian randomization (MR) analysis, utilizing genetic variation as an instrumental variable (IV), is a novel approach for investigating the causal relationship between hormones and diseases [ 11 , 12 ]. Numerous genome-wide association studies (GWAS) have identified multiple loci strongly linked to androgen and ovulatory dysfunction, offering the potential for MR analysis application [ 13 , 14 ]. In this study, we selected TT, BIOT, DHEAS, and SHBG as exposure factors as well as anovulation-related female infertility as disease outcomes. Using the two-sample MR approach, we comprehensively assessed the causal association between these androgen-related indices and ovulatory dysfunction. The results of this study will provide valuable evidence for the screening of anovulation-related female infertility and surveillance of ovulation disorders.

Supplementary Material 1 Supplementary Material 1