Dysfunctional epigenetic aging of the normal colon and colorectal cancer risk

This study included 334 tissue samples of normal colon mucosa collected at the University of Washington Medical Center (Seattle, WA, USA) by endoscopic biopsy from patients undergoing colonoscopies (age 19–85) [27] and by surgical resection from newly diagnosed CRC patients (age 28–89, stages I–IV) [28, 29], following the protocols approved by the Institutional Review Board. To avoid the potentially confounding effects of anatomic location, only the samples from the left colon were included in the study. Genome-wide DNA methylation levels were assessed using the Illumina Infinium HumanMethylation450 (HM450, N = 120, completed years 2012–2016) and Infinium MethylationEPIC (EPIC, N = 214, completed years 2017–2019) BeadChip arrays.

Risk group assignment was based on the subject’s personal history of adenomas or CRC, which is known to associate with the risk of developing CRC in the future [2]. We defined three risk groups: low, which was based on no concurrent adenomas; medium, which was based on non-advanced adenomas or advanced adenomas (defined as being an adenoma > 1 cm or having tubulovillous histology or high-grade dysplasia); and high, which was based on concurrent CRC. Table 1 summarizes the risk groups and characteristics of the study subjects. We adjusted for the clinical covariates, especially gender and age, and corrected batch effects in our analyses (Additional file 1: Figure S1).

DNA extraction and bisulfite conversion were performed as described previously [29]. In brief, genomic DNA samples were extracted from the fresh frozen normal colon mucosa tissue samples using the DNeasy Blood and Tissue Kit (Qiagen). Genomic DNA quantification was performed using the Quant-iT PicoGreen DNA assay kit (Life Technologies). DNA (500 ng) from each sample was bisulfite converted using the EZ DNA Methylation Kit (Zymo Research, Irvine, USA). The DNA samples were submitted to the Genomics Core at the Fred Hutchinson Cancer Research Center where they were processed and run on HM450 or EPIC arrays following the manufacturer’s instructions (Illumina, Inc.). The returned raw intensity (IDAT) files were then preprocessed and normalized as described below.

The raw IDAT files of the two methylation arrays were read into R with the minfi package separately [30]; the combineArrays function was utilized to combine the two arrays’ data together based on their common CpG sites. Then, the data was preprocessed with background and dye bias correction using the Noob method [31], followed by the functional normalization [32]. CpG probes that were SNP-associated, cross-reactive, located on sex chromosomes, and unreliably detected (> 10% of samples with detection p value > 0.01), with the exception of the epigenetic clock CpGs, were excluded from the analysis [33–35]. Methylation β value for each CpG site in each sample was calculated as M/(M + U + α), where M and U represent methylated and unmethylated signal intensities at the CpG site, respectively, and α is an arbitrary offset (usually 100) intended to stabilize β values where fluorescent intensities are low. NA values, if existed after the QC filtering, were imputed as the means of all non-NA values of the corresponding CpGs. The β values were transformed into M values as log₂(β/(1 – β)), and the batch effects were removed based on the M values using the Combat approach [36].

Of note, 64% of the samples were run on the EPIC array platform, while the rest were run on the HM450 array (Table 1 and Additional file 1: Figure S2). In our study, we first analyzed the data of the two array platforms separately and found comparable results with regard to the determined epigenetic age and acceleration and CRC risk; we also performed an inverse variance-based meta-analysis [37] to combine the testing statistics and p values of the two datasets and confirmed the results of using the combined dataset are similar to those obtained from the separate datasets (see results in the “Discussion” section). Therefore, we finally combined the EPIC and HM450 datasets to increase the sample size and gain more statistical power for the studies we conducted.

Cell-type heterogeneity may cause somatic DNA methylation variation between tissue samples and may be an important confounder in the study of DNAm and epigenetic age alterations in association with CRC risk in the normal colon tissues [18]. Therefore, we used EpiDISH [38], a reference-based algorithm for the inference of cell-type proportions in cell mixture samples, to estimate the fractions of epithelial cells, fibroblasts, and total immune cells in our samples.

The epigenetic ages of each sample were estimated using 4 popular epigenetic clocks, which were the Hannum clock, which relies on 71 CpGs identified in blood DNA samples [21]; the Horvath clock, which relies on 353 CpGs and is based on the analysis of DNA methylation from multiple tissue types [22]; the PhenoAge clock, which is based on 513 CpGs derived to measure phenotypic aging [23]; and the EpiTOC clock, which is derived from the analysis of 385 Polycomb group target promoter CpGs [24]. Note these epigenetic clocks were developed using data from the HM450 array. Although the EPIC array lacks some of the CpGs in the Hannum (6), Horvath (19), and EpiTOC (31) clocks due to the differences in the array design between the HM450 and EPIC arrays (Additional file 1: Figure S3), McEwen et al. have demonstrated that the missing clock CpGs on the EPIC array do not substantially affect the accuracy of the Hannum or Horvath age determination [39]. To verify this observation for all 4 clocks, we performed a sensitivity analysis on our HM450 data. We selected the common clock CpGs on both arrays to calculate the epigenetic ages for the HM450 samples and compared these results with their epigenetic ages derived from using all clock CpGs (see results in the “Discussion” section).

A linear regression model was used to describe the relationship between the epigenetic age and chronological age at the time of tissue collection. The deviation between epigenetic age and chronological age, also known as epigenetic age acceleration, was calculated for every sample based on the residuals of regressing the epigenetic ages on the chronological ages of all the samples, as described by McEwen et al. [39].

The correlation between epigenetic ages and chronological ages of the samples was calculated with the Pearson correlation coefficient. To investigate the change of this association between the different CRC risk groups, a linear regression model with interaction effect was adopted by considering the chronological age as a linear predictor and the risk status as a categorical predictor. To test the association of epigenetic clock with cancer risk, a multivariate linear regression was applied with the epigenetic age acceleration as the dependent variable and the cancer risk as the independent variable with adjustment for other covariates, such as gender and cell-type fractions. When considering other clinical variables that might be able to affect DNA methylation or aging, such as the BMI, smoking, and NSAID use (Table 1), a subset of 293 samples that had no missing data of these variables were analyzed. The Fligner-Killeen test was used to test the homogeneity of variances between different CRC risk groups.

We implemented the Surrogate Variable Analysis (sva) [40] on the DNAm data (M values) by setting a null model matrix (mod0 ~ age + gender) and a full model matrix (mod ~ risk + age + gender) to estimate the surrogate variables (SVs) that represented other latent confounding factors. Then, an epigenome-wide association study (EWAS) was performed to identify DNAm changes between different CRC risk groups using a multivariate linear regression model, where DNAm level of a CpG was the outcome, CRC risk status was the independent variable of interest, and age, gender, and SVs were the adjustment variables. The output included effect size (i.e., M value mean difference) and p value for each CpG. The false discovery rate (FDR)-adjusted p values were calculated for the multiple testing adjustment. FDR < 0.01 was used to determine the differentially methylated CpGs between different risk groups. One-way Fisher’s exact test was performed on a two-by-two contingency table, which contained the numbers of differentially methylated CpGs (from the EWAS analysis) in a clock and in a whole array as well as the numbers of total CpGs in the corresponding clock and array to test if the clock was enriched with more differential CpGs. Gene Ontology (GO) functional annotation for the genes close to the differentially methylated clock CpGs was analyzed using the online Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 [41].

To study epigenetic aging in the normal colon and its correlation with the risk for developing colorectal cancer, we conducted a genome-wide DNA methylation study on normal colon mucosa samples collected from patients assigned to the low, medium, and high CRC risk groups (N = 105, 128, and 101, respectively, see Table 1 for subject information). We calculated the Pearson correlation coefficient between the epigenetic age and chronological age for subjects in each risk group and for the combined set of samples. We found that the Horvath clock had the strongest and most significant correlation with the chronological age in all of the groups (r > 0.8, p < 1 × 10⁻³⁰), while the Hannum, PhenoAge, and EpiTOC clocks showed weaker correlations with chronological age, particularly in the high-risk group (Fig. 1). The linear regression with interaction effect did not reveal significant changes in the association of epigenetic age with chronological age between the risk groups (p > 0.05 in all tests). Our observations demonstrate that the Horvath clock is the most accurate clock for predicting the chronological age in the normal colon and that the epigenetic ages derived from the Hannum, PhenoAge, and EpiTOC clocks diverge from the chronological age of the samples.

Next, we analyzed the association between the epigenetic ages of the samples and their CRC risk status to determine whether the biological tissue age of the colon was associated with an increased risk of developing CRC. We initially assessed the epigenetic ages of the samples using the four different epigenetic clocks and found an older mean age in the medium-risk group for the Horvath clock but a younger mean age in the high-risk group for the Hannum, PhenoAge, and EpiTOC clocks compared to the low-risk group (Additional file 1: Figure S4). In addition, the analysis using the Fligner-Killeen test revealed that the Hannum (p = 7.5 × 10⁻⁴) and PhenoAge (p = 0.047) epigenetic ages of the high-risk samples had a significantly larger amount of variances compared to the low-risk samples. To adjust for the bias due to individual chronological age, we assessed epigenetic aging using a popular measure named epigenetic age acceleration, which was obtained from the residuals of regressing epigenetic ages of all the samples onto their chronological ages [39].

Cell-type heterogeneity may cause somatic DNA methylation variation between different groups of colon samples and may be an important confounder affecting epigenetic clocks in association with CRC risk [18]. Therefore, we estimated the fractions of epithelial cells, fibroblasts, and total immune cells within each sample using the EpiDISH algorithm [38]. We found that cell-type fractions were highly correlated with top PCs of the DNAm data (Additional file 1: Figure S5A), indicating cell-type heterogeneity had a considerable influence on the DNAm results of the samples. The percentage of fibroblasts was significantly higher (p = 2.6 × 10⁻¹¹), and the percentage of total immune cells was significantly lower (p = 2.8 × 10⁻⁹) in the high-risk group compared to the low-risk group (Additional file 1: Figure S5B). We also found that the estimated cell-type fractions were correlated with the Hannum, PhenoAge, and EpiTOC age estimates but were not correlated with the Horvath age estimate (Additional file 1: Figure S5C), perhaps because the Horvath clock was built on large-scale multi-tissue data and hence can adjust for the influence of cell-type heterogeneity intrinsically, while the other three clocks are more sensitive to changes in cell-type composition.

By taking the potential confounders into consideration, we used a multivariate linear regression with the adjustment for the gender and the estimated cell-type fractions to test the difference of epigenetic age acceleration between the CRC risk groups. We observed significant deceleration of PhenoAge (p = 0.0012) in the high-risk samples compared to the low-risk samples (Fig. 2). A similar phenomenon was also observed after additionally adjusting for other relevant clinical variables, such as BMI, smoking, and NSAID use, in the regression model using a subset of 293 samples that had sufficient clinical annotation for these variables (Additional file 1: Figure S6).

To validate our results and considering the lack of datasets from genome-wide DNAm analysis of normal colons from healthy individuals, we generated 2 HM450 datasets by (1) combing the raw IDAT files of our low-risk normal colon samples (N = 48, UWAS-Low) with the TCGA-COAD adjacent normal left colon samples (N = 9, TCGA-High) and (2) combing the raw IDAT files of our 48 low-risk normal colon samples with normal left colon samples from patients with CRC from the Australian Melbourne Collaborative Cohort Study (N = 14, MCCS-High) [42]. We repeated all the analyses using these datasets and observed epigenetic age deceleration in the TCGA-High group for the Hannum, PhenoAge, and EpiTOC clocks and in the MCCS-High group for the Hannum and EpiTOC clocks (Additional file 1: Figure S7, p < 0.05).

We wish to note that it might be also interesting to estimate the epigenetic age of matched cancer tissues from the high-risk group patients; hence, we assessed the CRC samples from a subset of the people with cancer (N = 13). We combined all the normal colon and CRC samples together and assigned them into four groups (low, medium, high, and CRC). We determined the epigenetic age in these samples and found that the Horvath clock was significantly decelerated while the PhenoAge and EpiTOC clocks were significantly accelerated in the CRC samples (Additional file 1: Figure S8).

In light of our observation of deceleration of the epigenetic clocks in the high-risk normal colon samples, we next assessed for epigenome-wide methylation changes in association with CRC risk status of the samples in order to determine if the risk-associated methylation changes in the normal colon were skewing the performance of the epigenetic clocks. We performed an EWAS analysis on the methylation data of all 334 samples to identify genome-wide DNA methylation changes that were associated with cancer risk by applying a multivariate linear regression to each CpG. Using a significance threshold of FDR < 0.01, we identified 14,947 differentially methylated CpGs in the high-risk group compared to the low-risk group (see Manhattan, QQ, and histogram plots of EWAS p values in Additional file 1: Figure S9A). We noticed that 5 of the Hannum clock CpGs, 18 of the Horvath clock CpGs, 20 of the PhenoAge clock CpGs, and 20 of the EpiTOC clock CpGs were differentially methylated in the high-risk group vs. the low-risk group (see volcano plots in Additional file 1: Figure S9B). None of the four clocks was significantly enriched with differentially methylated CpGs (Fisher’s exact test p value > 0.01). Gene Ontology (GO) functional annotation for the genes close to these differentially methylated clock CpGs indicated that they were significantly relevant to the biological process of “cardiac cell fate determination” (p value < 0.01). We further investigated the relationship between the methylation mean differences of the clock CpGs in the high-risk group derived from the EWAS and their coefficients (or weights) in the weighted sum-based epigenetic clock models (i.e., y = βX). We wish to note the EpiTOC clock is an average methylation model, where each clock CpG has the same coefficient that is 1/385. We multiplied the coefficient of each clock CpG by its EWAS mean difference to quantify its overall mean difference in terms of the epigenetic clock (see the scatter plots in Additional file 1: Figure S9C). Although some CpGs in the PhenoAge clock were significantly hypermethylated in the high-risk group, their negative clock coefficients made them contribute to the observed age deceleration. In contrast, the EpiTOC clock was directly affected by methylation changes of the clock CpGs.

The raw methylation IDAT files and processed data are available online at the GEO website (Accession number is GSE132804↗).