Mutation accumulation differentially impacts ageing in mammalian tissues

Medawar’s mutation accumulation (MA) hypothesis explains ageing by the declining force of natural selection with age: slightly deleterious germline mutations that are functional in old age are not effectively eliminated by selection and therefore lead to ageing-related phenotypes. Although widely cited, empirical support for the MA hypothesis, particularly molecular evidence, has remained limited. Here we test one of its predictions, that genes relatively highly expressed in old adults vs. young adults should be under weaker purifying selection than those relatively highly expressed in young adults. To do so, we combine 23 RNA-sequencing and 35 microarray gene expression datasets (including 9 tissues from 5 mammalian species) with protein and regulatory sequence conservation estimates across mammals. We identify age-related decrease in transcriptome conservation (ADICT) in four tissues, brain, liver, lung, and artery, but not in other tissues, most notably muscle and heart. ADICT is driven both by decreased expression of highly conserved genes and up-regulation of poorly conserved genes during ageing, in line with the MA hypothesis. Lowly conserved and up-regulated genes in ADICT-associated tissues have overlapping functional properties, particularly involving apoptosis and inflammation, with no evidence for a history of positive selection. Our results suggest that tissues vary in how evolution has shaped their ageing patterns. We find that in some tissues, genes up-regulated during ageing, possibly in response to accumulating cellular and histological damage, are under weaker purifying selection than other genes. We propose that accumulation of slightly deleterious substitutions in these genes may underlie their suboptimal regulation and activity during ageing, shaping senescent phenotypes such as inflammaging.


Introduction
To date, more than 300 theories have been postulated to explain senescence, i.e. age-related loss of function and increase in mortality rates 1 . The mutation accumulation (MA) hypothesis, an evolutionary explanation for ageing first developed by J.B.S. Haldane 2 and Peter Medawar 3 , is among the most influential of such theories. It states that negative selection against germ-line substitutions that exhibit harmful effects only during old age will be inefficient. Such substitutions can eventually fix through genetic drift, and thus contribute to observed senescent phenotypes 4 . The MA hypothesis generates several testable predictions. One expectation is that genetic variance in fitness-related traits, such as reproductive success or survival, will increase with age 2,5 . A second expectation is that inbreeding depression will also increase with age. In line with these predictions, several studies have shown age-related increase in genetic variance in fitness-related traits in either laboratory (Drosophila melanogaster: refs. 6,7 ) [but see refs. [8][9][10] or natural populations (Soay sheep, red deer: ref. 11 ). In humans, heritability of CpG methylation patterns was shown to increase with age for about 100 genome-wide loci also consistent with MA, although possible fitness consequences were not evaluated 12 . In an indirect test of the expectation for inbreeding depression, outbreeding was reported to reduce age-related increase in mortality in hermaphroditic snails 13 , again in line with MA. Finally, Rodríguez et al. 14 studied >2,500 human genetic variants linked to 120 genetic diseases and reported that variants associated with late-onset disease segregate at higher frequencies than those associated with early-onset disease, a third prediction under MA 14 .
Beyond those cited above, relatively few studies have used empirical data to test the MA hypothesis. In particular, the possibly variable contribution of MA in the ageing processes affecting different species and their different tissues has not yet been tackled. In addition, we have limited understanding of the nature and prevalence of late-expressed substitutions, with the exception of a few extreme cases such as the CAG repeat variants in the huntingtin gene that cause Huntington's disease 5 .
The role of MA in ageing therefore awaits testing through new approaches that encompass a larger number of traits, a wider array of species, different tissues, and that include molecular data. One such approach would be to take advantage of widely available transcriptome data, in particular genome-wide gene expression datasets that comprise adult individuals of varying age. Such transcriptome datasets have traditionally been used to identify functional processes affected by or underlying senescence, although they can also be used to test evolutionary theories, as we show here.
In previous work, we used prefrontal cortex transcriptome age-series from humans to test whether protein sequence conservation varies among genes that are highly expressed at different ages 15 . This analysis showed that are relatively highly expressed genes in young adults vs. old adults are evolutionarily more conserved than those that are relatively highly expressed genes in old adults vs.
young adults, which we call age-related decrease in conservation of the transcriptome (ADICT), consistent with the MA hypothesis. However, this work analysed only one brain region and did not distinguish between two distinct processes: (a) up-regulation of lowly conserved genes with age and, (b) down-regulation of highly conserved genes with age. Both processes could cause the ADICT effect, but only (a) would be predicted by MA.
To address these limitations, here we expanded our analysis to include 9 different tissue types and 5 mammalian species. First, we investigated the prevalence of the ADICT pattern across multiple mammalian ageing datasets, using estimates of protein and regulatory sequence conservation across mammals. Second, we determined whether genes up-regulated late in life show low evolutionary conservation, as predicted by MA. In other words, we tested whether slightly deleterious mutations are more likely to fix in genes that are more highly expressed in old age, such as genes that respond to age-associated tissue damage 16,17 .

Age-related decrease in conservation of the transcriptome
We collected published transcriptome age-series of young and old adults of 5 mammalian species, generated using RNA-sequencing or microarrays (Homo sapiens, Macaca mulatta, Macata fascicularis, Rattus norvegicus, Mus musculus; n = 58 datasets and 2041 unique samples in all). The datasets include different brain regions (humans, macaques, rats, and mice), muscle (humans, rats, and mice), artery (humans, macaques, and rats), heart (humans), skin (humans and mice), kidney (humans and mice), liver (humans and mice), lung (humans and mice), and spleen (mice). Across all analysed datasets (n = 9 to 116 individuals, mean = 35.2), human ages range from 16 to 106 years, macaque ages range from 4 to 28 years, rat ages range from 3 to 30 months, and mouse ages range from 8 to 130 weeks (Supplementary Table 1).
We studied congruence in age-related gene expression change across the 58 datasets. First, for each gene in each dataset we calculated the Spearman correlation coefficient between gene expression level and individual age (ρ EA ). We then compared datasets to estimate pairwise similarity in ρ EA values across common genes. ρ EA values were mostly (72% of comparisons) positively correlated across data sets, indicating that the same genes' expression levels were similarly affected by age ( Supplementary Fig. 1).
As measure of gene sequence conservation we used estimates of purifying selection on protein sequence (ω 0 ), calculated by Kryuchkova-Mostacci and Robinson-Rechavi (2015) 18 and estimated for the human or the mouse branch using the branch-site model 19 . ω 0 is the dN/dS ratio calculated for those sites determined to be under purifying selection, and thus is expected to be a direct measure of the strength of purifying selection on a gene. We further calculated an adjusted protein conservation metric (ω 0 * ) for each gene, factoring out the possible effects of GC content, CDS length, intron length, intron number, mean expression, median expression, maximum expression, tissue specificity, network connectivity, phyletic age, and number of paralogs, using a multiple regression model following Kryuchkova-Mostacci and Robinson-Rechavi (2015) 18 . The value -ω 0 * (ω 0 * multiplied by -1) represents the main protein sequence conservation metric we used in our analyses, where more positive values represent more conserved genes.
We then investigated the prevalence of ADICT in mammalian ageing. To do so, we first calculated the Spearman correlation coefficient between gene expression level and the protein sequence conservation metric (which we call ρ EC ) across all genes, for each individual in each dataset (Figs. 1a, 1b). Note that the conservation metric (-ω 0 * ) is a constant value per gene, while gene expression levels will differ among individuals. In general, a positive ρ EC is expected, such that more highly expressed genes tend to be more conserved in their protein sequence 20 , but its degree may vary among individuals depending on age. To test this idea, in each dataset we determined the correlation between individual ages and ρ EC (ρ AρEC ). Fig. 1c provides an example of such a pattern in one brain ageing dataset 21   In each dataset, we used two gene sets for testing ADICT: (a) genes that showed significant agerelated change in expression levels (at Spearman correlation test q-value < 0.10), and (b) all expressed genes. We conducted analyses using all expressed genes in order to avoid a reduction in statistical power in datasets with low sample sizes and to determine whether patterns that hold for strongly age-associated genes also apply across the entire transcriptome (Supplementary Table 2).
Note that in 21 of 58 datasets, and in particular smaller datasets, we could not identify a set of significant age-related genes at q < 0.10; for these datasets we only conducted the analysis using all expressed genes (see Methods).
In the brain, for 18 / 19 datasets with age-related genes, we found negative ρ AρEC values consistent with the notion of ADICT. ρ AρEC values for 17 of these 18 datasets were significant at nominal p < 0.05. When repeating this analysis with all expressed genes, 27 / 28 datasets had negative ρ AρEC values, 16 of which were nominally significant. Together, these results support a general trend of ADICT in the brain (Fig. 2). We also found significant negative ρ AρEC values in the majority of liver (3 / 4) and lung (3 / 3) datasets, and in 2 of 5 artery datasets. In contrast, we found no consistent ADICT pattern in various muscle types (n = 10, where we identified only one nominally significant case) or in heart, skin, spleen or kidney. We also note that in 6 / 10 muscle datasets, we did not detect significant age-related expression change in gene-by-gene analyses (Supplementary Table 2).
Overall, 26 / 58 of the datasets showed significant ADICT signatures across age-related genes after correction for multiple testing (q < 0.10), with this pattern driven by brain, liver, lung, and artery ( Fig. 2). ADICT signatures were also consistent with our findings using -ω 0 * when we used the mean PhastCons score as the measure of negative selection on protein coding sequences ( Supplementary Fig. 7).
To determine the robustness of this result with respect to our protein-coding sequence conservation metric, we repeated the analysis using ω 0 values without applying multiple regression, as well as ω values obtained from the Ensembl database for "one-to-one orthologs" between human-mouse, human-elephant, and human-cow (i.e. raw dN/dS values). We further tested whether the trend holds when we exclude (a) putatively positively selected genes (with ω >1 in our data), (b) immune system genes known to be generally fast evolving 19,22,23 , and (c) genes down-regulated with age in each dataset (ranging from n = 1,086 to 6,717). In addition, to exclude the possibility that ADICT signals are driven by gene expression changes involving only few functional processes (e.g. highly conserved developmental genes being down-regulated), we calculated ρ AρEC separately for genes in each of the largest GO Biological Process categories (n = 19, each with node size >1000 annotated genes) ( Supplementary Fig. 2). Negative ρ AρEC values were repeatedly detected in the same 25 datasets (excluding one muscle dataset showing ADICT), irrespective of the metric used, the gene sets, and GO categories involved (Supplementary Table 2 and Supplementary Fig. 2). ADICT thus appears to be a consistent trend in multiple tissues, although it is also not a universal pattern across all mammalian tissue types.

Fig. 2 | Age-dependent changes in transcriptome conservation
. The x-axis shows the Spearman correlation coefficient (ρ AρEC ) between individual age and expression-conservation correlations (ρ EC described in Fig. 1). The statistics are calculated separately for each dataset, and for significant agerelated genes in that dataset (light bars), as well as for all expressed genes (dark bars). Note that in 21 of 58 datasets (cases where the light bar is missing), no significant age-related gene set could be identified (at q < 0.10). The asterisks indicate nominal significance of the Spearman correlation test, (*): p ≤ 0.05, (**): p ≤ 0.01, (***): p ≤ 0.001. In the analysis using age-related genes, the 25 datasets showing nominal significance for ADICT remained significant at q < 0.10 after applying Benjamini-Hochberg correction (excluding one muscle dataset showing ADICT). In the analysis using all genes, the 21 / 22 datasets showing nominal significance for ADICT remained significant at q < 0.10 after applying Benjamini-Hochberg correction, with "Mm_Liver3" the only exception.

Distinct processes contribute to ADICT
We next investigated two non-exclusive scenarios that could lead to ADICT: (a) genes that show age-related increase in expression could be lowly conserved, consistent with MA, or (b) genes that show age-related decrease in expression could be highly conserved, relative to genes that show no change in expression. The latter scenario could occur if a set of highly conserved genes (e.g. synaptic genes) are down-regulated during the postnatal lifespan, as previously reported 15,24 , but would not provide direct evidence in support of MA.
To test which of these scenarios underlie ADICT, we compared the mean conservation metric among (a) genes that show increases in expression with age (i.e. up-regulation, ρ EA > 0.1, q < 0.1), and (b) genes that show decreases in expression with age (ρ EA < -0.1, q < 0.1), using (c) genes that show no age-related changes in expression level as a control. We repeated this analysis across the 25 datasets showing the ADICT signature. We found results consistent with both scenarios (Fig. 3): genes that show decreases in expression with age were more strongly conserved than genes with no change (n = 22 / 25; 14 with bootstrap support >95%). Genes up-regulated with age were also more weakly conserved than genes with no change, in nearly all cases (n = 23 / 25; 15 with bootstrap support >95%). The latter observation is in line with the MA hypothesis.
We further hypothesized that genes that more frequently exhibit increased expression with age across tissues would be less likely to be conserved. To test this hypothesis, we selected genes shared across the 25 ADICT datasets and counted how many times each gene was up-regulated with age.
As predicted, we observed a negative correlation between the number of cases in which a gene was up-regulated with age and its conservation metric (ρ = -0.17, p < 0.001) (Supplementary Fig. 8a).

Functional analysis of ADICT
To find functionally coherent gene sets that may contribute to ADICT patterns in brain, liver, lung, and artery, we conducted Gene Ontology (GO) analysis for the 3 GO domains (Biological Process, Cellular Component, Molecular Function). We calculated statistics for GO analyses by ranking genes according to both the conservation metric (-ω 0 * ) and expression-age correlations (ρ EA ) and investigating GO term enrichment in the 10% tails of the distributions. We separately analysed (a) genes that showed increased expression with age and low relative conservation (IELC, consistent with MA), and (b) genes that showed decreased expression with age and high conservation (DEHC). We sought shared GO categories enriched in either IELC genes or DEHC genes across all 25 datasets showing the ADICT signature. To determine the random expectation for shared GO categories, we randomly permuted ages of individuals in each dataset 1000 times, calculated ρ EA again, and repeated the gene ranking and GO analysis.
IELC genes were enriched in the same 24 GO Biological Process categories in all the 25 datasets (expected = 0; permutation test p < 0.001) (Supplementary Figs. 3 and 4, Supplementary Table 4).
These included categories related to apoptosis, inflammation, and the immune response, among others (see the REVIGO summary in Supplementary Fig. 4). In addition, four GO Cellular Component categories (expected = 0; p < 0.001) and one GO Molecular Function category (expected = 0; p = 0.022) were shared among IELC genes across the 25 datasets ( Supplementary   Fig. 3). Among DEHC genes, we found shared enrichment only in Cellular Component and Molecular Function categories (permutation test p < 0.05); significant gene sets included synapse and signaling related functions ( Supplementary Fig. 4 and Supplementary Table 4). Although functional analysis suggests DEHC as a potential factor for functional decline, we will only focus on IELC as the objective of this study is to test MA hypothesis.

Age-dependent effects on regulatory region conservation
Finally, we asked whether ADICT also extended to conservation of transcriptional regulatory regions. To test this possibility, we calculated the mean PhastCons score (a metric for conserved elements) from the UCSC database based on a 100-way vertebrate alignment 25 , for (a) +/-2000 bp around the transcription start site (TSS) and (b) the 3'-UTR. We then repeated our ADICT analysis by substituting these two regulatory conservation metrics for -ω 0 * . This again revealed a heterogeneous trend toward ADICT across tissues, with consistent ADICT trends in brain, liver, and lung ( Supplementary Fig. 5 and Supplementary Fig. 6).

Discussion
The MA hypothesis predicts that burden of slightly deleterious germline substitutions will increase with age due to the declining force of negative selection 3 . Our approach differs from earlier attempts to test this hypothesis [6][7][8][9][10][11]13,14 in two respects. First, instead of relying on intra-species variation to estimate mutational load, we used inter-species divergence, which may be statistically more powerful as it involves a larger number of substitutions. Second, we studied the mutational load on multiple tissues, thus taking into account the possibility that age-dependent germline mutational load may vary among tissues, depending on tissue-specific developmental patterns, mitotic capacity, and consequences for organism-level fitness.
We found the age-related decrease in transcriptome conservation (ADICT) pattern in datasets from brain, liver, lung, and artery, consistently across human, macaque, rat, and mouse. Among datasets from all four of these tissues, genes that increased expression with age showed low conservation (IELC). Furthermore, we found that a shared trend of up-regulation among all nine tested tissues predicted lower evolutionary conservation. These results are consistent with the MA hypothesis. In addition, our functional analysis suggests that processes known to underlie senescent phenotypes in multiple tissues, apoptosis and inflammation 17 , are particularly influenced by age-dependent mutational load. Nevertheless, the methodology depends on expression data and if the function of a gene is modulated through other mechanisms, such as post-translational modifications or alterations in the interaction partners, these will not be captured in our study.
Two notes of caution are warranted. First, ADICT and IELC were not consistently detected in muscle, heart, kidney, skin, and spleen datasets, and the reason is unclear. These cases could represent false negatives due to experimental artefacts, or they might reflect heterogeneity of the ageing process. It is notable that we observed ADICT and/or IELC in only one of the 10 muscle datasets. This result could be related to a weaker signature of ageing in muscle; indeed, we could only identify age-associated genes in 4 / 10 muscle datasets (at q < 0.10), in stark contrast to brain datasets (18 / 19). This latter result cannot be trivially explained by reduced statistical power in the muscle data sets (sample sizes are comparable with those of brain), but could be due to increased inter-individual variation 26 in muscle ageing. The MA hypothesis assumes late age-specific mutations; if age-specific expression patterns are absent in a tissue, we would not expect the contribution of the MA process to that tissue's ageing. This could explain variation among tissues in ADICT prevalence.
Second, the IELC propensity could be related to up-regulation of genes evolving under positive selection, such as immune-related genes, instead of genes under weak negative selection. Hence, low conservation might reflect the accumulation of beneficial substitutions rather than weakly deleterious ones. This scenario would be consistent with the antagonistic pleiotropy (AP) hypothesis, which argues that substitutions that are positively selected for their early life benefits may be harmful in late life 27 . However, (a) our analysis is based on an estimate of negative selection rather than raw ω, and thus is not expected to be affected by positive selection; (b) when we removed genes with ω > 1 or immune genes from our analysis, we still found the same ADICT and IELC signals (Supplementary Fig. 10); (c) we compared IELC genes and 370 genes identified to be under positive selection in humans through multiple genome scans 28 , and we did not find more overlap than expected compared to the background set of all genes we analysed (Fisher's exact test p = 0.3). While IELC may still partly derive from as-yet undetected positive selection in early life and represent a case of antagonistic pleiotropy, it remains at least as likely that deficient purifying selection and accumulation of slightly deleterious substitutions 29 , as predicted by the MA hypothesis, underlies the observed IELC signal.
Might the IELC phenomenon we detect contribute to physiological decline in ageing? Our finding that inflammation and apoptosis are shared functional characteristics of IELC genes in four tissues is telling, especially given increasing appreciation of the role of inflammaging, i.e. low-level inflammation observed in many ageing tissues 16,17,30 . There are multiple examples of how chronic inflammation can impair housekeeping functions, especially in the brain (e.g. refs. 17,31 ). Meanwhile, apoptosis is crucial for eliminating senescent cells during healthy ageing, and disruptions in apoptosis could lead to accumulation of dysfunctional cells over time. Conversely, apoptosis is also thought to have a role in neurodegenerative disease aetiology, such as Alzheimer's Disease aetiology, by causing neuronal loss 32 . Our results suggest that genes involved in cellular and tissue level damage response, such as those with roles in inflammation and apoptosis, are subject to weaker purifying selection than other genes, possibly due to their limited activity early in life. The resulting mutational load leads to suboptimal regulation and function during ageing in particular tissues, when these genes show elevated activity. The MA process may thus contribute to mammalian senescent phenotypes, although at varying levels in different tissues.
Probeset-to-gene conversion. Affymetrix probe set IDs were converted to Ensembl gene IDs using the Bioconductor "biomaRt" package 57 . We used the "useMart" function to select the dataset for the species of interest, and the "getBM" function to retrieve the Ensembl gene IDs corresponding to Affymetrix probe set IDs. We then followed two steps: (a) if one probe set corresponded to more than one Ensembl gene, we removed that probe set, (b) if >1 probe set corresponded to one Ensembl gene, we chose the probe set which had the maximum expression value across all samples in that dataset. This approach used information only from the highest expressed and best-measured transcript per gene in each dataset (in other words, we discarded information from more lowly expressed and possibly noisy transcripts in that dataset).
Age test and age-related gene sets. Genes showing age-related changes in expression levels were identified using the Spearman correlation test. We used the R "cor.test" function using "method='Spearman'" for calculating the age-expression correlation coefficient ρ EA . The p-values were corrected for multiple testing using the "p.adjust" function with the "Benjamini-Hochberg (BH)" method in R, yielding q-values as a measure of the false discovery rate. We used the nonparametric Spearman rank correlation test to overcome several problems related to conducting meta-analysis (e.g., each dataset displays unique and sometimes non-normal distributions; outliers can have large effects on data analysis). We used a q-value cutoff of q < 0.10, which is a commonly used threshold (e.g. refs. 15,58 ). Among 58 datasets, 21 had a low number of age-related genes (n < 50), so to limit type II error we did not include these datasets in analyses of age-related gene sets.
Gene set sizes for age-related genes and all detected genes for all 58 datasets are shown in Supplementary Table 2. Unsurprisingly, the number of age-related genes is partially affected by sample size (at Spearman correlation test ρ = 0.35, p = 0.03), but this does not influence our ability to observe ADICT (see below).
We then defined two further categories based on the expression-age Spearman correlation coefficient (ρ EA ): (a) genes that showed age-related increases, with ρ EA > 0.1 and q < 0.1; (b) genes that showed age-related decreases, with ρ EA < -0.1 and q < 0.1. In addition to the q-value to define our cutoffs, we also used the correlation coefficient (ρ EA ). Therefore genes with small effect size that can be identified in large datasets (i.e. with high power) but not in small datasets will not be included in our analysis. Genes with q > 0.10 were considered to have no change in expression level with age. Genes with q < 0.10 and |ρ EA | < 0.1 were discarded from further analysis.
ADICT. The ADICT pattern was calculated as the Spearman rank correlation between age and ρ EC in each dataset, using all genes in each dataset, and using age-related genes, if detected in that dataset. The Spearman p-values were corrected using the BH method as described above (across all datasets included in an analysis). Note that correlation between |ρ AρEC | and sample size across datasets were negative (ρ AρEC calculated for all genes: ρ = -0.48, p < 0.001; ρ AρEC calculated for agerelated genes: ρ = -0.64, p < 0.001). This is simply because finding large correlation coefficients is unlikely with large sample sizes. However, this pattern cannot explain why we observe a consistent trend for negative ρ AρEC values (i.e. ADICT) in some tissues: e.g. in brain 27 / 28 datasets show a negative ρ AρEC whereas only 4 / 10 muscle datasets show a negative ρ AρEC .
Protein sequence conservation metrics. We used several types of metrics to estimate negative selection pressure on protein coding sequences.
First, we used ω 0 , a statistic based on coding sequence alignments across mammalian species. ω 0 is estimated for the Homininae branch for human and the Murinae branch for mouse, using the branch-site model 22  This measure of ω 0 can vary among genes due to multiple factors that are not the focus of this study 18 . To disentangle the effects of such factors from the effect of protein sequence conservation per se, we used information on GC content, CDS length, intron length, intron number, mean expression, median expression, maximum expression, tissue specificity, network connectivity, phyletic age, and number of paralogs, which were directly obtained from the Supplemental Material of Kryuchkova-Mostacci and Robinson-Rechavi (2015) 18 . To remove the effect of these variables from ω 0 , we used the "lm" function in the R "stats" package to calculate the residuals (ω 0 * ) from a multiple regression model with ω 0 as the response variable and all other variables as predictors. The ω 0 * statistic was calculated separately for human and for mouse ω 0 values. We used the human ω 0 * data in analyses involving primate transcriptome datasets, and the mouse ω 0 * data in analyses involving rodent transcriptome datasets.
Second, we calculated conservation in protein coding regions between pairs of species separated by different evolutionary distances, using dN (nonsynonymous substitution rate) and dS (synonymous substitution rate) statistics downloaded from Ensembl Biomart (v.83) 60 . Here we used "one-to-one orthologs" between human-mouse, human-elephant, and human-cow, in order to identify whether evolutionary distance between species affects estimated levels of sequence conservation. Because dN/dS ratios measure both the strength of negative and of positive selection, we repeated our analysis only using genes with dN/dS < 1 (i.e., excluding the genes most likely to evolve under positive selection). In addition, GO categories and subcategories related to immune system genes ("GO:0002376"), which are known to be fast-evolving, were selected using the R "get" function.
We then repeated the analysis after discarding these genes.
Finally, we calculated the conservation of protein coding sequences using the PhastCons scores (phastcons100way) downloaded from the UCSC database 25  Using the BEDTools 61 software package we obtained a list of all PhastCons scores (phastcons100way) for the defined 3'UTR bases or promoter bases of each gene. We then calculated the mean PhastCons score value as a metric to represent that gene's 3'-UTR or promoter region conservation.
Bootstrapping. Bootstrapping was performed using the "sample" function in R, with "replacement=TRUE". We used bootstrapping to calculate 95% confidence intervals for the mean conservation metric among genes that showed (a) age-related increases in expression levels, (b) age-related decreases in expression levels, and (c) no age-related changes in expression levels. For each case we resampled genes 1000 times, and calculated the mean. To visually compare the conservation metric among datasets, we then subtracted the median for genes that showed no agerelated change, from genes that showed age-related increase or age-related decrease. The upper and lower 2.5% quantiles are plotted in Fig. 3.

Testing linear changes.
To determine whether the relationship between individual age and ρ EC (calculated across age-related genes) was linear across adulthood, we compared linear regression models and quadratic regression models for each dataset, with ρ EC as the response variable and age as the explanatory variable (using the R "lm" function).
Defining IELC and DEHC gene sets. We developed a non-parametric statistic, z, which captures the relative correlation coefficient between gene expression and age, and relative conservation level of a gene simultaneously: where x is the rank of a correlation coefficient between gene's expression level across all detected genes in a dataset and age, and y is the rank of that gene's conservation metric. Using squared values gives additional weight to differences between higher ranks. High values of z indicate genes that have relatively high expression and low conservation, whereas low values of z indicate genes that have relatively low expression and high conservation. After sorting z values, the top 10% of genes were included in the increasing expression and low conservation (IELC) gene set and the bottom 10% were included in the decreasing expression and high conservation (DEHC) gene set.

Gene Ontology analysis.
Here we sought to find functional groups associated with either IELC or DEHC patterns that were shared across datasets of a tissue, and across all datasets. We conducted Gene Ontology (GO) analyses for 3 GO domains -Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). For this, we (I) chose GO groups showing enrichment tendencies in each dataset, using liberal cutoffs, (ii) determined the overlap among chosen GO groups among datasets, (iii) tested significance of the overlaps using random permutations of individual age in each dataset. Specifically, in each dataset, we chose GO groups with an odds ratio >1, comparing either IELC or DEHC genes (10% tails of the z statistic described above) to the rest (90%). We preferred not to use a p-value cutoff and to use liberal odds ratio cutoff (>1) in order to avoid type II error and to ensure that datasets with different numbers of genes contributed equally to downstream analysis. We then counted the number of overlapping GO groups that were thus chosen (odds ratio > 1) across all 25 datasets, or among different datasets for the same tissues. Next, we randomized ages of individuals in each dataset by conducting 1000 permutations, calculated expression correlations with age, and repeated the GO analysis using these correlation values. We finally compared the number of GO groups that showed enrichment tendency (odds ratio > 1) in the random permutations, with the observed values.