Methylation QTLs are associated with coordinated changes in transcription factor binding, histone modifications, and gene expression levels

Abstract

DNA methylation is an important epigenetic regulator of gene expression. Recent studies have revealed widespread associations between genetic variation and methylation levels. However, the mechanistic links between genetic variation and methylation remain unclear. To begin addressing this gap, we collected methylation data at ∼300,000 loci in lymphoblastoid cell lines (LCLs) from 64 HapMap Yoruba individuals, and genome-wide bisulfite sequence data in ten of these individuals. We identified (at an FDR of 10%) 13,915 cis methylation QTLs (meQTLs)-i.e., CpG sites in which changes in DNA methylation are associated with genetic variation at proximal loci. We found that meQTLs are frequently associated with changes in methylation at multiple CpGs across regions of up to 3 kb. Interestingly, meQTLs are also frequently associated with variation in other properties of gene regulation, including histone modifications, DNase I accessibility, chromatin accessibility, and expression levels of nearby genes. These observations suggest that genetic variants may lead to coordinated molecular changes in all of these regulatory phenotypes. One plausible driver of coordinated changes in different regulatory mechanisms is variation in transcription factor (TF) binding. Indeed, we found that SNPs that change predicted TF binding affinities are significantly enriched for associations with DNA methylation at nearby CpGs.

Figure 1

A) QQ plot of –log10 p-values for testing the null of no association between methylation levels measured by all probes that passed our quality filters, and all SNPs within 3 kb of these probes. Data for SNPs within the candidate window are in black; negative control SNPs for which we chose a random 6 kb window elsewhere in the genome are in green; SNPs with the genotype labels permuted are in blue. B) Average methylation levels estimated using the bisulfite sequence data at meQTL probes, segregated by meQTL genotype. C) Histogram showing the distribution of distances between meQTL SNPs and the associated methylated sites in base pairs, for meQTLs where there is a single most likely causal site.

Figure 2

A) QQ plot of –log10 p-values for testing the null of no association between eQTL SNPs and methylation levels in sites within 3 kb. Positive correlations between expression and methylation levels are in red; Negative correlations are in blue, Data for random SNPs within the candidate window are in green; and data for a set of permuted genotype labels are in black. B) A plot of similar structure considering the associations of dsQTL SNPs and with methylation levels at sites within 3 kb. C) A plot of similar structure considering the QQ plots of associations between histone modification QTLs and methylation levels at sites within 3 kb.

Figure 3. Read counts segregated by meQTL genotype for multiple regulatory phenotypes.

The green line denotes the meQTL and the location of the probe measuring methylation data associated with the meQTL is identified by the black rectangle. The different colored data series indicate mean read depths segregated by genotype at the meQTL site: blue shows the homozygous genotype associated with low methylation level, orange shows the heterozygote, and purple the homozygous genotype associated with high methylation level. In this example, all of the regulatory phenotypes are negatively associated with DNA methylation levels.

Figure 4

A) Two-sided QQ-plots describing the effect of TF binding on DNA methylation. For each SNP in a predicted TF binding site we tested whether the SNP was associated with methylation at sites within 500 bp. Positive associations (upper right quadrant) indicate that the allele associated with increased PWM score for the TF in question is associated with increased methylation; negative associations (lower left quadrant) indicate that increased PWM score is associated with decreased methylation. We used a random set of SNPs in DNase I hypersensitive sites (DHSs) to indicate the expected baseline. When considering the control DHS SNPs, the direction of the effects was chosen randomly for the purpose of plotting. Panel B) additionally highlights four TFs that show particular strong association with changes in methylation levels. C) Two-sided QQ-plot of associations between Stat5 expression and DNA methylation at sites within 500 bp of Stat5 binding sites. D) QQ-plot of associations between ZNF274 expression and DNA methylation near ZNF274 binding sites. In both C and D, the grey shading indicates a region that would contain the data 95% of the time when the null hypothesis is true for all tests, obtained based on permutation of the expression data while holding the methylation data constant.