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

A conceptual difficulty in genetics of age-related phenotypes that make individuals vulnerable to disease in post-reproductive life is genetic heterogeneity attributed to an undefined role of evolution in establishing their molecular mechanisms. Here, we performed univariate and pleiotropic genome-wide meta-analyses of 20 age-related phenotypes leveraging longitudinal information in a sample of 33,431 individuals and dealing with the natural-selection–free genetic heterogeneity. We identified 142 non-proxy single nucleotide polymorphisms (SNPs) with phenotype-specific (18 SNPs) and pleiotropic (124 SNPs) associations at genome-wide level. Univariate meta-analysis identified two novel (11.1%) and replicated 16 SNPs whereas pleiotropic meta-analysis identified 115 novel (92.7%) and nine replicated SNPs. Pleiotropic associations for most novel (93.9%) and all replicated SNPs were strongly impacted by the natural-selection–free genetic heterogeneity in its unconventional form of antagonistic heterogeneity, implying antagonistic directions of genetic effects for directly correlated phenotypes. Our results show that the common genome-wide approach is well adapted to handle homogeneous univariate associations within Mendelian framework whereas most associations with age-related phenotypes are more complex and well beyond that framework. Dissecting the natural-selection–free genetic heterogeneity is critical for gaining insights into genetics of age-related phenotypes and has substantial and unexplored yet potential for improving efficiency of genome-wide analysis.


AGING
The problem becomes even more challenging for agerelated phenotypes, i.e., phenotypes that make individuals vulnerable to disease in post-reproductive life. A conceptual difficulty in genetics of age-related phenotypes is the natural-selection-free genetic heterogeneity attributed to an undefined role of evolution in establishing their molecular mechanisms [7,8]. This problem is complicated by recent changes in human life expectancy [9] and the fitness landscape [10][11][12][13]. Accordingly, in evolutionary biology, age-related phenotypes are viewed as the results of indirect mechanisms ("side effects") such as co-evolution with fast-evolving pathogens, mismatch with environments, reproductive success at the expense of health, trade-offs that leave every trait suboptimal, defenses and their special costs, etc. [7]. The concept of heritability in the evolutionary framework for age-related phenotypes becomes even more problematic because the estimates of heritability change as environment changes [5,6] and significant heritability does not imply that the same genetic variant carry the same risk in different population groups, even of the same ancestry [11].
The natural-selection-free genetic heterogeneity has not been addressed in mainstream GWAS. This heterogeneity implies that differences in genetic predisposition to age-related phenotypes across different population groups is biologically motivated. Accordingly, different, even antagonistic, effects of the same allele on the same phenotype in different population groups are biologically plausible [14,15]. Another challenge in the evolutionary framework is that genetic variants predisposing to a phenotype may not necessarily predispose to another, even causally related, phenotype [16,17] or such genetic variants can predispose to seemingly unrelated phenotypes [18][19][20].
Here, we examine genetic predisposition to age-related phenotypes following the concept of an undefined role of evolution in establishing their molecular mechanisms. We performed the univariate and pleiotropic GWAS meta-analyses of 20 age-related phenotypes in the sample of 33,431 Caucasians from five longitudinal studies. We identified 142 non-proxy (defined as linkage disequilibrium [LD] r 2 <70%) single nucleotide polymorphisms (SNPs) with genome-wide (GW) significance (p<p GW =5×10 -8 ). They include two novel and 16 previously reported SNPs attained GW significance in univariate meta-analysis of individual phenotypes, and 115 novel and nine previously reported SNPs attained GW significance in pleiotropic metaanalysis. We show that pleiotropic associations for most novel SNPs, 93.9% (108 of 115 SNPs), and for all replicated SNPs were strongly affected by the naturalselection-free genetic heterogeneity in a rarely recognized form of antagonistic heterogeneity, implying antagonistic directions of genetic effects for directly correlated phenotypes.

Study overview
We performed two-stage univariate and pleiotropic meta-analyses of genetic predisposition to 20 weaklyto-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 naturalselection-free genetic heterogeneity in genetic predisposition to age-related phenotypes (see the Introduction) by performing five meta-tests ( Fig. 1) (details in Methods).

AGING
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 10 (p-value) to the sample size [23] yielding 4-fold (=log 10 (9.4×10 -14 )×140,059/log 10 (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.  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.
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 1 and that P 1 is directly correlated with a phenotype P 2 , e.g. individuals with larger values of P 1 tend to have larger values of P 2 . 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 2 . In the framework of the undefined role of evolution in establishing molecular mechanisms of age-related Domain indicates the pleiotropic domain in which minimal GW significant p-values were attained in our pleiotropic meta-analysis examining: 12 quantitative markers (12M), 7 diseases and death (8DD), and all 20 phenotypes (20P) domains. P-values are given for all pleiotropic metatests for a given domain in a given Pathway (Fig. 1). Group: the type of an association defined in the "Antagonistic genetic heterogeneity in pleiotropic meta-analysis" section. P GRASP : p-values from the pleiotropic meta-analysis of results for the index SNPs reported in GRASP for phenotypes used in 12M or 20P domain reported in column "Domain". Empty cells: no index SNPs for phenotypes from the corresponding domain (shown in column "Domain") were reported in GRASP. N P and N G show the number of phenotypes that attained a minimal p-value (p<0.05) in the meta-or Fisher test in a given phenotypic domain in pathway 1 and Fisher test of GRASP results, respectively. P MOp , P MOb , P MFp , P FcFp , P OpFc , P ObFc , and P FpFc : p-values for meta-tests (subscripts) defined in Fig. 1. More details on the associations for these SNPs are given in Supplementary Tables 5 and 6. AGING 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 metaanalysis.  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 metastatistics for correlation among phenotypes ( -based omnibus test) and the effect statistics ( -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 omni-bus 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 -based omnibus tests (i.e., OpFc or MOp, Fig. 1 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   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 10 (p-value) trimmed at GW level -log 10 (5×10 -8 )=7.3 for better resolution.
AGING results in orthogonal directions of the bivariate vector of genetic effects β 1 and β 2 for phenotypes P1 and P2 and the vector of the correlation of these phenotypes and, accordingly, in opposite directions of vectors β 2 and P 2 . As seen from Fig. 3, the antagonistic heterogeneity implies antagonistic directions of genetic effects for directly correlated phenotypes (see Methods).
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).

DISCUSSION
Our univariate and pleiotropic meta-analyses of 20 agerelated phenotypes dealing with the natural-selectionfree genetic heterogeneity identified large number of non-proxy SNPs, 142 SNPs, with GW significance in a relatively modest sample of 33,431 individuals of Caucasian ancestry from five longitudinal studies ( Figure 5). Only 18 SNPs (12.7%) were identified in the univariate meta-analysis with most SNPs, 88.9% (16 of 18 SNPs), replicating previously reported associations, primarily with lipids. Two novel SNPs were associated with HC (rs6745983 in TMEM163 gene) and BG (rs10885409 in TCF7L2 gene). Two of 16 replicated SNPs, rs780094 (GCKR gene) and rs261332 (LIPC/LIPC-AS1 genes) associated with BG and TC, respectively, showed evidences for inter-cohort heterogeneity in the effect sizes. Accordingly, smaller p-values were attained by addressing this heterogeneity. The association of rs2155216 (BUD13-ZPR1/APOA5 gene locus) with TG was replicated with substantially higher (4-fold) efficiency in our analysis than in previous study [22] that was achieved, in part, by leveraging information on repeated measurements from longitudinal follow up.  Figure 5). Novel pleiotropic SNPs attained GW significance by combining associations with multiple phenotypes, which individually did not attain GW significance in univariate meta-analysis ( Figure 2). Therefore, pleiotropic meta-analysis has power to identify SNPs leveraging signals, which are often considered as noise in univariate GWAS. Most importantly, our analysis showed that the associations for most novel pleiotropic SNPs, 93.9% or 108 of 115 SNPs (and all 9 replicated SNPs) were strongly affected by the natural-selectionfree genetic heterogeneity in a rarely recognized form of antagonistic heterogeneity, implying antagonistic directions of genetic effects for directly correlated phenotypes. This heterogeneity has not been addressed in the currently prevailing GWAS.
Thus, our analysis provides compelling evidences that the traditional sample-size-centered GWAS approach is well adapted to handle homogeneous univariate associations (one SNP -one phenotype) within Mendelian framework whereas most associations with age-related phenotypes are more complex and well beyond that framework. This complexity substantially decreases the efficiency of the analysis of such phenotypes within the current GWAS framework. This conclusion leads to three major implications.
First, discovery of a large number of common genetic variants associated with different phenotypes in GWAS [25], raised concern known as a "missing heritability problem" [26]. The problem is that the estimated variance attributed to GWAS-identified common SNPs is only a fraction of the expected genetic variance from the estimates of heritability of a phenotype. This problem even leads to questioning whether genetics could be helpful for improving health care [27]. Our results show that in large part missing heritability problem for age-related phenotypes can be due to missing associations with complex genetic predisposition to these phenotypes.
Second, the traditional GWAS approach is built on the concept of a "true" or causal genetic effect (see the Introduction) whereas biologically plausible concept for age-related phenotypes is that based on the undefined role of evolution in establishing their molecular mechanisms. The letter implies that the same genetic variant can predispose differently to the same phenotype in different population groups, even of the same ancestry, e.g., throughout the life course and/or across generations. This concept particularly implies that "small" genetic effects (often expected in GWAS of age-related phenotypes) can be due to a complex superposition of large effects rather than small penetrance. Then, just relaying on large samples in GWAS of agerelated phenotypes without dissecting the naturalselection-free genetic heterogeneity is inefficient because increasing the sample size will also increase heterogeneity [28]. The current study supports this concept by highlighting strong impact of antagonistic heterogeneity, which is counter-intuitive in medical genetics but natural within the evolutionary framework. The antagonistic heterogeneity is also supported by the analyses of alleles from well-known apolipoprotein B and E genes [16,17].
Third, all pleiotropic SNPs attained GW significance in the two largest domains of 12 quantitative markers and all 20 phenotypes. These findings support the hypothesis that pleiotropic predisposition to age-related phenotypes can be driven by fundamental mechanisms associated with aging [8,18,19,29]. This attractive hypothesis in gerontology conceptualized as geroscience [28,30] is based on observations that agingrelated processes [31,32] are among the most important risk factors for geriatric diseases of distinct etiologies. The dominant role of antagonistic heterogeneity in AGING pleiotropic associations cautions against simplistic approaches in studies of pleiotropic effects on agerelated phenotypes and emphasizes the importance of personalized medicine which can potentially handle complexity of risk profiles on an individual basis [15,33,34]. Dissecting the role of antagonistic heterogeneity is particularly important as this provides a genetic basis for strategies focusing on anti-side-effect treatments in medical care.
Finally, our bioinformatics analysis of genes for pleiotropic SNPs identified three interlinked pathways, G-protein coupled receptor signaling, GPCR-mediated integration of enteroendocrine signaling exemplified by an L cell, and cAMP-mediated signaling. G-protein coupled receptors (GPCRs), members of the largest family of membrane proteins, are activated by a wide variety of external signals (e.g., light, odorants, hormones, neurotransmitters, pharmacologic agents, etc.) and modulate the activity of signaling pathways involved in most physiological processes. In particular, GPCRs partly mediate secretion of enteroendocrine peptide hormones glucagon-like peptide 1 (GLP-1) and peptide YY (PYY) by L-cells, relevant to food intake and blood glucose homeostasis [35,36]. GPCRs play a causative role in many human diseases and are the important drug targets [37,38]. Adenylyl cyclase, a membrane-associated enzyme, when activated by G proteins, catalyzes synthesis of important second messenger cyclic adenosine monophosphate (cAMP). cAMP is involved in diverse biological processes (e.g., glycogen metabolism, hormone regulation, gene expression, sensory signal transduction). The cAMP levels also partly mediate secretion of enteroendocrine peptide hormones linked to feeding/metabolic control [35] and are associated with immune function [39]. Top GO BPs (adenylate cyclase-activating G-protein coupled receptor signaling pathway and activation of adenylate cyclase activity) support an enrichment of BPs related to GPCRs function and cAMP synthesis. Thus, these results suggest an important role of G protein-dependent signaling and adenylyl cyclase activity necessary for the proper biological response of cells to hormones and other extracellular signals in pleiotropic effect on age-related phenotypes.

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 AGING 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-selectionfree heterogeneity in genetic predisposition to agerelated phenotypes (see below), we used four tests in pathway 1 and three tests in pathway 2 (Fig. 1).
In pathway 1, univariate (phenotype-specific) metaanalysis 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 Zscores [56]. In pathway 1a, the weighted average of the effect sizes was calculated as Σ /Σ with variance 1/Σ , where is the inverse variance of effect size in the cohort 1,5 for a given AGING 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 cohortspecific estimates. Accordingly, the Fisher test can indicate heterogeneity in genetic associations by providing smaller p-value than that from the fixedeffects 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][59][60].
For a certain SNP we have an estimated effect size and its standard error for the phenotype 1, in the cohort 1,5 , where is the number of phenotypes in study . The omnibus test statistic is constructed as where / is a z-score vector of associations of SNPs with phenotypes and 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 ( ), the test statistic follows a chi-squared distribution with degrees of freedom ~ , from which we obtained a combined p-value in the study . Accordingly, this statistic addresses the problem of multiple testing by increasing the number of degrees of freedom.

Correlation matrices
were estimated based on the effect statistics of simulated SNPs in association models with phenotypes (denoted ). 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 based on the 250 effect sizes of permuted SNPs with K phenotypes. Likewise, can be also constructed by evaluating the correlations between phenotypes [61] (denoted ). The correlation matrices, and , were constructed for pathways 1 and 2 ( Fig. 1) separately. Cohort-specific matrices for pathway 2 were constructed by evaluating correlations of the effect statistics in each cohort. For pathway 1, matrix 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, , were evaluated using phenotype measurements in each cohort separately and matrix for pathway 1, , 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 in twodimensional case, Either the or -based omnibus test (see above) penalize pleiotropic meta-statistics, if ̂ ̂ Σ , ̂ ̂ Σ 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 agerelated 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 ̂ ̂ Σ , ̂ ̂ Σ 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 pvalues 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 pvalues 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 2 . The I 2 can be interpreted as the percentage of the total variability in a set of effect sizes due to between-sample variability.
AGING Funding for SHARe genotyping was provided by NHLBI Contract N02-HL-64278. This manuscript was not prepared in collaboration with MESA investigators and does not necessarily reflect the opinions or views of MESA, or the NHLBI.
The Health and Retirement Study (HRS) genetic data is sponsored by the Genetics Resource with HRS April 21, 2010, version G Page 5 of 7 National Institute on Aging (grant numbers U01AG009740, RC2AG036495, and RC4AG039029) and was conducted by the University of Michigan. This manuscript was not prepared in collaboration with HRS investigators and does not necessarily reflect the opinions or views of HRS, or the NHLBI.

CONFLICTS OF INTEREST
The authors declare no conflicts of interest.

FUNDING
This research was supported by Grants No P01 AG043352 and R01AG047310 from the National Institute on Aging. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This manuscript was prepared using a limited access datasets obtained though dbGaP (accession numbers: phs000007.v22.p10 (FHS), phs000280.v3.p1 (ARIC), phs000209.v13.p3 (MESA), phs000287.v5.p1 (CHS), phs000428.v1.p1 (HRS)). Phenotypic HRS data are available publicly and through restricted access from http://hrsonline.isr.umich.edu/ index.php?p=data.