Distinct subtypes of polycystic ovary syndrome with novel genetic associations: An unsupervised, phenotypic clustering analysis

Summary

Background

Polycystic ovary syndrome (PCOS) is a common, complex genetic disorder affecting up to 15% of reproductive-age women worldwide, depending on the diagnostic criteria applied. These diagnostic criteria are based on expert opinion and have been the subject of considerable controversy. The phenotypic variation observed in PCOS is suggestive of an underlying genetic heterogeneity, but a recent meta-analysis of European ancestry PCOS cases found that the genetic architecture of PCOS defined by different diagnostic criteria was generally similar, suggesting that the criteria do not identify biologically distinct disease subtypes. We performed this study to test the hypothesis that there are biologically relevant subtypes of PCOS.

Methods and findings

Using biochemical and genotype data from a previously published PCOS genome-wide association study (GWAS), we investigated whether there were reproducible phenotypic subtypes of PCOS with subtype-specific genetic associations. Unsupervised hierarchical cluster analysis was performed on quantitative anthropometric, reproductive, and metabolic traits in a genotyped cohort of 893 PCOS cases (median and interquartile range [IQR]: age = 28 [25–32], body mass index [BMI] = 35.4 [28.2–41.5]). The clusters were replicated in an independent, ungenotyped cohort of 263 PCOS cases (median and IQR: age = 28 [24–33], BMI = 35.7 [28.4–42.3]). The clustering revealed 2 distinct PCOS subtypes: a “reproductive” group (21%–23%), characterized by higher luteinizing hormone (LH) and sex hormone binding globulin (SHBG) levels with relatively low BMI and insulin levels, and a “metabolic” group (37%–39%), characterized by higher BMI, glucose, and insulin levels with lower SHBG and LH levels. We performed a GWAS on the genotyped cohort, limiting the cases to either the reproductive or metabolic subtypes. We identified alleles in 4 loci that were associated with the reproductive subtype at genome-wide significance (PRDM2/KAZN, P = 2.2 × 10⁻¹⁰; IQCA1, P = 2.8 × 10⁻⁹; BMPR1B/UNC5C, P = 9.7 × 10⁻⁹; CDH10, P = 1.2 × 10⁻⁸) and one locus that was significantly associated with the metabolic subtype (KCNH7/FIGN, P = 1.0 × 10⁻⁸). We developed a predictive model to classify a separate, family-based cohort of 73 women with PCOS (median and IQR: age = 28 [25–33], BMI = 34.3 [27.8–42.3]) and found that the subtypes tended to cluster in families and that carriers of previously reported rare variants in DENND1A, a gene that regulates androgen biosynthesis, were significantly more likely to have the reproductive subtype of PCOS. Limitations of our study were that only PCOS cases of European ancestry diagnosed by National Institutes of Health (NIH) criteria were included, the sample sizes for the subtype GWAS were small, and the GWAS findings were not replicated.

Conclusions

In conclusion, we have found reproducible reproductive and metabolic subtypes of PCOS. Furthermore, these subtypes were associated with novel, to our knowledge, susceptibility loci. Our results suggest that these subtypes are biologically relevant because they appear to have distinct genetic architecture. This study demonstrates how phenotypic subtyping can be used to gain additional insights from GWAS data.

Keywords:

Genetic loci – Phenotypes – Molecular genetics – Genome-wide association studies – Human genetics – insulin – Polycystic ovary syndrome – Quantitative traits

Introduction

Understanding the genetic architecture of complex diseases is a central challenge in human genetics [1–3]. Often defined according to arbitrary diagnostic criteria, complex diseases can represent the phenotypic convergence of numerous genetic etiologies under the same clinical diagnosis [4–8]. Recent studies in type 2 diabetes (T2D) support the concept that there are disease subtypes with distinct genetic architecture [7,8]. Identifying and addressing genetic heterogeneity in complex diseases could increase power to detect causal variants and improve treatment efficacy [9].

Polycystic ovary syndrome (PCOS) is a highly heritable, complex genetic disorder affecting up to 15% of reproductive-age women worldwide, depending on the diagnostic criteria applied [10]. It is characterized by a variable constellation of reproductive and metabolic abnormalities [11–13]. It is the leading cause of anovulatory infertility and a major risk factor for T2D in young women [14]. Despite these substantial morbidities, the etiology (or etiologies) of PCOS remains unknown [15]. Accordingly, the commonly used diagnostic criteria for PCOS, the National Institutes of Health (NIH) criteria [16] and the Rotterdam criteria [17,18], are based on expert opinion rather than mechanistic insights and are designed to account for the diverse phenotypic presentations of PCOS. The NIH criteria require the presence of hyperandrogenism (HA) and chronic oligo/anovulation or ovulatory dysfunction (OD) [16]. The Rotterdam criteria include polycystic ovarian morphology (PCOM) and require the presence of at least 2 of these 3 key reproductive traits, resulting in 3 different affected phenotypes: HA and OD with or without PCOM, also known as NIH PCOS, as well as 2 additional non-NIH Rotterdam phenotypes, HA and PCOM and OD and PCOM.

Genome-wide association studies (GWAS) have considerably advanced our understanding of the pathophysiology of PCOS. These studies have implicated gonadotropin secretion [19] and action [20,21], androgen biosynthesis [20–22], metabolic regulation [22,23], and ovarian aging [23] in PCOS pathogenesis. A recent meta-analysis [22] of GWAS was the first study to investigate the genetic architecture of the diagnostic criteria. Only one of 14 PCOS susceptibility loci identified was significantly more strongly associated with the NIH phenotype compared to non-NIH Rotterdam phenotypes or to self-reported PCOS. These findings suggested that the genetic architecture of the phenotypes defined by the different PCOS diagnostic criteria were generally similar. Therefore, the current diagnostic criteria do not appear to identify genetically distinct disease subtypes.

It is possible to identify physiologically relevant complex disease subtypes through cluster analysis of phenotypic traits [8,24,25]. Indeed, there have been previous efforts to subtype PCOS using unsupervised cluster analysis of its hormonal and anthropometric traits [26–29]. However, there has been no validation that the resulting PCOS subtypes were biologically meaningful by testing their association with genetic variants, with other independent biomarkers, or with outcomes such as therapeutic responses. In this study, we sought to 1) identify phenotypic subtypes of PCOS using an unsupervised clustering approach on reproductive and metabolic quantitative traits from a large cohort of women with PCOS, 2) validate those subtypes in an independent cohort, and 3) test whether the subtypes thus identified were associated with distinct common genetic variants. As an additional validation, we investigated the association of the subtypes with rare genetic variants we recently identified in a family-based PCOS cohort [30].

Methods

Subjects

This study used biochemical and genotype data from our previously published PCOS GWAS, Hayes and Urbanek and colleagues [19], in which a discovery sample (Stage 1) of 984 unrelated PCOS cases and 2,964 population controls was studied, followed by a replication sample (Stage 2) of 1,799 PCOS cases and 1,231 phenotyped reproductively normal control women. All cases were of European ancestry. The present study began as an exploratory analysis to test the hypothesis that subtypes existed within the PCOS GWAS cohorts. Further analyses were performed once subtypes were identified. This study is reported according to the Strengthening the Reporting of Genetic Association Studies (STREGA) guideline (S1 Checklist). The study was approved by the Institutional Review Board of Northwestern University Feinberg School of Medicine, and each subject provided written informed consent prior to the study [19].

All PCOS cases were aged 13–45 years and were diagnosed according to the NIH criteria [10] of hyperandrogenism and chronic anovulation (8 or fewer menses per year), excluding specific disorders of the adrenal glands, ovaries, or pituitary gland [31]. Cases fulfilling the NIH criteria also meet the Rotterdam criteria for PCOS [10]. The GWAS cohorts included in the cluster analysis, the PCOS Family Study and Pregnancy in PCOS II (PPCOSII) study [19] (S1 Table), had complete data for the following traits: body mass index (BMI), testosterone (T), sex hormone binding globulin (SHBG), dehydroepiandrosterone sulfate (DHEAS), luteinizing hormone (LH), follicle-stimulating hormone (FSH), fasting insulin (Ins0), and fasting glucose (Glu0). Complete data for these quantitative traits were not available in the other GWAS cohorts because of differences in phenotyping protocols [19] (S1 Table). Two additional NIH PCOS cohorts with complete quantitative trait data were included in the present study. An ungenotyped cohort of 263 cases was used for clustering replication. A family-based whole-genome sequencing cohort of 73 PCOS cases was investigated to assess subtype clustering in families and for rare variant analysis [30].

Contraceptive steroids had been stopped at least 3 months prior to screening for the PCOS Family Study, ungenotyped, and whole-genome sequencing PCOS cohorts. Elevated T, non-SHBG bound T, and/or DHEAS levels were documented in all PCOS cases prior to enrollment in these cohorts. PPCOSII was a randomized clinical trial of letrozole versus clomiphene citrate for infertility in PCOS [32]. The PCOS cases in this study had contraceptive steroids discontinued at least 2 months prior to their baseline phenotyping visit. Since the PCOS women in this trial were seeking fertility, the majority were not on recent contraceptive steroid therapy. Therefore, it is unlikely that recent contraceptive steroid use altered T or SHBG levels in the PCOS cases included in the cluster analysis.

All subjects included in the cluster analysis were from US-based study sites. The GWAS Stage 2 replication included 2 cohorts from Europe in addition to US-based cohorts [19]. Neither European cohort was included in the cluster analysis because of incomplete quantitative trait data. We compared age and BMI in the cohorts included in the cluster analysis of cases with complete quantitative trait data versus cases from the same cohort not included because of missing data. There were no significant differences in these parameters, suggesting that the included cases were similar to those excluded because of missing data (S2 Table).

Population-based control DNA samples for the GWAS Stage 1 sample were obtained from the NUgene biobank [33] from women of European ancestry, aged 18–97 years. Control women in the Stage 2 sample were phenotyped reproductively normal women of European ancestry, aged 15–45 years, with regular menses and normal T levels, and who were not receiving contraceptive steroids for at least 3 months prior to study [19]. T, DHEAS, SHBG, LH, FSH, Glu0, and Ins0 levels were measured as previously reported [19].

Clustering

Clustering was performed in PCOS cases on 8 adjusted quantitative traits: BMI, T, DHEAS, Ins0, Glu0, SHBG, LH, and FSH. There were 893 combined cases from the GWAS samples with complete quantitative trait data available for clustering. Quantitative trait values were first log_e-normalized and adjusted for age and assay method, which varied according to the different study sites where samples were collected [19], using a linear regression. An inverse normal transformation was then applied for each trait to ensure equal scaling. The normalized trait residuals were clustered using unsupervised, agglomerative, hierarchical clustering according to a generalization of Ward’s minimum variance method [34,35] on Manhattan distances between trait values. Differences in adjusted, normalized trait values between subtypes were assessed using Kruskal–Wallis and unpaired Wilcoxon rank–sum tests corrected for multiple testing (Bonferroni). Cluster stability was assessed by computing the mean Jaccard coefficient from a repeated nonparametric bootstrap resampling (n = 1,000) of the dissimilarity matrix [36]. Jaccard coefficients below 0.5 indicate that a cluster does not capture any discernable pattern within the data, while a mean coefficient above 0.6 indicates that the cluster reflects a real pattern within the data [36]. Cluster reproducibility was further assessed by repeating the clustering procedure in an independent cohort of 263 PCOS cases.

Association testing

Stage 1 samples were genotyped using the Illumina OmniExpress (HumanOmniExpress-12v1_C; San Diego, CA, USA) single nucleotide polymorphism (SNP) array. Stage 2 samples were genotyped using the Metabochip [37] with added custom variant content based on ancillary studies and the discovery results [19]. The Stage 2 association replication in this study was therefore limited; many of the loci from Stage 1 were therefore not characterized in Stage 2. Low-quality genotypic data were removed as described previously [19]. SNPs were filtered according to minor allele frequency (MAF ≥ 0.01), Hardy–Weinberg equilibrium (P ≥ 1 × 10⁻⁶), call rate (≥0.99), minor allele count (MAC > 5), mendelian concordance, and duplicate sample concordance. Only autosomal SNPs were considered. Ancestry was evaluated using a principal component analysis (PCA) [38] on 76,602 linkage disequilibrium (LD)-pruned SNPs [19]. Samples with values >3 standard deviations from the median for either of the first 2 principal components (PCs) were excluded (34 in discovery; 37 in replication). Genotype data were phased using ShapeIT (v2.r790) [39] and then imputed to the 1000 Genomes reference panel (Phase3 v5) [40] using Minimac3 via the Michigan Imputation Server [41]. Imputed SNPs with an allelic r² below 0.8 were removed from analysis [42].

Association testing was performed separately for Stage 1 and Stage 2 samples. Of the 893 combined cases from both stages included in the clustering analysis, 555 were from the Stage 1 sample, and 338 were from the Stage 2 sample. In Stage 1, 2,964 normal controls were used, and 1,134 were used in Stage 2. Logistic regressions were performed using SNPTEST [43] for case–control status under an additive genetic model, adjusting for BMI and first 3 PCs of ancestry. P-values are reported as P₁ and P₂ for Stage 1 and Stage 2, respectively. Cases were limited to specific subtypes selected from clustering results. The betas and standard errors were combined across Stage 1 and Stage 2 samples for each subtype under a fixed meta-analysis model weighting each strata by sample size [44]. Association test outputs were aligned to the same reference alleles and weighted z-scores were computed for each SNP. The square roots of each sample size were used as the proportional weights. Meta-analysis P-values (P_meta) were adjusted for genomic inflation. Associations with P-values < 1.67 × 10⁻⁸ were considered statistically significant, based on the standard P < 5 × 10⁻⁸ used in conventional GWAS adjusted for the 3 independent association tests performed.

Chromatin interactions

Neighboring chromatin interactions were investigated in intergenic loci using high-throughput chromatin conformation capture (Hi-C) data from the 3DIV database [45]. Topologically associating domains (TADs) were identified using TopDom [46] with a window size of 20.

Identifying subtypes in PCOS families

Quantitative trait data from the affected women (n = 73) in the family-sequencing cohort [30] were adjusted and normalized as described above. Subtype classifiers were modeled on the adjusted trait values and cluster assignments from the genotyped clustering cohort. Several classification methods were compared using 10-fold cross-validation, including support vector machine (SVM), random forest (RF), Gaussian mixed model (GMM), and quadratic discriminant analysis (QDA) [47]. The classifier with the lowest error rate was then applied to the affected women in the family-sequencing cohort to identify subtypes of PCOS in the family data. Some of the probands from the family-based cohort were included in our previous GWAS [19]. Therefore, there was some sample overlap between the training and test data: of the 893 genotyped women used to identify the original subtype clusters, 47 were also probands in the family-based cohort. Differences between subtypes in the proportion of women with DENND1A rare variants were tested using the chi-square test of independence.

Results

PCOS subtypes

Clustering was first performed in a cohort of 893 genotyped PCOS cases (Table 1, S3 Table). The clustering revealed 2 distinct phenotypic subtypes: 1) a group (23%, 207/893) characterized by higher LH and SHBG levels with relatively low BMI and Ins0 levels, which we designated “reproductive,” and 2) a group (37%, 329/893) characterized by higher BMI and Glu0 and Ins0 levels with relatively low SHBG and LH levels, which we designated “metabolic” (Fig 1). The key traits distinguishing the reproductive and metabolic subtypes were BMI, insulin, SHBG, glucose, LH, and FSH, in order of importance according to relative unpaired Wilcoxon rank–sum test statistics (Fig 2). The remaining cases (40%, 357/893) demonstrated no distinguishable pattern regarding their relative phenotypic trait distributions and were designated “indeterminate” (S4 Table, S5 Table). The reproductive and metabolic subtypes clustered along opposite ends of the SHBG versus Ins0/BMI axis, which was highly correlated with the first PC of the adjusted quantitative traits (Fig 3). The reproductive subtype was the most stable cluster, with a mean bootstrapped Jaccard coefficient (γ¯C) of 0.61, followed by the metabolic subtype with γ¯C = 0.55. The indeterminate group did not appear to capture any discernable pattern within the data (γ¯C = 0.41) and was both overlapping and intermediate between the reproductive and metabolic subtypes on the SHBG versus Ins0/BMI axis.

The clustering procedure was then repeated in an independent, ungenotyped cohort of 263 NIH PCOS cases diagnosed according to the same criteria as the genotyped clustering cohort (Table 1). The clustering yielded similar results, with a comparable distribution of reproductive (26%, 68/263, γ¯C = 0.57), metabolic (39%, 104/263, γ¯C = 0.46), and indeterminate clusters (35%, 91/263, γ¯C = 0.40) (Fig 4).

Subtype genetic associations

Genome-wide association testing was performed for each of the 3 subtypes: reproductive, metabolic, and indeterminate (Table 2). We identified alleles in 4 novel, to our knowledge, loci that were associated with the reproductive PCOS subtype at genome-wide significance (chromosome [chr]1 p36.21 PRDM2/KAZN, P = 2.23 × 10⁻¹⁰; chr2 q37.3 IQCA1, P = 2.76 × 10⁻⁹; chr4 q22.3 BMPR1B/UNC5C, P = 9.71 × 10⁻⁹; chr5 p14.2–p14.1 CDH10, P = 1.17 × 10⁻⁸) and one novel, to our knowledge, locus that was significantly associated with the metabolic subtype (chr2 q24.2–q24.3 KCNH7/FIGN, P = 1.03 × 10⁻⁸). Association testing on the indeterminate subtype replicated the 11p14.1 FSHB locus from our original GWAS [19] (Table 3; Figs 5 and 6).

The strongest association signal with the reproductive subtype appeared in an intergenic region of 1p36.21 579 kb downstream of the PRMD2 gene and 194 kb upstream from the KAZN gene (Fig 7A). The lead SNP in the locus (rs78025940; odds ratio [OR] = 4.75, 2.82–7.98 95% confidence interval [CI], P₁ = 2.16 × 10⁻¹⁰, P_meta = 2.23 × 10⁻¹⁰) was imputed (r² = 0.91) in Stage 1 only. The SNP was not genotyped in Stage 2. The lead genotyped SNP in the locus (rs16850259) was also associated with the reproductive subtype with genome-wide significance (P_meta = 2.14 × 10⁻⁹) and was genotyped only in Stage 1 (OR = 5.57, 3.24–9.56 95% CI, P₁ = 2.08 × 10⁻⁹). In ovarian tissue, the SNPs appear to be centrally located within a large 2 Mb TAD stretching from the FHAD1 gene to upstream of the PDPN gene (Fig 8).

The 2q37.3 locus spanned a 50-kb region of strong LD overlapping the 5′ end and promoter region of the IQCA1 gene (Fig 7B). The SNP rs76182733 had the strongest association in this locus (P_meta = 2.76 × 10⁻⁹) with the reproductive subtype. The signal was genotyped only in Stage 1 (OR = 5.68, 3.00–10.78 95% CI, P₁ = 2.69 × 10⁻⁹) and was imputed with an imputation r² value of 0.84.

The 4q22.3 locus spanned a 500-kb region of LD, including the 3′ ends of both the BMPR1B and UNC5C genes (Fig 7C). The most strongly associated SNP (rs17023134; P_meta = 9.71 × 10⁻⁹) in the locus was within an intronic region of UNC5C and was associated with the reproductive subtype in the Stage 1 discovery sample with genome-wide significance (OR = 3.02, 2.06–4.42 95% CI, P₁ = 1.40 × 10⁻⁸) but was not significantly associated in the Stage 2 replication analysis (OR = 1.71, 0.98–2.99 95% CI, P₂ = 7.8 × 10⁻²). The SNP was imputed with an r² of 0.87 and 0.83 in the Stage 1 and Stage 2 analyses, respectively. The most strongly associated genotyped SNP in the locus (rs10516957) was also genome-wide significant (P_meta = 1.46 × 10⁻⁸) and was located in an intronic region of BMPR1B. The genotyped SNP was nominally associated with the reproductive subtype in both the Stage 1 (OR = 2.42, 1.66–3.52 95% CI, P₁ = 6.72 × 10⁻⁶) and Stage 2 (OR = 2.40, 1.51–3.82 95% CI, P₂ = 4.7 × 10⁻⁴) analyses with nearly identical effect sizes.

In the 5p14.2–p14.1 locus, 83 kb upstream of the CDH10 gene (Fig 7D), 2 adjacent SNPs (rs7735176, rs16893866) in perfect LD were equally associated with the reproductive subtype with genome-wide significance (P_meta = 1.17 × 10⁻⁸). The SNPs were imputed in Stage 1 (OR = 5.09, 2.62–9.86 95% CI, P₁ = 1.14 × 10⁻⁸) with an imputation r² of 0.93.

The single locus containing genome-wide significant associations with the metabolic subtype was in an intergenic region of 2q24.2–q24.3 roughly 200 kb downstream from FIGN and 500 kb upstream from KCNH7 (Fig 7E). The lead SNP, rs55762028, was imputed in Stage 1 only (OR = 1.86, 0.92–3.75 95% CI, P₁ = 9.17 × 10⁻⁹, P_meta = 1.03 × 10⁻⁸). In pancreatic tissue, the lead SNPs appear to be located terminally within a 1.3-Mb TAD encompassing the FIGN gene and reaching upstream to the GRB14 gene (Fig 9).

Association testing on the indeterminate cases replicated the genome-wide significant association in the 11p14.1 FSHB locus (Fig 7F) identified in our original GWAS (14). The lead SNP (rs10835638; P_meta = 4.94 × 10⁻¹²) and lead genotyped SNP (rs10835646; P_meta = 2.75 × 10⁻¹¹) in this locus differed from the index SNPs identified in our original GWAS (rs11031006) and in the PCOS meta-analysis (rs11031005), but both of the previously identified index SNPs were also associated with the indeterminate subgroup with genome-wide significance in this study (rs11031006: P_meta = 2.96 × 10⁻¹⁰; rs11031005: P_meta = 2.91 × 10⁻¹⁰) and are in LD with the lead SNP rs10835638 (r² = 0.59) [40]. The other significant signals from our original GWAS [19] were not reproduced in any of the subtype association tests performed in this study (Table 4).

Subtypes in PCOS families

The RF classifier yielded the lowest mean subtype misclassification rate (13.2%) compared to the SVM (13.6%), GMM (17.0%), and QDA (18.1%) models, according to 10-fold cross-validation of the genotyped clustering cohort. Affected women from the family-based cohort were classified accordingly using an RF model. Seventy-three daughters of the 83 affected women from the family-based cohort had complete quantitative trait data available for subtype classification. Seventeen (23.3%) were classified as having the reproductive subtype of PCOS, and 22 (30.1%) were classified as having the metabolic subtype. Of 14 subtyped sibling pairs, only 8 were concordantly classified (57.1%); however, there was only one instance of the reproductive subtype and metabolic subtype occurring within the same nuclear family because the remaining discordant pairs each featured one indeterminate member. The proportion of affected women with one or more of the previously identified [30] deleterious, rare variants in DENND1A varied by subtype. Women classified as having the reproductive subtype of PCOS were significantly more likely to carry one or more of the DENND1A rare variants compared to other women with PCOS (P = 0.03; Fig 10). The distribution of affected women and DENND1A rare variant carriers are shown relative to the adjusted quantitative trait PCs in Fig 11.

Discussion

It is becoming increasingly clear that common, complex traits such as T2D are a heterogeneous collection of disease subtypes [8,25,48,49]. There is emerging evidence that these subtypes have different genetic architecture [7,8,25]. Consistent with these concepts, we identified reproductive and metabolic subtypes of PCOS by unsupervised hierarchical cluster analysis of quantitative hormonal traits and BMI and found novel, to our knowledge, loci uniquely associated with these subtypes with substantially larger effect sizes than those associated with PCOS disease status in GWAS [19–23]. We also found that a significantly greater prevalence of women classified with the reproductive subtype of PCOS carried at least one of the previously reported deleterious DENND1A rare variants [30] compared with those with other PCOS subtypes. These findings suggest that these subtypes are both genetically distinct as well as more etiologically homogenous [9]. Our findings are in contrast to the recent PCOS GWAS meta-analysis [22] that found that only one of 14 loci was uniquely associated with the NIH phenotype compared to non-NIH Rotterdam phenotypes. These latter findings suggest that the NIH and Rotterdam diagnostic criteria do not identify biologically distinct subtypes of PCOS. There have been previous efforts to subtype PCOS using unsupervised clustering [26–29], but no subsequent investigation into the biologic relevance of the resulting subtypes using genetic association analyses.

The key traits driving the subtypes were BMI, insulin, SHBG, glucose, LH, and FSH levels. The reproductive subtype was characterized by higher LH and SHBG levels with lower BMI and blood glucose and insulin levels. The metabolic subtype was characterized by higher BMI and glucose and insulin levels with relatively low SHBG and LH levels. The remaining 40% of cases had no distinguishable cluster-wide characteristics, and the mean trait values were between those of the reproductive and metabolic subtypes. The relative trait distributions and results of the PCAs (Figs 2, 3 and 4B) showed the reproductive and metabolic subtypes as collections of subjects on opposite ends of a phenotypic spectrum with the remaining indeterminate subjects scattered between the two. Bootstrapping and clustering in an independent cohort revealed that the reproductive and metabolic subtypes were stable and reproducible. When the GWAS was repeated, different susceptibility loci were associated with the reproductive and metabolic subtypes, suggesting that they had distinct genetic architecture. The indeterminate PCOS cases were associated with the reported locus at FSHB, but the association signal was stronger than that of our original GWAS [19], suggesting that the indeterminate group was also more genetically homogenous after the reproductive and metabolic subtypes were removed from the analysis.

Two of the loci associated with the reproductive subtype implicate novel biologic pathways in PCOS pathogenesis. The association signal on chr1 appeared downstream of and within the same TAD as the PRDM2 gene (Figs 7A and 8), which is an estrogen receptor coactivator [50] that is highly expressed in the ovary [51] and pituitary gland [52]. In an independent rare variant association study in PCOS families, PRDM2 demonstrated the fifth strongest gene-level association with altered hormonal levels in PCOS families (P = 6.92 × 10⁻³) out of 339 genes tested [30]. PRDM2 binds with ligand-bound estrogen receptor alpha (ERα) to open chromatin at ERα target genes [50,53]. PRDM2 also binds with the retinoblastoma protein [54], which has been found to play an important role in follicular development in granulosa cells [55,56].

The reproductive subtype association in the 4q22.3 locus overlapped with the BMPR1B gene, which transcribes a type-I anti-Müllerian hormone (AMH) receptor highly expressed in granulosa cells and in gonadotropin-releasing hormone (GnRH) neurons [57] that regulates follicular development [58]. Bone morphogenetic protein receptor type IB (BMPR1B), also known as ALK6 (activin receptor-like kinase 6), heterodimerizes with the transforming growth factor beta (TGF-β) type-II receptors, including AMH receptor type 2 (AMHR2), and binds AMH and other BMP ligands to initialize TGF-β signaling via the Smad proteins 1, 5, and 8 [59]. BMPR1B has been found to mediate the AMH response in ovine granulosa cells [60], and BMPR1B-deficient mice are infertile and suffer from a variety of functional defects in the ovary [61,62]. One of the BMPR1B ligand genes, BMP6, had the third-strongest gene-level association with altered hormonal levels (P = 4.00 × 10⁻³) out of 339 genes tested in our rare variant association study in PCOS families [30]. Collectively, these results make BMPR1B a compelling candidate gene in PCOS pathogenesis. These findings also support our sequencing studies that have implicated pathogenic variants in the AMH signaling pathway in PCOS [63,64].

The nature of the potential involvement in PCOS is less clear for the other loci associated with the reproductive subtype. The 2q37.3 locus overlapped with the promoter region of the IQCA1 gene. Its function in humans is not well characterized, but IQCA1 is highly expressed in the pituitary gland [52]. The 5p14.2–p14.1 locus overlapped the promoter region of the CDH10 gene. CDH10 is almost exclusively expressed in the brain [51] and is putatively involved in synaptic adhesions, axon outgrowth, and guidance [65].

The lone significant association signal with the metabolic subtype was located in an intergenic region 200–280 kb downstream of the FIGN gene, 490–570 kb upstream of KCNH7. KCNH7 encodes a voltage-gated potassium channel (subfamily H member 7, alias ERG3 [early growth response protein 3]). KCNH7 is primarily expressed in the nervous system [66] but has been found in murine islet cells [67,68]. FIGN encodes fidgetin, a microtubule-severing enzyme most highly expressed in the pituitary gland and ovary [51]. A genetic variant in FIGN was found to reduce the risk of congenital heart disease in Han Chinese by modulating transmembrane folate transport [69,70]. The TAD encompassing the association signal in this locus includes FIGN and extends upstream to the GRB14 gene (Fig 8). GRB14 plays an important role in insulin receptor signaling [71,72] and has been associated with T2D in GWAS [73]. Given the various metabolic associations for the genes in this chromosomal region, it is plausible that causal variants in this locus could impact a combination of these genes.

Despite evidence linking neighboring genes to PCOS pathways in each of the aforementioned loci, it remains possible, of course, that other, more distant genes in LD underlie the association signals. Causal variants are often up to 2 Mb away from the associated SNP, not necessarily in the closest gene [74]. Fine-mapping and functional studies are needed in order to confirm the causal variants in each of these loci. In addition, the sample sizes for the subtype GWAS were small, some of the associations were based only on imputed SNPs in Stage 1, and a replication association study has not yet been performed. However, the aforementioned functional evidence for several of the loci—particularly for PRDM2 and BMPR1B—support the validity of their associations. Further, the fact that one of the genes associated with the reproductive subtype, PRDM2, was associated with PCOS quantitative traits in our family-based analysis [30] does represent a replication of this signal by an independent analytical approach. Nevertheless, our genetic association results should be considered preliminary.

The effect sizes of the subtype alleles, particularly those associated with the reproductive subtype (OR 3.02–5.68) (Table 3), were substantially greater than the effects (OR 0.70–1.51) observed for alleles associated with PCOS diagnosis in previous GWAS [19–23]. In general, there is an inverse relationship between allele frequency and effect size [1] because alleles with larger phenotypic effects are subject to purifying selection and, therefore, occur less frequently in the population [75,76]. Accordingly, in contrast to the common variants (effect allele frequency [EAF] > 0.05) associated with PCOS in previous GWAS [19–23], the alleles associated with the subtypes were all of low frequency (EAF 0.01–0.05; Table 3). However, given the limited cohort sizes in this study, the subtype association testing did not have adequate power to detect associations with more modest effect sizes, such as those from our previous GWAS [19]. It is also possible that the large effect sizes were somewhat inflated by the so-called “winner’s curse” [77,78], but they nonetheless suggest that the subtypes were more genetically homogeneous than PCOS diagnosis in general.

In applying a subtype classifier to our family-based cohort, we found 12 affected sibling pairs in which at least one of the daughters was classified with the reproductive or metabolic subtype. Six of these sibling pairs were classified with the same subtype. There was only one discordant pairing of the reproductive subtype with the metabolic subtype. This further suggests that the reproductive and metabolic subtypes are genetically distinct in their origins. The greater prevalence of DENND1A rare variant carriers observed in women with the reproductive subtype in the family-based cohort implicates this gene in the pathogenesis of this subtype. DENND1A is known to regulate androgen biosynthesis in the ovary [79,80]; therefore, we would expect DENND1A-mediated PCOS to be more closely associated with the reproductive subtype of PCOS. However, we did not find an association between any DENND1A alleles and the reproductive subtype in the subtype GWAS, perhaps because of allelic heterogeneity or our limited power to detect associations with more modest effect sizes.

We only studied women with PCOS as defined by the NIH diagnostic criteria. Future studies will investigate whether similar reproductive and metabolic clusters are present in non-NIH Rotterdam PCOS cases. In particular, it is possible that there will be no metabolic subtype in non-NIH Rotterdam PCOS cases because these phenotypes have minimal metabolic risk [81,82]. Indeed, in a previous effort to identify phenotypic subtypes in Rotterdam PCOS cases [29], the cluster that most closely resembled the reproductive subtype represented the largest proportion of PCOS women at 44%, of whom only 78% met the NIH criteria for PCOS, whereas the cluster that most closely resembled the metabolic subtype constituted only 12% of the total sample, but 98% met the NIH diagnostic criteria. Furthermore, trait distributions may vary among women with PCOS from different geographic locations, such as in some of the sites excluded from our analysis because of incomplete quantitative trait data. For example, European PCOS cases have a lower prevalence of obesity compared to US PCOS cases [83]. Because of the within-cohort normalization of quantitative traits prior to clustering, our method is well-suited for identifying subsets of cases within populations, but therefore, it may not be suitable for directly comparing subtype membership between populations.

Our clustering cohorts included only US-based women of European ancestry. It will be of considerable importance to investigate whether subtypes are present in women with PCOS of other ancestries and geographic regions. Women with PCOS of diverse races and ethnicities have similar reproductive and metabolic features [84–86]. However, there are differences in the severity of the metabolic defects due to differences in the prevalence of obesity [83], as well as to racial/ethnic differences in insulin sensitivity [87,88]. Furthermore, the susceptibility loci associated with subtypes in other ancestry groups may differ because the low frequency and large effect size of the variants associated with the reproductive subtype in our European cohort suggests these variants are of relatively recent origin and therefore may be population-specific [1,89,90].

While the bootstrapping and clustering in an independent cohort demonstrated that the subtypes were reproducible, the Jaccard scores were relatively modest, with only the reproductive subtype yielding a mean Jaccard coefficient γ¯C > 0.6. At least part of this outcome was likely due to the fact that all traits were fitted to a normal distribution using an inverse normal transformation prior to clustering. This transformation was done in order to prevent outliers from dictating cluster formations but also likely resulted in greater cluster overlap. Consequently, the metabolic and reproductive clusters we identified appear to represent opposite ends of a phenotypic spectrum with imperfect delineation. This spectrum, however, aligns with the known pathophysiology of PCOS and is bolstered by our genetic association findings. Our approach, therefore, appears to be a more reliable way of identifying subgroups of PCOS cases who have been noted in the literature [91] but have previously been defined using only a single trait like BMI [92–95] or by diagnostic criteria that do not reflect the genetic heterogeneity of the disorder [22]. Perhaps future studies that use clustering to identify reproductive and metabolic subtypes in PCOS can omit nondistinguishing traits such as DHEAS and T in an effort to reduce noise and improve subtype delineation and reproducibility.

Our study provides support for the hypothesis that PCOS is in fact a heterogeneous disorder with different underlying biological mechanisms. As a consequence, grouping women with PCOS under a single diagnosis may be counterproductive because distinct disease subtypes will likely benefit from different interventions.

In conclusion, using an unsupervised clustering approach featuring quantitative hormonal and anthropometric data, we identified reproductive and metabolic subtypes of PCOS that appeared to have distinct genetic architecture. The genomic loci that were significantly associated with either of these subtypes include a number of new, to our knowledge, highly plausible PCOS candidate genes. Moreover, our results demonstrate that precise phenotypic delineation, resulting in more homogeneous subsets of affected individuals, can be more powerful and informative than increases in sample size for genetic association studies. Our findings indicate that further study into the genetic heterogeneity of PCOS is warranted and could lead to a transformation in the way PCOS is classified, studied, and treated.

Supporting information

S1 Checklist [docx]
STREGA checklist.

S1 Table [docx]
GWAS cohorts used in cluster analysis.

S2 Table [xlsx]
Age and BMI distributions for subjects excluded from cluster analysis.

S3 Table [dpc]
Assays used to measure quantitative traits.

S4 Table [xlsx]
Quantitative traits in genotyped PCOS cohort by cluster.

S5 Table [txt]
Subtypes of Stage 1 GWAS samples.