DNA methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer

1The steroid hormone oestrogen activates the oestrogen receptor (ESR1) to mediate a variety of functions that are central to the normal development and maintenance of multiple tissues. The unique transcriptional response to oestrogen in each tissue-specific cell subtype is, in part, regulated by the epigenome. Differential DNA methylation and chromatin remodelling serve to dictate accessibility to functional, oestrogen-responsive regions of the genome, and thus define endocrine response. Inappropriate activation of the ESR1 signalling network in mammary epithelial cells initiates neoplastic transformation and drives ESR1-positive breast cancer. Patients with this disease commonly receive adjuvant endocrine therapy, which serves to inhibit ESR1 signalling. Although endocrine therapy reduces the risk of disease recurrence, a third of patients acquire drug resistance and experience disease relapse. Thus, endocrine sensitivity of both normal breast cells and breast cancer cells is dynamic, raising the hypothesis that global epigenetic reprogramming of oestrogen-responsive regions of the genome can modulate endocrine sensitivity and contributes to the onset of ESR1-positive breast cancer and the acquisition of endocrine resistance.

2While recent studies have provided excellent proof of principle that the DNA methylation profile of mammary epithelial cells is altered in early carcinogenesis, and further modified in cell models of endocrine-resistant breast cancer, they do not address how these changes could directly affect endocrine sensitivity. Here we identify DNA methylation as a key determinant of endocrine response in breast cancer. We show that differential DNA hypermethylation occurs predominantly at oestrogen-responsive enhancer, not promoter regions, and is associated with reduced ESR1 binding and decreased gene expression of key regulators of ESR1 activity. In addition, we demonstrate that the methylation status of these regulatory regions is associated with endocrine resistance in human disease, thus providing a novel mechanism by which endocrine response is abated in ESR1-positive breast cancers.

3To interrogate DNA methylation remodelling as a critical component of acquired endocrine resistance, we performed methylation profiling in duplicate using the Infinium HumanMethylation 450 beadchip, on ESR1-positive hormone sensitive MCF7 cells, and three different well-characterized endocrine-resistant MCF7-derived cell lines; tamoxifen-resistant (TAMR), fulvestrant-resistant (FASR) and oestrogen deprivation-resistant (MCF7X) cells. Density plots showing the correlation between the DNA methylation profile of parent MCF7 cells and individual endocrine-resistant cell lines indicate that the MCF7X and TAMR cells, which are both ESR1 positive, predominantly gained DNA methylation as indicated by the increased density of points above the trend line. In contrast, FASR cells, which are ESR1 negative, exhibited both hyper and hypomethylation events relative to parent MCF7 cells as indicated by a symmetrical density distribution (Fig. 1a–c). We first sought to identify the common differential DNA methylation events present in each of the three uniquely derived endocrine-resistant cell models by carrying out paired analyses (that is, each endocrine-resistant cell line versus MCF7 parent control) and overlapping the data (Fig. 1d). We found that across the individual resistant cell lines, 14,749 CpG probes were commonly hypermethylated (false discovery rate, FDR<0.01), whereas only 192 probes exhibited shared hypomethylation (FDR<0.01; Fig. 1d).

4To comprehensively characterize the functional genomic location of differential methylation observed in the endocrine-resistant cell models, we used ChromHMM segmentation of the MCF7 genome (previously described in Taberlay et al.; Fig. 1e). Strikingly, significant enrichment of commonly hypermethylated probes was exclusively observed in enhancer regions of the genome (n=3,932 probes, P<<0.0001; hypergeometric test; Fig. 1e). We next sought to determine whether the enhancer regions identified as being more heavily methylated in all endocrine resistance models were regulated by the ESR1 in the parental MCF7 cells. Using reprocessed, publically available MCF7 ESR1 (ref. ), GATA3 (ref. ) and FOXA1 ChIP-Seq data (two transcription factors closely associated with ESR1 activity), we found that enhancer-specific CpG-hypermethylated probes were enriched in ESR1-binding sites by approximately sixfold, FOXA1-binding sites by fivefold and GATA3-binding sites by eightfold (P<<0.0001; hypergeometric test; Fig. 2a). The greatest number of hypermethylated enhancer probes were found to overlap ESR1-binding sites (n=801), which represents ∼20% of all hypermethylated probes in enhancer regions. Significantly, 47% (379 out of 801) of the hypermethylated enhancer probes that were located within an ESR1-binding site were also located within a FOXA1 and/or GATA3-binding site (Fig. 2b), which is particularly noteworthy since these transcription factors cooperatively modulate ESR1-transcriptional networks by forming a functional enhanceosome.

5Having defined a subset of ESR1-binding sites that overlap enhancer regions which contain hypermethylated loci in multiple models of endocrine resistance (see Methods section; n=856 sites, Supplementary Data 1), we sought to determine whether DNA methylation affected the intensity of ESR1 binding at these sites. Using MCF7 and TAMR ESR1 ChIP data, we compared the change in ESR1 binding signal intensity at ESR1-enhancer sites that contained (a) hypermethylated probe(s) to that of all other ESR1-enhancer sites (Fig. 2c). At methylated ESR1-enhancer sites, there was a 2.29-log-fold reduction in ESR1 binding in TAMR compared with MCF7 cells. In contrast, at all other ESR1-enhancer-binding sites, there was a 0.52-log-fold reduction in ESR1 binding in TAMR compared with MCF7 cells. Thus, increased methylation at ESR1-enhancer sites is associated with reduction in ESR1 binding (P<<0.0001; t-test; Fig. 2c). Four illustrative examples show the loss of ESR1 binding in the TAMR cells at enhancer regions that are more heavily methylated in the endocrine-resistant versus the parent MCF7 (Fig. 2d). The examples include enhancer regions located within the gene body of death-associated protein 6 (DAXX), golgi to ER traffic protein 4 homologue (GET4; a member of the BAG6-UBL4A-GET4 DNA damage response/cell death complex), ESR1 itself and nuclear receptor co-repressor 2 (NCOR2; Fig. 2d).

6Since the vast majority of ESR1-enhancer-binding sites identified as hypermethylated in the endocrine-resistant cell lines compared with the parent MCF7 cells were intragenic (that is, 617 out of 856, 72% with at least partial overlap; Supplementary Data 1), we next sought to determine if the DNA methylation of these regions correlated with the expression of the genes in which they were located (or closest TSS if intergenic) in human breast cancer. Using RNA-seq and HM450 methylation data derived from TCGA breast cohort (n=459 patients), we determined that out of the 856 ESR1-enhancer-binding sites of interest, hypermethylation of 328 sites (that is, 38% of ESR1-enhancer sites) correlated with the reduced expression of the genes with which they were most closely associated (Spearman's correlation coefficient; P<0.001; Supplementary Data 2). The 328 ESR1-enhancer-binding sites represented 291 unique genes (including those presented in Fig. 2d; Supplementary Data 3). Gene set enrichment analysis revealed that these genes were over-represented in gene sets upregulated by ESR1 activation, downregulated in the acquisition of endocrine resistance and gene sets lowly expressed in basal versus luminal disease, thus suggesting that such genes were critical drivers of oestrogen-driven tumours (Supplementary Fig. 1a). Interestingly, using unsupervised clustering analysis, this gene set (n=291) stratifies ESR1-positive and ESR1-negative breast cancer patients (Supplementary Fig. 1b). Cumulatively, this indicates that the methylation events occurring throughout the acquisition of endocrine resistance are serving to facilitate an oestrogen-independent phenotype reflective of a breast cancer subtype that is refractory to endocrine therapy.

7We next sought to determine whether ESR1-enhancer hypermethylation was indicative of breast cancer subtype. We assessed the median methylation of all hypermethylated ESR1-enhancer probes (n=801) in TCGA normal (n=97), luminal A (n=301), luminal B (n=52) and ESR1-negative (n=105) patient HM450 data (Fig. 3a). In normal breast tissue (which is reported to be ∼7% ESR1 positive), the median methylation of the ESR1-enhancer sites was highest, while median DNA methylation was significantly reduced in luminal A disease (P<<0.0001; Mann–Whitney U-test), which is indicative of its endocrine-responsive state. Interestingly, median ESR1-enhancer methylation was greater in luminal B patients compared with luminal A patients (P=0.017; Mann–Whitney U-test), who are almost twice as likely to acquire endocrine resistance. In ESR1-negative disease, median methylation was higher than in luminal disease (versus luminal A, P<<0.0001; versus luminal B, P<<0.0001; Mann–Whitney U-test; Fig. 3a). A heatmap highlights the hypomethylated status of the ESR1-enhancer sites in luminal A disease relative to normal breast tissue and the other breast cancer subtypes (Fig. 3b). This trend is clearly illustrated at the DAXX enhancer region in which each CpG within the ESR1-binding site was hypomethylated in luminal A disease compared with normal tissue and luminal B and ESR1-negative cancer (Fig. 3c). Critically, no such variability was apparent at the DAXX promoter region (1,000 bp upstream and 100 bp downstream of the transcription start site; Fig. 3c), suggesting a significant regulatory effect of increased methylation at the enhancer locus.

8Given that ESR1-enhancer hypermethylation is prevalent in acquired endocrine resistance in vitro (Figs 1e and 2a-d) and in molecular subclassifications of breast cancer that are intrinsically less responsive to endocrine therapy (Fig. 3a–c), we next sought to determine the methylation status of a panel of these loci in ESR1-positive (luminal A) breast cancer samples from patients with different outcomes. Primary samples were sourced from patients who received endocrine therapy for 5 years and either experienced relapse-free survival (RFS; >14 years) or those who had relapsed (<6 years), defined as no relapse-free survival (n/RFS). Matched local relapse samples were also compared with the primary n/RFS patient samples. All patients received the same endocrine therapy (tamoxifen; anonymized patient data is given in Supplementary Data 4). Using a multiplex bisulphite-PCR resequencing methodology specifically devised for formalin-fixed, paraffin-embedded (FFPE)-derived DNA, the methylation of multiple CpG sites across a panel of nine oestrogen-responsive enhancer regions was interrogated (technical duplicate correlates for all amplicons investigated are shown in Supplementary Fig. 2). These enhancer regions included those located within DAXX, MSI2, NCOR2, RXRA and C8orf46 (Fig. 4a–e) and enhancer regions located within GATA3, ITPK1, ESR1 and GET4 (Supplementary Fig. 3a–d). The assay was repeated with DNA extracted from biological duplicates of the endocrine-resistant cell lines and the parent MCF7 cells to ensure its viability (Supplementary Fig. 4a–i; technical duplicate correlates for all amplicons investigated are shown in Supplementary Fig. 5). The average methylation levels detected at all enhancer loci were significantly higher in the recurrent tumours compared with the matched primary (n/RFS) tumours (DAXX, P<0.0001; ESR1, P<0.0005; RXRA, P<0.005; GET4, NCOR2, GATA3, MSI2, P<0.01; C8orf46, ITPK1, P<0.05; t-test), confirming that DNA methylation at ESR1-responsive enhancers is acquired in resistant disease (Fig. 4 and Supplementary Fig. 3). The difference in DNA methylation between RFS and n/RFS primary tumours was less considerable, although a statistically significant difference was observed for DAXX, P<0.0001; RXRA, P<0.01; C8orf46, P=0.01; NCOR2 and MSI2 (P<0.05; t-test) enhancer regions (Fig. 4).

9Our results support a model whereby ESR1-responsive enhancer DNA methylation is a fundamental unifying characteristic that defines endocrine sensitivity in breast cancer. Interestingly, previous studies interrogating DNA methylation changes in endocrine-resistant cell models have predominantly reported ESR1-regulated promoter methylation. Our study is the first to combine in depth MCF7 ChromHMM annotation and genome-wide methylation data from multiple resistance models to more comprehensively characterize global differential methylation across diverse genomic regions. We show for the first time that the methylation status of enhancers is associated with the inhibition of ESR1 binding in vitro and with the reduced expression of critical regulators and effectors of ESR1 activity in human disease. The identification of ESR1-responsive enhanceosome hypermethylation is both novel and considerably pertinent in the context of endocrine resistance, since genome-wide positional analyses defining the set of cis-regulatory elements that recruit ESR1 in breast cancer cells have revealed its predominant recruitment to enhancers as opposed to promoter regions. Enhancers are more common than promoters in the mammalian genome and can regulate gene transcription from tens to thousands of kilobases away by promoting communication with target promoters through chromatin looping. In our study, the majority of ESR1-regulated enhancer regions identified as hypermethylated in the resistant cells were located within gene bodies. Strikingly, hypermethylation of these enhancer regions was frequently correlated with reduced expression of the host gene, which is in line with recent studies that have shown that over half of all enhancer regions are located within a gene body and that the activation of these enhancers can indeed affect the transcription of the host gene. Examples of genes whose expression inversely correlated with ESR1-enhancer DNA methylation include DAXX and GET4, which have been previously associated with roles in apoptosis. It is conceivable that the loss of expression of genes associated with pro-apoptotic functions facilitates the progression of endocrine resistance by reducing the efficacy of apoptotic signalling pathways activated by endocrine therapies.

10Importantly, the ESR1-responsive enhancer hypermethylation events identified in the endocrine-resistant cell lines were also differentially methylated in endocrine-sensitive and endocrine-resistant breast cancer patient samples. Therefore, it is feasible that ESR1-responsive enhancer methylation status is reflective of endocrine dependence and could potentially be used to stratify patients as responders to endocrine therapy. For example, NCOR2, a gene whose expression has previously been associated with metastasis-free survival in 620 lymph node-negative patients with ESR1-positive breast cancer, was shown to negatively correlate with ESR1-enhancer methylation. In the present study, NCOR2 enhancer methylation was significantly higher in the poor (non-relapse-free) prognosis patients, compared with the good (relapse-free) prognosis primary luminal A breast cancer patients. Critically, however, in matched recurrent tumours, enhancer DNA methylation was further increased, supporting the hypothesis that the endocrine-resistant methylation profile is acquired, rather than pre-existing, limiting its potential prognostic value. Intriguingly, it could be a combination of both acquired and intrinsic methylation differences that give rise to endocrine-resistant disease. One possible explanation is that sparse, or ‘seeding' methylation at ESR1-responsisive enhancer sites in primary tumours could reflect a propensity to gain extensive methylation that spreads as resistance develops, which then becomes firmly established in recurrent disease (as discussed in ref. ). Further characterization of ESR1-responsive enhancer methylation in endocrine-resistant disease will hereafter be an important area of future investigation, as will be the assessment of its potential predictive and prognostic application in breast cancer.

11MCF7 breast cancer cells and the corresponding endocrine-resistant sub-cell lines were kindly given to our laboratory by Dr Julia Gee (Cardiff University, UK). In brief, MCF7 cells were maintained in RPMI-1640-based medium containing 5% (v/v) fetal calf serum (FCS). TAMR MCF7 cells were generated by the long-term culture of MCF7 cells in phenol red-free RPMI medium containing 5% charcoal stripped FCS and 4-OH-tamoxifen (1 × 10−7 M; TAM). FASR MCF7 cells were generated by the long-term culture of MCF7 cells in phenol red-free RPMI medium containing 5% charcoal stripped FCS and fulvestrant (1 × 10−7 M; FAS). Long-term oestrogen-deprived MCF7 (MCF7X) cells were generated by the long-term culture of MCF7 cells in phenol red-free RPMI medium containing 5% charcoal stripped FCS. Endocrine-resistant sub-lines were established and characterized following 6 months endocrine challenge/oestrogen deprivation exposure. All cell lines were authenticated by short-tandem repeat profiling (Cell Bank, Australia) and cultured for <6 months after authentication. Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue kit according to the manufacturer's instructions. HumanMethylation450K arrays were carried out by the Australian Genome Research Facility (AGRF; Melbourne, Australia).

12Two biological replicates per condition—MCF7, TAMR, MCF7X or FASR—were profiled on Illumina's HumanMethylation450K array. Raw HM450 data was preprocessed and background normalized with the Biconductor minfi package using preprocess Illumina(..., bg.correct=TRUE, normalize='controls', reference=1); resulting M values were used for statistical analyses and β values for heatmap visualizations and clustering. Differential methylation analysis of the preprocessed data was performed using the Bioconductor limma package.

13The ChromHMM segmentation of the MCF7 genome was obtained from Taberlay et al.. Enhancer (‘Enhancer' and ‘Enhancer+CTCF') and Promoter categories (‘Promoter', ‘Promoter+CTCF' and ‘Poised Promoter') were collapsed into a single ‘Enhancer' and ‘Promoter' state respectively for the purposes of our analysis. RefSeq transcript annotations were obtained from UCSC genome browser.

14ESR1 ChiP-seq data for ESR1 in MCF7 and TAMR was utilized in this study. Reads were mapped to genome build HG19 (GRCh37) with bowtie and mismatched (>3 mismatched bases), multiple mapping and duplicate reads were excluded from downstream analysis. ESR1 enrichment peaks were identified with the HOMER software suite using the findPeaks utility (-style factor -fragLength 200 -size 300 -F 0 -L 0 -C 0 -poisson 1e-06) on each experiment separately. We merged the resulting peaks to produce a ground set of 120,735 regions for subsequent analysis. Active ESR1 regions were identified in MCF7 by comparing the distribution of reads overlapping the ground set of ESR1 regions in the three MCF7 ESR1 experiments (GSM798423, GSM798424 and GSM798425) and MCF7 input experiment (GSM798440) with edegR. This yielded 54,265 active ESR1 regions in MCF7 (FDR<0.05). A similar strategy was applied to TAMR data to yield 49,511 ESR1 regions in TAMR cells. Regions of differential ESR1 binding were identified by comparing the distribution of sequence reads in MCF7 and TAMR across the ground set of ESR1 regions using edgeR and potential variation in copy number was accounted for using DiffBind. This analysis resulted in 24,711 regions with statistical significant gain (FDR 5%) and 32,343 regions with statistical significant loss (FDR 5%) of ESR1 binding in TAMR cells as compared to MCF7 cells. ESR1 peaks overlapping HM450 probes were assigned to the nearest RefSeq transcript (<20 kb distance) for the purposes of gene expression analysis. Raw MCF7 GATA3 and FOXA1 ChIP-Seq data was obtained from Theodorou et al. and Hurtado et al., respectively. Data were processed in the same manner as outlined for ESR1 ChIP-seq above.

15GSEA was performed against the Molecular Signatures Database v4.0 (MSigDB) C2 Collection. Enrichment was assessed by hypergeometric testing as implemented in the R stats package.

16DNA methylation analysis utilized clinical data available through the TCGA Breast Invasive Carcinoma cohort. Raw HM450 methylation data (level 1) were obtained from the TCGA data portal (normal samples=97, ESR1-positive tumours=353 and ESR1-negative tumours=105). ESR1-positive tumours were further divided into luminal A (lumA=301) and luminal B (lumB=52) populations using progesterone receptor (PR) expression, such that lumA were ESR1+/PR+ and lumB were ESR1+/PR−. Processed RNA-Seq expression data (level 3) were obtained from TCGA data portal (588 ESR1 positive tumours with 73 matched normals and 174 ESR1 negative samples with 19 matched normals).

17Bisulfite DNA conversions were performed using a manual protocol. For each conversion, ∼100 ng was bisulfite converted at a time. Conversion took place at 80 °C for 45 min in the presence of 0.3 M NaOH, 3.75 mM quinone and 2.32 M sodium metabisulfite, as per the method of Clark et al.. The multiplex bisulfite-PCR reaction was performed as follows. In brief, Promega HotStart GoTaq with Flexi buffer (M5005) was used with the following components at the indicated concentrations: 5 × green (1 × ), CES 5 ×, (0.5 ×, N.B. refer to ref. for CES recipe), MgCl2 (4.5 mM), dNTP's (200 μM each), primers (forward and reverse at 100 mM), Hot Start Taq (0.025 U μl−1), DNA (2 ng μl−1 ). All primers used are listed in Supplementary Data 5. Cycling conditions were: 94 °C, 5 min; 12 cycles of (95 °C, 20 s; 60, 1 min); 12 cycles of (94 °C, 20 s; 65 °C, 1 min 30 s); 65 °C, 3 min, 10 hold. Agencourt XP beads were using to clean-up and concentrate the multiplex reaction for subsequent barcoding (that is, addition of Illumina p5/p7 sequences and sample-specific DNA barcodes). The barcoding PCR used the following reagents at the indicated final concentrations in a 100-μl reaction: 1 × GoTaq Green Flexi buffer; 0.25 × CES; 4.5 mM MgCl2; 200 μM dNTPs; 0.05 U μl−1 HotStart Taq; 25 μl of pooled template after Agencourt XP bead clean-up; and 20 μl MiSeq (Fluidigm PN FLD-100-3771). Cycling conditions were: 94 °C, 5 min; 9 cycles of (97 °C, 15 s; 60 °C, 30 s; 72 °C, 2 min); 72 °C, 2 min; 6 °C, 5 min. MiSeq sequencing was performed used the MiSeq Reagent Kit v2, 300 cycle; PN MS-102-2002. Bioinformatic analysis started with adaptor trimming using Trim galore (options: --length 100). Mapping used the Bismark methylation mapping programme running Bowtie2 (ref. ) (options: --bowtie2 -N 1 -L 15 --bam -p 2 --score L,-0.6,-0.6 --non_directional; bismark_methylation_extractor -s -merge_non_CpG –comprehensive --cytosine_report). To reduce computational overhead, mapping took place against only those genomic regions which were being investigated, plus an additional 100 bp–1 kb of flanking sequence.

18FFPE breast cancer samples were obtained from the St George Hospital, Kogarah, Australia (Ethics approval reference from St George Hospital Human Research Ethics Committee is 96/84). The deidentified haematoxylin–eosin-stained sections were reviewed by a pathologist and representative tumour areas were marked and blocks were cored accordingly. Genomic DNA was extracted using the Qiagen AllPrep DNA/RNA FFPE kit according to the manufacturer's instructions.

19Accession codes: Cell line HumanMethylation450K array data is available online at GEO (GSE69118).

20How to cite this article: Stone, A. et al. DNA methylation of oestrogen-regulated enhancers defines endocrine sensitivity in breast cancer. Nat. Commun. 6:7758 doi: 10.1038/ncomms8758 (2015).