Effect of estrogen receptor α binding on functional DNA methylation in breast cancer

Abstract

Epigenetic modifications introduce an additional layer of regulation that drastically expands the instructional capability of the human genome. The regulatory consequences of DNA methylation is context dependent; it can induce, enhance, and suppress gene expression, or have no effect on gene regulation. Therefore, it is essential to account for the genomic location of its occurrence and the protein factors it associates with to improve our understanding of its function and effects. Here, we use ENCODE ChIP-seq and DNase I hypersensitivity data, along with large-scale breast cancer genomic data from The Cancer Genome Atlas (TCGA) to computationally dissect the intricacies of DNA methylation in regulation of cancer transcriptomes. In particular, we identified a relationship between estrogen receptor α (ERα) activity and DNA methylation patterning in breast cancer. We found compelling evidence that methylation status of DNA sequences at ERα binding sites is tightly coupled with ERα activity. Furthermore, we predicted several transcription factors including FOXA1, GATA1, and SUZ12 to be associated with breast cancer by examining the methylation status of their binding sites in breast cancer. Lastly, we determine that methylated CpGs highly correlated with gene expression are enriched in regions 1kb or more downstream of TSSs, suggesting more significant regulatory roles for CpGs distal to gene TSSs. Our study provides novel insights into the role of ERα in breast cancers.

Keywords: :

• ChIP-seq
• DNase I hypersensitivity
• breast cancer
• differential gene expression
• differential methylation
• estrogen receptor α
• transcription factor

Introduction

DNA methylation is a biochemical process that modifies the cytosine nucleotides in the context of CpG dinucleotides (CpGs) by the addition of a methyl group to the fifth carbon position. DNA methylation plays critical roles in many important biological processes including genomic imprinting,¹ X-chromosome inactivation,² transposable elements silencing,³ stem cell differentiation,⁴ embryonic development,⁵ and inflammation.⁶ In humans, DNA methylation patterns are precisely regulated to maintain a delicate balance between stability and plasticity. Alterations in DNA methylation have been demonstrated to interact with genetic events and to be involved in human carcinogenesis for nearly all cancer types.^7,8 A global shift in DNA hypomethylation in cancer cells has been reported, which is implicated in the development and progression of cancer.^9,10 More specifically, genome wide DNA methylation profiling has been performed in breast cancer to identify genes associated with tumorigenesis.^11-16 DNA methylation signatures or markers have been defined to classify breast cancer subtypes^17-19 and to predict prognostic outcomes, e.g., patient survival.^10,20

DNA methylation may affect gene expression by directly impacting the binding of transcription factors (TFs).²¹ It has been suspected that DNA methylation physically impedes the binding of transcription factors to their binding sites.^22-24 While this might be the case for most transcription factors, exceptions have been encountered in several studies. For example, Holler et al. showed that Sp1 is capable of binding DNA and activating transcription even when the binding site is methylated.²⁵ In addition, Guillaume et al. showed that a family of zinc finger proteins can bind methylated DNA and repress gene transcription.²⁶

Alternatively, DNA methylation may also regulate transcription by modifying local chromatin structure; however, the exact mechanisms by which this occurs are unclear.²⁷ There is convincing evidence demonstrating a linkage between DNA methylation and chromatin structure mediated by methylcytosine-binding proteins (MBPs).²⁸ A subset of these MBPs contain conserved methylcytosine-binding domains (MBD) that recognize and bind methylated cytosines and recruit additional chromatin remodeling factors such as histone deacetylases and histone methylases, leading to compacted inactive local chromatin structure. Another set of proteins that have been shown to function as methylcytosine-binding proteins contain SET- and Ring-associated (SRA) domains such as UHRF1.^29,30 Furthermore, Kaiso-like zinc finger motifs have been shown to bind single methylated CpGs.³¹ The large variety of protein motifs capable of binding methylated DNA is indicative of the complex interplay involving protein factors that couple DNA methylation to chromatin structure.

Although both of the above two mechanisms imply a repressive effect of DNA methylation on gene transcription, studies have shown a more complicated relationship between DNA methylation and gene expression.³² Generally, methylation in the immediate vicinity of the TSS blocks initiation, but methylation in the gene body does not block and might even stimulate transcription elongation.³² Thurman et al. examined the correlation between methylation levels at transcription factor binding sites (TFBS) and transcription factor abundance within DNase I hypersensitive (DHS) sites. They observed that 70% of transcription factors were negatively correlated with DNA methylation, whereas only a few transcription factors exhibited significant positive correlations. In general, CpG methylation within transcription factor binding sites is negatively correlated with the expression level of the corresponding transcription factor. Furthermore, they argue that a negative correlation between CpG methylation and transcription factor gene expression indicates that DNA methylation is a passive process, i.e., methylation fills in the voids left by vacating transcription factors.³³

In breast cancer, estrogen receptor (ESR1) activity status is a critical biomarker for subtype classification and is widely used to determine whether or not a patient should receive hormone therapy such as Tamoxifen treatment.³⁴ Consequently, DNA methylation predictors have been proposed as a clinical marker for ESR1 activity.³⁴ Differential methylation of ESR1 in breast carcinomas was first described by Piva et al. using the methylation–sensitive endonuclease HpaII.³⁵ More recently, several studies have applied high-throughput technologies including deep sequencing and microarrays to study DNA methylation at a genome-wide level. Li et al. identified 5 genes that were significantly differentially methylated between 12 ER+ and 12 ER- breast tumors using the Infinium Methylation Assay.³⁶ Similarly, Fackler et al. interrogated 27 578 CpG loci to deduce which genes were most associated with ER status in 103 breast tumor samples.¹⁴ Another approach used MethyLight to measure the methylation levels of 35 gene markers to classify 148 primary breast carcinomas.³⁷ All of these studies performed 2-dimensional unsupervised hierarchical cluster analysis to identify the most differentially methylated CpGs between ER+ and ER- breast cancer.^14,36,37 To understand the spatial distribution of aberrant CpG methylation, Ruike et al. used methyl-DNA immunoprecipitation followed by high-throughput sequencing to identify genomic regions in breast cancer cell lines that exhibit hyper- or hypomethylation.³⁸ However, these studies did not investigate the association of DNA methylation with breast cancer by interrogating >485 000 CpG sites at the level of specific transcription factor binding using methylation and gene expression data from 222 TCGA-derived breast tumor samples.¹¹

In this report, we conducted a detailed study of the relationship between the binding of a sequence-specific transcription factor and the methylation level at its corresponding binding sites, using the well-studied estrogen receptor in breast cancer as a model system. Our study revealed that methylation level of ESR1 binding sites is negatively correlated with ESR1 expression levels, and ESR1 binding sites tend to be methylated in ER- breast cancers. In addition, our results indicate that ESR1 exerts its effect on DNA methylation within its binding regions in a localized fashion. Based on this conjecture, we further predicted FOXA1 and GATA3 to be overactive in ER+ breast cancers. In addition, we determined CTBP2 and PRC2 family member SUZ12 to be positively associated with DNA methylation. Finally, we found that CpGs in DNase I hypersensitive regions are more likely to be negatively correlated with expression of corresponding genes, which is consistent with the findings that most transcription factors are trans-activating. This analysis bridges a comprehensive and high-resolution portrait of the breast cancer DNA methylome to the regulatory processes responsible for breast cancer classification. Specifically, by integrating ENCODE and TCGA data sets, we link DNA methylation to transcription factor binding to chromatin state; all of which are integral in determining a final gene expression output.

Results

Correlation of DNA methylation with ESR1 expression

An overview of our analysis strategy is provided in Figure 1. We focused on determining whether or not genomic features of CpG sites (those which are bound by ER α or other TFs, or located in DNase I hypersensitive sites) impact their methylation levels and their correlation with gene expression. To achieve high-resolution, this analysis was conducted by considering CpGs specifically located in TF binding sites and in DHS regions. In our first analysis, we operate under the assumption that ESR1 (gene that encodes ERα) expression is a proxy for ERα activity and correlated the DNA methylation level of all CpGs with ESR1 expression levels across all TCGA breast cancer samples stratified on ER status (see Fig. 2A for an example). On average, the Spearman correlation coefficient (SCC) between overall CpG methylation (across the whole genome) and ESR1 expression is -0.056. As part of the ENCODE blueprint, ChIP-seq data was generated for >100 TFs in various cell lines incubated under different treatments.³⁹ Using TCGA and ENCODE data sets, we defined a CpG set consisting of CpGs located in genomic regions not bound by ERα and determined that the average correlation between non-ERα binding CpGs and ESR1 expression was –0.083. Conversely, we correlated methylation of CpGs in genomic regions bound by ERα with ESR1 expression and obtained a striking correlation coefficient as extreme as –0.20 (Fig. 2B). Among CpGs not in ERα binding regions, 5.8% yield r > 0.4 and 2.8% yield r < –0.4 in their correlation with ESR1 expression (Fig. 2C). Contrastingly, in the case of CpGs in ERα binding regions, 0.55% yield r > 0.4 and 24% yield r < –0.4 (Fig. 2C, all SCCs are based on 222 samples; a correlation coefficient of r > 0.4 or r < –0.4 corresponds to a p-value of P < 3e-10). This suggests a strong enrichment of methylated CpGs in ERα binding sites that are negatively correlated with ESR1 expression. Moreover, for each CpG in an ERα binding region we calculated and compared its average DNA methylation level between ER+ and ER- breast cancer samples. As expected, the majority of CpGs in ERα binding regions exhibit lower average methylation levels in ER+ than in ER- samples (Fig. 2D).

Figure 1. The schematic diagram of our analysis. We combined DNA methylation and gene expression data in breast cancer samples from TCGA, and TF binding and DNase I hypersensitivity data from ENCODE. We identified the CpG sites with differential methylation levels between ER+ and ER- breast cancer samples, and examined the correlation of their methylation levels with expression of associated genes. Blue double arrows denote comparative analysis of regions of interest to outside regions.

Figure 2. Correlation between CpG methylation and ESR1 expression levels. (A) An example- correlation between methylation level of cg03387103 (corresponding to LETM1) and ESR1 expression level in all breast cancer samples. (B) CpGs in ER binding peaks have larger negative correlations with ESR1 expression in their methylation levels. (C) Fraction of CpGs highly correlated with ESR1 expression. There is a higher fraction of anti-correlated CpGs in ERα binding sites compared with non-ERα binding sites at ± 0.4 SCC cutoff. (D) CpGs in ER binding peaks have higher average methylation levels in ER- than ER+ samples. Each point is a CpG. SCC: Spearman correlation coefficient.

Distribution of differential DNA methylation between ER+ and ER

The DNA methylation levels of many CpGs in breast cancer samples are dependent on ER status. Some CpGs demonstrate higher methylation levels in ER+ than in ER- samples (Fig. 3A), while others show the opposite trend (Fig. 3A). We systematically investigated the distribution of CpGs with significant differential methylation levels between ER+ and ER- samples. Specifically, we calculated the position of each significant CpG relative to the transcription start site (TSS) of the gene it is associated with (Fig. 3B). As shown in Figure 3C, the distribution of significant CpGs is centered at the TSS of genes, suggesting that there is a greater number proximal to gene TSS. However, this provides no indication as to how probable any CpG selected at random from all CpGs will be significantly differentially methylated because there is an inherent enrichment of significant and non-significant CpGs vicinal to gene TSS. Therefore, by calculating the fraction of significant CpGs to the total number of CpGs at each genomic coordinate we account for the non-uniformity of CpGs distributed across genes. Consequently, we observe that CpGs nearby TSS have a lower likelihood (relative frequency) of exhibiting significant differential methylation than those in distal DNA regions (Fig. 3D). Overall, this result suggests that CpGs at locations distal from TSS might be equally or even more functionally relevant. Moreover, we observed that, at the same significance level, the fraction of hypermethylated CpGs (Fig. 3D, red) is higher than that of hypomethylated CpGs in ER+ samples (hypermethylated CpGs in ER- samples.) (Fig. 3D, green).

Figure 3. Distribution of CpGs with differential methylation levels between ER+ and ER- breast samples. (A) Examples- CpG may have higher methylation in ER+ (cg05846044) or ER- (cg05859267). (B) Relationship between differential methylation and CpG position relative to transcription start site (from –1500 upstream to 4500 downstream of TSS). (C) Distribution of CpGs with significant (P < 1e-6) differential methylation between ER+ and ER- samples. (D) Fraction of CpGs with significant differential methylation levels at different positions. The fraction is the ratio of the number of significant CpGs to the total number of CpGs in a DNA window. CpGs with significantly higher methylation levels in ER+ (red) and in ER- (green) samples are examined separately. The relative frequency of significantly differentially methylated CpGs increases as the distance from TSS increases.

Impact of ERα binding on DNA methylation

We next investigated the relationship between DNA methylation and transcription factor (TF) binding based on ENCODE ChIP-seq data. First, we considered the question: Are CpGs in ERα binding regions more likely to be differentially methylated between ER+ and ER- breast cancer samples? Based on the ERα binding peaks in T47d cell line, we defined an ERα-binding region CpG set and a non-ERα-binding CpG set as a control. The former consists of CpGs that fall precisely within an ERα binding peak. The latter contains CpGs that do not fall directly within ERα binding peaks, but do fall in gene regions that contain these binding peaks. Overall, CpGs included in the analysis are selected from regions within the genes that are bound by ERα (i.e., has a binding peak in their gene body or promoter).

We find that CpGs that are localized within ERα binding peaks exhibit lower DNA methylation in ER+ than in ER- samples (Fig. 4A), suggesting a negative correlation between ERα binding and site-specific DNA methylation. For example, if ER+ samples are analyzed at a significance level α = 1e-6, 31% of CpGs in ERα binding regions exhibit lower methylation levels whereas only 1.1% exhibit higher methylation levels (Table S1). In contrast, if CpGs from all genomic locations are considered, only 11.3% of CpGs have higher methylation levels and 6.7% of CpGs have lower methylation levels in ER+ compared with ER- samples. Likewise, similar fractions were observed with non-ERα binding CpGs. This trend remains stable when different significance thresholds are used (Fig. 4A).

Figure 4. CpGs in ER binding sites tend to have lower methylation levels in ER+ breast samples. (A) The fraction of CpGs with significant differential methylation levels between ER+ and ER- samples. Note that CpGs in ER binding regions tend to have higher methylation levels in ER- samples, while CpGs not in ER binding regions tend to have higher methylation levels in ER+ samples. Four different thresholds are used to determine differential methylated CpGs. (B) Distribution of t-scores (ER+ vs. ER-) of methylation levels for CpGs. Genes are divided into 3 classes based on their expression levels in ER+ vs. ER- samples: ER+ > ER- (red), ER+ < ER- (green), and ER+ = ER- (white). Distributions of CpGs associated with the three gene classes are shown separately. (C) CpGs in ER binding regions tend to have lower methylation levels in ER+ samples (lower t-scores) compared with thosenot in ER binding regions, which is the case for CpGs associated with all three gene classes.

Because promoter DNA methylation is generally negatively correlated with gene expression status, we compared expression levels of genes between ER+ and ER- samples. First, we defined 3 gene categories: upregulated, downregulated, and non-differentially expressed genes in ER+ vs. ER- samples. To quantify the difference in methylation levels between ER+ and ER- samples, we calculated the t-scores of β values for each gene category (ER+ vs. ER-). As shown in Figure 4B, CpGs in ER+ upregulated genes tend to have lower t-scores (i.e., ER+ is hypomethylated) as compared with CpGs in ER+ downregulated genes. In spite of this trend, the ERα binding CpGs demonstrated significantly lower methylation t-scores than non-ERα binding CpGs in all three of the gene categories (Fig. 4C).

To further investigate the impact of ERα binding on DNA methylation, we calculated the CpG methylation levels as a function of its distance to the center of an ERα binding peak (Fig. 5A). Strikingly, we find that CpGs closer to the center are more likely to have larger negative t-scores, namely, are more likely to have lower DNA methylation levels in ER+ than in ER- samples (Fig. 5B). Consistent with the function of SUZ12, CpGs closer to the center of SUZ12 binding peaks are more likely to have larger positive t-scores (Fig. 5C).

Figure 5. Relationship between differential methylation of CpGs and TF binding. (A) Binding of some TFs is correlated with reduced methylation level of CpGs in ER+ relative to ER- samples, while binding of others (SUZ12 and CTBP2) is correlated with increased methylation level. (B) CpGs proximal to ER binding center are more likely to have lower methylation levels in ER+ (smaller t-scores for ER+ vs. ER- comparison). (C) CpGs proximal to SUZ12 binding center are more likely to have higher methylation levels in ER+ (larger t-scores for ER+ vs. ER- comparison).

Taken together, our results indicate that the impact of ERα binding on DNA methylation is restricted to a local genomic region. The methylation level of ERα binding region CpGs is determined mainly by the ER status of samples (ER+ or ER-) rather than by the transcriptional status of genes (upregulated or downregulated in ER+).

Impact of other TF binding on DNA methylation

We next extended the analysis to other TF binding sites by defining a TFBS CpG set and a non-TFBS CpG set for all TF binding data from ENCODE, and compared methylation levels of TFBS CpG sets between ER+ and ER- samples. Although ERα seems to be the TF with the most significant impact on differential DNA methylation between ER+ and ER- samples, there are some other TFs that also exhibit influence (Fig. 5A). For instance, the CpGs located in FOXA1 and GATA3 TF binding sites are significantly hypomethylated in ER+ than in ER- as compared with the corresponding non-TFBS control CpGs. Consistently, these two TFs have been reported to function upstream of ER and mediate ER binding in breast cancer.^40,41

From the ENCODE data, four data sets containing ERα binding peaks were generated, which includes ERα binding analysis performed under treatment with two different steroid hormones (Gen1h and Estradia1h) in two cell lines responsive to primary steroid hormone treatment (T47d and Ecc1). Interestingly, the methylation difference (t-scores in ER+ vs. ER-) between ERα binding CpGs and non-ERα binding CpGs is much more obvious for peaks identified in T47d than those identified in Ecc1 (Fig. 5A). Given the fact that T47d is a breast epithelial-derived and Ecc1 is an endometrium epithelial-derived cell line, this result, as expected, likely indicates that T47d better reflects the ERα binding events in human breast cancer tissue than Ecc1.

Moreover, we also identified a number of TFs whereby the TFBS CpGs had larger t-scores than the non-TFBS CpGs, implying that binding of these TFs would enhance DNA methylation. One of these TFs is SUZ12, a component of the polycomb repressive complex 2 (PRC2), which catalyzes methylation of H3K7.⁴² Considering the chromatin-silencing role of PRC2,⁴³ it may not be surprising to observe the enhanced DNA methylation in the SUZ12 binding CpGs. Another example is CtBP2, which also show higher DNA methylation in its binding sites. Interestingly, CTBP2 has been reported to function as a transcriptional repressor.⁴⁴

Correlation between DNA methylation and gene expression

Depending on the genomic position and other factors, methylation of a CpG can be positively (Fig. 6A) or negatively (Fig. 6B) correlated with the expression levels of its associated genes. Therefore, for each CpG with a unique gene assignment, we calculated the Spearman correlation coefficient between its methylation level and the gene expression levels. Figure 6C shows the relationship between correlation and relative position (the distance from CpG to the TSS of its associated gene) of CpGs. As shown, there are significantly more instances of negative correlation than positive correlation and a large number of negative correlations occur at the DNA region proximal to the TSS. More clearly, the distributions of CpGs with r > 0.4 (the red line) or r < –0.4 (the green line) are shown in Figure 6D. As shown, most of the CpGs negatively correlated with expression are located in a DNA region upstream of TSS, whereas the CpGs with positive correlations exhibit two peaks, one in the gene body and the other in the promoter region. However, after taking into account the biased distribution of CpGs interrogated by the Illumina 450k DNA methylation array (the black line in Fig. 6D), the fraction of high correlation CpGs (r > 0.4 or r < –0.4) is maximal in the DNA region more than 1kb downstream of TSS (Fig. 6E) rather than in the TSS region.

Figure 6. Correlation of CpG methylation level with expression level of the associated genes. (A) Methylation level of cg06228260 is positively correlated with its associated gene PTPRN2. (B) Methylation level of cg01586506 is negatively correlated with its associated gene SOX10. (C) Relationship between methylation-expression correlation and CpG position relative to transcription start site (from –1500 upstream to 4500 downstream of TSS). (D) Distribution of CpGs with strong correlations in methylation with expression level of the associated genes. Positive correlation (red, r > 0.4) and negative correlation (green, r < –0.4) are examined separately. (E) Fraction of CpGs strongly correlated with expression of the associated genes at different positions. (F) CpGs in ER binding regions tend to have negative correlation in their methylation with expression of their associated genes.

In parallel to what we have done in differential DNA methylation analysis, we compared the correlations with gene expression levels between ERα binding region CpGs and non-ERα binding region CpGs. The results indicate that CpGs in ERα binding peaks are more likely to be negatively correlated with gene expression levels (Fig. 6F). Overall, 1.4% and 3.3% of all CpGs have positive (>0.4) and negative (<–0.4) correlations, respectively. Similarly, the fraction of high correlations is just slightly higher than the overall fraction in the non-ERα binding CpGs. However, 0.4% of ERα binding CpGs are positively correlated with gene expression and 7.7% are negatively correlated when the same significance threshold is maintained. This implies that CpGs specifically in ERα binding regions are more likely to exert their influence on gene expression regulation in breast cancers.

CpGs in DNase hypersensitive sites

TF binding data from ChIP-seq capture the binding events of a single TF in each experiment. The DNase I hypersensitivity data, however, identify the DNA regions enriched for all DNA regulatory elements.³³ Based on the DNase hypersensitivity data in T47d, we defined a DHS CpG set and a non-DHS CpG set. First we compared the fraction of differentially methylated CpGs (ER+ vs. ER-) in these two sets. Interestingly, we found that differentially methylated CpGs are depleted in DHS (Fig. 7A). This observation is consistent with the fact that DHS is enriched for all types of regulatory elements, among which only a small fraction (e.g., ER binding sites, FOXA1 binding sites) are correlated with ER status. Most regulatory elements should show similar activities between ER+ and ER- samples. Consequently, these elements should have similar DNA methylation states in ER+ and ER- samples.

Figure 7. Comparison of CpGs in and not in DNase hypersensitive sites. (A) CpGs in DHS and non-DHS exhibit no significant difference in differential methylation between ER+ and ER- breast cancer samples. (B) CpGs in DHS tend to have negative correlation in their methylation level with the expression of their associated genes.

We also compared DHS CpGs and non-DHS CpGs in the correlation of their methylation levels with gene expression. As shown in Figure 7B, CpGs with high negative correlations are enriched in DHS sites, whereas CpGs with high positive correlations are depleted in DHS sites. This suggests that CpGs involved in gene expression regulation is enriched in DHS regions. Most of such regulation might be mediated by the binding of positive regulators, e.g., transcription activators, which leads to reduced DNA methylation and thus a negative correlation with gene expression.

Discussion

In this study, we investigated the relationship of TF binding regions with DNA methylation using ERα binding activity as a model. We found that CpGs in ERα binding peaks were more likely to be hypomethylated in ER+ than in ER- breast cancer samples. Furthermore, methylation of these CpGs had a greater likelihood of being negatively correlated with gene expression. These results indicate that CpG methylation in distinct ERα binding sites may be dependent on ERα activity and that physical binding of ERα to its cognate DNA sequence has the potential to inhibit methylation of these CpGs. Moreover, we showed that such an effect was restricted to a local DNA region proximal to the center of ERα binding peaks (Fig. 5A). Lastly, by increasing the resolution of our analysis by considering CpGs in DHS regions, we observed that these regions harbor a large fraction of CpGs negatively correlated with gene expression. This result suggests that these CpGs have a higher probability of being functionally relevant since DHS regions are more accessible to protein regulators. Overall, this suggests a model whereby TF binding events impact the methylation status of local CpGs, and the final effect of DNA methylation on gene expression is determined by the overall output of each neighboring CpG’s methylation status. Additionally, the methylation status of CpGs might be determined by the binding of many different TFs cooperatively or competitively. Instead of acting as the readout of gene expression, DNA methylation may participate in transcriptional regulation of genes in a more active and delicate manner than has been expected; different CpGs independently read in binding signals of different TFs and integrate them.

Limited by the data source, we focused on the CpGs that were included in HumanMethylation450k array, which contained probes that mainly targeted CpGs in the transcribed region of genes or that were nearby gene TSSs. This most likely reflects an inherent genetic bias that may or may not be intensified by the array platform. After correcting for this bias, our results indicate that genomic coordinates localized more than 1kb downstream of gene TSS tend to have a higher fraction of significantly methylated CpGs. This challenges the notion that methylated CpGs proximal to gene TSSs are the major players in gene expression regulation. A recent paper by Aran et al. explored the relationship between the DNA methylation of distal regulatory sites and the dysfunctional regulation of cancer genes. They showed that hypomethylated enhancer sites correlated with upregulation of cancer-related genes and hypermethylated sites with downregulation. Moreover, the association between enhancer methylation and gene deregulation in cancer was significantly stronger than the association of promoter methylation with gene deregulation.⁴⁵ It would be interesting to investigate the effect on TF binding on methylation of CpGs located in enhancers.

Thurman et al. observed in ENCODE cell lines that the methylation levels of TF binding sites were correlated with the expression levels of the corresponding TFs, and proposed that DNA methylation might be a passive reflection of transcription factor binding, i.e., filling in the voids left by vacating transcription factors. Here we validated their observations in tumor samples from patients with breast cancer. We found that the methylation levels of ERα binding CpGs tended to be lower in ER+ than in ER- samples. Compared with the ER- samples, ER+ samples have significantly higher ERα activity. Additionally, we also confirmed that binding of some TFs (e.g., ERα and FOXA1) were associated with reduced methylation levels, while binding of other TFs (such as SUZ12 and CtBP2) were associated with enhancer methylation levels.

Overall, our study integrates multiple large-scale data sets from TCGA and ENCODE to construe the association of DNA methylation patterning with the underlying transcriptional machinery within specific regions of the genome, in particular TF-binding sites and DHS regions. We expand on prior studies by providing a high-resolution analysis that illustrates the potential mechanistic relationship between CpG methylation and TF binding and describe how it affects differential gene expression observed between ER+ and ER- breast carcinomas. More specifically, we are able to assess the CpG methylation patterning at specific binding sites and show how it can influence cancer phenotype via its interaction with transcription factors. To our knowledge, this is the most comprehensive analysis of DNA methylation in breast cancer since we have used data from the latest HM450K technology coupled with gene expression data from 222 TCGA primary breast carcinoma samples, along with ENCODE ChIP-seq TF profiles. By understanding how DNA methylation patterning affects the activity of specific transcription factors, we can better determine molecular characteristics of patients with tumors. For example, if we can dissect the “methylation code” of a transcription factor, we can use the information to understand the transcriptional aberrations implicated in tumor types. This may aid in the development of biomarkers and/or targeted therapy.

Materials and Methods

Data sets

The gene expression and DNA methylation data for breast cancer patients were downloaded from the TCGA (The Cancer Genome Atlas) project at https://tcga-data.nci.nih.gov/tcga/. Expression levels of genes were quantified using the two-channel Agilent microarrays. Methylation levels of CpG were measured using the HumanMethylation450 arrays, and represented as β values. The β value is a quantitative measure of DNA methylation levels of specific CpGs, which ranges from 0 for completely unmethylated to 1 for completely methylated.

The genome wide TF binding data were generated by ENCODE (The Encyclopedia of DNA Elements) project based on ChIP-seq experiments. We downloaded the binding peaks of TFs from UCSC Genome browser at http://genome.ucsc.edu/ENCODE/downloads.html. We used the binding peaks identified by the peak calling algorithm, PeakSeq.¹³ The data set contains TF binding data in a number of different cell lines, from which we only select the data in breast epithelial cell lines (MCF7 and T47D) for our analysis to achieve the best match with data from TCGA.

The DNase hypersensitivity data were generated by ENCODE project based on DNase-seq experiments and were downloaded from UCSC Genome browser. The data provide a complete list of DNA regions that are sensitive to DNase I treatment, also known as DNase hypersensitive sites. Again, to match TCGA data, we only selected DNase data obtained from MCF7 and T47d cell lines.

Differential DNA methylation between ER+ and ER- breast cancer samples

The DNA methylation data from TCGA contains methylation levels of 485 577 CpGs in 630 ER+ and 187 ER- breast cancer samples. Most of the CpGs can be associated with a gene based on their localization: in the transcribed region or proximal to the transcription start site of a gene. For each CpG, we compared its β values in ER+ with respect to ER- samples by using the Student t test. Given a significance cut-off (e.g., P < 0.001), we identified a hypermethylated CpG set and a hypo-methylated CpG set in ER+ samples with respect to ER- samples. We examined a number of different cut-off values for significance.

Differential gene expression between ER+ and ER- breast cancer samples

The gene expression data from TCGA contains expression levels of 17 814 human genes in 401 ER+ and 118 ER- breast cancer samples. We compared the expression levels of genes between ER+ and ER- samples using the Student i test. By setting a cut-off value of P < 0.001, we divided genes into three classes: upregulated in ER+, downregulated in ER+, and non-differentially expressed genes.

Relating CpGs with ER binding, other TF binding, and DNase I hypersensitive sites

Given the complete list of ERα binding peaks in a cell line (e.g., T47d), we can determine whether a CpG is located within an ER binding peak. We defined ER binding CpGs as those falling into an ERα binding peak. In general, a gene is associated with multiple CpGs in HumanMethylation450k array. In this study, we aim to investigate the local effect of ERα binding on DNA methylation, thus we defined non-ERα binding CpGs as those not in any ERα binding peak, but were associated a gene with at least one ERα binding CpG. Since both ERα binding CpGs and non-ERα binding CpGs are from genes associated with at least one ERα binding peak, for which we expect comparable global effect (i.e., effect of ERα binding on genes) of ERα binding. This enables us to investigate the local effect of ERα binding on DNA methylation by comparing ERα binding CpGs and non-ERα binding CpGs.

In a similar way, we defined the DHS-associated CpGs and non-DHS associated CpGs. Based on the binding data of other TFs, we defined TFBS associated and non-TFBS associated CpG sites separately for each TF with ChIP-seq data.

Correlation of DNA methylation with gene expression

Gene expression data and DNA methylation data are available for 222 of the TCGA breast cancer samples. In this data, we investigate the correlation between DNA methylation of CpGs with gene expression. For each CpG, we calculated the Spearman correlation coefficient of its β values with the expression levels of its associated gene across all the samples. The correlation between methylation level of a CpG and the expression level of ESR1 was calculated in the same way.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Acknowledgments

This work was supported by the American Cancer Society Research Grant, #IRG-82-003-27, and by the start-up funding package provided to C.C. by the Geisel School of Medicine at Dartmouth College.

10.4161/epi.27688