Estrogen receptor 1 chromatin profiling in human breast tumors reveals high inter-patient heterogeneity with enrichment of risk SNPs and enhancer activity at most-conserved regions

Abstract

Estrogen Receptor 1 (ESR1; also known as ERα, encoded by ESR1 gene) is the main driver and prime drug target in luminal breast cancer. ESR1 chromatin binding is extensively studied in cell lines and a limited number of human tumors, using consensi of peaks shared among samples. However, little is known about inter-tumor heterogeneity of ESR1 chromatin action, along with its biological implications. Here, we use a large set of ESR1 ChIP-seq data from 70 ESR1⁺ breast cancers to explore inter-patient heterogeneity in ESR1 DNA binding to reveal a striking inter-tumor heterogeneity of ESR1 action. Of note, commonly shared ESR1 sites show the highest estrogen-driven enhancer activity and are most engaged in long-range chromatin interactions. In addition, the most commonly shared ESR1-occupied enhancers are enriched for breast cancer risk SNP loci. We experimentally confirm SNVs to impact chromatin binding potential for ESR1 and its pioneer factor FOXA1. Finally, in the TCGA breast cancer cohort, we can confirm these variations to associate with differences in expression for the target gene. Cumulatively, we reveal a natural hierarchy of ESR1-chromatin interactions in breast cancers within a highly heterogeneous inter-tumor ESR1 landscape, with the most common shared regions being most active and affected by germline functional risk SNPs for breast cancer development.

Figure 1.

The largest inter-patient heterogeneity in ESR1 chromatin binding is found at putative enhancers. (A) Graphical representation of study design. ESR1 ChIP-seq on tumor samples from 30 male and 40 female breast cancer patients analyzed for the level of overlap and biological features. For sample details, see Supplemental Table S1. (B) Percentage of ESR1 peaks included or excluded in consensus, by varying the threshold of minimal overlap of peaks between female patients. (C) Genomic distribution of ESR1 consensus by varying threshold in females. (D) Percentage of distal and proximal regions retained by varying threshold for consensus in females. (E) ESR1 binding sites in the vicinity of FOXA1, showing the number of patients in which these peaks were called. Green lines represent enhancer regions; red line indicates promoter. Enhancer regions were coupled to FOXA1 on the basis of work by Corces et al. (2018). Gray lines represent peaks that were not coupled to FOXA1 on the basis of work by Corces et al. (2018), but these are shown for completeness as they were located in between peaks that were coupled to FOXA1.

Figure 2.

Characterization of enhancers ranked from commonly to less frequently bound by ESR1 shows distinct biological features. (A) A ranked overview of 74,438 distal ESR1 peaks showing in how many tumor samples each peak was found in a cohort of 40 female patients. Heatmap showing the average ESR1 ChIP-seq score at a specific peak for each sample. The bar plot (left) indicates the fraction of peaks found in each patient of the total peaks found. Clustering is based on the Pearson correlation at ESR1 peaks for the ESR1 ChIP-seq signal as defined in Supplemental Figure S2A. (B) Examples of ESR1 peaks that were peak-called in tumor samples in all 40 females (left), in 16 females (middle), and in only one female patient (right). (C) Examples of per-patient heatmaps of ESR1 signal, of peaks called in that female patient sample, ranked as in A. (D) For more commonly occurring and unique peaks, examples of the average intensity of ESR1 ChIP-seq signal in four female patients are shown. (E) Correlation plot of the total number of distal peaks in a patient sample (x-axis), versus the percentage of patient-unique peaks in that sample (y-axis).

Figure 3.

ESR1 heterogeneity in female patients is not explained by molecular and clinical features. (A) Cumulative percentage of ESR1 peaks shared among female patients within the five clusters defined in Supplemental Figure S2A. (B) Percentage of female patients sharing ESR1 within groups of patients based on outcome, PGR, and ERBB2 status. The global distribution over all the peaks is depicted by a dashed green line. (C) Cumulative percentages and heatmap of the overlaps between ranked female ESR1 patients and good/poor outcome-associated (Ross-Innes et al. 2012) or aromatase inhibitor (AI) response–associated ESR1 peaks (Jansen et al. 2013).

Figure 4.

ESR1 female peaks converge to redundant enhancers regulating estrogen response genes. (A) Heatmap shows the number of ESR1 peaks that are overlapping with a region associated to a gene (x-axis) (Corces et al. 2018) per each female patient (y-axis). Each gene is ranked by decreasing number of patients carrying ESR1 peaks associated to that specific gene. The number of patients sharing a gene is shown by the line above the heatmap. The global distribution of ESR1 peak conservation among samples is depicted by a black dashed line. Ranked genes are grouped in seven bins depending on the degree of coregulation among patients. For each bin, the statistically significantly enriched cancer hallmark gene sets are shown (bottom heatmap), and the bar plot on the bottom left shows the number of bins sharing a given hallmark. The left heatmap depicts the cancer hallmarks enriched in each patient; above this heatmap, a bar plot indicates the percentage of patients showing the enrichment of each hallmark. (B) Same heatmap as in A, but in this case, the gene associated is based on chromatin loops identified by ESR1 ChIA-PET in MCF-7 (Fullwood et al. 2009; The ENCODE Project Consortium 2012).

Figure 5.

Common ESR1 peaks are associated with stronger ERE motif, increased chromatin interactions, and higher enhancer activity. (A) The percentage of common, less common, and patient-unique ESR1 peaks in females that contain an estrogen response element (ERE). (B) The strength of those EREs as determined by HOMER, ranked from those in common to those in more patient-unique ESR1 peaks. Black dots represent outliers. (C) Aggregate region analyses (ARAs) showing the average Hi-C contacts (observed over expected scores) at ESR1 binding sites shared by an increasing number of patients from left to right. The matrices include a window of ±250 kb from the ESR1 peak centers. (D) Schematic overview of STARR-seq methodology. (E) Stacked bar plot showing the overlap between STARR-seq regions and MCF-7 ESR1 peaks (Ross-Innes et al. 2012) in bins of female patient STARR-seq shared regions. (F) Volcano plot of STARR-seq results in the cell line MCF-7 upon 6 h of 10 nM estradiol (E2) stimulation. (G) Distribution of enhancer activity as determined by STARR-seq upon 6 h of 10 nM estradiol stimulation, from common to more patient-unique peaks. Details on cutoffs for categories induced, not-induced, and inactive are described in the Methods section.

Figure 6.

ESR1⁺ breast cancer rSNPs are enriched at regions with low inter-patient heterogeneity in ESR1. (A) Manhattan plot of ESR1⁺ breast cancer risk SNPs (rSNPs) with genome-wide significance originating from Michailidou et al. (2017). Highlighted in orange are 318 rSNPs, for which the coordinates intersect with one of the 74,438 ESR1 peaks found among 40 female breast cancer patients. (B) The position of these 318 rSNPs in the ranked peaks introduced in Figure 2A. (C, top) Comparison (Fisher's exact test) of the percentage of ESR1 peaks with which coordinates overlap with at least one rSNP coordinate, for common and less common ESR1 peaks. (Bottom) Comparison (Fisher's exact test) of the percentage of bases, present in common or less common ESR1 peaks, that overlap with at least one rSNP coordinate. (D) Correlation between the P-value of rSNP (x-axis) and its position in the ranking of ESR1 peaks introduced in Figure 2A (y-axis). If multiple rSNPs overlapped the same ESR1 peak, the strongest P-value was used for analysis. (E) Overview of beta values corresponding to rSNPs with which a coordinate intersected an ERE. Negative beta values correspond with rSNPs that confer less risk to breast cancer, whereas positive beta values correspond to increased risk of ESR1⁺ breast cancer.

Figure 7.

rs9952980 affects SLC14A2 expression via reduced ESR1 binding by impacting ERE. (A) Snapshots of ESR1 peak intersecting the coordinate of rs9952980. The peak, positioned in an intron of SLC14A2, was found in 11 female patients. (B) ERE at this peak, in reference allele and rSNP format. (C) Predicted score of position weight matrix for WT and rSNP ERE, by SNP2TFBS (Kumar et al. 2017). (D) Using MCF-7 lysate, an immunoprecipitation (IP) was performed with 50-bp biotin-labeled oligos containing the WT or the rs9952980 variant of the ERE, followed by mass spectrometry. (E) ESR1 western blot (WB) of IP by 50-bp biotin-labeled oligos containing the reference allele or the rs9952980 variant of the ERE. (F) Snapshot of STARR-seq normalized signal at the SLC14A2 locus. (G) Luciferase reporter assay in MCF-7 cells stimulated or not by estradiol (E2) for the SLC14A2 locus enhancer activity containing or not the rs995290 variant. Bar plot represents the fold change of luciferase expression over the untreated empty vector condition. (H) TCGA gene expression of SLC14A, which rs9952980 is predicted to affect, by homozygous or heterozygous genotype.

Figure 8.

rs6420415 affects CDYL2 expression via reduced FOXA1 binding by impacting the forkhead motif. (A) Snapshots of ESR1 peak intersecting the coordinate of rs6420415. The peak, positioned in an intron of CDYL2, was found in four female patients. (B) Forkhead motif at this peak, in reference allele and rSNP format. (C) Predicted score of position weight matrix for reference allele and rSNP forkhead motif by SNP2TFBS (Kumar et al. 2017). (D) FOXA1 western blot of pulldown by 50-bp biotin-labeled oligos containing the WT or rs6420415 variant of the forkhead motif. (E) Distribution of reads in the ESR1 and FOXA1 ChIP-seq peak performed on tumor tissue from the same breast cancer patient, at the locus surrounding rs6420415. (F) Snapshot of STARR-seq normalized signal at the CDYL2 locus. (G) Luciferase reporter assay in MCF-7 cells stimulated or not by estradiol (E2) for the CDYL2 locus enhancer activity containing or not the rs6420415 variant. Bar plot represents the fold change of luciferase expression over the untreated empty vector condition. (H) TCGA gene expression of CDYL2, which rs6420415 is predicted to affect, by homozygous or heterozygous genotype for rs6420415.