Human longevity: 25 genetic loci associated in 389,166 UK biobank participants

GWAS of parental attained age

We examined associations with both parents’ combined attained ages (including ages of living parents), using Martingale residuals from Cox’s proportional hazards regression models, following the approach used by Joshi et al. [13]. In this analysis we found 10 loci associated at genome wide significance (p<5x10^-8) (Table 2; Supplementary Table 1 for details of all significant variants; full results available to download https://doi.org/10.6084/m9.figshare.5439382.v1).

Associations with APOE related rs429358 (p=5x10^-22) variant were confirmed, as were associations with the CHRNA3/5 (nicotinic acetylcholine receptor) locus (rs1317286, 33.2% minor allele frequency, p=1.2x10^-26). Associations were also found for six other loci already implicated in disease specific GWAS; the strongest association was with rs55730499 (8% frequency, p=1.7x10^-18) in an intron of LPA gene, associated with lipoprotein A and LDL cholesterol levels, plus CAD [14]. Variants in the 9p21 region, which includes CDKN2B-AS1 (ANRIL) and CDKN2A/B, were also associated (lead SNP rs1556516, 49.8% frequency, p=4.7x10^-16); this locus is involved in cellular senescence and associated with myocardial infarction and peripheral vascular disease [14]. Genetic variants across the 9p21.3 locus have previously also been associated with type-2 diabetes and cancers [14] and in longevity/survival analyses (see below), but only the CAD-associated loci were associated (p<5x10^-8) with parent’s lifespan (Figure 2; Supplementary Table 2).

Figure 1. Manhattan plots for GWAS of parent’s lifespan. (A) combined parents' attained age GWAS including living parents, (B) combined age at death of mother and father, (C) both parents reached top 10% of survival (mothers reached ≥90 years, father reached ≥87 years). Red points are variants (p<1x10^-5) previously associated with other traits. The blue line indicates p=5x10^-8. The Y-axis has been capped at 20 (p=1x10^-20) to aid in visualising the results.

Figure 2. LocusZoom plot of the CDKN2A/B/ANRIL gene cluster and known disease variants. Known disease-associated genetic variants have been highlighted in the 9p21.3 region, in addition to the lead SNP identified in our analysis. Association p-values indicate the strength of association with parents' attained age. The variants known to influence coronary artery disease (CAD) and related-traits are those most-associated with parents’ attained age. Variants previously associated with cancer are markedly less associated with parents’ attained age, with type-2 diabetes variants not associated (p>0.05). Details, including references, are included in Supplementary Table 2.

Associations with parents' attained age were also present for rs7137828 (48.3% frequency, p=3.4x10^-14) near ATXN2, coding for Ataxin-2 which is involved in endocytosis and is associated with blood traits and autoimmune disease [14]. SNPs in LD (R²=0.91) with rs7137828 include rs3184504 (48% frequency, p=4x10^-8), a missense variant in the SH2B3 gene, associated with blood pressure, maternal birth weight, TNF-α levels, colorectal cancer, Inflammatory Bowel Disease, LDL cholesterol levels and stroke [14]. Other known loci included rs28383322 (located between HLA-DRB1 and HLA-DQA1, p=5.3x10^-11), previously associated with atopic dermatitis; rs17514846 (located in an intron of FURIN, p=7.1x10^-10), previously associated with blood pressure; and rs602633 (located between CLESR2 and PSRC1, p=2.7x10^-8), previously associated with CAD and LDL [14].

Two associations were discovered for novel variants not previously linked in GWAS to other traits or diseases; rs7844965, located in an intron of EPHX2 (23.2% frequency, p=7.7x10^-9), and rs61978928, located in the 3’ UTR of gene PROX2 (31.2% frequency, p=2x10^-8).

A parents' attained age genetic risk score in independent cohorts

We calculated a genetic risk score (GRS) combining the 10 attained parental age associated variants, weighted by the estimated effect sizes. Genotyped samples of European descent from the Health and Retirement Study (HRS, n=12,940 in genetic analyses) and the Wisconsin Longitudinal Study (WLS, n=8,952 from 6,683 families) were too small to test individual variant associations: the sample size needed for 80% power to detect an association between the APOE variant and parents' lifespan at p<0.05 (R²=0.0005 in UK Biobank) is 15,598, but much larger samples were required for the smaller-effect variants discovered in UK Biobank. We therefore tested the parents’ attained age GRS; the sample size needed for an 80% chance of detecting an association of this magnitude (R²=0.002 in UKB) with p<0.05 was 3,924. In Cox’s proportional hazards survival models against both dead and alive parents attained ages in HRS (n=12,345, number of events=9,809), increasing GRS was associated with longer survival (mortality Hazard Ratio=0.961, 95% CIs=0.942 to 0.980, p=6x10^-5, adjusted for offspring age, sex, and genetic principal components 1-5). In the WLS (n=7,437, number of events=6,683, errors clustered by family), increasing GRS was associated with longer survival (mortality Hazard Ratio=0.963, 95% CIs=0.938 to 0.988, p=4x10^-3, adjusted for age at saliva collection, sex, and genetic principal components 1-5).

Functional implications and pathways analysis

We investigated effects of the parents' attained age SNPs on expression of genes in multiple tissues, such as whole blood [15], using the FUMA tools [16]. Expression of 53 genes were associated, including 8 of the 10 nearest genes to the variants identified above (Supplementary Table 5 and Supplementary Table 6).

In MAGENTA gene-set analysis, only two pathways were significantly enriched after correction for multiple testing, both related to nicotinic acetylcholine receptor binding and downstream events (Supplementary Table 4).

Sensitivity analyses of combined parents' attained age analyses

We generated QQ plots for all GWAS (Supplementary Figure 2) and observed some genomic inflation (Lambda GC for parents' attained age = 1.199), consistent with human longevity being highly polygenic, but perhaps also representing some effect of relatedness in inflating results. The BOLT-LMM software used for analysis accounts for relatedness in the methods, but we investigated further; the 10 variants associated with attained parental age were tested in Cox proportional hazards (CoxPH) models against parents’ attained age, to assess whether the associations observed with rank-normalized Martingale residuals were robust. These models also excluded related (≥3^rd-degree) participants, leaving a reduced sample size of n=326,864, including n=153,101 with censored data due to at least one parent still being alive. Five of the 10 loci remained genome wide significant (p<5x10^-8) in these models; APOE, CHRNA3, LPA, ANRIL, and ATAXN2/SH2B3 (Supplementary Table 3) with the remaining five directionally consistent and associated at p≤3.1*10^-6.

We performed sensitivity analyses only including parents who had died (Figure 1B; Table 2), as the ages of living parents at the UKB study baseline are clearly an incomplete measure of their longevity. This analysis included a smaller sample size of n=208,118 UKB participants with both parental ages of death known. The same five of the ten loci reported above remained genome-wide significant (APOE, CHRNA3, LPA, ANRIL, and SH2B3) with the remaining loci remaining nominally significant. Two additional variants not identified in the previous analysis emerged as candidates for follow-up: rs1627804 (32% frequency, p=4x10^-8), in an intron of BEND3, a transcriptional repressor which associates with the NoRC (nucleolar remodeling complex) complex and playing a key role in repressing ribosomal DNA transcription; and rs146254978 (2% frequency, p=4.6x10^-8), located in an intron of both FPGT and TNNI3K.

GWAS of both parents achieving longest 10% survival

We next tested binary associations with having both parents surviving to the longest 10% of lifespans, using sex-specific age-at-death cut-offs (father ≥87, mother ≥90, including living parents aged over these age cut-offs), to assess whether associations were attributable to longer lives rather than early deaths. This analysis confirmed associations (Figure 1C; Table 1) with variants in the APOE, CHRNA3, LPA and CDKN2B-AS1, with a suggestive association with ATXN2 (p=3.4x10^-6). Two further loci just reached significance in this analysis: rs28926173 (4.7% frequency, p=2.3x10^-8) in an intron of gene MC2R, and rs139137459 (0.4% frequency, p=2.7x10^-8) in an intron of USP2-AS1 (near USP2). The MC2R variant was associated at p=6.4x10^-6 with parents’ attained age, but the USP2-AS1 variant had a low population frequency (0.4%) and was not significant in the continuous trait analyses.

Genetic risk score associations with centenarian status and offspring health status in UKB

There were 691 UKB participants with a centenarian mother (the longest lived 0.1% of maternal deaths in the sample), and 747 with a father aged ≥98 years (the corresponding age cut-off of longest lived 0.1%); only five participants had both parents in the longest lived 0.1%. Although this relatively small group (0.3% of the full attained age sample) were included in the discovery analysis, we tested the GRS for association with having a centenarian parent (in the un-related UK Biobank participants, n=1,181), against controls whose (both) parents died before age 80 (n=66,280). Increasing attained age GRS was associated with increased odds of having a centenarian parent (OR per standard deviation of GRS=1.193, 95% CIs=1.130 to 1.260, p=3.0x10^-10).

GWAS of mothers' and fathers' age at death separately

We tested for parental gender differences by conducting separate GWAS for mothers' and fathers' attained ages (including living parents). Eleven additional genetic variants were identified, three only associated with mother’s age at death, eight associated with father’s age at death (Table 2). Results included variants previously associated with the autoimmune diseases psoriasis and idiopathic membranous nephropathy (for mothers), and an independent nicotine-dependence locus on chromosome 20 for fathers. Novel gene-longevity associations included genes MICA/B (a locus previously linked to auto-immune conditions rheumatoid arthritis and multiple sclerosis), TOX, ZW10, SEMA6D, EGLN2/CYP2A6, EXOC3L2/MARK4 and C20orf187.

Comparison with recently reported GWAS analyses of longevity

A recent large analysis by Joshi et al. [2], with up to 294,998 participants in their analysis of European parents' attained age, included the first third of UK Biobank genotyping data plus 26 other cohorts. For survival past age 40 years, associations were confirmed for APOE and CHRNA3/5, plus novel associations with overall longevity (deaths or survival to age 40 years and over) with HLA-DQA1/DRB1 and LPA, and with replication of CDKN2A/B and SH2B3. While these results overlap with our attained age analysis, they also reported an association with FOXO3A variants, but we found no significant associations for this locus (p>0.01) (see Supplementary Table 2 for details). No analyses of longer or extreme survival were reported.

McDaid et al. [6] undertook a Bayesian association scan of general lifespan (in the first third of UK Biobank), weighted for 58 disease related GWAS studies, and identified 16 individual variants at genome wide significance. Of these, four (LPA, CDKN2BAS, CHRNA5 and APOC1) reached conventional GWAS significance in our analyses, and CELSR2 only reached significance for mothers' age at death.

While the above studies focused mainly on overall lifespan, for aging research survival to extreme ages is of special interest. In a powerful study of extreme longevity, Fortney et al. used samples including 801 centenarians and 5406 subjects aged 90 plus. Disease GWAS results were used to weight genetic variants, using a false discovery rate of 10% [5], at which associations with APOE, SH2B3/ATXN2, ANRIL, LPA, and HLA were reported (although only APOE reached conventional genome wide significance). Fortney et al. also identified variants at KCNT2 and FADS1, which were not associated with parents’ attained age in our analysis, with the ABO SNP (rs514659) only just reaching nominal significance (p=0.015).

A previous meta-analysis comparing older individuals (age 80 or 90+) to younger controls reported a genome wide association at the 5q33.3 locus, but no significant associations were found for the other two loci (p>0.01).

Various other variants, including in mitochondrial DNA, have been suggested as associated with longevity. In our GWAS we tested associations with mitochondrial variants, plus variants on the X and Y chromosomes, using the directly genotyped microarray data: none were significantly associated (p>1x10^-5) with the studied phenotypes.

Disease-associated genetic variants; candidates for future work

We investigated genetic variants associated with parents' attained age (p<1x10^-5) that have previously been associated with other traits for promising candidates for future work and identified 12 independent signals (Table 3). These included SNPs associated with BMI (rs11126666 [17],), breast cancer (rs11552449 [18],), educational attainment (rs7593947 [19],), and kidney function (estimated Glomerular Filtration Rate, rs4014195 [20],), in addition to a variant in IL6R (rs6689306) associated with CAD.

Future work

Replication and extension of the novel associations in other large cohorts and in different ancestries is needed, especially for variants narrowly achieving genome wide significance. The associations identified should be examined in future studies of age-related traits and disease. More work is needed to characterize the variants identified for biological functions; this could include epigenetic and gene expression studies, such as analysis of allele-specific expression imbalance in heterozygote individuals to determine the gene of action for genetic variants [35]. The variants identified in our analysis may provide useful information for designing targeted prevention or precision treatments to promote human longevity, perhaps particularly related to targets in inflammatory pathways and cardiovascular disease.

Parental age at death and longevity phenotypes

Participants were asked the age at which their parents had died (or their current age if still alive). Analyses were performed separately on mother’s age at death and father’s age at death, and also on a combined phenotype. To reduce the effect of higher ages at death of mothers (compared to the fathers) we first z-transformed the mothers and fathers age at deaths before combining the z-scores into a single summed phenotype. Offspring of parents who died prematurely were excluded because the cause of death of the participant’s parents was not asked, so we could not exclude accidental deaths explicitly. To determine the age at death considered to be too young to be part of the normal distribution of age-related mortality, we used previously described methods to define a normal distribution of ages at death for mothers and fathers separately, as described previously [9]). In brief, non-linear least square regression models were used to fit a normal curve on age at death data from the modal (M) age at death to the oldest ages; the younger (left) half of the curve is then imposed, to avoid the severe skewing of the distribution caused by the excess of deaths below the modal age. The analysis identified the following early death cut-points (based on modal age at death minus 1 Standard Deviation): mothers 57 years, fathers 46 years.

In addition, combined parents’ age (including alive parents) were derived from Cox’s proportional hazards regression models using Martingale residuals, as described by Joshi et al. [13]. Events were defined as the both parents having died. We also defined a binary phenotype by determining the top 10% of age at death for mothers (≥90 years) and fathers (≥87 years) separately. Controls were participants where both parents died before the age of 80. We were able to increase the number of “long lived” parents by including alive participants beyond the cut-off.

UK Biobank

Between 2006 and 2010, 503,325 volunteers (aged 45 to 69 years old) were recruited from across the United Kingdom to the UK Biobank study [36]. Of these, 389,166 participants met the inclusion criteria for at least one analysis: participants with complete genetics data and either date of death of their parents or their current age if alive. Participants were excluded if they reported themselves as adopted, or whose parents died more than 1 standard deviation below the modal age of death (mothers <57 years, fathers <46 years). Several longevity phenotypes are defined below based on the age at death of the participant’s parents.

Genetic data was available on 488,377 UK Biobank participants after genotype calling and quality control performed centrally by the UK Biobank team [37]. We selected 451,447 participants identified as ‘white European’ through self-report and verified through principal components analysis based on genotypes. Briefly, principal components were generated in the 1000 Genomes Cohort using high-confidence SNPs to obtain their individual loadings. These loadings were then used to project all of the UK Biobank samples into the same principal component space and individuals were then clustered using principal components 1 to 4. Related individuals were identified through kinship analysis, although these participants were included in genome-wide analysis using BOLT-LMM (see below), those related to the third-degree or closer were excluded in sensitivity analyses.

Imputation of 39,235,157 genetic variants from the Haplotype Reference Consortium panel was performed using IMPUTE4 centrally by the UK Biobank team [37]. After filtering for variants with MAF ≥0.1%, missingness <1.5%, imputation quality >0.1 and with Hardy-Weinberg equilibrium (HWE) P>1x10^-6 within the white British participants 11,516,125 imputed autosomal variants were eligible for the analyses.

We also utilized data directly from the microarrays for variants on the X (n=19,381) and Y (n=284) chromosomes, and on the mitochondrial genome (n=135), which were unavailable in the imputed dataset.

Health and Retirement Study (HRS)

In HRS v2 data, phase 1-3 samples were combined (n=15,708) and SNPs were kept if no discordant genotype call among duplicate samples (2,315,518 SNPs). Related or duplicated participants and genotyping controls were removed, as well as those with a missing call rate equal to or greater than 2%. After the sample quality control, 15,454 samples were kept and 84% of them (n=12,940) were of European descent.

2,075,208 out of 2,315,518 SNPs passed the composite filters by the University of Washington (Hardy-Weinberg equilibrium p-value < 1e-4, missing call rate ≥ 2%, sex difference in allelic frequency ≥ 2% and sex difference in heterozygosity > 0.3) and had a minor allele frequency greater than zero. As recommended, SNPs with a minor allele frequency smaller than 2% further were removed. At the end, 1,437,608 SNPs were carried forward for SNP pruning. The purpose is to choose less correlated SNPs for principal components analysis.

The SNP pruning was performed following the procedure: a) consider a window of 50 SNPs, b) calculate LD between each pair of SNPs in the window, c) remove one of a pair of SNPs if the LD is greater than 0.5, d) shift the window 5 SNPs forward and repeat the procedure. 514,246 SNPs remained after the SNP pruning and were used to run principal components analysis. As a result, the top 10 principal components were used to adjust for ancestry in a regression model. Both SNP pruning and principal components analyses were performed in PLINK1.9

The samples in the genetic data were mapped to the phenotype data of HRS to retrieve parental longevity phenotypes. We used the phenotype data from the most recent wave or visit of HRS where parental data was missing if the parent died before the last visit. To recover the data, we sequentially looked back at the last visit until the parental data was available. Father’s and mother’s data therefore could be from different visits and the participant’s age varied with phenotypes. Mothers dying before 61 years and fathers dying before 46 years were excluded from the analysis, based on the same analysis of distribution of ages of death applied in UK Biobank. Parent’s age was set to missing if he or she was still alive and younger than the child’s age +10.

The genetic risk score analysis was conducted in a Cox regression model (against combined parents' attained age, alive or dead: derived by summing the z-transformed mothers and fathers' attained ages) with adjustment for age and sex of the participant and the top 10 principal components.

Wisconsin Longitudinal Study (WLS)

In the WLS saliva samples were collected by mail in 2006-7 and during in-person interviews in 2011 for those who did not mail back samples. Genotyping on the Illumina HumanOmniExpress array was completed at Johns Hopkins’ Center for Inherited Disease Research (CIDR) and data cleaning and imputation was performed in collaboration with the Genetic Analysis Center at the University of Washington. After quality control, there are a total of 9,012 respondents and 688,596 SNPs available for analysis. The sample has 4,601 singletons and 2,263 families with more than one member genotyped. Only 66 respondents are not of European descent.

Imputation of 33,549,349 genetic variants from the 1000 Genomes Project reference panel (Phase 3) was performed using IMPUTE2, with 686,143 study SNPs serving as the imputation basis. Of the 10 variants used to construct the GRS, two were directly genotyped and eight were imputed.

Parent’s age was constructed by subtracting the year of birth from the year of death. When the latter was not available, the year of the most recent interview was used instead and the parent was classified as not dead. Mothers dying before 56 years and fathers dying before 46 years were excluded from the analysis, based on the same modelling of distributions of ages of death applied in UK Biobank. Parent’s age was also set to missing if the disagreement between sibling reports exceeded 10 years. The analysis phenotype was the standardized sum of both parents' ages and was set to missing if either element was not available.

The genetic risk score analysis was conducted in a Cox regression with censored data, adjusting for age at saliva collection, sex of the participant and the top 5 principal components. Standard errors were clustered at the family level. Limiting the sample to the primary WLS respondents had a negligible effect on the results.

Power Calculations

R package `pwr` and function ‘pwr.f2.test’ was utilized for power calculations of regression analysis. The numerator degrees of freedom were set to 1, significance level to 0.05, and power to 0.8.

Genome-Wide Association Study

We used BOLT-LMM v2.2 to model the associations between imputed variants (dosages) and each phenotype [38], which uses a linear mixed effects model approach. We looked at the results for variants with imputation quality >0.1, HWE p-values >1x10^-6 and minor allele frequencies >0.1% in the white/British subset used for all analyses. For variants on the X, Y and mitochondrial chromosomes only in the directly-genotyped data we used Plink (v1.9) [39] in linear (additive) or binary (fisher) models, as appropriate, adjusted for the same covariates as above including the first 5 principal components from FlashPCA.

Pathways Analysis and Functional Annotation

We utilized FUMA (Functional Mapping and Annotation of Genome-Wide Association Studies) to investigate the genetic variants associated with parents’ attained age [16]. In particular, FUMA performs gene-set analysis (using MAGMA) to identify pathways enriched amongst the significant genes (weighted by the SNP-associations in proximity to them). In addition, FUMA searched eQTL databases to identify SNPs that significant alter the expression of genes in various tissues.