Shifting epigenetic contexts influence regulatory variation and disease risk

Defining chromatin accessibility change, its genomic context, and loci subject to change

To investigate epigenetic changes occurring over the course of post-natal development and aging, we focused on chromatin accessibility, as it reflects the regulatory potential of a genetic locus and can be considered a property of the epigenome which integrates a number of possible epigenetic phenomena (e.g., regulatory factor binding, chromatin remodelling, etc.) [39]. We thus consider regions with altered chromatin accessibility as being indicative of epigenetic modifications or ‘shifts’ in context. As a read-out of accessibility we analyzed DNase-I hypersensitivity datasets acquired from primary human tissue, and obtained fetal/adult sample pairs for eight distinct primary tissues (spleen, lung, muscle, stomach, kidney, brain, heart, and skin; see Supplementary Table 1) [6, 40]. For each tissue at each time-point, called accessible or “open” chromatin regions were consolidated across biological replicates, then further aggregated by tissue and stage (see Supplementary Methods).

We first identified chromatin regions exhibiting recurrent accessibility changes between fetal and adult samples across tissue types (see Supplementary Methods, Figure 1A, 1B and Supplementary Figures 1, 2). We define regions as ‘adult-biased’ if they exhibit increased differential accessibility in adults compared to in fetuses. Conversely, we define regions as ‘fetal-biased’ if they exhibited decreased differential accessibility in adults compared to in fetuses. These ‘pan-tissue’ altered regions were compared to those defined in individual tissue comparisons, showing substantial but not complete overlap (Figure 1C) – suggesting that our approach captures cross-tissue signals of broader developmental changes and not tissue-specific effects. We next explored possible signals of epigenetic aging occurring in the context of fetal to adult changes, by further dividing our adult tissue samples into ‘younger’ and ‘older’ age categories (Methods, Supplementary Figure 3). We then assessed accessibility change between young and old occurring within the ‘adult-biased’ and ‘fetal-biased’ regions defined above (Figure 1D). This approach identified regions for which young-old differences mirrored fetal-adult differences, as well as regions where aging changes appear to counter developmental patterns. We observed a tendency for shared directionality in gains or losses of accessibility; i.e., adult-biased regions tended to also have increased accessibility in older adult samples, while regions losing accessibility in adult samples (i.e., are fetal-biased) continued this trend in older samples (chi-sq test, p < 0.05, Supplementary Table 1). In this text, we therefore refer to regions with greater accessibility in older samples as ‘old-biased’ and regions with lower accessibility in older samples as ‘young-biased’. As described in the Supplementary Information, we considered histone mark and DNA methylation changes, key features of developmental [21, 26, 33] and aging epigenetic changes [3, 20] as additional means to validate the behavior of these region sets (see also Supplementary Table 1).

To gain insights into the roles these region sets have in transcriptional regulation, we next characterized the genomic distribution of our adult-biased and fetal-biased region sets using adult tissue epigenetic states [26] (Supplementary Methods). We found that our region sets preferentially fell within different epigenetic states (e.g., enhancers, heterochromatin) depending on the nature of their accessibility shift (e.g., adult-biased, old-biased), suggesting that these shifts may be associated with altered regulatory biology at different loci, and that the interaction between fetal and adult shifts as well as young and old-age shifts heavily favors developmental changes to accessibility (see Supplementary Figures 4, 5, Supplementary Information).

As accessible chromatin regions often serve to regulate gene expression [39] by acting as cis-regulatory sequences, we next sought to identify the potential role of our regions in regulatory changes occurring during development and aging. We did this by considering promoter-level accessibility (see Methods, Supplementary Figure 6), promoter-capture (Hi-C) interactions [41], and regulatory-domain annotations [42] for genes which may be subject to control by these regulatory regions. We found a general pattern for enrichment of immune-related gene sets with the adult-biased set, while development-related (e.g., cellular proliferation) terms were enriched with fetal-biased regions, patterns echoed when considering old-biased and young-biased region sets, respectively (see Supplementary Information, Supplementary Table 2).

We next incorporated tissue expression datasets looking for general gene expression trends between fetal and adult tissues (see Methods). We observe similar enrichment terms as well as significant overlaps with gene sets defined on the regulatory level (see Supplementary Information, Supplementary Table 2). Similarly, we utilized GTEx (gene tissue expression) datasets [43] to look for corresponding shifts in gene expression with age, similar to previous work [44] (see Methods)(Supplementary Table 2). While we did not observe significant overlaps between these aging-expression gene sets and those defined using aging-accessibility changes, we did see significant overlaps with the fetal/adult expression comparisons, along with enriched gene sets with relevance to aging biology (see Supplementary Information, Supplementary Table 2).

As development and aging are phenomena subjected to the actions of random and directed evolutionary forces [13, 45, 46], we next develop expectations for how these epigenetically-altered regions may have evolved over time.

Sequence evolution of epigenetically-altered regions

Development and aging are simultaneously very ancient and variable [45, 47] biological processes and are particularly divergent in key species [48]. Thus, it may be the case that both development and age-associated regions have been shaped by a mix of evolutionary forces acting to either maintain or modify genetic sequences (e.g., regulatory enhancers). To address this possibility, we examined patterns of sequence conservation in our epigenetic datasets using phyloP [49], a measure of nucleotide conservation and/or acceleration across species (Figure 2A). Across primates, we observed that fetal-biased regions tended to have greater sequence conservation than adult-biased regions, and furthermore that both sets differed significantly from those regions not defined as developmentally-altered (Supplementary Table 3). These patterns were similarly observed when comparing age-associated regions (Supplementary Table 3). These findings of conservation differences between sets suggests that the greater regulatory and developmental role associated with fetal-biased and young-biased regions (e.g., enriched for enhancer elements) exerts functional sequence constraint while adult-biased and old-biased regions (e.g., enriched for repressed segments) are less conserved across species.

Within this broader context of species diversity and evolution, humans and chimpanzees display marked and obvious differences in development and longevity [45, 50]. This relatively recent divergence is thought to be driven largely by non-coding changes to cis-regulation [51]. We therefore next looked for evidence of regulatory modifications to biological processes that may contribute to these human/chimp differences. To do this, we intersected our regions sets with sequences demonstrating significant divergence along the human lineage (e.g., ‘human accelerated regions’ [52]). We found that fetal-biased regions were enriched for signals of acceleration while adult-biased regions were depleted (Figure 2B). Similar patterns were seen for young versus old-biased regions (Figure 2B and Supplementary Table 3). We found a number of genes involved in development and aging processes with putative nearby regulatory elements intersecting accelerated regions, two examples of which are shown in Supplementary Figure 7 (see also Supplementary Information).

To gauge the evolutionary interaction between sequence constraint across species and within-species variation, we next assessed modern-day human diversity within region sets (Methods). We found that genetic diversity in fetal-biased regions was markedly reduced (i.e., constrained) compared to genomic backgrounds, as well as to intronic and promoter-TSS elements (Figure 2C and Supplementary Table 3). Conversely, adult-biased regions were enriched for sequence diversity, at the level of annotated repeat elements. These patterns were accentuated when examining young- and old-biased region datasets, and comparing region sets directly (Supplementary Table 3, Supplementary Information). Importantly, when we considered sequence diversity within other ape species, we also observed a decrease in fetal-biased, and young-biased sequence diversity (relative to adult-biased and old-biased, respectively). This latter finding further suggests that fetal-biased regions are associated with conserved regulatory function that discourages mutation and/or drift (Supplementary Figure 8 and Supplementary Table 3).

Overall, natural selection appears to have acted upon regions subject to accessibility shifts in development and aging, modifying some loci (i.e., accelerated divergence indicative of ancient positive selection) while protecting others (i.e., reduced variation indicative of more recent negative selection) (Figure 2D).

Importantly, selective forces, both positive and negative, manifest phenotypically through the effects of random genetic mutations, which act to modify gene regulatory networks to varying degrees. We next examine this relationship.

Epigenetic shifts in age-associated trait associations

In our above analyses, the observed signals of consistent inter- and intra-species conservation in regions most associated with early development (i.e., the fetal-biased set) follows with the expectation that variants negatively impacting early-life would be subject to stronger purifying selection [9, 53, 54]. Conversely, variants with later-manifesting effects, i.e., those within regions increasing in local accessibility with age (i.e., the adult-biased set), would be subjected to substantially weaker selection and may therefore be ‘tolerated’ [55, 56]. To test expectations of the possible deleterious effects of variants subject to accessibility change over development and aging, we utilized GWAS datasets available from the UK Biobank [57]. We extracted summary-statistics for a collection of 127 complex diseases/pathologies falling into aging-related categories [58], including metabolic disorders, cancers, cardiovascular disease, and musculoskeletal impairment (Supplementary Table 4, Methods). We similarly analyzed a set of developmental trait GWAS to act as a control for our fetal/adult accessibility comparisons, and finally considered longevity GWAS data (Supplementary Table 4, see Supplementary Information).

It has been suggested that the highly polygenic nature of complex traits and diseases reflects cumulative regulatory modification to a ‘core’ set of genes which functions most proximately in relevant biology [59]. Across age-associated diseases, this may reflect general aging processes, and regulatory variants impacting these would be expected to contribute to heritable aging-disease risk broadly. Given this rationale, we first considered the behavior of individual SNPs nearby accessibility-altered regions across diseases, and subsequently these behaviors at the gene-locus level (below). We aggregated per-SNP associations across diseases as a singular cross-set metric of risk association (Supplementary Methods). We confirmed that ClinVar variants, variants for which possible clinical significance have been described [60], tended to be more risk-associated by this metric, as we would expect (Supplementary Table 4). Additionally, across all diseases we individually performed enrichment tests for strongly-associated variants nearby our region sets, which corroborated the cross-disease results described below (Figure 3A, see Supplementary Information).

As variants in accessible non-coding regions likely have regulatory impacts generally [39], we confirmed that variants within or nearby regions not classified as strictly developmental nor age-altered tended to have greater association than non-accessible variants. However, this control set had significantly less association than variants nearby sets of both developmentally- and age-altered regions (i.e., fetal/adult-biased, and young/old-biased regions) (Supplementary Table 4).

Considering first developmental change, we found that variants in regions gaining nearby accessibility in adults (i.e., adult-biased) have greater association with disease than those in regions losing nearby accessibility (i.e., fetal-biased) (Figure 3B). Unexpectedly, when looking at aging accessibility changes, we observed that variants in regions gaining nearby accessibility in older-samples actually have lower cross-disease associations than those in regions becoming more accessible in younger samples (Figure 3B and Supplementary Table 4). Furthermore, we found that for intersections of development and age-altered regions that the increased disease association with adult-biased regions was abrogated when intersected with old-biased regions. The magnitude of region-set differences in disease associations was also greater in the young/old-biased comparisons (see Figure 3B, Supplementary Information).

Taken together, these results would suggest that those variants most accessible in younger adults stand to have the greatest impact (in terms of association p-value) on late-life disease risk – a finding that may have important implications for understanding the development of disease over adulthood (see Figure 3C, Discussion).

We next considered these disease-association patterns at the gene-locus level. Briefly, for a given disease we assign the most significant nearby SNP to all genes, and subsequently rank genes based on their assigned GWAS signal. Gene ranks are then aggregated across diseases, looking for genes consistently ranked higher across sets (Supplementary Methods). To confirm the behavior of this gene-ranking method, we compared the cross-set ranking of genes associated with homeostatic processes (based on GO annotations) to randomized gene sets, finding that these gene loci tend to harbor stronger genetic variants across a larger number of diseases than expected (compared to randomized sets), but not so for genes involved in reproductive organ development (see Supplementary Information, Supplementary Table 4).

We applied this method to the sets of development and age-associated genes we defined above and asked whether they tended to have more or less cross-disease GWAS signals than expected. Our sets defined by accessibility-region contacts supported our earlier findings on strong GWAS signals nearby development and age-altered regions – namely, loci of both adult-biased and young-biased gene sets were enriched for strong GWAS signals across diseases, while fetal-biased and old-biased gene sets were associated with relatively weaker GWAS signals (Supplementary Table 4). Sets defined by RNA-seq data showed more of a mix of enriched/depleted GWAS signals across developmental and age comparisons, reflecting the possibility that a mixture of genes increasing and decreasing expression over time may additively contribute to aging disease biology (see Supplementary Information).

Given our results, found at both genome-wide and gene-locus set levels, we finally sought to take an unbiased approach to identify relevant ‘core’ aging-related genes solely on the basis of aggregate GWAS signals (see Supplementary Methods, Supplementary Information). Overall, we had limited success in defining a set of genes with clear, pan-tissue biological relevance, suggesting that, if such a core does exist, that it may be too broad, or the per-locus signals too moderate, for our method to robustly detect. However, since our results suggest the importance of altered epigenetics in modifying GWAS associations, we performed a similar gene-prioritization analysis using variants occurring nearby altered-accessibility regions (Supplementary Methods). This yielded markedly different enrichments for terms relating to immune processes and gene regulation (see Supplementary Information, Supplementary Table 4). One particular set of genes, involved in histone deacetylation, has repeatedly been linked to aging and epigenetics [61, 62] and was identified using our set of young-age regions (see Figure 4). We explore this set in more detail below.

Sequence evolution and disease association

Our previous analyses found that patterns of inter- and intra-species sequence conservation depended on epigenetic status (i.e., degree of accessibility) of regulatory elements. Subsequently, we found that the risk association of variants across a number of age-associated diseases also varied based on accessibility change in the vicinity of the variant. Much work has been done on understanding the relationship between sequence conservation and disease risk [63–66]. For example, a transcription factor (TF) binding site may be subjected to negative selection to conserve its sequence and hence function. Mutations that occur within this site would more likely impact cis-regulatory biology, and therefore manifest an association with disease. If this disease impacts fitness, then over time, such mutations will be eliminated, so that genetically ‘older’ mutations are less prevalent [65, 67]. Given our interest in the evolution of development and aging processes, we wanted to investigate the role that epigenetics has on this disease-evolution relationship - and whether this holds with our data. By comparing the cross-trait associations of variants falling within and outside primate-conserved sequences (phastCons) [68] (Supplementary Methods), we found that variants within conserved sequences tend to have greater disease associations, along with younger estimated allele age (Supplementary Table 4). These patterns also hold true for phastCons sequences within age- and developmentally-altered regions (Supplementary Table 4), and follow with previously observed enrichments for GWAS associations of conserved, younger (allele age) variants [65, 67, 69].

We next considered primate conservation, estimated allele age, and cross-set association of variants, looking for the effects of nearby accessibility change on these metrics (Supplementary Methods). As an additional metric for predicted fitness consequences, particularly of non-coding variants, we also included per-bp LINSIGHT scores [70], which integrates data on chromatin accessibility, TF binding motifs, and comparative genomics.

First, we found that variants falling near fetal-biased regions were more conserved, younger, and had lower cross-set association, while variants near adult-biased regions behaved oppositely (Figure 3B and Supplementary Table 4). We also found that the predicted functional consequences associated with fetal-biased regions were greater than with adult-biased regions, despite the lower cross-set association with aging-associated diseases (Figure 3C, see Discussion). Variants falling near old-biased regions were less conserved, had older allele age, and had lower cross-set association than their young-biased counterparts (Supplementary Table 4). These old-biased regions were also associated with the lowest predicted functional consequences (in aggregate) of any set, while the set of young-biased regions had the second-highest average. To compare these behaviors with variants of annotated clinical significance we independently examined ClinVar variants, which while demonstrating increased cross-set association, tended to also be more conserved, younger, and have stronger predicted fitness consequences (Supplementary Table 4).

Collectively, our results indicate that variants stratified by nearby accessibility change violate the expected relationship between sequence conservation and disease association (behaviors instead observed for ClinVar variants). Namely, those regions exhibiting the highest sequence constraint (Figure 3C, left) do not also exhibit the strongest aging-disease associations, nor do those regions exhibiting the weakest constraint, as might be expected in a ‘mutation accumulation’ theory of aging [35]. However, when considering predicted functional consequences (LINSIGHT), which are not defined based on aging demographic data, this pattern is reversed (i.e., the most constrained set, fetal-biased regions, had the strongest predicted consequences despite weaker aging-disease associations). This unexpected behavior may have important implications for evolutionary models of late-onset complex disease genetics. Based on our results, we propose such a model suggesting the outsized impact of regulatory sequences active in early adulthood on genetic contributions to aging-associated disease risk (see Figure 3C, Discussion).

Our proposed model suggests that focusing on disease risk loci containing such putative regulatory sequences (i.e., young-biased regions), should implicate sets of genes involved in aging biological processes. Our gene-level GWAS analyses using young-biased regions identified genes involved with histone deacetylation as being more consistently associated with aging-disease GWAS signal, a pattern which was diminished when considering gene-level associations in the absence of this epigenetic information (Figure 4A), and when using other region sets (e.g., old-biased regions) (Figure 4B). Histone deacetylation enzymes have known impacts on epigenetic aging biology [7, 61, 62] and aging diseases [71]. Within our young-biased enriched gene set we identified SIRT6 and SIRT7 as having multiple variants falling nearby young-biased regions which contacted their respective gene promoters (Figure 4C). Both these enzymes have been associated with maintaining heterochromatin during aging [72–74]; SIRT7 decreases expression with age, and antagonizes hMSC epigenetic aging [73], while SIRT6 loss manifests an aging-like state [75]. It may be possible that decreased accessibility of regulatory regions controlling the expression of these genes are involved in decreases in sirtuin expression and heterochromatin [72].

Accessibility datasets

DNase-I hypersensitivity datasets were obtained from ENCODE [40] for eight different fetal and adult tissues (adrenal gland, brain, heart, lung, muscle, skin, spleen and stomach) (see Supplementary Table 1 for accessions and metadata). Raw data was processed as described in the Supplementary Methods, with called open-chromatin regions consolidated across replicates and tissues in order to define a final set of reproducible regions. This aggregated set of peaks was then used to assess both pan-tissue, as well as per-tissue, accessibility changes between fetal and adult datasets using the limma package (version 3.46) in R [81, 82]. Differentially-accessible regions were defined using a Benjamini-Hochberg FDR [83] cutoff of < 0.05.

Adult DNase samples were further stratified in order to define age-altered chromatin regions, splitting samples used in the above analysis into those individuals younger than 50 (‘young-adult’) and those older than 50 (‘old-adult’), this age representing a roughly equal split of sample numbers. Not all tissues used in the initial fetal/adult comparison were represented in these age-stratified sets – thus we restricted the tissue comparisons to brain, heart, lung, muscle and stomach tissues. A similar computational method as that used in defining fetal- and adult-biased regions was applied here (see Supplementary Methods). We compared accessibility changes between young and old-adult samples within those regions exhibiting fetal/adult biases, defining young-biased and old-biased regions (again, using an FDR cutoff of < 0.05).

Promoter accessibility change

All hg19 Refseq gene TSS were obtained from the UCSC Genome Browser [84] and padded 1kb up/downstream to define promoter regions. For each promoter region, DNase read coverage was compared between adult and fetal samples, with resulting data processed using a similar differential-accessibility method as that used above (see Supplementary Methods). Significantly differentially-accessible promoters were defined using an FDR cut-off of 0.05. As an additional, more stringent analysis, we also defined differentially-accessible promoters based on intersections with the above defined region sets (see Supplementary Methods).

Promoter capture datasets: Promoter-capture data was obtained from Jung et al., 2019 [41]; this dataset was generated from promoter-capture assays across a number of different tissues and cell-types. Given our pan-tissue approach, we considered all data (with the exception of OV2, as we excluded sex-specific tissues from all previous obtained datasets). To generate a set of genomic regions which show evidence of contacting gene promoters, we filtered interacting regions to those which contacted their respective promoters in at least two different tissues/cell-types. This moderate filter was used to exclude those regions for which interactions appear to be exclusive to one dataset, while allowing for regions that do not show such exclusivity.

Gene-set enrichment analyses: Gene sets generated in our analyses were tested for enrichment in different GO Biological Process terms using the ‘enrichGO’ function from the clusterProfiler [85] package version 3.16.1, with semantically-similar GO terms collapsed and significantly-enriched terms defined as adjusted p-value < 0.05.

ENCODE RNA-seq datasets: Processed per-gene quantification files, as generated by the ENCODE pipeline were obtained from the ENCODE web portal [40] (see Supplementary Table 2 for file accessions and metadata). Given the limited availability of adult tissue samples for use in differential-expression analysis, we instead defined a less-stringent method to identify broad changes in gene expression which demonstrate consistency across tissues (see Supplementary Methods).

GTEx RNA-seq processing: Processed RNA-seq quantification files were obtained from the GTEx web portal [43] for the following tissues (matching the above young/old-age accessibility comparison): brain (Brain - Cerebellum), heart (Heart – Left Ventricle), lung (Lung), muscle (Muscle - Skeletal) and stomach (Stomach). Similar to the processing performed in Benayoun et al [44], we applied quality filters to remove lowly-expressed and non-coding genes, and subsequently used the same definitions of ‘young-age’ and ‘old-age’ (as in the above analyses) to calculate differential expression using limma-voom (see Supplementary Methods).

Human sequence variation datasets: Variation data from the 1000 Genomes Project phase 3 (1KGP) [86] (n = 2504 individuals) in. vcf.gz format was obtained and intersected with our region sets using tabix [87] (version 1.9) to obtain variants occurring within these altered-accessibility regions. Common variants were defined using a minor allele frequency (MAF) threshold of >= 0.05. These sets of intersected variants were subsequently used to compare sequence variation across region sets, as well as comparing region-intersected variation to genomic backgrounds and feature sets (see Supplementary Methods).

GWAS summary statistics data: To define a set of aging-associated diseases for use in our analyses, we first used broadly-defined categories as described in Chang et al., 2019 [58]. This study described 92 age-related diseases grouped into broader disease categories based on analyses of large-scale demographic datasets. We took these classifications and manually extracted relevant GWAS phenotypes assessed by the UK Biobanks study [57], obtaining pre-processed summary statistics for these phenotypes provided by the Neale lab [77] (https://nealelab.github.io/UKBB_ldsc/downloads.html). These data were subsequently utilized across several bioinformatic analyses (see Supplementary Methods).

Additional computational methods, including implementations of statistical tests described in the Results, are described in detail in the Supplementary Methods document. The datasets supporting the findings of this study are publicly-available – accession codes and URLs are provided in the Supplementary Methods and Tables. Computational code for processing these datasets is available upon reasonable request from the Lead Contact.