Stochastic epigenetic mutations (DNA methylation) increase exponentially in human aging and correlate with X chromosome inactivation skewing in females

Overall DNA methylation is correlated with age and BMI

Methylation levels of 25.014 CpG sites were evaluated in whole blood from 170 subjects ranging from 3 to 106 years. A first descriptive evaluation of the methylation status was performed using the Principal Component Analysis (Supplemental Figure 1). This analysis excluded the presence of samples with aberrant epigenetic profiles and the presence of clusters of samples carrying peculiar epigenetic features. Subjects were then grouped in 5 age ranges (0-19 years, 20-39 years, 40-59 years, 60-79 years, 80-106 years) including at least 10 subjects in each age range.

As shown in Figure 1A, global DNA methylation β levels progressively decreased with age and were inversely correlated with aging (β = −0.36; confidence interval = −0.48, −0.22; p-value = 1.43 × 10⁻⁶). Moreover, a significant progressive decrease in overall DNA methylation levels was also associated with BMI, as shown by the density plot in Figure 1B. Indeed, DNA methylation was inversely correlated with BMI also after adjustment for age (β = −0.19642; confidence interval: −0.37, −0.02; p-value < 0.05).

Figure 1. DNA Methylation profile distribution. The density plot describes the mean methylation profile of samples grouped by age range (A) and by BMI range (B).

The number of epimutations increases exponentially during aging

The distribution and variability of methylation levels were studied for each one of the 25.014 CpG sites using Box-and-whiskers plots (Figure 2) as described in the materials and methods section. This allowed to obtain the total number of SEMs for each subject.

Figure 2. Schematic description of the SEM. For each probe, a Box-and-whiskers plot analysis identifies extreme outliers samples. For each sample epimutations are detected, counted and mapped considering their genomic position.

As shown in Figure 3A, considering age ranges, the mean number of SEMs increased exponentially during aging.

Figure 3. Correlations between epimutations and aging. (A) Exponential relation between age ranges and number of SEM. (B) Linear correlation between age ranges and log(SEM), red squares indicates mean log(SEM) value. (C) Linear correlation between age log(SEM) considering all samples independently. (D) Density plot describing age dependent shift in log(SEM).

SEMs data were not normally distributed and thus were logarithmic transformed before regression analysis. A strong correlation between age ranges and log(SEM) was evident (η= 0.92) (Figure 3B). This correlation was still significant when the subjects were not grouped in age ranges but considered individually (r = 0.63; confidence interval: 0.53, 0.71; p-value = 2.2 × 10⁻¹⁶)(Figure 3C).

The density plot in Figure 3D shows the age dependent shift in SEMs frequency depending on age range.

Multivariate analysis performed considering log(SEM), overall DNA methylation levels, BMI and aging confirmed that these findings were not influenced by the overall DNA methylation levels. Interestingly, no significant correlation between log(SEMs) and BMI even after adjustment for age and other covariates was detected.

The number of SEMs is correlated with skewing of X chromosome inactivation

Using the HUMARA method, skewing of X chromosome inactivation (XCI) was evaluated in a subset of 45 heterozygous women. Multivariate analysis indicated a significant correlation between log(SEM) and degree of XCI after adjustment for age (β = 0.41; confidence interval: 0.14, 0.68; p-value = 0.0053).

Subsequently, we performed a PATH analysis using the not imputed original dataset to test the complete model containing the variables skewing of XCI, age, log(SEMs) and overall CpG methylation in order to understand the relations between them. After adjusting for the number of epimutations, the marginal correlation between skewing of XCI and aging was not significant yet. Figure 4 shows the path with the respective standardized coefficients. Similar results were obtained also using and combining the five imputed datasets.

Figure 4. Path diagram with associated path coefficients (β). Arrows indicate the interrelationships tested in the analysis: red arrows indicate significant relationships (p value < 0.05), black arrows indicate non-significant relationships AVG beta = average DNA methylation β values.

Study population

We enrolled a total number of 178 subjects spanning from 3 to 106 years selected from Istituto Auxologico Italiano Biobank. In order to better estimate epimutations frequency we identified 5 age ranges (0-19 years, 20-39 years, 40-59 years, 60-79 years, 80-106 years), enrolling at least 10 subjects for each age range. An histogram of age and BMI distribution of subjects is shown in Supplemental Figure 2.

Subjects with imprinting disorders, emathological cancers, immune diseases were excluded from the study. The study protocol was approved by the Ethical Committee of the Istituto Auxologico Italiano.

DNA extraction and bisulphite treatment of the DNA

Genomic DNA was extracted from peripheral blood using the Wizard genomic DNA purification kit (PROMEGA, Madison WI, USA). Samples were processed for DNA isolation as follows: 3 to 5 mL of blood were lysed with lysis solution provided, digested with proteinase K in sodium dodecyl sulfate (SDS) buffer at 37°C for 1 hours, then DNA was extracted by salting-out, and resuspended in Tris-EDTA (TE) buffer. Sodium bisulphite conversion of DNA (500 ng) was performed by the EZ DNA Methylation-Gold Kit (Zymo Research Corporation, Orange, CA) according to the manufacturer's protocol.

Quality control and quantification of DNA were performed before and after bisulphite conversion. DNA was quantified with NanoDrop (NanoDrop Products Thermo Scientific Wilmington, DE) and quality was assessed by visualisation of genomic DNA on 1% agarose gel electrophoresis. Only DNA samples not fragmented were subsequently processed.

Genome-wide methylation analyses

Whole-genome methylation analysis of the DNA from 178 subjects was performed using the Illumina Infinium Methylation Platform (Illumina, San Diego, CA). 94 samples were analysed using the Infinium HumanMethylation 450K BeadChip while 84 samples were analysed using the 27K BeadChip. Those chips allow to assess the methylation level of 485,764 and 27,578 CpG sites over the entire genome, respectively.

The methylation profile was analysed according to the manufacturer's instructions using Illumina-supplied reagents and conditions. In brief, after bisulfite conversion, 250 ng of DNA were whole-genome amplified (WGA) and enzymatically fragmented. The bisulfite-converted WGA-DNA samples were purified and randomized again prior to hybridization to BeadChips. During hybridization, the WGA-DNA molecules annealed to locus-specific DNA oligomers linked to individual bead types, one designed against the unmethylated site and one against the methylated site. After hybridization, allele-specific single-base extension provided another level of specificity and incorporated a fluorescent label for detection. The level of methylation was determined at each locus by the intensity of the two possible fluorescent signals, specific for the methylated and unmethylated alleles. DNA methylation values, described as β-values, were recorded for each locus in each sample. Beta-value provides a continuous measure of levels of DNA methylation at a CpG site, ranging from 0 in the case of completely unmethylated sites to 1 in completely methylated sites.

HUMARA assay

X-inactivation pattern was determined in 45 female blood DNA samples using a modification of the methylation analysis of the HUMARA locus as described previously (26-GENTILINI). Briefly, 250 ng of DNA were digested at 37°C for 2 h with 10 U HpaII and 10 U HhaI (CELBIO) and a no-enzyme control digest was also setup for each sample. Digested and undigested DNAs were then amplified in duplicate PCRs using primers, amplifying the highly polymorphic CAG repeat region in Exon 1 of the AR gene at Xq12. The sequences of the primers used were as follows: forward; 5′-GCT GTG AAG GTT GCT GTT CCT CAT-3′ labelled with 5′-phosphoamidite dye, reverse; 5′-TCC AGA ATC TGT TCC AGA GCG TGC-3′. The samples were amplified for 35 cycles comprising of 15 s at 95°C, 30 s at 62°C and 30 s at 72°C with an initial denaturation at 95°C for 5 min. The PCR products were separated on an ABI 310 automated sequencer. The size of PCR product from each allele was analysed by Genescan software for the quantification by peak height.

XCI skewing was measured as Degree of Skewing (DS), which designates the percentage of the preferentially active allele. DS varies between 0% and 50%, where 0% indicates a random X inactivation pattern and 50% a completely skewed inactivation pattern.

Methylation data management, normalisation and quality control

Methylation raw data were generated using GenomeStudio software (Illumina, San Diego, CA)and were pre-processed using IMA R package [23]. Samples with low bisulphite conversion (BS) efficiency (BS control intensity values < 4000) were immediately excluded from the analysis as well as samples that failed the quality control analysis performed according the IMA package pipeline. Moreover, a principal component analysis was performed on methylation data in order to identify and remove samples with aberrant methylation profiles. A total of 8 samples were excluded from the study. Density plot of β-values showed a bimodal distribution with a shift between the positions of the peaks derived from type I and type II assay. Therefore, peak correction was performed using the IMA package, moreover data were normalised using background correction followed by quantile normalisation method. Furthermore, probes containing missing β-values, probes having <95% of samples with detection p-value <0.05 and probes on chromosome X and Y were also removed from the analysis.

We finally pooled data obtained from 27K and from 450K arrays and created a new data matrix containing 25.014 probes that were in common between the two array platforms and passed quality control steps.

Epimutations detection

The distribution and variability of methylation levels in our population were studied for each one of the 25.014 CpG sites using Box-and-whiskers plots in order to identify SEMs. For each probe, whenever the methylation level of one subject extremely differed from the rest of the population we considered the outlier sample as epimutated for that locus.

Thus, for each locus epimutated subjects were identified as the extreme outliers with methylation level exceeding three times interquartile ranges Q1-(3 × IQR) and Q3+(3 × IQR). Finally, all epimutated loci were annotated in a new data matrix that allowed to calculate, for each subject, the total amount of epimutations and their genomic position.

A schematic description of the analysis is shown in Figure 2.

The Box-and-whiskers plot analysis was conducted using boxplot function provided in the R car package and confirmed using the outlier function in the R outliers package.

Validation of the SEM analysis

In order to confirm the power of this analytical approach to detect epimutations we performed two separate tests after introducing positive controls.

1. We analysed 3 samples in duplicates and compared epimutations found in each of them. Results showed a mean correlation of 0.99 p< 0.01 among the experiments. It's important to underlie that the duplicate samples underwent independent bisulfite conversion reactions and this suggests that epimutations are not significantly influenced by bisulfite conversion errors.

2. We analysed 48 whole blood DNA samples obtained from subjects affected by imprinting diseases (BWS syndrome, Angelman Syndrome and Silver Russel syndrome) that attended diagnostic test at Istituto Auxologico Italiano. For these subjects we already had a medical report indicating the genomic position of their epigenetic alteration. For this control test, we used only data obtained from Illumina Infinium 450K arrays thus reducing the number of the control population (n=91). However, using a more informative array, we were able to better define all probes showing epigenetic alterations. Briefly, after the identification of the epimutated probes we performed a test for overrepresentation of epimutated probes inside each gene using the hypergeometric cumulative function. We considered the total number of probes (N_G = 485.764), the total number of epimutated probes in the selected sample (N_S) the number of probes in each gene (G_I) and the epimutated probes in each gene (S_I).

The analysis identified genes with enriched number of epimutated probes (bonferroni corrected p-value < 0.05) confirming the presence of the epigenetic alterations previously reported in the medical report.

Supplemental Figure 3 shows positive control test results.

All the R scripts employed as well as the set of methylation data generated will be made freely available to the scientific community.

Multiple regression and path analysis

The variables (age, BMI, overall methylation level, Skewing of XCI and log(SEM))were Z transformed and Z scores were used for the multiple regression and path analysis. Results of regression analyses were indicated as standardised regression coefficients β.

The univariate and multivariate linear regressions were conducted using the generalised linear mixed model function provided in the R base package.

The interrelationships among skewing of XCI and the other variables were examined by a path analysis model. A path diagram with associated path coefficients (β) was constructed as shown in Figure 4 based on previous findings [13] and theoretic rationales.

As age is generally considered a determinant of changes in XCI skewing profile and in overall DNA methylation, direct and indirect Paths from age to these components were tested. The Path Analysis was performed using the General Structural Equation Models R package.

Missing data imputation

Multiple imputation, producing five imputed datasets, was carried out to allow full use of all available data.

Imputation was performed using the bootstrap-based EMB algorithm included in the R package Amelia II. Five imputations were generated and models fitted to each imputed dataset. Model results were consolidated using Rubin's rules [43].