Abstract
While breast cancer patients with tumors that express estrogen receptor α (ER) generally respond well to hormone therapies that block ER activity, a significant number of patients relapse. Approximately 30% of these recurrences harbor activating mutations in the ligand binding domain (LBD) of ER, which have been shown to confer ligand-independent function. However, much is still unclear regarding the effect of mutant ER beyond its estrogen independence. To investigate the molecular effects of mutant ER, we developed multiple isogenic ER-mutant cell lines for the most common LBD mutations, Y537S and D538G. These mutations induced differential expression of thousands of genes, the majority of which were mutant allele specific and were not observed upon estrogen treatment of wild-type (WT) cells. These mutant-specific genes showed consistent differential expression across ER-mutant lines developed in other laboratories. WT cells with long-term estrogen exposure only exhibited some of these transcriptional changes, suggesting that mutant ER causes novel regulatory effects that are not simply due to constant activity. While ER mutations exhibited minor effects on ER genomic binding, with the exception of ligand independence, ER mutations conferred substantial differences in chromatin accessibility. Mutant ER was bound to approximately a quarter of mutant-enriched accessible regions that were enriched for other DNA binding factors, including FOXA1, CTCF, and OCT1. Overall, our findings indicate that mutant ER causes several consistent effects on gene expression, both indirectly and through constant activity.
This study demonstrates the multiple roles of mutant ER in breast cancer progression, including constant ER activity and secondary regulatory effects on gene expression and chromatin accessibility.
Introduction
Estrogen receptor α (ER, also known as ESR1) is a ligand-inducible nuclear hormone receptor that binds estrogens. ER is expressed in roughly 70% of breast cancers and plays a key role in the development and progression of these tumors (1). Because of ER's role in the growth of ER-positive (ER+) breast cancer, these tumors are generally treated with endocrine therapies, including aromatase inhibitors (AI), selective estrogen receptor modulators, and selective estrogen receptor degraders. While hormone therapies have been effective in preventing recurrence, approximately 20% of these cancers develop resistance to hormone therapies and will eventually recur (2, 3). Recent genome sequencing efforts have established mutations in the ligand binding domain (LBD) of ER as a common path to hormone therapy resistance that are found in approximately 30% of metastatic ER+ tumors (4–6).
ER mutations are rarely found in primary ER+ breast cancers, but have been observed often in metastatic tumors and in circulating cell-free DNA particularly after treatment with AIs (4–10). ER mutations occur almost exclusively in the receptor's LBD, which is responsible for binding to estrogens and recruiting cofactors, with the majority of mutations occurring at residues Y537 and D538. Crystal structures of the mutant LBD revealed that mutant ER adopts an active conformation in the absence of estrogen binding, and several studies have observed ligand-independent transcriptional regulation, growth, and proliferation of ER-mutant breast cancer cells both in vitro and in vivo (4–6, 10–16). Large-scale alterations in transcription have also been observed in cells expressing mutant ER (17–19). RNA-sequencing (RNA-seq) experiments identified hundreds of ER target genes that are differentially expressed in mutant ER cell lines in the absence of estrogens. A number of genes were also observed as potential “novel” targets of mutant ER. These genes are not regulated by wild-type (WT) ER, but are differentially regulated in mutant ER cells compared with WT regardless of estrogen treatment. Mutant ER also exhibits ligand-independent activity in regard to DNA binding. Recent chromatin immunoprecipitation sequencing (ChIP-seq) experiments showed that mutant ER is capable of binding at many ER binding sites (ERBS) in the absence of estrogens (12, 18, 20). In addition, mutant ER may differentially bind at certain sites compared with WT ER (18).
While much has been elucidated regarding the molecular consequences of ER mutations in breast cancer, what drives these effects is less clear. Multiple studies have reported that ER mutations alter the expression of a “novel” or mutant-specific set of genes; however, the consistency of these reported gene expression changes has not been evaluated. Establishing which gene expression changes are found across multiple mutant ER models will assist in identifying molecular or cellular pathways consistently altered when ER is mutated, and reveal potential factors driving mutant ER's transcriptional impact. In addition, how these mutant-specific transcriptional changes are achieved remains unclear. One model is that mutant-specific gene expression changes are caused by mutant ER's constant activity. The constitutive activity of mutant ER likely results in the regulation of genes not typically regulated by WT ER with short-term estrogen exposure. These effects could account for many of the mutant-specific expression changes. A second model is that mutations confer neomorphic properties to ER, which could induce new ER genomic binding sites or alter its ability to regulate expression. A third model is that ER mutations cause indirect changes in expression by changing the activity of other transcription factors and/or epigenetic regulators, which, in turn, impact expression. It is unclear which of these models contribute to the novel gene expression observed in breast cancer cells harboring ER mutations.
Here, we set out to answer these outstanding questions concerning ER mutations. Using a unique CRISPR/Cas9 strategy that allowed us to introduce the Y537S and D538G mutations and an epitope tag at ER's endogenous locus, we examined the effects of mutant ER on expression and found thousands of genes consistently impacted by ER mutations across models from different studies. We discovered that approximately half of the mutant ER–specific expression changes can be attributed to constant ER activity. Mutant-specific gene expression changes could be partially explained by alterations in ER genomic binding; however, chromatin accessibility alterations not involving ER were more commonly found nearby mutant-specific genes. Motif analysis of differentially accessible regions identified transcription factors bound at these loci and a role for CTCF and OCT1 in contributing to the unique transcriptional program. Together, our results suggest that ER mutations consistently impact the expression of thousands of novel genes, partially through constant ER activity and partially by indirectly altering regulatory elements through additional transcription factors.
Materials and Methods
Cell culture
T-47D cells were obtained from the ATCC and MCF7 cells were obtained from Dr. Jennifer Richer (University of Colorado Anschutz Medical Campus, Aurora, CO). WT and mutant T-47D cells were cultured in RPMI1640 Media (Thermo Fisher Scientific) supplemented with 10% FBS (Sigma-Aldrich) and 1% penicillin–streptomycin (Thermo Fisher Scientific; full media). MCF7 WT and mutant cells were cultured in Minimum Essential Media (MEM; Thermo Fisher Scientific) with 5% FBS, 1% penicillin–streptomycin, 1 × final concentration nonessential amino acids (Thermo Fisher Scientific), and 1 nmol/L final concentration Insulin (Sigma-Aldrich; full media). For all experiments, cells were kept at 37°C with 5% CO2. Five days before all experiments requiring an 17β-estradiol (E2) induction, cells were placed in hormone-depleted media. Hormone-depleted media consisted of phenol red-free MEM (Thermo Fisher Scientific), 5% charcoal-stripped FBS (Thermo Fisher Scientific), 1% penicillin–streptomycin, and nonessential amino acids and insulin as described above for MCF7 cells or phenol red-free RPMI (Thermo Fisher Scientific), 10% charcoal-stripped FBS, and 1% penicillin–streptomycin for T-47D cells. Hormone-depleted media were changed twice during each 5-day estrogen deprivation period to ensure complete removal of estrogens from the media. After 5 days of estrogen deprivation, cells were treated with DMSO as a vehicle control (Thermo Fisher Scientific) or 10 nmol/L E2 (Sigma-Aldrich) for 1 hour for ChIP-seq and Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) experiments and for 8 hours for qPCR and RNA-seq experiments.
Generation of mutant clones
Mutant clones were generated using the CRISPR epitope tagging ChIP-seq (CETCH-seq) method as described by Savic and colleagues (21). We followed the procedures outlined by Blanchard and colleagues (Supplementary Fig. S1A; ref. 22) using the same guide RNA/Cas9 vector and WT and D538G pFETCH vectors for mutant and WT ER clone generation. For this study, we generated a Y537S pFETCH vector using the same approach as in Blanchard and colleagues, substituting the D538G gBlock with a Y537S gBlock (Supplementary Table S1). Transfections were performed on cells in full media using either Lipofectamine (Thermo Fisher Scientific) or Fugene (Promega) transfection reagents, both of which yielded similar success rates. Cells were treated with 1 μmol/L SCR7 (Xcessbio) for 3 days posttransfection to block nonhomologous end joining. Seventy-two hours posttransfection, cells were treated with G418 (Thermo Fisher Scientific) at 300 μg/mL final concentration, which was applied every 2 days in conjunction with media changes. Single-cell clones were picked and validated using the same procedures as described by Blanchard and colleagues (22), including limiting dilution plating, colony picking, Sanger sequencing, and FLAG and ER immunoblots. In all, two clones for each genotype (WT, Y537S, and D538G) were validated for both T-47D and MCF7 cell lines. All clones were authenticated as matching the correct parental origin using short tandem repeat marker analysis and found to be Mycoplasma negative (IDEXX, T-47D clones, January 2020 and MCF7 clones, August 2017). Once established, early-passage lines were preserved and aliquots were grown for no more than eight passages before being discarded.
Proliferation assay
Cells from T-47D WT and mutant clones were cultured in hormone-depleted media for 3 days before plating for proliferation assays. Cells were plated in 96-well plates at approximately 5,000 cells per well for each WT or mutant clone. Proliferation was monitored for 48 hours on the IncuCyte Zoom Live Cell Imaging Platform (Sartorius) with 10× magnification images obtained at 2-hour intervals. Confluence was measured for three replicates for each cell line in each condition. Linear regression analysis was performed on the log2 confluency percentages that were divided by the initial confluency. Significant differences in proliferation rates were determined using ANOVA to compare the slopes of confluency over time.
Gene expression analysis
See Supplementary Materials and Methods for experimental and analytic details of quantitative PCR, RNA-seq, siRNA knockdown, and reverse phase protein array (RPPA).
Chromatin analysis
See Supplementary Materials and Methods for experimental and analytic details of ChIP-seq and ATAC-seq.
Statistical analysis and graphical packages
All statistical analyses were performed in R versions 3.5.2 or 3.5.3 (23), except P values for gene ontology and pathway enrichments, which were calculated by Illumina BaseSpace Correlation Engine, and P values for motif enrichments, which were calculated by MEME suite (24). P values and statistical tests used can be found throughout the text. Heatmaps for gene expression and differential chromatin accessibility were generated using the pheatmap package in R. Heatmaps display z-scores based on reads per million that align to each gene or region. ChIP-seq results were visualized using the deepTools package (25) to compare constant and mutant-enriched and -depleted ERBSs. Distances between mutant up- or downregulated genes and mutant-enriched or -depleted ERBSs or ATAC-seq sites were calculated using a Perl script. Distance plots were generated using the R plot function.
Data access
Raw and processed data are available at the Gene Expression Omnibus under accession GSE148279.
Results
Establishment of endogenously expressed FLAG-tagged ER-mutant models
Investigation of mutant ER's molecular and phenotypic consequences relies on the development of models that faithfully recapitulate mutant ER in tumors. Many previous studies investigating mutant ER's molecular and phenotypic consequences used ectopic expression of mutant ER (4–6, 10, 14–16, 18, 26). To characterize endogenous ER mutations, we created multiple isogenic clonal lines that heterozygously express FLAG-tagged mutant or WT (control) ER from the endogenous locus (Supplementary Fig. S1A; ref. 22). We developed two clones each for WT and the two most common ER LBD mutations (Y537S and D538G) in both T-47D and MCF7 breast cancer cell lines (Supplementary Fig. S1B). Engineered lines included a FLAG epitope tag at the C-terminus to allow for downstream analyses. The heterozygous expression of FLAG-tagged mutant or WT ER and the availability of multiple clones per genotype provided a robust system to investigate mutant ER's molecular effects.
Previous studies have shown estrogen-independent expression of known ER target genes in cells expressing mutant ER (5, 12, 14, 17). To determine whether our mutant clones demonstrated similar gene regulatory effects, we measured the expression of known ER target genes using quantitative PCR in cells treated with E2 or vehicle (DMSO; Supplementary Fig. S1C). In WT clones, ER-regulated genes showed significant changes in expression upon E2 treatment, as expected. In mutant ER clones, these genes were significantly differentially expressed compared with WT controls even with no E2 treatment, while the addition of E2 in some cases further increased the expression of these genes. These results indicate that our mutant clones have ligand-independent function similar to that observed in previous studies (4, 5, 12, 14, 17).
Expression levels of thousands of genes are consistently affected by ER mutations
To investigate whether mutant ER exhibits ligand-independent gene regulation on a genome-wide scale, we performed RNA-seq on WT and mutant clones grown in hormone-depleted media and treated with E2 or DMSO for 8 hours. Principal component analysis (PCA) of RNA-seq results indicated that T-47D–mutant clones, regardless of E2 treatment, differed distinctly from WT clones in their gene expression profiles (Fig. 1A). Further analysis of T-47D RNA-seq data identified a set of genes that were differentially regulated (Padjusted < 0.05) in the WT clones with the addition of E2 (E2-regulated genes). More than half (472/771) of the E2-regulated genes in the WT clones were similarly up- or downregulated in clones from one or both mutants independently of E2 (ligand-independent genes; Fig. 1B; Supplementary Table S2). A relatively small number of E2-regulated genes were differentially expressed in the opposite direction in ER-mutant cells (85 in D538G clones and 32 in Y537S clones). Known E2-regulated genes, including CISH, SUSD3, KCNK5, and PGR, were included in the ligand-independent gene set (example of qPCR-validated CISH in Fig. 1C). MCF7 mutants also distinctly clustered apart from MCF7 WT cells based on RNA-seq data as observed via PCA (Supplementary Fig. S2A). MCF7-mutant clones exhibited similar ligand-independent effects for both the Y537S and D538G mutations. In these clones, 30% of E2-regulated genes (85/289) were regulated in a ligand-independent fashion in one or both mutants (Supplementary Fig. S2B, example of CISH in Supplementary Fig. S2C). These results support previous data showing the ligand-independent function of mutant ER.
In addition to ligand-independent gene regulation, ER-mutant cells also exhibit differential expression of a large number of genes that are not differentially regulated upon an 8-hour E2 treatment in WT ER cells. In T-47D clones, more than 1,900 genes were up- (1,117) or downregulated (785) in both the Y537S and the D538G mutations (shared mutant-regulated genes; Fig. 1B, example of qPCR-validated CASC5 in Fig. 1C). The majority of mutant-regulated genes were allele specific, with 1,569 Y537S-specific and 4,311 D538G-specific up- or downregulated genes (Fig. 1B). Similar results were seen in our MCF7-mutant clones, with more than a thousand shared mutant-specific genes and thousands of genes regulated in an allele-specific fashion (Supplementary Fig. S2B, examples in Supplementary Fig. S2C). Overall, mutant-specific genes made-up the vast majority of differentially regulated genes in our mutant clones, indicating that a major result of ER mutation, in addition to ligand-independent regulation of WT ER target genes, is the differential regulation of a novel gene expression program. Of these mutant-specific genes, nearly 30% of D538G-regulated genes (upregulated genes, P = 5.1e−55 and downregulated genes, P = 1.8e−59) and 15% of Y537S differentially regulated genes (upregulated genes, P = 2.3e−28 and downregulated genes, P = 7.04e−17) were shared across both cell lines (T-47D and MCF7, examples in Supplementary Fig. S2D), suggesting that while some mutant-specific genes are consistently observed in multiple cell lines, many may be cell line specific.
To investigate enriched pathways in mutation-regulated genes, we performed gene ontology and pathway analysis and found that mutant-specific genes in T-47D cells were highly enriched for positive roles in cell cycle and cell division or protein synthesis and processing (Fig. 1D). Mutant-specific genes in MCF7 cells were enriched for genes involved in increased cellular migration and motility, but decreased cell division (Supplementary Fig. S2E). Using RPPA analysis, we validated that a number of gene expression changes that occurred in our T-47D–mutant clones also occurred at the protein level (Supplementary Fig. S3A). Many proteins that were significantly different between mutant and WT cells were also differentially regulated at the transcriptional level in our mutant clones according to our RNA-seq results (Supplementary Fig. S3B). These proteins included cyclins and cyclin-dependent kinases, as well as DNA repair proteins, which are key to successful cell-cycle progression. In addition to cell-cycle–related genes, we also validated increased expression of IGF1R at the protein level and decreased expression of CDH1, which codes for E-cadherin (Supplementary Fig. S3B). These proteins play important roles in cell growth and cell motility, and IGF1R has been shown to be activated by ER mutations (27, 28).
To explore phenotypes associated with the observed gene expression changes, we analyzed proliferation in T-47D cells. Using live cell imaging on the IncuCyte Zoom platform, we found that ER-mutant cells grow significantly faster in full media than WT lines, with the D538G lines exhibiting the fastest growth (Supplementary Fig. S4A). In hormone-depleted media, the effects on growth were much more pronounced, with all mutant clones growing more rapidly than WT cells (Supplementary Fig. S4B). The increased growth of ER-mutant cells may be related to E2F factors, which have been shown to mediate secondary effects of estrogens in breast cancer cells (29, 30). E2F factors, E2F1, E2F2, E2F7, E2F8, and E2F dimerizing partner, TFDP1, were significantly upregulated in ER-mutant T-47D cells. In addition, motif comparison of T-47D ER-mutant upregulated and downregulated gene promoters showed strong enrichment for E2F motifs in upregulated gene promoters (P = 1.4e−19 and 1.8e−20 for Y537S and D538G). These observations suggest that E2F factors may play a role in mediating the gene expression effects of ER mutations. The gene expression changes found in the MCF7 ER-mutant cells suggested a change in migratory potential. Consistent with this observation, previous studies have shown that ER-mutant MCF7 cells migrate faster in a scratch wound assay (14). We also found that our ER-mutant MCF7 cells were more capable of invading Matrigel (31). Together, these findings indicate that the gene expression changes observed in ER-mutant cells correspond to phenotypic changes. This expanded transcriptional impact may also affect other properties, such as anchorage-independent growth and metastatic outgrowth, as observed by Williams and colleagues (31).
We observed increased expression of ITGA6 (Supplementary Fig. S2D), which is also known as CD49f and is used as a marker of luminal progenitor cells (32, 33). This led us to question whether ER-mutant cells resemble luminal progenitor cells more than WT cells. We compared our mutant-specific up- and downregulated genes with genes found to be up- or downregulated in luminal progenitor cells in a previous study (33). For our T-47D-mutant–specific genes, we observed marginally significant overlaps with luminal progenitor genes. However, we found that our MCF7 mutant–specific upregulated genes exhibited highly significant overlaps with luminal progenitor increased genes (Supplementary Table S3). These results suggest that ER mutations may lead to a more stem-like luminal progenitor expression state, at least in some contexts.
Previous studies have described T-47D and MCF7 ER-mutant lines made using CRISPR or lentiviral techniques and have reported similar E2-independent and mutant-specific gene regulation (12, 17–19). In an effort to define which genes are consistently regulated by ER mutations, we compared RNA-seq data from the clones described above with RNA-seq data collected from ER-mutant lines described in two other studies (17, 18). A PCA plot of these multi-laboratory data from MCF7 lines clustered them by study more than by ER mutational status (Supplementary Fig. S5A). This is likely due to multiple factors, including differences in mutation introduction techniques, duration of hormone depletion, E2 induction time, and sequencing library generation protocol. However, within PCA study clusters, clones were segregated by mutational status and E2 treatment. Multivariate analysis of these data to account for variation between studies identified genes that were consistently differentially regulated by ER mutations across all datasets. More than 1,900 genes were up- (993) or downregulated (910) by both mutations across all MCF7 lines (Fig. 2A, examples in Fig. 2B; Supplementary Table S4). Thousands of genes were specific to each mutation (1,378 D538G-specific genes and 2,678 Y537S-specific genes). Genes that were consistently differentially expressed in mutant lines from all three studies were enriched for a variety of gene ontology terms shared between both mutations, including innate immune response, mitochondrial matrix, and extracellular matrix. Allele-specific terms were also observed, including cellular respiration and ribosomal subunit for the Y537S mutation and cell–cell adhesion and lipid biosynthesis for the D538G mutation (Fig. 2C). Although these terms differ somewhat from those identified in individual studies' lines, they represent underlying changes that consistently occur in response to the expression of ER mutations and, therefore, may provide direction for future studies regarding mutant ER action and effects.
We used WT and mutant T-47D RNA-seq data from the same three studies to perform an identical multivariate analysis as above, with the exception that one study lacked the D538G mutation in T-47D cells. The T-47D multivariate analysis yielded similar results to the MCF7 analysis, with more than 1,700 up- or downregulated genes shared between both mutations and thousands of additional genes differentially regulated by each mutant allele uniquely (Supplementary Fig. S5B, examples in Supplementary Fig. S5C). These genes were heavily enriched for cell-cycle–related genes (Supplementary Fig. S5D). Overall, the comparison of multiple lines from multiple individual groups provides evidence of consistent regulation of novel genes by ER mutations.
To determine whether the expression changes that we observed across studies can be seen in patient samples, we analyzed RNA-seq data from the MET500 study (34). By splitting metastatic breast cancer samples by ER mutation status, we identified 70 upregulated genes in ER-mutant tumors and 160 downregulated genes in ER-mutant tumors at an FDR of 5%. Comparing these gene lists with the multi-study genes described above, we observed significant overlaps of 16%–22% of upregulated genes and 10%–13% of downregulated genes (Supplementary Table S5). The overlaps were more significant between upregulated genes, with T-47D Y537S and MCF7 D538G showing the highest levels of enrichment. These results indicate that genes whose expression is impacted in well-controlled cell line experiments are enriched for differential expression in patient tumors based on ER mutation status.
Constant ER activity can explain some of mutant ER's impact on gene expression
The constant activity of ER brought about by LBD mutations could explain some of the mutant-specific gene expression changes observed in our RNA-seq analyses. Long-term constitutive ER activity may lead to the differential expression of genes that are in fact regulated by WT ER, but do not change expression with short-term E2 treatment. To determine how much of mutant ER's gene regulatory changes can be attributed to mutant ER's constitutive activity, we treated T-47D and MCF7 WT clones with E2 for 25 days. WT ER cells were collected at 5-day intervals and RNA was harvested for RNA-seq analysis. A PCA of RNA-seq data compared WT T-47D cells treated with long-term (5, 10, 15, 20, and 25 days) E2 with WT ER short-term (0 and 8 hours) E2 treatments and mutant T-47D cells. PC1, which accounted for 50% of the variance between samples, revealed a distinct separation between long-term E2-treated WT ER cells and both WT ER short-term E2 treatments and mutant ER cells. (Fig. 3A). PC2 grouped long-term E2 treatments apart from short-term E2 treatments and instead clustered them closer to Y537S- and D538G-mutant cells. Together, these observations indicate that while long-term ER activity may account for some of the differences between mutant and WT cells, it may not account for all mutant ER gene expression changes. Further analysis of RNA-seq data revealed that 54% of all mutant-regulated genes were similarly up- or downregulated by long-term E2 treatment (Padjusted < 0.05; Fig. 3B and C, example of TBCD in Fig. 3D). However, a considerable number of mutant ER–regulated genes were not differentially regulated with long-term E2 treatment in WT cells, but remained specific to one or both mutations (example of COPS2 in Fig. 3D). We performed gene ontology analysis on mutant-specific genes that were either differentially regulated or unaffected by prolonged E2 treatment in WT clones. Pathways regulated by these two groups of genes differed greatly and highlight processes affected by the constitutive activity of ER mutants or their potential novel function (Fig. 3E). Similar results were observed in MCF7 WT cells treated with long-term E2, with 38% of all mutant-specific genes being explained by constant ER activity (Supplementary Fig. S6A–S6D). These results suggest that while the constitutive activity of the mutant receptor may account for approximately half of the gene expression changes observed in mutant ER breast cancer cells, a large percentage of mutant-specific gene expression changes cannot be attributed to constant ER activity. These unexplained mutant-specific gene expression changes may instead result from additional properties of mutant ER, which allow it to directly or indirectly alter the expression of these genes.
ER mutations do not cause broad reprogramming of genomic binding
Mutant-specific gene expression changes could result from altered genomic binding by mutant ER compared with WT ER. To investigate this possibility, we performed ChIP-seq experiments on WT and mutant T-47D and MCF7 clones, grown for 5 days in hormone-depleted media followed by 1-hour E2 or DMSO treatments, using an antibody that recognizes the FLAG epitope tag. Results from these ChIP-seq experiments compared well with previously performed ER ChIP-seq experiments performed in another laboratory, with 68%–73% binding site overlaps observed between samples of the same ESR1 genotype (18). In WT T-47D clones, ER binding by WT ER was absent at thousands of genomic regions in the absence of E2, but was gained with the addition of E2 (Fig. 4A, example in Fig. 4B). ChIP-seq of mutant ER showed that, for both mutations, ER bound to the majority of WT ERBSs independent of E2 treatment and mutant ER binding at these sites increased with the addition of E2. The majority of ER-bound sites were identified in both mutants and in E2-treated WT clones (34,920 constant sites); only approximately 10% (3,837 sites) of ER-bound regions exhibited significant differential binding between WT and mutant, with the large majority of these sites being differentially bound in Y537S-mutant clones (Fig. 4C, example in Fig. 4D). In both constant and mutant-enriched or -depleted sites, the ER binding motif (ERE) was the most highly enriched motif.
To determine how much of the mutant-specific expression differences that we observed could be explained by ER binding, we analyzed the distance between mutant-specific genes and constant, mutant-enriched, or mutant-depleted ERBSs. Both Y537S- and D538G mutant–specific up- and downregulated genes were significantly closer to constant ERBS than were all other genes (background). Using a previously described distance threshold of 100 kb (35–37), we determined the number of mutant-specific genes that contained ERBS near their transcription start site (TSS). Of both Y537S and D538G mutant–specific genes, 55%–80% contained nearby constant ERBS, which was 13%–16% more than background genes. Mutant-specific up- and downregulated genes were also significantly closer to mutant-enriched and -depleted sites, respectively (Fig. 4E; Supplementary Table S6). Ten percent of Y537S up- and downregulated genes contained at least one nearby Y537S-enriched or -depleted ERBS (Supplementary Table S6). Enriched and depleted ERBSs were found to a much lesser extent in D538G-mutant clones; less than 1% of D538G-specific upregulated genes contained a nearby D538G-enriched ERBS and less than 2% of D538G-specific downregulated genes contained a nearby D538G-depleted ERBS. These results indicate that while constant ligand-independent ER genomic binding likely affects the expression of a portion of mutant-specific genes, differential ER binding, particularly in D538G-mutant cells, may only have a subtle effect on gene expression.
Similar results were seen in MCF7 WT and mutant clones where 36,606 sites were bound in WT cells treated with E2 that were also bound in mutant cells regardless of E2 treatment (Supplementary Fig. S7A–S7C). Differential ER binding in MCF7 clones also followed similar patterns to those observed in T-47D clones, with fewer D538G differentially bound ERBSs (1,161) than Y537S differentially bound ERBSs (3,020). Similar numbers of mutant-specific genes harboring an ERBS within 100 kb of the TSS were observed in MCF7 cells compared with T-47D cells (Supplementary Fig. S7D; Supplementary Table S6). Approximately 13% of Y537S-enriched and -depleted genes contained at least one nearby enriched or depleted Y537S ERBS, while only 1%–6% of D538G genes contain nearby D538G-enriched or -depleted ERBS. While ligand-independent genomic binding by mutant ER was apparent, the low percentage of ERBSs that exhibited significant differential binding suggests that changes in ER binding may not have a broad impact on the regulation of mutant-specific genes. However, mutant-specific genes were significantly closer to mutant-altered ERBS and, therefore, may account for a small portion (2%–14%) of mutant-specific differential gene expression.
ER mutations cause widespread changes in chromatin accessibility
Because ESR1 mutations do not appear to drastically reprogram ER genomic binding, we hypothesized that additional changes might occur at genomic loci bound by ER, or other transcription factors, that contribute to mutant ER's unique gene expression program. To identify genomic regions that could contribute to the observed mutant ER gene expression changes, we analyzed genome-wide chromatin accessibility using ATAC-seq in our WT and mutant T-47D clones with a 1-hour E2 or DMSO treatment. Minimal changes were observed in chromatin accessibility in response to E2 treatment, with only 169 regions exhibiting significantly increased accessibility and no regions of decreased accessibility in WT clones. However, in mutant ER cells, thousands of genomic regions exhibited increased (1,808 in Y537S clones and 1,563 in D538G clones; mutant-enriched sites; Fig. 5A, example in Fig. 5B) or decreased accessibility (2,060 in Y537S clones and 1,587 in D538G clones; mutant-depleted sites; Fig. 5A) when compared with WT clones. The changes in chromatin accessibility were largely mutation allele specific, with only 21.6% of mutant-enriched and 23.7% of mutant-depleted regions shared between both mutations.
We explored the relationship between changes in chromatin accessibility and mutant-specific gene expression by comparing ATAC-seq and RNA-seq results from our T-47D–mutant clones. We found that mutant-specific up- and downregulated genes were significantly enriched near regions of mutant-altered chromatin accessibility (Fig. 5C; Supplementary Fig. S8A). Of mutant-specific upregulated genes, 17% above background (all other genes) had at least one mutant-enriched region within 100 kb. Of mutant-specific downregulated genes, more than 15% above background contained at least one nearby mutant-depleted region (Supplementary Table S6). To illustrate the effect that altered chromatin accessibility could have on gene expression, two regions of increased chromatin accessibility near the mutant-specific CEP78 gene are shown. These regions correlate positively to the increased expression of the CEP78 gene and may function as regulatory elements controlling its expression (Fig. 5B). We next compared genes that had nearby mutant-altered ERBS with those harboring nearby mutant-altered accessible chromatin regions, and found that many Y537S-specific genes harbored both differential ER binding and accessible chromatin sites within 100 kb of their TSS (Supplementary Fig. S8B). However, chromatin accessibility changes were found much more often nearby mutant-specific genes than changes in ER binding, especially in the context of the D538G mutation. These findings suggest that while some mutant-specific gene expression changes may be a direct consequence of altered ER binding, more often these changes appear to be associated with regions of altered chromatin accessibility and potentially other transcription factors binding these regions.
To determine the extent to which ER could be directly involved at regions of altered chromatin accessibility, we overlapped ER-bound sites identified from our FLAG tag ChIP-seq experiments with all regions of differential chromatin accessibility. Of the mutant-enriched ATAC-seq regions, only 30%–50% overlapped with ERBS, while only 20%–30% of mutant-depleted ATAC-seq regions were bound by ER (Fig. 5D), indicating that many of these chromatin changes are not direct effects of mutant ER. To determine whether the results seen in our mutant cells were recapitulated in mutant lines described in previous studies, we performed ATAC-seq on T-47D WT and mutant D538G and Y537S cells first described by Bahreini and colleagues (17). We performed multivariate differential accessibility analysis on the combined ATAC-seq data from these sets of mutant ER models and identified thousands of genomic regions that showed consistent differential chromatin accessibility between mutant ER and WT lines (Supplementary Fig. S9A). ER binding again occurred at a minority (20%–30%) of the differentially accessible regions from the multivariate analysis (Supplementary Fig. S9B). These findings suggest that considerable changes in chromatin accessibility were consistently seen at thousands of sites in multiple independent ER-mutant models, and that they were largely mutation allele specific and mostly lacked ER binding.
Transcription factors associated with regions of differential accessibility in ER-mutant cells
To identify additional factors associated with mutant-enriched accessible chromatin, we performed motif enrichment analysis on mutant-enriched ATAC-seq sites, identified above, that were negative for ER binding (ERBS negative). We also excluded any sites that were proximal to gene TSS (within 2 kb) from our analysis due to the heavy overrepresentation of GC-rich regions that could bias motif analysis results. Motif analysis of mutant-enriched ATAC-seq sites resulted in a number of motifs that were significantly enriched in one or both of the mutant genotypes (Fig. 6A; Supplementary Fig. S9C). In D538G-enriched sites, the forkhead factor (FOX) motif was most highly enriched (P = 1.5e-31). Y537S-enriched sites were highly enriched for POU (OCT) factor (P = 8.0e-78) and CTCF (P = 3.3e-50) motifs. Because the FOX and OCT motifs can bind to many related transcription factors, we analyzed expression data and found that four FOX factors and three OCT factors had detectable mRNA (Fig. 6B); however, FOXA1 exhibited much higher expression than other FOX factors and OCT1 (POU2F1) was the highest expressed OCT factor. OCT1 was also more highly expressed in ER-mutant cells compared with WT. These observations led us to examine OCT1, FOXA1, and CTCF in more depth.
To investigate the potential role of these factors in mutant ER cells, we determined the extent to which these factors bind to mutant-enriched accessible chromatin. We performed ChIP-seq for OCT1 and CTCF in the T47-D WT and ER-mutant cells in hormone-depleted conditions. We also analyzed FOXA1 ChIP-seq data collected previously in T-47D WT cells (38). OCT1 binding was highly enriched in TSS-distal, mutant-enriched accessible chromatin regions in both the Y537S- (P = 2.88e-59; hypergeometric test) and D538G-mutant (P = 2.02e-85) clones. Significant OCT1 binding was found at 25% of both ERBS-negative and ERBS-positive accessible chromatin regions (Fig. 6C). CTCF showed significant enrichment in Y537S mutant–enriched, TSS-distal accessible chromatin (P = 1.59e-4) and overlapped with 35% of these sites, but was not enriched in D538G-enriched accessible regions, matching the motif pattern. CTCF was more prevalent at ERBS-negative than at ERBS-positive regions. FOXA1 was enriched in TSS-distal accessible chromatin regions in both Y537S- (P = 7.8e-52) and D538G-mutant (P = 3.13e-102) clones and was found at approximately half of both mutants' TSS-distal accessible regions, but was more evident at ERBS-positive regions than ERBS-negative regions. The enrichment for FOXA1 in ERBS-positive regions is not surprising, as it is a known pioneer factor for ER (39) and could be functioning in this role to increase chromatin accessibility in ER-mutant breast cancer cells. The binding of these factors at regions of mutant-enriched chromatin accessibility suggests that they could function in altering the expression of mutant-specific genes. Increased CTCF or OCT1 binding in mutant compared with WT cells near mutant-specific genes positively correlates with increased mutant-specific gene expression and suggests a role for these factors in altering gene expression in these cells. Overall, the chromatin accessibility analysis has led us to three transcription factors (OCT1, CTCF, and FOXA1) that bind to regions of the genome that became more accessible in ER-mutant cells, indicating that these factors may contribute to the unique gene expression patterns attributed to mutant ER.
To determine whether OCT1 and CTCF binding to regions of increased accessibility have a functional impact on gene expression, we performed OCT1- and CTCF-knockdown experiments. We selected four and six mutant-specific genes that harbored nearby regions of increased chromatin accessibility and increased CTCF or OCT1 occupancy, respectively. We performed siRNA-mediated knockdowns of OCT1 and CTCF, achieving approximately 50% protein knockdown of CTCF and 70% protein knockdown of OCT1 (Supplementary Fig. S10A). We found that upon OCT1 knockdown, the expression of four of the six genes showed a significant reduction in gene expression (Fig. 6D; Supplementary S10B). Two genes exhibited significantly reduced expression in the context of both mutations, while two genes had significant reduction in expression in the context of one mutation, although nonsignificant trends were observed for some genes. With CTCF knockdown, two genes exhibited significantly lower expression in the context of both mutations, while the other two genes showed no expression effects (Fig. 6E; Supplementary Fig. S10C). In the WT context, one gene was significantly reduced by CTCF knockdown and two genes were significantly reduced by OCT1 depletion. Loss of expression in both the WT and mutant settings indicates that these factors may have a general role in the expression of these genes regardless of ER genotype; however, we found that the effects of CTCF and OCT1 knockdown on expression were larger in ER-mutant cells. Overall, these results suggest that CTCF and OCT1 play some role in regulating mutant-specific genes, with OCT1 potentially playing a larger role.
Discussion
In breast cancer, ER mutations arise under the challenge of hormone therapy that aims to block estrogen production or disrupt ER activity through binding to the LBD. From previous studies, it is clear that mutations in ER's LBD confer ligand-independent activity to the receptor (4, 5, 10, 12, 14–16), which in large part explains the clinical observation of ER mutations in metastatic breast cancer. However, it has been reported that ER mutations cause additional gene expression changes that appear unrelated to ER's usual target genes (17–19). In this study, we aimed to confirm ER mutant's ability to regulate genes that are not normally E2 responsive. By creating isogenic models of the Y537S and D538G mutations in MCF7 and T-47D breast cancer cell lines, we found that thousands of non-E2–responsive genes changed expression because of the mutation, in addition to hundreds of E2-responsive genes that exhibited ligand-independent regulation in mutant cells. The genes that changed expression were related to increased growth and migration/invasion, phenotypes that were observed in the ER-mutant cells. The increased proliferation may be mediated through E2F transcription factors, based on expression and motif analysis. To determine whether mutant ER gene expression effects were specific to our clones or consistently observed across multiple models from different studies, we performed an aggregate gene expression analysis with a total of 16 lines for T-47D and 18 lines for MCF7 from three studies. When taking into account study-to-study variation, we identified thousands of genes that consistently changed expression because of ER mutation and did not change expression in response to 8–24 hours of E2 treatment. We observed that most genes had different responses to Y537S and D538G, confirming that allele-specific changes in gene expression were consistent across models from different studies. In addition, we found a significant overlap between the ER-mutant effects in cell lines and genes that were differentially expressed between ER-mutant and ER WT metastatic breast tumors.
One possible explanation for the effect of mutant ER on genes that are not normally E2 responsive is that constant ER activity could have different gene expression consequences than short-term E2 treatments. In an effort to understand the contribution of constant ER activity to the unique ER-mutant gene expression patterns, we treated WT cells with E2 for up to 25 days and performed RNA-seq at 5-day intervals. Prolonged exposure to E2 explained approximately half of the mutant-specific gene expression effects with Y537S being more similar to long-term E2 treatment than D538G reiterating the differences unique to each mutant allele. Our findings indicate that approximately half of the ER-mutant gene expression effects that appeared unrelated to normal ER regulation were in fact estrogen target genes that rely on long-term E2 exposure.
Another way to interpret our prolonged E2 treatment experiment is that half of ER-mutant regulated genes cannot be attributed to long-term E2 exposure and must involve novel function of the mutant receptor. We first analyzed genomic binding of mutant and WT ER taking advantage of the FLAG epitope tag that we introduced in the native ER locus. While we observed ligand-independent ER binding in mutant ER cells at the majority of ERBSs, we did not observe large-scale reprogramming of ER genomic binding, where approximately 10% of ER-bound loci were significantly altered by the mutations. These findings were in contrast to our recent studies on the D538G mutation in endometrial cancer cells that revealed more substantial changes in ER binding (46% of ER-bound sites were significantly affected; ref. 22). This indicates that the manner in which mutations influence ER could be unique to different cell types or that ER genomic binding is more robust to mutation in breast cancer cell lines. Both mutant-specific alterations in ER genomic binding (18) and lack of alterations (12) have been reported and our study supports the conclusion that mutations have a relatively small impact on ER genomic binding. While the overall picture of ER genomic binding does not dramatically change with mutation, we still found that approximately 10% of Y537S expression changes could be explained by significantly altered ER binding in ER-mutant cells. D538G changes in ER binding were minimal in both T-47D and MCF7 lines, suggesting that the Y537S mutation may have more of an impact on ER's genomic binding site selection. Looking at constant ERBS, we found a 10% enrichment near mutant-specific genes above background; however, this may be an underestimate of constant ERBS impact on mutant-specific expression due to the large number of constant ERBSs, which leads to a high likelihood of a gene having an ERBS nearby simply by chance.
Considering the modest changes in ER genomic binding, we looked for other regulatory regions and factors that could contribute to the unexpected mutant-specific gene expression patterns. By performing ATAC-seq in our WT and ER-mutant clones, we found that thousands of loci exhibited differential chromatin accessibility in ER-mutant cells. We were able to show consistency in chromatin accessibility changes by analyzing another set of ER-mutant clones from a different study (17). Some of these changes in chromatin could be a direct result of mutant ER binding; however, the majority of differentially accessible regions did not harbor ER binding. In addition, differentially accessible regions were at least twice as likely to be found near mutant-specific genes as differential ERBSs and could explain approximately 15% of mutant-specific up- and downregulated genes. These observations led us to the conclusion that other transcription factors contribute to the mutant-specific chromatin landscape and gene expression program.
Motif analysis of differential ATAC-seq sites uncovered OCT, forkhead, and CTCF motifs. Additional experiments confirmed that OCT1, FOXA1, and CTCF were bound at a significant fraction of differentially accessible regions, suggesting that they may be playing a role in increasing chromatin accessibility or possibly taking advantage of new open regions. FOXA1 is a critical factor for ER genomic binding in breast cancer (39) and previous reports have suggested that FOXA1 may not bind to mutant-specific ER-bound sites (18). However, our data reveal that FOXA1 could be playing a role in mutant-specific gene expression. OCT1 has been shown to both negatively and positively impact gene expression depending on its associating factors. It can both maintain a poised state of gene expression by facilitating the removal of repressive histone marks and drive a repressed gene expression state by failing to facilitate the removal of these marks (40). In its role as a transcription factor, OCT1 has been shown to contribute to increased cell proliferation, cellular stress response, altered metabolism, regulation of the immune response, and invasion and metastasis (41, 42). OCT1 also regulates somatic and cancer stem cell phenotypes in normal and tumor cells (43). OCT1's presence at mutant-enriched accessible chromatin and increased expression in ER-mutant cells suggest that it may be playing a regulatory role and contributing to the expression changes caused by mutant ER. CTCF was found specifically at the Y537S-enriched ATAC-seq sites. CTCF plays a key role in the 3D organization of the genome, as reviewed by Ong and colleagues (44), and the presence of CTCF at differential ATAC-seq regions indicates that 3D genome architecture could be altered in ER-mutant cells. Depletion of CTCF and OCT1 in ER-mutant cells revealed a role for these factors in the regulation of some ER mutant–specific genes. The fact that some mutant-specific genes did not change expression with OCT1 or CTCF knockdown, indicates that many factors may be contributing to the unique gene expression program observed in ER-mutant cells. Overall, our analyses of ER genomic binding and chromatin accessibility can explain approximately a quarter of mutant-specific gene expression effects and prolonged ER activity can explain an overlapping 50% of mutant-specific genes. This indicates that additional gene regulation mechanisms are contributing to the ER-mutant gene expression program, which may include changes in the 3D genome architecture, posttranscriptional changes, and differences at ERBSs, possibly through differential recruitment of cofactors.
Our findings show that the majority of gene expression patterns, as well as a small portion of ER genome binding sites are unique to either the Y537S or D538G mutations. In addition, we saw clear allele-specific changes in chromatin accessibility and accompanying transcription factors. Previous studies have also demonstrated allele-specific differences between ER mutations that appear in the structures and molecular effects of mutant ER, including differences in their ability to bind cofactors and in their response to hormone therapies (5, 6, 11, 13, 17, 18, 45). These unique differences between ER mutants have significant clinical relevance as ESR1 mutational status could be used to determine treatment strategies that would best serve the patient. This is particularly true as novel therapies are identified that more effectively treat breast cancers harboring specific ESR1 mutations. Overall, our study shows that ER mutations consistently impact the expression of thousands of genes that are not normally estrogen regulated and they do so partially through constant ER activity and partially through the use of novel regulatory elements.
Authors' Disclosures
S. Arnesen reports grants from Department of Defense during the conduct of the study. M.M. Williams reports grants from NIH during the conduct of the study. K.C. Berrett reports grants from Department of Defense during the conduct of the study. S. Oesterreich reports grants from NCI during the conduct of the study. J. Gertz reports grants from Department of Defense during the conduct of the study and Zenopharm outside the submitted work. No disclosures were reported by the other authors.
Authors' Contributions
S. Arnesen: Conceptualization, data curation, formal analysis, investigation, visualization, writing-original draft, writing-review and editing. Z. Blanchard: Investigation, writing-review and editing. M.M. Williams: Investigation, writing-review and editing. K.C. Berrett: Investigation. Z. Li: Formal analysis, writing-review and editing. S. Oesterreich: Resources, supervision, funding acquisition, project administration, writing-review and editing. J.K. Richer: Supervision, funding acquisition, project administration, writing-review and editing. J. Gertz: Conceptualization, formal analysis, supervision, funding acquisition, methodology, writing-original draft, project administration, writing-review and editing.
Acknowledgments
This work was supported by a DOD Breast Cancer Research Program Breakthrough Award (BC151357 to J. Gertz and J.K. Richer), a DOD Breast Cancer Research Program Expansion Award (BC181341 to J. Gertz), NIH R01CA221303 (to S. Oesterreich), and the Huntsman Cancer Institute. Research reported in this article utilized the High-Throughput Genomics Shared Resource at the University of Utah and was supported by NIH/NCI award P30 CA042014. We thank K.-T. Varley, D. Tantin, and D. Savic as well as Gertz, Varley, Richer, and Oesterreich laboratory members for their helpful comments on the study and the article.
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.