DNA methylation loss in late-replicating domains is linked to mitotic cell division

Abstract

DNA methylation loss occurs frequently in cancer genomes, primarily within lamina-associated, late-replicating regions termed partially methylated domains (PMDs). We profiled 39 diverse primary tumors and 8 matched adjacent tissues using whole-genome bisulfite sequencing (WGBS) and analyzed them alongside 343 additional human and 206 mouse WGBS datasets. We identified a local CpG sequence context associated with preferential hypomethylation in PMDs. Analysis of CpGs in this context ('solo-WCGWs') identified previously undetected PMD hypomethylation in almost all healthy tissue types. PMD hypomethylation increased with age, beginning during fetal development, and appeared to track the accumulation of cell divisions. In cancer, PMD hypomethylation depth correlated with somatic mutation density and cell cycle gene expression, consistent with its reflection of mitotic history and suggesting its application as a mitotic clock. We propose that late replication leads to lifelong progressive methylation loss, which acts as a biomarker for cellular aging and which may contribute to oncogenesis.

Figure 1. Solo-WCGW CpGs are prone to hypomethylation

(a) Each genomic CpG dinucleotide was placed into one of four CpG density categories (0, 1, 2, or 3+, depending on the number of additional CpGs within a +/− 35 bp window), and one of the three flanking nucleotide categories (SCGS, SCGW and WCGW, with “S” being C:G and “W” being A:T). Because CpGs are palindromic, WCGS and SCGW were combined. Each of the 4×3=12 possible contexts are shown as columns for CpGs within common HMDs (left) or common PMDs (right). In the illustrations, a star indicates the target CpGs, and solid circles indicate all neighboring CpGs within the window. The number of CpGs in each context is shown as a percentage of all genomic CpGs; for instance, the first column shows that 6% of all CpGs in the human genome are within HMDs, have 3+ flanking CpGs, and SCGS tetranucleotide context. (b) Violin plots show beta value distributions for CpGs in each context, for five human tissues (two normal colon tissues and three colon tumors) and two mouse tissues (one normal colon tissue and one colon tumor). Violin color indicates mean beta value. Columns shaded orange and green indicate the most hypomethylation-resistant and most hypomethylation-prone categories, respectively. (c) Average methylation values (non-overlapping 100-kb bins) across a 12-mb section of chr16p, for the human colon samples. Values were calculated using all CpGs (left), only hypomethylation resistant CpGs (orange, middle), or only Solo-WCGW CpGs (green, right). CpG islands were removed in all analyses.

Figure 2. Most PMDs are shared across cancer and normal tissues

(a) Average methylation values (non-overlapping 100-kb bins) for chr16p, shown for the core tumor/normal dataset. The “tumor” field indicates tumors (black) vs. adjacent normals, and “this study” field indicates samples that were newly sequenced as part of this study (black). Within both normal and tumor classes, tissue types are grouped and ordered by average methylation level of samples from the group. For instance, “endometrium” is the first normal group because it has the highest methylation among normal groups, and likewise for “GBM” among tumor groups. (b) Average methylation across all normal (upper) or tumor samples (lower), calculated for multiple window sizes from 10 kb to 10 mb (“multi-scale plot”). (c) SD across all normal or tumor samples as multi-scale plots. (d) 100-kb SD values for the all non-overlapping genomic bins, plotted for tumors (red histogram, X-axis) vs. normals (blue histogram, Y-axis). Bimodal peaks for each were identified via a Gaussian mixture model, and cutoffs dividing low and high SD values are indicated by dashed lines for each axis. A scatter cloud shows the correlation between SD values between the tumors and normals, indicating the percentage of 100-kb bins falling into each of the four quadrants as well as Spearman’s ρ. (e) Illustration of method used to rescale each sample’s methylation values based on genome-wide levels within a common set of PMDs (Online Methods). (f) Same data as panel (a), but using rescaled methylation values.

Figure 3. Most PMDs are shared across developmental lineages

(a) Average solo-WCGW methylation levels were plotted along chromosome 16p for 390 WGBS samples, organized into 6 groups: Germline and preimplantation embryo (GE). Post-implantation embryonic/fetal samples (FT), grouped first by embryonic vs. extra-embryonic, then by average methylation. Cell lines (CL). Post-natal non-blood normal tissue samples (PN). Post-natal blood-derived samples (PB). Primary tumors (TM). Within each of the 6 groups, samples were organized by cell type (labeled with color codes). Lamin B1 signal and replication timing of IMR90 lung fibroblast are shown below methylation heatmaps (bottom). (b) Mean methylation levels within each of the 5 major groups (excluding group GE), plotted as in Fig. 2b. (c) SD within each of the 5 major groups, plotted as in Fig. 2c. (d) SDs shown for the 100-kb scale alone. (e) Distribution of SD for all non-overlapping 100-kb genomic bins across all samples of the core tumor group (from panel (d)) are plotted on the Y-axis, compared to each of four major groups (FT, CL, PN, and PB), shown on the X-axis. Group GE is omitted due to lack of PMD structure.

Figure 4. Most PMDs are shared across developmental lineages in mouse

Average solo-WCGW methylation levels were plotted along a representative 30-mb region of chromosome 17 in mouse. 206 WGBS samples are organized into four groups: Embryonic Stem Cells (ESC); Germline and embryos (GE); Fetal tissues (FT); Postnatal tissues (PN); Grouping and ordering of samples were performed as described in Fig. 3. Lamin and replication timing are shown on the bottom of the heatmap. Lamin A DamID from wild type mouse ESCs were downloaded from GEO with accession GSE62683. Replication timing of day 9 differentiated ESCs were downloaded from GEO with accession GSE17983.

Figure 5. PMD hypomethylation emerges during embryonic development

(a) Multi-scale solo-WCGW average plots are shown for samples divided into seven developmental stages, as diagrammed in (b): paternal (I) and maternal (II) germ cells, implantation-related tissues (III), primordial germ cells (IV), embryonic soma (V), fetal soma (VI) and postnatal soma (VII). (c) Rank-based analysis of the 792 genomic 100-kb bins from chr16, comparing methylation ranks of the core tumors (Y-axis) to each developmental sample (X-axis), with each axis going from a rank of 1 (lowest methylation) to the rank of the highest methylation (excluding bins with missing value from either of the samples). Greater correlations (indicated by the Spearman’s correlation coefficient ρ) indicated stronger HMD/PMD structure.

Figure 6. PMD hypomethylation is associated with chronological age

(a) Multi-scale solo-WCGW average plots are shown for newborn CD4 T cell, 103-year-old CD4 T cell (GSE31438) and T cell prolymphocytic Leukemia (BLUEPRINT accession S016KWU1). (b–f) Summarization of average PMD hypomethylation in HM450-based samples, by averaging beta values for 6,214 solo-WCGW probes mapped to common PMDs (Online Methods). Peripheral Blood Mononuclear Cell (PBMC) in newborns and nonagenarians (left, from GSE30870, p=8.8e-5, one-way Wilcoxon Rank Sum test), and disease-free fetal and adult liver tissue (right, from GSE61278). Center lines of the box plots indicate median, and the lower and upper bounds indicate lower and upper quartiles. The lower and upper whiskers indicate smallest and largest methylation values. ** p <= 0.001 from Wilcoxon Rank Sum test. (c–f) HM450-based solo-WCGW averages vs. age for individual donors for several tissue types. N is the number of donors/samples, r is Pearson’s product moment correlation, b1 is the estimated rate of methylation loss, and p is the p-value based on Pearson correlation test. (c) Four fetal tissue types during three pre-natal time points (from GSE56515). (d) Sun-exposed and sun-protected dermis and epidermis (from GSE51954). (e) Sorted blood cells of the myeloid lineage (D1: GSE35069; D2: GSE56046). (f) Sorted blood cells of lymphoid lineage (D1: GSE35069; D3: GSE71955; D4: GSE59065).

Figure 7. PMD hypomethylation is linked to mitotic cell division in cancer

(a) PMD-HMD solo-WCGW methylation difference for 9,072 tumors from TCGA HM450 data. Each sample is ordered within cancer type by PMD-HMD difference, and cancer types are ordered by average PMD-HMD difference. (b) PMD methylation (X-axis) vs. somatic mutation density (Y-axis) for all 3,959 high purity TCGA cases (purity>=0.7), with Spearman’s ρ indicated. The blue line represents the regression line for all samples, while the red regression line excludes “hypermutator” samples (Online Methods). (c) Density of somatic LINE-1 insertions (violin plot elements) in non-overlapping 1-mb genomic bins (N=3,053), stratified by percent of bin overlapping common PMDs (only cases with whole-genome sequencing are included). (d) PMD methylation (X-axis) vs. LINE-1 insertion counts (Y-axis) for nine TCGA cancer types having substantial LINE-1 insertion counts. * (ρ < 0.05) and ** (ρ <= 0.01) indicate Spearman’s test significance. (e) The 10 most significantly enriched Gene Ontology (GO) terms for the 60 genes with the most strongly correlated expression vs. PMD hypomethylation in TCGA tumors, showing fold enrichment (grey) and false discovery rate (olive). (f) Gene Set Enrichment Analysis (GSEA) for 350 cell-cycle-dependent genes from Cyclebase, ranking all genes according to degree of expression vs. PMD hypomethylation correlation. (g) Normalized expression (Z-scores) of cell-cycle-dependent genes from Cyclebase (categorized by cell cycle phase) in 3,414 high purity TCGA tumor samples (purity >= 0.7), ordered by PMD-HMD methylation difference.

Figure 8. Replication timing and H3K36me3 contribute independently to methylation maintenance

(a) Multi-scale plot of chr16p showing similarity between solo-WCGW methylation and other chromatin marks in the IMR90 fibroblast cell line. (b) Average methylation level of all genomic solo-WCGWs in IMR90, stratified by (1) overlap with H3K36me3 peaks (left vs. right), (2) context relative to gene annotations (“Genic” vs. “Intergenic”), and (3) Repli-seq replication timing bin (red, yellow, light blue, dark blue). For Solo-WCGWs residing within +/− 10 kb of an annotated gene (Genic), meta-gene plots show methylation averages in relation to the Transcription Start Site (TSS) and the Transcription Termination Site (TTS). For all other Solo-WCGWs (Intergenic), each replication timing group is shown as a single violin plot. (c) The same representation of data plotted for the H1 hESC cell line (using Repli-chip data rather than Repli-seq). (d) Schematic summary, showing Solo-WCGW CpG methylation loss primarily determined by replication timing domain but locally protected by H3K36me3. (e) Schematic model illustrating DNMT1 processivity favoring dense CpGs and leading to incomplete re-methylation of Solo CpGs. (f) Schematic illustration of the “re-methylation timing model” where genomic regions synthesized earlier in S-phase (HMDs) spend more time exposed to methylation maintenance machinery and thus more complete methylation maintenance than PMDs. (g) Illustration of the relationship between major determinants of hypomethylation and 3D nuclear topology, with Lamina Associated Domains (LADs) occupying a distinct heterochromatic nuclear compartment.