DNA methylation patterns associate with genetic and gene expression variation in HapMap cell lines

Abstract

Background: DNA methylation is an essential epigenetic mechanism involved in gene regulation and disease, but little is known about the mechanisms underlying inter-individual variation in methylation profiles. Here we measured methylation levels at 22,290 CpG dinucleotides in lymphoblastoid cell lines from 77 HapMap Yoruba individuals, for which genome-wide gene expression and genotype data were also available.

Results: Association analyses of methylation levels with more than three million common single nucleotide polymorphisms (SNPs) identified 180 CpG-sites in 173 genes that were associated with nearby SNPs (putatively in cis, usually within 5 kb) at a false discovery rate of 10%. The most intriguing trans signal was obtained for SNP rs10876043 in the disco-interacting protein 2 homolog B gene (DIP2B, previously postulated to play a role in DNA methylation), that had a genome-wide significant association with the first principal component of patterns of methylation; however, we found only modest signal of trans-acting associations overall. As expected, we found significant negative correlations between promoter methylation and gene expression levels measured by RNA-sequencing across genes. Finally, there was a significant overlap of SNPs that were associated with both methylation and gene expression levels.

Conclusions: Our results demonstrate a strong genetic component to inter-individual variation in DNA methylation profiles. Furthermore, there was an enrichment of SNPs that affect both methylation and gene expression, providing evidence for shared mechanisms in a fraction of genes.

Figure 1

Distribution of methylation patterns across the genome. (a) Methylation patterns for CpG-sites on autosomes, X-chromosome, and in the vicinity of imprinted genes. Methylation values are plotted for 77 individuals at 21,289 autosomal CpG-sites (left), for 43 females at 997 CpG-sites on the X-chromosome (middle), and for 77 individuals at 153 CpG-sites in 33 imprinted genes (right). (b) Methylation levels with respect to the TSS (negative distances are upstream from the TSS), where the line represents running median levels in sliding windows of 300 bp. (c) Correlations in methylation levels for all pair-wise CpG-sites (black), and for CpG-sites where both probes are in the same CGI (red), or where at least one probe is outside of CGIs (blue). Lines indicate smoothed spline fits of the mean rank pairwise correlation between CpG-sites in 100 bp windows, weighted by the number of probe pairs. (d) Methylation levels inside and outside of annotation categories, including CpG Islands (CGIs) for probes within 100 bp of the TSS, and histone modifications and transcription factor (TF) binding sites for all probes (see Additional file 1).

Figure 2

DNA methylation is negatively correlated with gene expression. (a) Methylation levels are low in the top quartile of highly expressed genes (left), and high in the bottom quartile of lowly expressed genes (right), looking across 12,670 autosomal genes. (b) Methylation levels with respect to the TSS in sets of genes categorized by gene expression levels, from highest (red) to lowest (blue), using the quartiles of gene expression with respect to gene expression means, where fitted lines represent running median levels (see Figure 1b).

Figure 3

Cis methylation QTLs. (a) Quantile-quantile (QQ) plot describing the enrichment of association signal in cis compared to the permuted data (90% confidence band shaded). (b) The cis-meQTL SNPs were enriched for association signal at additional CpG-sites near to the CpG-site for which they are meQTLs. The 180 best-associated SNPs were tested for association to probes that fell within 2 kb (red), within 2 kb to 10 kb (purple), and within 10 kb to 50 kb (blue) of the original best-associated CpG-site. The majority (96%) of probes within 2 kb (red) were in the same CGI as the best-associated probe. (c) Spatial distribution of cis-meQTLs with respect to the CpG-site as estimated by the hierarchical model.

Figure 4

The overlap between meQTLs and eQTLs. (a) QQ-plot describing the eQTL association P-values in 180 cis-meQTL SNPs (red) and in eight samples of SNPs that match the cis-meQTL SNPs for minor allele frequency and distance-to-probe distributions (black). (b) Association signals in 508 FDR 10% eQTLs before and after regressing out gene-specific methylation. In black are 439 eQTLs that overlap across the two phenotypes, in red are 45 eQTLs present before methylation regressions, and in blue are 24 eQTLs present after regressing out methylation. The flat lines (green) correspond to the FDR 10% eQTL threshold.

Figure 5

C21orf56 gene region. (a), (b), (c) Genotype at rs8133082 is associated with methylation (cg07747299) and gene expression at C21orf56, plotted per individual colored according to genotype at rs8133082 (GG = black, GT = green, TT = red) for directly genotyped (circles) and imputed (triangles) data. (d) Gene expression levels at C21orf56 after regressing out methylation. (e) Gene expression at C21orf56 (+/-2 kb) genomic region on chromosome 21. Distance is measured on the reverse strand relative to C21orf56 TSS at 46,428,697 bp. Barplots show average gene expression reads per million in the subsets of individuals from each of the three rs8133082-genotype classes. Middle panel shows histone-modification peaks in the region from Encode LCL GM12878. Bottom panel shows the gene-structure of C21orf56, where exons are in bold and the gene is expressed from the reverse strand. Green points indicate the location of four HapMap SNPs (rs8133205, rs6518275, rs8133082, and rs8134519) associated at FDR of 10% with both methylation and gene expression, and Figure S11 in Additional file 1 shows association results for this region with SNPs from the 1,000 Genomes Project. (f) Bisulphite-sequencing results for eight rs8133082-homozygote individuals (4 GG black, 4 TT red) validates the genome-wide methylation assay at cg07747299 and shows the extent of methylation in the surrounding 411 bp region.