Rare genetic variation at transcription factor binding sites modulates local DNA methylation profiles

Abstract

Although DNA methylation is the best characterized epigenetic mark, the mechanism by which it is targeted to specific regions in the genome remains unclear. Recent studies have revealed that local DNA methylation profiles might be dictated by cis-regulatory DNA sequences that mainly operate via DNA-binding factors. Consistent with this finding, we have recently shown that disruption of CTCF-binding sites by rare single nucleotide variants (SNVs) can underlie cis-linked DNA methylation changes in patients with congenital anomalies. These data raise the hypothesis that rare genetic variation at transcription factor binding sites (TFBSs) might contribute to local DNA methylation patterning. In this work, by combining blood genome-wide DNA methylation profiles, whole genome sequencing-derived SNVs from 247 unrelated individuals along with 133 predicted TFBS motifs derived from ENCODE ChIP-Seq data, we observed an association between the disruption of binding sites for multiple TFs by rare SNVs and extreme DNA methylation values at both local and, to a lesser extent, distant CpGs. While the majority of these changes affected only single CpGs, 24% were associated with multiple outlier CpGs within ±1kb of the disrupted TFBS. Interestingly, disruption of functionally constrained sites within TF motifs lead to larger DNA methylation changes at nearby CpG sites. Altogether, these findings suggest that rare SNVs at TFBS negatively influence TF-DNA binding, which can lead to an altered local DNA methylation profile. Furthermore, subsequent integration of DNA methylation and RNA-Seq profiles from cardiac tissues enabled us to observe an association between rare SNV-directed DNA methylation and outlier expression of nearby genes. In conclusion, our findings not only provide insights into the effect of rare genetic variation at TFBS on shaping local DNA methylation and its consequences on genome regulation, but also provide a rationale to incorporate DNA methylation data to interpret the functional role of rare variants.

Fig 1. Effect of regulatory rare variants on shaping local DNA methylation profiles.

(A). Analysis pipeline of the study. Briefly, after identifying rare variants that overlap TFBSs, we extracted DNA methylation values corresponding to CpGs that fell within the targeted TFBSs and their extended flanks (±1 kb) present on the EPIC array. For each CpG site, DNA methylation values were transformed into ranks across all individuals (n = 247) in ascending order, ranging from 1 to 247. Thus, DNA methylation values rise through the ranks, with the minimum and the maximum DNA methylation values in our cohort represented by the lowest and highest ranks, respectively. CpGs showing both DNA methylation ranks present into the 5% tails of the distribution and DNA methylation difference of ≥0.05 (β-value) as compared to the mean of controls, i.e, non-carriers for the tested rare variant, were considered as differentially methylated (also referred as DNA methylation outlier). (B) Plot showing DNA methylation profile at chr2:158,183,226–158,185,240 (CTCF binding site ±1 kb). Red dots indicate DNA methylation values for individual carrying rare variant (chr2:158,184,228), while black dots and bars represents the mean and mean ±2 standard deviation of controls, respectively. Disrupted CTCF motif (chr2:158,184,226–158,184,240) is highlighted in light green within the CTCF ChIP-peak (black box). Position of rare variant is depicted by vertical gray dashed line. (C) Plot showing distribution of DNA methylation ranks across population for all tested TFBSs (n = 133). Red dots represent the frequency of the rank for individuals carrying the rare variant, while the black horizontal line is the frequency in controls. Peaks at gray shaded areas of the graph represent an excess of individuals with extreme DNA methylation values for a given CpG site.

Fig 2. TFBS associated with local DNA methylation profile.

Bars indicate P-values obtained for 46 TFBSs that were significantly enriched for outlier methylation marks when mutated by permutation analysis (see Methods), after Bonferroni-correction for the number of tested TFBSs (n = 133). For each factor, enrichment for differentially methylated CpGs between regulatory rare variants carriers and non-carriers is shown at the base of the bar in white. Significance threshold is indicated by the red-dashed vertical line.

Fig 3. Mutation of highly constrained positions in TFBS motifs are associated with larger changes in DNA methylation.

(A and B) Sequence logos showing SNVs within the consensus CTCF binding motif. While substitutions at degenerate positions (A) result in small changes to the position weight matrix score (ΔPWM), substitutions at highly conserved positions (B) cause large changes in ΔPWM. (C) Violin plots showing the distribution of ΔPWM scores for SNVs in TFBSs that are associated with different degrees of change in local DNA methylation. SNVs were stratified into eight different bins according to the degree of change in DNA methylation of associated CpGs (black and gray-filled plots corresponding to DNA methylation outlier and non-DNA methylation outlier, respectively). Red dots represent the median ΔPWM score for each bin, while the red line represents the smoothed median of three consecutive points. Above each violin is shown the number of associated CpGs per bin. (D) Violin plots showing the distribution of ΔPWM scores for SNVs that are associated versus those not associated with outlier methylation. The p-value is derived from the Wilcoxon matched-pairs signed rank test. Interquartile range (IQR) and median of the distribution are represented by boxes and white dots in the overlaid box plots, respectively. Whiskers represent upper/lower quartiles +/-1.5 IQR.

Fig 4. Mutations within TFBS are associated with larger effects on local DNA methylation.

Violin plots showing enrichments of DNA methylation outliers associated with SNVs occurring within or in the flanks of TFBSs. Overlaid box plots as described in Fig 3. P-value derived from paired t-test.

Fig 5. Long-range effects of regulatory rare variants on DNA methylation.

For each TFBS containing a rare variant, CpGs located up to 100 kb away were grouped into bins according to their distance. Dots indicate p-values obtained for each TFBS by permutation testing after Bonferroni-correction for the number of TFBSs tested (n = 133). TFBSs with Bonferroni corrected p<0.05 in each bin are colored red, all others in black.

Fig 6. DNA methylation outliers linked with disruption of TFBS are enriched for gene expression outliers.

(A) Bar plot showing the fraction of rare variants disrupting TFBSs that overlap promoters (TSS±2kb) of genes with outlier (black) and non-outlier expression (gray) in cardiac tissue. P-value derived from two-tailed Fisher’s exact test. (B) Plot showing DNA methylation profile at chr22:24,058,597–24,060,611 (UA4 binding site ±1 kb) genomic region. DNA methylation values are shown for an individual with a rare variant at chr22:24,059,610 (red dots), while black dots and bars represents the mean ±2 standard deviation of controls. Position of rare variant that disrupts the UA4 TFBS motif is depicted by vertical gray dashed line. Expression of the nearby GUSBP11 gene (C) is down-regulated in this individual carrying the rare variant (red dot) compared to controls. Expression of GUSBP11 in controls is represented by box plot, with gray shaded areas indicating outlier expression levels.