Proteome-wide Mendelian randomization and functional studies uncover therapeutic targets for polycystic ovarian syndrome

Study population and data source

The plasma pQTL data were retrieved from six large-scaled proteomic studies (Ferkingstad et al., 4,719 proteins¹¹; Sun et al., 2,994 proteins¹²; Sun et al., 2,923 proteins¹³; Suhre et al., 1,124 proteins¹⁴; Folkersen et al., 90 proteins¹⁵; and Zhang et al., 4,657 proteins¹⁶). Each study provided comprehensive summary-level statistics, and detailed characteristics of each study are available in Table S1.

For the primary MR analysis, summary statistics of PCOS were obtained from a meta-analysis of seven previously published PCOS GWASs, including 113,238 people (GWAS Catalog: GCST007089, n_{PCOS individual} = 10,074, n_Control = 103,164).¹⁷ For external validation, summary statistics of PCOS were obtained from the UK Biobank (UK Biobank: 20002_1350, n_{PCOS individual} = 436, n_Control = 193,717)¹⁸ and the FinnGen research project (FinnGen: E4_PCOS, R9 release, n_{PCOS individual} = 1,424, n_Control = 200,581).¹⁹ In cases where SNP or effect allele frequency information was missed in the GWAS summary statistics, we referred to the matched 1000 Genomes reference panel to complete the dataset.

Local study population recruited individuals diagnosed with PCOS or without PCOS from the Reproductive Medicine Center of Women’s Hospital, School of Medicine, Zhejiang University from September 2018 to June 2020. The clinical characteristics of the individuals are shown in Table S2. PCOS was diagnosed according to the revised Rotterdam criteria.²⁰ Women with any endocrine- or metabolism-related disorder were excluded, including congenital adrenal hyperplasia, androgen-secreting tumor, hypogonadotropic hypogonadism, thyroid dysfunction, Cushing syndrome, hyperprolactinemia, and premature ovarian failure. The control group comprised women without PCOS that attended the clinic owing to fallopian tube factors or male factors (including azoospermia, oligospermia, asthenozoospermia, and teratozoospermia) with regular menstrual menstruation, normal ovarian morphology, and no signs of hyperandrogenism and other systemic, gynecological or endocrine dysfunction. This study received approval from the Ethics Committee of Women’s Hospital, School of Medicine, Zhejiang University (no. IRB-20240097-R). The written informed consents were obtained from all participants.

Genetic instrument selection and primary mendelian randomization analysis

Genetic associations with plasma protein levels were extracted from the aforementioned six proteomic studies. Each protein was mapped to corresponding Uniprot ID using the bioDBnet database.²¹ Genomic coordinates for subsequent analyses were based on the human genome Build 37 (NCBI GRCh37) and processed using the “MungeSumstats” package.²² Instrumental variables (IVs) were selected based on the following criteria: (1) SNPs demonstrating genome-wide significant association with any protein were retained (p value <5 × 10⁻⁸), (2) SNPs located outside the major histocompatibility complex region were retained to avoid complex linkage disequilibrium (LD) structure, (3) LD clumping was conducted to identify independent pQTL (relative relaxed parameters [r2 < 0.1, clump_kb < 5,000] were used to enhance precision by incorporating more partial IVs^{9,23,24,25,26,27,28}), and (4) the R² and F-statistic (R² = 2 × EAF × [1-EAF] × beta²; F = R² × [N−2]/[1−R²]) were used to estimate the strength of IVs.²⁹ R² represented the proportion of variability in protein levels explained by each IV, while the F-statistic assesses the strength of the association between IVs and outcome. Only strong IVs with an F-statistic ≥10 were retained for each exposure of interest. In cases of duplicate proteins across studies, the protein with the highest sum of R² was selected. Additionally, IVs were classified as cis or trans pQTL based on the following criteria: cis pQTL was defined when the leading SNP was located within ±1 megabase (Mb) region of the transcription start site of the protein-coding gene, whereas trans pQTL was defined when the leading SNP laid outside of this region.¹² Ultimately, we identified a total of 74,693 cis pQTL and 2,590 unique plasma proteins as IVs for MR analysis.

We assessed the causality of circulating plasma proteins on PCOS by performing two-sample cis-MR using the “TwoSampleMR” package.⁸ In our study, the pQTL data served as IVs, and the GWAS of PCOS was utilized as the outcome trait data to investigate whether changes in plasma protein level have causal effects on PCOS. For proteins with only a single pQTL or with multiple pQTLs, we employed either the Wald ratio method or the inverse variance weighted MR (MR-IVW) to calculate MR estimates, respectively.³⁰ The Wald ratio estimates were calculated as the regression coefficient for genetic association with the outcome divided by the regression coefficient for genetic association with circulating protein levels. The IVW estimates were calculated as the average of instrument ratio coefficients weighted by the inverse variance. The estimates obtained from MR analyses approximated odds ratios for an increased risk of PCOS per standard deviation (SD) increase in normalized protein levels. To account for multiple testing, we used the false discovery rate (FDR) at ɑ = 0.05 based on Benjamini-Hochberg method to prioritize the results for further analysis. Simultaneously, the MR-Egger, weighted median, simple mode, and weighted mode methods were used as supplementary analyses to enhance the robustness of our findings.^31,32,33,34 There are partial sample overlaps between the proteomic study of Ferkingstad et al.¹¹ and PCOS study of GCST007089,¹⁷ as well as the proteomic study of Sun et al.¹³ and PCOS study of ukb_20002_1350.¹⁸ To assess the impact introduced by sample overlap, we evaluated the bias and type 1 error incidence for sample overlap using an online calculator.³⁵ Additionally, the MR-Pleiotropy Residual Sum and Outlier (MR-PRESSO) analysis detected and attempted to reduce horizontal pleiotropic SNPs and provided causal estimates after outlier correction.³⁶

Bayesian colocalization analysis

To assess whether associations between protein and PCOS risk were driven by the confounding of LD, we performed Bayesian colocalization analysis using the “coloc” package.³⁷ As described previously, Bayesian colocalization analysis assesses the probability that two traits share the same causal variant.³⁷ It provides the posterior probability for different hypotheses: (1) neither the protein nor PCOS is associated with any variants in the genomic locus (H0); (2) there is one causal variant associated with protein only (H1); (3) there is one causal variant associated with PCOS only (H2); (4) there are two distinct causal variants associated with protein and PCOS separately (H3); and (5) there is a shared causal variant associated with protein and PCOS (H4). For each protein, we included SNPs within ±500 kb of the pQTL. When a protein had multiple pQTLs, colocalization analysis was performed based on each pQTL, respectively. The pQTL with the strongest evidence for colocalization was shown. We used the coloc.abf algorithm with default parameters (p1 at 1 × 10⁻⁴, p2 at 1 × 10⁻⁴, and p12 at 1 × 10⁻⁵) and defined a gene as having evidence of colocalization if the posterior probability for shared causal variants (PH4) exceeded 80%.

SMR and heterogeneity in dependent instruments tests

The SMR and the HEIDI test occupied the ability to distinguish a pleiotropic model from a linkage model. We adopted the SMR and the HEIDI test to provide further investigation for the causal relationships identified through the primary MR analyses.³⁸ SMR analysis served as a supplementary method to examine for association between plasma protein level and PCOS, which may arise from a shared and potentially causal variant at a locus. The HEIDI test aimed to identify whether the observed association between proteins and PCOS was attributable to genetic LD of distinct causative variants.³⁹ A significance threshold of p value <0.05 was set for the SMR analysis. In the HEIDI test, a p value less than 0.05 showed that the association was driven by LD between two distinct genetic variants.

PPI network and druggability evaluation

We conducted PPI network analysis of potential targeted plasma proteins associated with PCOS risk (FDR-corrected p value <0.05 at MR discovery stage). Our aim was to explore the interactions between those proteins identified in MR analysis and the targets for current medications of PCOS. Firstly, we obtained drugs and their corresponding targets for PCOS from the Opentargets database.⁴⁰ Subsequently, we performed PPI network analysis using the Search Tool for the Retrieval of Interacting Genes (STRING) database.⁴¹ A minimum interaction score of 0.4 was set as the threshold for the PPI analysis.⁴² To assess the potential of identified proteins as potential therapeutic targets, we also extracted the druggable profile of candidate target proteins from the Drugbank.⁴³

Functional enrichment analysis for therapeutic targets for PCOS

To assess the potential role of identified therapeutic targets for PCOS, we classified their function using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. The potential therapeutic targets for PCOS in primary MR analysis (p value <0.05) were used as the input for the GO and KEGG enrichment analysis. Functional enrichment analysis relative to the overall GO pattern and the KEGG pathway was performed using “clusterprofiler” R package with Fisher’s exact test (FDR <0.05).⁴⁴

Sample collection

Ovarian GCs, follicular fluids, and peripheral blood samples were collected from both individuals with PCOS and controls. Follicular fluids were obtained from the largest follicle aspirated from each ovary, followed by centrifugation at 3,000 g for 10 min and subsequent storage at −80°C. GCs were isolated from the follicular fluid aspirated during oocyte retrieval by Ficoll density gradient centrifugation and then lysed in TRIzol for RNA extraction. Peripheral blood was collected on days 2～5 of a menstrual cycle. Venous blood was collected into heparin sodium tubes. The plasma was collected after centrifugation at 4°C and 3,000 g for 10 min and stored at −80°C for analysis.

Enzyme-linked immunosorbent assay

EGLN1 concentrations in human plasma and follicular fluid and mouse plasma were measured using enzyme-linked immunosorbent assay (ELISA) kits (MM-61907H1; MM-47886M1, Meimian, China). The ELISA procedure was conducted in accordance with the manufacturer’s recommendations.

Experimental animal care and ethics

Female Balb/c mice were housed in the laboratory animal center of Zhejiang University with free access to food and water and under a 12-h light/dark cycle. At 21 days of age, mice were randomly assigned to two groups as follows: the dehydroepiandrosterone (DHEA) group received subcutaneous injections of DHEA (HY-14650, USA) (6 mg/100 g body weight per day, dissolved in sesame oil) for 21 consecutive days to induce PCOS, while the control group received daily injections of the same volume of sesame oil. For roxadustat administration, mice were divided randomly into three groups from 21 days of age: the Oil + dimethyl sulfoxide (DMSO) group, the DHEA + DMSO group, and the DHEA + FG-4592 (FG) group. Roxadustat/FG-4592 (HY-13426, MCE) was dissolved in a mixed solvent containing 5% DMSO, 40% PEG-300, 5% Tween-80, and 50% N.S. Mice aged 21 days were intraperitoneally treated with roxadustat at a dosage of 10 mg/kg or an equal volume of vehicle once daily for 21 days as previously described.⁴⁵ The Animal Care and Use Committee of Zhejiang University of Medicine approved all care and experimental procedures for mice (ID: ZJU20240116).

Estrous cycle determination

From the 10^th day of modeling, vaginal smears were conducted daily at 09:00 for at least two estrous cycles. The stage of the estrus cycle was determined by microscopic examination of vaginal cytology using the Shorr method. Proestrus was characterized by round nucleated epithelial cells, estrus by keratinized squamous epithelial cells, metestrus by a mixture of epithelial cells and leukocytes, and diestrus by nucleated epithelial cells and a predominance of leukocytes.

Hematoxylin and eosin staining assay

The ovaries were initially fixed in 4% paraformaldehyde solution, then transitioned to 70% ethanol, and subsequently embedded in paraffin before being sectioned into 5 μm slices. Following deparaffinization, these sections were stained with hematoxylin and eosin (H&E). Observations, analysis, and photography of representative tissue sections were conducted under a light microscope.

Transmission electron microscopy

The ovaries were fixed with 2.5% glutaraldehyde solution in phosphate buffer overnight and washed with PBS three times for 10 min each. The ovaries were then immersed in 1% osmium tetroxide solution in phosphate buffer for 1 h and washed again with PBS three times for 10 min each. Thereafter, the ovaries were fixed in 2% uranyl acetate for 30 min; dehydrated with 50%, 70%, 90%, and finally 100% ethanol for 10 min each step; and rinsed twice with 100% acetone for 15 min each. Then, the ovaries were embedded in resin and sectioned. Observations and photography of the sections were carried out under an electron microscope (Thermo Scientific Talos L120C).

Protein isolation and western blot analysis

Cellular or tissue proteins were extracted by dissolving them in radioimmunoprecipitation assay (RIPA) lysis buffer containing protease and phosphatase inhibitors. Subsequently, total proteins were denatured in the 5× SDS-PAGE loading buffer by boiling at 95°C for 5 min. Equal amounts of proteins were separated using FuturePAGE 4%–20% gels (ET15412Gel, ACE Biotechnology) and then transferred onto nitrocellulose filter membranes. The membranes were blocked with 5% bovine serum albumin, incubated with primary antibodies overnight, and horseradish peroxidase-conjugated secondary antibodies at room temperature for an hour. Details regarding the antibodies and dilution ratio are listed in Table S3. The membrane was visualized using a BIO-RAD ChemiDoc MP Imaging System and analyzed with the Image Lab Software.

Quantitative polymerase chain reaction analysis

To measure the mRNA level of target genes, total RNA was extracted from cells or tissues using TRIzol reagent. Total RNA (1,000 ng) was reverse transcribed to complementary DNA (cDNA) using a DNase and reverse transcriptase reagent kit (R233, Vazyme). The quantitative polymerase chain reaction (qPCR) was carried out in a reaction mix containing SYBR Green (Q711-02, Vazyme) and cDNA on a StepOne Real-Time PCR System (Applied Biosystems) following the manufacturer’s instructions. To determine the relative expression levels of the target genes, the 2^−ΔΔCt method was used. Ct values were normalized to the RPS18 and Rps18 (Ribosomal Protein S18) genes, respectively. The primer sequences used in this study were listed in Table S4.

GCs extraction and culture

For the culture of primary GCs, mice were sacrificed 46 h after intraperitoneal injection of 5IU of pregnant mare serum gonadotropin, as previously described.⁴⁶ Ovarian antral follicles were punctured with a needle to release GCs. Cells were then cultured in DMEM/F12 supplemented with 5% fetal bovine serum and 1% penicillin and streptomycin at 37°C under a humidified atmosphere with 5% CO₂.

FerroOrange assay

Mouse GCs were seeded into a 12-well cell culture plate and incubated overnight for detection. Cells were washed in PBS once and stained in 1 μM Ferroorange (F374, DOJINGO, Japan) in DMEM/F12 for 30 min. Then, images were captured by a confocal microscope (Leica DMi8). The fluorescence intensity of each individual cell was quantified using ImageJ2 software (version 2.9.0).

Statistical analysis

Statistical analyses were conducted using GraphPad Prism 9 software. Each experiment was conducted with a minimum of three independent samples and repeated at least three times. All data are presented as the mean ± standard error of the mean (SEM) derived from the specific number of experiments. Prior to subsequent analyses, the Shapiro-Wilk normality test was utilized to assess data distribution. For two-group comparisons, statistical significance (p value <0.05) was determined using unpaired Student’s t tests for parametric data, and Mann-Whitney U test for non-parametric data. For multiple-group comparisons, statistical significance (p value <0.05) was evaluated by one-way analysis of variance (ANOVA) test for parametric data and Kruskal-Wallis test for non-parametric data.

Proteome-wide MR analysis identified 65 circulating proteins for PCOS

A total of 65 proteins exhibited significant associations with PCOS risk with FDR-corrected p value <0.05 using either the Wald ratio or IVW method (Figure 2; Table S5). Genetically predicted higher levels of 30 proteins were associated with an increased risk of PCOS, while the other 35 proteins were negatively associated with PCOS risk. In the replication stage, causal relationships for twelve proteins (FSHB, SIGLEC12, FUT10, NCR3LG1, NUCB2, C1QTNF9, EGLN1, MAPK9, SBSN, DBI, MRI1, and UROD) were successfully validated in independent datasets including the FinnGen or the UK Biobank dataset (p value <0.05) using the Wald ratio or IVW method (Figure 2; Table 1).^18,19 Out of the 65 proteins identified, the majority showed consistent associations between circulating protein levels and higher PCOS risk across various MR methods (including MR-Egger, weighted median, simple mode, and weighted mode^31,32,33,34) at both the discovery and replication stages, in comparison to the causalities obtained from the Wald ratio or IVW method (Table S6). In the MR-PRESSO analysis, there was insufficient evidence for horizontal pleiotropy in the association between the 65 circulating proteins and PCOS (global test p value >0.05; Table S7).

To assess potential bias introduced by sample overlap between the exposure and outcome, we performed evaluation of the type 1 error rate and biases as previous study described.³⁵ For each pair of causal association, the type 1 error rate was controlled under 0.05, and the biases due to sample overlap was less than 0.001 (Table S8). Hence, considerable weak instrument bias would not be expected, and the causalities between circulating proteins and PCOS were less likely to be biased by sample overlaps.

Colocalization analysis, SMR, and HEIDI tests supported the causality of circulating proteins with PCOS

Genetic colocalization analysis was conducted to mitigate potential confounding effect stemming from LD. Out of the 65 proteins identified as potential causal candidates through proteome-wide MR, 58 proteins exhibited evidence of genetic colocalization (PP4 > 80%) (Table S9). There was a substantial likelihood of a shared causal variant influencing both protein levels and PCOS risk.

To further assess the observed findings, we performed SMR and the HEIDI tests for the 65 proteins. Twenty-one of these proteins successfully passed the SMR test (p value <0.05) and aligned with the MR results in effect direction. Additionally, 63 proteins passed the HEIDI test (p value >0.05) (Table S10). Based on the above evidence, we categorized these proteins into three tiers. Four proteins (FSHB, SIGLEC12, FUT10, and NUCB2) successfully passed all tests and were classified into tier 1 (Table 1). Eight proteins that failed the colocalization analysis, SMR, or the HEIDI test (NCR3LG1, C1QTNF9, EGLN1, MAPK9, SBSN, DBI, MRI1, and UROD) were assigned to tier 2. The remaining 53 proteins, which failed in replication MR, were classified into tier 3.

PPI and druggability evaluation for candidate protein targets

To evaluate the druggability of candidate targets, the PPI network analysis revealed interactions between prioritized proteins (ADH7, TFRC, CTRB2, CPM, CXCL11, DBI, PER3, FSHB, HABP4, MAPK9, PSG5, CCN2, RET, FABP2, IDI2, and POSTN) and the targets of current medications for PCOS (Figure 3; Tables S11 and S12). Utilizing STRING database,⁴¹ several interactions were determined as the most reliable interactions (known interactions), including FSHR-FSHB, NDUFS2-POSTN, ESR1-CCN2, DPP4-CXCL11, MT-ND4-POSTN, NDUFS8-POSTN, MT-ND1-POSTN, MT-ND3-POSTN, and NDUFS3-POSTN (Table S12). Specifically, FSHB was found to interact with FSHR, the target of LH/CG receptor agonist (menotropins) and FSHR agonist (follitropin alfa). CCN2 was associated with ESR1, the target of estrogen receptor alpha modulators (tamoxifen, enclomiphene, and clomiphene) and estrogen receptor alpha agonist (ethinyl estradiol and estradiol). CXCL11 was associated with DPP4, the target of dipeptidyl peptidase IV inhibitor (sitagliptin, saxagliptin, and linagliptin). Additionally, STRING revealed that POSTN interacted with a series of mitochondrially encoded NADH dehydrogenase (MT-ND1, MT-ND3, and MT-ND4) and NADH dehydrogenase (ubiquinone) Fe-S Protein (NDUFS2, NDUFS3, and NDUFS8).⁴¹ MT-ND1, MT-ND3, MT-ND4, NDUFS3, and NDUFS8 are all targets of metformin.

To evaluate drug development activities of candidate targets, we searched for current medications targeting identified potential causal proteins in the Drugbank database⁴³ (Table S13). We identified 80 medications targeting 20 proteins that may modify the disease, including antioxidant agent, estrogenic steroid, or drugs treating hypoestrogenism, acne, and autoimmune diseases. Overall, therapeutics targeting 65 proteins requires further evaluation among which 20 are druggable, and 45 represent future breakthrough points.

Validation and functional enrichment analysis for therapeutic targets of PCOS

To assess the impact of tier 1 and 2 therapeutic targets for PCOS on human follicle development, we analyzed the mRNA expression changes in human GCs. In comparison to the control group, PCOS individuals exhibited decreasing mRNA levels of FUT10 (MIM: 616931), NCR3LG1 (MIM: 613714), NUCB2 (MIM: 608020), and DBI (MIM: 125950), as well as increasing mRNA levels of EGLN1 (MIM: 606425) and MAPK9 (MIM: 602896) in GCs (Figure 4A). However, there were no significant changes in mRNA levels of SIGLEC12 (MIM: 606094), C1QTNF9 (MIM: 614285), SBSN (MIM: 609969), MRI1 (MIM: 615105), and UROD (MIM: 613521) in GCs between control group and PCOS individuals (Figure 4A). Furthermore, we utilized a DHEA-induced PCOS mouse model to detect the expression of therapeutic targets for PCOS. Compared with female control mice, DHEA mice exhibited significantly reduced mRNA levels of Fut10, Nucb2, C1qtnf9, and Dbi, as well as elevated mRNA levels of Egln1, Mapk9, Mri1, and Urod (Figure 4B). To determine the effect of therapeutic targets on the occurrence of PCOS, all potential therapeutic targets for PCOS were subjected to GO enrichment and KEGG pathway analysis to explore their putative biological function (Figures S1–S4). The KEGG pathway analysis identified the hypoxia-inducible factor 1 pathway among the various pathways related to the therapeutic targets for PCOS (Figure S4). EGLN1, identified as a tier 2 therapeutic target protein for PCOS, is known as prolyl hydroxylase domain-containing protein. In the presence of oxygen, EGLN1 hydroxylates HIF1α, resulting in a binding site for a ubiquitin ligase complex that includes von Hippel-Lindau (VHL), leading to HIF1α degradation.⁴⁷ Based on the above analyses, there was substantial evidence supporting EGLN1 in the MR analysis and colocalization analysis, along with its consistent expression patterns observed in GCs of individuals with PCOS and ovaries of DHEA-induced PCOS mouse model.

Roxadustat improves the onset of PCOS through targeting the EGLN1-HIF1α-ferroptosis axis

In both the follicular fluids and plasma of PCOS individuals, EGLN1 protein levels were significantly elevated compared with control group (Figures 5A and 5B). DHEA-induced PCOS mice also exhibited significantly increased plasma levels of EGLN1 compared with oil mouse controls (Figure 5C). Consequently, we delved further into the role of EGLN1-HIF1α signaling in the pathogenesis of PCOS and the potential of Roxadustat (FG-4592, an approving EGLN1 inhibitor for treating renal anemia) in intervening in PCOS. Compared with the control group (Oil + DMSO), the DHEA-induced PCOS mouse model (DHEA + DMSO) displayed typical polycystic ovaries and notable estrous acyclicity. After intraperitoneal administration of roxadustat, the DHEA + FG group reversed the disrupted estrous cycle and morphological changes in the ovary (Figures 5D–5F). Roxadustat-induced stabilization of HIF1α has been found to mitigate acute kidney injury by decreasing ferroptosis.⁴⁸ Recent studies revealed that elevated miR-93-5p in GCs of PCOS individuals promoted ferroptosis, whereas metformin can inhibit ovarian ferroptosis through the SIRT3-AMPK-mTOR pathway to improve ovarian function in a PCOS mouse model.^49,50 Transmission electron microscopy (TEM) was utilized to observe mitochondria ultrastructure in mouse GC cells before and after roxadustat administration (Figure 5G). Compared to the Oil + DMSO mice, we observed increased mitochondrial membrane density, reduced mitochondrial crista structures, and broken mitochondrial outer membrane in DHEA + DMSO mice, which was characteristic morphological features of ferroptosis reported previously.⁵¹ Relatively normal mitochondria ultrastructure in mouse GC cells was observed in the DHEA + FG group compared to the DHEA + DMSO group (Figure 5G). We also used the FerroOrange probe to detect intracellular Fe²⁺ levels in GCs. We observed enhanced intracellular Fe²⁺ levels in GCs of the DHEA + DMSO group and lowered Fe²⁺ levels after roxadustat treatment in the DHEA + FG group (Figure 5H).

Next, we detected the expression of ferroptosis-related genes in ovarian tissue. Key genes involved in the positive regulation of ferroptosis, including Acsl4, Tfrc, Alox5, Alox12, Alox15, and Lpcat3, were markedly upregulated in the ovaries of the DHEA + DMSO group compared with those in the Oil + DMSO group (Figure 5I). Western blot analysis also showed higher protein accumulation of TFRC and ACSL4, as well as lower protein accumulation of GPX4 in the ovaries of the DHEA + DMSO group (Figure 5J). After administration with roxadustat, the expression of ferroptosis positive regulation genes (Acsl4, Tfrc, Alox5, Alox12, Alox15, Aloxe5, and Lpcat3) was lower than that in the DHEA + DMSO group (Figure 5I). Lower levels of TFRC and ACSL4, as well as higher level of GPX4, were also observed in the ovaries of the DHEA+FG group according to Western blot analysis (Figure 5J). Additionally, Roxadustat administration restored protein accumulation of HIF1α in the ovaries of the DHEA + FG group. Moreover, the mRNA expression of Fabp3 (fatty acid binding protein 3), a key HIF1α target genes responsible for fatty acid uptake and lipid storage,⁵² was restored by HIF1A activation (Figure 5K). Collectively, these findings confirmed that roxadustat treatment allows for HIF1α stabilization and antagonizes ferroptosis in the ovaries of the DHEA-induced mouse model.

Principal findings

The current study systematically investigated circulating proteins associated with PCOS risk, particularly focusing on the impact of EGLN1 on onset of PCOS using an MR-Bayesian model averaging causal analyses, biological pathway enrichment analysis, and experimental approaches. A comprehensive pipeline comprising analytical techniques revealed causal associations between the decreased level of plasma proteins (FSHB, SIGLEC12, FUT10, NCR3LG1, NUCB2, and C1QTNF9) and a higher risk of PCOS. Causal associations were found between the elevation of plasma proteins (EGLN1, MAPK9, SBSN, DBI, MRI1, and UROD) and a higher risk of PCOS. In PPI analyses, FSHB, POSTN, CCN2, and CXCL11 ranked among the top proteins interacting with the targets of current PCOS medications. A list of 80 current drugs and 20 druggable proteins further supported the potential for drug development programs targeting PCOS. Additionally, KEGG functional enrichment analysis for therapeutic targets revealed involvement of the HIF1α pathway in PCOS. Using a DHEA-induced PCOS mouse model, we explored the EGLN1-HIF1α-ferroptosis axis in inducing PCOS and investigated the therapeutic potential of the EGLN1-targeting roxadustat.

Concordance of observational and causal associations for identified proteins

One of the key strengths of this study lies in comparing evidence between observational and MR analyses for candidate proteins. Given the heritability of PCOS, the MR method has recently gained traction in assessing causal associations between PCOS and etiological risk factors.⁵³ Utilizing MR approaches, we have discovered drug targets having a causal association with an increased risk of PCOS. Several strategies were employed to enhance the biological credibility of these candidate proteins. Through our experimental approaches, we have partially confirmed the consistency between the expression levels of tier 1 and 2 targets (FUT10, NUCB2, EGLN1, and MAPK9) in ovaries and GCs, as well as the effect size of these targets in MR analyses. FUT10, a member of the α-1,3 fucosyltransferase family, mediates the α-1,3 fucosylation level of TNFR1, enhances its binding ability to TNF-α, and eventually accelerates the osteoarthritis process.⁵⁴ In the ventricular zone of the embryonic brain, FUT10 activity is necessary to maintain the self-renewal of neurosphere-forming neural stem cells.⁵⁵ While the regulation role of FUT10 in PCOS remains largely unknown, it is worth noting that core fucosylation catalyzed by fucosyltransferase 8 in GCs controls oocyte and follicular development via the FSH/FSHR pathway, loss of which contributes to poor ovary responders.⁵⁶ Further investigation into the associations of fucosylation catalyzed by FUT10 and PCOS occurrence is warranted. NUCB2, initially identified in the nucleus of the hypothalamus, plays a role in diet and energy homeostasis.⁵⁷ In a letrozole-induced rat PCOS model, plasma NUCB2 level was decreased compared to the control group, consistent with our findings.⁵⁸ MAPK9, a member of the MAP kinase family, phosphorylated BCL2L11 upon activation, leading to its dissociation from BECN1 and subsequently participating in RIPK1-mediated activation of autophagy.⁵⁹ Cultured GCs from women with PCOS exhibited increased levels of autophagy markers (ATG7, ATG5, BECN1, and LC3-II:LC3-I).⁶⁰ Moreover, the mRNA level of MAPK9 was markedly elevated in human PCOS ovarian tissue compared to normal ovarian tissue,⁶¹ suggesting its involvement in autophagy activation in PCOS.

FSHB encodes the beta subunit of follicle-stimulating hormone, which is secreted from the pituitary gland. A recent study has revealed that FSH promotes glutamine synthesis in GCs, thereby maintaining follicle wall integrity and inhibiting GC apoptosis, which is crucial for initiating ovulation.⁶² Inadequate FSH levels are strongly associated with anovulation in PCOS. This finding aligns with our MR analysis, which suggests that FSHB is a protective factor for PCOS. In a GWAS of PCOS, rs11031006, an independent signal reach genome-wide significance for association with PCOS, is located near FSHB. The PCOS-susceptibility allele at rs11031006 is also robustly associated with lower circulating FSH concentrations.^63,64 The allele rs11031006-G has been associated with higher serum FSH levels and a lower risk of PCOS.⁶⁵

Implications for therapeutics targeting EGLN1-HIF1α-ferroptosis axis inducing PCOS

EGLN1 occupies hydroxylase activity to allow ubiquitin-mediated proteasomal degradation of HIF1α via VHL binding under physiological conditions.⁶⁶ HIF1α is involved in glycolysis and lipid metabolism and regulates ferroptosis with varying effects among cell lines and animal models. In mouse models of diabetic renal tubular injury and testicular injury, HIF1α induces ferroptosis, resulting in damage to renal tubules, testicular Leydig, and Sertoli cells.^67,68 In a cancer cell model and model of liver fibrosis, HIF1α antagonizes ferroptosis by upregulating lipid metabolism-related genes, including FABP3 and FABP7.^69,70 In our DHEA-induced PCOS mouse model, we found that increased levels of ELGN1 were negatively associated with HIF1α stabilization, inducing ferroptosis of GCs in the ovary. Roxadustat, an EGLN1 inhibitor developed for treating anemia induced by chronic kidney disease or myelodysplastic syndromes, reversibly binds to and inhibits EGLN1.⁷¹ Roxadustat reverses the degradation of transcription factors in the HIF family under normal oxygen conditions. The role of roxadustat is cell background dependent. In some contexts, roxadustat amplifies HIF signaling and suppresses chemoresistant glioblastoma growth via inducing ferroptosis.⁷² Conversely, in other instances, HIF-1α stabilization by roxadustat has been found to mitigate acute kidney injury by decreasing ferroptosis.⁴⁸ Additionally, roxadustat can also serve as a protective agent to reduce organ fibrosis and inflammation, correct metabolic disorders, reduce oxidative stress, improve mitochondrial function, and normalize angiogenesis.⁷³ In our study, roxadustat improved the onset of PCOS by targeting the EGLN1-HIF1α-ferroptosis axis and reversed the disrupted estrous cycle and ovarian morphology. These findings suggest a considerable therapeutic potential for roxadustat in PCOS, providing an experimental basis before expanding its indication.

Strength and limitations

The strength of our investigation lies in the utilization of proteome-wide MR analyses to systematically estimate the causal effects of plasma proteins on PCOS. The MR design offers several advantages, including large-scale sample sizes, comprehensive proteome coverage, and the ability to minimize bias due to confounding and reverse causality. The replication MR analyses conducted across multiple datasets further validated the robustness of our findings. Additionally, colocalization analysis proved to be a powerful tool in revealing the pleiotropic effects of certain loci on multiple traits. Supplementary evidence from PPI networks, druggability evaluations, and experimental research provided insights into the potential pathogenic effect of druggable targets. Through this approach, we have discovered the HIF1α-ferroptosis axis as a relevant pathogenic pathway of PCOS in which the elevated level of EGLN1 plays a crucial role.

There are still several limitations in our analysis that merit consideration. Firstly, the plasma proteome data were obtained from European cohorts, which may limit the generalization of our conclusions to other populations.^74,75 Secondly, there is a small percentage of sample overlap, which may result in potential bias and type 1 error inflation due to weak instrument bias. The F statistic quantifies the strength of the association between IVs and exposures. In our preliminary genetic instrument selection, we retained only IVs with an F statistic ≥10 to avoid weak instrument bias. We also calculated the type 1 error rate and found the estimated bias from current sample overlap rate was less than 0.1%. Recent studies have demonstrated the reliability and validity of two-sample MR using overlapping samples, under the premise of low expected analytic bias and minimal type 1 error.^76,77,78,79 Therefore, the current degree of sample overlap does not significantly affect the robustness of causalities in our finding. Thirdly, in our experimental research, we detected expression levels of identified drug targets in GCs of PCOS but didn’t estimate their levels in other ovarian cells or organ tissues. Assessing protein levels from other ovarian cells or organ tissues in PCOS may provide further insights into pathogenesis and complication of PCOS. Fourthly, for candidate targets consistent across the observational validation and PCOS GWASs, our interpretation and hypothesis are limited by the lack of detailed understanding of pathogenesis and phenotypes of ovarian structure and function. Moreover, the MR analysis was based on the impact of circulating proteins, including cell-secreted and leaked proteins. The abundance of circulating proteins may differ from the abundance within cells and tissues; therefore, potential effects of cell- or tissue-specific protein abundance are not directly explored in our study.

Conclusions

In summary, this study identified numerous associations of circulating proteins with the risk of PCOS using integrated genetic approaches and experimental validation. FSHB, SIGLEC12, FUT10, NCR3LG1, NUCB2, C1QTNF9, EGLN1, MAPK9, SBSN, DBI, MRI1, and UROD were prioritized as potential drug targets for PCOS, pending verification in future trials. Building upon this, the EGLN1-HIF1α-ferroptosis axis was found to potentially mediate the pathogenesis of PCOS. Our findings underscore the value of population-scale genomic and proteomic analyses in identifying critical loci and pathways with causal implications for PCOS.