Authors: Sven E. Ojavee, Liza Darrous, Marion Patxot, Kristi Läll, Krista Fischer, Reedik Mägi, Zoltan Kutalik, Matthew R. Robinson
Categories: Article, age at natural menopause, age at onset, age-specific genetic effect, aging, common complex disease, disease etiology, genome-wide association study, interaction, significance testing, time to event
Source: American Journal of Human Genetics
There is currently little evidence that the genetic basis of human phenotype varies significantly across the lifespan. However, time-to-event phenotypes are understudied and can be thought of as reflecting an underlying hazard, which is unlikely to be constant through life when values take a broad range. Here, we find that 74% of 245 genome-wide significant genetic associations with age at natural menopause (ANM) in the UK Biobank show a form of age-specific effect. Nineteen of these replicated discoveries are identified only by our modeling framework, which determines the time dependency of DNA-variant age-at-onset associations without a significant multiple-testing burden. Across the range of early to late menopause, we find evidence for significantly different underlying biological pathways, changes in the signs of genetic correlations of ANM to health indicators and outcomes, and differences in inferred causal relationships. We find that DNA damage response processes only act to shape ovarian reserve and depletion for women of early ANM. Genetically mediated delays in ANM were associated with increased relative risk of breast cancer and leiomyoma at all ages and with high cholesterol and heart failure for late-ANM women. These findings suggest that a better understanding of the age dependency of genetic risk factor relationships among health indicators and outcomes is achievable through appropriate statistical modeling of large-scale biobank data.
Keywords: genome-wide association study, age at natural menopause, age-specific genetic effect, aging, age at onset, time to event, common complex disease, disease etiology, significance testing, interaction
Age at onset and time-to-event observations are among the most important traits of interest in cohort studies of age-related diseases as they are critical to gain insight into the genetics of disease development and progression.^1^^,^^2^ The underlying etiology of age-related outcomes likely reflects a variety of biological processes that are triggered at different stages of life, long before the onset of observable symptoms. As a result, the underlying genetic propensity for outcomes may vary with age and depend upon different sets of genetic risk factors at different time points, reflecting the range of underlying molecular mechanisms that shape the onset distribution. Therefore, identifying the genetic variants associated with onset at different stages of life will improve our understanding of disease progression.
Here, we seek to test the hypothesis that genetic propensity for age at onset is age specific by focusing on the most commonly experienced timing-related phenotype in the human population, age at natural menopause (ANM). Menopause is the permanent cessation of the menstrual cycle in women following the loss of ovarian function and occurs at an average age of 51 years, with 4% of the female population experiencing early menopause prior to age 45. Current evidence suggests that early menopause is associated with a risk for cardiovascular disease^3^ and osteoporosis,^4^ and late menopause is associated with a risk for breast cancer.^5^ Recent genomic studies find ∼50% of menopausal timing variation is attributable to genetic markers^6^ that are linked to regulation of DNA repair and immune function.^7^^,^^8^^,^^9^ However, previous analyses make strong assumptions that genetic effects are constant throughout life (Figure 1A). By modeling the quantitative genetic basis of ANM in a way that enables detection of the age at which genetic risk factors have the greatest influence, we report evidence for widespread age-specific genetic effects underlying population-level variation in ovarian aging in both the UK and Estonian Biobank data.
Figure 1 Statistical model description and previously unreported discoveries(A) The CAMP model enables a more flexible and accurate description of the SNP effect size by introducing a slope term. The linear change model enables three example questions to be (1) what is the interval at which there is a significant effect on the trait (“significant interval”)? (2) At which age is there strongest evidence for an effect (“strongest evidence for effect”)? (3) Is the slope (β1) significantly different from zero (“β1”)? Even though a more complex model can result in generally wider confidence intervals, it can still result in a more accurate representation of the effect size, often accompanied by higher statistical power. By estimating the effect-size change, it is also possible to accurately determine trends, which constant effect assumptions cannot capture. The lower two panels represent corresponding survival function estimates showing the probability that an individual has not experienced an event. We observe that under linear effect assumption, the survival functions can cross, which is not possible with the constant effect assumption.(B) 19 previously unreported discoveries for ANM from the CAMP model across the UK and Estonian Biobank data. The coral line indicates the genome-wide significance level of 5×10−8. Arrows indicate the positions of the previously unreported discoveries.
We present a marginal Cox age-specific mixed proportional hazards (CAMP) model paired with a unique significance-testing framework. Following the success of many genome-wide association studies (GWASs), more attention has been attributed to better characterizing SNP effects under different environmental conditions, leading to genotype-covariate analyses.^10^ For a continuous trait, one of the simplest ways to model this type of interaction would be to include a linear interaction term. For example, to estimate the impact of age on SNP effects, we could write the model as
where (for individual i and SNP j) yi is a continuous trait, μ is the intercept, xij is the SNP value, βj0 is the SNP effect at time t0, ti is the age when yij is measured, βj1 is the linear effect of age, and εi is the residual variance. However, for age-at-onset phenotypes, the model as specified in Equation 1 would not be identifiable because yi=ti. Previous studies have proposed to analyze age-specific effects by splitting timescales into non-overlapping intervals. Individuals who have the event in a future interval are treated as right censored, with individuals who have had an event in a previous interval excluded from the analysis. For example, this idea has been suggested by Joshi et al.,^11^ where time intervals of 40–75 and 75+ were used. Although it is correct to conduct the analysis in such a way, doing so requires defining intervals that could be seen as an arbitrary choice, with ill-defined intervals leading to an incomplete understanding of the effect-size distribution. Furthermore, this type of modeling will not scale well with the added number of intervals as each interval requires an additional parameter. Therefore, we propose a form of Cox PH model that allows specifying a functional shape for age-specific effects. To estimate the marginal effect of SNP j in chromosome k, the general form for CAMP model for each SNP j∈1,…,M is
where λ0(t) is the baseline hazard, i denotes the ith individual, xij is the standardized jth marker value, gik is the genetic predictor from all other chromosomes other than the SNP j is located at, ξj is its corresponding effect when estimating marker j, zi is the summarized covariate value, and δj is its corresponding effect when estimating marker j. βj(t) is the effect-size change function for SNP j and we define it in two ways. The effect-size function βj(t) should be defined such that its domain is the set of positive real numbers and it is (piecewise) differentiable. Here, we define the effect-size function as a linear function in age
where βj0 is the intercept term, effectively estimating the effect size at time t0, and βj1 is the slope showing how much the effect changes each year. It is possible to define any other parametric shape; for example, the exponential decay function could be a natural choice. However, in the first experiments on real data, the linear effect-size change gave a higher likelihood compared with exponential decay models. Hence, we decided to resort to the linear effect-size model that is easier to interpret, and the variance function can be represented without Taylor expansion-based approximations. Given the effect-size function definition, we can calculate the variance at each time point t as
where Var(βj0), Var(βj1), and Cov(βj0,βj1) can be estimated from the Hessian of the Cox model. We use the linear effect-change definition from Equation 3 to test whether there exists a change in the effect size across the lifespan. As including the genetic values from other chromosomes shares properties with mixed modeling, we are going to refer to this model as the CAMP model.
At each time point t, we can define the test statistic function χj2(t) for SNP j as the square of the ratio of effect size and the standard error of the effect size
All that remains is to estimate gik, the genetic predictor from all other chromosomes other than the SNP j. To do this, we use a BayesW model,^1^ which assumes that for an individual i, the age at onset of a disease yi has Weibull distribution, with a re-parameterization of the model to represent the mean and the variance of the logarithm of the phenotype as
where Xφ is a standardized genotype matrix containing SNPs allocated to group φ, μ is an intercept, βφ is the vector of SNP effects in group φ, Zi is additional covariates (such as sex or genetic principal components [PCs]), δ is the additional covariate effect estimates, and α is the Weibull shape parameter. For each group, we assume that βφ are distributed according to a mixture of Gaussian components with mixture-specific proportions πφ and mixture variances σφ12,…,σφL2 and a Dirac delta at zero, which induces
where L is the number of mixture components. We estimated the hyperparameters such as genetic variance and prior inclusion probability by grouping markers into “MAF-LD” bins, as recent theory suggests this yields improved estimation.^12^ We used 20 MAF-LD groups that were defined as minor-allele frequency (MAF) quintiles, where we then split each quintile into quartiles by linkage disequilibrium (LD) score. The cutoff points for creating the MAF quintiles were 0.006, 0.013, 0.039, and 0.172; the cutoff points for creating LD score quartiles were 2.11, 3.08, and 4.51 for the first; 3.20, 4.71, and 6.84 for the second; 4.70, 6.89, and 9.94 for the third; 7.65, 11.01, and 15.70 for the fourth; and 10.75, 15.10, and 21.14 for the fifth MAF quintile, exactly as in the age-at-menopause analysis by Ojavee et al.^1^ The posterior mean BayesW model estimates of βφ are then used to create gik, the genetic predictor from all other chromosomes other than the SNP j. This gives a two-step leave-one-chromosome-out (LOCO) approach, where first a BayesW model is used to estimate the genetic predictor and then a marginal age-specific CAMP model is used for the second step. Next, we discuss how we can conduct significance testing in the second step in an efficient manner while ensuring that the type I error is bounded below the fixed threshold α even with a more complex model.
We demonstrate that to limit type I error rate below α and given the null hypothesis of no effect at SNP j (βj0=0, βj1=0), it is sufficient to compare the test statistic χj2(t∗) with the χdf=22 1−α quantiles at any time point t∗.
We will naturally assume that Varˆ(βˆj0)>0 and Varˆ(βˆj1)>0. As χj2(t) is twice differentiable, it is possible to find its respective first and second derivative. This will give us two extreme points t∗ at which the ddtχj2(t∗)=0. Then, the χ2 score has a local maximum if d2dt2χj2(t∗)<0 or the χ2 score has a local minimum weakest if d2dt2χj2(t∗)>0. It can be shown that the function χj2(t) has two extreme points located at
Using the second derivative, it can be shown that χj2(t) will always have a (global) maximum at t2∗ and a (global) minimum (χ2(t1∗)=0) at t1∗. We find that in the limiting cases, the test statistic is testing the significance of the slope limt→∞χj2(t)=(βj1)2Var(βj1) and limt→−∞χj2(t)=(βj1)2Var(βj1). As the domain of χj2(t) is the set of positive real numbers and there are no breakpoints in the function, then χj2(t) is bounded within interval [χj2(t1∗),χj2(t2∗)].
We are especially interested in the distribution of the maximum possible χ2 statistic χj2(t2∗) under the null hypothesis that both βj0=0 and βj1=0.
Under the null hypothesis that both βj0=0 and βj1=0, the chi-squared statistic evaluated at the maximum point t2∗ χj2(t2∗) follows a χdf=22
We define r:=Covˆ(βˆj0,βˆj1)Varˆ(βˆj0)Varˆ(βˆj1) and we express χj2(t2∗) such that it would be a sum of two uncorrelated random variables.
Under the null hypothesis of βj0=0 and βj1=0, we know that βˆj0Varˆ(βˆj0)∼N(0,1) and βˆj1Varˆ(βˆj1)∼N(0,1) and therefore
The last result implies that under the null of βj0=0 and βj1=0, Equation 12 is a sum of two uncorrelated standard Gaussian random variables squared, which means that χj2(t2∗) is from the chi-squared distribution with degrees of freedom of 2.
This naturally gives us a rule for hypothesis testing at time t2∗. If the test fails to disprove the null hypothesis at time t2∗, it will fail to disprove the null hypothesis at any possible t. If the test accepts the alternative hypothesis, it means that there must exist an interval (or at least one point) at which the variable has an effect on the phenotype.
Furthermore, we can show that the quantiles of χdf=22 distribution result in a stringent-enough test at any time point.
Suppose that we have estimated effect sizes βj0 and βj1 from a linear effect change model βj(t)=βj0+βj1(t−t0) and that the null hypothesis of βj0=0 and βj1=0 holds. Then, for every time point t, the probability of type I error (α) is bounded when using the χ22 distribution 1−α quantile as a critical value.
To prove the lemma, we need to demonstrate that 1−α quantile of χ22 distribution (q1−α) is greater than 1−α quantile of χ2(t) at any time point t under the null hypothesis βj0=0 and βj1=0. Suppose that the maximum test statistic value is achieved at t2∗ with value χj2(t2∗).
We suppose, in contradiction, that under the null hypothesis, there exists some time point t∗≠t2∗ at which the distribution of χj2(t∗) would have a higher 1−α quantile value q˜1−α than the 1−α quantile of χdf=22 distribution q1−α:
Given this, we can write the following inequalities,
where the first inequality follows from the contradiction and the second inequality from the fact that χj2(t∗) is the maximum possible value of the test statistic. The inequalities result in a contradiction, which therefore proves the lemma.
An important corollary of this result is that we can use χdf=22 quantiles to do statistical testing at any time point, and doing tests at (many) different time points will not increase type I error. For example, we can simultaneously test the significance of the slope (corresponds to limt→∞χj2(t)) and significance at t2∗ (using χj2(t2∗)) at the 1−α-quantile of χdf=22 distribution while limiting the type I error at α as given the effect and variance estimates.
This project uses UK Biobank data under project 35520. UK Biobank genotypic and phenotypic data are available through a formal request at http://www.ukbiobank.ac.uk. The UK Biobank has ethics approval from the North West Multi-centre Research Ethics Committee (MREC). We first restricted our analysis to a sample of European-ancestry UK Biobank (UKB) individuals. To infer ancestry, we used both self-reported ethnic background (UKB field 21000-0) and genetic ethnicity (UKB field 22006-0) and selected coding 1 in both cases. We projected the 488,377 genotyped participants onto the first two genotypic PCs calculated from 2,504 individuals of the 1,000 Genomes project. Using the obtained PC loadings, we then assigned each participant to the closest 1,000 Genomes project population, selecting individuals with PC1 projection < absolute value 4 and PC2 projection < absolute value 3. Samples were also excluded based on UKB quality control procedures with individuals removed of (1) extreme heterozygosity and missing genotype outliers; (2) a genetically inferred gender that did not match the self-reported gender; (3) putative sex chromosome aneuploidy; (4) exclusion from kinship inference; and (5) withdrawn consent. We used genotype probabilities from version 3 of the imputed autosomal genotype data provided by the UKB to hard call the genotypes for variants with an imputation quality score above 0.3. The hard-call threshold was 0.1, setting the genotypes with probability ≤0.9 as missing. From the good-quality markers (with missingness less than 5% and p value for Hardy-Weinberg test larger than 10−6, as determined in the set of unrelated Europeans), we selected those with MAF > 0.0002 and rs identifier in the set of European-ancestry participants. We then took the overlap with the Estonian Biobank data described below to give a final set of 8.7 million SNPs using both autosomal chromosomes and the X chromosome. This provides a set of high-quality SNP markers present across both discovery and prediction datasets.
We created the phenotypic data of ANM similarly to Ojavee et al.^1^ We used UKB field 3,581 to obtain the time, if available, and excluded from the analysis (1) women who had reported having and later not having had menopause or vice versa, (2) women who said they had menopause but with no record of the time of menopause (UKB field 2,724), (3) women who have had a hysterectomy or the information about this is missing (UKB field 3,591), and (4) women whose menopause is before age 33 or after 65. Within the UKB data, there were a total of 173,424 unrelated (only one person kept from second-degree or closer relative pairs) European-ancestry women, out of which 125,697 had experienced menopause and 47,727 had not had menopause based on data field 2,724. For computational convenience when conducting the joint BayesW analysis, we created an additional subset of markers by removing markers in very high LD through the selection of the highest MAF marker from any set of markers with LD R2≥0.8 within a 1-Mb window. These filters resulted in a dataset with 173,424 individuals and 2,174,071 markers for the first-step estimation of the LOCO genetic predictors, and then in the second step CAMP model we analyzed 8.7 million SNPs using both autosomal chromosomes and the X chromosome.
To replicate the findings, we used the Estonian Biobank with 70,082 women (22,740 with menopause and 47,342 without menopause). For access to be granted to the Estonian Biobank genotypic and corresponding phenotypic data, a preliminary application must be presented to the oversight committee, who must first approve the project. Ethics permission must then be obtained from the Estonian Committee on Bioethics and Human Research, and finally a full project must be submitted and approved by the Estonian Biobank. This project was granted ethics approval by the Estonian Committee on Bioethics and Human Research (https://genomics.ut.ee/en/content/estonian-biobank). Similar to the UKB data, we used the age range of 33–65 and thus excluded women with age of menopause outside this interval. In the total Estonian Biobank data, there were 195,432 individuals genotyped on the Illumina Global Screening Array (GSA), which were imputed to an Estonian reference created from the whole-genome sequence data of 2,244 participants.^13^ From 11,130,313 markers with imputation quality scores >0.3, we selected SNPs that overlapped with those selected in the UKB as described above by using the same SNP sets for the first and second steps of the analyses.
We investigated the three questions as proposed in Figure 1: (1) we checked the interval at which there is a significant effect, (2) we tested the significance at ages at which the variants had the most evidence for an effect, and (3) we tested for the existence of an age-specific effect.
We evaluated the test statistic function using Equation 5 at ages 41, 43, 45, 47, 49, 51, 53, and 55 and calculated the p values using the χdf=22 distribution quantiles, as suggested by our theory above. On these results, at each age, we applied the clumping procedure (using plink 1.9^14^) with a window size of 1 Mb, LD threshold of r2=0.05, p value threshold for index SNPs of p=5×10−8, and no p value threshold for other SNPs belonging to a clump of an index SNP. To detect clumps with an independent signal, we applied the COJO procedure^15^ implemented in GCTA software^16^ with a window size of 1Mb, and the SNPs were considered independent if the p value in the joint model was less than 5×10−8. The independent index SNPs from the COJO analysis were then replicated at each age in the Estonian Biobank data with replication defined as a p value lower than 0.05 and the same effect-size estimate sign as in the discovery analysis. To check the period during which a SNP has a significant effect (Figure 2A), we checked whether the same SNP or a SNP in the same clump also has an effect in the consecutive grid point. Specifically, we took all the significant independent and replicated SNPs at ages 41 and 43 and then checked whether the index SNPs of age 41 mapped directly to an index SNP at age 43 or a clump of an index SNP. Then, we compared the ages 43 and 45 and iteratively so forth until 55.
Figure 2 Age distribution of significant effects and effect-size change(A) We evaluated the effect-size and standard error estimates for every SNP at each age on a grid from 41 to 55 and counted the number of significant hits replicated in the Estonian Biobank. The significant focal hits were mapped to consecutive ages, summarizing the count since when the effects were significant.(B) Effect sizes (log hazard ratio, HR) for 245 significant effects; the majority of variants have a larger absolute effect size at age 41 than at age 55, and only 9 variants have an increasing effect size. We observe that the model manages to capture the effect-size change for many variants.(C) Classification of the menopause-associated variants by age-specific evidence by testing whether the slope parameter is equal to zero. Variants with weak evidence have a p value lower than p<0.05, moderate evidence requires p<5×10−8, and strong evidence requires the slope to be significant also in the replication dataset.
We tested the significance at ages at which the variants had the most evidence for an effect to understand the total number of significant effects and possibly identify previously undiscovered loci. Furthermore, we defined a period of interest between ages 45 and 52 so the significance would only be evaluated during this period. That was done to avoid over-interpretation of the linear effect at uncommon high or low ages. Therefore, using Equation 10, we first calculated the age t2∗ at which the test statistic achieves the highest value, and secondly, we evaluated the chi-squared function (Equation 5) at that age t2∗. If the time point t2∗ was outside of the interval (45,52), then we instead evaluated the function at the time point where the function χj2(t) achieved its maximum within the interval (either age 45 or 52). We compared the chi-squared statistics with the χdf=22 distribution to calculate the p values. We applied a similar procedure to the previous case. First, we applied the clumping procedure (using plink 1.9^14^) with a window size of 1 Mb, LD threshold of r2=0.05, p value threshold for index SNPs of p<5×10−8, and no p value threshold for other SNPs belonging to a clump of an index SNP. Secondly, to detect clumps with an independent signal, we applied the COJO procedure^15^ implemented in GCTA software^16^ with a window size of 1 Mb, and the SNPs were considered independent if the p value in the joint model was less than 5×10−8. We then checked the independent and significant SNPs from the COJO analysis for previous association signals. We removed all the markers that had a correlation of r2>0.1 with a marker that had been previously found associated with ANM using the GWAS catalog (published until April 2022) and LDtrait tool with the British in England and Scotland population. Furthermore, we specifically compared our candidate set for the significant SNPs reported by Ruth et al.,^9^ removing the markers reported by them and those with a correlation r2>0.1. Finally, we checked our candidate set with the Phenoscanner database^17^^,^^18^ to find any previous associations with variants of interest or variants in LD. The index SNPs not removed by the three filters were then replicated at age with maximum evidence in the period (45,52) in the Estonian Biobank data, with replication defined as a p value lower than 0.05 and the same effect estimate sign as in the discovery analysis.
To verify whether there is a genome-wide significant age-specific effect for every SNP j, we checked whether the slope parameter is significantly different from 0
where βj1 is estimated from the model specified in Equation 2. The chi-squared statistic was calculated as the squared ratio of the slope size estimate and standard error estimate. As this quantity naturally corresponds to limt→∞χj2(t) where χj2(t) is defined as in Equation 5, we compare the test statistic again with the χdf=22 distribution quantiles to get the p values. Similarly to two previous cases, we first applied the clumping procedure (using plink 1.9^14^) with a window size of 1 Mb, LD threshold of r2=0.05, p value threshold for index SNPs of p<5×10−8, and no p value threshold for other SNPs belonging to a clump of an index SNP. Second, to detect clumps with an independent signal, we applied the COJO procedure^15^ implemented in GCTA software^16^ with a window size of 1 Mb, and the SNPs were considered independent if the p value in the joint model was less than 5×10−8. The independent index SNPs from the COJO analysis were then replicated in the Estonian Biobank data, with replication defined as a p value lower than 0.05 and the sign of the slope the same as in the discovery analysis. We classified the SNPs with different levels of evidence of age-specific effects. SNPs with a p value below the nominal significance threshold (p<0.05) are said to have at least weak evidence for an age-specific effect; variants with a p value below the genome-wide significance threshold (p<5×10−8) are said to have at least moderate evidence for an age-specific effect; variants with moderate evidence that also replicate in the Estonian Biobank are considered to have strong evidence for an age-specific effect. The variants that do not fall under these three categories are said to have no evidence for age-specific effects.
We used recently presented Downstreamer software^19^ to identify genes connected to our association study results through gene expression and to identify enriched pathways. We calculated test statistics using Equation 5 at ages 41, 43, 45, 47, 49, 51, 53, and 55 and calculated the p values using the χdf=22 distribution quantiles, as suggested by our theory above.
Downstreamer implements a strategy that accounts for LD structure and chromosomal organization, operating in two steps. In the first step, gene-level prioritization scores were calculated for each age group’s summary statistics and a null distribution. This aggregates p values per variant into a p value per gene while accounting for local LD structure. GWAS gene p values were calculated for all 20,327 protein-coding genes (Ensembl release version 75). The gene p values were then converted to gene Z scores for use in subsequent analysis. To account for the long-range effects of haplotype structure, which results in genes getting similar gene Z scores, a generalized least-squares (GLS) regression model is used, which requires a gene-gene correlation matrix. This correlation matrix is calculated by first simulating 10,000 random phenotypes by drawing phenotypes from a normal distribution and then associating them to the genotypes of the 1000 Genomes phase 3 non-Finnish European samples. The GWAS gene Z scores for each of the 10,000 simulated GWAS signals alongside the Pearson correlations between the GWAS gene Z scores are then calculated. Correlations between simulated GWAS gene Z scores reflect the underlying LD patterns and chromosomal organization of genes. An additional 10,000 GWASs were simulated to empirically determine enrichment p values, and, finally, an additional 100 simulations were used to estimate the false discovery rate (FDR).
In the second step, the gene-level prioritization scores were associated with the co-regulation matrix and pathway annotations. We used a previously generated co-regulation matrix that is based on a large multi-tissue gene network.^19^ Publicly available RNA-seq samples were downloaded from the European Nucleotide Archive (https://www.ebi.ac.uk/ena) containing 56,435 genes and 31,499 samples covering a wide range of human cell types and tissues. 165 leading PCs representing 50% of the variation were selected. For protein-coding genes, centered and scaled eigenvectors for these 165 components (mean = 0, SD = 1) were calculated. The co-regulation matrix is then defined as the Pearson correlation between the genes from the scaled eigenvector matrix, with diagonal zero and Pearson r values converted to Z scores. To identify pathway and disease enrichments, the Human Phenotype Ontology (HPO), KEGG, Reactome and Gene Ontology (GO) Biological Process, Cellular Component, and Molecular Function databases were used. A Z score per pathway or term per gene is calculated, giving how much each gene contributes to these gene sets. We collapsed correlated genes in parallel with the GWAS step to ensure compatibility with the GWAS gene Z score and scaled all pathway Z scores to zero mean and unit variance.
Real and simulated GWAS Z scores were rank-based inverse-normal transformed. A linear model was used to correct for gene length, as longer genes will typically harbor more SNPs. Genes with a Pearson correlation r≥0.8 in the 10,000 GWAS permutations were treated as one gene. A GLS regression is used to associate the GWAS gene Z scores to the pathway Z scores and co-regulation Z scores, with β=(XTω−1X)−1XTω−1y, where β is the estimated effect size of a pathway, term, or gene from the co-regulation matrix; ω is the gene-gene correlation matrix; X is the design matrix of real GWAS Z scores; and y is the vector of gene Z scores per pathway, term, or gene from the co-regulation matrix. β for the 10,000 random GWASs were estimated in the same way and used to estimate the empirical p value for β. These two analyses result in pathway enrichments and core gene prioritizations, respectively. The gene-gene correlation matrix derived from the 10,000 permutations is used as a measure of conditional covariance of the error term (ω) in the GLS to account for the relationships between genes due to LD and proximity. We present full results for each age in Table S2. Pathway enrichment analysis for age 41, Table S3. Pathway enrichment analysis for age 43, Table S4. Pathway enrichment analysis for age 45, Table S5. Pathway enrichment analysis for age 47, Table S6. Pathway enrichment analysis for age 49, Table S7. Pathway enrichment analysis for age 51, Table S8. Pathway enrichment analysis for age 53, Table S9. Pathway enrichment analysis for age 55. In the main figures, we present the results significant at both FDR correction and Bonferroni correction offered by the Downstreamer software.
We used LD score regression to calculate genetic correlations among the test statistics generated using Equation 5 at ages 41, 43, 45, 47, 49, 51, 53, 55, and among these ages and other phenotypes using publicly available GWAS summary data. We used UKB results for 100 phenotypes released by Neale group and the Global Biobank Meta-analysis Initiative consortium. The significance threshold of 0.05 was corrected by the total number of tests (800). We present only estimates in Figure 4 where the correlation was significant in at least one age point.
Figure 4 Age-specific genetic correlations and causality of ANM and health-related outcomes(A and B) We evaluated the effect-size and standard error estimates for every SNP at each age on a grid from 41 to 55, and from the resulting summary statistics we tested for genetic correlations among (A) age groups and (B) among observed ANM values at the age groups and 100 other health-related indicators and outcomes using LD score regression. In (B), we present correlations for outcomes with evidence of a significant non-zero genetic correlation at one age group or more.(C) In (C), we present results from inverse-weighted Mendelian randomization that estimates the potential causal relationship between ANM and outcomes where a significant genetic correlation was found in (B) across ages. Black boxes depict ages for which significant causal estimates were found. Major depression is excluded from (C) as there are no instrumental variables available for this analysis. Full results from a range of Mendelian randomization models are given in Table S10.
We calculated the causal effect estimates that ANM at different ages has on various traits (Table S11) using Mendelian randomization (MR), a statistical method that utilizes the randomized inheritance of genetic variations in the population to estimate the potential causal effect a modifiable risk factor or exposure has on a health-related outcome of interest.^20^^,^^21^ The genetic variants used as instrumental variables (IVs) for our exposure were selected to have a genome-wide significant association with the exposure (p<5×10−8) and were then pruned using LD distance to ensure that they were independent. This was done using the “ld_clump” function of the “ieugwasr” R package^22^ with default settings (clump_kb = 10,000, clump_r2 = 0.001, clump_p = 0.99, pop = “EUR”). After the IVs of our exposure were selected, their association effects were then obtained for each of our outcome traits of interest. A single-sided t test was carried out to check whether the IVs had a stronger association with the outcome than with the exposure and were subsequently removed if so (for violating the MR assumptions). The two sets of association effects were then harmonized and used to calculate the causal effect estimates using the inverse-variance weighted method found in the “TwoSampleMR” R package.^23^ This analysis was repeated for each varying age of our exposure. It is important to note that in the case of educational attainment as an outcome, there were few exposure IVs that overlapped with the outcome genetic variants, especially as the age increased; hence, in the presence of a single IV, a Wald ratio was used to calculate the MR causal effect estimate. Moreover, when the trait was of a case-control nature, the effective sample size was calculated using the following (4^∗^cases^∗^controls)/(cases+controls).^24^
We show in the methods that CAMP determines the time dependency of marker age-at-onset associations without a significant multiple-testing burden. We applied the CAMP model to 173,424 unrelated observations of self-reported ANM in the UKB data (125,697 reported events and 47,727 right-censored observations; 8,747,951 SNPs) and 70,082 observations in the Estonian Biobank (22,740 reported events and 47,342 censored observations, Figure S1). We find 312 ANM associations in the UKB, of which 226 replicate previous studies^7^^,^^8^^,^^9^ and 19 are previously unreported and replicate for the first time within the Estonian Biobank (Table 1; Figures 1B and S2). In addition, we find 67 associations that have not previously been reported, but they did not replicate in the Estonian Biobank. Nevertheless, 46 out of 67 previously unreported associations show consistency with signs in the discovery and replication datasets (Fisher’s exact test, p=0.007), suggesting that a larger replication dataset could lead to further replications. To test for novelty, we test for significance at time points where the evidence is the highest within the intervals 45 to 52 (Figure 1A, see methods). For significant SNPs detected in the UKB, the age distribution of maximum association evidence is concentrated between the ages of 43 and 51 (Figure S3). That is different compared with the maximum association evidence age distribution for all SNPs, which has thicker tails with nearly four times higher standard deviation even if the distributions have similar centers (median age of 51 and 49 for all SNPs and significant SNPs, respectively) (Figure S3). In conclusion, the CAMP approach yields an increase in previously unreported 8% of the replicated marker associations are previously unreported, a somewhat expected increase resulting from our approach leveraging potentially existing age specificity.
For quantifying the existence of age-specific effects, we first test the null hypothesis of whether the slope term βj1 is different from 0, with a rejection of the null implying the existence of a time-varying genetic effect. Of the 245 UKB associations (226 previously reported + 19 previously unreported associations), we find that 72% (176) show at least nominally significant (p<0.05) age-specific effects within UK women (Figure 2C). 37 of these 176 associations, representing 15% of all associations, have a slope with a genome-wide significant p value (p<5×10−8), constituting a more stringent criterion (strong and moderate evidence in Figure 2C). Checking the significance of the slope for all SNPs (not only genome-wide significantly associated to ANM), we find 63 regions that exhibit a significant slope term in the UKB (Tables 2 and S1), and we replicate the age-specific effects for 20 regions in the Estonian Biobank (Table 2), yielding a replication rate of 32% (Figure 2C). These 20 variants have stronger effect sizes earlier in life that mostly decay toward zero after age 50 (Figure S5), making them early-ANM specific. Although 139 variants do not pass the threshold for genome-wide significance, they still indicate that for many regions previously identified as menopause associated, the assumption of constant effect size (assumption of proportional hazards at the SNP) is generally invalid. Indeed, the 43 UKB-discovered variants with significant slope terms that did not replicate in the Estonian Biobank had effect-size directions that were broadly concordant across studies (Fisher’s exact test p=0.051, Table S1).
Second, we observe that the number of regions affecting ANM changes considerably with the peak number of ANM-affecting regions observed at age 49 (Figure 2A). Moreover, we find that the period during which a particular region can significantly impact ANM varies considerably with only half of the significant associations at age 47 also significant at age 41. In general, we observe that effects tend to become insignificant with increasing age, with the drop in significance occurring at age 53, so that by age 55 only 8 loci have a genome-wide significant effect on ANM (Figure 2A). A similar result can be seen if we observe the distribution of ages when the evidence for the menopause effect is the strongest (Figure S3), as very few significant SNPs achieve the strongest association after age 51.
In contrast with most associations discovered at age 49, the general trend across 245 significant SNPs is that the effect size estimates shrink toward zero (Figure 2B). That might imply that the increase in the number of discoveries in the period 41 to 49 is instead due to the reduction in the standard error, and with a higher sample size, it could be possible to detect more associations already at earlier ages. Interestingly, only 9 of the significant SNPs have a larger absolute effect size at age 55 than at 41. That is in line with many previous results reporting a reduction in relative genetic risks with the increase in age.^25^^,^^26^ Finally, we observe that there exists a stark difference between the effect-size profiles of significant and insignificant effects (Figure S4) with a much narrower effect-size distribution for the non-significant SNPs. Meanwhile, menopause-associated variants stand out as their effect size can change greatly across the period of interest. We find it important to stress that the interpretation of the effects is done using (log−) hazard ratios as the underlying model is the CAMP model. Specifically, in the context of the CAMP model, hazard ratios need to be interpreted at each age separately as the hazard ratio changes with age, for example as seen in Figure 2.
Our analysis differs in one other key way from previous ANM genetic association studies. Here, we do not censor women who were placed on hormone replacement therapy (HRT). In survival models, declaring HRT individuals as censored makes the modeling assumption that age-at-HRT start and ANM are independent. They are clearly not, as for women on HRT there is a correlation between the age of HRT and ANM of 0.58 within the UKB. For women who were placed on HRT prior to the recorded date of ANM, this correlation is stronger at 0.69. A Cox Proportional Hazards model for ANM including a categorical covariate of whether a woman was given HRT prior to menopause (1 if on HRT prior to menopause, 0 otherwise) shows that censoring for HRT prior to ANM would significantly censor for earlier menopause (HR = 0.95, p=6.08×10−15). Thus, censoring for HRT is not the optimal modeling choice, and additionally, it results in the loss of 34,031 observations, reducing power. Nevertheless, we conduct a sensitivity analysis of our estimated effect sizes with and without adding HRT as a time-varying covariate to the CAMP model at the 312 top loci (both replicated and unreplicated regions) identified within our study (Figure S6). We find very strong concordance of effect sizes across loci (Figure S6), highlighting that in practice these different modeling choices have no detectable impact on the leading SNP association findings.
Finally, we highlight notable examples of the previously unreported replicated associations with significant slope terms. For example, we note that 63,595,648 (in GRCh37 coordinates), which is downstream of SPINDOC, is increasingly associated with menopause genetic risk as age increases, with the highest effect size being at later age groups. Another example is 63,795,628, which is upstream of USP3, where the menopause association disappears with increasing age. Both of these associations were previously suggestively associated with ANM, but they pass the significance threshold in the UKB, and they replicate in the Estonian Biobank, when our proposed model was used.
We conduct a number of follow-up analyses to support our age-specific association results. First, we test for significant enrichment of the summary statistics generated by our approach for each age group. For all categories showing significant enrichment after Bonferroni multiple testing correction, we find that their significance does not hold across all age groups (Figure 3; Table S2. Pathway enrichment analysis for age 41, Table S3. Pathway enrichment analysis for age 43, Table S4. Pathway enrichment analysis for age 45, Table S5. Pathway enrichment analysis for age 47, Table S6. Pathway enrichment analysis for age 49, Table S7. Pathway enrichment analysis for age 51, Table S8. Pathway enrichment analysis for age 53, Table S9. Pathway enrichment analysis for age 55).
Figure 3 Age-specific enrichment of genetic associations across multiple genomics resourcesWe evaluated the significance of every SNP at each age on a grid from 41 to 55 and from the resulting summary statistics we tested for enrichment across multiple genomics resources. Circle circumference gives the −log10 p value and the color gives the enrichment Z score calculated from the Downstreamer software. GTEx tissue-specific expression and GO terms are given on the y axis for annotations with genome-wide significance after multiple testing correction at one or more age groups. Full results are given in Table S2. Pathway enrichment analysis for age 41, Table S3. Pathway enrichment analysis for age 43, Table S4. Pathway enrichment analysis for age 45, Table S5. Pathway enrichment analysis for age 47, Table S6. Pathway enrichment analysis for age 49, Table S7. Pathway enrichment analysis for age 51, Table S8. Pathway enrichment analysis for age 53, Table S9. Pathway enrichment analysis for age 55.
Effect sizes for ANM between the ages of 41–49 were enriched in genes differentially expressed in the uterus, thyroid, prostate, ovary, fallopian tubes, and cervix within the Genotype-Tissue Expression (GTEx) consortium data (Figure 3). Additionally, we find enrichment between the ages of 41 and 49 for KEGG pathway NOTCH signaling associated with cell proliferation and death, the GO terms for an intrinsic pathway for apoptosis, and BH3-only proteins (Figure 3). In contrast, associations with variation in ANM for individuals older than 51 were all enriched for genes with differential expression in several brain regions within the GTEx data with no evidence for enrichment in reproductive tissues (Figure 3). These results suggest that genetic effects may differ across the age range.
Our next follow-up analysis used LD score regression, where we find that genetic correlations across ages are significantly less than 1 (Figure 4A). Genetic correlations of ANM and other phenotypes were also largely age dependent (Figure 4B). Note here that effect-size estimates for ANM are calculated on the menopause hazard scale, and thus a positive correlation estimated by LD score regression would refer to a high hazard of ANM (earlier ANM) corresponding to high trait values; in other words, the observed value of ANM and the trait are in fact negatively correlated. Thus, to ease interpretation, we flip the sign of the estimated correlation to display the genetic correlation of the observed values of ANM and each trait.
Between 41 and 49 years, we find a significant positive genetic correlation of observed ANM values with age at last birth and age at first birth (Figure 4B), implying a genetic relationship between later reproduction and later ANM within this age range.
Of significant note, genetic propensity for breast cancer was significantly associated with a later ANM (Figure 4B) before age 51. This supports previous evidence, where genetically mediated delays in ANM were found to increase the relative risks of several hormone-sensitive cancers.^8^ Additionally, evidence linking exposure to high levels of estrogen hormones with an increased risk of breast cancer is supported by a significant positive genetic correlation of ANM and oestradiol levels for women of early menopause, implying a genetic propensity for high estrogen levels associated with a genetic propensity for a later ANM (Figure 4B). Furthermore, we find that a high genetic risk for leiomyoma is consistently associated with later ANM (Figure 4B). Together with our enrichment results presented above (Figure 3) showing early ANM genetic associations are enriched for genes differentially expressed in female reproductive organs, oocytes, and DNA damage-repair mechanisms, our findings suggest that at the genetic level, breast cancer risk, hormone levels, and ANM are linked prior to age 50.
We also find significant genetic correlations implying that genetically mediated later ANM is correlated with a lower genetic predisposition to hypothyroidism, stroke, major depression, blood albumin levels, and obesity prior to age 51 (Figure 4B). For women at later ages, we find positive genetic correlations of ANM with cholesterol, low-density lipoprotein (LDL), obesity, and heart failure, implying that later ANM is correlated with an increased genetic predisposition for these metabolic-associated health measures (Figure 4B).
We find a significant positive genetic correlation of ANM values with both educational attainment and fluid intelligence between 41 and 49 years (Figure 4B), implying a genetic relationship between later reproduction, higher education, and later ANM for women before the age of 50. Interestingly, these genetic correlations also significantly change in sign for women whose ANM occurred after age 53, with a significant negative genetic correlation of ANM with educational attainment and fluid intelligence (Figure 4B), implying delayed reproduction and high educational attainment are associated with reproductive senescence post-age 50.
In a further follow-up analysis, we used MR, which utilizes the randomized inheritance of genetic variations in the population, to estimate the potential causal effect a modifiable risk factor or exposure has on a health-related outcome of interest. We used menopause at different ages as an exposure in five different MR methods (weighted median, inverse-variance weighted, simple mode, weighted mode, and MR-Egger; see Table S10) found in the “TwoSampleMR” R package. Note here again that effect-size estimates for ANM are calculated on the menopause hazard scale, and thus, to ease interpretation, we flip the sign of the estimated potential causal effect to give values on the observed ANM scale. When repeating the analysis for each varying age of our exposure, we find changes in the magnitude of the potential causal effect with age for educational attainment, leiomyoma, oestradiol, and neutrophil count (Figure 4C).
Taken together, we find that the majority of ANM genetic associations display some form of age specificity in their effects. In turn, that translates into the associations being differentially enriched in different biological pathways across ages, which then leads to different genetic associations of ANM and other health indicators and outcomes depending upon the timing of ANM, with different potential statistical causal relationships.
We find evidence that prolonged and delayed reproduction are genetically associated with reproductive senescence post-age 50 as the genetic correlations significantly change in the sign for women whose ANM occurred after age 53, with a significant negative genetic correlation of ANM with age at last birth and age at first birth (Figure 4B). Similarly, we find genetic correlations between ANM and educational attainment or fluid intelligence significantly turn negative for women after age 53. The latter patterns of changing genetic correlation may simply reflect changes in schooling opportunities or system. Furthermore, both age at first and last birth show a similar, but stronger, pattern, and it may be the underlying factor acting as a confounder in the education-ANM relationship.
Complementing the results from genetic correlations, our enrichment analysis results (Figure 3) show early ANM genetic associations being enriched for genes differentially expressed in female reproductive organs, oocytes, and DNA damage-repair mechanisms. Hence, our findings suggest that at the genetic level, breast cancer risk, hormone levels, and ANM are linked prior to age 50. The patterns observed in the MR analyses largely reflect those of the genetic correlations described above, but here we find little evidence for a causal relationship between ANM and breast cancer, heart failure, age at first or last birth, or hypothyroidism (Figure 4C). That implies that genetic correlation estimates likely reflect reverse causation or the presence of heritable confounders of the trait pairs.
Our enrichment analysis findings support a link between DNA damage-repair genes and repair and surveillance for the development of oocytes for early-ANM women. The size of the initial oocyte pool at birth, along with the rate of atresia, influences the age at which the oocyte pool is depleted. The meiosis that occurs in oocytes necessitates programmed double-stranded breaks (DSBs) that must be repaired through the homologous recombination pathway, with oocytes that do not properly repair DSBs after this first phase of meiosis undergoing apoptosis. Here, early-ANM-associated common variants are enriched at loci, harboring genes involved in the DNA repair and replication checkpoint processes, such as RNA polymerase II, histone methyltransferase complex, and histone acetylation (Figure 3). One-carbon metabolism has the ability to regulate the estrus cycle and modulate the initiation of reproductive senescence through the loss of methyl-donor production needed to properly maintain the epigenome. Our results support the existence of this mechanism as early-menopause associations are enriched in pathways associated with the hypothalamic-pituitary-gonadal (HPG) axis and with methylation in the nucleosome, with later-menopausal genetic associations showing no evidence of enrichment in these pathways (Figure 3). In humans, it has been suggested that postmenopausal women exhibit accelerated aging compared with premenopausal women of the same biological age.^27^ However, the cause-effect relationship between epigenetic changes and reproductive senescence remains unclear, and our results imply early-ANM women may have a methylation pattern associated with one-carbon metabolism that differs from the general population. Generally, our follow-up analyses support previous studies,^9^^,^^28^ but we demonstrate that almost all underlying pathways associated with variation in ANM act in an age-specific manner.
There are several important caveats to our study. First, we have assumed that the effect size can only change linearly with age, whereas in reality, they could consist of more complicated patterns that could be captured with piece-wise exponential models. However, introducing many more parameters on a genome-wide scale would lead to a high multiple-testing burden, potentially hampering the capability to detect the actual signal. Furthermore, especially for traits with a moderate range of values (90% of the observed ANM happens between ages 45 to 55, Figure S1), introducing many parameters could lead to overfitting. Therefore, we find that, especially in the context of traits such as ANM, assuming a linear effect change is a suitable compromise between the added value of learning new information about effect change and limiting the model complexity without damaging the statistical power. Nevertheless, the analyses presented here represent a first step, and we encourage specifying different functional forms for the effect size, preferably for traits with a broader range of values or a reasonable prior guess. Similarly, we have refrained from drawing conclusions about the causal effect change tendencies (e.g., linear vs. non-linear) as analyses claiming statistically significant non-linearity would be under-powered.
Second, the current implementation of the model is not computationally efficient, and to handle the computational burden we have utilized the computational resources of two universities to produce these results. Our objective was simply to highlight the existence of changing genetic relationships between phenotypes and health outcomes across the lifespan. Although it is possible to make marginal analyses embarrassingly parallel, it is inherently time consuming to fit a CAMP model. Scaling the inference requires new research into novel algorithms for computationally heavy high-dimensional statistical problems of this kind.
Finally, our study focused only on European ancestry individuals in the UK and Estonian Biobanks, and future analyses must take into account populations with more diverse ancestries to get a fuller picture of the genetic architecture of ANM across the globe. This requires research into statistical models that are capable of learning both shared and unique age-dependent effect sizes across populations, and it requires large-scale data to be collected from worldwide populations.
In summary, we propose an analysis approach for GWASs of age-at-onset phenotypes using a two-stage mixed linear-association model, where marker effect sizes are estimated using a CAMP model. Our approach provides a better understanding of the genetic basis of ANM and applies to any form of time-to-event phenotype.
This project was funded by an SNSF Eccellenza grant to M.R.R. (PCEGP3-181181) and by core funding from the Institute of Science and Technology Austria. K.L. and R.M. were supported by the Estonian Research Council grant 1911. Estonian Biobank computations were performed in the High-Performance Computing Center, University of Tartu. We thank Triin Laisk for her valuable insights and comments that helped greatly. We would like to acknowledge the participants and investigators of UK Biobank and Estonian Biobank studies. This project uses UK Biobank data under project number 35520.
S.E.O. and M.R.R. conceived and designed the study. S.E.O. and Z.K. conceived and derived the significance testing. S.E.O., L.D., M.P., and M.R.R. conducted the analysis. K.L., K.F., and R.M. provided study oversight and contributed data. S.E.O. and M.R.R. wrote the paper. All authors approved the final manuscript prior to submission.
M.R.R. receives research funding from Boehringer Ingelheim. S.E.O. is an employee of MSD at the time of the submission; contribution to the research occurred during the affiliation at the University of Lausanne.
Published: August 4, 2023
Sven E. Ojavee, Email: sven.ojavee@gmail.com.
Matthew R. Robinson, Email: matthew.robinson@ist.ac.at.
Age-specific summary statistic estimates are released publicly on Dryad: https://doi.org/10.5061/dryad.nvx0k6dx5. The BayesW model was executed with the software Hydra, with full open source code available at https://github.com/medical-genomics-group/hydra.^29^ The scripts used to execute CAMP model are available at https://github.com/svenojavee/CAMP. R version 4.2.1 is available at https://www.r-project.org/.
Age-specific summary statistic estimates are released publicly on Dryad: https://doi.org/10.5061/dryad.nvx0k6dx5. The BayesW model was executed with the software Hydra, with full open source code available at https://github.com/medical-genomics-group/hydra.^29^ The scripts used to execute CAMP model are available at https://github.com/svenojavee/CAMP. R version 4.2.1 is available at https://www.r-project.org/.