Research into the connection between aging and cancer is being fueled by advances in molecular profiling of age-related processes such as chronic inflammation, accumulation of somatic DNA/mtDNA mutations, and epigenetic changes in tissues in which cancers arise [1]. Two distinct, albeit not entirely independent, concepts have emerged recently that relate changes in DNA methylation to biological tissue age. The first is based on the discovery (of sets) of CpG dinucleotides (CpGs) in the genome that are subject to age-dependent, possibly complex changes in methylation levels that, in aggregate, correlate strongly with chronological age [2–4]. We refer to these types of CpGs as clock-CpGs. A second and simpler concept is based on the observation of gradual age-related changes in methylation levels at specific CpG sites or CpG-rich regions, a process commonly referred to as epigenetic or methylomic drift [5–11]. For example, some CpG islands show very low methylation levels early in life but are known to become gradually methylated over time as a result of sporadic de novo methylation events during DNA replication. We identify these as drift CpGs. It is worth pointing out that data supporting these concepts come mainly from cross-sectional studies that include individuals of different age. In contrast, individual-level (longitudinal) drift, unless studied directly in select individuals over time as we have collected for this study, is typically inferred from population drift.

In a recent study, we used a combination of cross-sectional and longitudinally collected biopsy samples to identify a set of highly correlated CpGs in premalignant Barrett’s esophagus (BE) tissue that undergo differential epigenetic drift relative to normal squamous (NS) tissue [12]. This study arrived at a set of 67 drift CpGs that show significant age-related methylation differences between NS and BE and was used as an epigenetic clock to estimate the time of BE onset. Unlike previous epigenetic clocks that were constructed to predict the age of an individual, the BE tissue-specific clock model was designed to infer unknown tissue ages. Because BE is essentially asymptomatic, it is usually not known how long a patient has lived with BE. However, the time a patient has lived with BE may be considered a risk factor since older BE tissue has had more time for cancer to evolve compared to younger BE tissue which is more likely to be free of neoplastic changes [13, 14]. Although the earlier study was the first to develop an epigenetic clock for BE tissue age, we did not evaluate the full scope of epigenetic drift occurring in BE following its formation, nor its dynamics or functional consequences at the genomic level.

The aims of this study are to characterize individual heterogeneity and genome-wide patterns of age-related epigenetic drift in tissue samples from BE and EAC patients, to develop a mechanistic understanding of methylation dynamics during tissue aging, including the degree of cooperativity and possible presence of nonlinear feedback within CpG islands during the temporal evolution of epigenetic drift, and to explore the impact of advanced drift on gene expression in EACs for which we have both gene expression and methylation data.

To accomplish the aims of this study, it was necessary to analyze methylation drift patterns from an extensive collection of NS, BE, and EAC biopsy tissue samples, including 173 tissue samples from 100 nonprogressed and progressed (nonfamilial) BE patients, along with methylation and gene expression data in 87 EAC from The Cancer Genome Atlas (TCGA) [15] and in 47 EAC + 4 BE samples previously analyzed by Krause et al. [16].

Using data from the HM450 methylation arrays (see “Methods”), we aimed to better characterize the full extent of epigenetic drift occurring in BE and the temporal dynamics of epigenetic drift and to explore the impact of advanced drift on gene expression in EAC for which we have both gene expression and methylation data. Methylation levels are measured either in terms of a β value (methylation fraction) or M value (logit₂(β)) as indicated in the text.

Genomic scope of drift

Out of 146,029 hypomethylated CpG probes in normal squamous (NS) tissue, we identified 18,013 (12%) probes that have significant positive correlation and 560 (0.4%) that have a significant negative correlation (q < 0.01) with the mean differential drift (relative to NS levels) of 67 previously validated drift CpGs in 64 BE samples from patients without a diagnosis of dysplasia or cancer (see “Methods”). In contrast, out of 133,857 CpG probes that were hypermethylated in NS tissue, only 795 (0.6%) probes correlated positively and 3402 (2.5%) probes correlated negatively with the mean methylation drift levels of our 67 probe reference clock (Fig. 1). Thus, significant differential drift in BE involves thousands of CpGs, occurring predominantly in hypomethylated regions that are associated with CpG islands, affecting 4024 (24%) of the 16,984 hypomethylated CpG islands in NS tissue. In contrast, we found only 7% of the identified drift CpG probes to be “open sea,” i.e., isolated in the genome [17], compared to about 10% in the hypomethylated normal background.

The majority (63%) of islands that include drift CpGs are associated with gene promoter regions, i.e., they involve a transcription start site (TSS200 or TSS1500), while only 11% of islands that undergo drift overlap with the gene body, compared to 73% (TSS-associated) and 10% (body-associated) CpG islands on the HM450 array, respectively. In contrast, the relative abundance of intergenic CpG islands is significantly higher among CpG islands that undergo drift compared to the fraction of intergenic CpG islands found on the array (24 vs. 16%, see Fig. 1).

Out of 16,832 island-based (BE-specific) drift CpGs we identified, 1317 unique CpG islands with 5 or more drift CpGs per island (comprising a total of 11,425 drift CpGs). Figure 2 presents a karyograph of 4 autosomes indicating the genomic locations and mean β values for these islands across the 64 BE samples (Additional file 1: Figure S4 for all 22 autosomes). Figure 3 shows a heatmap of the island-level mean β values of the 1317 drift-associated CpG islands for the first 10 NS tissue samples, and all 64 BE samples used to identify CpG probes undergoing drift. For this map, both CpG islands and tissue samples were ordered by their respective mean values. As expected, all 10 normal control tissue samples show no island-level drift. In contrast, we see significant heterogeneity in mean methylation levels of these CpG islands ranging from < 20 to > 80% methylation across the 64 BE samples. With the notable exception of a group of samples that have undergone minimal drift, most BE samples show bimodal patterns of drift where some islands appear to linger at low levels and others show advanced drift. We later use the following categorization for the observed drift patterns in BE and EAC: unimodal low drift (group L), bimodal high drift with a major mode β > 50% (group H) and the remaining bimodal intermediate drift (group I).

Pairwise correlations between island-associated CpGs

CpG islands are considered functional genomic units that may exert transcriptional control by their collective state of methylation rather than through individual CpG sites. To demonstrate this collective behavior in an island-level DNA methylation, we evaluated the pairwise correlations between all island CpGs that are hypomethylated in NS tissue. In general, for static islands (that do not show significant drift), pairwise correlations are moderate (< 0.5) across the span of an island and exhibit anti-correlations near and beyond the island boundaries (Fig. 4). In contrast, island CpGs that drift have stronger pairwise correlations reflecting a collective response of these CpGs to drift consistent with processive DNA methylation maintenance [18, 19].

Figure 4 also shows that the pairwise correlations decay with genomic distance and, for drift CpGs, extend further into the island shelves than static CpGs. Islands that show tissue age-related drift are also significantly larger (in terms of genomic length) than static islands. The mean sizes of the static vs. drift CpG islands were 0.9 vs 1.1 kb, respectively (p = 2 × 10⁻¹⁰, two-sided t test).

Bimodal nature of epigenetic drift in BE and EAC

To see whether methylomic drift is uniformly distributed within our BE and EAC samples, we examined β value distributions for a subset of island-associated drift CpGs with a minimum of 5 detected drift CpGs per island (11,425 drift CpGs in total) for 64 BE and 24 EAC samples from the BETRNet (see “Methods”). Consistent with the patterns seen in Fig. 3, these individual-level distributions show signatures that fall into the arbitrary three types: with unimodal distributions showing low or no drift (group L), distinctly bimodal distributions with intermediate drift (both modes β < 0.5, group I), or bimodal with a major mode near or above β = 0.5 and a minor mode at lower levels (group H). (See Fig. 5a, for an aggregated view of the samples in these groups). While the distributions are similar for BE and EAC, EAC show more advanced drift in the third group (bimodal high) which may be attributed to EAC patients being on average older than the BE patients (68 vs 62 years, respectively), or to the fact that EAC undergoes more frequent stem cell divisions thereby increasing replication-coupled de novo methylation, or to the possibility that BE arises earlier in patients with EAC compared to patients who have not progressed to dysplasia or EAC. We found similar unimodal/bimodal drift signatures in 87 EAC from TCGA and in a combined set of 19 BE and 47 EAC tissue samples provided by Krause et al. [16] (GEO accession number: GSE72874).

Evidence for a threshold effect leading to bimodal drift

The drift patterns shown in Fig. 5a for BE and EAC samples suggest nonlinear drift dynamics in BE tissues. Specifically, the presence of a persistent mode in the drift distribution at low levels (β < 0.2) is indicative of a threshold below which drift is suppressed but advances rapidly once the mean level is surmounted. To validate that epigenetic drift occurs in our longitudinal samples (20 patients with 2 biopsies each separated by at least 3–4 years), we determined for each individual at two time points the number of drift CpGs that remained below (n₁₁), respectively, the number that remained above (n₂₂), and the number of drift CpGs that had crossed the threshold from low to high at β = 0.2 (n₁₂) and, vice versa, from high to low (n₂₁). The results, including the % fraction of drift CpGs advancing, n₁₂/(n₁₂ + n₁₁), and the % fraction retarding between the two time points, n₂₁/(n₂₁ + n₂₂), are listed in Additional file 2: Table S3 and shown as annual rates in Fig. 6. For 19/20 patients, we detect greater methylation flow from sub-threshold levels to higher levels at the second (later) biopsy compared to flow in the opposite direction.

We note that these findings are surprisingly consistent with the unimodal-to-bimodal epigenetic drift predictions made by Sontag et al. [8] who proposed a mathematical model that included a nonlinear relationship between de novo methylation and the ambient level of methylation present in a region of CpGs. To demonstrate that such a model results in unimodal-to-bimodal drift transitions over time, we explicitly simulated sporadic de novo methylation on an island of 50–100 CpGs, independently in 1000 cells, mimicking crudely the cell population in the tissue samples. Each CpG was assumed to be in a binary state (0/1 of being (un)methylated), and the states of the CpGs initially (at time t = 0) were sampled from a binomial distribution with probability 0.06 which equals the mean methylation level in our NS tissue samples. The CpG states were then propagated stochastically with a rate (probability per time step) of becoming methylated that increases 100-fold from a background of 10⁻⁴ to 10⁻² when the mean level of methylation on the island crosses a threshold of β = 0.2.

Without a mathematical exploration of this Markov model, but straightforward in silico experimentation with the baseline distribution of methylation rates (specifically, a gamma distribution with mean 10⁻⁴ and variance 4 × 10⁻⁸) and threshold value, our simulations show that this simple model generates methylation density trajectories that typically bifurcate and strikingly resemble the observed drift signatures in our samples. Figure 5b shows a typical density trajectory for a region of 50 CpGs, an arrayed population of 1000 cells, every 100 time steps, for a total duration of 3000 time steps. Although our model differs in functional form from the model described in Sontag et al., it shares important features, including a suppression of de novo methylation at low levels and a nonlinear acceleration as the ambient (regional) level of methylation increases. In contrast, models that do not include this ambient methylation feedback on the local (site-specific) rate of methylation do not, in general, lead to bifurcations in the main (initial) mode of the evolving drift pattern, but still exhibit a weak bimodality as shown in Additional file 1: Figure S2. Additional file 1: Figure S3 further illustrates the stochastic nonlinear behavior of our model via simulated time course trajectories of the mean methylation levels for 10 CpG islands that share an identical drift rate distribution across their CpGs.

Here we take a closer look at how differential epigenetic drift is organized in BE-associated genomes, and its scope and association with gene expression, motivating further investigation of its role in neoplastic progression in BE. To do so, we first surveyed the array-based DNA methylome for significant correlations with the mean drift measured by 67 drift CpGs previously identified by our group to estimate BE dwell time, i.e., the time a patient has lived with BE [12, 13]. Following this study, we targeted CpGs that are hypomethylated in NS tissue but are subject to differential drift in BE tissue caused by accelerated age-related de novo methylation. While NS tissue may not be the tissue of origin for BE, the similarity of methylation levels at drift-associated CpGs between NS and other normal tissues, such as fundus (see Additional file 1: Figure S1), justifies the use of NS as a normal reference tissue to identify differential drift in BE. Our previous study did not reveal the full extent of this differential drift due to highly restrictive pre-filtering.

Our genome-wide “drift survey” revealed that, at the island level, > 24% of CpG islands undergo methylomic drift and are predominantly promoter-associated (i.e., overlap transcription start sites (TSS)). To investigate whether epigenetic drift occurs on isolated CpG sites or is a nonlocal phenomenon at the CpG island level, we evaluated correlations of methylation between pairs of CpGs (across BE samples) using all island-associated CpG probes (available on the HM450 platform) as a function of genomic distance between the probes (Fig. 4). Our results confirm the prevailing view that CpG islands essentially exert epigenetic control by their collective methylation state rather than through specific CpG sites [25, 26]. Importantly, we found evidence that drift does not evolve uniformly in BE and EAC but appears to be governed by a nonlinear, threshold-like stochastic methylation process which depends nonlocally on the methylation status of other island CpGs. Simulations using a stochastic model, which reflects these dynamics at the island level, show characteristic transitions from unimodal to bimodal drift similar to what we observe in our data. Although other models may provide similar fits to the observed drift distributions, this model has its origins in earlier work aimed at understanding the stable, somatic inheritance of methylation imprints [8] and predicts epigenetic drift as a series of sporadic de novo methylation events at the island level. Our nonlinear feedback model for methylomic drift suggests that the various drift distributions we see in our tissue samples may simply be attributed to tissue age itself (i.e., at what points in time the tissue samples were obtained during the dynamic process of methylomic evolution). Furthermore, analyses of consecutive biopsies in the same patient separated by several years further confirmed that epigenetic drift, as defined in this study, involves the sporadic departure from normal (hypomethylated) levels to higher levels as the tissue ages. Taken together, these findings suggest that epigenetic drift in BE advances non-uniformly by departing from unimodal (low-drift) distributions of methylation and gradually bifurcating into bimodal distributions over time. Similar unimodal and bimodal methylation distributions are observed in EAC samples although the bimodality appears more pronounced in EAC.

To investigate potential functional consequences of epigenetic drift, we compared gene expression in BE and EAC samples showing no (or low) drift to gene expression in samples that show definite drift β > 0.2. This comparison revealed statistically significant differences in gene expression between the two sample groups that are predominantly repressive involving several checkpoint and tumor suppressor-like genes, in particular CHFR (checkpoint with forkhead and ring finger domains), a mitotic stress checkpoint gene that has been observed to undergo promoter-associated hypermethylation in colon, gastric, and esophageal cancers and is associated with chromosomal instability [27, 28]. Submitting the 200 differentially repressed genes in the TCGA EAC samples to a statistical over-representation test (Additional file 2: Table S2) further revealed an unexpected high number of KRAB domain zinc finger genes (greater than sixfold enrichment using DAVID) that are subject to epigenetic drift and transcriptional repression possibly compromising their KAP1(TRIM28)-mediated repressive function. This finding is intriguing because KRAB domain ZNF also target endogenous retroviruses and transposable elements.

Finally, comparison of island-level drift with gene expression in NS and BE tissue samples from the Krause study [16] revealed that the majority of genomic loci undergoing epigenetic drift in BE are transcriptionally silent, consistent with the notion of neutral (clock-like) drift. However, the majority of differentially expressed genes associated with CpG islands that exhibit advanced drift are repressed in EAC when methylation levels increase beyond a threshold of approximately 20%. These findings support the hypothesis that neoplasia, such as dysplastic BE and EAC, may develop in response to epigenetically driven selective pressure exerted on gene expression as methylation levels (on CpG islands associated with gene promoters) advance via random drift beyond a critical, repressive threshold.

Our results are consistent with the hypothesis that epigenetic drift heralds the onset of (epi)genomic instability via bifurcations (as seen in Fig. 5) that associate with the transcriptional repression of important regulatory genes [29–32]. Thus, under this hypothesis, epigenetic drift not only defines tissue aging (i.e., provides a molecular clock) but also “throttles” the expression and function of developmental genes forcing transitions in tissue characteristics that better cope with the erosive and damaging milieu in BE. Further studies of whether changes in methylomic drift simply reflect transcriptional changes during neoplastic progression or induce such changes are therefore of critical importance to better understand mechanisms that drive age-related cancer evolution.