Aromatase inhibitors are the major first-line treatment of estrogen receptor–positive breast cancer, but resistance to treatment is common. To date, no biomarkers have been validated clinically to guide subsequent therapy in these patients. In this study, we mapped the genome-wide chromatin-binding profiles of estrogen receptor α (ERα), along with the epigenetic modifications H3K4me3 and H3K27me3, that are responsible for determining gene transcription (n = 12). Differential binding patterns of ERα, H3K4me3, and H3K27me3 were enriched between patients with good or poor outcomes after aromatase inhibition. ERα and H3K27me3 patterns were validated in an additional independent set of breast cancer cases (n = 10). We coupled these patterns to array-based proximal gene expression and progression-free survival data derived from a further independent cohort of 72 aromatase inhibitor–treated patients. Through this approach, we determined that the ERα and H3K27me3 profiles predicted the treatment outcomes for first-line aromatase inhibitors. In contrast, the H3K4me3 pattern identified was not similarly informative. The classification potential of these genes was only partially preserved in a cohort of 101 patients who received first-line tamoxifen treatment, suggesting some treatment selectivity in patient classification. Cancer Res; 73(22); 6632–41. ©2013 AACR.

Breast cancer is the most frequently diagnosed malignancy among women worldwide, with annually around 1.4 million new cases and half a million patients who die from the disease each year (1). Seventy-five percent of all breast tumors are of the luminal subtype and tumor cell proliferation is thought to depend on activity of the estrogen receptor α (ERα). Inhibition of ERα by endocrine therapy is therefore a major treatment modality for these tumors, either by tamoxifen or the state-of-the-art aromatase inhibitors. Although many studies have focused on defining predictive markers for tamoxifen resistance, relatively little is known about the molecular determinants of aromatase inhibitor response. Such knowledge is essential as intrinsic and acquired resistance to treatment is common (2, 3). Whole-genome sequencing analyses on breast tumor samples revealed a set of 18 genes to be mutated in patients with breast cancer that correlated with differential survival upon aromatase inhibitor treatment (4), including PI3K, TP53, MAP3K1, and GATA3. These data were supported by other reports, indicating differential gene expression pathways to be enriched in poor versus good outcome patients upon aromatase inhibitor treatment (5, 6).

ERα–binding profile assessment in breast cancer cell lines has greatly increased our knowledge of hormonal receptor action. ERα rarely binds promoters and most estrogen receptor/chromatin interactions occur at distal enhancers (7), which are involved in chromatin loop structures to regulate gene expression (8). ERα/chromatin interactions require the functional involvement of other proteins, including FOXA1. FOXA1 is one of the key luminal-defining transcription factors (9, 10) and acts as pioneer factor for ERα by untangling chromatin structures, which enables ERα to bind its targets as shown in cell lines (11) and breast tumor samples (12). Other pioneer factors have been described for ERα function, including AP-2γ (13) and PBX1 (14), which function in synergy with FOXA1 in enabling ERα/chromatin interactions and subsequent activity. These reports tightly link ERα action with chromatin structure.

The epigenetic regulation of chromatin structure is a highly complex interplay between multiple histone modifiers and their downstream histone modifications, each with their intrinsically distinct patterns and functions (15). Two of the best-studied histone marks are the repressive trimethylation on lysine 27 of histone 3 (H3K27me3; ref. 16) and the activating trimethylation on lysine 4 on histone 3 (H3K4me3; ref. 17). H3K4me3 is found enriched at promoter regions (18, 19), whereas H3K27me3 can be found over large genomic regions spanning one or more epigenetically silenced genes (20, 21), which is directly regulated by the EZH2-driven polycomb repressive complex PRC2 (22).

To test for a possible interplay between differential ERα binding and epigenetic regulation of gene activity, we determined the genome-wide chromatin-binding patterns of ERα, H3K4me3, and H3K27me3 in primary breast tumor samples. Following relapse, patients were treated with an aromatase inhibitor for metastatic disease and time to progression (TTP) was correlated with altered ERα/H3K4me3/H3K27me3 binding profiles to identify distinct patterns that could hallmark aromatase inhibitor resistance.

Patients, tumor samples, and processing

The Erasmus University Medical Center (EMC; Rotterdam, the Netherlands), the Netherlands Cancer Institute (Amsterdam, the Netherlands), and the Translational Cancer Research Unit (Saint Augustinus Hospital, Antwerpen, Belgium) participated in this study. A detailed description of the tumor samples has been described previously (23) and can be found in the Supplementary Methods.

In addition, previously described dataseries from patients with breast cancer, receiving neoadjuvant letrozole treatment (6) or tamoxifen for metastatic disease (24) were applied.

Chromatin immunoprecipitations

Chromatin immunoprecipitations (ChIP) were performed as described before (25) with minor adjustments. Firstly, as input material, tumor samples were cryosectioned (30× 30 μm sections) before further processing for ChIP-seq. For each ChIP, 10 μg of antibody was used, and 100 μL of Protein A magnetic beads (Invitrogen). ERα (SC-543; Santa Cruz), H3K4me3 (ab8580; Abcam), and H3K27me3 (07–449; Millipore) were used as antibodies. Primer sequences for quantitative PCR (qPCR) analyses are in Supplementary Table S5.

Solexa sequencing and enrichment analysis

ChIP DNA was amplified as described (25). Sequences were generated by the Illumina Hiseq 2000 genome analyzer (using 50 bp reads), and aligned to the Human Reference Genome (assembly hg19, February 2009). Enriched regions of the genome were identified by comparing the ChIP samples to mixed input using the MACS peak caller (26) version 1.3.7.1. Details on the number of reads obtained, the percentage of reads aligned, and the number of peaks called can be found in Supplementary Table S3. Bioinformatic analyses are described in the Supplementary Methods.

RNA isolation and mRNA expression analysis

Total tumor RNA was isolated as described previously (27). mRNA quality was assessed by quantitative real-time PCR (qRT-PCR) and bioanalyzer. Amplification, labeling, and hybridization of samples to 44k mRNA oligonucleotide-arrays (Agilent Technologies) were performed as described (27). Samples were hybridized against a breast cancer reference pool consisting of RNA from 105 primary breast tumors. A detailed description of mRNA expression analyses and progression-free survival analyses can be found in the Supplementary Methods.

Data access

All genomic data are deposited at the NCBI GEO, with accession numbers GSE40867 (ChIP-seq data) and GSE41994 (expression data).

Genome-wide binding patterns of ERα, H3K4me3, and H3K27me3 in primary breast tumor specimens

Fresh frozen primary tumor specimens from a cohort of 84 patients with breast cancer were tested. These patients received aromatase inhibitor treatment of metastatic disease, i.e., anastrozole, letrozole, or exemestane (Supplementary Table S2). TTP was taken as an endpoint (Fig. 1A). Poor outcome patients were defined as patients with a TTP < 12 months, whereas good outcome patients were defined as patients with a TTP > 24 months (Supplementary Tables S1 and S2). To determine the differential patterns in the ERα/chromatin–binding landscape and the epigenetic modifications H3K4me3 and H3K27me3 as well as their possible correlations with outcome after aromatase inhibitor treatment, five good outcome tumor samples and seven poor outcome tumors were randomly selected as “discovery set.” Remaining samples from the entire cohort were used as “validation sets,” as will be discussed later. The tumor cell percentage was consistently high (>70%), and all tumors were ERα and progesterone receptor (PR) positive, whereas negative for HER2. Other clinicopathologic parameters are shown in Supplementary Table S2.

Figure 1.

Genome-wide DNA-binding patterns of ERα, H3K4me3, and H3K27me3 in primary human breast tumors. A, Kaplan–Meier survival curve for the cohort of breast cancer patients who received aromatase inhibitors for metastatic disease. TTP < 12 months (red) is defined as poor outcome, whereas TTP > 24 months (green) is defined as good outcome. B, genome browser snapshot of ERα (red), H3K4me3 (blue), and H3K27me3 (green) ChIP-seq data from the same tumor sample. Genomic coordinates are indicated. Tag count is shown on the y-axis. C, Venn diagram, illustrating shared and unique chromatin-binding events between ERα, H3K4me3, and H3K27me3 from the same human breast tumor specimen. D, enriched motifs for the ERα-, H3K4me3-, and H3K27me3–binding events, as highlighted in B and C. E, genome browser snapshot examples of binding events for ERα (red), H3K4me3 (blue), and H3K27me3 (green) for all the tested tumor samples. F, genomic distributions of binding events of ERα, H3K4me3, and H3K27me3 in poor (red) and good outcome (blue) tumors, illustrating the positioning of the binding events related to the most proximal gene. Bars, SD values of all tested tumor samples.

Figure 1.

Genome-wide DNA-binding patterns of ERα, H3K4me3, and H3K27me3 in primary human breast tumors. A, Kaplan–Meier survival curve for the cohort of breast cancer patients who received aromatase inhibitors for metastatic disease. TTP < 12 months (red) is defined as poor outcome, whereas TTP > 24 months (green) is defined as good outcome. B, genome browser snapshot of ERα (red), H3K4me3 (blue), and H3K27me3 (green) ChIP-seq data from the same tumor sample. Genomic coordinates are indicated. Tag count is shown on the y-axis. C, Venn diagram, illustrating shared and unique chromatin-binding events between ERα, H3K4me3, and H3K27me3 from the same human breast tumor specimen. D, enriched motifs for the ERα-, H3K4me3-, and H3K27me3–binding events, as highlighted in B and C. E, genome browser snapshot examples of binding events for ERα (red), H3K4me3 (blue), and H3K27me3 (green) for all the tested tumor samples. F, genomic distributions of binding events of ERα, H3K4me3, and H3K27me3 in poor (red) and good outcome (blue) tumors, illustrating the positioning of the binding events related to the most proximal gene. Bars, SD values of all tested tumor samples.

Close modal

For all 12 tumors from the discovery set, fresh frozen samples were cryosectioned, and chromatin was isolated for ChIP. For each tumor sample, the ChIP for ERα, H3K4me3, and H3K27me3 was performed and isolated DNA fragments were analyzed by high-throughput sequencing (ChIP-seq; ref. 25). The number of aligned reads, unique reads, and number of peaks for each ChIP-seq sample are shown in Supplementary Table S3. Clear and distinct peaks were observed for each condition, as exemplified for one tumor sample (good outcome tumor #5; Fig. 1B). In this example specimen, 13,575 binding events for ERα, 19,012 binding events for H3K4me3, and 33,661 binding events for H3K27me3 were found. Although the vast majority of binding events found from this tissue sample were unique among the markers, overlap was found between ERα and H3K4me3 as well as between H3K4me3 and H3K27me3 (Fig. 1C; for all tumor samples, see Supplementary Fig. S1). Even though the total numbers of binding events greatly varied among tumors, these relative distributions of the three ChIP conditions were consistently found for all tumor samples tested (Supplementary Fig. S1). For the exemplified tumor sample, motif analysis was performed for each ChIP condition, Fig. 1D. As expected, ERα–binding events were enriched for ESR1 motifs, but also for its designated pioneer factors FOXA1 (7) and TFAP2 (13). For H3K4me3 ChIP on this tumor sample, enriched motifs were found for the promoter-selective transcription factors GTF2I (28), ZIC1 (29), and E2F1:TFDP2 (30), whereas only motifs for the mitochondrial transcription factor A (TFAM) were observed for H3K27me3-bound regions.

For all tested tumor samples, example binding events for ERα, H3K4me3, and H3K27me3 are shown (Fig. 1E), illustrating clear and high-quality data for all samples. For one tumor sample, no data could be generated for the H3K4me3 ChIP (poor outcome sample #2). Most ERα–binding events were found at distal enhancers and introns (Fig. 1F), consistent with previous reports in cell lines (7) and breast tumors (12). H3K4me3 was more markedly enriched at promoters, 3′-UTRs and exons, this in contrast to H3K27me3. These distributions were consistent in all analyzed tumors studied. The binding site distributions of all peaks related to the most proximal gene were found not to differ between good and poor outcome tumors.

Distinct genome-wide binding patterns of ERα, H3K4me3, and H3K27me3 correlate with patient survival after aromatase inhibitor treatment

Next, we aimed to determine whether the chromatin-binding patterns of ERα-, H3K4me3-, and H3K27me3–binding events would deviate between patients with a good versus a poor outcome upon aromatase inhibitor treatment. This was achieved through differential binding analysis (DBA; refs. 12, 31), directly comparing good outcome and poor outcome tumor ChIP-seq data. DBA normalizes the sequencing data for each run over the effective library size (reads in peaks) after input background subtraction, and determines relative enrichment of raw sequence reads over two distinct subgroups of tumor samples, in this case good versus poor outcome. Called peaks that were found in at least two tumor samples were considered for DBA to minimize noise, resulting in 14,232 peaks for ERα, 22,587 peaks for H3K4me3, and 35,602 peaks for H3K27me3 (Supplementary Fig. S2). On the basis of these peaks, the raw read counts for all tumor samples were checked for differential binding intensities between the two patient subgroups, resulting in lists of 222 (ERα), 66 (H3K4me3), and 351 (H3K27me3) peaks, with a false discovery rate of < 0.1 (Fig. 2A). Tumors from patients with the same outcome clustered together for ERα and H3K27me3 and, although less pronounced, for the H3K4me3 signals (Fig. 2B). This class separability was also observed in a principal component analysis (Supplementary Fig. S3). The read count and peak caller score did not bias patient stratification (Supplementary Fig. S4).

Figure 2.

Distinct chromatin-binding patterns of ERα and H3K27me3 in tumor samples with differential aromatase inhibitor response. A, example genomic regions with differential ERα (red), H3K4me3 (blue), and H3K27me3 (green) binding events. Genomic coordinates are indicated. Tag count is shown for each position. B, cross-correlation analysis for ERα (red), H3K4me3 (blue), and H3K27me3 (green) on the basis of the differential-bound profiles. Heatmap intensities and counts are depicted in the top left corner for each factor. C, heatmap visualization of the ERα (red), H3K4me3 (blue), and H3K27me3 (green) peak intensities that were differentially enriched in good and poor outcome tumors. Heatmap intensities and counts are depicted in the top left corner for each factor. D, pie chart depicting the relative number of peaks that were differentially enriched in good and poor outcome tumors. E, genomic locations of the differentially enriched peaks of ERα, H3K4me3, and H3K27me3 in good and poor outcome, related to the most proximal genes. F, top motif enrichment for the unique binding events in good and poor outcome tumors for all three markers. No enriched motifs were found for the H3K4me3 peaks that were only observed in the good outcome patients.

Figure 2.

Distinct chromatin-binding patterns of ERα and H3K27me3 in tumor samples with differential aromatase inhibitor response. A, example genomic regions with differential ERα (red), H3K4me3 (blue), and H3K27me3 (green) binding events. Genomic coordinates are indicated. Tag count is shown for each position. B, cross-correlation analysis for ERα (red), H3K4me3 (blue), and H3K27me3 (green) on the basis of the differential-bound profiles. Heatmap intensities and counts are depicted in the top left corner for each factor. C, heatmap visualization of the ERα (red), H3K4me3 (blue), and H3K27me3 (green) peak intensities that were differentially enriched in good and poor outcome tumors. Heatmap intensities and counts are depicted in the top left corner for each factor. D, pie chart depicting the relative number of peaks that were differentially enriched in good and poor outcome tumors. E, genomic locations of the differentially enriched peaks of ERα, H3K4me3, and H3K27me3 in good and poor outcome, related to the most proximal genes. F, top motif enrichment for the unique binding events in good and poor outcome tumors for all three markers. No enriched motifs were found for the H3K4me3 peaks that were only observed in the good outcome patients.

Close modal

Peak regions that were significantly differentially enriched between the poor and good outcome patients are shown in a heatmap for all three conditions (Fig. 2C), with a clear distinction for ERα and H3K27me3, whereas case mixing was observed for H3K4me3 signals. The number of poor outcome peaks greatly outweighs the number of good outcome events for both histone marks, but not for ERα (Fig. 2D). A list of all differentially bound regions between good and poor outcome for ERα, H3K4me3, and H3K27me3 is shown in Supplementary Table S4.

The binding events that were differentially enriched between the two patient subgroups were analyzed for their localization relative to the most proximal genes (Fig. 2E). No clear differences in genomic locations of the ERα peaks were observed between good and poor outcome patients. H3K4me3 poor outcome peaks were found enriched at 3′-UTR regions and downstream thereof, whereas the 3′-UTR signals for H3K27me3-binding events were selectively lost in poor outcome tumor samples.

The ERα-, H3K4me3-, and H3K27me3–binding events were analyzed for DNA motifs differentially enriched in the good and poor outcome patient subgroups (Fig. 2F). Although only ESR1 motifs were found in good outcome patients, TFAP2A and TCF4 motifs were enriched next to ESR1 for the poor outcome patients. For H3K27me3, good outcome–enriched binding events were enriched for TERF1 and PPARG:RXRA motifs, whereas for poor outcome–enriched binding events, motifs were found for ERG and SIX1/SIX3. No distinct enriched motifs were found for the good outcome–enriched binding event for H3K4me3, but the poor outcome sites were statistically enriched for LHX2/GBX2, ARID5B, PAX7, and HOXB8 motifs.

Mutual exclusivity of differentially enriched chromatin-binding events between the patient subgroups

A large variation was found in the number (Supplementary Table S2) and overlap of chromatin-binding events between samples as exemplified for three specimens (Fig. 3A; for all samples Supplementary Fig. S2). This is consistent with a previous report that studied ERα ChIP-seq on breast tumor specimens (12). Because ERα/chromatin interactions are dictated by and have a facilitating effect on histone accessibility (11), we determined if altered ERα–binding patterns are accompanied by changes in the epigenetic profile of H3K4me3 and H3K27me3. When ERα–binding events observed in the poor outcome tumor sample were absent in the good outcome tumor, these alterations are not accompanied by a loss or gain of proximal histone marks, as exemplified in Fig. 3B. In line with these data, differences of either histone mark between good and poor outcome patients were virtually mutually exclusive and not shared with changed ERα–binding patterns (Fig. 3C). Consequently, a limited overlap was observed for the genes that were proximal to the altered ERα-, H3K4me3-, and H3K27me3–binding events (Fig. 3D). The altered binding events were mapped over all chromosomes and no clear bias towards distinct regions or chromosomes was observed (Fig. 3E).

Figure 3.

Mutual exclusivity of altered binding events. A, Venn diagram, showing shared and unique binding events of ERα (red), H3K4me3 (blue), and H3K27me3 (green) from three example tumor samples. ERα and H3K27me3–binding patterns are highly heterogeneous between tumors, in contrast to H3K4me3-binding events. B, genome browser snapshot from a good outcome (top) and a poor outcome (bottom) tumor sample. Even though the ERα–binding sites between the two tumor samples were altered, no change was found for the present H3K4me3 (blue) and absent H3K27me3 (green) signals. C, Venn diagram, showing the shared and unique differentially bound binding events for ERα, H3K4me3, and H3K27me3, comparing good and poor outcome patients. D, as in C, but now analyzing the genes proximal to the altered binding events. All peaks within a gene body or 20 kb upstream from the transcription start site (TSS) were considered as proximal. E, visualization of the genome-wide distribution of the altered binding events of ERα (red), H3K4me3 (blue), and H3K27me3 (green) depicted over all chromosomes.

Figure 3.

Mutual exclusivity of altered binding events. A, Venn diagram, showing shared and unique binding events of ERα (red), H3K4me3 (blue), and H3K27me3 (green) from three example tumor samples. ERα and H3K27me3–binding patterns are highly heterogeneous between tumors, in contrast to H3K4me3-binding events. B, genome browser snapshot from a good outcome (top) and a poor outcome (bottom) tumor sample. Even though the ERα–binding sites between the two tumor samples were altered, no change was found for the present H3K4me3 (blue) and absent H3K27me3 (green) signals. C, Venn diagram, showing the shared and unique differentially bound binding events for ERα, H3K4me3, and H3K27me3, comparing good and poor outcome patients. D, as in C, but now analyzing the genes proximal to the altered binding events. All peaks within a gene body or 20 kb upstream from the transcription start site (TSS) were considered as proximal. E, visualization of the genome-wide distribution of the altered binding events of ERα (red), H3K4me3 (blue), and H3K27me3 (green) depicted over all chromosomes.

Close modal

Validation, integration with gene expression, and clinical outcome

Next, differential enriched binding patterns between good and poor outcome patients were validated in an independent group of 10 patients (4 good outcome and 6 poor outcome) using qPCR, Fig. 4A (for clinicopathologic parameters, see Supplementary Table S2). For each of the ChIP conditions, four to six primer pairs were designed, detecting randomly picked regions enriched in the good or poor outcome patients from the discovery set. ChIP efficiency can vary among samples due to tumor cell percentage, expression levels, and experimental variations, making inter-sample normalization an essential step in the DBA for the ChIP-seq pipeline. To implement inter-sample normalization in qPCR, the ratios of average good outcome over poor outcome–enriched binding site intensities were calculated. For ERα and H3K27me3, qPCR ratio separated poor outcome patients from good outcome patients, whereas differential enrichment could not be confirmed for H3K4me3.

Figure 4.

ChIP-qPCR validation and correlation with survival and treatment selectivity. A, heatmap, illustrating qPCR-ChIP–based enrichment of regions enriched for patients with a good or poor outcome after aromatase inhibitor treatment, normalized over negative control CCND1 primers. Average signals for the good outcome and poor outcome sites were calculated and ratio was determined, as also visualized in a bar plot (bottom). Separate primer sets were used for the altered ERα (red), H3K4me3 (blue), and H3K27me3 (green) binding events. B, Kaplan–Meier survival curves of aromatase inhibitor–treated patients, using the proximal genes from the differentially enriched binding patterns of ERα, H3K4me3, or H3K27me3, as well as the established Oncotype DX, PAM50, and TAM78 as classifiers. TTP is shown with time expressed in months. C, heatmap visualization of patient classification. Gene expression data from the Miller dataset were applied, where 54 patients received neoadjuvant letrozole treatment. Patients were stratified in responders and nonresponders, and pre- and posttreated samples from the same patients were separately analyzed. Adherence to the gene classifier is visualized in a heatmap, where green indicates a good outcome signature and red a poor outcome signature. D, identical analyses as depicted in B, but now applying the aromatase inhibitor ChIP-seq–based classifiers on a cohort of breast cancer patients who received tamoxifen for metastatic disease. E, hazard rates (HR), including 95% CI values, for the ChIP-seq–based classifiers (ERα, H3K4me3, or H3K27me3) as compared with established Oncotype DX, PAM50, and TAM78 classifiers. Both the aromatase inhibitor cohort (red) and tamoxifen cohort (blue) of patients treated for metastatic disease are shown.

Figure 4.

ChIP-qPCR validation and correlation with survival and treatment selectivity. A, heatmap, illustrating qPCR-ChIP–based enrichment of regions enriched for patients with a good or poor outcome after aromatase inhibitor treatment, normalized over negative control CCND1 primers. Average signals for the good outcome and poor outcome sites were calculated and ratio was determined, as also visualized in a bar plot (bottom). Separate primer sets were used for the altered ERα (red), H3K4me3 (blue), and H3K27me3 (green) binding events. B, Kaplan–Meier survival curves of aromatase inhibitor–treated patients, using the proximal genes from the differentially enriched binding patterns of ERα, H3K4me3, or H3K27me3, as well as the established Oncotype DX, PAM50, and TAM78 as classifiers. TTP is shown with time expressed in months. C, heatmap visualization of patient classification. Gene expression data from the Miller dataset were applied, where 54 patients received neoadjuvant letrozole treatment. Patients were stratified in responders and nonresponders, and pre- and posttreated samples from the same patients were separately analyzed. Adherence to the gene classifier is visualized in a heatmap, where green indicates a good outcome signature and red a poor outcome signature. D, identical analyses as depicted in B, but now applying the aromatase inhibitor ChIP-seq–based classifiers on a cohort of breast cancer patients who received tamoxifen for metastatic disease. E, hazard rates (HR), including 95% CI values, for the ChIP-seq–based classifiers (ERα, H3K4me3, or H3K27me3) as compared with established Oncotype DX, PAM50, and TAM78 classifiers. Both the aromatase inhibitor cohort (red) and tamoxifen cohort (blue) of patients treated for metastatic disease are shown.

Close modal

After qPCR validation, all ChIP-seq–identified differential binding patterns were coupled to genes, based on proximity (gene body plus 20 kb upstream from the transcription start-site), enriching for ERα–binding sites involved in gene regulation (8). For the good outcome binding sites, this resulted in 84 (ERα), 19 (H3K4me3), and 29 (H3K27me3) genes. For poor outcome sites, 99 (ERα), 22 (H3K4me3), and 158 (H3K27me3) genes were found (Supplementary Table S6).

Gene Ontology (GO) analyses were performed on the proximal gene sets to identify plausible functional regulatory networks (Supplementary Table S7). Genes between the different classifiers only marginally overlapped (Fig. 3D), whereas pathways between ERα, H3K4me3, and H3K27me3-based classifiers showed no overlap. The ERα–based classifier was enriched for metabolic processes, whereas nucleoside transport, chemotaxis, and angiogenesis were enriched GO terms for H3K4me3. For the H3K27me3 classifier, developmental processes were strongly enriched. Comparing our data with pathways identified in mutational analysis of aromatase inhibitor–treated patients with breast cancer (4) showed shared processes involved in cell adhesion, cell cycle, chemotaxis, developmental processes, immune responses, metabolism, signal transduction, and transcriptional regulation (Supplementary Table S8).

Next, array-based mRNA expression levels of all proximal genes for ERα, H3K4me3, and H3K27me3 sites were tested for a correlation with time to progression (TTP) in a second independent set of 72 tumors (validation set), of aromatase inhibitor–treated patients with metastatic disease (Fig. 4B, clinicopathologic parameters Supplementary Table S2). Expression of proximal genes for differential ERα (183 genes), H3K4me3 (41 genes), and H3K27me3 (187 genes) binding was used to classify tumors as “poor” and “good outcome.” Classifications for altered ERα (P = 0.0016; HR, 2.42; 95% CI, 1.40–4.19) and H3K27me3 (P = 0.0001; HR,3.19; 95% CI, 1.83–5.57) binding events correlated with a differential TTP, in contrast to H3K4me3 (P = 0.386; HR, 1.38; 95% CI, 0.76–2.49).

Previously reported classifiers, PAM50 (P = 0.0093; HR, 2.11; 95% CI, 1.20–3.71; ref. 32), Oncotype DX (P = 0.0256; HR, 1.87; 95% CI, 1.08–3.23; ref. 33), and the TAM78 Rotterdam classifier (P = 0.0151; HR, 1.98; 95% CI, 1.14–3.42; ref. 34), classified patients in this cohort, performing equally well as the ERα–proximal gene classifier (Fig. 4E). The H3K27me3-proximal gene classifier showed a higher HR, even though the 95% CI did overlap (Fig. 4E). ERα and H3K27me3 ChIP-seq–based classifiers remained significant after multivariate correction analyses (Supplementary Table S9). The classifier genes are listed in Supplementary Table S6. As a second expression-based validation, a cohort of neoadjuvant letrozole treatment patients (N = 54) was used, analyzing samples before and after treatment, where patients were stratified in responsive and non-responsive groups (Fig. 4C; refs. 5, 6). PAM50, Oncotype DX, and TAM78 as well as our ChIP-seq classifiers successfully identified patients with differential outcome.

ERα can affect gene regulation by long-range genomic chromatin-loop interactions, as shown by ERα ChIA-PET (8). Therefore, these published long-range interactions were also considered in our study. In addition, binding sites were analyzed for any transcription factor enrichment, DNAse hypersensitivity, and H3K27Ac, representing active enhancers (35), using ENCODE datasets (Supplementary Table S10). Genes were selected that either had a binding event at a promoter region, or chromatin looping toward the promoter was found, generating three separate lists of 92 (ERα), 38 (H3K4me3), and 150 (H3K27me3) genes. Overlap between the “proximity-based” and “ChIA-PET–based” gene sets was high for H3K4me3 (76%) and H3K27me3 (78%), but considerably lower for ERα (49%; Supplementary Fig. S5). Analogous to our initial classifiers, ChIA-PET–based classifiers for ERα (P = 0.0041; HR, 2.45; 95% CI, 1.41–4.26) and H3K27me3 (P = 0.0002; HR, 3.02; 95% CI, 1.73–5.27) could classify patients, in contrast to H3K4me3 (P = 0.606; HR, 1.30; 95% CI, 0.71–2.37; Supplementary Fig. S6A and Supplementary Table S9). Also for this “ChIA-PET–based” classifier, samples from the neoadjuvant treated letrozole cohort were analyzed (6), identifying patients with differential response (Supplementary Fig. S6B). Conserved genes between both gene classifiers performed equally as the separate classifications (ERα: P = 0.05; HR, 1.96, 95% CI, 1.10–3.49; H3K27me3: P = 0.0001; HR, 3.68, CI 95%, 2.09–6.50; H3K4me3: P = 0.3944; HR, 1.30, 95% CI, 0.711–2.37; Supplementary Fig. S6C and S6D).

A third cohort of patients was analyzed, receiving first-line tamoxifen for metastatic disease instead of aromatase inhibitors (Fig. 4D; ref. 27). The metastatic setup of this cohort enabled a direct comparison between the patient sets. No significant correlation with TTP after tamoxifen was found for the H3K4me3 and H3K27me3 ChIP-seq–based classifiers (H3K27me3: P = 0.0838, HR, 1.53, 95% CI, 0.95–2.47; and H3K4me3: P = 0.0602, HR, 1.556, 95% CI, 0.98–2.47). However, significance was reached for the ERα–based classifier (P = 0.0426, HR, 1.63, 95% CI, 1.02–2.63). These data indicate that while the ERα ChIP-seq–based classifier may identify breast cancer patients with a poor outcome, irrespective of the type of endocrine treatment, the H3K27me3 classification suggests aromatase inhibitor–treatment selectivity in patients with metastatic breast cancer (Fig. 4E).

Estrogen receptor genomics has greatly increased our knowledge of hormone receptor functioning in breast cancer cell lines and tumors. The number of ERα/chromatin–binding events greatly exceeds the amount of estradiol-responsive genes (11, 36), suggesting a high level of complexity in ERα–mediated gene regulation (8). ERα/chromatin–binding events have been found to correlate with survival of patients with breast cancer (12), but survival in this heterogeneous study also correlated with traditional pathologic parameters, including PR and HER2, potentially hampering further clinical interpretation of the data. Therefore, we selected ERα+, PR+, and HER2− tumors from a homogeneous cohort of breast cancer patients whose metastasis were all treated with aromatase inhibitors.

Altered ERα–binding patterns between patient groups were not accompanied with altered epigenetic profiles of H3K4me3 and H3K27me3. Our data suggest that any dynamic behavior of ERα uses the accessible regions in the genome that are imprinted and readily accessible in a static epigenetic landscape. Still, H3K27me3-binding events enable the identification of patients with a poor outcome after aromatase inhibitor treatment, in contrast to H3K4me3. The presented data suggest that aromatase inhibitor resistance is accompanied by a specific gain of polycomb-mediated gene repression at distinct sites. Moreover, it suggests an ERα–independent mechanism of therapy resistance in ERα–positive tumors. Our data could be validated in two other cohorts of breast cancer patients, using two different technologic approaches, namely ChIP-qPCR and gene expression analysis.

Because genomic patterns of ERα and H3K27me3 were indicative for patient survival, the enriched motifs may provide clues for transcription factors involved in treatment outcome. For ERα, TFAP2 and TCF4 motifs were selectively enriched for the poor outcome patients. AP-2 can directly guide ERα/chromatin interactions (13), promotes breast cancer cell proliferation (37), and correlates with poor outcome (38). TCF4 enhances breast cancer cell invasion (39) and binds ERα (40), providing a level of cross-control between estradiol and wnt pathways (40).

For H3K27me3, “good outcome” sites were enriched for TERF1 and PPARG:RXRA motifs. TERF1 is a component of the telomere nucleoprotein complex. SNPs in TERF1 have been tested for breast cancer susceptibility and prognosis, but no correlations were found (41). PPARG:RXRA ligands can trigger breast cancer cell apoptosis (42) and PPARγ activation blocks breast cancer cell invasion (43) and induces terminal cell differentiation (44). For poor outcome–enriched H3K27me3 regions, motifs were found for ERG and SIX1/SIX3. ERG and ERα mutually repress each other's activities (45). SIX1 is expressed in breast cancer, stimulating tumor cell proliferation (46), inducing genomic instability and malignant transformation (47), correlating with poor prognosis (48). Collectively, these data highlight possible transcriptional mechanisms that may form the basis for aromatase inhibitor response.

The ERα ChIP-seq classifier was applicable for aromatase inhibitor- and tamoxifen–treated patients with breast cancer with metastatic disease. Because aromatase inhibitors and tamoxifen both affect the functionality of ERα, the genomic downstream signatures for treatment outcome could overlap as well. Because ERα is targeted by these endocrine agents, the ChIP-seq approach could aid in removing noise from expression analyses to exclusively monitor genes that are directly affected by this hormone receptor.

Testing the aromatase inhibitor ChIP-seq classifiers in tamoxifen-treated tumors showed no significant difference in TTP for H3K27me3. Aromatase inhibitor–treated patients cohorts are relatively rare, and our cohort is on metastatic disease. To compare between endocrine treatments, data from a metastatic cohort of tamoxifen-treated patients were used. Performing these ERα ChIP-seq experiments and epigenetic assessments in the adjuvant setting would enable a direct comparison with any of the large (adjuvant-treated) breast cancer patient cohorts (34, 49, 50) for more extensive in silico validations.

Tumor-intrinsic plasticity of ERα and H3K27me3 can be a hallmark of endocrine therapy resistance in breast cancer and may ultimately be applicable to guide endocrine treatment selection for patients with breast cancer.

I. Simon is employed as Senior Director, R&D in Agendia and has commercial research grant from the same. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.P.H.M. Jansen, L. Dirix, E.M.J.J. Berns, W. Zwart

Development of methodology: M.P.H.M. Jansen, W. Zwart

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.P.H.M. Jansen, E.A. Reijm, I. Simon, R. Kerkhoven, A. Velds, S. van Laere, L. Dirix, X. Alexi, J.A. Foekens, S.C. Linn, W. Zwart

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.P.H.M. Jansen, T. Knijnenburg, M. Droog, A. Velds, L. Dirix, X. Alexi, L. Wessels, S.C. Linn, E.M.J.J. Berns, W. Zwart

Writing, review, and/or revision of the manuscript: M.P.H.M. Jansen, T. Knijnenburg, E.A. Reijm, S. van Laere, L. Dirix, X. Alexi, L. Wessels, S.C. Linn, E.M.J.J. Berns, W. Zwart

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.P.H.M. Jansen, I. Simon, M. Droog, W. Zwart

Study supervision: M.P.H.M. Jansen, L. Dirix, E.M.J.J. Berns, W. Zwart

The authors thank Shan Baban and Marja Nieuwland for sample processing, Kirsten Ruigrok-Ritstier, Paul Roepman, Anieta Sieuwerts, Vanja de Weerd, Anne van Galen, Maxime Look, Marion Meijer-van Gelder, Marleen Kok, Karin Beelen, Sander Canisius, and Hilde Wuyts for their contribution and technical support, Michael Hauptmann for help with statistics, and Jason Carroll (Cambridge, UK) for critically reading the article and valuable suggestions. The authors also thank the surgeons, pathologists, and medical oncologists of the St. Clara Hospital, Ikazia Hospital, St. Franciscus Gasthuis (Rotterdam), and Ruwaard van Putten Hospital (Spijkenisse) for the supply of tumor tissues and/or clinical follow-up data collection.

This work was supported by grants from the KWF Dutch Cancer Society, A Sisters Hope and Top Institute Pharma, grants T3-108 and T3-502.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Ferlay
J
,
Shin
HR
,
Bray
F
,
Forman
D
,
Mathers
C
,
Parkin
DM
. 
Estimates of worldwide burden of cancer in 2008: GLOBOCAN 2008
.
Int J Cancer
2010
;
127
:
2893
917
.
2.
Baselga
J
,
Campone
M
,
Piccart
M
,
Burris
HA
 III
,
Rugo
HS
,
Sahmoud
T
, et al
Everolimus in postmenopausal hormone-receptor-positive advanced breast cancer
.
N Engl J Med
2012
;
366
:
520
9
.
3.
Bertelli
G
. 
Sequencing of aromatase inhibitors
.
Br J Cancer
2005
;
93
:
S6
9
.
4.
Ellis
MJ
,
Ding
L
,
Shen
D
,
Luo
J
,
Suman
VJ
,
Wallis
JW
, et al
Whole-genome analysis informs breast cancer response to aromatase inhibition
.
Nature
2012
;
486
:
353
60
.
5.
Miller
WR
,
Larionov
A
. 
Changes in expression of oestrogen regulated and proliferation genes with neoadjuvant treatment highlight heterogeneity of clinical resistance to the aromatase inhibitor, letrozole
.
Breast Cancer Res
2010
;
12
:
R52
.
6.
Miller
WR
,
Larionov
A
,
Renshaw
L
,
Anderson
TJ
,
Walker
JR
,
Krause
A
, et al
Gene expression profiles differentiating between breast cancers clinically responsive or resistant to letrozole
.
J Clin Oncol
2009
;
27
:
1382
7
.
7.
Carroll
JS
,
Liu
XS
,
Brodsky
AS
,
Li
W
,
Meyer
CA
,
Szary
AJ
, et al
Chromosome-wide mapping of estrogen receptor binding reveals long-range regulation requiring the forkhead protein FoxA1
.
Cell
2005
;
122
:
33
43
.
8.
Fullwood
MJ
,
Liu
MH
,
Pan
YF
,
Liu
J
,
Xu
H
,
Hohamed
YB
, et al
An oestrogen-receptor-alpha-bound human chromatin interactome
.
Nature
2009
;
462
:
58
64
.
9.
Perou
CM
,
Sorlie
T
,
Eisen
MB
,
van de Rijn
M
,
Jeffrey
SS
,
Rees
CA
, et al
Molecular portraits of human breast tumours
.
Nature
2000
;
406
:
747
52
.
10.
Sorlie
T
,
Perou
CM
,
Tibshirani
R
,
Aas
T
,
Geisler
S
,
Johnsen
H
, et al
Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications
.
Proc Natl Acad Sci U S A
2001
;
98
:
10869
74
.
11.
Hurtado
A
,
Holmes
KA
,
Ross-Innes
CS
,
Schmidt
D
,
Carroll
JS
. 
FOXA1 is a key determinant of estrogen receptor function and endocrine response
.
Nat Genet
2011
;
43
:
27
33
.
12.
Ross-Innes
CS
,
Stark
R
,
Teschendorff
AE
,
Holmes
KA
,
Ali
HR
,
Dunning
MJ
, et al
Differential oestrogen receptor binding is associated with clinical outcome in breast cancer
.
Nature
2012
;
481
:
389
93
.
13.
Tan
SK
,
Lin
ZH
,
Chang
CW
,
Varang
V
,
Chnq
KR
,
Pan
YF
, et al
AP-2gamma regulates oestrogen receptor-mediated long-range chromatin interaction and gene transcription
.
EMBO J
2011
;
30
:
2569
81
.
14.
Magnani
L
,
Ballantyne
EB
,
Zhang
X
,
Lupien
M
. 
PBX1 genomic pioneer function drives ERalpha signaling underlying progression in breast cancer
.
PLoS Genet
2011
;
7
:
e1002368
.
15.
Ernst
J
,
Kheradpour
P
,
Mikkelsen
TS
,
Shoresh
N
,
Ward
LD
,
Epstein
CB
, et al
Mapping and analysis of chromatin state dynamics in nine human cell types
.
Nature
2011
;
473
:
43
9
.
16.
Agger
K
,
Cloos
PA
,
Christensen
J
,
Pasini
D
,
Rose
S
,
Rappsilber
J
, et al
UTX and JMJD3 are histone H3K27 demethylases involved in HOX gene regulation and development
.
Nature
2007
;
449
:
731
4
.
17.
Bernstein
BE
,
Humphrey
EL
,
Erlich
RL
,
Scheider
R
,
Bouman
P
,
Liu
JS
, et al
Methylation of histone H3 Lys 4 in coding regions of active genes
.
Proc Natl Acad Sci U S A
2002
;
99
:
8695
700
.
18.
Guenther
MG
,
Levine
SS
,
Boyer
LA
,
Jaenisch
R
,
Young
RA
. 
A chromatin landmark and transcription initiation at most promoters in human cells
.
Cell
2007
;
130
:
77
88
.
19.
Kirmizis
A
,
Santos-Rosa
H
,
Penkett
CJ
,
Singer
MA
,
Vermeulen
M
,
Mann
M
, et al
Arginine methylation at histone H3R2 controls deposition of H3K4 trimethylation
.
Nature
2007
;
449
:
928
32
.
20.
Mikkelsen
TS
,
Ku
M
,
Jaffe
DB
,
Issac
B
,
Lieberman
E
,
Giannoukos
G
, et al
Genome-wide maps of chromatin state in pluripotent and lineage-committed cells
.
Nature
2007
;
448
:
553
60
.
21.
Barski
A
,
Cuddapah
S
,
Cui
K
,
Roh
TY
,
Schones
DE
,
Wang
Z
, et al
High-resolution profiling of histone methylations in the human genome
.
Cell
2007
;
129
:
823
37
.
22.
Hansen
KH
,
Bracken
AP
,
Pasini
D
,
Dietrich
N
,
Gehani
SS
,
Monrad
A
, et al
A model for transmission of the H3K27me3 epigenetic mark
.
Nature Cell Biol
2008
;
10
:
1291
300
.
23.
Ramirez-Ardila
DE
,
Helmijr
JC
,
Look
MP
,
Lurkin
I
,
Ruigrok-Ritstier
K
,
van Laere
S
, et al
Hotspot mutations in PIK3CA associate with first-line treatment outcome for aromatase inhibitors but not for tamoxifen
.
Breast Cancer Res Treat
2013
;
139
:
39
49
.
24.
Kok
M
,
Koornstra
RH
,
Mook
S
,
Hauptman
M
,
Fles
R
,
Jansen
MP
, et al
Additional value of the 70-gene signature and levels of ER and PR for the prediction of outcome in tamoxifen-treated ER-positive breast cancer
.
Breast
2012
;
21
:
769
78
.
25.
Schmidt
D
,
Wilson
MD
,
Spyrou
C
,
Brown
GD
,
Hadfield
J
,
Odom
DT
. 
ChIP-seq: using high-throughput sequencing to discover protein-DNA interactions
.
Methods
2009
;
48
:
240
8
.
26.
Zhang
Y
,
Liu
T
,
Meyer
CA
,
Eeckhoute
J
,
Johnson
DS
,
Bernstein
BE
, et al
Model-based analysis of ChIP-Seq (MACS)
.
Genome Biol
2008
;
9
:
R137
.
27.
Kok
M
,
Linn
SC
,
Van Laar
RK
,
Jansen
MP
,
van den Berg
TM
,
Delahaye
LJ
, et al
Comparison of gene expression profiles predicting progression in breast cancer patients treated with tamoxifen
.
Breast Cancer Res Treat
2009
;
113
:
275
83
.
28.
Makeyev
AV
,
Enkhmandakh
B
,
Hong
SH
,
Joshi
P
,
Shin
DG
,
Bayarsaihan
D
. 
Diversity and complexity in chromatin recognition by TFII-I transcription factors in pluripotent embryonic stem cells and embryonic tissues
.
PLoS ONE
2012
;
7
:
e44443
.
29.
Salero
E
,
Perez-Sen
R
,
Aruga
J
,
Gimenez
C
,
Zafra
F
. 
Transcription factors Zic1 and Zic2 bind and transactivate the apolipoprotein E gene promoter
.
J Biol Chem
2001
;
276
:
1881
8
.
30.
Liu
W
,
Tanasa
B
,
Tyurina
OV
,
Zhou
TY
,
Gassmann
R
,
Liu
WT
, et al
PHF8 mediates histone H4 lysine 20 demethylation events involved in cell cycle progression
.
Nature
2010
;
466
:
508
12
.
31.
Sharma
NL
,
Massie
CE
,
Ramos-Montoya
A
,
Zecchini
V
,
Scott
HE
,
Lamb
AD
, et al
The androgen receptor induces a distinct transcriptional program in castration-resistant prostate cancer in man
.
Cancer Cell
2013
;
23
:
35
47
.
32.
Parker
JS
,
Mullins
M
,
Cheang
MC
,
Leung
S
,
Voduc
D
,
Vickery
T
, et al
Supervised risk predictor of breast cancer based on intrinsic subtypes
.
J Clin Oncol
2009
;
27
:
1160
7
.
33.
Cobleigh
MA
,
Tabesh
B
,
Bitterman
P
,
Baker
J
,
Cronin
M
,
Liu
ML
, et al
Tumor gene expression and prognosis in breast cancer patients with 10 or more positive lymph nodes
.
Clin Cancer Res
2005
;
11
:
8623
31
.
34.
Wang
Y
,
Klijn
JG
,
Zhang
Y
,
Sieuwerts
AM
,
Look
MP
,
Yang
F
, et al
Gene-expression profiles to predict distant metastasis of lymph-node-negative primary breast cancer
.
Lancet
2005
;
365
:
671
9
.
35.
Creyghton
MP
,
Cheng
AW
,
Welstead
GG
,
Kooistra
T
,
Carey
BW
,
Steine
EJ
, et al
Histone H3K27ac separates active from poised enhancers and predicts developmental state
.
Proc Natl Acad Sci U S A
2010
;
107
:
21931
6
.
36.
Zwart
W
,
Theodorou
V
,
Kok
M
,
Canisius
S
,
Linn
S
,
Carroll
JS
. 
Oestrogen receptor-co-factor-chromatin specificity in the transcriptional regulation of breast cancer
.
EMBO J
2011
;
30
:
4764
76
.
37.
Williams
CM
,
Scibetta
AG
,
Friedrich
JK
,
Canosa
M
,
Berlato
C
,
Moss
CH
, et al
AP-2gamma promotes proliferation in breast tumour cells by direct repression of the CDKN1A gene
.
EMBO J
2009
;
28
:
3591
601
.
38.
Gee
JM
,
Eloranta
JJ
,
Ibbitt
JC
,
Robertson
JF
,
Ellis
IO
,
Williams
T
, et al
Overexpression of TFAP2C in invasive breast cancer correlates with a poorer response to anti-hormone therapy and reduced patient survival
.
J Pathol
2009
;
217
:
32
41
.
39.
Ravindranath
A
,
Yuen
HF
,
Chan
KK
,
Grills
C
,
Fennell
DA
,
Lappin
TR
, et al
Wnt-beta-catenin-Tcf-4 signalling-modulated invasiveness is dependent on osteopontin expression in breast cancer
.
Br J Cancer
2011
;
105
:
542
51
.
40.
El-Tanani
M
,
Fernig
DG
,
Barraclough
R
,
Green
C
,
Rudland
P
. 
Differential modulation of transcriptional activity of estrogen receptors by direct protein-protein interactions with the T cell factor family of transcription factors
.
J Biol Chem
2001
;
276
:
41675
82
.
41.
Varadi
V
,
Brendle
A
,
Brandt
A
,
Johansson
R
,
Enquist
K
,
Henriksson
R
, et al
Polymorphisms in telomere-associated genes, breast cancer susceptibility and prognosis
.
Eur J Cancer
2009
;
45
:
3008
16
.
42.
Bonofiglio
D
,
Cione
E
,
Qi
H
,
Pingitore
A
,
Perri
M
,
Catalano
S
, et al
Combined low doses of PPARgamma and RXR ligands trigger an intrinsic apoptotic pathway in human breast cancer cells
.
Am J Pathol
2009
;
175
:
1270
80
.
43.
Liu
H
,
Zang
C
,
Fenner
MH
,
Possinger
K
,
Elstner
E
. 
PPARgamma ligands and ATRA inhibit the invasion of human breast cancer cells in vitro
.
Breast Cancer Res Treat
2003
;
79
:
63
74
.
44.
Mueller
E
,
Sarraf
P
,
Tontonoz
P
,
Evans
RM
,
Martin
KJ
,
Zhang
M
, et al
Terminal differentiation of human breast cancer through PPAR gamma
.
Mol Cell
1998
;
1
:
465
70
.
45.
Vlaeminck-Guillem
V
,
Vanacker
JM
,
Verger
A
,
Tomavo
J
,
Stehelin
D
,
Laudet
V
, et al
Mutual repression of transcriptional activation between the ETS-related factor ERG and estrogen receptor
.
Oncogene
2003
;
22
:
8072
84
.
46.
Coletta
RD
,
Christensen
K
,
Reichenberger
KJ
,
Lamb
J
,
Micomonaco
D
,
Huang
L
, et al
The Six1 homeoprotein stimulates tumorigenesis by reactivation of cyclin A1
.
Proc Natl Acad Sci U S A
2004
;
101
:
6478
83
.
47.
Coletta
RD
,
Christensen
KL
,
Micalizzi
DS
,
Jedlicka
P
,
Varella-Garcia
M
,
Ford
HL
. 
Six1 overexpression in mammary cells induces genomic instability and is sufficient for malignant transformation
.
Cancer Res
2008
;
68
:
2204
13
.
48.
Iwanaga
R
,
Wang
CA
,
Micalizzi
DS
,
Harrell
JC
,
Jedlicka
P
,
Sartorius
CA
, et al
Expression of Six1 in luminal breast cancers predicts poor prognosis and promotes increases in tumor initiating cells by activation of extracellular signal-regulated kinase and transforming growth factor-beta signaling pathways
.
Breast Cancer Res
2012
;
14
:
R100
.
49.
Loi
S
,
Haibe-Kains
B
,
Desmedt
C
,
Lallemand
F
,
Tutt
AM
,
Gillet
C
, et al
Definition of clinically distinct molecular subtypes in estrogen receptor-positive breast carcinomas through genomic grade
.
J Clin Oncol
2007
;
25
:
1239
46
.
50.
Buffa
FM
,
Camps
C
,
Winchester
L
,
Snell
CE
,
Gee
HE
,
Sheldon
H
, et al
microRNA-associated progression pathways and potential therapeutic targets identified by integrated mRNA and microRNA expression profiling in breast cancer
.
Cancer Res
2011
;
71
:
5635
45
.

Supplementary data