A methylome-wide study of aging using massively parallel sequencing of the methyl-CpG-enriched genomic fraction from blood in over 700 subjects

Abstract

The central importance of epigenetics to the aging process is increasingly being recognized. Here we perform a methylome-wide association study (MWAS) of aging in whole blood DNA from 718 individuals, aged 25-92 years (mean = 55). We sequenced the methyl-CpG-enriched genomic DNA fraction, averaging 67.3 million reads per subject, to obtain methylation measurements for the ∼27 million autosomal CpGs in the human genome. Following extensive quality control, we adaptively combined methylation measures for neighboring, highly-correlated CpGs into 4 344 016 CpG blocks with which we performed association testing. Eleven age-associated differentially methylated regions (DMRs) passed Bonferroni correction (P-value < 1.15 × 10(-8)). Top findings replicated in an independent sample set of 558 subjects using pyrosequencing of bisulfite-converted DNA (min P-value < 10(-30)). To examine biological themes, we selected 70 DMRs with false discovery rate of <0.1. Of these, 42 showed hypomethylation and 28 showed hypermethylation with age. Hypermethylated DMRs were more likely to overlap with CpG islands and shores. Hypomethylated DMRs were more likely to be in regions associated with polycomb/regulatory proteins (e.g. EZH2) or histone modifications H3K27ac, H3K4m1, H3K4m2, H3K4m3 and H3K9ac. Among genes implicated by the top DMRs were protocadherins, homeobox genes, MAPKs and ryanodine receptors. Several of our DMRs are at genes with potential relevance for age-related disease. This study successfully demonstrates the application of next-generation sequencing to MWAS, by interrogating a large proportion of the methylome and returning potentially novel age DMRs, in addition to replicating several loci implicated in previous studies using microarrays.

Figure 1.

Flow diagram of data processing and quality control (QC). MBD-seq, methyl-CpG binding domain (MBD) protein-enriched sequencing; PCA, principal component analysis; MWAS, methylome-wide association study.

Figure 2.

Quantile–quantile (QQ) plot of observed CpG block association P-values (y-axis) against expected P-values (x-axis). The negative logarithm in base 10 of the association P-value is plotted. The red line depicts the expectation under the null hypothesis. The deviation of many points above the 95% confidence interval (blue lines) in the right upper corner points to a considerable number of significant findings.

Figure 3.

Circular plot with MWAS results and pathways/networks. Each dark-red spot in outermost track with gray background is −log10(MWAS P-value). Genomic features associated with these sites are shown by repeats (dark-purple), exons (green) and introns (orange), respectively, in the concentric circular tracks between the P-values and chromosome number. Of the 68 unique genes associated (±20 kb flaking region) with the top 70 DMRs (q-value < 0.1), the most relevant 45 genes for pathways and networks are shown. Links between genes indicate co-occurrence in pathways and protein–protein interactions networks (Table 3). Lines connecting genes in deep-purple represent protein–protein interactions. Genes co-occurring in reference pathways are linked as follows: c-Jun kinases (light blue), Downstream TCR signaling (deep blue), p38 MAPK signaling pathway (light green), MAPK signaling pathway (green), protein processing in endoplasmic reticulum (pink) and amyotrophic lateral sclerosis, ALS (red).

Figure 4.

Bioinformatics plot of genomic features for the top sites. The horizontal axis shows the proportion of features (vertical axis) overlapping with the top sites (q-values < 0.1). Three groups of points are plotted: red circles corresponding to all 70 top DMRs, blue squares for hypermethylated age-DMRs and green triangles for hypomethylated age-DMRs. Where enrichment compared with the genome-wide average is significant (P < 0.05), points are shown as solid color blocks. In panel A, ‘Exon’, ‘Intron’ and ‘Gene’ designate overlap with RefSeq genes; ‘CGI’ denotes overlap with a CpG island; ‘Shore’ is ±2 kb flanking a CGI; ‘Upstream 2 kb/8kb’ indicates within 2 or 8 kb upstream of transcription start site; ‘DNase cluster’ indicates a genomic region hypersensitive to DNaseI; ‘Conserv’ indicates regions of high conservation across eutherian mammals; ‘TFBS’ indicates conserved motifs for transcription factor-binding sites in humans and rodents. In panel B, ‘CTCF’ is CCCTC-binding factor protein; ‘EZH2’ is enhancer of zeste homolog 2, a histone-lysine N-methyltransferase; ‘H2AZ’ is H2A histone family, member Z, whereas the remaining categories indicate histone modifications described using standard nomenclature.