Genome-wide causal mediation analysis identifies genetic loci associated with uterine fibroids mediated by age at menarche.

In this study, we recruited individuals in two TWB cohorts with 13 899 females in TWB1 and 55 653 females in TWB2. Two distinct sets of individuals participated in the two cohorts, except for 730 individuals who were enrolled in both cohorts. Following quality control processes and removal of samples with missing data, 13 437 and 53 003 individuals were finally included in the TWB1 and TWB2 cohorts, respectively. DNA samples from both cohorts were genotyped using two different customized chips. The TWB1 custom array, which was released in April 2013, was designed for Taiwan’s Han Chinese using the Axiom Genome-Wide Array Plate System. Based on the experience with TWB1, the TWB2 custom array was released in August 2018 to cover SNPs specific to the Taiwanese population. Briefly, the TWB1 and TWB2 arrays included 653 291 and 714 431 SNPs, respectively. Approximately 100 000 of SNPs overlapped between the TWB1 and TWB2 arrays. Considering that only 15.4% of TWB1 SNPs and 14.6% of TWB2 SNPs were in overlap regions, the study subjects genotyped by two different platforms were not be used in a replication analysis but in a complementing GWAS. The demographic and clinical characteristics of the TWB1 and TWB2 cohort participants are summarized in Supplementary Table SII , showing that the overall prevalence of uterine fibroids in the Taiwanese population is ∼21% (TWB1: 19.6%, TWB2: 21.4%). This table shows that subjects in TWB1 were slightly older, had an earlier at age menarche, and had lower education levels than subjects in TWB2. Lifestyle characteristics (i.e. marital status, alcohol drinking, smoking and exercise) of the TWB1 and TWB2 cohorts are highly similar, indicating the homogeneity of risk factors in these two populations. Therefore, we treated the TWB1 and TWB2 cohort samples as two sets of complementing samples from an identical population with similar lifestyles. SNPs identified in one cohort were regarded as complement SNPs to the other cohort. Association analysis and causal mediation analysis were subsequently conducted to identify genetic loci associated with uterine fibroids mediated by age at menarche. We adopted an imputation method for predicting unobserved genotypes. To minimize the imputation bias, we performed genotype imputation only for the SNPs we identified in the following initial step of association analysis rather than for the original set of SNPs. Prior to causal mediation analysis, the following three analyses were conducted to confirm the plausibility of the mediation model between SNPs and uterine fibroids: (i) examining the effect of age at menarche on uterine fibroids, (ii) screening potential risk SNPs for age at menarche and (iii) screening potential risk SNPs for uterine fibroids. In the first analysis, we evaluated the effect size of age at menarche on uterine fibroids in the respective cohorts using logistic regression models adjusting for independent variables (i.e. education, marriage, alcohol drinking, smoking, age, exercise, and BMI). The main results of logistic regression are presented in Supplementary Table SIII . The results of logistic regression show that later age at menarche is significantly associated with a decreased risk of uterine fibroids (adjusted odds ratios (OR adj ) = 0.909 with P  <   0.05 in both cohorts). In addition, we applied G-estimation, a semiparametric approach ( Dukes and Vansteelandt, 2018 ), to verify the causal effect of abnormal timing of menarche (16 years old) on uterine fibroids. Education and birth year are considered in G-estimation as potential confounders in the relationship between age at menarche and uterine fibroids. While birth year itself should not directly affect the age at menarche, it is regarded as a valid surrogate of exposure to the macro-environment, such as early-life events, health and nutritional intake, which are risk factors for early and late age at menarche. The remaining independent variables, including marriage, alcohol drinking, smoking, exercise and BMI, were measured after menarche, and they would not confound the causal interpretation. More specifically, these can be treated as potential mediators between age at menarche and uterine fibroids. Since controlling possible mediators in analysis results in a biased result for causal inference, we performed G-estimation for estimating causal effects without adjusting for marriage, alcohol drinking, smoking, exercise and BMI. As a result, the causal odds ratios were 0.769 with P  =   0.029 in TWB1 and 0.672 with P  <   0.001 in TWB2. These findings are consistent with those of previous epidemiological studies that have shown elevated risk for developing uterine fibroids in women with early age at menarche ( Schwartz, 2001 ; Laughlin et al. , 2010 ; Velez Edwards et al. , 2013 ). In the second analysis, a GWAS was performed to identify potential risk SNPs associated with age at menarche. An additive linear regression model adjusting for independent variables was applied to GWAS for age at menarche. The independent variables in this analysis were education, marriage, alcohol drinking, smoking, age, exercise and BMI. The normality assumption of age at menarche was checked using a quantile–quantile plot (Q–Q plot) and a density plot for huge sample sizes (shown in Supplementary Fig. S1 ). Although the Q–Q plots do perfectly follow the line y  =  x and show light-tailed distributions, a violation of the normal distribution assumption in a linear model does not impact the regression estimates in studies with large sample sizes ( Schmidt and Finan, 2018 ). Given that the TWB1 and TWB2 cohort DNA samples were genotyped using two different arrays, the aforementioned analyses and results of the respective arrays were reported separately. In total, 17 and 142 SNPs were associated with age at menarche in the TWB1 and TWB2 cohorts, respectively, at a false discovery rate (called q -value) of 5% (corresponding to P  =   8.9 × 10 −7 in TWB1 and P  =   1.4 × 10 −5 in TWB2) ( Fig. 3a and c ). Complete SNP lists are provided in Supplementary Table SIV . These SNP markers were considered in the following mediation analysis. To assess the plausibility of the identified SNPs, we further highlighted SNPs that showed stronger associations with age at menarche by setting a q -value of <1 × 10 −3 (shown in Table I ) and subsequently examined the consistency between the most strongly associated SNPs and previous findings in the literature. Overall, the SNPs in Table I clustered in 14 loci: INADL , EXTL2P1 , RIPK1 , KRT19 , CASC17 , RXRG , GPR45 , LINC00577 , LIN28B , ZNF483 , GAB2 , CHEK1 , BRCA2 and RP11-359N5.1 . Among these 14 identified loci associated with age at menarche in the Taiwanese population, five ( RXRG , GPR45 , LIN28B , ZNF483 and GAB2 ) were previously reported for different populations ( Elks et al. , 2010 ; Delahanty et al. , 2013 ; Tanikawa et al. , 2013 ; Day et al. , 2017 ). Manhattan plot showing the -log 10 ( P values) for individual variant association with age at menarche and uterine fibroids. Results are from genome-wide association studies in TWB1 ( a, b ) and TWB2 ( c, d ). Age at menarche is shown in blue dots in a, c and uterine fibroids is shown in green dots in b, d. The red lines indicate -log 10 ( P values) at false discovery rate of 5%: P  =   8.93 × 10 −7 (a: age at menarche in TWB1), P  =   7.71 × 10 −7 (b: uterine fibroids in TWB1), P  =   1.39 × 10 −5 (c: age at menarche in TWB2) and P  =   4.20 × 10 −6 (d: uterine fibroids in TWB2). TWB, Taiwan Biobank. SNPs associated with age at menarche reaching the genome-wide significant threshold at a q -value of <1 × 10 −3 . β adj , adjusted beta coefficient; MAF, minor allele frequency; P , P value; SNPs, single-nucleotide polymorphisms; TWB, Taiwan biobank. In the third analysis, we identified potential risk SNPs associated with uterine fibroids. We adopted an additive logistic regression model by adjusting for independent variables in the GWAS for age at menarche and added age at menarche as an independent variable. The Manhattan plots of the GWAS results for uterine fibroids are shown in Fig. 3b and d . At a q -value of <0.05 (corresponding to P  =   7.7 × 10 −7 in TWB1 and P  =   4.2 × 10 −6 in TWB2), we identified 10 and 50 uterine fibroid-related SNPs in the TWB1 and TWB2 cohorts, respectively ( Supplementary Table SV shows the complete list of SNPs). Furthermore, we selected the SNPs with the strongest association with uterine fibroids at a q -value of <1 × 10 −3 ( Table II ). Notably, we could not identify any significant association with uterine fibroids under this strict threshold in the TWB1 cohort. Consequently, the SNPs in Table II represent 13 loci: LINC00339 , CDC42 , WNT4 , GREB1 , REXO4 , OBFC1 , SLK , COL17A1 , SFR1 , BET1L , MST4 , RAP2C and RAP2C-AS1 . Seven of the 13 loci we identified ( WNT4 ( Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Gallagher et al. , 2019 ; Sakai et al. , 2020 ), GREB1 ( Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Gallagher et al. , 2019 ; Sakai et al. , 2020 ), SLK ( Sakai et al. , 2020 ), CDC42 ( Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Gallagher et al. , 2019 ), OBFC1 ( Cha et al. , 2011 ; Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Gallagher et al. , 2019 ), BET1L ( Cha et al. , 2011 ; Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Gallagher et al. , 2019 ) and RAP2C ( Välimäki et al. , 2018 ; Gallagher et al. , 2019 )) were previously associated with uterine fibroids. In addition, LINC00339 (long intergenic non-protein coding RNA 339) and COL17A1 were previously associated with the endometrium. SNPs associated with uterine fibroids reaching the genome-wide significant threshold at a q -value of <1 × 10 −3 . No SNP was found with significant association with uterine fibroids under the strict threshold at a q -value of <1 × 10 −3 in TWB1 cohort. MAF, minor allele frequency; OR adj , adjusted odds ratio; P , P value; SNPs, single-nucleotide polymorphisms; TWB, Taiwan biobank. Two sets of risk SNPs associated with age at menarche and uterine fibroids were identified using GWAS. Next, we examined which SNPs in these sets play a role in the causal mechanism between age at menarche and uterine fibroids. Based on the finding that age at menarche is a risk factor for uterine fibroids, we treated age at menarche as the mediator of the association of genetic factors with uterine fibroids. The causal diagram is presented in Fig. 1 . The methodology of causal mediation analysis enables us to assess the effect of exposure on the outcome, which is mediated through the mediator ( Imai et al. , 2010 ; Valeri and VanderWeele, 2013 ; VanderWeele, 2016 ). In causal mediation analysis, satisfying exchangeability is one of the key assumptions to ensure that the result has a causal interpretation. The exchangeability assumption requires a comprehensive collection of common causes in exposure–mediator, exposure–outcome and mediator–outcome relationships to avoid confounding effects. Therefore, we included education and birth year as confounders in the analysis. The effect of an individual SNP on uterine fibroids is the TE. As mentioned above, TE could be separated into NDE and NIE. For easier understanding, we used DE and indirect effect (IE) to represent NDE and NIE, respectively, throughout this study. DE represent the effects of SNPs unmediated by age at menarche on uterine fibroids, whereas IE represent the effects of SNPs mediated by age at menarche on uterine fibroids. In other words, TE is given by DE and IE. For decomposition, we used the R package ‘mediation’ ( Tingley et al. , 2014 ), which uses a model-based approach to estimate causal mediation effects based on the counterfactual model ( Rubin, 1974 ; Little and Rubin, 2000 ). Accordingly, we could estimate the risk difference scale for each SNP associated with either age at menarche or uterine fibroids. In the TWB1 cohort, 27 (17 from the aforementioned GWAS analysis for age at menarche + 10 from the aforementioned GWAS analysis for uterine fibroids) SNPs in were considered in the following mediation analysis, while there were 192 SNPs (142 + 50) identified in the TWB2 cohort. Notably, the menarche-associated SNPs and fibroid-associated SNPs identified in this study are generally disjoint in both cohorts (only six common SNPs were identified). We found 162 SNPs (6 common SNPs + 12 TWB1-array specific SNPs + 144 TWB2-array specific SNPs) associated with both uterine fibroids and age at menarche ( Fig. 2 ). The association between these 162 SNPs and uterine fibroids was significantly mediated by age at menarche ( P  <   0.05) ( Supplementary Table SVI ). We used rs809302 to explain the findings of TE, DE, and IE. Our analysis showed that, overall, there were approximately four excess cases of uterine fibroids per 100 subjects in the group with rs809302, compared to those of the group without this SNP (TE = 0.0397 in Supplementary Table SVI ). Thereafter, we dissected the effect of this SNP on uterine fibroids into DE and IE. The positive DE value indicates that subjects with this SNP are more likely to have uterine fibroids irrespective of their age at menarche. The positive IE value indicates that subjects with this SNP are more likely to have menarche at earlier ages, leading to an increased risk of developing uterine fibroids. In the case of rs809302, a DE of 3.92% and an IE of 0.05% were found, indicating that among the four excess cases of uterine fibroids per 100 subjects, 3.92 could be attributed to DE and only 0.05 could be attributed to IE. We were able to calculate the proportion of the DE of risk SNPs on uterine fibroids and the age at menarche-mediated effect of risk SNPs. In this case, the mediation effect accounted for only 1.3% (0.05/4) of the mechanism. We followed the above criteria to classify our identified SNPs into four groups according to DE and IE on uterine fibroids: (i) positive DE and positive IE ‘both-harmful’ group, (ii) positive DE and negative IE ‘mediator-protective’ group, (iii) negative DE and positive IE ‘mediator-harmful’ group and (iv) negative DE and negative IE ‘both-protective’ group ( Fig. 4a and c ). SNPs in the ‘both-harmful’ group increased the risk of developing uterine fibroids. rs750648715 ( KRT19 ) (DE = 1.74%, P  <   0.001; IE = 0.24%, P  =   0.013) and rs59037454 ( AL162725.2 ) (DE = 1.54%, P  <   0.001; IE = 0.17%, P  <   0.001) were reported in TWB1 with relatively large effect size (i.e. SNPs located in the right-upper quadrant of Fig. 4a ). In TWB2, rs809302 ( SLK ) was identified as the most representative risk SNP (DE = 3.92%, P  <   0.001; IE = 0.05%, P  =   0.067) (i.e. SNPs located in the right-upper quadrant of Fig. 4c ). In contrast, the ‘both-protective’ group includes SNPs that are associated with a reduced risk of developing uterine fibroids. The top two SNPs in TWB1 were rs7755818 ( RIPK1 ) (DE = −1.5%, P  =   0.02; IE = −0.16%, P  <   0.001) and rs72869451 ( CASC17 ) (DE = −1.2%, P  <   0.001; IE = −0.22%, P  <   0.001) (i.e. SNPs located in the left-lower quadrant of Fig. 4a ). In TWB2, the top three SNPs of the ‘both-protective’ group were rs587776903 ( ZNF644 ) (DE = −8.81%, P  =   0.040; IE = −0.69%, P  =   0.013), rs587784142 ( NSD1 ) (DE = −9.5%, P  =   0.053; IE = −1.12%, P  <   0.001) and rs73069835 ( RP11-761N21.1 ) (DE = −6.42%, P  =   0.08; IE = −1.04%, P  <   0.001) (i.e. SNPs located in the left-lower quadrant of Fig. 4c ). The ‘mediator-harmful’ and ‘mediator-protective’ groups contained SNPs that displayed different directions in terms of DE and IE on uterine fibroids. The ‘mediator-harmful’ group contained SNPs that were associated with reduced fibroid risk; however, they were associated with early menarche, which subsequently increased the probability of developing uterine fibroids. In TWB1, there was no SNP in the ‘mediator-harmful’ group (i.e. SNPs located in the left-upper quadrant of Fig. 4a ). In TWB2, although some SNPs were classified into this group, the effects were relatively small (i.e. SNPs located in the left-upper quadrant of Fig. 4c ). The ‘mediator-protective’ group was composed of SNPs that were associated with increased fibroid risk; however, they were associated with late menarche, which may reduce the harmful effect of an SNP on uterine fibroids (i.e. SNPs located in the right-lower quadrants of Fig. 4a and c ). rs314272 ( LIN28B ) was reported in both TWB1 and TWB2 in the ‘mediator-protective’ SNP group. In addition, rs80357190 ( BRCA1 ) (DE = 1.56%, P  =   0.867; IE = −1.81%, P  <   0.001) and rs371721345 ( HLA-DOB ) (DE = 1.7%, P  =   0.967; IE = −2.70%, P  <   0.001) were detected in TWB2. Scatter plots. Left: scatter plot between direct and indirect effects ( a ) for TWB 1 and ( c ) for TWB2. X -axis and Y -axis represent the estimates of the direct and indirect effects separately. The color indicates the absolute value of proportion mediated (|PM|), which is defined as the ratio of the indirect effect to the sum of direct effect and indirect effect. The dashed line represents y + x = 0 . Right: scatter plot between age at menarche and the predicted risk difference for top SNPs ( b ) for TWB1 and ( d ) for TWB2. The X-axis represents age at menarche ranging from 7 to 25 years old. Y -axis represents the risk difference of developing uterine fibroids at a specific age at menarche. TWB, Taiwan Biobank. We considered the aforementioned SNPs as high-impact factors for uterine fibroids. Thereafter, we evaluated their effects on uterine fibroids at a specific age at menarche using CDE. Unlike DE, which is used to explore the natural mechanism of exposure on the outcome without intervening in the mediator, CDE is typically used to assess how exposure affects the outcome when the mediator had a specific value ( VanderWeele, 2011 , 2013 ). The results of CDE analysis on the risk difference scale are illustrated in Fig. 4b and d . The CDE analysis revealed that the absolute value of the risk change is decreased along with the age of menarche for all SNPs that we identified. Such a finding implies that the impact of a genetic factor on the risk of uterine fibroids is generally decreased with later menarche ( Supplementary Table SIII ). Assessing the risk change for each SNP under different ages at menarche provides insights into the dynamic relationship between age at menarche and uterine fibroids. For example, as shown in Fig. 4d , the risk of developing uterine fibroids increased by ∼5% for women with rs809302 ( SLK ) if they experienced early menarche (≤10 years old). In addition, our analysis further reported that the estimated risk changes for three SNPs (i.e. rs73069835 ( RP11-761N21.1 ), rs587776903 ( ZNF644 ) and rs587784142 ( NSD1 )) always presented a negative association with uterine fibroids regardless of early or late age of menarche, indicating that these three SNPs would inhibit the occurrence of uterine fibroids. We developed a hierarchical Bayesian logistic regression model with Polya-Gamma sampling to estimate the SNP heritability of uterine fibroids. As a result, we estimated the SNP heritability to be 6.9% in the TWB1 cohort. This indicates that 6.9% of the variance in uterine fibroids within the TWB1 cohort is due to the 18 identified SNPs (=6 common SNPs + 12 TWB1-array specific SNPs). In the TWB2 cohort, we identified 150 SNPs (=6 common SNPs + 144 TWB2-array specific SNPs), and the estimate of SNP heritability was 30.1%. It has been known that many reproductive traits in women have high heritability ( McGrath et al. , 2021 ). For uterine fibroids, previous estimates of heritability in the European population were 54.8% in a familial aggregation study ( Luoto et al. , 2000 ) and 69% in a twin study ( Snieder et al. , 1998 ). Note that the SNP heritability estimated in this study is the proportion of variance in fibroid risk in the population explained by genetic variation when the SNP effects on uterine fibroids were mediated by age at menarche. If we assumed the SNP heritability for uterine fibroids as 54.8% (like the estimate in the European population), the result showed that the percentages of SNP heritability attributable to age at menarche account for 12.6% of heritability (=6.9/54.8) in TWB1 and 54.9% of heritability (=30.1/54.8) in TWB2. In total, mediation analysis reported 162 distinct SNPs representing 77 transcripts. To confirm the plausibility of our analysis, we additionally evaluated the predictive ability of our genomic strategy based on these 77 transcripts using microarray data for uterine leiomyosarcomas, available in the gene expression omnibus database ( GSE64763 ) ( Barlin et al. , 2015 ). The dataset had 79 samples: 29 from normal myometrium, 25 from fibroids and 25 from uterine leiomyosarcomas, which are rare mesenchymal malignant tumors. In the uterine leiomyosarcoma dataset, RNA was hybridized to Affymetrix U133A 2.0 transcription microarrays, and a total of 22 218 transcripts were included. Among the 77 transcripts we identified, 41 were included in the microarray data. Accordingly, we examined the predictive ability of these 41 transcripts for uterine leiomyosarcomas and fibroids. A heatmap for the gene expression profiles of these transcripts generated using the R function ‘heatmap’ is presented in Fig. 5 . This function further incorporates the ‘hclust’ function in R for clustering. In this evaluation, we used complete-linkage clustering, one of the agglomerative methods of hierarchical clustering. If we separated the samples into two groups based on the information of the tree, the accuracy is ∼80% in distinguishing the case group (i.e. uterine leiomyosarcomas and fibroids) from the control group (i.e. normal myometrium). This is significantly higher than the results of random separation. Moreover, sensitivity and specificity were ∼76% and 86%, respectively. The high value of accuracy implies that the identified SNPs have a superior ability to predict the development of uterine leiomyosarcomas and fibroids. Heatmap of 41 transcripts obtained from causal mediation analysis in the relation between age at menarche and uterine fibroids. The gene expression profiles are microarray data from a previous study for uterine leiomyosarcomas. Expression values for each transcript (column) are normalized across all samples (row). Both column and row clustering were applied. To identify potential biological pathways that influence uterine fibroids by the age at menarche, we performed enrichment analysis of the 77 transcripts in which the 162 SNPs were found. Enrichment analysis has been widely used to examine whether identified transcripts are overrepresented in an a priori defined set of genes. Using ConsensusPathDB web ( http://cpdb.molgen.mpg.de/ ), at a P of <1 × 10 −3 , 13 biological pathways were identified ( Fig. 6 ). The most significant pathway was that regulating transferase activity ( P  =   6.64 × 10 −5 ). These results cover a wide and nonspecific range of biological pathways, such as regulation of molecular processes, regulation of biological processes or cellular response to stimuli. Gene Ontology (GO) enrichment analysis of 77 transcripts containing SNPs with significantly age at menarche-mediated effects on uterine fibroids. The X -axis represents -log 10 ( P value), and Y -axis represents the ratio of the identified transcript number to the total gene number in a GO term. The color corresponds to the category. The blue line indicates -log 10 ( P value) = 3. SNPs, single-nucleotide polymorphisms. Furthermore, we used GeneMANIA, an online network integration tool for predicting gene interactions, to analyze the intergenic interactions between the above-mentioned 77 transcripts ( Franz et al. , 2018 ). As the long noncoding RNA (lncRNA) genes or long intergenic noncoding RNA (lincRNA) transcripts (e.g. LIN28B-AS1 , AL162725.2 , CASC17 , LINC00339 , AC096559.1 , AC019330.1 , RP11-761N21.1 , ADAMTS9-AS2 , RP11-384F7.1 , RP11-96H19.1 , RP11-219B17.1 , RP11-429P3.5 and RP11-142G1.1 ) were absent from the database, we excluded them from this bioinformatic analysis. In addition, AC079447.1 , UGT3A2 , C2orf74 and SKOR2 were excluded from this analysis because they were isolated and did not interact with the other transcripts. Consequently, 60 identified transcripts were included in the resulting network ( Fig. 7 ). Most of these genes corresponded to physical interactions, followed by predicted interactions. Intergenic interactions comprised physical interactions (48.92%), predicted interactions (22.23%), co-expression (16.16%), co-localization (8.24%), genetic interactions (3.44%), pathway (0.65%) and shared protein domains (0.37%). GeneMANIA interaction network between the 77 transcripts identified with a causal effect on uterine fibroids mediated by age at menarche. The transcripts determined in the present study are cross-shaded. Edge thickness represents the interaction strength. The color represents the transcript interaction type as described in the figure legend.

The Taiwan Biobank (TWB) is an ongoing community- and hospital-based cohort aiming to enroll 200 000 individuals from the general Taiwanese population between 30 and 70 years old. The inclusion criteria were individuals without cancer diagnosis at the time of enrollment. A pilot study was conducted to demonstrate the feasibility of this program with 7956 subjects enrolled; among them, 6008 (75.5%) were retained in the formal cohort. The formal cohort began enrolling Taiwanese study participants since 2012, with extensive phenotypic data collected from 148 291 individuals as of May 2021. Over 99% of participants are Han Chinese. Supplementary Table SI illustrates the distributions of the self-reported parental origins of the TWB participants. Furthermore, a recent study revealed the genetic population structure of the Han Chinese populations in Taiwan based on the TWB data ( Chen et al. , 2016 ). They reported that 79.9% of Taiwanese Han Chinese were of southern Han Chinese ancestry, implying that the Taiwanese population is a homogeneous racial/ethnic group. The effective number of females with whole-genome genotyping data is 66 440. The variables for analysis for all participants included age, BMI, self-reported age at menarche, lifestyle (i.e. alcohol drinking, smoking, regular exercise), socioeconomic status (e.g. education) and marital status. Among the enrolled females, 21% reported having a medical diagnosis of uterine fibroids through the following question: ‘Have you ever been diagnosed with uterine fibroids by a physician?’. The study was approved by the Central Regional Research Ethics Committee, China Medical University, Taiwan (CRREC-109-190). A total of 653 291 single-nucleotide polymorphisms (SNPs) were genotyped in TWB1 samples using the Axiom Genome-Wide TWB 1.0 Array Plate, which was designed based on the Thermo Fisher Axiom Genome-Wide CHB Array with customized contents. Furthermore, 714 431 SNPs were genotyped in replication cohort (TWB2) samples using the Axiom Genome-Wide TWB 2.0 Array Plate. We removed 6318 and 27 992 SNPs from TWB1 and TWB2, respectively. Consequently, 646 973 TWB1 SNPs (99% genotyping rate) and 686 439 TWB2 SNPs (96% genotyping rate) passed quality control and were assessed in our genome-wide association study (GWAS) and mediation analysis. A total of 99 939 SNPs overlapped between the TWB1 and TWB2 arrays. There were 547 034 TWB1 array-specific SNPs and 586 500 TWB2 array-specific SNPs. Imputation of missing genotypes was conducted after the initial association analyses. For each individual, an imputed SNP was specified according to the probabilities of the three possible genotypes obtained from the observed population. The genotype with the maximum probability was imputed as the genotype for a missing SNP. Causal mediation analysis is a typical technique for assessing mediation when a causal effect is confirmed ( Robins and Greenland, 1992 ; Pearl, 2001 ). In this study, we treated age at menarche as a potential mediator in the pathway from SNPs to uterine fibroids. We first screened potential causal SNPs for age at menarche and uterine fibroids, separately, using causal inference. In causal inference, the population level total effect (TE) (or called an average causal effect) is widely used to assess the causal effect of exposure on disease. Since the TE is generally unidentifiable, we require exchangeability, consistency, and positivity assumptions to identify the TE from observational data ( Greenland and Robins, 1986 ). The key concept of the exchangeability assumption is to completely collect the common causes of exposure and disease to avoid confounding effects. Under these assumptions, the estimator of TE is given in the Supplementary Materials and methods . TE can be estimated by several statistical methods, such as marginal structural model, inverse-probability-weighting, standardization, and G-estimation. In this study, we adopted standardization, a regression-based approach, to identify potential causal SNPs with a significant TE. In the regression-based approach, we considered a Gaussian model for age at menarche and a logistic model for uterine fibroids. Given a confirmed causal SNP of uterine fibroids, we used causal mediation analysis to its mediation effect by age at menarche. Natural direct effect (NDE) and natural indirect effect (NIE) decomposed from TE are the measures to quantify the importance of a particular mediator in the mechanism. In our study, NDE captures the effect of an SNP on uterine fibroids through pathways that do not involve age at menarche, whereas NIE captures the age at menarche-mediated effect of an SNP on uterine fibroids. Four assumptions are necessary for identifying NDE and NIE: (i) no unmeasured exposure–outcome confounder, (ii) no unmeasured mediator–outcome confounder, (iii) no unmeasured exposure–mediator confounder and (iv) no mediator–outcome confounder affected by the exposure ( Avin et al. , 2005 ; VanderWeele and Vansteelandt, 2009 ). Under these four assumptions, the empirical expressions of NDE and NIE are given in the Supplementary Materials and methods . The conventional approach for the estimation of causal mediation analysis is a regression-based method, which entails a two-step procedure. The first step is to separately estimate the distributions of the outcome and mediator. In this analysis, the outcome (uterine fibroids) was modeled using a logistic regression model, and the mediator (age at menarche) was modeled using a linear model. In the second and final step, given the estimated distributions, Monte Carlo integration was adopted to approximate NDE and NIE. Proportion mediated (PM) is defined as the ratio of NIE to TE and is used to assess the degree to which TE operates through the mediator. The assessment of PM provides insight into the role of different pathways. In particular, if NIE and TE are both positive and negative, a large PM value implies that an SNP causes uterine fibroids primarily by age at menarche than by other factors. Controlled direct effect (CDE) is an alternative measure to assess the direct effect (DE) where a hypothetical intervention controls the mediator to a particular value ( Robins and Greenland, 1992 ; Pearl, 2001 ; VanderWeele, 2011 ). Two properties of CDE are crucial for mediation analysis. First, the assumptions required for CDE identification are considerably weaker than those of other methods. Second, CDE is of greater relevance for policy-making, as it corresponds to the effect of exposure on the outcome after the intervention to fix the mediator to the status or value of interest. We performed CDE analysis to estimate DEs on uterine fibroids at different ages at menarche. SNP heritability is the fraction of the phenotypic variance explained by a specific set of SNPs. Assessing the relative effect of genetics versus environment for uterine fibroids can facilitate us to understand the causal mechanism of disease etiology. For the SNP heritability estimation, we use a hierarchical Bayesian logistic regression to model the relationship between the outcome and SNPs. Given proper prior distributions, we developed an algorithm for parameter estimation by using Gibbs Sampling and the Metropolis–Hastings algorithm. In this study, the number of iterations is 10 000, and the length of burn-in is 0.8 times the number of iterations. We obtained the posterior distribution of the parameters and then performed an analysis using the posterior mean. The proposed algorithm is detailed in the Supplementary Materials and methods .

Using causal mediation analysis, we identified 162 SNPs in 77 transcripts associated with uterine fibroids mediated by age at menarche via four different causal mechanisms: ‘both-harmful’ (52 SNPs), ‘both-protective’ (34 SNPs), ‘mediator-harmful’ (22 SNPs) and ‘mediator-protective’ (54 SNPs). Specifically, in TWB1, we found 18 SNPs with a causal effect on uterine fibroids mediated by age at menarche (6 SNPs in the ‘both-harmful’ group, 3 in the ‘both-protective’ group, none in the ‘mediator-harmful’ group and 9 in the ‘mediator-protective’ group). In TWB2, 150 SNPs were found using causal mediation analysis (46 SNPs in the ‘both-harmful’ group, 31 in the ‘both-protective’ group, 22 in the ‘mediator-harmful’ group and 51 in the ‘mediator-protective’ group). The proportion at which age at menarche mediates the relationship between individual SNPs and uterine fibroids could be up to ±41%. This result underscores the importance of genetic variation that affects the occurrence of uterine fibroids mediated by age at menarche. For example, LIN28B and GAB2 have been frequently associated with the timing of puberty in GWAS ( Ong et al. , 2009 ; Day et al. , 2017 ; Ponomarenko et al. , 2020 ). In our mediation analysis, eight SNPs located in LIN28B (i.e. rs314266, rs314272, rs160593, rs314268, rs1744206, rs364663, rs9377684 and rs9391264) were further characterized as ‘mediator-protective’ as they delay age at menarche and subsequently protect women from uterine fibroids. A similar conclusion was reached for eight SNPs located in GAB2 (i.e. rs1385600, rs2248407, rs1981405, rs7101429, rs7107174, rs10899469, rs2373115 and rs7937277). Furthermore, our results showed that WNT4 , a plausible predisposition gene to uterine fibroids ( Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Gallagher et al. , 2019 ; Sakai et al. , 2020 ), belonged to the ‘mediator-harmful’ group. Although subjects with SNP rs3765350 located in WNT4 should have a lower risk of uterine fibroids, they may have increased risk owing to early menarche. Moreover, we obtained suggestive evidence that some risk SNPs directly lead to uterine fibroids regardless of age at menarche. As mentioned above, rs809302 ( SLK ) increases the risk of developing uterine fibroids, with 3.92 excess cases of uterine fibroids per 100 subjects through a mechanism other than age at menarche ( P  <   0.001), whereas the increased risk of rs809302 mediated by age at menarche accounts for only 0.05 cases per 100 subjects ( P  =   0.067). As previously reported, SLK was associated with uterine fibroids in a Japanese population ( Sakai et al. , 2020 ). However, that study did not explore the underlying causal mechanism among SLK , age at menarche and uterine fibroids. In contrast to SLK SNPs, other gene SNPs cause uterine fibroids, primarily mediated by age at menarche. For instance, rs371721345 in the intron of the HLA-DOB gene is associated with a 2.70% decreased risk ( P  <   0.001) in the occurrence of uterine fibroids mediated by age at menarche. However, its DE is not significant (DE = 1.7%, P  =   0.967). In addition to estimating DE and IE for each SNP, we also assessed the overall heritability of risk SNPs. We compared the heritability reported in the literature and concluded that 12.6% of overall heritability in TWB1 and 54.9% of overall heritability in TWB2 could be explained by the mediatic mechanism of age at menarche. The pronounced difference in the heritability of uterine fibroids between TWB1 and TWB2 results from an almost 9-fold increase in the SNP number (i.e. 18 SNPs were identified in TWB1 and 150 SNPs in TWB2). In addition, DNA samples in TWB1 and TWB2 cohorts were genotyped using two different SNP arrays (TWB1 array and TWB2 array). The total SNPs in TWB1 and TWB2 were 653 291 and 714 431, respectively, and only 99 939 SNPs overlapped between the TWB1 and TWB2 arrays. The two estimates of heritability cannot be compared with each other. Both results are biologically meaningful to represent the degree of variation and can complement each other. Moreover, the TWB2 array was designed based on the experience of TWB1 use to further cover more SNPs specific in the Taiwan population. Thus, a higher estimate of heritability in the TWB2 cohort was reasonable and expected. The assessment of the heritability due to mediation can help shed insight into how the overall genetic factors affect uterine fibroids mediated by age at menarche. In addition to the heritability, we evaluated the overall predictive ability of the 77 transcripts using an additional gene expression dataset, and the accuracy was reported as 80%. The high accuracy is impressive because these SNPs are primarily identified as their mediation effects are significant. It reveals the importance of age at menarche in the development of uterine fibroids caused by genetic factors. In contrast, the gene set enrichment analysis reported that no specific biological pathway is identified. While the corresponding biological pathway is unclear, the network analysis by GeneMANIA enables us to see that 77 transcripts are physically interacted with each other. Exploring the biological pathway of the transcripts identified from the mediation analysis in different datasets warrants future research. Recently, a related study has identified 14 uterine fibroid-associated loci associated with age at menarche ( Ponomarenko et al. , 2020 ). Seven out of these 14 SNPs were also included in our TWB1 and TWB2 arrays. Among these seven SNPs, three (i.e. rs2090409, rs6438424 ( LSAMP ), and rs7759938) were identified as noteworthy in our causal mediation analysis. However, the interpretations for these risk SNPs differ in two respects. First, Ponomarenko et al. (2020) aimed to assess genetic associations rather than to evaluate causality as we did. Causal inference is the primary goal for scientific research. Therefore, compared with previous research, which only identified uterine fibroid-associated loci, our study, in addition, revealed the causal ‘both-harmful’ mechanism underlying the effects of the first two SNPs and the ‘mediator-harmful’ mechanism underlying the effects of the third SNP. Second, the effects of SNPs determined using causal mediation analysis can explain how SNPs affect uterine fibroids when age at menarche is considered a mediator. Furthermore, the previously used approach for investigating SNPs associated with age at menarche and uterine fibroids lacks the ability to reveal the detailed underlying mechanism. In contrast, our study fills this research gap and elucidates the mechanism linking SNPs, age at menarche, and uterine fibroids by estimating the proportion of the SNP-associated risk mediated by age at menarche (i.e. age at menarche mediates 18–34% of the effects of these three SNPs on uterine fibroids). It is important to note that the present study identified more SNPs (162 SNPs) than the conventional GWAS method did for uterine fibroids. A conventional GWAS method can only identify SNPs with non-zero TEs (TE = IE + DE) and lacks the ability to detect an SNP with negative DE and positive IE (or positive DE and negative IE). Namely, the conventional method is typically restricted to studying the SNPs belonging to the ‘both-harmful’ or ‘both-protective’ groups. In contrast, the method of mediation analysis can completely explore the aforementioned four groups, and it would be expected to find a more comprehensive set of SNPs by using causal mediation analysis. Among our 162 menarche-mediated SNPs, over 45% of SNPs were identified in the new ‘mediator-harmful’ or ‘mediator-protective’ group. It reveals that our mediation analysis is more powerful than a GWAS analysis in terms of investigating the genetic mechanism of developing uterine fibroids. While our 162 SNPs significantly induced the development of uterine fibroids mediated by menarche in the TWB cohorts, future research should focus on a replication analysis for validation and to avoid false positives. The present study further reveals an alternative and precise way to investigate genetic pleiotropy between uterine fibroids and age at menarche. Pleiotropy occurs when a single genetic variant affects multiple traits, and these variates can be the so-called pleiotropic genes or pleiotropic SNPs. The conventional methods of pleiotropic gene/SNP identification usually rely on test statistics for association to determine the existence of a pleiotropic association between two traits. These methods can identify three types of pleiotropy: biological pleiotropy, mediated pleiotropy and spurious pleiotropy ( Solovieff et al. , 2013 ), but lack the ability to distinguish them. For example, Sakai et al. (2020) employed pleiotropic analysis based on GWAS results to identify nine uterine fibroid-associated loci associated endometriosis and other malignant tumors, and these pleiotropic effects would be induced by any type of pleiotropy ( Sakai et al. , 2020 ). In contrast, causal mediation analysis used in this study emphasizes the mediated pleiotropy which exists when one phenotype is itself causally affects a second phenotype ( Solovieff et al. , 2013 ). The present study is the first to assess the mediated pleiotropy for the relationship between age at menarche and uterine fibroids. Specifically, the estimated SNP heritability of uterine fibroids attributable to age at menarche in this study can be interpreted as its overall mediated pleiotropic effect. Inferring the mediated pleiotropic effect, instead of the effect corresponding to all types of pleiotropy, reveals more details about how genetic variants affect the causally related phenotypes. A potential limitation of this study is that it relied upon self-reported age at menarche and uterine fibroid information, which could be subject to the validity of self-reports conducted many years after these events and the consistency between self-reports and medical records. For the first issue, two cohort studies in the USA ( Must et al. , 2002 ) and the UK ( Cooper et al. , 2006 ) indicated that the correlation between original and recalled age at menarche remained high, with correlation coefficients of 0.66–0.79. A previous systematic review by Stewart et al. (2017) indicated that the source of clinical data did not substantially impact the prevalence of uterine fibroids (e.g. medical record review or self-reported). Therefore, self-reporting years later remains a practical means for collecting data on menarche and uterine fibroids. For the second issue, universal access to medical care requires medical diagnosis information. In Taiwan, national health insurance is a universal and mandatory program with coverage to 99.9% of the entire population. Although no study analyzed the consistency between self-reports and medical records for uterine fibroids in Taiwan, one Taiwanese study found good consistency for hypertension (Kappa coefficient = 0.80) and diabetes (Kappa coefficient = 0.67) among a sample of Taiwanese aged 58 and over ( Chiu et al. , 2018 ). Their sub-group analysis indicated potential under-reporting of illness among those perceived poor health conditions ( Chiu et al. , 2018 ). Since the TWB dataset did not include individuals with a cancer diagnosis at the time of enrollment, the impact of under-reporting of uterine fibroids is less in our study. Additionally, the rate of reporting a diagnosis of uterine fibroids was 21% in our study population, which fell within the rates of medical diagnosis based on national health insurance data reported in previous Taiwanese studies (i.e. 10.5% in 25 years and over (24 281 samples) ( Tseng et al. , 2017 ), 13.4% in 30–45 years old females (3 840 samples) ( Tsai et al. , 2017 ) and 50% in 40–48 years old females (5 842 samples) ( Yeh et al. , 2014 )). Therefore, data from TWB can still represent the population of eligible females. Beyond the scope of this study, a separate study could help us understand the rate of consistency in Taiwan by entailing the collection of self-reported uterine fibroids that were medically confirmed in the national health insurance database and vice versa. In conclusion, the present study is the first to adopt causal mediation analysis to investigate the genetic risk between age at menarche and uterine fibroids. We analyzed data obtained from two cohorts of 69 552 females in the population of Taiwan and reported 162 SNPs with effects mediated by age at menarche. Overall, we explicitly identified the role of these SNPs in uterine fibroids and age at menarche by comparing the DE and IE sizes. Our study demonstrates a novel approach for assessing the underlying mechanism linking SNPs, age at menarche and uterine fibroids and estimating the proportion of SNP-associated risk mediated by age at menarche. Since the risk factors for uterine fibroids not only include age at menarche, future work is suggested to explore the underlying mechanisms of multiple time-varying mediators using causal mediation analysis.

Uterine fibroids, also known as uterine myomas or leiomyomas, are common benign tumors that develop from uterine smooth muscle ( Bulun, 2013 ). Most women with uterine fibroids are asymptomatic ( Donnez and Dolmans, 2016 ). However, if accompanying symptoms occur, they include excessive menstrual bleeding, abdominal pain, infertility and anemia ( Donnez and Dolmans, 2016 ), which are associated with increased medical costs and work loss ( Hartmann et al. , 2006 ), substantially affecting the patient’s quality of life ( Soliman et al. , 2015 ). The incidence and prevalence of uterine fibroids vary widely across epidemiological studies, a prevalence ranging from 4.5% to 68.6%, depending on the diagnostic method used and population studied ( Stewart et al. , 2017 ). The Global Burden of Disease Study estimated that uterine fibroids affected 226 million females in 2019, of whom 18% were of child-bearing age (20–34 years old) and 80% were of working age (25–54 years old) ( Global Burden of Disease Collaborative Network, 2020 ). Studies estimated that 40–63% of uterine fibroids in a given population are attributed to chromosomal abnormalities or genetic factors ( Luoto et al. , 2000 ; Commandeur et al. , 2015 ). Over 30 genetic loci associated with uterine fibroids have been identified in different ethnic populations ( Cha et al. , 2011 ; Eggert et al. , 2012 ; Edwards et al. , 2013 ; Zhang et al. , 2015 ; Hellwege et al. , 2017 ; Liu et al. , 2018 ; Rafnar et al. , 2018 ; Välimäki et al. , 2018 ; Edwards et al. , 2019 ; Ponomarenko et al. , 2020 ; Sakai et al. , 2020 ). Several genetic variations in or nearby these identified loci were also candidate loci directly associated either with early age at menarche (e.g. LIN28B ) ( Ong et al. , 2009 ; Ponomarenko et al. , 2020 ) or with regulatory potential on traits related to menarche (e.g. ESR1 and FSHB ) ( Ponomarenko et al. , 2019 , 2020 ), which are risk factors for developing uterine fibroids ( Styer and Rueda, 2016 ; Wise and Laughlin-Tommaso, 2016 ; Pavone et al. , 2018 ). Age at menarche also has a high proportion attributable to genetic effects (45%) ( Snieder et al. , 1998 ). We anticipated that the direction of minor allele associations with age at menarche and uterine fibroids as well as the relative direct contribution of genetic variants to the pathophysiology of uterine fibroids and the relative contribution mediated by age at menarche might differ. To the best of our knowledge, this is the first study to investigate whether and to what extent age at menarche mediates the association between genetic variants and uterine fibroids ( Fig. 1 ). The detailed phases of our genome-wide causal mediation analysis are provided in Fig. 2 . This study identified new genetic associations, providing insights into the proportion of the risk for uterine fibroids that is attributable to genetic effects mediated by age at menarche. Causal diagram indicating possible relationships of SNPs with age at menarche and uterine fibroids that we assumed and analyzed in our mediation analysis.   represents the effects of SNPs unmediated by age at menarche on uterine fibroids, classified as direct effects. represents the effects of SNPs mediated by age at menarche on uterine fibroids, classified as indirect effects. represents the effects of factors that influence both age at menarche and uterine fibroids, classified as confounding effects. SNPs, single-nucleotide polymorphisms. Schematic diagram indicating processes of data acquisition from two cohorts, data quality control, GWAS, causal mediation analysis, data validation, number of identified SNPs and three main outputs. GWAS, genome-wide association studies; SNPs, single-nucleotide polymorphisms.

Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file.