Genome-wide association meta-analysis identifies new endometriosis risk loci

We conducted a genome-wide association (GWA) meta-analysis of 4,604 endometriosis cases and 9,393 controls of Japanese1 and European2 ancestry. We show that rs12700667 on chromosome 7p15.2, previously found in Europeans, replicates in Japanese (P = 3.6 × 10−3), and confirm association of rs7521902 on 1p36.12 near WNT4. In addition, we establish association of rs13394619 in GREB1 on 2p25.1 and identify a novel locus on 12q22 near VEZT (rs10859871). Excluding European cases with minimal or unknown severity, we identified additional novel loci on 2p14 (rs4141819), 6p22.3 (rs7739264) and 9p21.3 (rs1537377). All seven SNP effects were replicated in an independent cohort and produced P < 5 × 10−8 in a combined analysis. Finally, we found a significant overlap in polygenic risk for endometriosis between the European and Japanese GWA cohorts (P = 8.8 × 10−11), indicating that many weakly associated SNPs represent true endometriosis risk loci and risk prediction and future targeted disease therapy may be transferred across these populations.

Encouraged by the WNT4 association and with accumulating evidence for many complex traits that the number of discovered variants is strongly correlated with experimental sample size 6 , we sought to increase the ratio of controls to cases in the Australian GWA cohort and to perform a formal meta-analysis of the Australian (QIMR), UK (OX) and Japanese (BBJ) GWA data.
To increase the power of the Australian GWA dataset we matched the existing QIMR cases and controls 2 on ancestry to individuals from the Hunter Community Study (HCS) 7 . After stringent quality control (QC), the combined QIMRHCS GWA cohort consisted of 2,262 endometriosis cases and 2,924 controls, increasing the number of controls by 1,054 and the Australian effective sample size by 24%. We also performed more stringent QC incorporating the OX dataset, resulting in a revised OX GWA cohort of 919 endometriosis cases and 5,151 controls. All cases in the QIMRHCS and OX studies have surgically confirmed endometriosis and disease stage from surgical records using the rAFS classification system 8 , subjects are grouped into stage A (stage I or II disease or some ovarian disease with a few adhesions; n = 1,680, 52.8%) or stage B (stage III or IV disease; n = 1,357, 42.7%), or unknown (n = 144, 4.5%). Details of the final GWA and independent replication case-control cohorts are summarized in Table 1 and a schematic of our study design is provided in Fig. 1.
Given the substantially greater genetic loading of moderate to severe (Stage B) endometriosis (rAFS stage III or IV disease) compared to minimal (Stage A) endometriosis (rAFS stage I or II disease) 2 , a secondary analysis was performed for the SNPs reaching genome-wide suggestive association, where the association results from QIMRHCS and OX Stage B cases versus controls, were meta-analyzed with the BBJ association results (stage information not available).
Interestingly, the most associated genotyped SNP at 9p21.3 (rs1537377) is 55 kb centromeric to the genome-wide significant SNP reported in the original BBJ GWA 1 (rs10965235) located in the CDKN2BAS gene, and 49 kb 3 to the transcription end site of CDKN2BAS. SNP rs10965235 is monomorphic in Caucasian populations and we investigated the independence of rs10965235 and rs1537377 in the BBJ GWA data. Firstly, in the BBJ GWA data, alleles of rs10965235 and rs1537377 are very weakly correlated, with linkage disequilibrium (LD) metrics of r 2 = 0.028 and D = 0.461. Secondly, the allelic association P values for rs10965235 and rs1537377 are P = 1.6 × 10 −4 and P = 1.8 × 10 −2 , respectively. After conditioning on rs10965235, weak residual association remains at rs1537377 (P = 9.0 × 10 −2 ). Consequently, the data suggest there may be two independent genetic risk factors near the CDKN2BAS locus on 9p21.3. CDKN2BAS is a long noncoding RNA adjacent to and transcribed from the opposite strand to CDKN2B (p15), CDKN2A (p16) and ARF (p14). Loss of heterozygosity of CDKN2A and hypermethylation of the CDKN2A promoter have been reported in endometriosis 10,11 .
To further validate the seven SNPs implicated by the meta-analysis, we carried out a replication study using a cohort of 1,044 cases and 4,017 controls obtained from the BioBank Japan independent of the BBJ GWA cohort. As shown in the forest plots of risk allele effects estimated using all cases versus controls (Fig. 3), the effects (ORs) were in the same direction for all seven implicated SNPs across the GWA and replication cohorts. With the exception of rs12700667, which was previously replicated (P = 1.2 × 10 −3 ) in 2,392 cases and 2,271 controls from the US 2 , and rs4141819 (with a marginal P = 5.1× 10 −2 ), all SNPs were replicated at the nominal P < 0.05 threshold ( Table 2). All seven SNPs surpass the conventional genome-wide significant threshold of P ≤ 5 × 10 −8 after combined analysis of the GWA and replication cases and controls ( Table 2). A conservative adjustment of the rs4141819 total P values (P all = 8.5 × 10 −8 ; P stageB = 4.1 × 10 −8 ) for performing two independent GWA studies (all and stage B endometriosis cases versus controls) would produce P > 5 × 10 −8 (P all-adjusted = 1.7 × 10 −7 ; P stageB-adjusted = 8.2 × 10 −8 ). However, the accurately imputed (Rsq > 0.95) SNP rs2861694 (P stageB = 7.9 × 10 −9 ), in strong LD with rs4141819 (r 2 = 0.981, D = 1.0; and r 2 = 0.867, D = 1.0, in the 379 European and 286 Asian 1000 Genomes reference samples, respectively), would remain genome-wide significant (P stageB-adjusted = 1.6 × 10 −8 ).  Fig. 17d), reveals a significant preponderance of small P values <10 −3 , suggesting many of these nominally significant SNPs likely represent true signals 12 . To further examine the shared genetic risk across our European and Japanese populations we performed polygenic prediction analysis 13 to evaluate whether the aggregate effects of many variants of small effect in the BBJ GWA cohort, could predict affection status in the European GWA cohorts. The BBJ-derived risk scores significantly predicted affection status in the QIMRHCS (R 2 = 0.0064; P = 6.9 × 10 −7 ), OX (R 2 = 0.0057; P = 9.6 × 10 −6 ) and combined QIMRHCS+OX all endometriosis case-control sets (R 2 = 0.0054; P = 8.8 × 10 −11 ). For the individual and combined QIMRHCS and OX case-control sets, the variance explained peaked in the SNP sets with BBJ GWA P < 0.1, using all GWA meta-analysis SNPs (Fig. 4a) and after excluding all SNPs within ±2500 kb of the seven implicated SNPs listed in Table 1 (Fig.  4b). Analogously, performing the prediction in reverse, the QIMRHCS+OX-derived risk scores significantly predicted affection status in the BBJ case-control set (R 2 = 0.0106; P = 3.3 × 10 −6 ) (Supplementary Fig. 18 and Supplementary Note).
A gene-based GWA analysis using VEGAS 14 , which accounts for gene size and LD between SNPs, revealed 1,184 genes with a combined P ≤ 0.05 and the top three ranked genes associated with endometriosis to be WNT4 on 1p36.12 (P = 5.0 × 10 −9 ), VEZT on 12q22 (P = 5.7 × 10 −7 ) and GREB1 on 2p25.1 (P = 2.5 × 10 −5 ) (Supplementary Table 3). In addition to having genome-wide significant SNPs near these three genes, the WNT4 and VEZT genes easily surpass our conservative gene-based significant association threshold of P ≤ 2.85 × 10 −6 (calculated as P = 0.05 / 17,538 independent genes). WNT4 encodes for wingless-type MMTV integration site family, member 4 and is important for the development of the female reproductive tract 15 and steroidogenesis 16 . VEZT encodes vezatin, an adherens junction transmembrane protein that is down regulated in gastric cancer 17 . GREB1 encodes growth regulation by estrogen in breast cancer 1, an early response gene in the estrogen regulation pathway involved in hormone dependent breast cancer cell growth 18  In conclusion, given their high gene-based ranking, proximity to genome-wide significant SNPs, known pathophysiology and reported gene expression (Supplementary Note and Supplementary Fig. 19), the WNT4, VEZT and GREB1 genes are strong targets for further studies aimed at understanding the molecular pathogenesis of endometriosis. Our results also suggest that a considerable number of SNPs nominally implicated (e.g. P < 0.1) in the European and Japanese GWA cohorts represent true endometriosis risk loci. Moreover, the significant overlap in common polygenic risk for endometriosis indicates genetic risk prediction and future targeted disease therapy may be transferred across these populations.

GWA samples and phenotyping
Initially, 2,351 surgically-confirmed endometriosis cases were drawn from women recruited by The Queensland Institute of Medical Research (QIMR) study 19  Ethics Committees, and the Oxford regional multi-centre and local research ethics committees. Informed consent was obtained from all participants prior to testing 2 .
All Japanese GWA case and control samples were obtained from the BioBank Japan (BBJ) at the Institute of Medical Science, the University of Tokyo. Quality control procedures for the OX genotype data resulted in the removal of SNPs with a genotype call rate < 0.99 and/or heterozygosity < 0.31 or > 0.33. Genome-wide IBS was estimated for each pair of individuals and one individual from each duplicate or related pair (IBS > 0.82) was removed. Genotype data were combined with CEU, CHB&JPT and YRI genotype data from HapMap 3 and individuals of non Northern European ancestry were identified using EIGENSOFT and subsequently removed. SNPs with a genotype call rate < 0.95 were removed, and this threshold was increased to 0.99 for SNPs with MAF < 0.05. In addition, SNPs showing a significant a) deviation from HWE (P < 1 × 10 −6 ); b) difference in call rate between 58BC and NBS control groups (P < 1 × 10 −4 ); c) difference in allele/ genotype frequency between control groups (P < 1 × 10 −4 ); d) difference in call rate between cases and controls (P < 1 × 10 −4 ) and e) a MAF < 0.01 were removed. 2 The BBJ cases and controls were genotyped using the Illumina HumanHap550v3 Genotyping BeadChip. Quality control included sample call rate ≥ 0.98, identity-by-state to exclude close relatedness samples and principal component analysis to exclude non-Asian samples. We also performed SNP quality control (call rate of ≥ 0.99 in both cases and controls and Hardy-Weinberg equilibrium test P ≥ 1.0 × 10 −6 in controls); 460,945 SNPs on all chromosomes passed the quality control filters and were further analyzed.

GWA meta-analysis
For SNPs passing QC, tests of allelic association (--assoc) were performed using PLINK 23 in the separate QIMRHCS, OX and BBJ GWA datasets. The primary meta-analysis of all endometriosis cases versus controls in the QIMRHCS, OX and BBJ GWA data was performed using a fixed-effect (inverse variance-weighted) model, where the effect size estimates, or -coefficients, are weighted by their estimated standard errors, utilizing the GWAMA software 24 .
The threshold of 7.2 × 10 −8 for GWA studies of dense SNPs and resequence data 25 proposed by Dudbridge and Gusnanto 26 was utilized to indicate genome-wide significant association, while SNPs with P ≤ 10 −5 were considered to show a suggestive association [as used in the online 'Catalog of Published Genome-Wide Association Studies'].
Also, given the substantially greater genetic loading of moderate to severe (stage B) endometriosis (rAFS stage III or IV disease) compared to minimal (stage A) endometriosis (rAFS stage I or II disease) 2 , a secondary analysis was performed for suggestive SNPs (P ≤ 10 −5 ); where the association results from QIMRHCS and OX stage B cases versus controls, were meta-analyzed with the BBJ association results. As previously demonstrated 2 , the exclusion of minimal endometriosis cases has the potential to enrich true genetic risk effects, even taking into account the reduced sample size.

Consistency of allelic effects across studies was examined utilizing the Cochran's Q test 27 .
Between-study (effect) heterogeneity was indicated by Q statistic P values < 0.1 28 . Metaanalysis of SNPs associated with fixed-effect P ≤ 10 −5 and showing evidence of effect heterogeneity were also analyzed using the recently developed Han and Eskin's random effects model (RE2) implemented in the Metasoft software 29 . In contrast to the conventional DerSimonian-Laird random effects (RE) model 30 , the RE2 model increases power under heterogeneity 29 .

Genotype imputation analysis
In order to assess the impact of variants not present on the Illumina BeadChips and better define the associated regions, we imputed genotypes within ±2500 kb of the most significant genotyped SNP using the full reference panel from the 1000 Genomes project Interim Phase I Haplotypes (2010-11 data freeze, 2011-06 haplotypes). Imputation was performed separately for the QIMRHCS, OX and BBJ GWA datasets with only the overlapping genotyped SNPs within ±2500kb of the most significant genotyped SNP, using the MaCH and minimac programs 31,32 and following the two-step approach outlined in the online 'Minimac: 1000 Genomes Imputation Cookbook'. Analysis of imputed genotype dosage scores was performed using mach2dat 31,32 and PLINK. The quality of imputation was assessed by means of the Rsq statistic. Results for poorly imputed SNPs, defined as having an Rsq < 0.3, were subsequently removed. The results from association analysis of imputed data in the QIMRHCS, OX and BBJ datasets were then combined via meta-analysis of thecoefficients weighted by their estimated standard errors using GWAMA.

Replication samples and genotyping
Independent of the BBJ GWA case-control cohort, a total of 1,044 cases and 4,017 controls were obtained from the BioBank Japan and utilized for replication. We note that 653 of these 1,044 cases were also utilized in a small GWA study (Adachi et al. 2010) of 696 cases and 825 controls 9 . To utilize all available association data for rs13394619 maximally, given there is incomplete overlap between the Adachi and our replication cases and zero overlap between the controls, we worked with the published results for rs13394619 in Adachi et al

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts (2010) and the results from comparing our non-overlapping 391 replication cases to our 4,017 replication controls.
The seven SNPs (rs7521902, rs13394619, rs4141819, rs7739264, rs12700667, rs1537377 and rs10859871) reaching genome-wide significance in the GWA meta-analysis were genotyped in the independent Japanese replication cohort using the multiplex PCR-based Invader assay (Third Wave Technologies), as previously described 1 .

Replication and total association analyses
Tests of allelic association were performed using PLINK in the independent Japanese replication cohort. Because only the associations in the same direction would be considered as replicated, one-sided P values were obtained by halving the standard (two-sided) PLINK P values. To determine the total evidence for association, the one-sided replication P values were meta-analyzed with the QIMRHCS, OX, BBJ [and Adachi 9 500K (290 cases and 262 controls) and 6.0 (406 cases and 563 controls) for rs13394619] GWA P values using METAL 33 . The P values observed in each case-control cohort were converted into a signed Z-score. Z-scores for each allele were combined across samples in a weighted sum, with weights proportional to the square-root of the sample size for each cohort 34 . Given that our cohorts have unequal numbers of cases and controls, we utilized the effective sample size, where N eff = 4 / (1 / N cases + 1 / N controls ) 33 . We also performed meta-analysis of thecoefficients weighted by their estimated standard errors using GWAMA to estimate the overall odds ratio and 95% CI for the genome-wide significant SNPs.

Polygenic prediction
The aim of the prediction analysis was to evaluate the aggregate effects of many variants of small effect. We summarized variation across nominally associated loci into quantitative scores and related the scores to disease state in independent samples. Although variants of small effect (e.g., genotype relative risk of 1.05) are unlikely to achieve even nominal significance, increasing proportions of "true" effects will be detected at increasingly liberal P value thresholds, e.g. P < 0.1 (i.e., ~10% of all SNPs), P < 0.2, etc. Using such thresholds, we defined large sets of "allele specific scores" in the "discovery" sample of the Japanese BioBank (BBJ) endometriosis case-control set (1,423 cases, 1,318 controls) to generate risk scores for individuals in the "target" sample of the QIMRHCS (2,262 cases, 2,924 controls), OX (919 cases, 5,151 controls) and combined European (QIMRHCS+OX) endometriosis case-control sets (3,181 cases, 8,075 controls). The term risk score is used instead of risk, as it is impossible to differentiate the minority of true risk alleles from the non-associated variants. In the discovery sample, we selected sets of allele specific scores for SNPs with the following levels of significance; P < 0.01, P < 0.05, P < 0.1, P < 0.2, P < 0.3, P < 0.4, P < 0.5, P < 0.6, P < 0.7, P < 0.8, P < 0.9, P < 1.0. For each individual in the target sample, we calculated the number of score alleles that they possessed, each weighted by the log odds ratio from the discovery sample. To assess whether the aggregate scores reflect endometriosis risk, we tested for a higher mean score in cases compared to controls. Logistic regression was used to assess the relationship between target sample disease status and aggregate risk score. Nagelkerke's pseudo R 2 was used to assess the variance explained.
Prediction was performed using all 407,632 SNPs overlapping the QIMRHCS, OX and BBJ GWA datasets, and after excluding the 6,163 SNPs within ±2500 kb of the seven implicated SNPs listed in Table 1. We also performed the predictions in reverse, using QIMRHCS +OX-derived risk scores to predict affection status in the BBJ case-control set.

Gene-based association analysis
Gene-based approaches can be more powerful than traditional individual-SNP-based approaches in the presence of allelic heterogeneity. Therefore, utilizing the QIMRHCS, OX Europe PMC Funders Author Manuscripts and BBJ GWA data, we performed a genome-wide gene-based association study using VEGAS 14 . Briefly, for the 407,632 overlapping SNPs, the P values from the European GWA study (i.e., FE meta-analysis of QIMRHCS and OX GWA data) and the P values from the Japanese (BBJ) GWA study were analyzed separately using VEGAS. The VEGAS test incorporates evidence for association from all SNPs across a gene and accounts for gene size (number of SNPs) and LD between SNPs by using simulations from the multivariate normal distribution. The resulting European and Japanese gene-based P values were meta-analyzed using Stouffer's Z-score combined p-value method 34 . A total of 17,538 genes (including 50 kb 5 and 3 of their transcription start and end site, respectively 14 ) contained association results for ≥1 SNP, so a Bonferroni adjusted significance threshold of P ≤ 2.85 × 10 −6 (0.05 / 17,538) was utilized to indicate genome-wide gene-based significant association.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Acknowledgments
We acknowledge with appreciation all the women who participated in the QIMR, OX and BBJ studies. We thank Endometriosis Associations for supporting the study recruitment. We also thank the many hospital directors and staff, gynecologists, general practitioners and pathology services in Australia, the UK and the United States who provided assistance with confirmation of diagnoses. We thank Sullivan and Nicolaides Pathology and the Queensland Medical Laboratory Pathology for pro bono collection and delivery of blood samples and other pathology services for assistance with blood collection. The Hunter Community Study team thanks the men and women of the Hunter region who participated in the study.   Table 2 Summary of the GWA and replication study results for the seven genome-wide significant loci