Polycomb repressive complex 2 epigenomic signature defines age-associated hypermethylation and gene expression changes

Abstract

Although age-associated gene expression and methylation changes have been reported throughout the literature, the unifying epigenomic principles of aging remain poorly understood. Recent explosion in availability and resolution of functional/regulatory genome annotation data (epigenomic data), such as that provided by the ENCODE and Roadmap Epigenomics projects, provides an opportunity for the identification of epigenomic mechanisms potentially altered by age-associated differentially methylated regions (aDMRs) and regulatory signatures in the promoters of age-associated genes (aGENs). In this study we found that aDMRs and aGENs identified in multiple independent studies share a common Polycomb Repressive Complex 2 signature marked by EZH2, SUZ12, CTCF binding sites, repressive H3K27me3, and activating H3K4me1 histone modification marks, and a “poised promoter” chromatin state. This signature is depleted in RNA Polymerase II-associated transcription factor binding sites, activating H3K79me2, H3K36me3, H3K27ac marks, and an “active promoter” chromatin state. The PRC2 signature was shown to be generally stable across cell types. When considering the directionality of methylation changes, we found the PRC2 signature to be associated with aDMRs hypermethylated with age, while hypomethylated aDMRs were associated with enhancers. In contrast, aGENs were associated with the PRC2 signature independently of the directionality of gene expression changes. In this study we demonstrate that the PRC2 signature is the common epigenomic context of genomic regions associated with hypermethylation and gene expression changes in aging.

Keywords:

• aging
• epigenetics
• epigenomics
• polycomb
• PRC2
• ENCODE
• GenomeRunner
• methylation

Abbreviations

aDMR	=	age-associated differentially methylated region
aGEN	=	promoter of an age-associated gene
ENCODE	=	Encyclopedia of DNA elements
PRC2	=	Polycomb repressive complex 2
TFBS	=	transcription factor binding site

Introduction

Aging affects all aspects of cellular regulation and maintenance, ultimately leading to detrimental physiological effects, impacting the well being of the whole organism. Several reports have identified age-associated genes (aGENs) whose gene expression patterns change with and predictive of age.^1-4 In the last decade, age-associated epigenetic changes attracted significant attention,^5,6 with DNA methylation being the best studied.^1,7-17 In human subjects, age-related changes in DNA methylation have been detected mostly in whole blood,^1,14,17 in purified subsets of blood,^11,14 and in the brain.^18,19 The relationship of age-related DNA methylation regions (aDMRs) with CpG islands,⁹ “bivalent domains,”¹⁴ and gene-centric regions, such as promoters, exonic, intronic and intergenic regions,¹¹ Alu, LINE-1, and other repetitive elements^11,20,21 has been extensively investigated. Yet, the unifying regulatory epigenomic mechanisms that may be targeted to neutralize and, potentially, reverse the detrimental effect of aging remain poorly understood.

Recent years have seen rapid growth of genomic organization and functional annotation data. The Encyclopedia of DNA Elements (ENCODE)²² and the NIH Roadmap Epigenomics²³ projects have aimed to identify the functional elements in the human genome across the major cell types and tissues. We refer to these functional and regulatory regions as epigenomic elements, i.e., genomic properties other than the DNA sequence that describes functions, properties, or experimental values associated with specific genomic regions.^24-27 Although the definition of “epigenomics” is still debated,²⁸ we feel a more inclusive definition of epigenomic data is better suited to convey the concepts of this study.

Epigenomic marks are generally tissue- and cell-type specific,^23,29,30 suggesting that age-dependent epigenomic changes may also be cell- or tissue-specific.³¹ However, several studies have shown that methylation changes can be defined independently of sex, tissue type, cellular composition, differentiation status, and array platform.^{7,10,12,14,17} In this study we used the ENCODE and Roadmap Epigenomics data to investigate potential cell- and tissue-type specificity of epigenomic signatures of aDMRs and aGENs.

It is generally thought that DNA methylation decreases with age.^10,11,20 Yet, some promoter-specific CpG dinucleotides have been shown to become hypermethylated with age,^9,10,14,17 which is a trend also associated with cancer progression.³² To investigate whether the directionality of age-associated changes is associated with similar epigenomic signature(s), we took advantage of published lists of aDMRs hyper- or hypo-methylated with age, as well as genes that increase or decrease their expression with age.

aDMRs have been linked to the promoters of key developmental genes, suggesting the former may directly affect gene expression.^14,16-18 However, methylation and gene expression changes are not necessarily correlated.^1,15,16 We investigated whether similar epigenomic signature(s) are shared by aDMRs that are associated vs. not associated with gene expression changes.

We found epigenomic marks from human embryonic stem cell lines most frequently and most significantly associated with aDMRs and aGENs, concordant with the hypothesis that age-associated methylation and gene expression changes are likely most relevant to developmental and differentiation processes. Both aDMRs and aGENs shared a similar combination of repressive (H3K27me3) and activating (H3K4me1) histone modification marks, as well as binding of Polycomb Repressive Complex 2 (PRC2) co-factors. Consequently, aDMRs and aGENs were enriched in chromatin regions labeled as “poised promoters.” On the contrary, binding sites of RNA Polymerase II co-transcription activators were absent in aDMRs and aGENs, concordant with the depletion of the “active promoter” chromatin state. This epigenomic signature was generally stable across cell lines, with the exception of monocyte-specific aDMRs¹⁵ and aGENs identified by de Magalhães et al.³³ in a meta-analysis. Promoters of genes that increased/decreased their expression with age were associated with the same PRC2 signature. In contrast, aDMRs that were hyper- and hypo-methylated with age were associated with different epigenomic contexts, with hypomethylated aDMRs being enriched in enhancers. Prioritizing all genes by the presence of PRC2 signature identified that age-associated processes most often affect organ formation and molecular communications among cells.

Results

Age-associated methylated regions and genes are largely non-overlapping

Sets of age-associated methylated regions (aDMRs) and promoters of age-associated genes (referred hereafter as aGENs for brevity) were obtained from supplementary material of 8 studies (Table 1). The majority of the studies used the Illumina Infinium HumanMethylation450 BeadChip³⁴ to identify aDMRs, with only one study using the Illumina HumanMethylation27 BeadChip.¹⁴ Genomic coordinates of aDMRs and aGENs were extracted from the platform-specific annotation files (Methods). Where information was available, aDMRs and aGENs were separated into sets positively and negatively correlated with age by methylation and expression, respectively. The aDMRs and aGENs from different studies were largely non-overlapping (Tables S1 and S2, Fig. S1, the fact also discussed in¹⁶), with only one aDMR being identified in 7 studies and 335 out of 16,854 aDMRs detected by 3 studies. These results suggest that different studies may point toward study-specific age-related epigenomic associations.

Table 1 Datasets used in the current study

Abbreviation	Brief description
Age-associated differentially methylated regions (aDMRs)
Hannum_aDMRs_primary ^1	71 methylation markers predictive of age, primary cohort, whole blood
Hannum_aDMRs_all ^1	89 methylation markers predictive of age, primary and validation cohorts, whole blood
Alisch_aDMRs_all; ^7 Alisch_aDMRs_all_pos; Alisch_aDMRs_all_neg	2,078 regions correlated with age by methylation in peripheral blood; 479 of them showed positive correlation with age by methylation, 1,599 were found to be negatively correlated with age
Rakyan_aDMRs_all ^14	131 regions correlated with age by methylation in whole blood, and show the same directional age-associated DNA methylation change in CD4+ T-cells and CD14+ monocytes.
Reynolds_aDMRs_tcells; ^15 Reynolds_aDMRs_tcells_pos; Reynolds_aDMRs_tcells_neg	2,595 regions correlated with age by methylation in CD4+ T-cells. Two,049 of them showed positive correlation with age by methylation, 546 were found to be negatively correlated with age
Reynolds_aDMRs_monocytes; ^15 Reynolds_aDMRs_monocytes_pos; Reynolds_aDMRs_monocytes_neg	2,259 regions correlated with age by methylation on CD14+ monocytes. 468 of them showed positive correlation with age by methylation, 1,791 were found to be negatively correlated with age
Steegenga_aDMRs_and_gene ^16	726 methylation markers correlated with age- and gene expression changes in peripheral blood
Steegenga_aDMRs_not_gene ^16	4,552 methylation markers correlated with age- but not with gene expression changes in peripheral blood
Steegenga_aDMRs_meta ^16	7,477 age-associated methylation markers identified in multiple studies
Heyn_aDMRs ^11	3,205 age-associated methylation markers differentially methylated in cord blood of newborns and CD4+ T cells of centenarian; 1,219 of them showed hypermethylation with age, 1,986 were found to be hypomethylated with age
Age-associated genes (aGENs)
Hannum_genes_associated ^1	326 genes correlated with age by expression, whole blood, data not related to methylation data
Hannum_genes_predictive; ^1 Hannum_genes_predictive_pos; Hannum_genes_predictive_neg	55 genes predictive of age based by their expression, whole blood, data not related to methylation data. Twenty-four of these genes showed positive correlation with age by expression, 31 were found to be negatively correlated with age
Genes_LongevityMap ^47	752 manually curated list of age-associated genes
deMagalhaes_genes_all; ^33 deMagalhaes_genes_pos; deMagalhaes_genes_neg	73 age-associated genes. 56 genes overexpressed with age; 17 genes underexpressed with age.

Epigenomic marks from human embryonic cell lines are most frequently and significantly associated with age-associated regions

The ENCODE and the Roadmap Epigenomic projects provide a wealth of cell type-specific genome annotation data. A total of 1,951 transcription factor binding site (TFBS) datasets (277 transcription factors, 91 cell lines), 722 histone mark data sets (41 marks, 69 cell lines), 136 chromatin state datasets (15 states, 9 cell lines) provided by the ENCODE project, and 979 histone marks data sets (127 cell- and tissue-types) provided by the Roadmap Epigenomics project were used in the analyses (Table S3).

To identify cell types with epigenomic marks most frequently and significantly enriched in aDMRs and aGENs, cell type-specific enrichment frequencies were compared with what can be observed by chance using Fisher's exact test (see Methods, Cell type enrichment analysis). Epigenomic marks from H7es cells, the undifferentiated human embryonic stem cell line, were most frequently and most significantly associated with aDMRs and aGENs (P = 3.94E-12). Other cell lines showing frequent epigenomic associations include K562 leukemia (P = 1.53E-06), HelaS3 cervical carcinoma (P = 3.10E-06), and H1hesc embryonic stem cell lines (P = 1.37E-05, Table S4A). These observations were confirmed by the cell type enrichment analysis of 127 cell- and tissue type-specific epigenomic datasets provided by the Roadmap Epigenomics project, with epigenomic marks from “ESC.H1,” H1 embryonic stem cell line, being most frequently and significantly associated with aDMRs and aGENs (P = 2.60E-12, Table S4B).

The epigenomic similarity analysis among cell types (see Methods, Epigenomic similarity among cell types) identified H1hesc and K562 cell lines as most epigenomically distinct from the other cell types (Fig. 1). The epigenomic similarity analysis of 127 cell- and tissue types used by the Roadmap Epigenomics project further highlighted a group of embryonic cell lines epigenomically distinct from the other cell types (Fig. S2). In summary, these results suggest preferential association of aDMRs and aGENs with epigenomic marks in cells with high developmental potential, such as embryonic stem cells.

Figure 1. Epigenomic similarity among the ENCODE cell lines. We computed the Spearman correlation coefficient between cell type-specific epigenomic enrichment profiles (Methods). The resulting correlation matrix was clustered using Euclidean/average clustering metrics, and visualized with darker blue/red gradient representing weaker/stronger epigenomic similarity, respectively. Each cell shows the numerical value of the corresponding Spearman correlation coefficient.

aDMRs and aGENs share a similar repressive epigenomic signature

To identify whether aDMRs and aGENs from multiple studies were enriched in similar epigenomic signatures, their epigenomic similarity was assessed across all data sets (see Methods, Epigenomic similarity among aDMRs and aGENs, and³⁵). Briefly, each set of aDMRs and aGENs was tested for enrichment in histone modification marks, TFBSs, and chromatin states obtained from H1hesc human embryonic cell line. The epigenomic data from H7es cell line was not used, as it provides insufficient data to evaluate epigenomic similarity, covering only 3 histone modification marks, H3K27ME3, H3K4ME3, and H3K36ME3. The aDMR- and aGEN set-specific enrichment P-values (epigenomic enrichment profiles) were collected, −log₁₀-transformed, and correlated with each other using Spearman correlation coefficient (Fig. 2). Notably, the aGENs identified by a meta-analysis³³ were the most epigenomically distinct. While these results suggest that the epigenomic signatures may be driven by study-specific conditions, the aDMRs and aGENs from different studies showed on average high epigenomic similarity (median Spearman correlation coefficient 0.45, median P -value < 1.00E-8). There results were confirmed by the epigenomic similarity analysis using the ENCODE and Roadmap Epigenomics data from all cell lines (Fig. S3A and S3B), showing even higher overall similarity among the aDMRs and aGENs from multiple studies (median Spearman correlation coefficient ranging from 0.51 to 0.55, median P -value < 1.00E-8). These results indicate that the largely non-overlapping sets of aDMRs and aGENs share similar epigenomic enrichments, suggesting the presence of a common epigenomic signature associated with them.

Figure 2. Epigenomic similarity among aDMRs and aGENs using the ENCODE data, H1hesc cell line. We computed the Spearman correlation coefficient between aDMR-and aGEN-specific epigenomic enrichment profiles (Methods). The resulting correlation matrix was clustered using Euclidean/average clustering metrics, and visualized with darker blue/red gradient representing weaker/stronger epigenomic similarity, respectively. Each cell shows the numerical value of the corresponding Spearman correlation coefficient.

Age-associated regions are concurrently enriched in repressive (H3K27me3, EZH2) and activating (H3K4me1, H3K4me2) histone modification marks

The H3K27me3 repressive mark was the most frequent and most significantly enriched in the aDMRs and aGENs (Fig. 3, Fig. S4, Table S5A and S5B). It was accompanied by enrichment in histone methyltransferase EZH2, a subunit of Polycomb Repressive Complex 2 (PRC2), known to trimethylate H3K27me3 and repress transcription.^17,36,37 Binding sites of another PRC2 subunit, SUZ12, and its interacting partner CTCF, along with the H2az histone modification mark known to co-localize with PRC2 at developmentally silenced genes³⁶ were also noted. These results suggest potential association of PRC2 with age-related methylation and gene expression changes.

Figure 3. Histone modification marks (ENCODE) enriched/depleted in aDMRs and aGENs. Darker blue/red gradient highlights depleted/enriched associations, respectively. Red/green bar gradient defines frequency of epigenomic marks enriched/depleted in aDMRs and aGENs, respectively.

Despite being enriched in repressive histone modification marks and TFBSs, the aDMRs and aGENs were also enriched in activating H3K4me1, H3K4me2, and H3K4me3 histone modification marks, known to surround transcription start sites.³⁸ These activating marks have been shown to co-localize with the H3K27me3 repressive mark in “bivalent domains” in embryonic stem cells,³⁹ confirming the functional role of aDMRs in orchestrating developmental regulation through such domains.¹⁴ In summary, these results suggest that binding of the PRC2 complex may be involved in regulation of aGENs poised for active transcription, and that its regulatory effect may be altered by the presence of aDMRs.

H3K27ac, H3K79me3, and H3K36me3 marks of active transcription are depleted in age-associated regions

H3K27ac, a mark of active promoters and enhancers, H3K36me3, a mark elevated predominantly in transcribed regions of active genes and peaked at the 3′ end of genes,⁴⁰ and H3K79me3, a promoter mark, were all depleted in the aDMRs and aGENs (Fig. 3, Fig. S4, Table S5A and S5B). Depletion of these marks in aGENs was the most pronounced, as compared with aDMRs. These findings are consistent with the observation that H3K27ac is enriched in enhancers poised for active transcription and depleted in active enhancers,⁴¹ further supporting the association of aDMRs and aGENs with “bivalent domains.”

Repressive transcription factor binding sites are enriched in age-associated regions

EZH2 and SUZ12 transcriptional repressors and histone modifiers were identified among the most frequent binding sites enriched in aDMRs and aGENs (Fig. 4, Table S5C). Two other transcriptional repressors, SIN3A and CtBP2, were also highly enriched. Another enriched transcription factor, E2F6, which plays a crucial role in controlling cell cycle and tumor suppressor proteins, also possesses repressive activity upon E2F-dependent transcription. Enrichment of aDMRs and aGENs in Rad21, a double-strand-break repair protein and a part of cohesin complex interacting with PRC2 complex and participating in controlling transcription of silenced and active genes⁴² was also noted. These results emphasize enrichment of aDMRs and aGENs in binding sites of other repressive regulators known to interact with PRC2.

Figure 4. Transcription factor binding sites (ENCODE) enriched/depleted in aDMRs and aGENs. Darker blue/red gradient highlights depleted/enriched associations, respectively. Red/green bar gradient defines frequency of epigenomic marks enriched/depleted in aDMRs and aGENs, respectively.

RNA Polymerase II and other transcription initiation factors are depleted in age-associated regions

Concordant with the presence of the PRC2 signature, aDMRs and aGENs were depleted in RNA Polymerase II binding sites, and sites of other proteins necessary for RNA Polymerase II-dependent transcription, such as TAF1, TAF7, TBP, CREB1 (Fig. 4, Table S5C), and a group of multi-functional transcription factors, such as the activators/repressors SP1 and SP2 and development and oncogenesis regulators SRF, SIX5, JunD, c-Myc. JARID1A, an H3K4 demethylase required for haematopoietic stem cell self-renewal,^43,44 showed marginal enrichment in some, but an overall strong depletion, in the majority of aDMRs and aGENs. These results further strengthen the connection of aDMRs and aGENs with the binding of Polycomb repressive complex 2 in the absence of activating transcription factors.

“Poised promoters” but not “active promoters” are enriched, and “heterochromatin” chromatin states are depleted in age-associated regions

Epigenomic enrichment analysis of chromatin states identified the “poised promoter” state as the most frequently enriched in aDMRs and aGENs, followed by the “strong enhancer” state (Fig. 5, Table S5D). Consistent with the depletion of RNA Polymerase II and activating TFBSs, the “active promoter” chromatin state was among the most depleted. A group of monocyte-specific,¹⁵ pediatric-teenager,⁷ pediatric-centenarian,¹¹ and meta aDMRs¹⁶ showed a weaker enrichment in “poised promoter” but a stronger depletion in “active promoter” chromatin states. The “transcriptional elongation” state was uniformly depleted in all sets of aDMRs and aGENs. However, the “heterochromatin” state was also depleted, suggesting the presence of open chromatin around aDMRs and aGENs. These results further confirm the association of aDMRs and aGENs with poised promoters and enhancer regions controlled by repressive transcription factors and the PRC2 complex.

Figure 5. Chromatin states (ENCODE) enriched/depleted in aDMRs and aGENs. Darker blue/red gradient highlights depleted/enriched associations, respectively. Red/green bar gradient defines frequency of epigenomic marks enriched/depleted in aDMRs and aGENs, respectively.

Repressive signature is associated with aDMRs and aGENs positively correlated with age, while enhancers are preferentially associated with aDMRs hypomethylated with age

Some studies provide the directionality of age-association changes, that is, whether aDMRs are hyper/hypomethylated with age, and whether aGENs are positively/negatively correlated with age by expression. This directionality of changes is referred to hereafter as “positive/negative.” Sets of positive and negative aDMRs and aGENs were clustered by their epigenomic similarity to identify epigenomically distinct groups. This clustering, as expected, separated positive and negative aDMRs into two distinct groups, with the exception of monocyte-specific positive aDMRs¹⁵ (Fig. 6, Fig. S5). Notably, positive and negative aGENs did not show epigenomic differences, sharing the same epigenomic signature with aDMRs hypermethylated with age. These results suggest that aDMRs hyper- or hypo-methylated with age may be associated with differential epigenomic signatures, whereas the same epigenomic mechanisms may drive increasing or decreasing gene expression changes of aGENs.

Figure 6. Epigenomic similarity among aDMRs hyper/hypomethylated with age, and aGENs showing increasing/decreasing gene expression with age (“pos/neg” postfixes). Pos/Neg aDMR-and aGEN epigenomic enrichment profiles obtained using H1hesc epigenomic data from the ENCODE project (see Methods) were correlated to each other using the Spearman correlation coefficient metric. The resulting correlation matrix was clustered using Euclidean/average clustering metrics, and visualized with darker blue/red gradient representing weaker/stronger epigenomic similarity, respectively. Each cell shows the numerical value of the corresponding Spearman correlation coefficient. Red/green bars define groups of aDMRs and aGENs showing high inter-group similarity.

In order to identify differential epigenomic signatures between the groups of aDMRs hyper- or hypomethylated with age, we compared them using differential epigenomic analysis (see Methods). The repressive H3K27me3 mark was not significantly associated with negative aDMRs, while being highly enriched in positive aDMRs. On the contrary, H3K36me3, a mark of transcribed regions, was more significantly depleted in positive aDMRs (Table S6). Several epigenomic marks were identified as differentially enriched, that is, significantly enriched in positive aDMRs and aGENs while depleted in negative aDMRs and aGENs, or vice versa. Among them were repressor EZH2, histone deacetylase HDAC2, lysine-specific demethylase JMJD2A (aka KDM4A),⁴⁵ and a component of methyltransferase complex, RBBP5⁴⁶ (Table S6). These results highlight clear differences in histone code and histone modifiers associated with positive and negative aDMRs.

Twenty-three transcription factors also showed differential directionality and significance of enrichment. The repressive EZH2, SUZ12, CtBP2, E2F6, Znf143, SIN3A, YY1, and Max TFBSs were significantly enriched in positive aDMRs, and depleted in negative aDMRs. The transcription-promoting TAF1, TBP and RNA Polymerase II factors, on the other hand, were significantly depleted in the negative aDMRs. Comparison of chromatin states identified “strong enhancer” and “weak enhancer” as enriched in negative aDMRs. Positive aDMRs, on the other hand, were significantly depleted in “weak transcription” chromatin states (Supplementary Table S6). These results emphasize the PRC2 epigenomic signature as associated with aDMRs hypermethylated with age, while linking hypomethylated aDMRs with enhancers.

Study-specific epigenomic associations are generally stable across cell lines

In order to identify potential cell type-specific enrichments in each study, sets of aDMRs and aGENs were analyzed individually (Supplementary Results S1). Different sets of aDMRs and aGENs emphasized the aforementioned histone modification marks differentially (Supplementary Results S1A, D), but consistently enriched/depleted across cell types. For example, aDMRs predictive of age¹ did not show depletions, while being enriched in EZH2, CTCF, H3K27me3 and other repressive and activating histone marks. Other aDMRs and aGENs showed consistent enrichment and depletion patterns, with enrichments in H1hesc epigenomic data among the most significant. The pediatric-teenager aDMRs⁷ showed the strongest association with the H2az histone mark. The monocyte-specific aDMRs,¹⁵ meta-aDMRs¹⁶ and pediatric-centenarian aDMRs¹¹ were the most prominently enriched in H3K4me1 across the majority of cell types. These results suggest that the age-associated histone code may be similarly regulated in any cell type.

Similarity of the enrichment patterns across cell types was also supported by the analysis of TFBSs (Supplementary Results S1B). As described in the previous section, a group of aDMRs, such as pediatric-teenager aDMRs,⁷ and T-cell-specific aDMRs,¹⁵ showed depletion in RNA Polymerase II related TFBSs, while the others^14-16 also showed enrichment in EZH2, SUZ12, CtBP2, RAD21, Znf143, SIN3A, and E2F6. Again, enrichments in H1hesc epigenomic data were most pronounced. In contrast with other sets of aDMRs and aGENs, meta-aGENs from³³ showed enrichment in RNA Polymerase II, TAF1, TBP and other transcription promoting transcription factors. This observation placed TFBSs in a category of regulatory marks most diversely associated with aDMRs and aGENs.

Cell type-specific chromatin states analysis showed consistent enrichment of aDMRs and aGENs in “poised promoter,” “strong enhancer” and “repressed” chromatin states, and depletion in “heterochromatin,” “transcription elongation” and “active promoter” states (Supplementary Results S1C). Several sets of aDMRs and aGENs,^1,14,47 including T-cell-specific aDMRs,¹⁵ were strongly enriched in the “poised promoter” state, while the others,^7,11,16,33 including monocyte-specific aDMRs,¹⁵ were predominantly enriched in “strong enhancer” state. Although meta-aGENs³³ were enriched in “strong enhancer” and depleted in “heterochromatin” chromatin state, comparable with other aGENs and aDMRs, they were also enriched in the “active promoter” states. In summary, each set of aDMRs and aGENs emphasized different aspects of the repressive PRC2 signature that was generally similar across multiple cell types.

Genes with an age-associated epigenomic signature are enriched in organ morphogenesis, cell signaling, extracellular matrix and brain-related functions

In order to identify genes showing the identified age-associated epigenomic signature in their promoters, independent of cell- or tissue specificity, we performed reverse epigenomic analysis of the promoters of all genes (Table S7). Briefly, for each gene, a promoter-specific score was calculated using a set of promoter-specific epigenomic marks weighted by their average significance identified from the enrichment analysis of aDMRs and aGENs (Methods). This analysis was designed to prioritize the promoters of genes having, i.e., H3K27me3 and other repressive marks, TFBSs, and the “poised promoter” chromatin state, and under-emphasize the promoters with activating H3K36me3, H3K27ac, RNA Polymerase II-related TFBSs and the “active promoter” chromatin state. Therefore, it was not surprising to identify “H3K27me3 bound,” “Suz12 targets” and similar gene signatures from MSigDB (P adj. = 5.52E-64/9.86E-48, respectively, Table S8) via functional enrichment analysis on the top 500 genes, prioritized by the age-associated epigenomic signature. 71 of these genes were identified as targets of the hsa-miR-1275 miRNA (P adj. = 9.260E-3), which participates in PRC2-mediated silencing and is regulated by the level of the H3K27me3 repressive mark,⁴⁸ supporting the notion these genes share the age-associated repressive epigenomic signature.

Transmembrane protein 61 (TMEM61) and solute carrier family 10, member 4 (SLC10A4) were among the top genes potentially regulated by the age-associated epigenomic signature (Table S7), suggesting involvement of cell signaling processes. Functional enrichment analysis further identified “gated channel activity” (P adj. = 4.36E-5) molecular function and “intrinsic component of plasma membrane” (P adj. = 3.61E-11) cellular components as the most significant gene ontologies (Table S9). Suggestive of age-related developmental processes, “organ morphogenesis” and “cell-cell signaling” were the most significant biological processes (P adj. = 1.16E-10/1.73E-10, respectively). Of note were enrichments in “extracellular matrix” cellular component functional category (P adj. = 1.41E-8), and “extracellular matrix organization” Reactome canonical pathway (P adj. = 3.97E-5) – processes increasingly recognized in the context of aging.⁴⁹ A separate group of functional enrichments included the “synaptic transmission” biological process (P adj. = 3.56E-8) and the “neuroactive ligand-receptor interaction” KEGG canonical pathway (P adj. = 1.79E-4, Table S10), which suggest involvement of genes regulated by the age-associated epigenomic signature in brain-related processes. This observation was further supported by enrichment in a gene set involved in neurotransmission and neurodevelopment from⁵⁰ (P adj. = 3.968E-5). There results outline core biological processes that may be affected by the PRC2 epigenomic signature in aging.

Discussion

This study identified the epigenomic signature of Polycomb Repressive Complex 2 binding (EZH2, SUZ12, CTCF) to be associated with aDMRs and aGENs across multiple studies. This finding is in line with the notion that PRC2 plays a direct role in modulating longevity and stress resistance.³⁷ The association of “poised promoters,” marked by the repressive H3K27me3, H3K9me3 and activating H3K4me1, H3K4me2, H3K4me3 histone modification marks, is consistent with previous findings of “bivalent domains,”^{11,14,15,23,51} and further establishes PRC2 as a repressor of regulatory regions poised for active transcription. We found clear differences between hyper- and hypo-methylated aDMRs, but not between aGENs positively and negatively correlated with age, linking hypermethylated aDMRs with the PRC2 signature and hypomethylated aDMRs with enhancers. Although epigenomic marks from the human embryonic stem cell line were most strongly associated with aDMRs and aGENs, the PRC2 epigenomic signature was generally consistent across the cell lines. These findings highlight PRC2 as a key regulator of age-related processes and strengthen the notion that it may be used as a target for increasing longevity³⁷ and, perhaps, reversing detrimental effects of aging.

Our selection of normal human embryonic stem cells was guided by the following:1) epigenomic data derived from these cell lines showed the most frequent and most significant enrichments in the age-associated regions, confirmed by the ENCODE and Roadmap Epigenomics data; 2) these cells are pluripotent, i.e., able to differentiate into multiple lineages, thus being relevant to investigation of the epigenomic context of development and aging; 3) epigenomic enrichments across multiple cell lines were generally similar; and 4) DNA methylation changes in genes required for stem cell differentiation have been demonstrated by others.^17,52 Notably, the other cell lines showing frequent enrichments in age-associated regions resembled embryonic cell lines in their developmental properties. For example, the multipotent K562, a highly undifferentiated lymphoblastoid cell line, is able to develop into different blood cell progenitors,⁵³ resembling the developmental potential of haematopoietic stem cells. Furthermore, monoclonal antibodies against K562 cells also recognize human haematopoietic pluripotent stem cells.⁵⁴ Our findings of concordant age-related epigenomic associations across multiple cell types confirm previous studies that reported a high correlation of the age-associated loci in multiple tissues and cell types.^7,14,17

The consistency of the epigenomic signature shared by different sets of aDMRs and aGENs cannot be attributed solely to their overlapping genomic locations. In fact, aDMRs from different studies were largely non-overlapping (Tables S1 and S2), and the number of aDMRs in each study differed. It is important to note that some sets of aGENs were defined independently of aDMRs,^1,33,47 yet, they showed the same epigenomic enrichments as the aDMRs. Note that age-associated regions do not have to be overlapping in order to show similar enrichment patterns, due to the fact that histone modification regions and chromatin states span large stretches of the genome and can overlap with multiple aDMRs and aGENs. If some aDMRs were missed in one study, due to tissue-specific, technical or data processing differences, the enrichment analysis would still identify similar enrichment patterns if they were driven by the same biological effect. Notably, both aDMRs associated and not associated with proximal gene expression changes¹⁶ possessed a similar PRC2 epigenomic signature, suggesting that some aDMRs may exert their regulatory effect via long range interactions, or indirectly affect gene expression by enhancing the effect of other aDMRs.

In human embryonic cell lines, EZH2 and PRC2 have been shown to repress key developmental genes poised for activation during ES cells differentiation,^55,56 regulating pluripotency and differentiation through accumulation of H3K27me3 on chromatin.⁵⁷ However, Polycomb proteins can also facilitate productive gene expression required for establishing the identity and for differentiation of stem cells,^42,58 reviewed in⁵⁹. This activation is accompanied by replacement of H327me3 by H3K27ac⁴¹ and increasing levels of H2az, H3K4me1, and RNA Polymerase II transcription.⁵¹ The predominant enrichment of the H2az histone modification mark in pediatric-teenager aDMRs⁷ suggests that some tissues and developmental stages utilize additional mechanisms together with PRC2 to facilitate active developmental processes. Our study cautions against labeling the PRC2 epigenomic signature as repressive only in context of age-associated changes and warrants further investigation of the role of PRC2 in aging on cell-, tissue- and developmental stage-specific level.

The directionality of age-associated changes deserves further attention. While some studies report predominant DNA hypermethylation¹⁴ or hypomethylation^21,52 with age, little is known about epigenomic mechanisms regulating the increase or decrease of gene expression. Our observation that the directionality of gene expression is not associated with differential epigenomic signatures highlights the dual role of PRC2 as both a repressor and a facilitator of age-associated gene expression changes. Our findings that PRC2 is strongly associated with hypermethylated aDMRs further suggest that PRC2-driven regulation may be affected by the increase of methylation with age.

Although we observed enrichment of PRC2 binding co-factors, and depletion in transcription activation mark H3K36me3 in aDMRs hypermethylated with age, we also noted the enrichment of demethylase JMJD2A and methyltransferase RBBP5. The role of Jumonji domain-containing protein 2A (JMJD2A), which targets methylated H3K36me3 and prevents senescence by modulating p53 and Rb pathways,⁶⁰ has been well studied in cancer.⁴⁵ Its potential role in aging has also been described.⁶¹ The presence of JMJD2A and RBBP5 in the absence of H3K36me3 in aDMRs hypermethylated with age suggests activity of p53 and Rb pathways at early developmental stages and their silencing followed by increased methylation over time.

A set of aGENs identified by a meta-analysis of gene expression profiles from mice, rats and humans³³ was markedly different in its enrichment in RNA Polymerase II and other transcription promoting factors from the overall repressive signature of PRC2 binding. This may be due to pre-selection of genes that consistently change their expression across organisms, hence, emphasizing their RNA Polymerase II-driven transcription. This observation suggests that meta-analysis may emphasize actively transcribed genes while eliminating some important age-related associations.

To extend our finding on gene-centric level, while considering potential limitations of gene expression meta-analysis, we reverse engineered the promoters of all genes in the human genome for presence of the PRC2 signature. In addition to “cell-cell signaling” and “extracellular matrix” processes, increasingly recognized in the context of aging,⁴⁹ we identified several development-related biological processes, ranging from general “organ morphogenesis” and “tissue development” to more specific “sensory organ development” and “neurogenesis” (Table S8). It is important to link these observations with our findings that the PRC2 signature is associated with aDMRs hypermethylated with age. Increase in methylation has traditionally been viewed as leading to decrease in gene expression,^29,32 suggesting these biological processes being downregulated with age. Indeed, prenatal hypomethylation following increased hypermethylation with age has been observed in normal brain development,¹⁹ further suggesting that genes enriched in “neurogenesis,” “synaptic transmission” and other brain-related processes identified in our study as being regulated by the PRC2 signature also undergo repression with age through hypermethylation.

The breadth of cell- and tissue specificity of the ENCODE and the Roadmap Epigenomics data allows unprecedented insights into the epigenomics of aging. More work is required in order to build a comprehensive map of common and cell- and tissue-specific epigenomic signatures enriched in aDMRs and aGENs across the whole human body, such as the Illumina Human BodyMap 2.0⁶² and the Digital Aging Atlas⁶³ projects covering gene expression changes in different tissues. As more age-associated genomic regions are identified, such as the age-associated miRNAs, long non-coding RNAs, epigenome quantitative trait loci from age-oriented epigenome-wide association studies, the methods described in this study will allow an even more comprehensive understanding of the genomics and epigenomics of aging. Last, but not least, identification of hypomethylated aDMRs associated with enhancers, generally regulating gene expression via long-range interactions,⁶⁴ underscores the importance of long-distance interactions of the age-associated genomic regions in the 3D structure of the genome. This understanding will help find a systemic approach to amending the detrimental effects of age in each cell in the human body.

Methods

Data sources and pre-processing

Genomic coordinates of all probes on the Illumina Infinium HumanMethylation27 and HumanMethylation450 BeadChips were extracted from the corresponding probe annotation files (GEO accession numbers GPL8490 and GPL13534, respectively). Unless otherwise specified, hg19 human genome assembly coordinate system was used. The coordinates of the HumanMethylation27 chip were converted from hg18 to hg19 genome assembly using the liftOver tool through the UCSC genome browser.⁶⁵ See https://github.com/mdozmorov/Aging/tree/master/Illumina for details.

Genomic coordinates of aDMRs were obtained from the lists of probes from supplementary material of each study by cross-matching them to the corresponding probe annotation file. Genomic coordinates of aGENs were obtained from the lists of gene names from supplementary material of each study, and genomic coordinates of their promoters were extracted using the refGene table via the UCSC genome browser MySQL database (accessed on January 10, 2015). A promoter was defined as the 2,000 bp region upstream of a gene's transcription start sites. See https://github.com/mdozmorov/Aging/tree/master/data for details.

The human hg19 genome annotation data from the ENCODE project was obtained through the UCSC genome database⁶⁵ (accessed on 10/12/2014), converted to the BED format and stored in a local SQLite database. Three groups of epigenomic data were used: histone modifications, transcription factor binding sites experimentally obtained by ChIP-seq, and chromatin state segmentation by Hidden Markov Models (Table S3). The additional histone modification data from the Roadmap Epigenomics covering 127 cell- and tissue-types⁶⁶ was downloaded from https://www.broadinstitute.org/˜anshul/projects/roadmap/peaks/consolidated/broadPeak/ (accessed on 03/16/2015).

Cell type enrichment analysis

Cell type enrichment analysis was performed to identify cell type specificity of the epigenomic marks showing most frequent and most significant enrichments in the aDMRs and aGENs. Cell type-specific and a total number of epigenomic datasets were compared with cell type-specific and a total number of epigenomic data sets identified as significantly enriched (P-values < 1.00E-10) in aDMRs and aGENs. The P-values of cell types showing most frequent and most significant enrichments that could happen by chance were calculated using Fisher's exact test. See https://github.com/mdozmorov/Aging/blob/master/R.Aging/cell_type_enrichment.R for more details.

Epigenomic enrichment analysis

The epigenomic enrichment analysis was performed using GenomeRunner.^35,67,68 Briefly, the enrichment analysis evaluates whether a set of aDMRs or aGENs is statistically significantly co-localizes with epigenomic datasets. The two-tailed Fisher's exact test was used to calculate enrichments, or depletions, of such co-localization. A ‘background’, or ‘universe’ of all chip-specific CpGs, or the promoters of all genes, was used to evaluate epigenomic enrichments that occur by chance. When testing a set of aDMRs or aGENs for enrichment in multiple epigenomic data sets, the enrichment P-values were corrected for multiple testing using a False Discovery Rate (FDR) calculation. The P-values are then −log₁₀-transformed, and a “−“ sign is added to denote depleted associations. This transformation converts the P-values into a linear scale, where larger/smaller numbers denote more statistically significant enriched/depleted associations, respectively. The following analyses use such transformed P-values, unless explicitly stated otherwise.

Epigenomic similarity analysis

Epigenomic similarity among aDMRs and aGENs

When several sets of aDMRs and aGENs are analyzed, a natural question may be asked – how similar is the epigenomic context associated with these aDMRs and aGENs? To answer this, the epigenomic similarity analysis³⁵ that compares the epigenomic enrichment profiles among the sets of aDMRs and aGENs was designed. An epigenomic enrichment profile is defined as a vector of P-values obtained by testing a set of aDMRs or aGENs for enrichment in several epigenomic datasets. If two sets of aDMRs are co-localized in similar epigenomic data sets, their epigenomic enrichment profiles will be similar. On the other hand, if the 2 sets of aDMRs are enriched in different epigenomic datasets, their epigenomic enrichment profiles will be different.

Epigenomic similarity among cell types

The same approach is used to evaluate epigenomic similarity among the cell types. Instead of collecting aDMR-or aGEN-specific epigenomic enrichment profiles, the epigenomic enrichment profiles for each cell type were assembled and compared using Spearman correlation coefficients.

Differential epigenomic enrichment analysis

If a group of aDMRs and aGENs is epigenomically different from another group, e.g., Figure 6, epigenomic data sets differentially enriched in aDMRs and aGENs between those groups can be identified. For each epigenomic mark, the distributions of −log₁₀-transformed P-values are compared between the groups using a moderated t-statistics, as implemented in the limma R package.⁶⁹ See https://github.com/mdozmorov/Aging/blob/master/R.Aging/GR_ALL.Rmd for more details.

Clustering and Visualization

The –log₁₀-transformed P-values are also used for visualization purposes. For the epigenomic enrichment analysis, the transformed P-values are assembled into an n × m matrix, where n (rows) is the number of aDMRs and aGENs and m (columns) is the number of epigenomic marks. To allow the user to focus on the most significant enrichments, this matrix is filtered to remove epigenomic marks not showing any significant enrichment with at least one set of aDMRs or aGENs. The resulting matrix is then clustered using “Euclidean” distance to measure dissimilarity between rows and columns, and the “average” agglomeration method. The clustered matrix is then visualized using Matlab-like color gradient, with darker blue/red gradient indicating under/overrepresented enrichments, respectively.

A similar visualization strategy is used to visualize the results of the epigenomic similarity analysis. The n × n matrix of Spearman correlation coefficients is clustered and visualized using Matlab-like color gradient. See https://github.com/mdozmorov/Aging/blob/master/R.Aging/GR_ALL.Rmd for more details.

Reverse epigenomic analysis

Genomic coordinates of the promoters of all genes defined in RefSeq database were obtained from the UCSC genome browser database⁶⁵ (accessed on 10/12/2014). The promoters of each gene were annotated for the presence of 26 histone modification marks, 76 TFBSs and 15 chromatin states. Each annotation was up-weighted by the significance of an average −log₁₀-transformed enrichment P-value of a given epigenomic category. Note that if an epigenomic category was identified as depleted, the absence of it was up-weighted. A summary score of such up-weighted annotations was calculated for each gene (see https://github.com/mdozmorov/Aging/blob/master/R.Aging/annotation_enrichment.R for more details).

Functional enrichment analysis

Gene-centric functional and canonical pathway enrichment analysis was performed using ToppGene Suite.⁷⁰

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

The author thanks the members of the Epigenome Data Analysis and Coordination Center and the NIH Roadmap Epigenomic project who have processed and made publicly available the epigenomic data used in this study. The author is very grateful to Yana Bromberg, Jonathan Wren, Cory Giles, Amy Olex for their review and critique of the manuscript.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.