Tissue-specific impacts of aging and genetics on gene expression patterns in humans

Age is the primary risk factor for many common human diseases. Here, we quantify the relative contributions of genetics and aging to gene expression patterns across 27 tissues from 948 humans. We show that the predictive power of expression quantitative trait loci is impacted by age in many tissues. Jointly modelling the contributions of age and genetics to transcript level variation we find expression heritability (h2) is consistent among tissues while the contribution of aging varies by >20-fold with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${R}_{{{{{{{{\rm{age}}}}}}}}}^{2} \; > \;{h}^{2}$$\end{document}Rage2>h2 in 5 tissues. We find that while the force of purifying selection is stronger on genes expressed early versus late in life (Medawar’s hypothesis), several highly proliferative tissues exhibit the opposite pattern. These non-Medawarian tissues exhibit high rates of cancer and age-of-expression-associated somatic mutations. In contrast, genes under genetic control are under relaxed constraint. Together, we demonstrate the distinct roles of aging and genetics on expression phenotypes.

This study previously demonstrated a reduction in eQTL 138 heritability with age supporting our results. We confirmed 139 using our approach that eQTLs were less predictive of gene 140 expression in 80, compared to 70 year old's ( Fig.S5, S6). 141 These results suggest that the predictive power of eQTLs 142 is impacted by the sample age across the vast majority of   Examples of gene expression binned by genotype and age for four genes in which eQTL p-values change in whole blood and heart atrial appendage. Example genes in whole blood show lower p-values within young individuals and lower p-values within old individuals in heart atrial appendage. NFRKB and HSD17B12 show positive effect size while MIF-AS1 and ZC3H3 show negative effect size. Center line of the boxplot indicates median, box limit indicates first and third quartiles and points indicate outliers. C) QQ plots of eQTL p-values (plotted as -log(P)) for old (red) and young (blue) individuals from a linear model correlating expression with the lead SNP for each gene in all 27 tissues (Table S2). P-value for significance differences in eQTL p-value distributions are obtained from two-sided Welch's t-test.    associated with genetic-associated variance (e.g. Fig. 4A).

308
We found more age-associated pathway enrichment even in 309 tissues for which the average age-associated contribution to 310 gene expression was low (e.g. Pancreas, Fig. S25). This  Table S6). 400 We found that the per-gene age of expression (β age ) was sig- tissues. 410 We also explored the relationship between gene expression 411 heritability and constraint. Across all tissues h 2 was signif-412 icantly negatively correlated with pLI (P-value < 10 −3 lin-413 ear model two-sided t-test, Fig. S35, S36, S37). While 414 this trend was consistent across tissues, intriguingly it was 415 strongest in heart tissues. Thus, on average genes in which 416 the variation in expression levels is heritable tend to be un-417 der significantly less functional constraint (Fig. 5F). In con-418 trast however, on average R 2 age increases as function of pLI 419 (Fig. 5F), highlighting the increased constraint of many of 420 the genes that exhibit age-associated changes in gene expres- terns of aging will exhibit substantial cell-type specificity. 456 We also present a novel approach to jointly model the impact highlight that genes with eQTLs tend to be subject to less 506 evolutionary constraint, and thus potentially less biologically 507 important, in contrast to genes with age-associated gene ex-508 pression patterns which exhibit increased constraint.

509
The critical role of aging as a risk factor for many common analyses we split the GTEx data for each tissue into two age 530 groups, "young" and "old," based on the median age of indi-531 viduals in the full dataset, which was 55 (Fig. S1) (Fig. S38).
where P i is the distribution for individual i and H is the Shan-586 non entropy function: JSD is known to be a robust metric that is less sensitive to 588 noise when calculating distance compared to traditional met-589 rics such as Euclidean distance and correlation. It has been 590 shown that JSD metrics and other approaches yield similar 591 results but that JSD is more robust to outliers (12). The square 592 root of the raw JSD value follows the triangle inequality, en-593 abling us to treat it as a distance metric.

594
Slope of JSD distance versus age. In addition to compar-595 ing JSD between the two age groups defined above, "young" 596 and "old", we also binned all GTEx individuals into 6 age 597 groups, from 20 to 80 years old with an increment of 10 years. 598 We then computed pairwise distance and average age for each   (F DR<0.2). We compare this metric to the per-tissue 2-bin 638 JSD (Fig. S39) and 6-bin JSD slope (Fig. S13). 646 Where β i,g,t is the coefficient or effect size for SNP X i in (4) The minimization problem contains both the error of our 658 model predictions (Y j − β 0 − X T j β) 2 and a regularization 659 term λ( 1−α 2 ||β|| 2 2 + α||β|| 1 ) to prevent model overfitting.

660
The elastic net regularization term incorporates both L1 661 (||β|| 1 )) and L2 (||β|| 2 2 ) penalties. Following PrediXcan, Y for a gene g and tissue t is modeled as: Where A is the normalized age of an individual. Coefficients 680 were fit using elastic net regularization, as above, which sets 681 coefficients for non-informative predictors to zero. The sign 682 of the fitted age coefficient (β age,g,t ), when nonzero, reflects 683 whether the gene in that tissue is expressed more in young 684 (negative coefficient) or old (positive coefficient) individuals. 685 We also evaluated the training set R 2 using the fit model coef-686 ficients separately for genetics (across all SNPs in the model) 687 and age: We also tested whether the age-related gene expression rela-689 tionship was sex-specific by rerunning the joint model with 690 an additional age-sex interaction term as follows: Where β age * sex,g,t is the additional model weight for the age- ance explained by each term in the model. We compared the 696 R 2 age between the models including or excluding the age-sex 697 interaction term (Fig. S24). We also compared the tissue-

705
Tissue specificity of age and genetic associations. 706 We evaluated the variability of age and genetic associations 707 across tissues using a measure of tissue specificity for age 708 and genetic R 2 (36). We measured the tissue-specificity of a 709 gene g's variance explained R 2 g using the following metric: Where n is the total number of tissues, R 2 g,t is the variance  and repeated the analysis between pLI and β age (Fig. S40).

737
We also analyzed the evolutionary constraint metric dN/dS 738 (21) and its tissue-specific relationship with β age by deter-