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 exist1-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 interactions4-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 prognosis7
, 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 aggressive8,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, MACS10
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 (). Peak calling details can be found in Supplementary figure 4
. An example of a core ER binding event is shown in . 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 () 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) (). The genes near (within 20kb: an optimal window between ER binding events and target genes11
) 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’ ( shows the results based on one study12
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 (more ...)
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 (). 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 () and 599 ER binding regions with more ER binding in the good outcome patients, when compared to the poor/met patients (). 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 (). 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 (more ...)
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
. 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 patients12
, 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 dataset12
and with the expected opposite directionality (). 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 pathways13
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 ( and Supplementary figure 15
). The majority (98.9%) of the core ER binding events that occurred in most primary tumours () 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 (). Examples of differentially bound regions are shown in . Using the differential ER binding events, the cell line classification can be visualised in principal component analysis () and in hierarchical clustering (). Enriched motif analysis revealed ERE and FoxA1 motifs in regions showing increased ER binding in tamoxifen resistant cell lines (), which are the same motifs observed in the poor outcome ER binding events in primary tumours (). GATA motifs were enriched in ER binding events depleted during acquisition of drug resistance (), 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 (more ...)
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 ) 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 () 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 response14,15
and have been shown to influence ER binding potential and gene expression profiles14
. 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 (). 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 (more ...)
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 (). Since FoxA1 is a pioneer factor required for ER-chromatin interactions4
and FoxA1 motifs were enriched in both the mitogen induced ER binding events (Supplementary figure 16
) and the tumours from the poor outcome patients (), 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 (). 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 (), a level of concordance that mirrors the ~50% overlap observed in wild type cells6
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 (). We found that ~87% of the metastases retained ER expression and that FoxA1 expression occurred in ~95% of the metastases (). 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.