The Gene-Regulatory Footprint of Aging Highlights Conserved Central Regulators.

Many genes and pathways have been linked to aging, yet our understanding of underlying molecular mechanisms is still lacking. Here, we measure changes in the transcriptome, histone modifications, and DNA methylome in three metabolic tissues of adult and aged mice. Transcriptome and methylome changes dominate the liver aging footprint, whereas heart and muscle globally increase chromatin accessibility, especially in aging pathways. In mouse and human data from multiple tissues and regulatory layers, age-related transcription factor expression changes and binding site enrichment converge on putative aging modulators, including ZIC1, CXXC1, HMGA1, MECP2, SREBF1, SREBF2, ETS2, ZBTB7A, and ZNF518B. Using Mendelian randomization, we establish possible epidemiological links between expression of some of these transcription factors or their targets, including CXXC1, ZNF518B, and BBC3, and longevity. We conclude that conserved modulators are at the core of the molecular footprint of aging, and variation in tissue-specific expression of some may affect human longevity.

Introduction

Aging is a multifactorial process characterized by the gradual decline in vitality and is accompanied by increased susceptibility to a wide range of pathologies, including cancers and neurodegenerative, cardiovascular, metabolic, muscular, and infectious diseases (Kenyon, 2010). Quantifying, analyzing, and understanding these complex processes are critical to modulate and manage its negative ramifications to increase health span.

Systems-level characterization of aging at the levels of the transcriptome and proteome have identified common pathways and signatures across species (Kenyon, 2010) as well as molecular phenomena, including genomic instability, epigenetic alterations, loss of proteostasis, and mitochondrial dysfunction (Kenyon, 2010; Riera et al., 2016). By interfering with such mechanisms, such as cytosolic or mitochondrial proteostasis, lifespan and/or health span can be prolonged in animal models (Chondrogianni et al., 2015; Houtkooper et al., 2013; Sorrentino et al., 2017; Zhang et al., 2016). However, some aging-driven changes are tissue specific, suggesting that different tissues of the same organism age differently (Schumacher et al., 2008; Ori et al., 2015). In addition, although most studies on aging signatures rely on gene expression, epigenetic alterations are a major nexus of genome stability and transcriptional control (Benayoun et al., 2015; Johnson et al., 2012; Pal and Tyler, 2016). DNA methylation levels at some CpG sites can accurately predict an individual’s chronological age (Horvath, 2013; Horvath et al., 2016), yet the functional consequence of these changes is not clear. There are also many age-associated epigenetic changes at the level of histone modifications and composition, yet their relevance in mammalian lifespan or health span remains to be established (Booth and Brunet, 2016; Morris et al., 2018; Tvardovskiy et al., 2017).

It is therefore not clear whether aging in different tissues of different organisms shares a common denominator in terms of the gene-regulatory drivers. Because of the inherent complexity of the aging process, answering this question necessitates a multi-layer and multi-tissue analysis as well as integration of data from different species. Here, we characterized the age-related epigenetic and gene expression changes in three mouse tissues (liver, heart, and quadriceps muscle). Although there are many tissue-specific molecular differences, common biological processes are affected across these layers, which are modulated by a central set of transcription factors (TFs). We extend these findings to large human population datasets, in which we identify the same gene-regulatory drivers. Finally, we establish epidemiologically relevant genetic links between specific TFs and their targets with human longevity through Mendelian randomization (MR). Thus, our integrative investigation of the gene-regulatory footprint of aging identifies regulatory drivers that are shared by different molecular layers, tissues, and species, with some of them potentially explaining genetic variation in human longevity.

Results

Multi-layer Characterization of Aging Reveals Tissue-Specific Gene Regulatory Differences

Liver, heart, and quadriceps muscle were harvested, and their transcriptome, DNA methylome, and histone modification profiles were measured in adult (6 months) and old (24 months) male C57BL/6J mice, an age range comparable with 20–80 human years (Dutta and Sengupta, 2016), through genome-wide profiling of five different omic layers: the transcriptome, the DNA methylome, and three histone modifications (STAR Methods; Figure 1A). We quantified both positive (H3K27ac and H3K4me3) and negative (H3K27me3) histone marks using chromatin immunoprecipitation (ChIP) sequencing. We performed peak-based and gene-based differential binding analyses (Figure S1; STAR Methods). To simplify integration of all the layers, we chose to use the gene-based results, in which we quantified the reads in windows spanning 5 kb around gene transcription start sites (TSSs) and performed differential binding. Similarly, we performed differential methylation analyses on CpG sites and aggregated results by gene in two ways: upstream of the gene (promoter) and in the gene body (GB).

Figure 1

Tissue-Specific Multiomic Footprint of Aging

(A) General outline of the study, starting with multi-layer profiling experiments of mouse tissues followed by the identification of gene regulatory modifiers of the molecular footprint of aging. Results in the mouse were then replicated in human population data. Finally, using Mendelian randomization, genetic, tissue expression, and longevity trait data were integrated to identify candidate aging regulators of epidemiological significance.

(B) Top panel: schematic of differential analysis comparison across layers. For each measured layer, a p value threshold of 0.05 is used to define genes that go up or down. Then, the resulting set is overlapped with all other sets across layers to obtain pairwise intersections represented in the graph as circles of different sizes in the scatterplot. Bottom panel: liver, heart, and quadriceps summary of differential analysis results showing the number of genes that are either up or down upon aging for each of the layers (gray circles; size reflects number of genes). The intersection between different layers is computed and listed, as well as its significance, which is indicated by the color code (Wang et al., 2015).

More liver genes are affected in the transcriptome and methylome levels than in the histone epigenome (false discovery rate = 10%; Figure 1B; Table S1). In contrast, the heart and quadriceps aging footprints are dominated by changes in the methylome and histone modification profiles (Figure 1B). There is more DNA hyper- than hypomethylation in liver and quadriceps but not in the heart (Figure 1B).

As for the relationship between the different layers, there is a significant overlap between differentially expressed (DE) and differentially methylated (DM) GBs and promoters in the liver (Figure 1B). In addition, the overlap between DM in the GB and gene expression is directional, with more genes having increased expression having increased rather than decreased GB and promoter methylation. Compared with the GB, change in promoter methylation has a smaller overlap with expression. In heart and muscle, we observe fewer overlaps between gene expression changes and the other epigenetic layers, likely because of the smaller number of DE genes (Figure 1B). In all tissues, overlaps between DM and differential histone modifications are directional, in that many genes with increased GB methylation also have increased H3K27ac in their promoters and vice versa. In addition, the promoter DM overlaps well with that of the GB, even though these signals are traditionally regarded as having opposing effects on transcription (Jones, 2012). From our observations in these three metabolic tissues, we conclude that some omic layers are affected by age in a tissue-specific manner. Hence, the conclusions based on one tissue may not apply to another, even when exploring fundamental aspects of gene regulation and chromatin structure.

Aging Affects Similar Biological Processes at Distinct Regulatory Layers in Different Tissues

In order to gain a higher level understanding of the molecular pathways affected by aging, we performed gene set enrichment analysis (GSEA) using Gene Ontology (GO) terms of biological processes (Figure 2A). The enrichments mirror the differential analyses in that the liver transcriptome has many enrichments, whereas quadriceps and heart have enrichments mainly at the histone modification level. Within each layer, we identified many significant terms that have been previously linked to aging, including an upregulation of the immune response and downregulation of telomere-related genes, development, protein quality control, mitochondrial processes, and RNA processing (Figure 2B; Table S2). We observed no functional enrichments for changes in H3K4me3 in our data. In heart and muscle, genes with increased H3K27ac and decreased H3K27me3 are enriched for many overlapping biological processes. This emerging pattern of GO enrichments points to possible functional implications of open chromatin in aging.

Figure 2

Integrated Overview of Gene Set Enrichment Results across Biological Layers and Tissues

(A) Schematic of gene set enrichment analysis (GSEA) and integration across different layers. For each omic layer, the aging-driven log fold change is used to order genes. Then GSEA is performed using Gene Ontology biological process terms. Gene sets representing GO terms are represented as circles that are ordered on the basis of their significance and direction of enrichment, with green and red circles representing gene sets that increase and decrease with age, respectively. The size of the circle reflects the number of genes in the gene set. The overlaps between the core enrichment genes of each gene set in each layer and all other layers are computed. Links are drawn between sets of genes that have more than 20% overlap in their core enriched genes and are colored differently for each layer pair to facilitate visual discrimination.

(B) Hive plot representation of the multi-layer enrichment results in the liver, heart, and quadriceps muscle. Each circle represents a Biological Process term from the Gene Ontology consortium that has an enrichment q value less than 0.1. Purple boxes list some representative gene sets that were manually selected, most of which have already been linked to aging. The signed −log10(q), which is mapped to the color, is defined as the sign of the GSEA normalized enrichment score, multiplied by −log10(q).

We explored commonalities in gene set enrichments across layers by calculating cross-layer overlaps of genes that are driving the enrichments (Figure 2B). H3K27ac and H3K27me3 share most of their enrichments in the heart, and to a lower extent in the muscle, but in an opposite direction as expected. As for the methylome, we observe enrichment in the GBs only in liver, in which we see hypermethylation of genes regulating ion transport and interleukin regulation. Furthermore, the genes driving these enrichments are upregulated at the transcriptome level, suggesting that aging leads to their GB hypermethylation and increased expression.

Interestingly, some of the same biological processes that are enriched in the liver transcriptome are also enriched at the histone modification level in heart and quadriceps, suggesting that the same biological processes may be affected by age at different gene-regulatory layers. For example, the biological process enrichment scores in H3K27ac in the heart positively correlate with those of the liver transcriptome. This correlation is even stronger than most pairwise correlations within the liver (Figure S2).

Taken together, our multi-layer and multi-tissue analysis reveals that aging affects different facets of gene regulation in a tissue-dependent manner. Whereas liver has the strongest gene expression and DNA methylation effect, heart and muscle show strong epigenetic alterations, mainly toward a gain in activating and loss of repressive marks. Even though different tissues react at different layers, the aging footprint converges on common biological processes.

The Epigenetic Footprint of Aging Is Tissue Specific

The epigenetic footprint of aging is strikingly tissue specific. For instance, liver and quadriceps have more hyper- than hypomethylated CpG sites (Figure 3A). These changes are concentrated around the TSS of genes, peaking downstream of the TSS, likely in the GB (Figure 3B). Muscle DM has a similar trend, with more diffused localization of hypermethylated sites around the TSS. The heart, on the other hand, has similar amounts of hyper- and hypomethylated CpG sites (Figures 3A and 3B). Our liver DM data correlate best with the weights associated with the CpG sites in two published epigenetic clocks, a general and a liver-specific clock (Horvath et al., 2016; Stubbs et al., 2017; Wang et al., 2017; Figure 3C; Figure S3A). The reason behind these tissue-specific correlations with the general clock is not clear and may stem from different factors. It may be that the liver’s epigenetic landscape in the ages we are measuring is indeed the most affected by age. Conversely, this can be due to the different ages at which the clocks have been derived (newborn to 41 weeks old) and the presence of more liver samples than samples from other tissues.

Figure 3

The Epigenetic Footprint of Aging Is Tissue Specific

(A) Histogram of CpG sites with a methylation difference greater than 10% in the three tissues. Aging-driven hypo- and hyper-methylated sites are colored differently. Note the greater number of hypermethylated CpG in liver and quadriceps.

(B) Metaplot of the locations of hypo- and hypermethylated CpG sites with respect to the TSS of genes showing that age-related hypermethylation is enriched around the TSS.

(C) Correlation between the published general epigenetic clock CpG site weights (Stubbs et al., 2017) and the methylation differences in this study (Pearson correlations: for liver, r = 0.49, n = 162, p = 5.1e−11; for heart, r = 0.096, n = 159, p = 0.23; for quadriceps, r = 0.16, n = 160, p = 0.047). The mouse general epigenetic clock is a model in which methylation of a selected set of CpG sites predicts chronological age. The weight of each CpG site determines its relative contribution to the prediction.

(D) Number of differentially bound peaks as a function of distance from the closest TSS and false discovery rate (FDR). Bars pointing downward represent peaks whose binding decreases with age. The dashed vertical red lines demarcate the 10 kb region used in the TSS-based differential binding analyses.

In terms of histone modification changes, there are more increased than decreased H3K27ac peaks upon aging. These peaks are predominantly overlapping or within 5 kb of a TSS (Figure 3D). H3K27me3, being a broad peak, shows a more diffuse pattern with respect to TSS. However, we observe an abundance of decreased H3K27me3 peaks around the TSS in the heart and quadriceps. As for H3K4me3, we observe the greatest changes in the quadriceps, with a reduction in this mark around and downstream of the TSS.

It is not clear whether these tissue-specific differences are due to specific changes in histone modifications at certain peaks or more global changes in histone composition. To address this, we used a dataset in which histone H3 post-translational modifications were quantified using mass spectrometry in different tissues of mice from different age groups (Tvardovskiy et al., 2017; Figure S3B). Although total H3K27ac and H3K27me3 are negatively correlated in adult livers, this relationship is abolished or even reversed in old livers. Interestingly, the relative abundance of heart H3K27ac increases in old age, in contrast to the H3K27me3, which remains stable (Figure S3B). This tissue specificity is in line with the observations from our sequencing-based results and demonstrates how different tissues exhibit qualitatively distinct aging-driven epigenetic changes summarized by increased DNA methylation around promoters in liver and quadriceps and global differences in H3K27ac and H3K27me3 in heart and muscle.

TF Motif Enrichment across Layers and Tissues Identifies Candidate Central Aging Regulators

To pinpoint molecular players that may drive the gene-regulatory footprint of aging, we performed TF motif differential enrichment analysis in each of the layers using the HOCOMOCO-v10 motif database (STAR Methods; Figure 4A; Table S3). The H3K27ac and H3K27me3 TF motif enrichments show opposing signals in heart but not in quadriceps and liver (Figure 4C; Figure S4B). In addition, heart transcriptome enrichments are related negatively to that of H3K27me3 and positively to that of H3K27ac, meaning that TF motifs enriched in genes that increase in expression are also enriched in genes that have an increase in H3K27ac and a decrease in H3K27me3. In the liver, however, the transcriptome enrichments relate positively with the H3K27me3 repressive mark, which is unexpected. The observed relationship between these different layers is therefore tissue dependent. In addition, the enrichments based on changes in H3K4me3 show opposite relationships with other layers in liver versus heart and quadriceps

Figure 4

Enrichment of TF Motifs across Different Layers and Tissues Identifies Candidate Regulators of Age-Driven Transcriptional Changes

(A) Schematic of the TF motif differential enrichment analysis. For each layer, 10 kb sequences centered around the TSS of genes that are up or down (FDR = 10%) are scanned for differential TF enrichment. In the diagram, TFx represents one of the 426 TF motifs from the mouse HOCOMOCO-v10 database (Kulakovskiy et al., 2013). TF motifs are represented as circles ordered by differential enrichment on each axis, with size and color indicating differential enrichment. Links connect the same TF across layers and are drawn only when a TF has a differential enrichment greater than 10 on the log scale. Links are colored differently for better visual discrimination.

(B) Detailed view of the differential enrichment of some top DNA-binding motifs across layers in liver and quadriceps.

(C) Same results as in (B) represented in a hive layout to show relationships between transcription factor enrichments across layers. To highlight the most significant TFs, links are drawn between pairs of layers only when the motif had a differential enrichment greater than 10 in at least one of the two assays. Links are colored differently for each layer pair for visual discrimination.

We used a ranking scheme to obtain the top TFs with increased or decreased age-dependent enrichment (Figure S4; STAR Methods). For example, ZIC1, ZIC2, ZIC3, KLF6, PLAG1, and IKZF1 are generally associated with genes that increase in expression or open chromatin with age, whereas HMGA1, MECP2, CXXC1, SRY, MSX2, TBP, and CDX1 are associated with genes that decrease in expression or open chromatin with age. Three of the most significant and pervasive motifs are those of three members of the zinc finger of the cerebellum (Zic) family, which have similar binding motifs and are enriched in upregulated genes in all tissues as well as genes with increases in H3K27ac in the heart and decreases in H3K27me3 in the heart and quadriceps (Figures 4B and 4C). On the opposite side, the TF motif of HMGA1 is associated with genes that are silenced or decrease in expression with age (Figures 4B and 4C). Collectively, these results raise the possibility that although aging affects the measured layers in a tissue-specific manner, the molecular landscape of aging may be modulated by common regulators.

Orthologs of Mouse Age-Related TFs Have an Age-Related Expression Pattern in Human Populations

We then analyzed their expression levels in large population-based human datasets in which age is available. Particularly, we used the Genotype Tissue Expression (GTEx) dataset (GTEx Consortium, 2013), the Human Liver Cohort (Schadt et al., 2008), and a compendium of skeletal muscle datasets (Su et al., 2015). We selected liver, left ventricle of the heart, and skeletal muscle from GTEx and performed differential expression analysis to estimate the effect of age on each gene (Figure 5A; Table S1). Interestingly, the genes encoding many of the age-related TFs are significantly affected by age, and the directionality of this effect (i.e., whether the TF increases in expression with age) reflects the TF’s enrichment direction in the mouse. For instance, ZIC1 significantly increases in expression with age in human liver and muscle and is known to be an activator, and its motif is enriched in genes that increase in expression in mice. We found that ZIC1 expression, but not that of ZIC2 or ZIC3, correlates with age in multiple tissues and cohorts (Figure 5B; Figure S5A). HMGA1, on the other hand, decreases in expression with age, which may explain why its motif is enriched in genes whose expression decreases with age. Finally, we identified ZNF518B as the TF with the best correlation with age in the Human Liver Cohort, which is followed by ZIC1, the second best correlating gene. ZNF518B expression significantly decreases with age in all datasets and tissues (Figure 5). Unfortunately, the binding motif of ZNF518B is unknown, and therefore it is impossible to estimate its enrichment like the other motifs. Collectively, on the basis of transcriptome data from human populations, we observe that the expression of many of the mouse aging-related TFs is age dependent.

Figure 5

Expression Levels of the Human Orthologs of Age-Related TFs Identified in the Mouse Are Significantly Affected by Age

(A) Volcano plot of the age effect on the expression of genes in liver, heart, and skeletal muscle from the GTEx database. The log2 values for fold change represent the change in expression per year. The labeled genes correspond to TFs, which have on average the highest rank in terms of absolute value of differential enrichment across mouse tissues and layers. In addition, ZNF518B is labeled because it is consistently downregulated with age, even if it could not be included in the mouse TF enrichment analysis (see text). Genes that are differentially expressed in the human and mouse at 0.05 and 0.1 FDR thresholds, respectively, are marked in black, and their symbols are listed in the boxes. The number of significant genes scales with the sample sizes of the different tissues (153, 272, and 491 for liver, heart, and muscle, respectively). See STAR Methods for details.

(B) mRNA Expression level of ZIC1 and ZNF518B as a function of age in human data available in the GTEx Consortium, the Human Liver Cohort (Schadt et al., 2008), and a large meta-analysis of human skeletal muscles (Su et al., 2015).

Identified TFs Are Conserved across Aging Studies and Species

Having established that a set of core TFs may drive aging in different tissues in the mouse and that the expression of some of these TFs is significantly affected by age in humans, we explored whether the same TF enrichments can be replicated in another mouse study and in human GTEx data. For that, we applied the same TF enrichment strategy on human GTEx data as well as on an independent RNA sequencing experiment from the Brunet laboratory consisting of different C57BL/6J mouse tissues collected at 3, 12, and 29 months (Benayoun et al., 2019). We identify very similar TF motif enrichments, indicating the conservation of the gene-regulatory landscape changes across independent mouse studies as well as in human populations (Figure 6A; Table S3).

Figure 6

TFs and Their Targets Are Conserved across Aging Studies and Species and Are Genetically Linked to Aging in Large Human Cohorts

(A) We compared our data, using the same differential enrichment pipeline, with another extant mouse RNA sequencing dataset (Benayoun et al., 2019) as well as human GTEx data, revealing similar TF enrichments. The ten most differentially enriched TFs, defined as the average rank of the absolute value of their enrichment in the different layers, are labeled for clarity (STAR Methods). Links connect the same TF in the different layers to highlight the conservation across different studies and tissues.

(B) Two-sample Mendelian randomization results taking the top ten TFs and their top ten targets in three human tissues as exposures and four aging-related GWAS from the UK Biobank as outcomes (STAR Methods). The graph shows only TFs and targets that have at least one significant MR regression (random-effects inverse variance-weighted regression p < 0.05). The color of the MR link represents the direction of the MR effect size and its width represents the log10 of its p value. TFs and their predicted targets are connected by a dotted line. The node sizes and colors represent the fold change and p value attributed to age.

(C and D) Examples of positive and negative MR associations with age-related TFs or their targets. MR scatter and forest plots for ZNF518B and BBC3 in the liver versus mother’s age at death.

We generated an updated list of age-related TFs that are consistently enriched in different mouse studies and human populations (Figure 6A; STAR Methods). The majority of TFs share the same direction of enrichment in the different layers and tissues. Namely, the ZIC family motifs rank highly in all analyses, further pointing to a possible central role in aging. Other TFs, such as CXXC1 and MECP2, differ in directionality, where their enrichment can be highly negative or positive, but always at an extreme. Other top TFs include SREBF1, SREBF2, ZBT7A, ETS2, and IKZF1, all with possible links to aging. The sterol regulatory element proteins SREBF1 and SREBF2 regulate lipid homeostasis and therefore may have wide-ranging effects in the aging process (Shao and Espenshade, 2012). ZBT7A (aka Pokemon) is a transcriptional repressor for a wide range of metabolic genes, with implications in cancer (Liu et al., 2014). The gene expression of this ZBTB7A repressor decreases with age (Figure 5A), and its motif is enriched in genes that increase in expression (Figure 6A). ETS2 levels in the heart explain longevity variation in rats through activation of necrosis (Sheydina et al., 2012). IKZF1 (aka IKAROS) is associated with chromatin remodeling, and its targets are associated with neurodegenerative disease (Li et al., 2014). Taken together, our cross-species approach highlights conserved central gene regulatory players in aging and may provide more insights into the progression of this complex phenomenon.

Genetic Evidence for Involvement of Aging Regulators and Their Targets in Human Longevity

We next sought to find genetic and epidemiological evidence for the involvement of these putative players in human longevity. For that, we performed two sample MR (Hemani et al., 2018) taking the ten top (Figure 6A) TF expression levels in human tissues from GTEx as exposures and longevity traits from the UK Biobank (Sudlow et al., 2015). In addition to the top TFs, we also used a list of top predicted target genes (see STAR Methods) in the analyses. To obtain genome-wide eQTL associations for each gene, we performed a modified pipeline from GTEx in which age was taken as a covariate prior to removing unwanted sources of variation (STAR Methods). For each gene in each tissue, we took all associations with nominal p values less than 1e−5 as input for MR. As for the outcomes, we used available outcomes related to age: the mother’s, father’s, and parents’ ages at death (continuous traits) and whether either of the parents’ survival is in the top 1% survival compared with all other parents (binary trait).

The MR analysis yielded 11 significant associations (MR p < 0.05), with the strongest being in the liver (Figure 6B; Figure S6; Table S4). For TFs, ZNF518B and CXXC1, liver expression shows a positive MR effect on mother’s age at death; that is, individuals with genetic propensity to express higher levels of these TFs in the liver may live longer or age more slowly (Figure 6C; Figure S6). However, it is noteworthy that liver ZNF518B expression decreases with age, whereas CXXC1 expression increases (Figure 5A). Combining these findings, we hypothesize that ZNF518B may have an anti-aging function, while CXXC1 may have be pro-aging effect. As for the top TF targets, the strongest link is also in the liver, with the BBC3/PUMA gene. BBC3/PUMA is a predicted target of ZIC1 and CXXC1 and has a negative MR effect direction; that is, a genetic predisposition for an increase in expression correlates with a decrease in mother’s age at death (Figures 6B and 6D). In the heart and muscle, we uncover relatively weaker links between the outcomes and EBF1 in the muscle and ETS2, JAG2, MDK, and NPDC1 in the heart. Taken together, MR analysis further identifies some candidate aging-related genes and TFs that may drive differences in longevity in human populations.

Discussion

Our main question was whether common regulatory players underlie the seemingly tissue- and species-specific molecular footprint of aging. Using our data as well as external datasets, we show that although the molecular footprint of aging evolves differently across tissues, striking similarities emerge in terms of affected pathways and underlying regulators. For instance, the liver’s aging footprint is dominated by changes in transcriptome and DNA methylome. In contrast, transcriptomes of heart and quadriceps are relatively stable but have marked changes in histone modification profiles around genes. Despite all these differences, similar pathways are affected in these distinct layers.

One possible explanation of this tissue-specific histone modification differences lies in these tissues’ replicative potential and their distinct rates of age-dependent replacement of histone subunits H3.1/2 with H3.3 (Tvardovskiy et al., 2017). Replication-dependent histones H3.1/2 are gradually replaced by replication-independent histone H3.3, which may alter the epigenetic landscape of aged cells. The liver H3 pool is dominated by the H3.3 isoform by the age of 10 months, with most differences occurring between 5 and 10 months. However, heart changes are more gradual and occur between 5 and 24 months (Tvardovskiy et al., 2017). The time points in our study (6 and 24 months) therefore correspond to different states of histone isoform abundances in these tissues. By 6 months of age, the liver’s histone H3 pool is dominated by H3.3 and therefore does not increase substantially by 24 months (Tvardovskiy et al., 2017). However, heart, and likely muscle, have different age-dependent H3 isoform compositions. The global increase in active marks in heart and muscle may be due to the increase of H3.3, an imprint known to be associated with active chromatin marks (McKittrick et al., 2004; Tvardovskiy et al., 2017).

The striking similarity in TF enrichment between different mouse and human tissues implies that there may be a common and perhaps restricted set of TFs underlying the aging footprint across tissues and species. The ZIC1 motif is highly enriched across multiple tissues and gene-regulatory layers. ZIC1 increases with age in many peripheral human tissues. Although its relationship or possible implication in aging has not been studied extensively, it has been shown that its brown adipose tissue (BAT) expression increases with age and body mass index, concurrently with a decrease in BAT activity (Nascimento et al., 2018). In addition, ZIC1 and ZIC2 transactivate apolipoprotein E (APOE) expression (Salero et al., 2001), one of the strongest human longevity determinants (Bien-Ly et al., 2012; Corder et al., 1993; Gottschalk et al., 2016; McDaid et al., 2017). The facts that APOE increases with age (Figure S5B) and that higher APOE levels correlate with negative outcomes in age-related diseases such as Alzheimer’s disease (Gottschalk et al., 2016) render ZIC1 a prime candidate driver of gene regulatory changes associated with aging.

We also identify other TFs, such as HMGA1, TBP, and CXXC1, as candidate regulators of the aging process. HMGA1 has been linked to mitochondrial function, repair, and maintenance (Dement et al., 2007; Li et al., 2018; Mao et al., 2009) and is implicated in promoting senescence-associated heterochromatic foci, which are associated with transcriptional repression (Narita et al., 2006). In addition, it has recently been shown to promote the senescence-associated secretory phenotype (SASP) through its effect on NAD+ metabolism (Nacarelli et al., 2019). The TATA box binding protein (TBP) motif is enriched in genes that decrease with age. Although this TF has not been directly linked to aging, it can harbor variations in polyglutamine repeats, which may be relevant in age-related processes such as neuro-muscular degenerative disease (Huang et al., 2015; Reid et al., 2004; Wu et al., 2005). CXXC1, or Cfp1, is a member of the Setd1 H3K4 methyltransferase complex and binds non-methylated DNA of transcriptionally permissive promoters. Given the trend for hypermethylation with age, CXXC1 binding to many promoters may be affected, which may lead to differences in H3K4 methylation. CXXC1 may therefore be an important link between the different molecular layers, which merits further mechanistic investigation, especially given that its motif’s enrichment varies in direction in different tissues, suggesting a complex context-dependent relationship with aging.

The expression levels of some of these TFs vary with age in different species, and in some cases, we can find genetic evidence of their possible involvement in determining variation in human aging and longevity. MR is a powerful tool to infer causality in human GWAS data. However, as the number of participants in GTEx is small compared with the UK Biobank studies (n = 150–491 in GTEx versus up to 500,000 in UK Biobank), the results must be treated with caution and may not necessarily imply causality. On an exploratory level, MR helps us focus on interesting aging candidates that are not only implicated in the process but also explain variation in human longevity. Although the MR outcomes in our study concern the parents of the studied individuals, parental longevity has been used successfully as a quantitative trait reflecting rates of aging (McDaid et al., 2017; Pilling et al., 2016, 2017). Our combined analyses point to CXXC1 and ZNF518B as not only mediators of aging-related genes but also as genes whose genetically modifiable expression levels may explain variation in human longevity. Importantly, analyses focused on TF expression and therefore do not take TF activity into account, which is regulated at many other levels. Hence, the lack of MR evidence for certain TFs does not diminish their relevance.

One of the top predicted targets of ZIC1 and CXXC1, BBC3 (BCL2 binding component 3; also known as PUMA), is a pro-apoptotic member of the Bcl-2 protein family that induces mitochondrial membrane permeabilization through p53-dependent and p53-independent signals (Han et al., 2001; Yu et al., 2001). There has been a recent surge in interest in senolytic therapies to combat human age-associated disease (van Deursen, 2019). One such strategy is the inhibition of anti-apoptotic Bcl-2 family genes, which results in the clearance of senescent cells and improved health span in mice (Zhu et al., 2016). From that perspective, an increase in BBC3 expression is expected to be beneficial, yet our data point to a negative relationship with longevity. Thus, the increased expression of BBC3 throughout life may affect other aspects of metabolism. For instance, BBC3/PUMA is highly expressed in hepatocellular carcinoma and reduces mitochondrial pyruvate uptake and oxidative phosphorylation, effectively driving the cancer metabolic switch (Kim et al., 2019).

In conclusion, integrated multi-layer, multi-tissue, and multi-species study serves as a hypothesis-generating resource of candidate aging regulators. Future genetic and experimental studies that bootstrap from these findings can provide a more mechanistic understanding of aging or age-related diseases and pave the way for new strategies to cope with their negative health consequences.

STAR★Methods

Key Resources Table

Resource Availability

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Johan Auwerx (admin.auwerx@epfl.ch).

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

Datasets were deposited in NCBI GEO: GSE120290) for transcriptome data, NCBI SRA: SRP162353) for DNA methylation data, NCBI SRA: SRP162386) for histone modification ChIP-seq data.

Essential scripts used in this study are available at Mendeley Data: https://doi.org/10.17632/s5p638kbws.1.

Experimental Model and Subject Details

Animals and Tissue Collection

The tissues used in this study come from an existing biobank from a previous study (Houtkooper et al., 2011). Mouse husbandry and tissue collection was performed at the Ecole Polytechnique Fédérale de Lausanne (EPFL) Center of PhenoGenomics as described previously. In that study, all animal experiments were performed according to Swiss ethical guidelines and approved by the local animal experimentation committee of the Canton de Vaud under license 2172. Briefly, male C57BL/6J mice of either mature adults (13 weeks old) or old (93 weeks old) were purchased from Janvier (St. Berthevin, France). Mice were group housed and received standard chow (#2018, containing 18% protein, 50% carbohydrate and 6.0% fat; Harlan Laboratories, Madison, WI, USA). After reaching the age of 24 (~6 months) or 103 weeks (~24 months), mice were sacrificed after overnight fasting. Tissues were frozen in liquid nitrogen for biochemical and molecular analyses. As we could not do all extractions required for the omic analysis on the same tissue, different mice were used for transcriptome, DNA methylation, and histone modification analysis. The three ChIP-sequencing analyses were performed on samples from the same individuals.

Method Details

Transcriptome Analysis

Total RNA was isolated using Trizol (Life Technologies) and purified using the RNeasy Mini Kit (QIAGEN). It was then assessed for degradation using an Advanced Analytical Agilent Fragment Analyzer. For each condition, the 3 best quality samples were selected for liver and muscle and 4 for heart and taken further for microarray analysis. Microarray analysis was performed using the Affymetrix MoGene 1.0ST for liver and quadriceps and the similar but slightly enhanced Affymetrix MTA 1.0 (an array with probes that target splice junctions that is also known as Clariom D) for heart. Microarray data were normalized using the rma function from the R Oligo package (Carvalho and Irizarry, 2010). Differential expression was performed using the limma package (Smyth, 2005).

DNA Methylation Analysis

We isolated genomic DNA from heart, quadriceps, and liver in three replicates using the QIAamp DNA Mini Kit (QIAGEN). To generate sequencing library, 3 μg of DNA was digested with the methyl insensitive enzyme MspI (NEB) at the CCGG site and then purified with QIAquick PCR Purification Kit (QIAGEN). The purified DNA was mixed with End Repair Mix, and incubated at 20°C for 30 min. The end-repaired DNA was purified and mixed with A-Tailing Mix following by incubation at 37°C for 30 min. The purified Adenylate 3′Ends DNA, Methylated Adaptor and Ligation Mix were combined and the ligation reaction was incubated at 20°C for 15 min. We ran a 2.5% agarose gel to select a narrow size-range of 150bp-400bp and purify the gel with QIAquick Gel Extraction kit (QIAGEN). The Methylation-Gold kit (ZYMO) was used for bisulfite conversion, followed by PCR and target fragment recovery from 2.5% agarose gel.

We performed library quantification and quality control using the Agilent 2100 bioanalyzer instrument (Agilent DNA 1000 Reagents. Paired end sequencing was then performed using the Illumina Truseq (TruSeq PE Cluster Kit V3–cBot–HS,Illumina) on the HiSeq 2000 System (TruSeq SBS KIT-HS V3,Illumina), with read length of 100bp to obtain a total output of 3Gb.

We aligned the sequencing reads and called methylation with Bismark v0.19.0 (Krueger and Andrews, 2011) and Bowtie2 v2.3.4.3 (Langmead and Salzberg, 2012), using default parameters and the Ensembl mm10 genome. We used methylKit v1.8 (Akalin et al., 2012) for differential methylation analyses. To obtain gene body or promoter-level summaries of methylation differences, CpG sites with a q value inferior to 0.2 were summarized by each feature and a mean difference in methylation was calculated, as well as a combined p- and q-values using the sum of logs method or Fisher’s method.

Histone Modification Analysis

We performed ChIP-sequencing for 2 positive histone marks (H3K4me3 - ab8580, Abcam; H3K27ac - ab4729, Abcam) and a negative histone mark (H3K27me3 with ab07-449, Millipore) in three replicates using the Auto iDeal ChIP-seq kit (Diagenode). DNA from ChIP was used to prepare barcoded libraries employing the Ion Xpress Plus Fragment Library Kit and Ion Express Barcode Kit from Life Technologies. After library preparation, the barcoded libraries were screened on an Agilent High Sensitivity DNA chip for size distribution. 1μl of each barcoded library were then pooled and sequenced on an Ion Torrent PGM 314 chip. The read counts from the 314 chip were then used to prepare a final equalized pool. The final library pool was quantified by real-time PCR, used to prepare beads, and sequenced using a P1 chip on the Ion Torrent Proton sequencer. The total output is 20 million reads for libraries from point-source modifications of H3K4me3 and H3K27ac, and 40 million for libraries from H3K27me3, which is a broad-source modification and requires more reads.

To remove low quality reads or noise at the 3′ end, we trimmed reads with a cutoff of 250 bp using FASTQ/A Trimmer (hannonlab.cshl.edu/fastx_toolkit, v0.0.13.2). We used Rsubread v1.32.4 (Liao et al., 2013) to align trimmed reads to mouse reference genome mm10 with default parameters. For each tissue, we used the R package GreyListChIP v1.10 to generate a tissue-specific gray list using the merged BAM alignments (Carroll et al., 2014).

For the TSS region-based analyses, we used the csaw package v1.16.1 (Lun and Smyth, 2016) to perform differential binding, excluding reads in gray lists and with a mapping quality score below 20, on regions ± 5kb around gene transcription starting sites (TSS). For the peak-based analysis, we used epic2 for peak calling (Stovner and Sætrom, 2019), followed by diffbind for differential analysis. There is a general agreement between the two methods (Figure S1).

Quantification and Statistical Analysis

Differential Analyses and Gene Set Enrichment Analysis (GSEA)

We used the clusterProfiler R package (Yu et al., 2012) to conduct GSEA analysis on gene ontology biological process terms. We used a minimum gene set size of 10, a maximum gene set size of 500, and performed 10000 permutations. For the transcriptome data, we used a gene list that is ordered by log2(Fold Changes) from the differential expression analysis. In the methylation analysis, we used the mean difference in methylation that is calculated by summarizing all the CpG sites that have q-value less than 0.2 within the gene body or promoter of a gene. The gene list is then ordered based on that value. For the histone modification analyses, we used the log2(Fold Change) values from the differential analyses of reads around TSS sites of genes (see above). Pairwise comparisons between layers was performed by overlapping the genes that are up or down with an FDR adjusted p value of 0.1. The R package SuperExactTest was used to calculate the significance of the overlap.

Differential Enrichment of TF Binding Sites

For enrichments, we used the R PWMenrich package in conjunction with the mouse HOCOMOCO-v10 motifs from the motifDB package. We first constructed a lognormal background distribution using the sequences of ± 5kb region around the transcription starting site (TSS) of all genes. For differential enrichment analysis in each tissue and assay, we selected genes that pass a 10% FDR threshold and separated them into two groups depending on the direction of their aging driven change. For the RRBS upstream data, we took the sequences 5kB, then we calculated the differential enrichment of each TFBS as the difference in enrichment between the genes whose measurement increases with age versus those that decrease. For human data, we used the differential expression results for age (see below) and took the 1000 up and downregulated genes. The reason behind this is the large discrepancy in sample sizes between the three tissues (153, 272, 491 for liver, heart, and muscle respectively), which affects significance levels. We used the HOCOMOCO-v10 motifs for humans. To define the top transcription factors in the mouse data alone, we first took the rank of the absolute value of enrichment and then calculated the mean per tested motif, combining all tissues and layers. To define top TFs across species, we used the same strategy on the enrichment results from the Auwerx, Brunet (Benayoun et al., 2019), and GTEx transcriptome data alone. Since the number of motifs in the human and mouse HOCOMOCO databases is not equal, we further normalized this rank. For defining top targets of top TFs, we defined a potential target as any gene that has a background corrected p value lower than 0.05, then we counted the frequency of these genes (how many times they are assigned to an enriched TF) in the mouse transcriptome data as well as in the human data. We retained mouse targets that appear at least in one layer in all three tissues. We retained human targets in at least 2 of the three tissues. The intersection of these two lists yielded a final list of 23 targets that were used for Mendelian randomization. While the HOCOMOCO-v10 for SREBF1/2 is SRBP1/2, the main figures and panels use SREBF1/2.

Analysis of External Datasets

GTEx Data

We used the GTEx v7 gene-level transcript per million data along with genotypes and phenotypes of the subjects (dbGAP approved request #10143 - AgingX). The GTEx consortium provide the covariates used for their eQTL analyses and they contain known (sex, genotyping principal components, and sequencing platform) and unknown factors estimated using probabilistic estimation of expression residuals (Stegle et al., 2012). The PEER factors may remove the effect of age, and in fact, the first PEER factor significantly correlates with age. We therefore re-estimated PEER factors using age in addition to the existing known covariates (sex, platform, and the three genotyping principal components). Then we performed differential expression analysis using voom and limma in order to obtain the coefficients of age in a model that include all the other covariates. For re-calculating genome-wide eQTLs, we used the same pipeline as in GTEx v7 with two exceptions. (1) we used our own estimated covariates that take age into account and (2) we used a linear model in R for the eQTL calculation. Associations with p values < 1e-5 were kept for the subsequent Mendelian Randomization analysis.

Mendelian Randomization

We used the TwoSampleMR package version 0.4.22 (Hemani et al., 2018). We used four outcomes that were included in the package’s database: Top 1% survival (Pilling et al., 2016) (id:1091), Parents’ age at death (id:1094), Father’s age at death (id:UKB-b:11303), Mother’s age at death (id:UKB-b:12687). We performed clumping to prune SNPs with high LD (r2 cutoff of 0.001) and followed the standard pipeline for MR analysis using Inverse Variance Weighted regression. Scatter and forest plots were obtained using functions from the same package.

Other Data

For the Human Liver Cohort dataset, we used Synapse.org to download the curated expression and phenotype data from the syn88644 dataset after obtaining appropriate permissions from the lead author (Schadt et al., 2008). The skeletal muscle data were obtained from ArrayExpress accession E-MTAB-1788 (Su et al., 2015). For the RNA-seq dataset from the Brunet lab, we obtained the supplemental data and code from the published article (Benayoun et al., 2019).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Rabbit polyclonal to H3K4me3	Abcam	ab8580, RRID: AB_306649
Rabbit polyclonal to H3K27ac3	Abcam	ab4729, RRID: AB_2118291
Rabbit polyclonal to H3K27me3	Millipore	ab07-449, RRID:AB_310624
Chemicals, Peptides, and Recombinant Proteins
TRIzol	Life technologies	15596026
EZ DNA Methylation-Gold	Zymo Research	D5005
16% Formaldehyde, methanol free	Perbio Science	28908
RNase cocktail	Ambion	AM2286
Critical Commercial Assays
Auto iDeal ChIP-seq kit for Histones	Diagenode	C01010171
GeneChip MoGene 1.0ST Array	Affymetrix	901171
GeneChip MTA 1.0	Affymetrix	902512
RNeasy Mini Kit	QIAGEN	74106
QIAamp DNA Mini Kit	QIAGEN	51304
QIAquick PCR Purification Kit	QIAGEN	28104
QIAquick Gel Extraction Kit	QIAGEN	28706
TruSeq PE Cluster Kit	Illumina	V3–cBot–HS
Qubit dsDNA HS assay kit	Life technologies	Q32854
Ion Xpress Plus Fragment Library Kit	Life Technologies	4471269
Ion Express Barcode Adapters 33-48 Kit	Life Technologies	4474518
Ion Torrent PGM 314 chip	Life Technologies	4482261
Ion Torrent PI Chip	Life Technologies	A26771
Deposited Data
Transcriptome data	NCBI GEO	GSE120290
DNA methylation data (RRBS)	NCBI SRA	SRP162353
Histone modification ChIP-seq data	NCBI SRA	SRP162386
Analysis code	Mendeley Data	https://doi.org/10.17632/s5p638kbws.1
Experimental Models: Organisms/Strains
Mouse: C57BL/6J	Janvier	NA
Software and Algorithms
FASTQ/A Trimmer v0.0.13.2	NA	http://hannonlab.cshl.edu/fastx_toolkit/
Rsubread v1.32.4	(Liao et al., 2019)	https://bioconductor.org/packages/release/bioc/html/Rsubread.html
Csaw v1.16.1	(Lun and Smyth, 2016)	https://bioconductor.org/packages/release/bioc/html/csaw.html
Genomation v1.14.0	(Akalin et al., 2015)	https://bioconductor.org/packages/release/bioc/html/genomation.html
clusterProfiler v3.10.0	(Yu et al., 2012)	https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
TwoSampleMR v0.4.22	(Hemani et al., 2018)	https://github.com/MRCIEU/TwoSampleMR
epic2	(Stovner and Sætrom, 2019)	https://github.com/biocore-ntnu/epic2
Other
GTEx data	NCBI dbGaP approved request #10143 – AgingX (GTEx Consortium, 2013)	phs000424.v7.p2
Human Liver Cohort data	https://www.synapse.org (Schadt et al., 2008)	syn88644
Human skeletal muscle expression data	ArrayExpress (Su et al., 2015)	E-MTAB-1788
External Mouse RNaseq data	(Benayoun et al., 2019)	https://github.com/BenayounLaboratory/Mouse_Aging_Epigenomics_2018