Genetic predisposition to uterine leiomyoma is determined by loci for genitourinary development and genome stability.

Uterine leiomyomas (ULs), also known as fibroids or myomas, are benign smooth muscle tumors of the uterine wall. They are extremely common; approximately 70% of women develop ULs before menopause ( Stewart et al., 2017 ). The symptoms, occurring in one fifth of women, include excessive menstrual bleeding, abdominal pain and pregnancy complications ( Stewart et al., 2017 ). In most cases, durable treatment options are invasive ( Stewart, 2015 ). ULs cause a substantial human and economic burden, and the annual cost of treating these tumors has been approximated to be as high as $34 billion in the United States, higher than the combined cost of treating breast and colon cancer ( Cardozo et al., 2012 ). Earlier studies have indicated strong genetic influence in UL susceptibility based on linkage ( Gross, 2000 ), population disparity ( Wise et al., 2012 ) and twin studies ( Luoto et al., 2000 ). The most striking UL predisposing condition thus far characterized is hereditary leiomyomatosis and renal cell cancer (HLRCC) syndrome, caused by high-penetrance germline mutations in the Fumarate hydratase ( FH ) gene ( Multiple Leiomyoma Consortium et al., 2002 ; Launonen et al., 2001 ). Genome-wide association studies (GWAS) have proposed several low-penetrance risk loci but few unambiguous predisposing genes have emerged. Cha et al. reported loci in chromosome regions 10q24.33, 11p15.5 and 22q13.1 based on a Japanese patient cohort ( Cha et al., 2011 ). The 11p15.5 locus - near the Bet1 golgi vesicular membrane trafficking protein like ( BET1L ) gene - was later replicated in Caucasian ancestry ( Edwards et al., 2013a ). The 22q13.1 locus has been replicated in Caucasian, American and Saudi Arabian populations suggesting trinucleotide repeat containing 6B ( TNRC6B ) as a possible target gene ( Edwards et al., 2013a ; Aissani et al., 2015 ; Bondagji et al., 2017 ). Further UL predisposition loci have been suggested at 1q42.2 and 2q32.2 by Zhang et al ( Zhang et al., 2015 ). and, at 3p21.31, 10p11.21 and 17q25.3 by Eggert et al ( Eggert et al., 2012 ). A recent work reported cytohesin 4 ( CYTH 4) at 22q13.1 as a novel candidate locus in African ancestry ( Hellwege et al., 2017 ). While multiple loci and genes have been implicated through these valuable studies it is not straightforward to connect any of them mechanistically to UL development. Most ULs show somatic site-specific mutations at exons 1 and 2 of the mediator complex subunit 12 ( MED12 ) gene ( Mäkinen et al., 2011 ; Heinonen et al., 2014 ). These observations together with further scrutiny of driver mutations, chromosomal aberrations, gene expression, and clinicopathological characteristics have led to identification of at least three mutually exclusive UL subtypes; MED12 mutant, Fumarate Hydratase ( FH ) deficient, as well as High Mobility Group AT-Hook 2 ( HMGA2 ) overexpressing lesions ( Mehine et al., 2016 ). Here we report the most powerful GWAS on uterine leiomyoma to date, and novel genome-wide significant UL susceptibility loci with plausible adjacent predisposition genes. These genes associate UL genesis to two distinct biological mechanisms: Genome stability related processes are implicated by genes Tumor Protein P53 ( TP53 ) and ATM Serine/Threonine Kinase ( ATM ) together with the telomere maintenance genes Telomerase Reverse Transcriptase ( TERT ) , Telomerase RNA Component ( TERC ) and STN1-CST Complex Subunit ( OBFC1 ). The other prominent group is genes relevant for genitourinary development, specifically Wnt Family Member 4 ( WNT4 ), Wilms Tumor 1 ( WT1 ), Spalt Like Transcription Factor 1 ( SALL1 ) , Estrogen Receptor 1 ( ESR1 or ERα ) , Growth Regulation By Estrogen In Breast Cancer 1 ( GREB1 ) , Forkhead Box O1 ( FOXO1 ) , Doublesex and Mab-3 Related Transcription Factor 1 ( DMRT1 ) and CD44 Molecule ( CD44 ). Our analysis of the X chromosome identifies a risk allele near MED12 that drives UL tumorigenesis towards somatic MED12 mutations. We report altogether 22 genome-wide significant susceptibility loci and compile them into a polygenic risk score. The UL association is then replicated in six independent cohorts of different ethnic origins: individuals of African origin are characterized by the highest risk load. Finally, we investigate the risk alleles’ association to clinical features, molecular UL subtypes, telomere length, gene expression and DNA methylation.

Figure 1 provides an outline of this study. At discovery stage 1,428 SNPs emerging from 22 distinct genetic loci passed the genome-wide significance level of 5 × 10 −8 . Figure 2 displays a Manhattan plot of these associations (15,453 UL cases and 392,628 controls; linear mixed model). Two of the significant loci (359/1,428 SNPs) were found on the X chromosome. After linkage disequilibrium (LD; r 2  ≤0.3) pruning the significant SNPs, a total of 50 LD-independent associations remained: the resulting SNPs are given in Appendix 1—table 2 , and the lead SNPs are summarized in Table 1 . GRS, genomic risk score. NFBC, Northern Finland Birth Cohort. ( A ), Manhattan plot of the UK Biobank GWAS on 15,453 UL cases and 392,628 controls. On Y-axis, logarithm transformed association values, and on X-axis, autosomes and the X chromosome. The blue horizontal line denotes genome-wide significance (p=5 × 10 −8 ). Gene symbols shown for reference. ( B ), MED12 region in more detail. ENCODE tracks (details in Supplementary Methods) are shown for reference. ( C ), The risk allele near MED12 (rs5937008) is observed with a significant increase in number of MED12 -mutation-positive tumors (p=0.009; negative binomial regression; RR = 1.23 per risk allele; n = 457 Helsinki cohort patients). The numbers for B allele frequency (B Freq), odds-ratio (OR, where B is the effect allele) and association (P) are based on the UKBB cohort (15,453 UL cases). Gene symbols are shown for reference. The genomic coordinates follow hg19 and dbSNP build 147. All genome-wide significant (p<5 × 10 −8 ) loci and their highest-association SNP are shown. * Previously implicated predisposition to ULs. # Previous associations to endometriosis, lung adenocarcinoma, glioma or telomere length; see literature in Appendix 1—table 11 . Appendix 1—figure 1 displays the regional structure of each locus together with flanking association values, linkage disequilibrium (LD) and genome annotation. Annotation tracks are included for tissue-specific data on open chromatin, topologically associating domains (TAD) and other regulatory features (details in Supplementary Methods). A polygenic risk score ( Abraham and Inouye, 2015 ) was compiled based on the discovery stage associations. After LD pruning (r 2  ≤0.3) the discovery-stage SNPs, 50 SNPs from the 22 distinct loci passed for the initial genomic risk score (GRS; Appendix 1—table 4 ). The SNP weights were based on UKBB log-odds. We applied this initial GRS of 50 SNPs to the Helsinki cohort and identified a significant association to the UL phenotype (p=8.3 × 10 −10 ; adjusted p=1.1 × 10 −8 ; one-tailed Wilcoxon rank-sum; W = 1.69 × 10 6 ; 457 cases and 8899 female controls). The second stage GWAS combined the UKBB and Helsinki cohorts for a meta-analysis approach. The genome-wide statistics revealed rs117245733, at 13q14.11, as the only SNP with a suggestive (p<10 −5 ) association in both the UKBB (OR = 1.26; p=4.2 × 10 −9 ) and Helsinki (OR = 1.82; p=8.1 × 10 −6 ) cohorts. Figure 3 shows the regional structure and combined association (fixed effect model p=3.1 × 10 −12 ) at the locus: the SNP resides on a gene poor region, at a conserved element that displays activity in uterus-specific H3K27ac and DNaseI data (see ENCODE track details in Supplementary Methods). The SNP is independent of the group of associations at FOXO1 (r 2  = 0.0; Figure 3 ). ( A ) meta-analysis P-values and the genomic context at the locus. Gene symbols and ENCODE tracks (details in Supplementary Methods) are shown for reference; coordinates follow hg19. ( B ) Hi-C, TADs and CpG methylation around the locus with a 1 Mb flank. The needle plot shows the meQTL associations (dashed lines at 10% FDR; green line denotes the SNP; gray ticks denote all CpGs tested; blue needle for positive coefficient, red for negative coefficient) for tumors (above x-axis; n AA  = 53, n AB  = 3) and normals (below x-axis; n AA  = 33, n AB  = 2). ( C ) UCSC genome browser tracks related to conservation and regulation at the locus. The meta-analysis identified altogether 112 genome-wide significant SNPs not seen in the discovery stage: seven of those were LD-independent (r 2  ≤0.3; Appendix 1-table 3 ) and their UKBB log-odds weights were appended to the initial GRS model. The final GRS model of 57 SNPs and their UKBB-based weights is given in Appendix 1—table 4 . Supplementary file 1 gives further details on the meta-analysis results and heterogeneity statistics. The third stage replicated the observations in NFBC and in five different ethnic groups. In NFBC, the SNP identified in the stage two meta-analysis, rs117245733 at 13q14.11, was replicated (p=0.034; linear mixed model; OR = 1.50; 95% CI 1.03 – 2.19). Additional analysis of all 57 SNPs did not reveal other associations: Supplementary file 2 gives further details on the meta-analysis results and heterogeneity statistics. The association between the GRS and UL phenotype was significant (p=1.1 ×10 −5 ; Wilcoxon rank-sum; adjusted p=1.1 × 10 −4 ; one-tailed; W = 4.7 × 10 5 ) in NFBC. These case-control distributions of GRS are displayed in Figure 4 . On top, GRS association to MED12 mutation status. The rest show GRS association to the UL phenotype in six independent replication cohorts. Associations ( P ) and test statistics ( W ) are from Wilcoxon rank-sum tests. Only females were included as the control samples. The X-axes show the GRS distributions for each phenotype. UL susceptibility is known to vary by ancestry ( Wise et al., 2012 ). Five different ethnic groups - African, Caribbean, Irish, Indian and ‘other white’ background - were available from the UKBB cohort. A total of 2,212 UL cases and 21,054 female controls could be utilized for replication ( Appendix 1—table 1 ). Supplementary file 3 includes all the 57 SNPs and their summary statistics in these five cohorts, together with heterogeneity estimates. Due to the small cohort sizes, none of the single-SNP associations passed genome-wide significance. The GRS model replicated with a significant phenotype association in all five ethnicities (Appendix 1-Table 6). A summary of test statistics, GRS distributions and the numbers of cases and controls for each population is given in Figure 4 . A more detailed summary of the GRS model and receiver operating characteristic (ROC) curve of each cohort are given in Appendix 1—figure 5 . The self-reported ‘Black African’ (mean GRS 4.83) had an outstanding risk-load compared to Caucasian (self-reported ‘White Irish’; mean GRS 4.04) background ( Figure 4 ; Wilcoxon rank-sum p<10 −15 ). As expected ( Wise et al., 2012 ), the African ethnicity displayed an increased prevalence (19%) compared to the Irish (6%). Assuming that the observed GRS weights have a linear relationship to the true risk, the GRS difference between African and Irish ancestries explains 9.0% of the increased prevalence in the African population. Similar population-specific GRSs could be estimated for the seven populations in the gnomAD database ( Appendix 1—table 7 ). Appendix 1—figure 6 shows an overview of the GRS for each of the populations. African ancestry has been shown to carry a two-to-three times higher prevalence when compared to Caucasian ancestry ( Wise et al., 2012 ). Based on the gnomAD frequencies, the increased GRS of African ancestry explains between 8 – 16% of this population difference. The number of ULs per patient had a significant positive association to GRS (negative binomial regression p=0.001; adjusted p=0.0032; rate ratio 1.25; 95% CI 1.09 – 1.43 for one-unit increase in GRS; Appendix 1—figure 7 ). No association was found between GRS and age at hysterectomy ( Appendix 1—table 6 ). Testing the 57 GRS SNPs separately did not reveal any associations that pass FDR ( Appendix 1—table 5 ). Our UL set of 1481 lesions included 1159 (78%) mutation-positive and 322 mutation-negative tumors. The occurrence of mutant tumors did not distribute evenly among the 457 patients. In total 221 (48%) and 123 (27%) patients had all their tumors identified as either MED12 -mutation-positive or -negative, respectively, suggesting that genetic or environmental factors contribute to the preferred UL type in affected individuals, as previously observed ( Mäkinen et al., 2011 ). Indeed, the 221 mutation positive patients were found to have a significantly higher GRS (Wilcoxon rank-sum p=7.9 × 10 −4 ; adjusted p=0.0032; two-sided; W = 1.6 × 10 4 ). This difference in GRS distributions is visualized in Figure 4 . Comparison against the population controls (n = 8899 females) revealed that the above-mentioned patient groups differ by their effect size: the MED12 -mutation-positive (221) subset of patients had an odds ratio of 2.28 for one-unit increase in GRS (95% CI 1.80 – 2.88) compared to the controls, while the mutation-negative (123) subset had an odds ratio of 1.20 (95% CI 0.88 – 1.66). Thus, the majority of the compiled case-control association signal had arisen from the MED12 -mutation-positive subset of the patients. The number of MED12 -mutation-positive tumors per patient had a significant positive association to GRS (p=3.2 × 10 −4 ; adjusted p=0.002; negative binomial model rate ratio 1.43; 95% CI 1.13 – 3.83 for one-unit increase in GRS; Appendix 1—figure 8 ). No association between the number of MED12 -mutation-negative tumors and GRS was found (adjusted p=0.053; Appendix 1—figure 8 ). The GWAS signal near MED12 was inspected for any associations to somatic MED12 mutations. Strikingly, the risk allele (rs5937008) did significantly increase the number of MED12 -mutation-positive tumors (p=0.0087; negative binomial model rate ratio 1.23; 95% CI 1.05 – 1.44). Among our 457 patients, the median number of MED12 -mutation-positive tumors increased from one to two for the risk allele carriers. The risk locus and its effect on the number of MED12 -mutation-positive tumors is visualized in Figure 2 . An additional analysis of each of the 57 GRS SNPs did not reveal any further associations ( Appendix 1—table 5 ). All the genome-wide significant SNPs from UKBB and the meta-analysis stage (altogether 1,540 SNPs) were tested with a permutation based approach. In total 34 and 24 genes passed the local permutation significance threshold (p<0.05) for tumor and matched myometrium data, respectively ( Appendix 1—table 8 ). Among the hits in tumors were WNT4 (p=0.01; permutation test) and CDC42 (p=0.03) at 1 p, TNRC6B (p=0.02) at 22q, FOXO1 (p=0.03) at 13q, and DMRT1 (p=0.04) at 9 p. None of the local associations passed a genome-wide FDR of 10%. No significant association was observed between the risk allele and MED12 expression (rs5936989; Appendix 1—figure 11 ). The full list of eQTL statistics can be found from Supplementary file 4 . Our analysis of the 57 GRS SNPs revealed altogether 17,030 (9,466 in tumors and 7564 in matched myometrium) cis methylation quantitative trait loci (cis-meQTL) with nominal p<0.05. Of these, 145 passed a 10% FDR. Of the plausible predisposition genes, FOXO1 , TERT and WNT4 showed significant meQTL associations ( Appendix 1—table 9 ). All the cis-meQTLs and annotation for their genomic context are in Supplementary file 5 . The UL predisposition loci at TERT, TERC and OBFC1 were examined for an effect on telomere length. Overall the telomere length was significantly shorter in tumors than in adjacent matched myometrium (p=0.01; Kruskal-Wallis), as previously reported ( Rogalla et al., 1995 ; Bonatz et al., 1998 ). One of the risk alleles at TERT (rs2736100) was significantly associated with shorter telomere length (p=0.01; Kruskal-Wallis) ( Appendix 1—figure 12 ). Adjusting for the patient age did not explain away the association. The association was not seen in myometrium. The other two LD-independent SNPs at TERT , rs72709458 and rs2853676, or the SNPs at TERC (rs10936600) and OBFC1 (rs1265164) did not show association to telomere length (p=0.24, p=0.57, p=0.07 and p=0.48, respectively; Kruskal-Wallis). The combined effect of TERT (rs72709458, rs2736100, rs2853676), TERC (rs10936600) and OBFC1 (rs1265164) had a negative trend with telomere length (p=0.055; linear model 95% CI −408.5 – 4.7 per one risk allele; see Appendix 1—figure 13 ). In whole genome sequencing data, no association was detected between genotype and the number of somatic structural variants. The DEPICT framework ( Pers et al., 2015 ) was ran using the genome-wide significant SNPs from the UKBB cohort, in total 1,069/1,428 autosomal SNPs. The resulting target gene prioritization, pathway enrichment and tissue enrichment results are given in Supplementary file 6 . The analysis did not reveal any significant enrichments with the exception of one pathway related to induced stress. ATM was the highest ranking target gene, and uterus/myometrium were among the highest ranking tissue types. Previous UL association studies ( Cha et al., 2011 ; Zhang et al., 2015 ; Eggert et al., 2012 ; Hellwege et al., 2017 ) have reported altogether seven genome-wide significant UL susceptibility loci. Two out of the seven loci - that is, 22q13.1 (at TNRC6B ) and 11p15.5 (at BET1L ) - replicated in UKBB using 15,453 cases and 392,628 controls. Cha et al ( Cha et al., 2011 ). highlight OBFC1 (at 10q24.33) as a candidate gene and, while the SNP that they reported does not replicate in UKBB, the OBFC1 region is identified in our discovery stage (rs1265164; Table 1 ). See Appendix 1—table 10 for a summary of these results.

The UK Biobank genotype-phenotype data revealed 22 novel predisposition loci for UL, most of them in close proximity to highly plausible predisposition genes. The combined UL risk of these loci was replicated in a subsequent analysis of the polygenic risk score (GRS) in six independent cohorts from different ethnic backgrounds. Our multi-ethnic replication implies that the discovered loci are indeed involved in UL development, and the early UL association studies have likely been underpowered to detect them. Three previously reported loci, at OBFC1 ( Cha et al., 2011 ), TNRC6B ( Cha et al., 2011 ; Edwards et al., 2013a ; Aissani et al., 2015 ; Bondagji et al., 2017 ) and BET1L ( Cha et al., 2011 ; Edwards et al., 2013b ), were also validated, however, the mechanistic connection to UL development remains obscure for the latter two. Though simple association is not sufficient to formally prove causality, 14 out of the 22 risk loci harbor plausible predisposition genes. These genes can be divided into two groups: TERT, TERC, OBFC1 (all involved in telomere length), ATM and TP53 guard stability of the genome. ESR1, GREB1, WT1, MED12, WNT4, FOXO1, DMRT1, SALL1, and CD44 play a role in genitourinary development. Estrogen is a well-known inducer of UL growth ( Borahay et al., 2015 ). The top association at 6q25.2 (rs58415480) resides within intron 107 of Spectrin Repeat Containing Nuclear Envelope Protein 1 ( SYNE1 ), 130 kb downstream of ESR1 , the latter being the only gene that resides completely within the topologically associating domain (TAD; Appendix 1—figure 1 ). While the role of estrogen in leiomyomagenesis has been firmly established, this is the first genetic evidence to this end. The lead SNP at 2 p resides in the third exon of the gene GREB1. GREB1 is an essential regulatory factor of ESR1 ( Mohammed et al., 2013 ). WT1 , WNT4 and FOXO1 are central factors in uterine development and in the preparation for pregnancy (decidualization) in endometrium ( Biason-Lauber and Konrad, 2008 ; Hill, 2018 ; Kaya Okur et al., 2016 ; Tamura et al., 2017 ). Perturbations in their function are known to have neoplastic potential. The strongest association at 11p13 (rs10835889) is 40 kb downstream of the closest gene WT1 , at a region with enhancer activity ( Appendix 1—figure 1 ). WT1 is a transcription factor that acts as both a tumor suppressor and an oncogene ( Yang et al., 2007 ). The lead SNP at 1p36.12 (rs2235529) resides at the second intron of WNT4 . The risk allele is associated with suggestive upregulation of WNT4 ( Figure 5C ). WNT4 is known to be overexpressed in uterine leiomyomas with MED12 mutations ( Markowski et al., 2012 ), and knock-down of MED12 in UL cells reduces WNT4 expression ( Al-Hendy et al., 2017 ). The risk locus in 1p36.12 was also associated with several meQTLs suggesting that methylation may have a role in WNT4 regulation ( Figure 5 ). WNT4 encodes a signaling protein that has a crucial role in sex-determination ( Vainio et al., 1999 ), and the WNT signaling pathway has a well-established role in various malignancies such as breast and ovarian cancer ( Peltoketo et al., 2004 ). Of note, recent GWAS on gestational duration suggested that binding of the estrogen receptor at WNT4 is altered by rs3820282 (r 2  = 0.92 with our lead SNP) ( Zhang et al., 2017 ). Both WNT4 and FOXO1 are decidualization markers regulated by ESR1 ( Kaya Okur et al., 2016 ). Though these considerations support WNT4 as a candidate predisposition gene at this locus, the near-by CDC42 has been shown to play a role in uterine pathology, in particular endometriosis ( Powell et al., 2016 ), and should not be overlooked in further work. ( A ), Hi-C, TADs and CpG methylation around the locus with an 1 Mb flank. The needle plot shows the meQTL associations (dashed lines at 10% FDR; green lines denote the two SNPs, rs2235529 and rs2092315; gray ticks denote all CpGs tested; blue needle for positive coefficient, red for negative coefficient) for tumors (above x-axis; n AA  = 40, n AB  = 15, n BB  = 1 for rs2235529 and n AA  = 32, n AB  = 23 for rs2092315) and normals (below x-axis; n AA  = 23, n AB  = 12 and n AA  = 17, n AB  = 17, n BB  = 1). ( B ), Methylation differences in tumors (n = 56) at CpG chr1:22456326 by SNP rs2235529. ( C ), WNT4 expression differences in tumors (n = 41) stratified by the rs12042083 genotype. B is the risk allele, and the P-value is corrected for local multiple testing (permutation based test). Also MED12 has been implicated in uterine development in a mouse model ( Wang et al., 2017a ). DMRT1 is a transcription factor associated with male sex-development ( Lindeman et al., 2015 ). CD44 is a plausible fibroid stem cell marker ( Mas et al., 2015 ). Mutations in SALL1 and a deletion at the GWAS signal have been associated with Townes-Brocks syndrome, a condition associated with kidney malformations ( Stevens and May, 2016 ). Thus genes involved in genitourinary development are strikingly associated with UL predisposition. ATM , TP53, TERT , TERC and OBFC1 could be involved in uterine neoplasia predisposition through genetic instability and telomere maintenance. The lead SNP at 11q (rs141379009) resides in the 22nd intron of ATM , and the SNP at 17 p in the 3′-untranslated region of TP53. ATM and TP53 are involved in DNA damage response ( Guleria and Chandna, 2016 ), and they are among the relatively few genes that have been found to be recurrently mutated in leiomyosarcoma ( Lee et al., 2017 ). TERT and TERC encode subunits of the telomerase enzyme, which guards chromosomal stability by elongating telomeres ( Blasco, 2005 ). In addition OBFC1 has been associated with telomere maintenance ( Lee et al., 2013 ). TERT is expressed in germ cells as well as in many types of cancers ( Blasco, 2005 ). The neoplasia predisposing effect of the risk alleles at the TERT locus (rs72709458; rs2736100; rs2853676) has been overwhelmingly documented ( Appendix 1—table 11 ). Previous studies have reported contradicting observations on the effect of rs2736100 on telomere length ( Liu et al., 2014 ; ENGAGE Consortium Telomere Group et al., 2014 ; Lan et al., 2013 ; Melin et al., 2012 ; Choi et al., 2015 ). ULs have been shown to display shortened telomeres ( Rogalla et al., 1995 ; Bonatz et al., 1998 ), potentially provoking chromosomal instability as the lengths of chromosome telomeres are diminished. In our patient cohort, the risk allele at TERT (rs2736100) is significantly associated with shorter telomere length ( Appendix 1—figure 12 ), whereas the combined effect of SNPs at TERT , TERC and OBFC1 did not reach statistical significance. GRS associated merely with a susceptibility to the most common UL subtype, MED12 mutation positive tumors. Indeed it has been known that MED12 -mutation-positive tumors do not distribute randomly among patients ( Mäkinen et al., 2011 ), and our data provide at least a partial explanation to this intriguing finding. An outstanding susceptibility locus was identified 250 kb upstream of MED12 : our in-house patient cohort - together with a mutation-screening of their 1481 tumors - revealed that the risk allele could facilitate selection of somatic MED12 mutations. It may be that environmental factors contribute more significantly to genesis of MED12 wild-type lesions. In our recent study this tumor type was associated with a history of pelvic inflammatory disease, and thus infectious agents could be one underlying factor ( Heinonen et al., 2017 ). Obviously, also the power of GWAS to detect genetic associations to rare UL subtypes – such as the HMGA2 overexpressing or FH deficient subtypes – is reduced. This work highlights several new genetic cornerstones of UL formation, highlights genitourinary development and maintenance of genomic stability as key processes associated with it, and represents another step towards a much-improved understanding of its molecular basis. The proposed risk score can stratify the female population to low and high-risk quartiles that differ by two-fold in their UL risk. The population-specific risk score was inflated towards the African and Caribbean cohorts, which connects the predisposition loci to the excess UL prevalence in these ethnicities. While the increased risk appears minor on an individual level, the population-level burden to women’s health arising from these risk loci is highly significant considering the incidence of the condition. Together with the recent progress in molecular tumor characterization and subclassification, the identification of the genetic components of UL predisposition should pave the way towards more sophisticated prevention and management strategies for these extremely common tumors. The risk SNP with the most immediate potential value is that at estrogen receptor alpha, and our findings should fuel much further work on the interplay between individual germline genetics, endogenous and exogenous hormonal exposure, and occurrence and growth rate of UL.

Figure 1 provides an outline of the four stages that were implemented. The discovery stage was conducted with UK Biobank resources (UKBB; project #32506; accessed on April 10, 2018). The resource included pre-imputed genotypes (version 3; March 2018) for a total of 487,409 samples (486,757 samples for the X chromosome) and 96 million SNPs. The background information on the imputation and data quality control (QC) can be found through the UKBB documentation ( www.ukbiobank.ac.uk ). The UL cases were identified on the basis of both the self-reported uterine leiomyoma (UL) phenotype (UKBB data-field 20002: Non-cancer illness code 1351) and International Classification of Diseases (ICD) codes (data-fields 41202 – 41205: Main and secondary diagnosis for ICD10 code D25 and ICD9 code 218). These phenotype data resulted in a total of 20,106 UL cases prior to any sample/genotype QC. Sample QC was based on the UKBB annotation as follows. In total 409,692 samples passed the initial QC on ethnic grouping (UKBB data-field 22006): self-identified as ‘White British’, and similar genetic ancestry based on a principal component analysis (PCA) of the genotypes. Further sample QC excluded excess kinship (field 22021; 408,797 samples passed), sex-chromosome aneuploidy (field 22019; 408,241) and inconsistent gender (fields 31 and 22001, and one male with self-reported ULs; 408,081). In total 15,453 UL cases and 392,628 population-matched controls (205,157 females and 187,471 males) passed all these criteria. Raw genotype calls (UKBB version 2; Affymetrix UK BiLEVE Axiom, or Affymetrix UKBB Axiom array) were available for 805,426 SNPs: after filtering out low genotyping rate (<95%), Hardy-Weinberg equilibrium (p<10 −10 ) and minor allele frequency (MAF) <0.001, the remaining 611,887 autosomal genotypes were used to train the mixed model for association testing. Imputed SNPs with MAF <0.001 and imputation score (INFO) <0.3 were excluded. Further SNPs were excluded due to imputation panel differences between cohorts, and the remaining 8.3 million SNPs (Haplotype Reference Consortium, HRC1.1 panel) were tested for case-control association with BoltLMM (version 2.3.2) ( Loh et al., 2015 ). The default linear, infinitesimal mixed model was used to adjust for any underlying population structure. The model included categorical covariates for the 22 UK Biobank assessment centres and two genotyping arrays. The second stage meta-analysis utilized the genome-wide summary statistics from UKBB and the Helsinki cohort of 457 UL cases and 15,943 controls. Details on the Helsinki cohort’s imputation, sample and genotype QC are given in the Supplementary Methods. A total of 8.3 million SNPs passed imputation QC and were utilized in the meta-analysis with PLINK (version 1.90b3i) ( Chang et al., 2015 ). The SNPs were tested for association in six independent cohorts: Northern Finland Birth Cohort (NFBC) and five non-overlapping subsets of UKBB. In addition to the single-SNP association tests, a polygenic risk score ( Abraham and Inouye, 2015 Abraham and Inouye, 2015 ) was compiled as follows. The genomic risk score (GRS) was computed as a sum over SNP dosages weighted by their observed log-odds: LD pruning (r 2  ≤0.3) was applied in the order of UKBB association, and the remaining, genome-wide significant SNPs were chosen for the GRS. The log-odds weights were taken from the UKBB statistics (i.e. logarithm of the dosage-based ORs). The resulting GRS model was evaluated using R (3.3.1) and the packages PredictABEL (1.2 – 2) and MASS (7.3 – 45). The Northern Finland Birth Cohort (NFBC) had in total 459 UL cases and 4943 controls; details of the imputation, sample and genotype QC are given in the Supplementary Methods. Five non-overlapping, self-reported population-strata were available from UKBB (data-field 21000) and could be utilized as an independent replication: the five self-reported ancestries were ‘Black African’, ‘Black Caribbean’, ‘Indian’, ‘White Irish’ and ‘Other white background’. Sample QC excluded excess kinship (field 22021), sex-chromosome aneuploidy (field 22019) and inconsistent gender (fields 31 and 22001). The numbers of cases and controls that passed the sample QC can be found in Figure 1 . A summary of background variables is given in Appendix 1—table 1 . These five sample subsets did not overlap with the discovery GWAS individuals. A collection of ancestry-informative genotypes was utilized to assess the genetic homogeneity of each of the self-reported ancestry (details in Supplementary Methods). Our in-house patient and tumor data were investigated regarding the identified risk loci. All tumors of ≥1 cm diameter had been harvested and stored fresh-frozen (details in Supplementary Methods). MED12 mutations were screened by Sanger sequencing the MED12 exons 1 and 2 and their flanking sequences (60 bp) from all uterine leiomyoma and matching normal myometrium samples ( Mäkinen et al., 2011 ; Heinonen et al., 2014 ). The resulting sequence graphs were inspected manually and with Mutation Surveyor software (Softgenetics, State College, PA). Clinical patient data was available for the number of ULs, menopause status, parity, body mass index (BMI) and age at hysterectomy ( Appendix 1—figure 2 ). This study was conducted in accordance with the Declaration of Helsinki and approved by the Finnish National Supervisory Authority for Welfare and Health, National Institute for Health and Welfare (THL/151/5.05.00/2017), and the Ethics Committee of the Hospital District of Helsinki and Uusimaa (HUS/177/13/03/03/2016). For the cis expression quantitative trait loci (cis-eQTL) analysis, genes with less than six reads in over 80% of the samples were filtered out. The between sample normalization was done with Relative Log Expression (RLE) normalization and each gene was inverse normal transformed. The eQTL analysis was run with FastQTL (version 2.184) ( Ongen et al., 2016 ) separately for 60 tumors and 56 patient-matched unaffected, adjacent myometrium samples using permutation approach. The permutation parameter was set to ‘1000 10000’. Sequencing batch was used as a covariate. The cis-region was set to be 2 Mb. FDR correction was applied for tumors and matched myometrium separately. DNA methylation was studied in 56 tumors and 36 matched myometrium samples. The methylation calls were analyzed with bsseq (version 1.12.2) ( Hansen et al., 2012 ). Only the methylation in CpG context was considered. Every locus was required to have the coverage of ≥2 in at least 90% of samples. The association between methylation and genotype was studied with MatrixEQTL (version 2.1.1) using a linear regression model ( Shabalin, 2012 ). The LD-independent (r 2  ≤0.3) SNPs from the discovery stage ( Appendix 1—table 2 ) and meta-analysis ( Appendix 1—table 3 ) were considered (altogether 57 SNPs). The SNPs with MAF <0.05 in the methylation samples were filtered out. This resulted in 44 SNPs in tumors and 45 SNPs in matched myometrium. Cis methylation quantitative trait loci (cis-meQTL) was determined to be within 1 Mb flank from the SNP of interest. To annotate the CpGs with genomic context, the overlap between UCSC’s gene track (hg19) and known CpG islands was studied. As the role of promoter methylation is well known, promoter methylation was studied in addition to gene body methylation. Core promoter was defined as a region −2 kb and +1 kb from the transcription start site. The methylation analysis was performed separately for tumors and matched normal myometrium to study whether the changes in methylation could be observed in both tissues. The whole genome sequenced (WGS) samples, in total 71 tumors (48 Illumina, 23 Complete Genomics) and 51 matched myometrium samples (28 Illumina and 23 Complete Genomics), were prepared following Illumina and Complete Genomics protocols and processed as described previously ( Mehine et al., 2013 ). Structural variation was defined as a structural rearrangement (e.g deletion, inversion or translocation) not detectable in matched normal myometrium. Structural variation was detected as described in Mehine et al. ( Mehine et al., 2013 ) The mean telomere length was estimated for Illumina samples using Computel (version 0.3) ( Nersisyan and Arakelyan, 2015 ) with the default settings. Clonally related tumors were excluded from the analysis by randomly sampling one tumor to represent each clonally related tumor group. Clonally related tumors had identical changes in driver genes and shared at least a subset of somatic copy-number changes and/or copy neutral loss of heterozygosity (see Mehine et al. ( Mehine et al., 2015 ) for further details in identification of the clonally related tumors). Kruskal-Wallis test was used to assess the telomere length differences between tumors and matched myometrium as well as between genotypes. Linear model was used to calculate the association between number of risk alleles and telomere length. Pathway enrichment of all genome-wide significant SNPs was tested with DEPICT (version 1 release 194) following the default settings ( Pers et al., 2015 ). The tool is designed to integrate multiple GWAS loci for in silico target gene prioritization, pathway enrichment and tissue-specific expression profiling. In short, the DEPICT framework combines phenotype-free co-expression networks, predefined pathways and protein-protein interaction networks in order to reveal functionally connected genes among the multiple risk loci. The tool is restricted to autosomal SNPs. Meta-analysis was implemented with an inverse-variance weighted, fixed effect model. Associations between the risk alleles and other variables were tested assuming an additive genotype model unless otherwise noted. The DHARMa (0.1.5) package was applied to evaluate the goodness-of-fit of the binomial and negative binomial models. The contribution of GRS to prevalence was estimated by [(E a /P i -1)/(P a /P i -1)], where E a  = P i *GRS a /GRS i assumes a linear relationship between GRS and the true risk, and P x and GRS x are the population-specific prevalence and mean GRS, respectively. All statistical tests were two-tailed unless otherwise noted. Summary statistics were collected from each of the study stages and are available as Appendix 1—table 2 and Supplementary file 1 – 3 . For each SNP, we report its allele frequency, effect size estimates and association based on the default linear, infinitesimal mixed model. For the meta-analysis stages, we also report the Cochrane’s Q statistic and I 2 heterogeneity index in addition to the fixed-effects meta-analysis association and effect size (random-effects meta-analysis is included for reference). BoltLMM reported lambda (λGC) 1.055, 1.045 and 1.016 for UKBB, Helsinki and NFBC, respectively. The X chromosome associations were processed separately and included only the female controls (λGC 1.052, 1.049 and 1.005 for UKBB, Helsinki and NFBC, respectively). For GWAS, p<5 × 10 −8 was reported as significant. The GRS association tests ( Appendix 1—table 6 ) were controlled for family-wise error rate (FWER) and reported significant for Holm-Bonferroni adjusted p<0.05. Large families of association tests were controlled for false discovery rate (FDR; Benjamini-Hochberg method) and noted significant at FDR < 10%. In the six telomere length association tests and the two structural variation association tests, p<0.05 was considered statistically significant.