Differential oestrogen receptor binding is associated with clinical outcome in breast cancer.

Oestrogen receptor-α (ER) is the defining and driving transcription factor in the majority of breast cancers and its target genes dictate cell growth and endocrine response, yet genomic understanding of ER function has been restricted to model systems. Here we map genome-wide ER-binding events, by chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq), in primary breast cancers from patients with different clinical outcomes and in distant ER-positive metastases. We find that drug-resistant cancers still recruit ER to the chromatin, but that ER binding is a dynamic process, with the acquisition of unique ER-binding regions in tumours from patients that are likely to relapse. The acquired ER regulatory regions associated with poor clinical outcome observed in primary tumours reveal gene signatures that predict clinical outcome in ER-positive disease exclusively. We find that the differential ER-binding programme observed in tumours from patients with poor outcome is not due to the selection of a rare subpopulation of cells, but is due to the FOXA1-mediated reprogramming of ER binding on a rapid timescale. The parallel redistribution of ER and FOXA1 binding events in drug-resistant cellular contexts is supported by histological co-expression of ER and FOXA1 in metastatic samples. By establishing transcription-factor mapping in primary tumour material, we show that there is plasticity in ER-binding capacity, with distinct combinations of cis-regulatory elements linked with the different clinical outcomes.

Recent technological advances have allowed mapping of oestrogen receptor-α (ER) binding events, with the goal of discovering the cis-regulatory elements and factors involved in mediating ER binding and transcription. Several genome-wide maps of ER in breast cancer cell line models exist[1-3], all showing that most ER binding events occur at distal cis-regulatory elements. Forkhead motifs are enriched within the regions bound by ER binding and numerous studies have identified the forkhead protein FoxA1 as an important pioneer factor for ER− chromatin interactions[4-6]. However, the ER mapping studies have been restricted to breast cancer cell lines, mostly the MCF-7 cell line. We sought to interrogate ER binding events, for the first time, in primary frozen breast cancer samples, to determine if ER binding is dynamic and if specific cis-regulatory elements can distinguish tumours from patients with distinct clinical outcomes.

ER ChIP-sequencing (ChIP-seq) was performed in eight ER+, PR+, HER2− primary breast tumours, representative of tumours from patients with better prognosis[7], a conclusion supported by the available long-term clinical follow-up (Supplementary figure 1). Also included were seven primary breast tumours from patients with a poor outcome (ER+ PR− HER2− or ER+ PR+ HER+), since PR− or HER2+ tumours are more likely to be aggressive[8,9]. As expected, the poor outcome patients who had long-term clinical follow-up died of breast cancer (Supplementary figure 1). Furthermore, three ER+ distant metastatic samples from women with breast cancer were included. The metastatic locations and sample preparation can be found in Supplementary figure 1. As a control, we included two breast cancer samples that were ER− (ER-α negative), but expressed high transcript levels of ER-β.

ER ChIP-seq was conducted and ER binding peaks were called using two different algorithms, MACS[10] and SWEMBL (http://www.ebi.ac.uk/~swilder/SWEMBL/), to minimise peak caller bias. The number of sequencing reads and ER binding events for each tumour is shown in Supplementary figure 2. ER binding could be mapped in all tumours, but total peak intensity and the number of identified binding events differed. Three tumours were split into two sections and ER ChIP-seq was conducted on the separate sections. We found very good concordance when comparing different sections of the same tumour (R2 = 0.954) suggesting that tumour heterogeneity did not substantially influence the ER binding signal obtained from a sample (Supplementary figure 3).

We initially assessed whether a conserved set of breast cancer ER binding events could be identified. We found a core set of 484 ER binding events that were identified in at least 75% of all the tumours, but not in either of the ER− tumours (Figure 1A). Peak calling details can be found in Supplementary figure 4. An example of a core ER binding event is shown in Figure 1B. This demonstrates that ER binding to chromatin still occurs even in tumours that are unlikely to respond to antioestrogen therapies, implying that drug resistance is not due to loss of ER binding to DNA. The average ER binding signal intensity was highest in the metastatic samples and lowest in patients with good outcome tumours, a phenomenon observed both within the 484 core ER binding regions (Figure 1C) and globally (Supplementary figure 4). These data suggest that there is an acquisition of binding signal intensity in tumours that progress towards a poorer prognosis and ultimately metastasise. The only DNA motif found enriched in the core ER binding events was an oestrogen responsive element (ERE) (Figure 1D). The genes near (within 20kb: an optimal window between ER binding events and target genes[11]) the 484 core ER binding events exhibited elevated expression in the ER+ tumours used for ChIP-seq, as compared to all other genes (data not shown) and were higher in ER+ tumours relative to ER− tumours in nine independent datasets (Supplementary figure 5). The genes are provided in Supplementary figure 6 and include classic ER target genes such as TFF1, GREB1 and RARα. A gene predictor was generated based on genes near the core ER binding events. Patients were stratified and the tumours with the highest ‘risk index’ had a poor clinical outcome when compared to the tumours with the lowest ‘risk index’ (Figure 1E shows the results based on one study[12] and additional datasets are shown in Supplementary figure 6. Only ER+ patients were considered). These conserved cis-regulatory elements and their putative target genes may be the elements that contribute to tumourigenesis and are maintained regardless of the clinical outcome of the breast cancer patient. In contrast to the primary breast cancers and metastases, we mapped ER binding in three normal human mammary glands and two normal human liver samples and found limited numbers of ER binding events, with almost no concordance in ER binding between individuals (data not shown).

Figure 1

A subset of ER binding events is conserved in primary breast tumours and distant metastases. A. Heatmap showing binding peak intensity of 484 core ER binding events that are common to primary breast tumours and distant metastases. The window represents −/+ 5kb regions from the centre of the binding events. B. Example of an ER binding event at the RAR α locus, that is present in at least 75% of primary ER+ breast tumours and metastases, but not the ER− tumours. C. Normalised average signal intensity of all core ER binding events. Also included are the ER− tumours. D. Motif analysis revealed the enrichment of oestrogen responsive elements (ERE). E. Genes within 20kb of the core ER binding events were used to generate a gene predictor that was tested in independent datasets for predictive value. Tumours were stratified according to expression of this gene signature and the top one third and bottom one third of tumours were compared. Results from one dataset is shown[12] and additional datasets are provided in Supplementary figure 6.

We sought to determine if differential ER binding events could discriminate the patients with good outcome (ER+ PR+ HER2− tumours), from patients with poor outcome or metastases (we described the combined set as poor/met tumours). After normalisation of the data to account for global differences in ER binding, Differential Binding Analysis (DBA) was used to identify ER binding events that were statistically enriched in one category or the other. This resulted in a set of ER binding events that could discriminate between the two groups when using principal component analysis (Figure 2A). In total, DBA revealed 1,192 genomic regions that had significantly more ER binding in the poor/met group, compared to the good outcome patients (Figure 2B) and 599 ER binding regions with more ER binding in the good outcome patients, when compared to the poor/met patients (Figure 2B). The clustering of the tumours based on the 1,791 differential ER binding events can be visualised in Supplementary figure 7. These findings suggest that there are specific and re-occurring cis-regulatory elements that are occupied by ER in breast cancers, but that these are different in tumours that respond, versus those that relapse and metastasise. Analysis of enriched DNA motifs identified the presence of ERE and FoxA1 motifs in the differential poor outcome ER binding events and ERE motifs in the good outcome ER binding events (Figure 2C). Correlation of the poor outcome ER binding events with known processes revealed an association with endocrine resistance and luminal B status (Supplementary figure 8).

Figure 2

ER binding profiles can discriminate between tumours from patients with different clinical outcomes. A. Principal component analysis of the 1,791 ER binding events that can discriminate between the patients with good outcome tumours and those with poor/met tumours. Included are the replicates from three tumours, which are highlighted. B. Box plot representing distribution of normalised read counts in differential ER binding events that are statistically enriched in either the patients with good outcome tumours (599 ER binding events) or the patients with poor outcome tumours and the metastases (1,192 ER binding events). The samples were pre-normalised. ** denotes p < 1 ×10−10. P-values were calculated using the Wilcoxon rank sum test. C. Enriched motifs in the poor/met tumours and the good outcome ER binding events. D. Genes within 20kb of the differentially bound ER binding event were used to generate a gene predictor. Tumours were stratified according to expression of this gene signature and the top one third and bottom one third of tumours were compared. Results from one dataset is shown[12] and additional datasets are provided in Supplementary figure 12.

To investigate if the genes near the differential ER binding events were potentially functional in breast cancer, we analysed genes within a 20 kb window around the 1,192 poor/met and 599 good outcome ER binding events. Using a training set, we generated a gene expression predictor for each of the good and poor outcome gene lists. The probability calculation and comparisons between the good and poor outcome genes is shown in Supplementary figures 9 and 10. Within the poor outcome gene list was the oncogene ErbB2 (all genes are shown in Supplementary figure 11). As expected, genes in the poor outcome predictor were preferentially up-regulated in poor outcome patients, while those in the good outcome predictor were preferentially down-regulated (Supplementary figure 10). We next tested the predictors in an independent large cohort of breast cancer patients[12], only considering ER+ tumours. Using distant metastases free survival as an endpoint, both gene sets predicted outcome (p = 3 − 10−5 for good and 3 × 10−8 for poor outcome genes) in this dataset[12] and with the expected opposite directionality (Figure 2D). The gene predictors were associated with survival in additional datasets (Supplementary figure 12) and were largely independent of histopathological factors (Supplementary figure 13). We tested 1,000 randomisations from the entire list of genes and determined that the probability that a random set of genes would yield an equally robust predictor of clinical outcome was p = 0.004. Furthermore, the good and poor gene predictors had no predictive power in four cohorts of ER negative patients (Supplementary figure 14). This suggests that the increased ER binding at distinct cis-regulatory elements is functionally and biologically relevant, resulting in altered gene expression profiles that contribute to differences in drug response and overall survival.

To validate the findings made in the tumours, we explored the possibility that ER binding events were acquired in cell line models of endocrine resistance. ER binding was mapped by ChIP-seq in three commonly used tamoxifen-responsive, ER+ breast cancer cell lines (MCF-7, T-47D and ZR75-1), and two tamoxifen-resistant, ER+, breast cancer cell lines, namely a tamoxifen-resistant MCF-7 derivative (TAM-R)[13] and BT-474 cells that are ER+ and contain the ERBB2 amplification (ER+ HER2+). Similar to BT474 cells, TAM-R cells have elevated ERBB2 pathways[13] and both represent cellular systems where increased growth factor signalling results in endocrine resistance. For all five cell lines, ER ChIP-seq was performed in at least duplicate, in asynchronous cells to recapitulate the situation observed in primary tumours (Supplementary figure 2).

Almost seven thousand (6,920) ER binding events were identified in all replicates of all cell lines (Figure 3A and Supplementary figure 15). The majority (98.9%) of the core ER binding events that occurred in most primary tumours (Figure 1A) overlapped with the cell line core ER binding events. DBA identified 8,188 ER binding events with significantly stronger binding affinity in the tamoxifen-resistant cell lines and 5,713 ER binding events that were stronger in the tamoxifen-responsive cell lines (Figure 3A). Examples of differentially bound regions are shown in Figure 3B. Using the differential ER binding events, the cell line classification can be visualised in principal component analysis (Figure 3C) and in hierarchical clustering (Figure 3D). Enriched motif analysis revealed ERE and FoxA1 motifs in regions showing increased ER binding in tamoxifen resistant cell lines (Figure 3E), which are the same motifs observed in the poor outcome ER binding events in primary tumours (Figure 2C). GATA motifs were enriched in ER binding events depleted during acquisition of drug resistance (Figure 3E), possibly due to competition between FoxA1 and GATA3, another prominent breast cancer transcription factor.

Figure 3

Identification of tamoxifen-resistant ER binding profile. A. Heatmap representing ER binding events found in all tamoxifen-responsive and tamoxifen-resistant cell lines, or those enriched in either sensitive or resistant cell lines. The window represents −/+ 5kb regions from the centre of the binding events. B. Examples of ER binding events that are unique to either tamoxifen-responsive or tamoxifen-resistant cell lines. C. Principal component analysis of differential ER binding events observed in sensitive (blue) or resistant (red) breast cancer cells. D. Hierarchical clustering of the top 1,500 differentially bound ER binding events (FDR < 0.006). E. Enriched motifs within the ER binding events that discriminate between tamoxifen-sensitive and resistant cell lines.

We hypothesised that the ER binding events induced in the tamoxifen-resistant breast cancer cell lines would be the same regions that were enriched in the poor/met clinical samples. However, 79.8% of the 1,192 ER binding events enriched in the poor/met samples (from Figure 2B) overlap with ER binding events in wild type MCF-7 cells, suggesting that the cell line models are closer to the tumours and metastases from poor outcome patients. In support of this, the 599 good outcome ER binding event observed in primary tumours (Figure 2B) overlap poorly with ER binding observed in MCF-7 cells (30.2% versus 79.8% for the poor outcome). Interestingly, MCF-7 cells (plus T-47D and ZR75-1 cells) are derived from the pleural effusion of metastatic breast cancer patients, but were established prior to tamoxifen use in the clinic. We hypothesise that MCF-7, ZR75-1 and T-47D cell lines possess an intermediate ER binding profile with the acquisition of additional ER binding regions required for resistance to antioestrogen treatment.

The differences in ER binding between sensitive and resistant contexts may be due to selection and expansion of a resistant subpopulation, or may be due to reprogramming of ER binding following specific stimuli. Growth factor pathways have long been implicated in modulating endocrine response[14,15] and have been shown to influence ER binding potential and gene expression profiles[14]. We identified various stimuli (EGF, IL-6, TNF-α and IGF-I) shown to induce increased cellular invasion and drug resistance and treated asynchronous MCF-7 cells with control or the cocktail of mitogens for 90 minutes (Figure 4A). Duplicate ER ChIP-seq replicates were performed (Supplementary figure 2).

Figure 4

ER and FoxA1 binding is dynamic and their expression correlates in metastases. A. Morphological changes in MCF-7 cells treated with a cocktail of mitogens (EGF, IGF-1, IL-6 and TNF-α) for 90 minutes. B. Principal component analysis of differentially bound ER binding events can discriminate between the control and mitogenic cocktail–treated samples in the first (and second) component. C. Heatmap showing induction and loss of ER and FoxA1 binding after 90 minutes of treatment with the mitogenic cocktail. Also shown is the normalised average signal intensity of all ER and FoxA1 binding events within the ‘gained’ ER binding events. D. Metastatic samples from different sites were stained by immunohistochemistry for ER and FoxA1. Example of ER and FoxA1 expression in a breast cancer bone metastasis. E. Graphical representation of ER and FoxA1 co-expression in metastases. Location of metastases is provided, as is the Allred score.

Differential binding analysis identified 6,089 ER binding regions that were differentially enriched (≥ 4 fold change difference, FDR < 0.1) following a 90 min treatment with the cocktail. These mitogen-induced ER binding differences could be visualised using principal component analysis (Figure 4B). Since FoxA1 is a pioneer factor required for ER-chromatin interactions[4] and FoxA1 motifs were enriched in both the mitogen induced ER binding events (Supplementary figure 16) and the tumours from the poor outcome patients (Figure 2C), we assessed whether the rapid, reprogrammed ER binding occurred at regions pre-determined by FoxA1. We repeated the mitogen treatment, but mapped FoxA1 binding by ChIP-seq and found that ~25% of the reprogrammed ER binding events (1,515) occur at regions that are already bound by FoxA1, prior to mitogen treatment (Figure 4C). A substantial proportion (37.6%) of the other reprogrammed ER binding events occur at regions where FoxA1 binding is also induced by mitogens. As such, ~53% of mitogen-induced ER binding events occur at regions pre-bound by FoxA1 or at regions that also acquire FoxA1 binding (Figure 4C), a level of concordance that mirrors the ~50% overlap observed in wild type cells[6].

To determine if FoxA1 expression was present in ER+ distant metastases, we obtained 24 metastatic samples (bone, brain and liver) from ER+ breast cancer patients and performed immunohistochemistry for ER and FoxA1 (Figure 4D). We found that ~87% of the metastases retained ER expression and that FoxA1 expression occurred in ~95% of the metastases (Figure 4E). Importantly, the concordance between ER and FoxA1 was high (R2 = 0.585), regardless of the site of metastasis. Therefore, the co-expression of ER and FoxA1 in distant metastases supports our conclusions that FoxA1 mediates ER reprogramming.

By mapping ER binding in clinical samples, we provide a first glimpse of the primary regulatory regions that contribute to differences within ER+ breast cancers, rather than secondary events such as gene expression profiles. Our findings suggest that there is plasticity in ER binding, with distinct ER binding profiles associated with clinical outcome. These differential ER binding profiles appear to be mediated by FoxA1. A remaining question is what dictates differential FoxA1 and subsequently ER binding. Possibilities include changes in the genomic landscape, alterations in co-factor levels or changes in FoxA1 structure and function, potentially by post-translational modifications. By establishing transcription factor mapping in primary samples, we find that differential ER binding patterns govern gene expression programs and are associated with clinical outcome in ER+ cancer.

Methods summary

MCF-7, ZR75-1, T-47D and BT-474 human cell lines were obtained from ATCC and grown in the relevant media. TAM-R cells[13] were a kind gift from Dr Iain Hutcheson and Prof. Robert Nicholson (Cardiff). The ER+ breast cancer tumours were obtained from the Nottingham Tenovus primary breast cancer series, Addenbrooke’s Hospital and Imperial College Healthcare NHS Trust, London, UK with appropriate ethical approval from the repositories. The malignant pericardial effusion and the two distant metastases were obtained from Imperial College Healthcare NHS Trust, London, UK. For ChIP in the tumours and metastases, the frozen sample was cut into smaller pieces prior to ChIP, which was then performed as previously described[16]. For the malignant pericardial effusion, epithelial cells were first enriched using Dynabeads conjugated with Epcam[17]. For ChIPs from cell line material, proliferating cells were cross-linked and processed for ChIP as previously described[16]. The antibodies used were anti-ER (sc-543) from Santa Cruz Biotechnologies and anti-FoxA1 (ab5089) from Abcam. Sequences generated by the Illumina Genome Analyzer were processed by the Illumina analysis pipeline version 1.6.1, and aligned to the Human Reference Genome (assembly hg18, NCBI Build 36.1, March 2008) using BWA version 0.5.5[18]. Differential binding analysis was performed using the DiffBind package[19]. For immunohistochemical analyses, ER staining was conducted using the 6F11/2 mouse monoclonal antibody (Novocastra, Leica Microsystems, Bucks, UK) and FoxA1 staining was conducted using a rabbit polyclonal antibody (ab23738) from Abcam. An Allred scoring system was used to assess staining accounting for both staining intensity and the proportion of cells stained.