Analysis of the machinery and intermediates of the 5hmC-mediated DNA demethylation pathway in aging on samples from the MARK-AGE Study

Study population

The sample of study consisted of PBMC from 188 volunteers enrolled in eight European countries (i.e. Austria, Belgium, Finland, Germany, Greece, Italy, The Netherlands, Poland) covering the age range between 34 and 74 years. The population was stratified into three age groups (i.e. 34-48; 49-65; 66-74) sharing comparable sample size and gender composition. However, female representation was slightly larger than males in the 34-48 age class. In all age classes the mean value of the body mass index (BMI) was largely below 30 kg/m², which is the threshold for obesity. Moreover, in agreement with literature data, the BMI appeared to increase with age indicating that the population analyzed was representative of a physiological aging process. Another feature given in Table 1 includes self-reported smoking status of the population studied, which was almost similar between age groups. About the half of the subjects had never smoked while the remaining part consisted of current or former smokers (Table 1).

Characterization of data values

The parameters analyzed were tested for normal distribution and most of them failed to pass the Kolmogorov-Smirnov and Shapiro-Wilk normality tests (data not shown). Due to skewness and high kurtosis of not-transformed values, data were log-transformed but, also in this case, parameters did not show normal distribution. As such, both parametric and non-parametric tests were applied to verify the robustness of correlations and the results of group comparisons. The generalized linear model (GLM) analysis was performed to investigate the influence of covariates and multiple confounding factors, including recruitment center, gender and the lymphocyte to monocyte (lympho/mono) ratio. The latter was taken into account as PBMC are a mixed cell population, mainly characterized by lymphocytes and monocytes, and the relative composition of each cell type has been demonstrated to affect methylation studies [10,23]. Furthermore, to rule out the possibility of misleading interpretation of data due to potential batch effects, the Partek Genomic Suite 6.6 ANOVA tool was run on log-transformed values adopting two diverse correction procedures: (1) by using default settings; (2) by attempting to retain age and gender differences [10].

Age-dependent expression of TET and TDG genes in PBMC

The assessment of TET1, TET2, TET3 and TDG genes mRNA expression was performed by RT-qPCR. Spearman’s correlation analysis yielded a highly significant negative linear association between TET1 expression and age for both not-transformed and log-transformed values, which was even more pronounced after batch effect correction. Pearson’s correlation analysis also supported a significant linear decline of TET1 transcript with age (Fig. 1, upper panels). No association with age was observed for TET2 expression (Fig. 2), whose values showed high variance. Accordingly, cluster analysis performed on TET2 data identified two subgroups separated around the 75th percentile value (Supplementary Figure 1A and B), where the one with higher TET2 levels showed a median value six times greater than the other subgroup (data not shown). These two subgroups of TET2 expression data were correlated separately to age but, even in this case, no association was obtained (Supplementary Figure 1C and D). Regarding TET3, both parametric and non-parametric analyses revealed a negative correlation with age, which was even more evident after the removal of batch effects (Fig. 3, upper panels). The relationship of TET1 and TET3 genes with age was further tested against gender and center by bootstrapped regression analysis (Supplementary Table 1 and 2). Notably, the negative correlation of TET1 and TET3 with age was retained in all conditions, indicating that even if an influence of gender and geographical origin could exist, age affects their expression almost independently. The same analyses confirmed the lack of association with age of both TET2 and TET2 expression clusters independently of gender and recruitment center (data not shown).

Figure 1. Age-related changes of TET1 mRNA levels in PBMC. Upper panels show scatter plots representing the linear correlation between TET1 mRNA levels and age in PBMC calculated from (A1) non-transformed TET1 data, (B1) log-transformed TET1 data, (C1) batch-corrected TET1 data, (D1) batch-corrected TET1 data retaining age and gender differences. Parametric (Pearson r) and non-parametric (Spearman’s ρ) correlation coefficients and statistical significance are given above each graph. Lower panels show bar graphs reporting the expression levels of TET1 gene in three different age classes calculated from (A2) non-transformed TET1 data, (B2) log-transformed TET1 data, (C2) batch-corrected TET1 data, (D2) batch-corrected TET1 data retaining age and gender differences. Boxplots show the median, the interquartile range (boxes) and the 5–95% data range (whisker caps). Comparisons between groups were performed by the Kruskal-Wallis test followed by post-hoc Bonferroni test (*P < 0.05; **P < 0.01). (y)= years.

Figure 2. Age-related changes of TET2 mRNA levels in PBMC. Upper panels show scatter plots representing the linear correlation between TET2 mRNA levels and age in PBMC calculated from (A1) non-transformed TET2 data, (B1) log-transformed TET2 data, (C1) batch-corrected TET2 data, (D1) batch-corrected TET2 data retaining age and gender differences. Parametric (Pearson r) and non-parametric (Spearman’s ρ) correlation coefficients and statistical significance are given above each graph. Lower panels show bar graphs reporting the expression levels of TET2 gene in three different age classes calculated from (A2) non-transformed TET2 data, (B2) log-transformed TET2 data, (C2) batch-corrected TET2 data, (D2) batch-corrected TET2 data retaining age and gender differences. Boxplots show the median, the interquartile range (boxes) and the 5–95% data range (whisker caps). Comparisons between groups were performed by the Kruskal-Wallis test followed by post-hoc Bonferroni test. (y)= years.

Figure 3. Age-related changes of TET3 mRNA levels in PBMC. Upper panels show scatter plots representing the linear correlation between TET3 mRNA levels and age in PBMC calculated from (A1) non-transformed TET1 data, (B1) log-transformed TET3 data, (C1) batch-corrected TET3 data, (D1) batch-corrected TET3 data retaining age and gender differences. Parametric (Pearson r) and non-parametric (Spearman’s ρ) correlation coefficients and statistical significance are given above each graph. Lower panels show bar graphs reporting the expression levels of TET1 gene in three different age classes calculated from (A2) non-transformed TET3 data, (B2) log-transformed TET3 data, (C2) batch-corrected TET3 data, (D2) batch-corrected TET3 data retaining age and gender differences. Boxplots show the median, the interquartile range (boxes) and the 5–95% data range (whisker caps). Comparisons between groups were performed by the Kruskal-Wallis test followed by post-hoc Bonferroni test (*P < 0.05; **P < 0.01). (y)= years.

Additional evidence of a link between aging and TETs gene expression was obtained by stratifying samples into three age classes. The age group including the younger individuals (34-48y) showed significantly higher expression of both TET1 and TET3 compared to the 66-74y and the 49-65y age groups, respectively. This was even more evident after correction of batch effects when preserving gender and age differences (Fig. 1 and 3, lower panels). Thereafter, more detailed analyses were adopted to determine the influence of covariates and multiple confounding factors including recruitment center, gender and lympho/mono ratio on the age group differences. In particular, GLM analysis revealed that the TET1 expression difference between age groups resisted the confounding factors in the case of batch effect corrections (Supplementary Table 3). Due to the notable impact of leukocyte composition of PBMC on TET1 expression, bootstrapping-coupled regression analysis was performed for TET1 and age stratifying for the lympho/mono ratio. This analysis definitively demonstrated that age independently affected TET1 expression (Supplementary Table 4). Similar approaches confirmed that the expression of TET2 was not associated with age but rather affected by lympho/mono ratio (Fig. 2, lower panels, and Supplementary Table 5). Regarding TET3 expression, it was sensitive to all tested variables, especially to gender, but none of them affected the observed differences between age groups (Supplementary Table 6).

The analysis of TDG expression revealed a slight negative correlation with age (Fig. 4, upper panels), which was also maintained after batch effect corrections. Regression analysis with bootstrapping supported the decline of TDG with increasing age (Supplementary Table 7), which probably was too moderate for being endorsed by the analysis by age classes (Fig. 4, lower panels). GLM analysis uncovered an influence of the recruitment center on TDG expression (Supplementary Table 8).

Figure 4. Age-related changes of TDG mRNA levels in PBMC. Upper panels show scatter plots representing the linear correlation between TDG mRNA levels and age in PBMC calculated from (A1) non-transformed TDG data, (B1) log-transformed TDG data, (C1) batch-corrected TDG data, (D1) batch-corrected TDG data retaining age and gender differences. Parametric (Pearson r) and non-parametric (Spearman’s ρ) correlation coefficients and statistical significance are given above each graph. Lower panels show bar graphs reporting the expression levels of TDG gene in three different age classes calculated from (A2) non-transformed TDG data, (B2) log-transformed TDG data, (C2) batch-corrected TDG data, (D2) batch-corrected TDG data retaining age and gender differences. Boxplots show the median, the interquartile range (boxes) and the 5–95% data range (whisker caps). Comparisons between groups were performed by the Kruskal-Wallis test followed by post-hoc Bonferroni test. (y)= years.

DNA methylation analysis of TET1 and TDG genes

Focal DNA hypermethylation events are part of age-associated changes of DNA methylation patterns and, if occurring on transcriptional regulatory regions, they may lead to down-regulation of gene expression [8]. DNA methylation-mediated control of the genes showing age-related decline has only been demonstrated for TET1 and TDG. In particular, hypermethylation of TET1 and TDG CGIs but also of the TET1 CGI 3'-shore has been associated with transcriptional repression [24–27]. Based on this evidence, MassARRAY EpiTYPER platform was used to determine the DNA methylation profile of TET1 and TDG genes in young (34-41y) and old (69-74y) individuals. The analysis of the TET1 CGI revealed a slight but significant hypermethylation of some CpGs in the elderly group (Fig. 5). On the contrary, no other CpGs, neither in the TET1 CGI 3'-shore nor in the TDG CGI, showed any statistically significant change of the DNA methylation levels (Supplementary Fig. 2).

Figure 5. DNA methylation profile of TET1 CGI in aging. Graphs represent the CpGs of TET1 CGI analyzed by the epiTYPER assay and show the difference in DNA methylation between the groups of young (34-41) and old (69-74) individuals. Statistical significance was obtained by the Mann-Whitney test (*P < 0.05; **P < 0.01; ***P < 0.001). n(34-41y)=36; n(69-74y)=34. (y)= years.

Content of 5hmC and its derivatives in PBMC in aging

Analysis of global content of 5hmC, performed by dot blot assay with specific antibodies to 5hmC, showed that its levels decreased linearly with age as demonstrated by Pearson’s and Spearman’s correlations for both non-transformed and log-transformed data. This result was even more evident after correction of batch effects (Fig. 6, upper panels). Regression analysis performed by using bootstrap resampling stratified for gender and recruitment center supported the solid negative association between age and 5hmC content (Supplementary Table 9). Confirmatory results were also obtained when samples were stratified in age classes. In fact, the group of younger individuals (34-48y) was significantly different from the remaining two groups (49-65y and 66-74y) (Fig. 6, lower panels). GLM analysis confirmed the negative association of 5hmC with age and also revealed a potential impact of recruitment center with no consequences on the age-dependent reduction of 5hmC though (Supplementary Table 10).

Figure 6. Age-related changes of 5hmC levels in PBMC. Upper panels show scatter plots representing the linear correlation between 5hmc levels and age in PBMC calculated from (A1) non-transformed 5hmC data, (B1) log-transformed 5hmC data, (C1) batch-corrected 5hmC data, (D1) batch-corrected 5hmC data retaining age and gender differences. Parametric (Pearson r) and non-parametric (Spearman’s ρ) correlation coefficients and statistical significance are given above each graph. Lower panels show bar graphs reporting the levels of 5hmC in three different age classes calculated from (A2) non-transformed 5hmC data, (B2) log-transformed 5hmC data, (C2) batch-corrected 5hmC data, (D2) batch-corrected 5hmC data retaining age and gender differences. Boxplots show the median, the interquartile range (boxes) and the 5–95% data range (whisker caps). Comparisons between groups were performed by the Kruskal-Wallis test followed by the post-hoc Bonferroni test (*P < 0.05; **P < 0.01). (E) Representative dot blot performed on DNA from 20 individuals by using anti-5hmC antibody and methylene blue (MB) staining to control DNA loading. (y)= years.

As TET enzymes are also responsible of the conversion of 5hmC into 5fC and then 5caC, which are considered DNA demethylation intermediates, their levels were also determined by dot-blot analysis. Due to the low genomic content of 5fC and 5caC and the limited amount of DNA available, DNA from multiple samples was pooled into 9 groups of increasing age. This analysis confirmed the decreasing trend of 5hmC with aging and revealed an accumulation of 5caC in older ages. By contrast, the levels of 5fC appeared to be comparable between groups (Fig. 7). Notably, correlation analyses showed a significant negative association between 5hmC and 5caC levels (Pearson r: -0.771; p: 0.015; Spearman ρ: -0.767; p: 0.021).

Figure 7. Age-related changes of 5hmC, 5fC and 5caC levels in PBMC. (A) The graph shows the amount of 5hmC, 5fC and 5caC determined by dot-blot assay on pooled DNA samples obtained by grouping individuals into nine different age classes. n(34-36y)=16; n(37-39y)=17 ; n(40-44y)=16; n(45-49y)=22; n(50-53y)=16; n(54-59y)=16; n(60-66y)=31; n(67-70y)=27; n(71-74y)=27. Analysis was performed with 5hmC, 5fC and 5caC specific antibodies. Methylene blue (MB) staining was used to monitor DNA loading. (B) Bar graphs show the densitometric quantification of 5hmC, 5fC and 5caC signal after normalization for loading by MB staining, shown as mean ± S.E.M. of three different technical replicates. (C) Linear dsDNA containing unmodified (C), full methylated (5mC) or hydroxymethylated (5hmC) cytosines were used as specificity control of the anti-5hmC antibody. DNA derived from HEK293T overexpressing TET1 catalytic domain (OE TET1) was used as positive control for anti-5fC and anti-5caC antibodies. (y)= years.

Correlation between components of DNA demethylation machinery, 5hmC levels and other epigenetic enzymes

A putative correlation between the expression of TET1-3 and TDG and the levels of 5hmC was investigated. A slight significant positive correlation was only detected by non-parametric test between TET1 expression and 5hmC levels after batch effect correction (Table 2). Further correlation analyses were then performed to trace back a possible co-regulative network involving TET genes, TDG and other epigenetic factors functionally related to DNA methylation/demethylation pathways. In particular, data from the same subset of PBMC samples were available for DNA methyltransferases DNMT1 and DNMT3B and for poly(ADP-ribose) polymerases PARP1 and PARP2. The most evident associations were obtained for TET1, which exhibited significant positive correlations with TDG, DNMT1, DNMT3B, along with PARP1 and PARP2. A strong correlation was also found between TDG and PARP2. Notably, correlations between these variables were significant after both parametric and non-parametric tests, with or without batch effect corrections. On the contrary, other associations, such as the one between TET3 and TET2 or PARP1, or between this latter and TDG, were not confirmed by either test or they only emerged after the removal of batch effects (Table 2 and 3).

Analysis of clinical variables associated with the bimodal distribution of TET2 mRNA levels

Apart from data on TET2 gene expression, all other parameters investigated displayed a homogenous distribution of values and a similar variance. In fact, TET2 values split the samples into two main subgroups, whose classification in high and low TET2 expression classes was not ascribable to demographic characteristics including gender, age or the country of origin (Supplementary Figure 1 and Supplementary Table 5). The availability of a large set of clinical parameters for the individuals enrolled prompted the investigation of the possible sources underlying the high variability of TET2 expression. For this purpose, a decision tree analysis was performed, which identified the serum levels of alanine aminotransferase (ALT) as the most important factor influencing TET2. Minor subgroups were determined on the basis of serum glucose and urine creatinine concentrations. In particular, individuals showing high TET2 expression also had elevated ALT and intermediate glycaemia while about one third of them showed lower levels of urine creatinine concentrations (Fig. 8). However, as the subdivision of TET2 values based on serum glucose levels fell within the normoglycemic range and the low creatinine subgroup concerned very few samples, only the relationship of TET2 with ALT levels was examined further. In particular, GLM analysis, with gender, age, lympho/mono ratio and recruitment center as covariates, definitively confirmed the positive association between ALT levels and TET2 expression (Table 4).

Figure 8. Identification of major variables that affect TET2 expression by decision tree analysis. Decision tree analysis was performed to identify potential variables responsible of TET2 gene bimodal distribution. Apart from demographic characteristics, several clinical biochemistry parameters were included.

Study population, recruitment, data and blood collection

PBMC samples used in the present work derived from donors aged 34-74 years that have been recruited in the context of the MARK-AGE project [21,22]. The PBMC isolation procedure, details of the recruitment procedure and of the collection of anthropometric, clinical and demographic data have been published [68,69].

Briefly, PBMCs deriving from donors of eight European countries were isolated from EDTA-whole blood, obtained by phlebotomy after an overnight fasting, by discontinuous density gradient centrifugation in Percoll and subsequently cryopreserved, stored in liquid nitrogen and then shipped to the MARK-AGE Biobank located at the University of Hohenheim, Stuttgart, Germany. From the Biobank, coded samples were subsequently sent to the Sapienza University of Rome on dry ice, where they were stored in liquid nitrogen.

RNA extraction and cDNA synthesis

Samples were thawed by incubation at 37°C, followed by dropwise addition of RPMI containing 10% FCS to a final dilution of 1:20. Cells were collected by centrifugation and processed for RNA extraction by using RNeasy Mini Kit (Qiagen) according to the manufacturer's instructions and subjected to DNase I digestion using RNase-free DNase (Qiagen). RNA concentration, purity and integrity were evaluated as previously described [70]. Reverse transcription was carried out using the SuperScript VILO cDNA Synthesis Kit (Invitrogen) on equal amounts of total RNA (0.5 μg).

Real-time quantitative RT-PCR

The expression of TET1, TET2, TET3 and TDG was determined by quantitative PCR using Taqman Gene Expression Assays (Applied Biosystems) following the manufacturer's protocol on the iCycler IQ detection system (Bio-Rad). Each set of primers and probe showed an efficiency of 90–100%. Assays were performed in duplicate with cDNA equivalent to 30 ng of reverse transcribed RNA. Gene expression analysis was performed by the relative calibrator normalized quantification method using the expression level of the β-glucuronidase gene (GUSB) as reference [70]. An inter-run calibration sample was used in all plates to correct for the technical variance between the different runs and to compare results from different plates. The calibrator consisted of cDNA prepared from HEK293T cells. Taqman Gene Expression Assays IDs for each set of primers and probe were as follows: Hs00286756_m1 (TET1); Hs00758658_m1 (TET2); Hs00379125_m1 (TET3); Hs00702322_s1 (TDG); Hs99999908_m1 (GUSB).

Dot blot assay

DNA was extracted from PBMC with DNeasy Blood & Tissue Kit (QIAGEN), samples denatured in 0.4 M NaOH, 10 mM EDTA at 95°C for 10 min, and then neutralized by adding an equal volume of ice-cold 4 M ammonium acetate (pH 7.0). Two-fold dilutions of DNA samples were spotted on a nitrocellulose membrane Hybond-N+ (Amersham Biosciences) in an assembled Bio-Dot apparatus (Bio-Rad). After vacuum-associated blotting procedure, the membrane was washed with 2X SSC buffer, air-dried and blocked with 5% non-fat milk followed by incubation with primary polyclonal anti-5hmC, 5fC, or 5caC antibodies (Active Motif). By the use of HRP-conjugated secondary antibodies, signals were visualized by chemiluminescence (Amersham ECL Western Blotting detection reagents). A total of 250 ng of DNA was used for 5hmC detection, while 750 ng of DNA for 5fC and 5caC. The same amount of DNA from multiple samples was pooled into different groups of increasing age for the analysis of 5fC and 5caC levels. To correct for technical variance between replicates and to compare results from different experiments, a control DNA deriving from HEK293T cells was used in all assays. To monitor DNA loading, filters were stained with 0.02% methylene blue (MB) in 0.3 M sodium acetate (pH 5.2). Densitometric analysis was performed by Quantity One Software (Bio-Rad Laboratories) according to manufacturer’s instructions. To evaluate the specificity of antibodies, 5mC and 5hmC DNA Standard Set (Zymo Research) was used together with DNA derived from HEK293T cells overexpressing the catalytic domain of TET1 (OE TET1).

EpiTYPER assay for quantitative DNA methylation analysis

The EpiTYPER assay (Sequenom) was used to quantitatively assess the DNA methylation state of TET1 and TDG CGIs and TET1 3’-shore. DNA (1 μg) was bisulfite-converted using the EZ-96 DNA Methylation Kit (Zymo Research) with the following modifications: Incubation in CT buffer was performed for 21 cycles of 15 min at 55°C and 30 sec at 95°C and elution of bisulfite-treated DNA was performed in 100 μl of water. PCR was performed on 10 ng of bisulfite-treated DNA using specific primers. The TET1 CGIs and TET1 3’-shore amplicons mapped in chr10:70,320,251-70,320,466 and chr10:70,321,719-70,322,204 (GRCh37/hg19) regions, respectively. TDG amplicon mapped in chr12:104,359,220-104,359,701 (GRCh37/hg19) region.

Bisulfite specific primers were the following:

TDG_Forward: aggaagagagGGTTGGTAGTATTTAGATAGTGGTTGG
TDG_Reverse: cagtaatacgactcactatagggagaaggctAAAACCCAAAATAACACAATAACCTC

TET1 CGI_Forward: aggaagagagGGTTTTTAGTTTTAAGTTTGTATTAGTTTT

TET1 CGI_Reverse: cagtaatacgactcactataGGGAGAAGGCTATCATACAACCCTACCTACCTCTCC

TET1 3’shore_Forward: aggaagagagTTTAAGTTTTTTGATTTTTGGTTTG
TET1 3’shore_Reverse. cagtaatacgactcactatagggagaaggctCTCTTAAAATACCTCTTCCCCTCC

Statistical analysis

Distribution of TET1, TET2, TET3, TDG and 5hmC variables was investigated Kolmogorov-Smirnov and Shapiro-Wilk normality tests. Identification of potential critical variables that can affect TET1, TET2, TET3, TDG mRNAs and 5hmC levels in PBMCs was performed by non-parametric tests (Kruskal-Wallis test or Mann-Whitney U test for two group comparisons) as well as by Generalized Linear Models (GLM). Pairwise comparisons (adjusted for multiple comparisons by Dunn's and Bonferroni’s methods for the Kruskal-Wallis and the GLM tests, respectively) were used to identify significant differences between percentile groups of each categorized variables.

Influence of age and the relative impact of confounding variables on analyzed parameters, was investigated by linear regression and also by stratified bootstrap sampling (1000 bootstrap samples). GLM were used to investigate the influence of confounding variables (tested as categorized and continuous variables) on age-related changes of analyzed parameters. The Bonferroni adjustment for multiple comparisons was used to identify differences between subgroups.

The identification of the major variables affecting TET2 mRNA in PBMCs was performed using decision tree analysis with exhaustive "Chi-squared Automatic Interaction Detector" (CHAID) algorithms. To determine the best split at any node, the CHAID algorithm chooses the predictor variable with the smallest adjusted p-value, i.e., the predictor variable that will yield the most significant split. If the smallest (Bonferroni) adjusted p-value for any predictor is greater than some alpha-to-split value, then no further splits are performed, and the respective node is a terminal node. The process repeats recursively until one of the stopping rules is triggered. In growing the tree the following stopping rules were used: minimum terminal parental node size of 20 cases, minimum terminal child node size of 10 cases and alpha = 0.05 for splitting nodes. The convergence criteria for the Exhaustive CHAID were: epsilon = 0.001 and 100 as the maximum number of iterations before stopping the process.

The potential influence of batch effects was investigated by correcting batch effects with the Partek Genomic Suite 6.6 implemented with the ANOVA tool (Partek Incorporated) using log-transformed variables. We considered each real-time PCR plate (containing 8 samples with study groups not evenly distributed across plates) as a batch. We performed adjustments for batches by applying the following strategies: (i) simple removal of batch effect (including PCR plates as batches in the ANOVA tool); (ii) removal of batch effects attempting to retain differences in age and gender (identified as variable of interest in the ANOVA tool). The same analysis was conducted using log-transformed variables of 5hmC. We considered each dot-blot filter (containing 19 samples with study groups not evenly distributed across plates) as a batch.

All statistical analysis (excluding analysis of batch effects) was carried out using SPSS software (SPSS Inc., Chicago, IL; Version 22.0).