Exome-Wide Association Study Identifies FN3KRP and PGP as New Candidate Longevity Genes.

Despite enormous research efforts, the genetic component of longevity has remained largely elusive. The investigation of common variants, mainly located in intronic or regulatory regions, has yielded only little new information on the heritability of the phenotype. Here, we performed a chip-based exome-wide association study investigating 62 488 common and rare coding variants in 1248 German long-lived individuals, including 599 centenarians and 6941 younger controls (age < 60 years). In a single-variant analysis, we observed an exome-wide significant association between rs1046896 in the gene fructosamine-3-kinase-related-protein (FN3KRP) and longevity. Noteworthy, we found the longevity allele C of rs1046896 to be associated with an increased FN3KRP expression in whole blood; a database look-up confirmed this effect for various other human tissues. A gene-based analysis, in which potential cumulative effects of common and rare variants were considered, yielded the gene phosphoglycolate phosphatase (PGP) as another potential longevity gene, though no single variant in PGP reached the discovery p-value (1 × 10E-04). Furthermore, we validated the previously reported longevity locus cyclin-dependent kinase inhibitor 2B antisense RNA 1 (CDKN2B-AS1). Replication of our results in a French longevity cohort was only successful for rs1063192 in CDKN2B-AS1. In conclusion, we identified 2 new potential candidate longevity genes, FN3KRP and PGP which may influence the phenotype through their role in metabolic processes, that is, the reverse glycation of proteins (FN3KRP) and the control of glycerol-3-phosphate levels (PGP).

Human longevity is a complex phenotype influenced by both genetic and environmental factors, which furthermore interact, for example, via epigenetic changes (1). Heritability estimates for longevity range from 12% (2) to 30% (3). However, in the oldest-old, the genetic influence appears to be even higher and it has been suggested to be as high as 48% (4). Candidate studies have yielded many potential longevity associations, but most of them have not been replicated in independent investigations. To date, variation in only 4 loci has been confirmed to influence longevity across populations: APOE (5,6), FOXO3 (7,8), the 5q33.3 locus (9,10), and CDKN2B-AS1 (11–13). The first 3 reached genome-wide significance in large genome-wide association studies (GWAS) (9,14). However, the identified variants together explain only a small proportion of the longevity heritability. Thus, novel approaches are needed to identify additional loci involved in the phenotype.

Longevity studies have mostly been conducted following the common disease/common variant hypothesis, which is based on the assumption that the probability of becoming long-lived depends on a small number of single-nucleotide variants (SNVs) that occur at high frequency in all populations. Yet, it has been estimated that common variants explain only 5%‐10% of the heritability of complex traits (15). Therefore, rare variants (minor allele frequency [MAF] < 0.01) may play an important role in the genetic regulation of longevity, but they are often not covered by genome-wide genotyping arrays or imputations and are therefore less studied (14). Common variants often reside in regulatory/intronic regions (16), whereas rare variants are more likely to be present in exons. Therefore, an effective way to target such rare variants is through exome-based genotyping or through whole genome or exome sequencing. Some studies have used these methods to investigate coding region variants associated with longevity (eg, (17–20)). These approaches yielded important results, not only in terms of new longevity variants in, for instance, CLEC3B and HLA-DQB1 (17,19), but also provided novel analysis methods that could extend the discovery spectrum of loci and genes influencing complex traits, like the combined analysis of multiple variants (11).

Here, we performed an exome-wide association study in a large German longevity cohort that comprised more than 1200 long-lived individuals (LLI) and 6941 younger controls using the Illumina Infinium HumanExome BeadChip that covers both rare and common variants to identify novel variants/genes associated with longevity.

Method

Study Populations

In the following, a brief overview of the populations used in this study is presented. For a more extensive description, please refer to the Supplementary Material. Our discovery cohort comprised 1248 German LLI (male/female ratio: approximately 1/3; age range: 94–110 years; mean age: 99 years) and 6941 younger controls (age < 60 years). For replication, we investigated a Danish data set with 1002 cases (male/female ratio: 1/3; age range: 90.0–102.5 years; mean age: 97.4 years) and 738 controls (mean age: 66.3 years), and a French cohort with 1264 LLI (male/female ratio ~ 1/4.5, age range 91–115+ years; mean age 102.4 years;) and 1830 younger subjects (age < 62 years). All individuals belonged to the top 1% of the survivors of their respective birth cohorts.

Exome Chip Genotype Calling and Quality Control in the German Study Population

The samples were genotyped on the Infinium HumanExome-12v BeadChip (244,770 SNVs) (Illumina Inc., San Diego, CA; Supplementary Figure 1). Quality control was done with the software Plink 1.9 (21). In total, 224 samples were removed because they had failed one or more of the following inclusion criteria: concordant sex information, missing genotype < 3%, heterozygosity rate greater or lower than ± 4SD from the mean, and no relatedness of individuals. Relatedness was estimated using the identity by descent (IBD) metric (22). In case of relatives (IBD > 0.1875; halfway between the expected IBD for third- and second-degree relative (22)), only one individual was included in the analysis. SNVs were excluded if the missing rate was too high (> 3%) or, for common SNVs in the control sample, if they deviated from the Hardy–Weinberg equilibrium (HWE) (p < .0001). Due to the lack of sufficient statistical power, SNVs with an extremely low minor allele frequency (MAF < 0.003) were removed prior to the single-variant analysis. This resulted in 62 488 remaining variants.

Population substructure was evaluated with the principal component analysis (PCA) using a common set of independent markers (HapMap3 ancestry set from 4 ethnic populations; 16 782 SNVs). The principal components (PC) were calculated using PLINK 1.9 (21). The first 5 PCs were selected for further association analysis to adjust for population stratification. Additionally, stratification outliers were identified based on the local outlier factor (LOF > 1.7) (23) and excluded to mitigate population stratification.

Association Analysis in the German Sample

Single-variant association analysis was performed using the logistic regression test implemented in PLINK 1.9 (21) assuming an additive genetic model and additional influence variables, namely the 5 PCs and sex (Supplementary Figure 1). Candidate longevity SNVs with a discovery p-value < 1 × 10E−04 were selected for replication. Relaxation of the significance threshold allows the identification of longevity SNVs that usually have small effects on the phenotype (24). To identify additional association signals and to test for independency of the newly identified SNVs from the effects of the known longevity-associated locus TOMM40/APOE/APOC1 (9,11,25), a conditional association test was performed using logistic regression in PLINK 1.9 (21), adjusting for the SNVs rs2075650, rs4420638, and rs769449. The Sanger imputation service (http://www.sanger.ac.uk/science/tools/sanger-imputation-service) and the 1000 Genomes phase I v.3 reference panel were used to enrich the pool of common SNVs.

In addition to the single-variant association testing, a gene-based analysis was performed. All SNVs that passed quality control (112,977 SNVs) were used to generate the gene set (Supplementary Figure 1). Potential cumulative effects of rare and common variants with longevity were tested using both burden and non-burden (ie, SKAT) approaches (application of the RC-SKAT algorithm from the R-package SKAT (26)). Both approaches were considered because burden tests were shown to perform better when multiple variants in a region are causal and influence the phenotype in the same direction (27), while non-burden tests, like SKAT (28), are more advantageous when SNVs in the region interact or show opposing directions of effect (29). The overall effect of rare and common variants in a gene was evaluated based on the adaptive sum test (26) in combination with either burden or SKAT. Concordantly, the gene-based Bonferroni-corrected p-value threshold was based on the number of gene sets tested: p < 3.3 × 10E−06 (significance threshold of 0.05/14 790 genes = the number of genes with SNVs that remained after quality control).

Epistatic interaction between pairs of candidate SNVs was assessed using multifactor dimensionality reduction (MDR) analysis (30) correcting for confounders such as sex and the known longevity-associated locus TOMM40/APOE/APOC1 (9,11,25). The entropy-based clustering algorithm used by MDR calculates case–control ratios for each of the possible multilocus genotypes; therefore, if a genotype combination is more abundant in cases than in controls, it is considered as a high-level interaction. The MDR interaction model describes the percentage of entropy (information gain) by each SNV or SNV interaction. Positive values of entropy indicate synergistic or nonadditive interactions, while negative entropy values indicate redundancy between the markers or lack of synergy between the interacting markers. The model (ie, combination of SNVs) that appeared most consistently among replicates in the training balance accuracy level was considered the best model. Significance was calculated via permutation tests using 1000 permutations and a significance level of α = 0.05. MDR was implemented employing the open-source MDR software package version 3.0.2 (https://www.mybiosoftware.com/mdr-2-0-multifactor-dimensionality-reduction.html) (31).

The functional implications of the longevity variants were assessed using our previously published gene expression data from whole blood samples of 55 LLI and 73 control individuals (independent cohort from Germany) (32) as well as publicly available data from the Blood eQTL browser (https://genenetwork.nl/bloodeqtlbrowser/) (33) and the GTEx eQTL (https://www.gtexportal.org/home/; accessed April 5, 2019).

Genotyping in the Replication Cohorts

The Danish LLI were genotyped using the Illumina HumanOmniExpress Array (Illumina Inc.). Preimputation quality control included filtering of SNVs on genotype call rate <95%, HWE p < 10E-04, MAF < 1% and filtering of individuals on sample call rate <95%, relatedness, and sex mismatch. After imputation to the 1000 Genomes phase I v.3 reference panel, genotype probabilities were converted to hard-called genotypes in PLINK 1.9 (using a cutoff of 90%) (21). For the Danish controls, data for the SNVs rs1063192, rs1046896, and rs13119846 were extracted from quality-controlled genotype data (as above) created using the Illumina Infinium PsychArray (Illumina Inc.).

French individuals were genotyped by TaqMan (Thermo Fisher Inc., Waltham, MA) on a 7900HT Fast Real-time PCR System (Thermo Fisher Scientific Inc., Waltham, MA). Association analysis for the Danish and French cohort was performed with logistic regression accounting for sex as covariate using PLINK 1.9 (21).

Availability of Data and Materials

All German samples and information on their corresponding phenotypes were obtained from the PopGen Biobank (Schleswig-Holstein, Germany). The data can be accessed through a Material Data Access Form (http://www.uksh.de/p2n/Information+for+Researchers.html).

Results

Single-Variant Association Analysis Reveals an Association of rs1046896 in FN3KRP With Longevity

To identify new common and rare genetic variants associated with longevity with moderate to high effect size, a chip-based exome-wide association study was conducted using the Illumina Infinium HumanExome BeadChip. This array covers rare and common variants in a ratio of approximately 8:1 (34). In total, 1248 German LLI and 6941 younger controls were included in the study (Supplementary Figure 1). The single-variant analysis was performed based on 62 488 nonimputed SNVs, 1212 LLI, and 6762 younger controls who remained after quality control (Supplementary Figures 2 and 3). The single-variant association approach for rare and common variants yielded 11 candidate SNVs (p < 1 × 10E−04; Table 1). The best association signal and the only one reaching exome-wide significance (except for the TOMM40/APOE/APOC1 region, see below) was obtained for rs1046896 in the gene FN3KRP (MAF = 0.32, p = 7.40 × 10E−07; Table 1; Figure 1; Supplementary Table 1). In addition, we observed a longevity association for rs1063192 in the CDKN2B-AS1 region (p = 2.99 × 10E−05; Table 1; Figure 1). Apart from these SNVs, we identified rs2075650, rs4420638, and rs769449 in the TOMM40/APOE/APOC1 region that is well known to be negatively associated with longevity (5,9,11,25).

Table 1.

Association Statistics for the 11 Longevity-Associated SNVs Identified by the Single-Variant Association Approach in the Whole German Study Population

			MAFa			Basic Association Test		Conditional Analysis
SNV	Gene	Chr	LLI	C	MA	ORb		ORb
						[95% CI]c	p d	[95% CI]c	p d
rs769449	APOE	19	0.056	0.109	A	0.48 [0.40–0.58]	7.77E−15	-	-
rs4420638	APOC1	19	0.109	0.169	G	0.60 [0.52–0.69]	3.55E−13	-	-
rs2075650	TOMM40	19	0.109	0.147	G	0.70 [0.61–0.80]	3.51E−07	-	-
rs1046896	FN3KRP	17	0.276	0.324	T	0.78 [0.71–0.86]	7.40E−07	0.77 [0.70–0.85]	2.32E−07
rs55882518	NOTCH3	19	0.013	0.005	T	2.69 [1.73–4.18]	1.07E−05	2.77 [1.78–4.30]	6.23E−06
rs13119846	TMEM131L	4	0.486	0.438	C	1.22 [1.11–1.33]	1.73E−05	1.21 [1.11–1.32]	2.45E−05
rs1063192	CDKN2B-AS1	9	0.482	0.439	G	1.21 [1.11–1.32]	2.99E−05	1.22 [1.12–1.33]	1.08E−05
rs184214819	SPZ1	5	0.009	0.003	A	3.01 [1.77–5.11]	4.72E−05	2.79 [1.64–4.74]	1.52E−04
rs200956599	SKOR1	15	0.014	0.006	T	2.34 [1.55–3.54]	5.24E−05	2.29 [1.52–3.47]	8.65E−05
rs63750412	GRN	17	0.008	0.002	T	3.57 [1.93–6.61]	5.18E−05	3.34 [1.80–6.18]	1.29E−04
rs1790706	DSC2	18	0.159	0.189	A	0.79 [0.70–0.89]	9.59E−05	0.78 [0.69–0.88]	3.35E−05

Note: APOC1 = apolipoprotein C1; APOE = apolipoprotein E; CDKN2B-AS1 = cyclin-dependent kinase inhibitor 2B antisense RNA 1; DSC2 = desmocollin 2; FN3KRP = fructosamine 3 kinase-related protein; GRN = granulin precursor; NOTCH3 = notch 3; SKOR1 = SKI family transcriptional corepressor 1; SPZ1 = spermatogenic leucine zipper 1; TMEM131L = transmembrane 131 like; TOMM40 = translocase of outer mitochondrial membrane 40.

C = controls; Chr = chromosome; LLI = long-lived individuals.

aMinor allele frequency, MAF; the definition of the minor allele (MA) is based on controls.

bOdds ratio for longevity, OR; based on the MA in controls.

c95% confidence interval, 95% CI; CI for the OR.

dAllelic p-values, calculated from logistic regression.

Figure 1.

(A) Circular Manhattan plot summarizing the findings from the single-variant analysis. The inner plot represents the basic association results, the outer plot the association results after conditioning on the longevity-associated locus TOMM40/APOE/APOC1. The y-axis shows the -log(p-value), while the dotted line depicts the p-value threshold (1 × 10E−05). SNVs with p < 1 × 10E−05 are shown as dots. (B,C). Regional plots for FN3KRP and CDKN2B-AS1 candidate variants.

The effects of the candidate variants were investigated for independency of the TOMM40/APOE/APOC1 locus by a conditional association test. The results of the conditional analysis confirmed the independency of the longevity association of rs1046896 (FN3KRP) as well as of the variants rs55882518 (NOTCH3), rs1063192 (CDKN2B-AS1 region), rs1319846 (TMEM131L), and rs1790706 (DSC2) (Table 1). Additionally, an association analysis with the centenarian subpopulation (n = 599 individuals ≥100 years) was performed. The association analysis yielded higher ORs for rs1046896 (FN3KRP), rs55882518 (NOTCH3), rs1063192 (CDKN2B-AS1 region), and rs184214819 (SPZ1) compared with the analysis using the whole study population (Supplementary Table 2). This effect has been reported before (35) and is explained by the larger genetic influence with increasing age. However, the centenarian subset comprised substantially fewer individuals and apart from rs1046896 (FN3KRP), the longevity associations did not reach exome-wide significance (Supplementary Table 2).

Since the top hit, rs1046896, was located in a 3’UTR and might therefore affect gene expression in an allele-dependent manner, we investigated whether this SNV or SNVs in high LD (r2 > 0.8 based on the CEU subpopulation from the 1000 Genomes Project Phase 3; Supplementary Table 3) influence local or distant gene expression. In our previously published whole blood transcriptome data from 55 LLI and 73 control individuals (independent cohort from Germany) (32), we observed significant cis-eQTL associations of rs1046896 (and high LD SNVs) with the expression of FN3KRP. FN3KRP gene expression was higher in the presence of the rs1046896 longevity allele C (major allele) (Figure 2A). Using publicly available data from the Blood eQTL browser (https://genenetwork.nl/bloodeqtlbrowser/) (33) and the GTEx eQTL database (https://www.gtexportal.org/home/; accessed April 5, 2019), we confirmed the allele-dependent expression of FN3KRP in several tissues (eg, brain regions, testis, pancreas; Supplementary Table 4; Figure 2B). For the expression of FN3KRP, we observed a C-allele dose effect (TTFN3KRP expression even when they were homozygous for the T allele (Figure 2C).

Figure 2.

rs1046896-C allele-dose influence on FN3KRP mRNA expression. (A) Boxplots of FN3KRP expression according to C allele-dose shown separately for males (M) and females (F). (B) cis-eQTL effect of rs1046896-C on FN3KRP expression in several human tissues based on the GTEx Portal (GTEx Portal V8 Release, https://www.gtexportal.org). Normalized effect size (NES) is defined as the slope of the linear regression and is calculated as the effect of the allele rs1046896-C relative to rs1046896-T. eQTL p-value was calculated by adaptative permutation tests and adjusted using the false discovery rate. (C) Expression of FN3KRP for each individual as a function of age.

The Longevity Association of rs1063192 (CDKN2B-AS1) Replicates With Borderline Significance in an Independent Cohort

We aimed for replication of the association results in 2 independent cohorts. Sample sizes of the Danish and French cohorts (Danish: 1002 LLI, 738 younger controls; French: 1264 LLI, 1830 younger controls) limited the replication approach to the common variants rs1046896 (FN3KRP), rs1063192 (CDKN2B-AS1), and rs1319846 (TMEM131L). The association of rs1063192 (CDKN2B-AS1) reached borderline significance in the French (p = .056 after Bonferroni correction (3 tests), odds ratio [OR] = 1.14; Table 2), but not in the Danish (p = 1.00 after Bonferroni correction (3 tests), OR = 1.04; Table 2). The longevity associations of the other SNVs, including rs1046896 in FN3KRP, did not replicate in the French or Danish (Table 2). Association tests performed with the centenarian subsets of all 3 cohorts yielded similar results (data not shown).

Table 2.

Single-Variant Replication and Meta-analysis Statistics for Candidate SNVs in the French and Danish Populations

					Single-Variant Association Test			Meta-analysis
SNV (Gene)	Sample	MAFa		MA	ORb [95% CI]c	p d	pcorr d	p(R) e	OR(R)f	Q g	I h
		LLI	C
rs1046896 (FN3KRP)	German	0.276	0.324	T	0.78 [0.71–0.86]	7.40E−07		5.54E−01	0.94	1.00E−04	89.35
	French	0.362	0.347		1.06 [0.95–1.19]	2.46E−01	7.38E−01
	Danish	0.278	0.288		1.00 [0.85–1.17]	9.94E−01	1.00
rs1063192 (CDKN2B-AS1)	German	0.482	0.439	G	1.21 [1.11–1.32]	2.99E−05		1.74E−03	1.14	1.85E−01	40.72
	French	0.416	0.385		1.14 [1.02–1.26]	1.88E−02	5.62E−02
	Danish	0.485	0.476		1.04 [0.90–1.20]	6.36E−01	1.00
rs13119846 (TMEM131L)	German	0.486	0.438	C	1.22 [1.11–1.33]	1.73E−05		4.10E−01	1.07	1.20E−03	85.22
	French	0.452	0.464		0.93 [0.84–1.04]	2.31E−01	6.92E−01
	Danish	0.456	0.444		1.08 [0.94–1.24]	3.05E−01	9.15E−01

Notes: Listed are rs-numbers, annotated gene name, chromosome, allele frequencies in cases and controls the minor allele, odds ratios with 95% confidence intervals and allelic p-values for each study population. The effective size of the German population was 1248 LLI and 6762 younger controls; for the French 1270 LLI and 1824 younger controls, and for the Danish 1002 LLI and 738 younger controls. CDKN2B-AS1 = cyclin-dependent kinase inhibitor 2B antisense RNA 1; FN3KRP = fructosamine 3 kinase related protein; TMEM131L = transmembrane 131 like.

C = controls; LLI = long-lived individuals.

aMinor allele frequency, MAF; the definition of the minor allele (MA) is based on controls.

bOdds ratio for longevity, OR; based on the MA in controls.

c95% confidence interval, 95% CI; CI for the OR.

dAllelic p-values, calculated from logistic regression; Pcorr, corrected P-value using Bonferroni (corrected for 3 tests in the French and Danish study populations).

eRandom-effects meta-analysis p-value; fRandom-effects OR estimate; gp-value for Cochrane’s Q statistic, Q; hI2 heterogeneity index (0–100).

In a meta-analysis, we observed large inconsistency (I2) of the genetic effects across the 3 studies (I2 > 75%) for rs1046896 (FN3KRP) and rs1319846 (TMEM131L), but moderate I2 (25% < I2 <75%) for rs1063192 (CDKN2B-AS1) (Table 2). The between-population heterogeneity of the genetic effects for rs1063192 (CDKN2B-AS1) was estimated as 40.72% (Table 2). To account for potential heterogeneity biases, we used the random effects summary OR (OR(R)). The strongest evidence for an association with longevity was observed for rs1063192 (CDKN2B-AS1) (OR(R) = 1.14, p = .00174; Table 2).

The different ancestry of the individuals may have contributed to the heterogeneity of the results from the 3 cohorts. While the heterogeneity within the German sample was low and corrected for by PCA (see Methods section), we did not have equivalent genetic information on the French and Danish cohorts to correct for any influence of genetic heterogeneity.

Gene-Based Analysis Reveals PGP as a Potential New Longevity Locus and Strengthens the FN3KRP Association

In addition to the single-variant association approach, we assessed the cumulative effects of common and rare variants within one genomic region. In the gene-based association test, 13 genes (apart from TOMM40/APOE/APOC1) were identified with an enriched burden of rare and common variants (p < 1 × 10E−04). However, only PGP survived Bonferroni correction in both the SKATO and the burden test (Table 3). Of the 13 genes, 4 (FN3KRP, GRN, SKOR1, and SPZ1) had already been observed in the single-variant association analysis (Table 3). An association analysis using the centenarian subset only, yielded significant associations (p < 3 × 10E−06) for 7 genes (apart from APOE) (Supplementary Table 5). Five genes (TMEM14A, IRAK1BP1, ACPP, PLXNB1, and GAR1) were identified in the centenarian subpopulation only, whereas FN3KRP and HMHA1 overlapped with the results of the gene-based analysis using the entire German study population. Epistatic interactions among SNV candidates (Table 3) were assessed through MDR analysis. This analysis confirmed that rs1046896 in FN3KRP had the largest univariate effect (p < .01) and that the CC genotype represented the favored allele combination for longevity. For the 2-locus interaction model, we observed a significant SNV-SNV interaction between rs1046896-FN3KRP and rs2074442-HMHA1, a gene associated with the immune response. For the 3-locus interaction model, we observed an interaction between the gene rs2074442-HMHA1 and the 2 transcription regulators rs2943549-HNF4G and rs4097-SET9 (Supplementary Figure 4).

Table 3.

Association Statistics for the 16 Longevity-Associated Genes Identified by the Gene-Based Association Approach in the German Study Population

				SNVs
Gene	Chr	P_skato	P_burden	All	Tested	Rare	Common	NCBI rs identification number
APOE	19	3.25E−15	3.25E−15	2	2	1	1	rs769449, rs769452
APOC1	19	2.59E−11	1.77E−02	3	3	0	3	rs439401, rs445925, rs4420638
TOMM40	19	1.35E−06	4.34E−02	3	3	1	2	rs157580, rs2075650, rs142412517
PGP	16	2.50E−06	8.90E−07	3	3	2	1	rs200526199, rs116977380, rs200615324
OTOL1	3	2.56E−04	7.11E−06	5	5	3	2	rs199791179, rs149127996, rs199985412, rs3921595, rs202021352
FN3KRP	17	9.19E−06	3.75E−02	5	5	3	2	rs138953335, rs61743692, rs144986629, rs142718764, rs1046896
SETD9	5	1.19E−05	6.95E−05	6	6	4	2	rs2257505, rs149334074, rs40497, rs141692637, rs150526244, rs146260337
RPS6KB1	17	1.56E−05	1.56E−05	2	2	1	1	rs201316437, rs1051424
GRN	17	1.59E−05	8.28E−03	4	4	4	0	rs63750723, rs63750043, rs63750541, rs63750412
PSG7	19	2.64E−05	4.19E−04	2	2	2	0	rs199532805, rs112354282
SKOR1	15	3.38E−05	9.20E−05	2	2	2	0	rs200956599, rs143419968
HNF4G	8	3.95E−05	2.17E−02	4	4	3	1	rs2943549, rs201625743, rs138897994, rs148532560
ASB17	1	4.42E−05	4.22E−05	3	3	1	2	rs149522654, rs11161887, rs3795251
SPZ1	5	8.09E−05	2.26E−04	5	5	2	3	rs1862136, rs139471643, rs184214819, rs200249535, rs35337118
BFSP1	20	1.25E−03	6.00E−05	6	6	5	1	rs145703098, rs140116733, rs6080719, rs147718368, rs143865632, rs142092768
HMHA1	19	1.54E−03	3.49E−05	8	8	5	3	rs1801284, rs2074442, rs36084354, rs142614852, rs150294461, rs139988914, rs61734935, rs139251906

Notes: Chromosome (Chr); P_skato and P_burden indicate the p-value calculated using SKAT-CommonRare function or the SKAT-O method, respectively, from R-package SKAT. All: number of SNVs genotyped for each gene set; Tested: all rare and common variants tested; Rare: number of rare variants tested. Common: number of common variants tested. NCBI: The National Center for Biotechnology Information. APOE = apolipoprotein E; APOC1 = apolipoprotein C1; ASB17 = ankyrin repeat and SOCS box containing 17; BFSP1 = beaded filament structural protein 1; FN3KRP = fructosamine 3 kinase-related protein; GRN = granulin; HMHA1 = Rho GTPase activating protein 45; HNF4G = hepatocyte nuclear factor 4 gamma; OTOL1 = otolin 1; PGP = phosphoglycolate phosphatase; PSG7 = pregnancy specific beta-1-glycoprotein 7; RPS6KB1 = ribosomal protein S6 kinase B1; SETD9 = SET domain containing 9; SKOR1 = SKI family transcriptional corepressor 1; SPZ1 = spermatogenic leucine.

Discussion

In our German longevity cohort comprising more than 1200 LLI (incl. ~600 centenarians), we screened 62 488 common and rare exonic SNVs for an association with longevity. First, we found rs1046896 in the 3’UTR region of FN3KRP to be associated with exome-wide significance. Second, we confirmed the CDKN2B-AS1 region as a longevity locus and third, we identified PGP as a new candidate gene.

Here, we suggest rs1046896 in FN3KRP as a novel longevity variant though the association only reached exome-wide significance in Germans. Single-association analysis (Figure 1B) revealed that rs1046896 was surrounded by several other signals that were in high LD with rs1046896 (r2 > 0.8) and that displayed significant cis-eQTL associations (Supplementary Table 4). Despite the fact that rs1046896 exhibited the strongest effect size, the gene-based analysis (Table 3, Supplementary Table 5) showed that the FN3KRP association was likely driven by more than just one SNV (5 SNVs in FN3KRP led to the significant association in the gene-based analysis). Additionally, the meta-analysis exposed a very high between-population heterogeneity for rs1046896 (Table 2), suggesting that the SNV could be in different LD with the causal polymorphism in the 3 populations studied. Therefore, our findings indicate that rs1046896 was not the only or main causative variant. This hypothesis could explain the heterogeneity with our replication results. Nevertheless, the longevity association with FN3KRP needs to be confirmed in larger or additional samples. An in-depth fine-mapping of the FN3KRP region would be the next necessary step to identify the true causative variant/s. Furthermore, we cannot exclude that other mechanisms (eg, epigenetics) may be involved in the regulation of FN3KRP expression. For instance, our observation suggests that a high FN3KRP expression seems to be linked to human longevity independent of rs1046896.

FN3KRP belongs to a gene family with an important role in the reversal of the nonenzymatic glycation of proteins (36). Glycation adversely affects protein functioning which leads, among other effects, to arterial stiffening (37) and to an accumulation of damaged proteins. This in turn promotes aging (36) and the development of age-related diseases (38). Both FN3KRP and fructosamine-3-kinase (FN3K), which shows 65% sequence similarity with FN3KRP, seem to protect proteins from nonenzymatic glycation and stop the formation of certain advanced glycation end products (36). Strikingly, rs1046896-T, which was depleted in the LLI in our study (ORcond. = 0.77; ORcond. centenarians = 0.70), was identified as a risk locus for glycated hemoglobin (HbA1c), a critical nonenzymatic glycation product used to monitor and diagnose diabetes (39). Our findings, together with the substantial role of FN3KRP in cell maintenance and viability, support it as a promising candidate that may facilitate healthy aging and longevity.

In the German as well as the French samples, we observed a statistically significant association for rs1063192 in the CDKN2B-AS1 region. This particular SNV was previously suggested as part of a signature of exceptional longevity (11). Other CDKN2B-AS1 SNVs, but not rs1063192, were described in (meta-)GWAS to be associated with longevity in various populations of European ancestry (12,13,40). Thus, our results support the relevance of this region, and maybe of rs1063192 (or one of its LD SNVs), for the phenotype. The G-allele of rs1063192 was enriched in our LLI relative to the controls (Table 2). Interestingly, this allele has been reported to be a protective variant in glaucoma, a classical age-related disease (41) that is also characterized by an increased burden of advanced glycation end products (38). Remarkably, we could not detect the rs1063192 association in the Danish sample. The considerably smaller sample size (Danish: 1740 individuals; French: 3094 individuals; German 8189 individuals) as well as the older age of the Danish controls (<80.2 years; mean age: 66.3 years) compared with the German (<60 years) or French (<62 years) controls might have contributed to the failed replication. The lack of standardized criteria for the selection of suitable controls, together with varying phenotype definitions of longevity, has been suggested as one cause for discrepant association results (11,40,42).

In the single-variant analysis, we found exome-wide significant associations only for the 3 variants in the TOMM40/APOE/APOC1 region (rs2075650, rs4420638, rs769449) and rs1046896 (FN3KRP). Apart from these, 7 additional SNVs showed a suggestive association with longevity (p < 2 × 10E−04, Table 1). Suggestively significant variants may still reflect true association signals with biological relevance, for example, it was shown that 8 of the 10 top not genome-wide significant CHARGE SNVs from a meta-analysis of GWA studies of longevity (25) corresponded to mouse life-span QTL (43). Of the 11 SNVs, for which we report a significant or suggestive longevity association, 4 are rare variants (MAF < 0.01, Table 1). In a recent study, in which the exomes of 100 LLI were sequenced, no rare protein-altering SNVs were observed to be enriched in the genomes of the LLI compared to younger controls (44). However, our results suggest that rare coding variants may very well be drivers of longevity.

Studies on other complex traits have shown that gene-based tests, in which the effects of all common and rare SNVs within a gene locus are considered jointly, can be more biologically relevant than single-variant association approaches. When we examined the cumulative effect of common and rare variants within each gene, we discovered PGP (cumulative effect of 3 variants: 2 rare and 1 common) as another novel longevity-associated locus (p = 8.90 × 10E−07). Through its function as a glycerol-3-phosphate (Gro3P) phosphatase, PGP controls the levels of Gro3P that is an important metabolite formed during glycolysis (45). The availability of Gro3P is crucial for the regulation of both glucose and fat metabolism and eventually determines the generation of signaling and regulatory molecules which further affect many biological processes, for example, insulin secretion and sensitivity, inflammation, fat synthesis and storage, and (cancer) cell proliferation (46). Therefore, PGP may aid in the detoxification of excess nutrient/fuel supplies, thereby preventing metabolic stress and eventually facilitating longevity.

In the gene-based tests, apart from PGP and TOMM40/APOE/APOC1, we identified 12 additional potential longevity genes (p < 1 × 10E−04), including otolin 1 (OTOL1, p-value < 1.61 × 10E−06 in the burden test). One intronic SNV in OTOL1, rs1425609, has already been reported in the context of longevity (47). The fact that the gene-based test also yielded a longevity association for OTOL1, which was first identified in a single-variant analysis (47), not only validates the gene-based approach, but also points towards the added value of analyzing both common and rare SNVs jointly for explaining the missing longevity heritability. Next to OTOL1, the 16 genes included the ribosomal protein S6 kinase B1 (RPS6KB1), a downstream effector of the nutrient-responsive mTOR (mechanistic target of rapamycin kinase) whose inhibition was shown to promote longevity in yeast, worms, and flies (48). Furthermore, a subset of 5 genes with suggestive statistical significance was involved in cell proliferation and cell growth (Supplementary Table 5). Interestingly, the gene FN3KRP reached a p-value of 3.08 × 10E−06 and survived multiple testing correction in the centenarian subset (Supplementary Table 5). This supports the FN3KRP-longevity association as a true-positive signal.

An accumulation of common and rare variants has already been observed in genes associated with other complex traits like type 2 diabetes (49) and Alzheimer’s disease (50). With respect to extreme aging, the genes LYST, MDN1, and RBMXL1 were found to harbor an increased burden of rare coding variants in LLI versus younger controls; however, with nominal significance only (44). To the best of our knowledge, the joint effect of common and rare variants on human longevity has not been investigated yet in a cohort of comparable size to ours. Recently, sets of variants with low frequency (MAF < 0.05) were analyzed for an association with longevity in 530 East Asian nonagenarians and centenarians. More than 100 genes reached nominal significance in that study. However, exome-wide significance was not met by any of the genes (17) and we could not identify any overlap between their and our results.

Conclusion

Our study contributes to the genetic framework of longevity with 2 new potential candidate genes, FN3KRP and PGP, which were identified by single-variant and gene-based analyses, respectively. The 2 genes likely influence longevity through their role in metabolic processes, that is, the reverse glycation of proteins (FN3KRP) and control of Gro3P levels (PGP). However, with respect to FN3KRP, the variant that we report here (rs1046896) is unlikely to be the main causative variant, considering the high between-population heterogeneity values in the meta-analysis of the German, Danish, and French populations. Future fine-mapping studies are warranted to identify the true functional variant(s). With the combination of analysis methods and, in particular, the investigation of cumulative effects of common and rare variants within one genetic region, we are a small step closer to accounting for the missing heritability of human longevity.

Funding

G.G.T. was funded by the Cluster of Excellence “Inflammation at Interfaces” (EXC 306). This work was supported by the German Federal Ministry of Education and Research (BMBF) within the framework of the e:Med research and funding concept (SysInflame grant 01ZX1606A). The biobank PopGen and the PopGen 2.0 network are financed by the BMBF (grant 01EY1103). The CEPH centenarian cohort, maintained at the CEPH Biobank, Paris, France [BIORESOURCES], was supported by the “Ministère de l’Enseignement supérieur et de la Recherche.” The SU.VI.MAX cohort was funded by the French Institut National de la Santé et de la Recherche Médicale, the Institut National de la Recherche Agronomique, the Université Paris 13. The genotyping of French samples was supported by the Commissariat à l’Energie Atomique-Centre National de Recherche en Génomique Humaine. The Danish replication study was supported by The National Program for Research Infrastructure 2007 (grant no. 09-063256), the Danish Agency for Science Technology and Innovation, the Velux Foundation, and the U.S. National Institute of Health (P01 AG08761). Genotyping of the Danish controls was conducted by the SNP&SEQ Technology Platform, Science for Life Laboratory, Uppsala, Sweden (http://snpseq.medsci.uu.se/genotyping/snp-services/) and supported by NIH R01 AG037985 (Pedersen). The KORA research platform (KORA, Cooperative Research in the Region of Augsburg) was initiated and financed by the Helmholtz Zentrum München – German Research Center for Environmental Health which is funded by the BMBF and the State of Bavaria. Furthermore, KORA research was supported within the Munich Center of Health Sciences (MC Health), Ludwig-Maximilians-Universität, as part of LMU innovativ.

Author Contributions

G.G.T., J.D., F.F., A.F., and A.N. designed research; J.D., A.F., and A.N. supervised the project; W.L., A.F., F.F., D.E., S.S., K.S., M.M.-N., A.P., H.N., and P.H. were involved in recruitment of German study subjects and assembling of phenotypic data; F.F., D.E., and A.N. organized chip genotyping of German long-lived individuals; M.N., H.B., S.C., P.G., L.C., J.-F.D., and K.C. performed replication experiments; G.G.T. analyzed the data and together with A.C. performed the statistical analysis; G.G.T., J.D., and A.N. interpreted the data and wrote the manuscript; all authors performed critical revision and approved the final version of the manuscript.

Conflict of Interest

None declared.

Click here for additional data file.