A genomic atlas of systemic interindividual epigenetic variation in humans

Abstract

Background: DNA methylation is thought to be an important determinant of human phenotypic variation, but its inherent cell type specificity has impeded progress on this question. At exceptional genomic regions, interindividual variation in DNA methylation occurs systemically. Like genetic variants, systemic interindividual epigenetic variants are stable, can influence phenotype, and can be assessed in any easily biopsiable DNA sample. We describe an unbiased screen for human genomic regions at which interindividual variation in DNA methylation is not tissue-specific.

Results: For each of 10 donors from the NIH Genotype-Tissue Expression (GTEx) program, CpG methylation is measured by deep whole-genome bisulfite sequencing of genomic DNA from tissues representing the three germ layer lineages: thyroid (endoderm), heart (mesoderm), and brain (ectoderm). We develop a computational algorithm to identify genomic regions at which interindividual variation in DNA methylation is consistent across all three lineages. This approach identifies 9926 correlated regions of systemic interindividual variation (CoRSIVs). These regions, comprising just 0.1% of the human genome, are inter-correlated over long genomic distances, associated with transposable elements and subtelomeric regions, conserved across diverse human ethnic groups, sensitive to periconceptional environment, and associated with genes implicated in a broad range of human disorders and phenotypes. CoRSIV methylation in one tissue can predict expression of associated genes in other tissues.

Conclusions: In addition to charting a previously unexplored molecular level of human individuality, this atlas of human CoRSIVs provides a resource for future population-based investigations into how interindividual epigenetic variation modulates risk of disease.

Fig. 1

Strategy for identifying correlated regions of systemic interindividual variation (CoRSIVs). a The tissues analyzed represent the three germ layer lineages; 10 Caucasian GTEx donors were studied, yielding 30 methylomes. b Initial unsupervised clustering of whole-genome bisulfite sequencing data; considering all informative bins, they cluster by tissue. c Example of a CoRSIV identified at the PM20D1 promoter. The blue triangle shows a region of correlated methylation comprising thirteen 100-bp bins; the three scatter plots illustrate its high inter-tissue correlation. d Plots of individual methylation at the PM20D1 CoRSIV illustrate systemic interindividual variation. Genotype data at rs708727 (bottom panel) indicate strong mQTL at the locus. e Scatter plot of interindividual methylation range vs. number of CpGs per CoRSIV, for all 39,424 CoRSIVs initially identified. Subsequent analyses focus on the 9926 CoRSIVs with ≥ 5 CpGs/CoRSIV and IIR ≥ 20 (shaded area). f Unlike genome-wide bins, the 9926 CoRSIV bins cluster by individual (B, H, T—brain, heart, thyroid). Box plots on right show that the 10 donors show no individual differences in average methylation across all the CoRSIVs. g An illustrative region from the CoRSIV plot of chr19. Inset shows example of annotation of a CoRSIV (chr19_8726) at ZNF714

Fig. 2

Distribution and characteristics of 9926 human CoRSIVs. a Scatter plot of CoRSIV size vs. CpGs per CoRSIV. b Circos plot of human autosomes indicates regions of high CoRSIV density. c Individual methylation levels across examples of CoRSIVs from each of the two most high-density regions: (left) HLA-C (chromosome 6) and (right) DUX4L35 (chromosome 20). d Compared to all informative bins, CoRSIVs are more than twofold enriched in sub-telomeric regions (χ² test P < 10^− 300). e Compared to control regions or tissue-specific differentially methylated regions (tDMRs), CoRSIVs are enriched for repetitive elements and depleted for CpG islands (CGI) and transcription binding sites (TFbs) (χ² test P < 10^− 8 for all comparisons). f Relative to control and tDMR regions, CoRSIVs are under-represented at transcription start sites (TSS), within gene bodies, and at transcription end sites (TES) (χ² test P < 10^− 16 for all comparisons). g ChromHMM epigenomic features significantly enriched or depleted in CoRSIVs

Fig. 3

Evaluating associations between CoRSIV DNA methylation and neighboring genetic variation. a CoRSIV-specific P value distributions from mQTL analyses at the two most CoRSIV-dense genomic regions. At the MHC locus on chromosome 6, most CoRSIVs show mQTL; none of 10 informative CoRSIVs in the pericentromeric region on chromosome 20 show mQTL. b Summary of pair-wise analyses. Scatter plot shows correlation of interindividual differences in CoRSIV methylation vs. interindividual differences in average CoRSIV-associated SNV genotype for one pair of GTEx donors (RVPV and X261). Distribution of all such pairwise comparisons (right); average R² = 0.61. c Heat map of published mQTL effects at 819 CoRSIVs overlapping probes on the HM450 array. Approximately 40% (white) show no evidence of mQTL. Dotted line shows location of probes within the PM20D1 CoRSIV. d Histogram of P values for the spearman correlation between methylation and read coverage within each CoRSIV. Two insets above the histogram show two scatter plots drawn from opposite ends of the histogram. The top left inset shows a strong association (P < 0.05) between average methylation in chr3_2398 CoRSIV and individual-average read depth in this region. The top right inset shows weaker association (P > 0.85) between average methylation in chr6_24299 CoRSIV and individual-average read depth in this region. e Decay curves of correlation in methylation and LD for CoRSIVs and control regions. f Histogram of CoRSIV-specific P values for the Spearman correlation analyses of methylation decay vs. LD decay. Locations of two CoRSIVs with (ITGB2) and without a significant association (AC07910.1) are indicated. g The CoRSIV at ITGB2 yielded a significant association (P < 0.05) between methylation decay and LD decay, consistent with a haplotype-effect on methylation. The scatter plot on the left (top, blue) shows all methylation R² for pairs of bins in which one bin is within the CoRSIV and other bin is within ± 20 kb. The scatter plot on the left (bottom, red) shows LD for all pairs of SNVs. For each pair one SNV is within the CoRSIV and the other is within 20kb. On the right, the upper diagonal shows the CoRSIV (blue) in this region, and the lower diagonal shows the LD pattern in this region (1000 Genomes CEU). h In contrast, the CoRSIV at AC07910.1 showed no association between methylation decay and LD decay (P > 0.9), suggesting pure epigenetic variation

Fig. 4

CoRSIVs show influence of periconceptional environment, correlate with gene expression, and are associated with human disease. a Seasonal variation of blood methylation at CoRSIVs, controls, and tDMRs in 233 Gambian 2-year olds [13]. Data represent CpGs showing significant seasonal variation (FDR < 20%, see “Methods”). Predicted methylation maxima for CpGs within CoRSIVs occur in conceptions at the peak of the Gambian rainy season (July–September); minima fall within the Gambian dry season (Jan–April). Seasonal patterns at control regions and tDMRs are less pronounced (Fisher’s exact test P > 0.8). b Results of analysis linking gene expression in adipose tissue, skin, and lymphoblastoid cell lines (LCL) vs. methylation in adipose tissue [26] for 645 gene-associated CoRSIVs that are informative on the HM450 array. Venn diagrams show that for most genes that show a significant association (Spearman P < 0.05) between methylation and expression in adipose tissue, methylation in adipose tissue is also associated with expression in skin or LCL. c Examples of expression vs. CoRSIV methylation data in all three tissues at the MRI1 promoter (left), the CNDP2 gene body (middle), and the 3′ end of SSNA1 (right). d Summary of an automated PubMed literature-search using PubTator. Shown are MESH code labels corresponding to the top 100 human diseases linked to CoRSIV-associated genes