Strong impact of natural-selection–free heterogeneity in genetics of age-related phenotypes

Study overview

We performed two-stage univariate and pleiotropic meta-analyses of genetic predisposition to 20 weakly-to-modestly correlated age-related phenotypes (Supplementary Figure 1). The data were drawn from five longitudinal studies for 33,431 men and women combined, who identified themselves as of Caucasian ancestry. In stage 1, we first performed simplified univariate GWAS of each phenotype in each cohort separately using plink software [21]. Then, these results were prioritized and top 1,000 promising SNPs were selected for more comprehensive analysis leveraging information on repeated measurements for quantitative markers and timing of the risk outcomes (Tables 1 and Supplementary Table 1). In stage 2, we performed univariate meta-analysis to combine statistics for these 1,000 selected SNPs across cohorts and pleiotropic meta-analysis to combine such statistics across phenotypes. These analyses addressed the natural-selection–free genetic heterogeneity in genetic predisposition to age-related phenotypes (see the Introduction) by performing five meta-tests (Fig. 1) (details in Methods).

Figure 1. Scheme of univariate and pleiotropic meta-analyses in stage 2. (A) Statistics from stage 1 univariate GWAS of 20 phenotypes in five cohorts denoted PiCj, $\mathrm{i}\in \left(\stackrel{-}{\mathrm{1,20}}\right)$ and $\mathrm{j}\in \left(\stackrel{-}{\mathrm{1,5}}\right)$. (B) Univariate statistics from the meta-analysis across cohorts using the fixed-effects meta-test (PiM) and Fisher test (PiFc). (C) Statistics from the pleiotropic meta-analysis across phenotypes in cohort j for: (i) omnibus tests with correlation matrix for phenotypes ${\mathrm{\Sigma }}_{\mathrm{j}}^{\mathrm{P}}$ (OpCj) and effect statistics ${\mathrm{\Sigma }}_{\mathrm{j}}^{\mathrm{B}}$ (ObCj) and (ii) Fisher test (FpCj). (D) Meta-statistics from Fisher test across cohorts for the results in (C) (OpFc, ObFc, and FpFc). (E) Meta-statistics across phenotypes for the results in (B) from (i) meta-test across cohorts and omnibus test across phenotypes with correlation matrix for phenotypes ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{P}}$ (MOp), (ii) meta-test across cohorts and omnibus test across phenotypes with correlation matrix for effect statistics ${\mathrm{\Sigma }}_{\mathrm{m}}^{\mathrm{B}}$ (MOb), (iii) meta-test across cohorts and Fisher test across phenotypes (MFp), and (iv) Fisher test across cohorts and Fisher test across phenotypes (FcFp). Pathway 1: meta-analysis combining statistics across cohorts (pathway 1a) and pleiotropic meta-analysis across phenotypes (pathway 1b). Pathway 2: meta-analysis combining statistics across phenotypes (pathway 2a) and cohorts (pathway 2b).

Univariate meta-analysis

Using the Fisher test and a fixed-effect meta-test (Fig. 1, pathway 1a), we identified two novel (see Methods) non-proxy SNPs (defined as LD r²<70%, actual LD for the reported SNPs in the same locus ranged r²=1% to 64%, Supplementary Table 2) associated at GW level with HC (rs6745983, β_meta=-0.077, p_meta=4.2×10^-8) and BG (rs10885409, β_meta=0.672, p_meta=1.9×10^-8) (Supplementary Table 3). We also replicated (p<p_GW) 24 associations, primarily with lipids (21 associations), for 16 SNPs in or nearby (within ±100 Kb flanking region for the index SNP) 11 protein coding genes (Supplementary Table 3). Seven associations were replicated in the LPL-gene locus. P-values were smaller for the meta-test than for the Fisher test (p_meta<p_Fisher) for all associations except for the associations of rs780094 with BG (p_meta=2.7×10^-8 vs. p_Fisher =2.4×10^-11) and rs261332 with TC (p_meta=5.1×10^-8 vs. p_Fisher =2.7×10^-9). The Fisher test benefited the meta-analysis because it disregarded strong heterogeneity in the effect sizes across cohorts measured by the heterogeneity coefficient I², i.e., I²=87.8% for rs780094 and I²=82.2% for rs261332, as supported by significant inverse correlation of the ratio log₁₀(p_meta)/log₁₀(p_Fisher) with I², r_Pearson=-0.634, p=5.0×10^-4.

Significance of the association of rs2155216 with TG in our modest sample of 33,431 individuals (p_meta=9.4×10^-14) was about the same as previously reported (p_grasp=1.6×10^-14) in a sample of 140,059 individuals [22]. These p-values indicate improvement of the efficiency of the analysis that can be quantified by the ratio of the log₁₀(p-value) to the sample size [23] yielding 4-fold (=log₁₀(9.4×10^-14)×140,059/log₁₀(1.6×10^-14)/33,431) larger efficiency in our study than in [22]. This improvement was partly achieved by leveraging longitudinal information on repeated measurements that resulted in smaller standard errors (Supplementary Figure 2).

Pleiotropic meta-analysis

We used four tests in pathway 1 and three tests in pathway 2 (Fig. 1) to examine pleiotropy in the domains of 12 quantitative markers, 8 diseases and death, and all 20 phenotypes (Table 1). Pleiotropic associations were defined as at least one test showing GW significance (see Methods).

This analysis identified 115 novel pleiotropic SNPs in or nearby (within ±100 Kb flanking region for the index SNP) 84 protein coding genes (Table 2) and replicated nine SNPs (Supplementary Table 4) with p_pleio<p_GW by combining associations with multiple phenotypes that individually did not attain GW significance in univariate meta-analysis (p_uni>p_GW) (Supplementary Table 5). SNPs were identified as novel if they did not attain GW significance in our pleiotropic meta-analysis of the results collected in GRASP for the index SNPs or for flanking SNPs, if no results for 20 selected phenotypes (Figure 1) were reported in GRASP for the index SNPs (see Methods). Specifically, of 115 SNPs, 31 were not reported in GRASP for the selected 20 phenotypes. For them, we identified multiple SNPs within ±1Mb flanking region in GRASP. For the flanking SNPs, which were in the strongest LD to the index SNPs, no GW significant pleiotropic associations were identified (Table 3). Nine flanking SNPs showed GW significant pleiotropic effects, which were independent of the pleiotropic effects for the index SNPs. Of 115 novel SNPs, 29 SNPs attained smaller p-value in pathway 1 and the remaining 86 SNPs in pathway 2.

All SNPs attained GW significant associations either in the domain of 12 markers (76 of 115 novel SNPs and seven of nine replicated SNPs) or all 20 phenotypes (39 novel and two replicated SNPs). No associations with p<p_GW were identified in the disease/death domain (Supplementary Table 6). The 115 novel and nine replicated SNPs were associated with up to 10 phenotypes at p<0.05 (Table 2 and Supplementary Table 4). Figure 2 shows 39 novel pleiotropic SNPs attained GW significance in the domain of 20 phenotypes.

Figure 2. Heat map of phenotype-specific associations for selected pleiotropic SNPs. Data are for 39 novel SNPs with pleiotropic associations in the domains of all 20 phenotypes (20P) from Table 2. Numbers in parentheses are SNP IDs in Table 2. Phenotypes are defined in Table 1. FF, Op, and Ob denote pleiotropic meta-tests either from pathway 1 or 2 based on Fisher test, omnibus test with correlation matrix for phenotypes, and omnibus test with correlation matrix for effect statistics (Fig. 1), respectively. MFp denotes pleiotropic meta-tests from pathway 1 (Fig. 1). Colors code -log₁₀(p-value) trimmed at GW level -log₁₀(5×10^-8)=7.3 for better resolution.

Because p_pleio<p_GW<p_uni for pleiotropic SNPs, this inequality automatically validates pleiotropy for 29 SNPs in pathway 1, as it implies that pleiotropic statistics improved to attain p_pleio<p_GW by pooling contributions from multiple phenotypes with p_uni>p_GW (Supplementary Table 5). Likewise, GW significant pleiotropic associations in meta-analysis in pathway 2 were automatically validated for 37 SNPs with p_pleio>p_GW in each individual cohort (i.e., in pathway 2a) as it implies that pleiotropic statistics improved to attain p_pleio<p_GW in meta-analysis of individual cohorts. For the remaining 49 SNPs in pathway 2, we observed p_pleio<p_GW in at least one cohort. For five of these 49 SNPs, pleiotropic associations did not attain Bonferroni adjusted significance level, p<0.05/(N_cohorts – 1), in at least one additional cohort, whereas they attained that level for the other 44 SNPs (Supplementary Table 7).

Antagonistic genetic heterogeneity in pleiotropic meta-analysis

Let us assume that an allele is associated with a higher risk of a phenotype P₁ and that P₁ is directly correlated with a phenotype P₂, e.g. individuals with larger values of P₁ tend to have larger values of P₂. An implicit expectation in medical genetics for partly correlated phenotypes, especially when they are etiologically related, is that the same allele will be also associated with a higher risk of a phenotype P₂. In the framework of the undefined role of evolution in establishing molecular mechanisms of age-related phenotypes (see Introduction), this logic may or may not hold because given allele may not be evolutionary selected in favor or against these phenotypes. A simple approach to examine this logic is to compare the results from different tests used in our pleiotropic meta-analysis. Indeed, the Fisher test in pathways 1b and 2a (Fig. 1) combined p-values across phenotypes assuming that they were from independent associations whereas the two omnibus tests adjusted the pleiotropic meta-statistics for correlation among phenotypes (${\mathbf{\Sigma }}^{\mathit{P}}$-based omnibus test) and the effect statistics (${\mathbf{\Sigma }}^{\mathit{B}}$-based omnibus test) (see Methods). Accordingly, the differences in p-values from the Fisher and omnibus tests reflect the impacts of heterogeneity and/or correlation in genetic associations. Below, we characterize these impacts. We used an ad-hoc cut-off for the difference in p-values between these tests of ≥2 orders of magnitude to characterize a strong impact.

Given a common assumption that the Fisher pleiotropic meta-test provides inflated p-values for correlated phenotypes, correlation-adjusted estimates in the omnibus tests are expected p_Omnibus>p_Fisher. Contrary to this expectation, for most novel SNPs, 93.9% (108 of 115 SNPs) and all nine replicated SNPs, we observe an opposite inequality, i.e., p_Omnibus<p_Fisher by ≥2 orders of magnitude (Table 2 and Supplementary Table 4A). An unconventional set of 108 novel SNPs includes 27 SNPs from HP group in pathway 1, and 54, 5, and 22 SNPs from HP, HB, and HP/HB groups, respectively, in pathway 2 (Table 2). The HP group (27+54=81 SNPs) was characterized by attaining GW significance after the ${\mathbf{\Sigma }}^{\mathit{P}}$-based omnibus tests (i.e., OpFc or MOp, Fig. 1) and substantially larger (≥2 orders of magnitude) p-values in the Fisher pleiotropic meta-tests, i.e., p_MOp<p_MFp for pathway 1 and p_OpFc<p_FpFc for pathway 2. The HB group (5 SNPs) resembled the HP group except for having p_ObFc<p_FpFc. The HP/HB group (22 SNPs) was characterized by attaining GW significance after the ${\mathbf{\Sigma }}^{\mathit{P}}$- and ${\mathbf{\Sigma }}^{\mathit{B}}$-based omnibus tests (p_OpFc, p_ObFc<p_FpFc) and substantially larger (≥2 orders of magnitude) p-values in the Fisher pleiotropic meta-tests (p_FpFc).

These results show that majority of novel SNPs with pleiotropic associations, 93.9%, is characterized by a strong impact of unconventional (natural-selection–free) genetic heterogeneity rather than commonly expected correlation. This form of the natural-selection–free genetic heterogeneity, called antagonistic heterogeneity, is schematically illustrated in Fig. 3 for two partly correlated phenotypes P1 and P2. This heterogeneity results in orthogonal directions of the bivariate vector of genetic effects β₁ and β₂ for phenotypes P1 and P2 and the vector of the correlation of these phenotypes and, accordingly, in opposite directions of vectors β₂ and P₂. As seen from Fig. 3, the antagonistic heterogeneity implies antagonistic directions of genetic effects for directly correlated phenotypes (see Methods).

Figure 3. Schematic illustration of antagonistic genetic heterogeneity in the associations with two partly correlated age-related phenotypes P1 and P2. Small dots represent a sample of carriers of an effect allele A (red color) and those who do not carry this allele (no A; green color). Ellipse shows correlation of P1 and P2 in this sample (r=0.6). Red color denotes vector of correlation of P1 and P2 (thick diagonal vector) and its projections on phenotypes, i.e. P₁ (horizontal) and P₂ (vertical). Black vectors β₁ and β₂ denote the effects in the associations of allele A with P1 and P2. Sum of β₁ and β₂ represents bivariate vector (thick line) of the effects.

Pathway 2 provides a natural opportunity to validate the antagonistic heterogeneity in different cohorts. Our analysis shows that this heterogeneity, characterized by p_OpFc, p_ObFc<p_FpFc with a difference of ≥0.2 orders of magnitude for associations with suggestive significance (p_OpFc, p_ObFc <0.1 or p_FpFc<0.1), was replicated for 77 of 81 novel pleiotropic SNPs in two or more cohorts (Supplementary Table 8).

Lastly, the M group (7 SNPs) included 2 SNPs from pathway 1 and 5 SNPs from pathway 2. It was characterized by attaining GW significance in several tests and/or minor (<2 orders of magnitude) differences between p-values in OpFc and FpFc tests (pathway 2) and MOp and MFp (pathway 1).

Bioinformatics analysis

We performed biological pathway (IPA, www.qiagenbioinformatics.com) and gene ontology (GO) biological processes (BPs) (DAVID [24]) enrichment analysis for 96 protein coding genes mapped from 108 pleiotropic SNPs (Table 2). We excluded genes for five non-validated SNPs (Table 2, ID 1, 24, 25, 49, and 95) and for two SNPs (Table 2, ID 53, 55) from the Major Histocompatibility Complex (MHC). The GO analysis (by DAVID 6.8 GO category “GO_direct”) identified 10 BPs with enrichment for genes at p<0.05 (Fisher’s exact test) (Supplementary Table 9). Two specific terms, adenylate cyclase-activating G-protein coupled receptor signaling pathway and activation of adenylate cyclase activity, were enriched at p<10^-4 and p<10^-3, respectively. The IPA analysis identified 18 pathways (Supplementary Table 10) enriched for genes at p<0.05 (Fisher’s exact test). The strongest enrichment was observed for G-protein coupled receptor signaling (p<10^-4), GPCR-mediated integration of enteroendocrine signaling exemplified by an L cell and cAMP-mediated signaling (p<10^-3) (Figure 4) that is consistent with enrichment of GO BPs.

Figure 4. Enrichment of pathways in the Ingenuity Pathway Analysis (IPA) bioinformatics tool. Blue bars (upper x-axis) show –log₁₀(p-value) for the top five IPA pathways. Proportion of genes from the identified sets to those in the IPA pathways (orange symbols and line) is shown on the lower x-axis. Numerical estimates are given Supplementary Table 10.

Study cohorts

Data were drawn from the Atherosclerosis Risk in Communities (ARIC) study [40,41], the Cardiovascular Health Study (CHS) [42], the Multi-Ethnic Study of Atherosclerosis (MESA) [43], the Framingham Heart Study (FHS) [44-46], and the Health and Retirement Study (HRS) [47] for individuals who identified themselves as of Caucasian ancestry.

Phenotypes

The analyses focused on 20 phenotypes (12 quantitative markers, 7 diseases, and death) listed in Table 1. Given longitudinal design of the studied cohorts, the analyses leveraged repeated measurements of quantitative markers during follow up and timing information on onsets of diseases or death (Supplementary Table 1). These not strongly correlated phenotypes (Supplementary Figure 1) were available in majority of the selected cohorts. All studies except HRS collected information on diseases and death in population samples during follow-up. The HRS assessed information on date of death from the National Death Index Cause of Death file. Age at onset of diseases in HRS was assessed from the linked Medicare service use files, which include enrollment information and the diagnoses made (International Classification of Disease-revision 9, Clinical Modification) during episodes of care paid for by the Medicare system.

Genotyping

SNPs were available from Affymetrix 1M chip in ARIC and MESA, Illumina CVDSNP55v1_A chip (~50K SNPs) in FHS and CHS cohorts, Illumina HumanCNV370v1 chip (370K SNPs) in CHS, Affymetrix 500K in FHS, and Illumina HumanOmni 2.5 Quad chip (~2.5M SNPs) in HRS. SNPs were included in the analyses after quality control in each study (call rate>95%, Hardy-Weinberg disequilibrium p<10^-6, minor allele frequency [MAF] >2%). In case of marginally smaller MAF for prioritized SNPs in a specific cohort compared to the MAF cut off, these SNPs were used regardless of the MAF cut off. To facilitate cross-platform comparisons, we selected directly genotyped (index) SNPs using all available arrays for each study. Non-genotyped SNPs were imputed (IMPUTE2 [48],) according to the 1000 Genomes Project Phase I integrated variant set release (SHAPEIT2) in the NCBI build 37 (hg19) coordinate. Only SNPs with high imputation quality (info>0.8) were retained for the analyses.

Mapping to genes

SNPs were mapped to genes using variant effect predictor from Ensembl and NCBI SNP database (assembly GRCh38.p7). If an index SNP was not within protein coding gene, the closest gene(s) within ±100 Kb flanking region was (were) assigned. Multiple genes were selected if they were at about the same distance up- and downstream from the index SNP or the index SNP was within the region of overlapping genes.

Selecting SNPs in univariate GWAS in stage 1

GWAS was performed for each phenotype in each cohort separately using plink software [21]. In this analysis, we used quantitative markers measured at baseline, and diseases and death as binary outcomes. An additive genetic model with minor allele as an effect allele was adopted in all analyses throughout this article. We computed pleiotropic p-value by combining p-values for individual markers and binary outcomes for SNPs attained nominal significance only (p<0.05) using Fisher’s test [49]. This computationally efficient but overly-simplified approach was used for prioritizing of promising SNPs. Because this procedure provides inflated p-values, we selected top 1,000 SNPs (with smallest p-values) for all downstream analyses, rather than selecting SNPs based on specific cut off for pleiotropic p-values.

Leveraging longitudinal information in stage 1 analysis

The analysis of the selected 1,000 SNPs were enhanced by leveraging longitudinal information. For quantitative markers, we used all measurements available during follow up of the same individuals. Information on longitudinal measurements has multiple advantages including potential gain in statistical power in the analyses [50]. To correct for repeated-measurements (all cohorts) and familial (FHS) correlations in the analyses of quantitative markers, we mostly used the linear mixed effects model (lme4 package in R [51]). Measurements of BG, BMI, HDL-C, HR, TC, and TG were natural-log-transformed to offset potential bias due to skewness of their frequency distributions. They were multiplied by 100 for better resolution. Measurements of creatinine, DBP, FVC, HC and SBP were used in their natural scale as no significant skewness was observed. Given gamma-like frequency distributions of CRP, a generalized linear mixed model (glmmPQL package in R) with a gamma function and log-link was used. We evaluated the associations for SNPs given the measurements of quantitative markers for individuals of a given age at each examination with available measurements.

Longitudinal information on time to events was implemented in the Cox proportional hazards mixed effects model (coxme package in R). This model addressed familial relatedness. In the FHS, we used both prospective and retrospective onsets. The use of retrospective onsets in a failure-type model is justified by Prentice, Breslow [52]. These analyses provide estimates of the effects in a given population. Time variable in these analyses was the age at onset of an event or at right censoring.

The models were adjusted for: (all studies) age and sex; (ARIC, CHS, and MESA) field center; (HRS) HRS cohorts, and (FHS) whether the DNA samples had been subject to whole-genome amplification [53]. No adjustment for principal components was performed as argued in Ref [54].

Meta-analyses in stage 2

After stage 1, for each SNP we have a table with the association statistics for up to 20 phenotypes in 5 cohorts. These statistics were combined along two possible pathways: (i) first across studies and then across phenotypes and (ii) first across phenotypes and then across studies. To deal with the natural-selection–free heterogeneity in genetic predisposition to age-related phenotypes (see below), we used four tests in pathway 1 and three tests in pathway 2 (Fig. 1).

In pathway 1, univariate (phenotype-specific) meta-analysis combining the stage 1 statistics for the selected SNPs across cohorts (pathway 1a) was performed using the Fisher test and the traditional GWAS fixed-effects meta-test. Pleiotropic meta-analysis (pathway 1b) was performed by combining the univariate meta-statistics for the same SNPs across phenotypes. In pathway 1b, we used the Fisher test and two omnibus tests. The latter were designed to address correlations between the effect statistics and phenotypes. In pathway 2, we performed first pleiotropic meta-analysis by combining the stage 1 univariate statistics across phenotypes in each cohort separately (pathway 2a) using the Fisher test and the two omnibus tests as in pathway 1b. Then, we combined these pleiotropic meta-statistics across cohorts using the Fisher test (pathway 2b) (details on all tests are below).

Fixed-effects meta-test

We adopted a conventional GWAS meta-test using a fixed effects model with inverse-variance weighting (METAL software [55]). This test combines the estimated effect sizes across cohorts for each phenotype. Combining effect sizes is feasible because each phenotype was harmonized to be on the same scale and unit across cohorts. This meta-test could account for the directions of effects in different cohorts, and is more powerful than those tests that combines p-values or Z-scores [56]. In pathway 1a, the weighted average of the effect sizes was calculated as ${\mathrm{\Sigma }}_j\left({\widehat{w}}_j{\widehat{\beta }}_j\right)/{\mathrm{\Sigma }}_j\left({\widehat{w}}_j\right)$ with variance 1/${\mathrm{\Sigma }}_j\left({\widehat{w}}_j\right)$, where ${\widehat{w}}_j$ is the inverse variance of effect size ${\widehat{\beta }}_j$in the cohort $j\in \left(\stackrel{-}{\mathrm{1,5}}\right)$ for a given phenotype and SNP. Wald test was then used to obtain p-value.

Fisher’s test

The Fisher method [49] combines p-values assuming that they are from independent tests. In pathway 1a, this test combines p-values across cohorts. It has power to reject the “null” hypothesis of no pooled effect regardless of the effect sizes and directions in the cohort-specific estimates. Accordingly, the Fisher test can indicate heterogeneity in genetic associations by providing smaller p-value than that from the fixed-effects meta-test. In pathways 1b and 2a, the Fisher test combines p-values across phenotypes. This test is often used for pleiotropic meta-analysis of modestly correlated phenotypes [57]. Because the Fisher method combines p-values from multiple tests, it addresses the problem of multiple testing by increasing the number of degrees of freedom.

Omnibus tests

The statistics from tests of the same SNPs with different phenotypes may or may not be independent. It is then argued that tests penalizing for correlation of such statistics should be used to deflate the Fisher test estimates. A commonly adopted test in this case is an omnibus test [58-60].

For a certain SNP we have an estimated effect size ${\widehat{\beta }}_{ij}$ and its standard error ${\widehat{\sigma }}_{ij}$ for the phenotype $i\in \left({\stackrel{-}{1,K}}_j\right)$ in the cohort $j\in \left(\stackrel{-}{\mathrm{1,5}}\right)$, where $K_j$ is the number of phenotypes in study $j$. The omnibus test statistic is constructed as ${\widehat{\mathit{z}}\text{'}}_j{\mathbf{\Sigma }}_j^{-1}{\widehat{\mathit{z}}}_j$ where ${\widehat{\mathit{z}}}_j={\widehat{\mathit{\beta }}}_j/{\widehat{\mathit{\sigma }}}_j$ is a z-score vector of associations of SNPs with phenotypes and ${\mathbf{\Sigma }}_j$ is the correlation matrix of the z-scores to be estimated [58,59]. Accordingly, this test takes into account the correlation of the effect statistics for different phenotypes. Under the null hypothesis (${\mathit{\beta }}_{\mathit{j}}=0$), the test statistic follows a chi-squared distribution with $K_j$ degrees of freedom

${\widehat{\mathit{z}}\mathrm{\text{'}}}_j{\mathbf{\Sigma }}_j^{-1}{\widehat{\mathit{z}}}_j~{\chi }_{K_j}^2\mathrm{ },$

from which we obtained a combined p-value $p_j$ in the study $j$. Accordingly, this statistic addresses the problem of multiple testing by increasing the number of degrees of freedom.

Correlation matrices $\mathbf{\Sigma }$ were estimated based on the effect statistics of simulated SNPs in association models with phenotypes (denoted $\mathbf{\Sigma }\equiv {\mathbf{\Sigma }}^{\mathit{B}}$). We randomly permuted the dosage data for alleles of a given SNP for 250 times, used each permuted SNP in association models with each of K phenotypes, and estimated ${\mathbf{\Sigma }}^{\mathit{B}}$ based on the 250 effect sizes of permuted SNPs with K phenotypes. Likewise, $\mathbf{\Sigma }$ can be also constructed by evaluating the correlations between phenotypes [61] (denoted $\mathbf{\Sigma }\equiv {\mathbf{\Sigma }}^{\mathit{P}}$). The correlation matrices, ${\mathbf{\Sigma }}^{\mathit{B}}$ and ${\mathbf{\Sigma }}^{\mathit{P}}$, were constructed for pathways 1 and 2 (Fig. 1) separately. Cohort-specific matrices ${\mathbf{\Sigma }}_{\mathit{j}}^{\mathit{B}}$ for pathway 2 were constructed by evaluating correlations of the effect statistics in each cohort. For pathway 1, matrix ${\mathbf{\Sigma }}_{\mathit{m}}^{\mathit{B}}$ was constructed by evaluating correlations of the effect statistics from a fixed effect meta-test of all cohorts for each permutation. The phenotype-based cohort-specific matrices for pathway 2, ${\mathbf{\Sigma }}_j^P$, were evaluated using phenotype measurements in each cohort separately and matrix for pathway 1, ${\mathbf{\Sigma }}_m^P$, was evaluated using phenotype measurements in the combined data from all cohorts. All correlation matrices were constructed using averaged values for quantitative traits measured longitudinally at different visits.

Correlation and genetic heterogeneity

Omnibus tests are commonly used to penalize correlation between the effect statistics or phenotypes in pleiotropic meta-analysis. This can be explicitly illustrated by the test statistics ${\widehat{\mathit{z}}\text{'}}_j{\mathbf{\Sigma }}_j^{-1}{\widehat{\mathit{z}}}_j$ in two-dimensional case,

${\widehat{\mathit{z}}\mathrm{\text{'}}}_j{\mathbf{\Sigma }}_j^{-1}{\widehat{\mathit{z}}}_j=\left[\left({\widehat{z}}_1^2{\mathrm{\Sigma }}_{22}-{{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{21}\right)+\left({\widehat{z}}_2^2{\mathrm{\Sigma }}_{11}-{{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{12}\right)\right]/\mathrm{d}\mathrm{e}\mathrm{t}\left(\mathbf{\Sigma }\right)\mathrm{ }$.(1)

Either the ${\mathbf{\Sigma }}^{\mathit{B}}$ or ${\mathbf{\Sigma }}^{\mathit{P}}$-based omnibus test (see above) penalize pleiotropic meta-statistics, if ${{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{21},{{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{12}>0$.

The undefined role of evolution in establishing molecular mechanisms of age-related phenotypes is the source of genetic heterogeneity, which is driven by the lack, rather than the presence, of the evolutionary pressure in favor or against these phenotypes. The natural-selection–free genetic heterogeneity suggests that mechanisms driving associations of genes with age-related phenotypes and correlation between these phenotypes are, generally, of different origins. In case of antagonistic genetic heterogeneity (see “Antagonistic genetic heterogeneity in pleiotropic meta-analysis” section), omnibus test provides smaller p-values compared to the Fisher test as can be seen from (1) because ${{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{21},{{\widehat{z}}_1{\widehat{z}}_2\mathrm{\Sigma }}_{12}<0$ in this case.

Significance and novelty

The fixed-effect, Fisher, and two omnibus meta-tests used in our analyses have power to identify associations, adjusted and unadjusted for correlation and heterogeneity in genetic predisposition to age-related phenotypes (see above). Fisher and two omnibus tests used in pleiotropic meta-analysis address the problem of multiple testing by increasing the number of degrees of freedom. Accordingly, GW level (p=5×10^-8) attained in either of these tests was used as a cut off for the significance for a given SNP. The difference between p-values from these tests for a given SNP was used to characterize the impact of correlation and heterogeneity on the association.

SNPs were considered as novel if they attained GW significance in: (i) our univariate meta-analysis but were not reported in GRASP catalog [62] at p≤5×10^-8 or (ii) our pleiotropic meta-analysis but not in our pleiotropic analysis of the results collected in GRASP. For univariate meta-analysis, we used evidences for the same (index) SNP in our study and GRASP. For pleiotropic meta-analysis, we selected associations from GRASP for the index SNPs with 20 phenotypes used in our analysis (Table 1). If an index SNP was not available in GRASP for the selected phenotypes, we selected SNPs within ±1Mb flanking region. Then we performed pleiotropic meta-analysis by applying the Fisher method to these GRASP results to present evidence for pleiotropy in prior studies. We used the Fisher statistics p_grasp because the effect sizes were not reported in GRASP. We used a flat p-value, p=0.4, to penalize the Fisher statistics for phenotypes with p-values not reported in GRASP. For flanking SNPs, we reported associations for SNPs having: (i) the strongest LD with the index SNP and (ii) the smallest p_grasp for SNPs within ±1Mb flanking region.

Heterogeneity coefficient

We used METAL software [55] to evaluate the heterogeneity coefficient I². The I² can be interpreted as the percentage of the total variability in a set of effect sizes due to between-sample variability.