Links between DNA methylation and nucleosome occupancy in the human genome

Abstract

Background: DNA methylation is an epigenetic modification that is enriched in heterochromatin but depleted at active promoters and enhancers. However, the debate on whether or not DNA methylation is a reliable indicator of high nucleosome occupancy has not been settled. For example, the methylation levels of DNA flanking CTCF sites are higher in linker DNA than in nucleosomal DNA, while other studies have shown that the nucleosome core is the preferred site of methylation. In this study, we make progress toward understanding these conflicting phenomena by implementing a bioinformatics approach that combines MNase-seq and NOMe-seq data and by comprehensively profiling DNA methylation and nucleosome occupancy throughout the human genome.

Results: The results demonstrated that increasing methylated CpG density is correlated with nucleosome occupancy in the total genome and within nearly all subgenomic regions. Features with elevated methylated CpG density such as exons, SINE-Alu sequences, H3K36-trimethylated peaks, and methylated CpG islands are among the highest nucleosome occupied elements in the genome, while some of the lowest occupancies are displayed by unmethylated CpG islands and unmethylated transcription factor binding sites. Additionally, outside of CpG islands, the density of CpGs within nucleosomes was shown to be important for the nucleosomal location of DNA methylation with low CpG frequencies favoring linker methylation and high CpG frequencies favoring core particle methylation. Prominent exceptions to the correlations between methylated CpG density and nucleosome occupancy include CpG islands marked by H3K27me3 and CpG-poor heterochromatin marked by H3K9me3, and these modifications, along with DNA methylation, distinguish the major silencing mechanisms of the human epigenome.

Conclusions: Thus, the relationship between DNA methylation and nucleosome occupancy is influenced by the density of methylated CpG dinucleotides and by other epigenomic components in chromatin.

Keywords: CpG; DNA methylation; Epigenetics; MNase-seq; NOMe-seq; Nucleosome.

Fig. 1

Genome-wide nucleosomal DNA sequence and methylation patterns in leukocytes and IMR90 cells. a For the in vitro MNase-seq data from leukocytes, average occurrences of select dinucleotides were computed from forward and reverse complement sequences aligned to nucleosome midpoints. Using BS-seq and MNase-seq data from leukocytes, average occurrences of CpGs (b), methylated and unmethylated CpGs (c, d) and mCpG/CpG fractions (e, f) were computed for both the in vitro and in vivo libraries from forward and reverse complement sequences aligned to nucleosome midpoints. Using MNase-seq and NOMe-seq data from IMR90 cells, average mHCG/HCG fractions (g, h) and uGCH/GCH fractions (i, j) were computed from forward and reverse complement sequences aligned to MNase-seq derived nucleosome midpoints. Nucleosomes within CpG Islands were excluded from this analysis

Fig. 2

Average DNA methylation levels and nucleosome occupancy as a function of nucleosome core CpG frequency. Nucleosome sublibraries were generated for the in vitro and in vivo leukocyte and IMR90 MNase-seq libraries based on the number of CpG occurrences between positions −61 to +61 relative to nucleosome midpoints. Using BS-seq and MNase-seq data from leukocytes, average mCpG/CpG fractions were computed from forward and reverse complement sequences aligned to nucleosome midpoints for each in vitro (a) and in vivo (b) nucleosome sublibrary. Using MNase-seq and NOMe-seq data from IMR90 cells, average mHCG/HCG fractions (c) and uGCH/GCH fractions (d) were computed from forward and reverse complement sequences aligned to MNase-seq derived nucleosome midpoints for each nucleosome sublibrary. Nucleosomes within CpG Islands were excluded from this analysis

Fig. 3

Effects of increasing methylated CpG density in the nucleosome core on nucleosome occupancy and the ratio of methylated CpG density in core versus linker DNA. Boxplots display the distribution of nucleosome occupancy as a function of the number of CpGs (a), the number of methylated CpGs (c), the number of unmethylated CpGs (d), G + C content (g), and G + C-CpG content (i) for all nucleosomes in the genome outside of CpG islands. Boxplots also show the distribution of DNA methylation levels (b), G + C content (f), and G + C-CpG content (h) as a function of the number of CpGs. e Boxplots show the distribution of the ratios of methylated CpG density per base pair in the core versus linker DNA as a function of the number of CpGs. See “Methods” section for more details

Fig. 4

Frequency profiles of CpGs, DNA methylation levels, and nucleosome occupancy surrounding nucleosomes positioned within exons, introns, and SINE-Alu elements. Using MNase-seq and NOMe-seq data from IMR90 cells, average occurrences of CpGs (a), mHCG/HCG fractions (b) and uGCH/GCH fractions (c) were computed from forward and reverse complement sequences aligned to MNase-seq derived nucleosome midpoints

Fig. 5

Analysis of nucleosome occupancy and mCpG density in differentially marked chromatin across the genome. Using MNase-seq, NOMe-seq, and Roadmap ChIP-seq data from IMR90 cells, average occurrences of CpGs (a), mHCG/HCG fractions (b), and uGCH/GCH fractions (c) were computed from forward and reverse complement sequences aligned to MNase-seq derived nucleosome midpoints within several histone modification peaks indicated in the key. The 12 histone modifications (and variant H2A.Z) in the key are ordered by decreasing average DNA methylation at the nucleosome midpoint. d Using data between positions −61 to +61 relative to the nucleosome midpoint, the average nucleosome occupancy for each histone modification was plotted against the corresponding average number of methylated CpGs in the central 123 base pairs of the nucleosome. e–g Nucleosome occupancy (e) and mCpGs density (f) were weighted by ChIP-seq Z-scores for both single and paired histone modifications (see “Methods” section). g Plots of nucleosome occupancy and mCpG density from the heatmaps in e and f show that H3K9me3-modified chromatin (blue dots) deviates from the correlation observed in the genome-wide data

Fig. 6

Characterization of chromatin at CpG islands. a Boxplots display the distribution of nucleosome occupancy as a function of the number of methylated CpGs. b, c Nucleosome with midpoints within CpG islands was divided into two groups, non-TSS and TSS. For each nucleosome, nucleosome occupancy and DNA methylation values were plotted in the 2D color-coded scatterplots. Using NOMe-seq and Roadmap ChIP-seq data from IMR90 cells, nucleosome occupancy, DNA methylation levels, and Z-scores from 12 histone modifications and the histone variant H2A.Z were aligned to TSSs that overlapped CpG Islands (d) and to the centers of CpG islands located in gene bodies and intergenic regions (e). Heatmap values were computed in 101 21-bp bins surrounding the TSSs and non-TSS CpG island centers, and sites (heatmap rows) were sorted by nucleosome occupancy

Fig. 7

Characterization of chromatin at unmethylated and methylated conserved transcription factor binding sites. All conserved transcription factor binding sites (TFBSs) from the UCSC genome browser were divided into unmethylated and methylated sites (see “Methods” section). a Using NOMe-seq, BS-seq, and ENCODE ChIP-seq data from HCT116 cells, nucleosome occupancy, DNA methylation levels, and DNase-seq and H3K27ac signals were aligned to unmethylated and methylated conserved TFBSs. b Average uGCH/GCH fractions were also computed from forward and reverse complement sequences aligned to these unmethylated and methylated conserved TFBSs. c Similar to a, nucleosome occupancy, DNA methylation levels, and DNase-seq, H3K27ac, and SP1 signals were aligned to unmethylated and methylated SP1 conserved TFBSs. d Using ENCODE ChIP-seq data for 10 transcription factors, including SP1, average signals were plotted with respect to their corresponding unmethylated and methylated conserved TFBSs. In the heatmaps, average signals were computed in 51 21-bp bins surrounding each conserved site. Scales for DNase-seq, H3K27ac, and the transcription factors represent occupancy levels in reads per million