Human body epigenome maps reveal noncanonical DNA methylation variation

Abstract

Understanding the diversity of human tissues is fundamental to disease and requires linking genetic information, which is identical in most of an individual's cells, with epigenetic mechanisms that could have tissue-specific roles. Surveys of DNA methylation in human tissues have established a complex landscape including both tissue-specific and invariant methylation patterns. Here we report high coverage methylomes that catalogue cytosine methylation in all contexts for the major human organ systems, integrated with matched transcriptomes and genomic sequence. By combining these diverse data types with each individuals' phased genome, we identified widespread tissue-specific differential CG methylation (mCG), partially methylated domains, allele-specific methylation and transcription, and the unexpected presence of non-CG methylation (mCH) in almost all human tissues. mCH correlated with tissue-specific functions, and using this mark, we made novel predictions of genes that escape X-chromosome inactivation in specific tissues. Overall, DNA methylation in several genomic contexts varies substantially among human tissues.

Extended Data Figure 1. Identification of differentially methylation regions (DMRs) and Multidimensional Scaling Analysis

a, Line plot showing the fraction of differentially methylated CG sites (DMSs, dynamic CGs) out of all CG sites under various methylation difference cutoffs. The methylation difference of a CG site is defined in Ziller et al. b, A plot of the first two principal components from the methylation level multi-dimensional scaling. Tissues are shaded by the organ group they belong to as in Figure 1c and 1d. c-d, Bar charts of the cumulative amount of variance explained by the first N principal components from the multi-dimensional scaling performed on the methylation levels of all DMRs (c) and the expression levels of all differentially expressed genes (d). e, A representative example of enriched GO biological process terms based on the most hypomethylated DMRs from LV-1. f, A representative example of enriched mouse phenotype terms based on the most hypomethylated DMRs from LV-1.

Extended Data Figure 2. DMRs and their correlation with transcription

a, A browser screenshot of an example DMR downstream of the TSS. b, Expression level of the BIN1 gene which contains the DMR in (a). c, The percentages of hypomethylated intragenic DMRs in each class of genomic features. c-h, Histone modification profiles of five categories of uiDMRs.

Extended Data Figure 3. Classification of uiDMR histone profiles and uiDMR properties

a, heatmap of the histone modification profiles for the five types of uiDMRs. The profiles were plotted for each mark across the DMR and the 5kb upstream and downstream and the colors of each cell indicate the input normalized ChIP-seq RPKM. The colors on the left indicate the group of each profile assigned by k-means clustering (red, weak enhancer; orange, promoter-proximal; green, transcribed; blue, unmarked; black poised enhancer). b, A pie chart of the distribution of uiDMRs across the classes defined by k-means clustering.

Extended Data Figure 4. Classification of promoter histone profiles

a, A heatmap of the histone modification profiles across strong (rows labeled with red) and unmarked (rows labeled with orange) promoters. The profiles were plotted for each mark across the promoter and the 5kb upstream and downstream and the colors of each cell indicate the input normalized ChIP-seq RPKM. b-c, The aggregate profiles for strong and unmarked promoters (b) and (c), respectively. d, The distribution of the Spearman correlation coefficients between the methylation level of different types of hypomethylated intragenic DMRs and the expression of the nearest gene. Notches indicate a confidence interval estimated from 1,000 bootstrap samples.

Extended Data Figure 5. uiDMR fetal DNase I profiles

DNase I profiles of various fetal tissues corresponding to the tissues presented in this study. The samples are arranged columnwise by age, and row-wise by fetal tissue. The uiDMR – unmarked line represents the DNase I profile of uiDMRs without histone modifications. The DMR – enhancer line represents the DNase I profile of DMRs that overlapped an enhancer in a matched tissue in this study (indicated in the row label in parentheses). The shuffled line represents the DNase I profile of uiDMRs randomly shuffled across the genome.

Extended Data Figure 6. PMD Features

a, A browser screenshot (see Figure 1 for description) of an example PMD found in IMR90, PLA, PA-2, and PA-3. RV-1 is included as a representative sample without PMDs. b, The distribution of sizes of PMDs in various samples. c, A heatmap representation of the overlap between various sets of PMDs. The denominator of the fraction of overlap is determined by the sample on the y-axis. d-e, ChIP-seq profiles of the PMD regions defined in PA-2 (c) and IMR90 (d) after shuffling.

Extended Data Figure 7. DNMT expression across tissues

a-d, Bar plots of the expression (measured in log10 FPKMs) of DNMT1 (a), DNMT3A (b), DNMT3B (c), and DNMT3L (d) across various samples.

Extended Data Figure 8. mCH distribution and correlation

a, A browser screenshot (see Figure 1 for description) of an example region with non-CG methylation (mCH). Purple and pink ticks are methylated CHG and CHH sites, respectively (H = A, C, or T). Ticks on the forward strand are projected upward from the dotted line and ticks on the reverse strand are projected downward. b, The distribution of methylation levels at mCH sites across all samples with a discernible TNCAC motif. Only mCH sites with at least 10 reads and a significant amount of methylation were considered. c, Boxplots of the expression values across different quantiles of CAC gene body methylation (Gene body mCAC). d, Scatterplot of mCAG vs. mCAC inside gene bodies. e, Bar plot of the correlation of mCAG and mCAC inside gene bodies (blue) and the theoretical maximal correlation (red) if mCAC and mCAG are perfectly correlated. f-h, The methylation levels of C (upper panel), CG (middle panel) and CH (lower panel) across the read positions for PO-2 (red line) and EG-3 (blue line). Vertical lines indicate the position (10^th base from the beginning) where trimming was applied. i, mCH motif from PO-2 with the first 10 bases of each read trimmed. j, mCH motif from PO-2 without trimming. k, mCH motif from EG-3 with the first 10 bases of each read trimmed l, mCH motif from EG-3 without trimming. The height of each letter represents its information content (i.e., prevalence).

Extended Data Figure 9. X chromosome inactivation

a, Distributions of promoter CG methylation (mCG) levels (mCG/CG), gene body non-CG methylation (mCH) levels (mCH/CH), gene body mCG levels and promoter mCH levels in genes previously reported to express from only one allele (inactivated) or biallelically (escapee)⁶³. Black ticks show median, and bars indicate 25-75^th percentile range. Genes more prone to escaping inactivation have lower promoter mCG, higher gene body mCH, higher gene body mCG and higher promoter mCH in females. b-e, Discriminability analysis using b, gender-specific gene-body mCH, c, promoter mCG, d, promoter mCH and e, gene body mCG to predict the escapee status of X-linked gene, respectively. Among them, gene body mCH is the most predictive feature of chromosome X inactivation escapees.

Extended Data Figure 10. Allele-specific Methylation and Expression

a, An example of allele-specific methylation (ASM). Reads that contain a heterozygous SNP (red box) are separated by allele. The number of methylated (reads containing Cs) and unmethylated (reads containing Ts) at adjacent CG sites (black boxes) and tested for differential methylation. b, Fraction of allele-specific expressed (ASE) genes (blue) and bi-allelically expressed genes (grey) that have at least one ASM event within a certain distance. Bi-allelically expressed genes were defined as genes that were covered by at least 10 reads and whose p-values given by binomial test for allelic expression were greater than 0.2 (i.e. no significance). c, Fraction of ASE genes that were linked to matched ASM event(s) (blue) and matched ASM events with their locations shuffled (grey). b-c are aggregated results using samples from triplicate tissues.

Figure 1. The methylomes and transcriptomes of human tissues

a, The tissues analyzed in this study. Samples are denoted by the two letter code in parentheses followed by an individual ID. b, Browser screenshot of an example DMR. The top track contains gene models. The following four tracks contain green blocks indicating the location of super enhancers, enhancers, and hypomethylated DMRs in aorta, respectively. The remaining tracks display methylation data from each sample. Gold ticks are CG sites with heights proportional to their methylation level. Ticks on the forward and reverse strand are projected upward and downward from the dotted line, respectively. c-d, Hierarchical clustering of DMR methylation levels (c) and expression levels of differentially expressed genes (d). Colors indicate organ systems each sample belongs to.

Figure 2. DNA methylation and its relationship with gene expression

a, The mean Spearman correlation coefficient at various distances between the methylation level of autosomal DMRs and the expression of the nearest gene. These correlations are shown for DMRs: overlapping genes (Genebody), overlapping enhancers (Enhancer), overlapping promoters or CpG islands (CGIs) or CGI shores (Promoter, CGI, CGI shore), not overlapping genes (Intergenic) and all remaining DMRs (Undefined). b, Heatmap showing each motif’s tissue-specific methylation preference. The tissues are colored according to Fig 1c., and the ordering is listed at the bottom of the figure. The bar plot at the end of the panel shows the number of times the motif was present in the 20 motif models. c, The number of base pairs covered by PMDs in all samples. d, The distribution of expression inside and outside of PA-2 PMDs across various samples. Notches indicate a confidence interval estimated from 1,000 bootstrap samples. Each PMD boxplot consists of 3,627 genes and each non-PMD boxplot consists of 22,907 genes. e-f, Histone modification profiles in and around PMDs in PA-2 (e) and IMR90 (f).

Figure 3. mCH is prevalent in human tissues

a, The fraction of methylated cytosines in the CH context by sample. b-d, Representative mCH motifs from embryonic, (H1; b), tissue (LI-11; c), and brain (NRN; d) samples. The height of each letter represents its information content. e, A heatmap of genic mCAS patterns normalized to the flanking region. Each gene was assigned to one of twenty clusters, which is indicated by the number and tick marks on the y-axis. The tick marks on the x-axis indicate the upstream, transcription start, transcription end, and downstream segments of each gene. The boxes around various patterns highlight regions referenced in the main text. f, Bar plot of the ratio of the genome-wide mCAC to mCAG in various samples.

Figure 4. Allele-specific Methylation and Expression

a, Browser screenshot of the increase in female mCH for a gene known to escape X chromosome inactivation (MED14). Sample names are colored by gender (male, black; female, red). b, Ratio of mCH level in female vs. male samples across genes with a significant difference in at least one sample. Cells boxed in black denote samples with a statistically significant difference between females and males. c, The number of ASM and ASE sites across the triplicated tissues. The top row depicts ASM events (left) and ASE events (right) which are allele-specific in all tissues (black), are variable across tissues (white), or do not possess enough data to tell (grey). The bottom row depicts the distribution of variable sites from the top row that vary by individual (white), tissue (black), or neither (grey).