1. Cancer Biology
  2. Cell Biology

Heat shock factor 1 (HSF1) cooperates with estrogen receptor α (ERα) in the regulation of estrogen action in breast cancer cells

  1. Natalia Vydra  Is a corresponding author
  2. Patryk Janus
  3. Paweł Kus
  4. Tomasz Stokowy
  5. Katarzyna Mrowiec
  6. Agnieszka Toma-Jonik
  7. Aleksandra Krzywon
  8. Alexander Jorge Cortez
  9. Bartosz Wojtas
  10. Bartłomiej Gielniewski
  11. Roman Jaksik
  12. Marek Kimmel
  13. Wieslawa Widlak  Is a corresponding author
  1. Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Poland
  2. Department of Systems Biology and Engineering, Silesian University of Technology, Poland
  3. Department of Clinical Science, University of Bergen, Norway
  4. Laboratory of Molecular Neurobiology, Nencki Institute of Experimental Biology, Polish Academy of Sciences, Poland
  5. Departments of Statistics and Bioengineering, Rice University, United States
https://doi.org/10.7554/eLife.69843
Share this article
Cite this article
  1. Natalia Vydra
  2. Patryk Janus
  3. Paweł Kus
  4. Tomasz Stokowy
  5. Katarzyna Mrowiec
  6. Agnieszka Toma-Jonik
  7. Aleksandra Krzywon
  8. Alexander Jorge Cortez
  9. Bartosz Wojtas
  10. Bartłomiej Gielniewski
  11. Roman Jaksik
  12. Marek Kimmel
  13. Wieslawa Widlak
(2021)
Heat shock factor 1 (HSF1) cooperates with estrogen receptor α (ERα) in the regulation of estrogen action in breast cancer cells
eLife 10:e69843.
https://doi.org/10.7554/eLife.69843
  2. Download BibTeX
  3. Download .RIS

Abstract

Heat shock factor 1 (HSF1), a key regulator of transcriptional responses to proteotoxic stress, was linked to estrogen (E2) signaling through estrogen receptor α (ERα). We found that an HSF1 deficiency may decrease ERα level, attenuate the mitogenic action of E2, counteract E2-stimulated cell scattering, and reduce adhesion to collagens and cell motility in ER-positive breast cancer cells. The stimulatory effect of E2 on the transcriptome is largely weaker in HSF1-deficient cells, in part due to the higher basal expression of E2-dependent genes, which correlates with the enhanced binding of unliganded ERα to chromatin in such cells. HSF1 and ERα can cooperate directly in E2-stimulated regulation of transcription, and HSF1 potentiates the action of ERα through a mechanism involving chromatin reorganization. Furthermore, HSF1 deficiency may increase the sensitivity to hormonal therapy (4-hydroxytamoxifen) or CDK4/6 inhibitors (palbociclib). Analyses of data from The Cancer Genome Atlas database indicate that HSF1 increases the transcriptome disparity in ER-positive breast cancer and can enhance the genomic action of ERα. Moreover, only in ER-positive cancers an elevated HSF1 level is associated with metastatic disease.

Editor's evaluation

The authors present an interesting genomics approach to understanding the role of heat shock factor 1 (HSF1) in breast cancer cells. They show that HSF1 indirectly interacts with estrogen receptor α (ERα) by regulating the transcription of HSP90, which is essential for normal folding and function of the receptor. They also show that HSF1 and ERα tether within the genome to enhance the transcription of a subset of genes associated with disease progression. Finally, they show the relevance to the breast tumors through comparing their data to publicly available data.

https://doi.org/10.7554/eLife.69843.sa0

eLife digest

About 70% of breast cancers rely on supplies of a hormone called estrogen – which is the main hormone responsible for female physical characteristics – to grow. Breast cancer cells that are sensitive to estrogen possess proteins known as estrogen receptors and are classified as estrogen-receptor positive. When estrogen interacts with its receptor in a cancer cell, it stimulates the cell to grow and migrate to other parts of the body. Therefore, therapies that decrease the amount of estrogen the body produces, or inhibit the receptor itself, are widely used to treat patients with estrogen receptor-positive breast cancers.

When estrogen interacts with an estrogen receptor known as ERα it can also activate a protein called HSF1, which helps cells to survive under stress. In turn, HSF1 regulates several other proteins that are necessary for ERα and other estrogen receptors to work properly. Previous studies have suggested that high levels of HSF1 may worsen the outcomes for patients with estrogen receptor-positive breast cancers, but it remains unclear how HSF1 acts in breast cancer cells.

Vydra, Janus, Kuś et al. used genetics and bioinformatics approaches to study HSF1 in human breast cancer cells. The experiments revealed that breast cancer cells with lower levels of HSF1 also had lower levels of ERα and responded less well to estrogen than cells with higher levels of HSF1. Further experiments suggested that in the absence of estrogen, HSF1 helps to keep ERα inactive. However, when estrogen is present, HSF1 cooperates with ERα and enhances its activity to help cells grow and migrate. Vydra, Janus, Kuś et al. also found that cells with higher levels of HSF1 were less sensitive to two drug therapies that are commonly used to treat estrogen receptor-positive breast cancers.

These findings reveal that the effect HSF1 has on ERα activity depends on the presence of estrogen. Therefore, cancer therapies that decrease the amount of estrogen a patient produces may have a different effect on estrogen receptor-positive tumors with high HSF1 levels than tumors with low HSF1 levels.

Introduction

Breast cancer is the most common malignancy in women worldwide. Four clinically relevant molecular types are distinguished based on the expression of estrogen receptors (ERs) and HER2 (ERBB2). Among them, luminal adenocarcinomas, characterized by the expression of estrogen receptors, constitute about 70% of all breast cancer cases. There are two classical nuclear estrogen receptors, ERα and ERβ (encoded by ESR1 and ESR2 genes, respectively), and structurally different GPR30 (GPER1), which is a member of the rhodopsin-like family of the G protein-coupled and seven-transmembrane receptors. ERα expression is most common in breast cancer, and its evaluation is the basis for determining the ER status. The activity of estrogen receptors is modulated by steroid hormones, mainly estrogens, which are synthesized from cholesterol via androgens in the reaction catalyzed by aromatase (Fuentes and Silveyra, 2019). According to epidemiological and experimental data, estrogens alongside the mutations in BRCA1 and BRCA2, CHEK2, TP53, STK11 (LKB1), PIK3CA, PTEN, and other genes, are key etiological factors of breast cancer development (Yaşar et al., 2017; Verigos and Magklara, 2015). The mechanism of estrogen-stimulated breast carcinogenesis is not clear. According to the widely accepted hypothesis, estrogens acting through ERα stimulate cell proliferation and can support the growth of cells harboring mutations that then accumulate, ultimately resulting in cancer. Another hypothesis suggests the ERα-independent action of estrogens via their metabolites, which can exert genotoxic effects, contributing to cancer development (Yager and Davidson, 2006; Pescatori et al., 2021). Nevertheless, hormonal therapies targeting either estrogen production (i.e., aromatase inhibitors) or the hormone receptor itself such as selective ER modulators (SERMs; i.e., tamoxifen) and selective ER degraders (SERDs; i.e., fulvestrant) are widely used to block the mitogenic action of estrogens in patients with ER-positive breast cancer (Renoir, 2012; Farcas et al., 2021), contributing to the decline in mortality from breast cancer in recent decades (Iwase et al., 2021).

Previously, we have found that the major female sex hormone 17β-estradiol (E2) stimulates activation of heat shock factor 1 (HSF1) in estrogen-dependent breast cancer cells via MAPK signaling (Vydra et al., 2019). HSF1 is a well-known regulator of response to cellular stress induced by various environmental stimuli. It mainly regulates the expression of the heat shock proteins (HSPs), which function as molecular chaperones and regulate protein homeostasis (Ran et al., 2007). HSF1-regulated chaperones control, among others, the activity of estrogen receptors (Echeverria and Picard, 2010). ERs remain in an inactive state trapped in multimolecular chaperone complexes organized around HSP90, containing p23 (PTGES3), and immunophilins (FKBP4 or FKPB5) (Segnitz and Gehring, 1995). Upon binding to E2, ERs dissociate from the chaperone complexes and become competent to dimerize and regulate the transcription. ERs bind DNA directly to the estrogen-response elements (EREs), or indirectly, via tethering factors, and promote the transcription at either nearby promoters or through chromatin loops from distal enhancers. The dynamic action of ERs, which enables the adaptation of cancer cells and impacts the clinical outcome, relies on many transcriptional coactivators and corepressors (Heldring et al., 2007; Renoir, 2012; Farcas et al., 2021). HSP90 is essential for ERα hormone binding (Fliss et al., 2000), dimer formation (Powell et al., 2010), and binding to the EREs (Inano et al., 1994). Also, the passage of the ER to the cell membrane requires association with the HSP27 (HSPB1) oligomers in the cytoplasm (Razandi et al., 2010). More than 20 chaperones and co-chaperones associated with ERα in human cells have been identified through a quantitative proteomic approach (Dhamad et al., 2016), but their specific contribution in the receptor action still needs to be investigated. Moreover, HSF1 is involved in the regulation of a plethora of non-HSP genes, which support oncogenic processes: cell cycle regulation, signaling, metabolism, adhesion, and translation (Mendillo et al., 2012). A high level of HSF1 expression was found in cancer cell lines and many human tumors (Vydra et al., 2014; De Thonel et al., 2011) and was shown to be associated with the increased mortality of ER-positive breast cancer patients (Santagata et al., 2011; Gökmen-Polar and Badve, 2016).

E2-activated HSF1 is transcriptionally potent and takes part in the regulation of several genes essential for breast cancer cell growth (Vydra et al., 2019). Furthermore, HSF1-regulated chaperones are necessary for ERα proper function. Thus, a hypothetical positive feedback loop between E2/ERα and HSF1 signaling may exist, which putatively supports the growth of estrogen-dependent tumors. Here, to study the cooperation of HSF1 and ERα in estrogen signaling and the influence of HSF1 on E2-stimulated transcription and cell growth and mobility, we created novel experimental models based on HSF1-deficient cells and performed an in-depth bioinformatics analysis of the relevant genomics data. We also compared the influence of HSF1 on ER-positive and ER-negative breast cancers transcriptomes from The Cancer Genome Atlas (TCGA) database.

Results

HSF1 deficiency reduces the estrogen-stimulated proliferation and mobility of ERα-positive MCF7 cells

To study the contribution of HSF1 in E2 signaling, we established MCF7 cell lines with reduced HSF1 expression. Firstly, we tested a few HSF1-targeting shRNAs (Figure 1—figure supplement 1A). Then, the most potent variant that reduced HSF1 level about 10-fold (termed afterward shHSF1) was chosen for further studies. Although the heat shock response was significantly reduced, the expression of HSP genes (HSPA1A, HSPH1, HSPB1, and HSPB8) was still induced after this HSF1 knockdown (Figure 1—figure supplement 1B). Thus, we additionally created MCF7 variants with HSF1 functional knockout using the CRISPR/Cas9 gene targeting approach (clones arisen from two individual cells termed KO#1 and KO#2 afterward). Then, considering the slight differences between clones, we created an additional experimental model: six new individual HSF1-negative (HSF1−) and six HSF1-positive (HSF1+) MCF7 clones obtained using the DNA-free CRISPR/Cas9 system (which was more effective) were pooled before analyses. The complete elimination of HSF1 (Figure 1A, Figure 1—figure supplement 1A) was connected with a substantial loss of inducibility of HSP genes (Figure 1—figure supplement 1B) and proteins (HSP105/HSPH1, HSP90, HSP70/HSPA1) following hyperthermia (Figure 1B). The ability of cells to form colonies in the clonogenic assay was reduced in all MCF7 experimental models of HSF1 depletion (using shRNA and sgRNA; Figure 1C, Figure 1—figure supplement 1C). Moreover, the population size of ALDH-positive (stem/progenitor) cells correlated with the HSF1 level and was reduced in HSF1-deficient cells (Figure 1—figure supplement 1D). Also, the increased contribution of cells in the G1 phase was associated with the HSF1 knockout (Figure 1—figure supplement 1E). HSF1 knockdown did not affect the proliferation rate, while the functional HSF1 knockout led to a slight reduction in the proliferation rate under standard conditions (this effect was not visible under less favorable growing conditions, i.e., in phenol red-free 5% dextran-activated charcoal-stripped fetal bovine serum (FBS); Figure 1D, Figure 1—figure supplement 1F). To check if HSF1 deficiency would affect the growth of another ERα-positive cell line, we modified T47D cells using the CRISPR/Cas9 method (Figure 1—figure supplement 2A). Under standard conditions, we did not observe differences between HSF1+ and HSF1− T47D cells in the proliferation and clonogenic assay (not shown). Unlike MCF7 cells, HSF1− T47D cells grew slightly faster than HSF1+ cells, but this difference was not statistically significant and no differences were observed in the cell cycle (Figure 1—figure supplement 2B and C).

Figure 1 with 2 supplements see all
Download asset Open asset
Effect of HSF1 depletion on MCF7 cell growth and migration.

(A) HSF1 level and (B) heat shock response assessed by western blot in unmodified cells (WT) and in cells obtained using DNA-free CRISPR/Cas9 system: HSF1+ (six clones with the normal HSF1 level were pooled) and HSF1− (six HSF1-negative clones were pooled). Actin (ACTB) was used as a protein loading control. Heat shock: 43°C/1 hr + recovery 37°C/6 hr. (C) The number of colonies formed in the clonogenic assay: representative images of single-cell clones stained with crystal violet and their quantification (mean ± SD, n = 4). (D) Growth curves of untreated (Ctr) and E2-stimulated cells in phenol red-free media with 5% or 10% charcoal-stripped FBS (assessed using crystal violet staining). Mean and standard deviation from three independent experiments (each in three technical replicates) are shown. (E) F-actin staining in cells treated with E2 (10 nM for 14 days), then seeded for 24 hr on fibronectin-coated slides. Arrowheads, ameboid-like cells; arrows, mesenchymal-like cells; scale bar, 50 μm. (F) The number of cells after F-actin staining was counted in 10 random fields and single cells (ameboid-like and mesenchymal-like) were calculated as a percent of all cells. (G) Cell adhesion to collagens and vitronectin analyzed after E2 treatment (10 nM for 14 days); adhesion to BSA serves as a negative control (n = 4) (H) The number of migrating cells assessed by Boyden chamber assay after E2 treatment (10 nM for 14 days) (n = 3, each in three technical replicates). Boxplots represent the median, upper and lower quartiles, maximum and minimum; ***p<0.0001, **p<0.001, *p<0.05 (significance of differences versus the corresponding control – next to the curve/box/bar or between cell variants). See Figure 1—figure supplement 1 and Figure 1—figure supplement 2 for an extended characteristic of other HSF1-deficient MCF7 and T47D cell models.

We have previously demonstrated that HSF1 was activated after E2 treatment of ERα-positive cells and it was able to bind to the regulatory sequences of several target genes, which correlated with the upregulation of their transcription (Vydra et al., 2019). Since most of these genes code for proteins involved in E2 signaling, we expected that HSF1 downregulation could affect E2-dependent processes, especially cell proliferation. Therefore, we compared E2-stimulated proliferation of HSF1-proficient and HSF1-deficient MCF7 cells in all experimental models. HSF1 deficiency resulted in weaker growth stimulation by E2 (than in the corresponding control cells), but a statistically significant difference was not observed in all experimental conditions/cell variants (Figure 1D, Figure 1—figure supplement 1G). However, E2-stimulated proliferation was not significantly reduced in HSF1-deficient T47D cells (Figure 1—figure supplement 2B). These results indicate that HSF1 may influence the growth of ER-positive breast cancer cells, unstimulated and stimulated by estrogen, although the effect also depends on other factors (differences between cells, culture conditions).

We then searched for differences between modified cells in response to longer E2 treatment. We noticed that stimulation of HSF1-proficient MCF7 cells with E2 for 7–14 days resulted in cell-cell dissociation, the acquisition of an ameboid- or mesenchymal-like morphology (Figure 1E and F, Figure 1—figure supplement 1I), and enhanced adhesion to collagens (I and IV) but reduced to vitronectin (Figure 1G; adhesion to fibronectin, laminin, and tenascin was not affected; not shown). These changes enabled cells to migrate faster (Figure 1H, Figure 1—figure supplement 1H). HSF1 deficiency counteracted cell scattering after E2 stimulation (Figure 1E and F, Figure 1—figure supplement 1I ). This was associated with the reduced adhesion to collagens and cell motility (Figure 1G and H, Figure 1—figure supplement 1H). It is noteworthy that T47D cells differed from MCF7 cells in response to E2 treatment for 14 days, especially in acquired cell morphology. Amoeboid-like morphology was dominant among the scattered MCF7 cells, while mesenchymal-like morphology was dominant in T47D cells (Figure 1E and F, Figure 1—figure supplement 2D and E). Also, adhesion to collagens was not affected by E2 in T47D cells (Figure 1—figure supplement 2F). Nevertheless, E2 treatment enhanced migration of HSF1-proficient but not HSF1-deficient T47D cells (Figure 1—figure supplement 2G).

Transcriptional response to estrogen is inhibited in HSF1-deficient cells

In a search for the mechanism responsible for a distinct response to estrogen in ER-positive cells with different levels of HSF1, we analyzed global gene expression profiles by RNA-seq in all MCF7 cell variants. At control conditions (no E2 stimulation), we found relatively few genes differentially expressed in HSF1-proficient and HSF1-deficient cells that were common for different models of HSF1 downregulation. These included mainly known HSF1 targets (e.g., HSPH1, HSPE1, HSPD1, HSP90AA1) slightly repressed in HSF1-deficient cells. Analyzing the response to E2, we initially compared cell variants from different models: with the normal level of HSF1 (WT, SCR, and MIX) and HSF1-deficient cells (shHSF1, KO#1, and KO#2) (Supplementary file 1, sheet 1). We found 50 genes similarly regulated by E2 (47 upregulated and 3 downregulated) in all HSF1-proficient MCF7 cell variants (Figure 2—figure supplement 1A and B). On the other hand, only 13 genes were similarly upregulated after E2 stimulation in HSF1-deficient MCF7 cell variants (Figure 2—figure supplement 1A and C). The gene set enrichment analyses indicated that HSF1 deficiency negatively affected the processes activated by estrogen (the early and late estrogen response; Figure 2—figure supplement 1E). Moreover, though almost all genes upregulated by E2 in HSF1-proficient cells were also upregulated in HSF1-deficient cells (except NAPRT), the degree of their activation (measured as a fold change E2 versus Ctr) was usually weaker in the latter cells (Figure 2—figure supplement 1F), which indicated that the transcriptional response to estrogen was inhibited in the lack of HSF1. Interestingly, however, several E2-dependent genes revealed slightly higher basal expression (without E2 stimulation) in HSF1-deficient cells (Figure 2—figure supplement 1G), which suggested that in the absence of E2, HSF1 could be involved in the suppression of these genes.

Considering differences between KO#1 and KO#2 HSF1 knockout clones derived from individual cells (Figure 2—figure supplement 1D), we performed an additional transcriptomic analysis using a putatively more representative MCF7 cell model obtained by DNA-free CRISPR/Cas9 method (heterogeneous populations of HSF1+ and HSF1− cells) (Supplementary file 1, sheet 2). The analysis showed that 3715 genes significantly changed the expression (2336 upregulated and 1479 downregulated) in HSF1+ cells after E2 stimulation. On the other hand, only 2969 genes (1818 upregulated and 1151 downregulated) changed the expression in HSF1− cells (Figure 2A). Thus, approximately 20% of genes responding to E2 treatment in HSF1+ cells did not respond similarly in HSF1− cells. Moreover, among genes up- or downregulated in both cell variants, approximately 68% or 81%, respectively, responded less effectively (fold change E2 versus Ctr) in HSF1− cells than HSF1+ cells (Figure 2A, bottom panel). The gene set enrichment analyses revealed the slight differences in the early and late estrogen response pathways (Hallmark gene sets M5906 and M5907) but also in genes defining epithelial-mesenchymal transition (M5930). Interestingly, the expression of genes from these pathways already differentiated untreated HSF1− and HSF1+ cells. Genes encoding cell cycle-related targets of E2F transcription factors (M5925), involved in the G2/M checkpoint (M5901) as well as ECM proteoglycans (M27219) and collagen formation (M631), also discriminated HSF1− and HSF1+ cells (Figure 2—figure supplement 2A). The analysis confirmed that the transcriptional response to estrogen was inhibited in the lack of HSF1. In addition, signaling pathways related to proliferation, migration, and collagen adhesion were identified as primarily affected, which was consistent with the results of functional tests. Among E2-responding genes that were common for all MCF7 cell models, 46 were upregulated and 2 were downregulated (Figure 2B). Though a fraction of genes with higher basal expression in HSF1− cells than in HSF1+ cells (potentially repressed by HSF1) was smaller compared to other models of HSF1 deficiency (Figure 2—figure supplement 1), it remained relevant (Figure 2C).

Figure 2 with 2 supplements see all
Download asset Open asset
The deficiency of HSF1 reduces a transcriptional response to estrogen (E2) in ER-positive MCF7 cells.

(A) Overlap of genes stimulated or repressed after the E2 treatment (RNA-seq analyses) in HSF1+ and HSF1− cells (model created as described in Figure 1). The bottom panel compares the degree of response to E2 of overlapping genes. (B) Heatmap with hierarchical clustering of normalized read counts from RNA-seq (row z-score) for selected genes (identified as similarly responding in all MCF7 cell models; see Figure 2—figure supplement 1) stimulated or repressed after the E2 treatment. (C) The response to E2 stimulation presented as a mean fold change E2 versus Ctr (dots; scale on the top) as well as changes in the expression level between Ctr and E2 (arrows begin at the level of mean expression in untreated cells and end at the level of mean expression in treated cells; normalized RNA-seq read counts with a scale on the bottom). Genes are sorted according to the hierarchical clustering shown in the heatmap. Upregulation, fold change >1.0; downregulation, fold change <1.0. Statistically significant differences between untreated HSF1+ and HSF1− cells are marked: ***p<0.0001, **p<0.001, *p<0.05. (D) Nascent RNA gene expression analyses by RT-qPCR in wild-type (WT), HSF1+, and HSF1− cells. The upper panel shows E2-stimulated changes (E2 versus Ctr fold change; E2 treatment: 10 nM, 4 hr), bottom panel shows basal expression level represented as fold differences between untreated wild-type control (WT), HSF1+, and HSF1− cells. Corresponding total RNA analyses are shown in Figure 2—figure supplement 2B. ***p<0.0001, **p<0.001, *p<0.05 (significance of differences versus the corresponding control – above the bar, or between cell variants). (E) Analyses at the protein level (western blot) after 48 hr treatment with E2. HSPA8 was used as a protein loading control. The graph below shows the results of densitometric analyses (n = 3).

To validate the RNA-seq results, we selected 13 estrogen-induced genes for RT-qPCR analyses using nascent RNA (Figure 2D). In the case of nine genes, the degree of activation was substantially lower in HSF1− than in HSF1+ cells. When the basal expression in E2-untreated cells was compared, there were 12 genes expressed at a higher level in untreated HSF1− cells in comparison to HSF1+ cells (Figure 2D). Additional RT-qPCR analyses using total RNA showed that of the 15 genes tested 12 were less activated after E2 treatment in HSF1− than in HSF1+ cells. When the basal expression in E2-untreated cells was compared, six genes were expressed at a significantly higher level and one at a lower level in HSF1− than HSF1+ cells (Figure 2—figure supplement 2B). Therefore, although the response to E2 was highly variable (differences were observed between cell models), RT-qPCR-based validation generally confirmed differences between HSF1-proficient and HSF1-deficient MCF7 cells revealed by the RNA-seq. These changes at the transcriptional level might have direct functional consequences in HSF1− cells (reduced level of E2-stimulated lcnRNAs, e.g., LINC01016) but also were connected with the reduced protein level of E2-stimulated genes (HSPB8, PHLDA1, and EGR3 are shown as examples; Figure 2E).

HSF1 influences the binding of ERα to chromatin

To further study the influence of HSF1 on estrogen signaling, we analyzed ERα binding to chromatin in HSF1-proficient and HSF1-deficient MCF7 cells. We performed ChIP-seq analyses using the first functional knockout model (KO#2 and MIX cells) and validation by ChIP-qPCR using the model obtained by the DNA-free CRISPR/Cas9 system. A list of all ERα-binding sites detected by ChIP-seq in unstimulated cells and after 30 or 60 min of E2 treatment is presented in Supplementary file 2. These analyses revealed that in unstimulated cells ERα binding was more efficient (more binding sites and increased number of tags per peak) in HSF1-deficient cell variant (KO#2) than in the corresponding HSF1-proficient control (MIX cells) (Figure 3A and B) (it is worth noting that the MIX cell variant was also different from wild-type cells, indicating that the genome organization was affected by the CRISPR/Cas9 procedure itself, possibly due to off-targets). ERα target sequences in IGFBP4 or GREB1 are examples of such increased binding efficiency in unstimulated HSF1-deficient cells (Figure 3D). Estrogen treatment for 30 or 60 min resulted in enhanced ERα binding in all cell variants. However, fold enrichment (E2 versus Ctr) was lower in HSF1-deficient cells than in HSF1-proficient cells (Figure 3C). Moreover, the number of detected peaks in the E2-treated HSF1-deficient cells was only slightly higher than in unstimulated cells (Figure 3A) and enhanced ERα binding was primarily manifested in sites already existing in unstimulated cells (Figure 3C and D). We additionally searched for ERα- binding preferences in HSF1-proficient and HSF1-deficient cells. After estrogen treatment, ERβ (ESR2) and ERα (ESR1) motifs were centrally enriched in ERα-binding regions in all cell variants (Figure 3—figure supplement 1). Moreover, in untreated cells, the motif for PBX1 (not centrally enriched in peak regions), which is a pioneer factor known to bind to the chromatin before ERα recruitment (Magnani et al., 2011), was identified by MEME-ChIP analysis in all cell variants (not shown). This indicates that ERα chromatin-binding preferences were not substantially changed in HSF1-deficient cells.

Figure 3 with 2 supplements see all
Download asset Open asset
HSF1 deficiency influences the binding of ERα to chromatin in ER-positive MCF7 cells.

(A) Number of peaks and peak size distribution (number of tags per peak), (B) heatmap visualization of ERα ChIP-seq data (versus input), and (C) binding enrichment (fold enrichment E2 versus Ctr) after E2 stimulation (10 nM for 30 or 60 min) in HSF1-deficient cells (KO#2) and corresponding control (MIX, a combination of control clones arisen from single cells following CRISPR/Cas9 gene targeting). Heatmaps depict all ERα-binding events centered on the peak region within a 3 kb window around the peak. Peaks in each sample were ranked on intensity. (D) Examples of ERα peaks identified in ChIP-seq analyses, normalized by scaling factor using bamCoverage tool, and visualized by the IGV browser in unstimulated cells (Ctr) and after E2 treatment (10 nM, 30 min). The scale for each sample is shown in the left corner. Line plots show the E2/Ctr ratio obtained using the bamCoverage tool. (E) Comparison of ERα-binding efficiency (by ChIP-qPCR; % of input) in selected sequences in untreated (Ctr; upper panel) and after E2 stimulation (10 nM, 30 min; bottom panel) HSF1+ and HSF1− MCF7 cells (model created as described in Figure 1). ***p<0.0001, **p<0.001, *p<0.05 (significance of differences versus the corresponding control – above the bar, or between cell variants). (F) Western blot analysis of ERα level and its phosphorylated form (pS118) after E2 treatment (1, 10, and 100 nM for 30 or 60 min) in HSF1+ and HSF1− cells. Actin (ACTB) was used as a protein loading control. Graphs below show the results of densitometric analyses (n = 4). *p<0.05 (significance of differences versus the corresponding control – above the bar, or between cell variants, versus the same treatment).

Validation of ChIP-seq results revealed that in the case of IGFBP4 and GREB1 (i.e., sequences highly enriched with ERα after E2 stimulation) the binding efficiency of ERα was higher in unstimulated HSF1− cells than in the corresponding HSF1+ cells. On the other hand, although estrogen treatment strongly induced ERα binding, this induction was considerably lower in HSF1− cells (Figure 3E). Therefore, we confirmed that in this experimental system the deficiency of HSF1 may result in enhanced binding of unliganded ERα (in particular at strongly responsive ERα-binding sites) and weaker subsequent enrichment of ERα binding upon estrogen stimulation. However, other patterns of the response to E2 treatment are also possible, especially in sequences that were weakly enriched in ERα after stimulation, as exemplified by AMZ1, SDK2, SMPD3, and SMTNL2 (Figure 3D and E). Observed differences in response to E2 between cells with different levels of HSF1 may result from altered expression of ERα in HSF1-deficient cells. We found that although the kinetics of ERα activation (as assessed by S118 phosphorylation) in response to E2 treatment was similar in HSF1+ and HSF1− MCF7 cells, ERα and pS118 ERα levels were lower in HSF1− cells (Figure 3F).

ERα is known to be kept in an inactive state by HSP90 (Pratt and Toft, 1997), in particular by HSP90AA1 (Dhamad et al., 2016), that is, the HSF1 transcriptional target. Thus, looking for a reason for the decreased ERα level and its dysregulated binding to DNA in HSF1-deficient cells, we focused on ERα and HSP90 interactions. Analyses of the proximity of both proteins by PLA revealed that the number of ERα /HSP90 complexes decreased after estrogen treatment in HSF1+ MCF7 cells (Figure 3—figure supplement 2A). This indicates that liganded (and transcriptionally active) ERα is indeed released from the inhibitory complex with HSP90. HSP90AA1 expression was substantially reduced in HSF1-deficient cells (RNA-seq analyses), which correlated with the reduced HSP90 protein level (Figure 3—figure supplement 2B). Also, the ERα level was considerably decreased in most HSF1-deficient cell variants (except KO#1 cells; Figure 3—figure supplement 2C), especially in cells cultured in phenol-free media (Figure 3F). Therefore, we hypothesized that the number of ERα/HSP90 complexes could be reduced in HSF1-deficient cells, which would result in enhanced basal transcriptional activity of ERα in untreated cells. However, we observed an increased number of such complexes both in untreated and E2-stimulated HSF1− cells when compared to HSF1+ cells (Figure 3—figure supplement 2A). This indicates that the response to estrogen could be dysregulated in HSF1-deficient cells, also at the level of ERα/HSP90 interactions, in a mechanism not related directly to the HSP90 and ERα downregulation mediated by the HSF1 deficiency.

HSF1 can cooperate with ERα in chromatin binding and participate in the spatial organization of chromatin loops

Since estrogen-activated HSF1 was shown to bind to chromatin, we compared the binding patterns of ERα and HSF1 in wild-type MCF7 cells (using our ChIP-seq data deposited in the NCBI GEO database; accession no. GSE137558; Vydra et al., 2019). Although in untreated cells (Ctr) there were 1535 and 2248 annotated peaks for ERα and HSF1 respectively (compared to the input), only a few (below 50) binding sites with overlapped peaks for both transcription factors were identified. Moreover, these common binding regions were characterized by a small number of tags (smaller in the case of ERα) (Figure 4A; Supplementary file 3, sheet 1). On the other hand, the search for ERα and HSF1 common binding regions created after estrogen treatment (E2 versus Ctr) returned more than 200 peaks (Supplementary file 3, sheet 2). They represented only a small fraction of the total number of ERα-binding sites (~2.6% from 8320 peaks; in the case of HSF1, this represents 35% of 571 peaks) (Figure 4B). Numbers of tags per peak and fold enrichment increased after E2 stimulation for both factors, yet more for ERα than HSF1 binding in such regions (Figure 4C). These results suggest that although there is a significant overlap between two sets of peaks (p-value=0.0099, ChIPpeakAnno, peakPermTest), the cobinding of both factors in the same DNA region may not be critical in the regulation of the ERα transcriptional activity. Instead, we postulate that HSF1 may influence the organization of the chromatin loops created after estrogen stimulation. When we combined ERα and HSF1 ChIP-seq peaks with data from chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) performed by Fullwood et al., 2009, it was evident that the HSF1-binding sites mapped to ERα-interacting loci (ERα anchor regions) (Figure 4D, Figure 4—figure supplement 1) even if actual ERα binding was not detected in the same locus (examples of such anchors in FAM102A, HSPB8, PRKCE, and WWC1 regulatory sequences are shown in Figure 4D). HSF1 peaks unrelated to ERα anchoring were also existing (Figure 4—figure supplement 1B). Further analyses of the spatial organization of chromatin by chromosome conformation capture (3C) technique revealed that some interactions between different ERα anchor regions were dependent on the presence of HSF1. This is exemplified by HSPB8 and WWC1 loci analyzed in HSF1-proficient and HSF1-deficient cells (Figure 4E), which confirms the role of HSF1 in the formation of ERα-mediated chromatin loops.

Figure 4 with 2 supplements see all
Download asset Open asset
HSF1 may cooperate with ERα in DNA binding and take a part in chromatin organization.

(A) Overlapped HSF1 and ERα ChIP-seq peaks in untreated wild-type MCF7 cells – peak size distribution (number of tags per peak). (B) The number of overlapped ERα and HSF1 peaks identified after E2 stimulation in wild-type MCF7 cells. (C) Overlapped HSF1 and ERα ChIP-seq peaks in wild-type MCF7 after E2 stimulation – peak size distribution (number of tags per peak) and fold enrichment. (D) Examples of ERα and HSF1 peaks identified by MACS in ChIP-seq analyses in wild-type MCF7 cells after E2 treatment and corresponding ChIA-PET interactions (Fullwood et al., 2009) downloaded from ENCODE database and visualized by the IGV browser. The red bar shows the ERα anchor region (interacting loci), the red line – the intermediate genomic span between the two anchors forming a putative loop; the scale for each sample is shown in the left corner. (E) ERα-mediated chromatin interactions analyzed by chromosome conformation capture (3C) technique in HSPB8 and WWC1 loci. The scheme represents ERα anchor regions (red bars), HSF1-binding sites (blue arrows), and forward (F) and reverse (R) primers around subsequent HindIII cleavage sites. A model of chromatin loops resulting from interactions between ERα anchor regions is also illustrated above. Interactions between selected DNA regions were analyzed by PCR in untreated and E2-stimulated HSF1+ and HSF1− cells. (F) Interactions between ERα and HSF1 assessed by Proximity Ligation Assay (PLA) (red spots) in HSF1+ and HSF1− MCF7 cells after E2 treatment. DNA was stained with DAPI. Scale bar, 20 μm. Representative nuclei are enlarged. The number of spots per nucleus is shown in boxplots (which represent the median, upper and lower quartiles, maximum and minimum). ***p<0.0001, *p<0.05. E2, 10 nM for 60 min (or 30 min for 3C).

Though the cobinding of HSF1 and ERα to DNA was rare and relatively weak, particularly in untreated cells, the proximity of both factors was easily detected. In general, both transcription factors colocalized in the nucleus when assessed by immunofluorescence (Figure 4—figure supplement 2A). Thus, PLA spots indicating putative HSF1/ERα interactions were mainly located in the nucleus and their number increased after E2 treatment (Figure 4F, Figure 4—figure supplement 2B). However, large diversity was observed between individual cells, which suggests that also HSF1 binding to DNA may be differentiated at the single-cell level. Nevertheless, we concluded that the proximity of HSF1 and ERα putatively reflecting their interactions frequently happens in the cell nucleus.

The PLA results showed that interactions between ERα and HSF1 are possible, while the binding patterns observed in ChIP-seq combined with the ChIA-PET results suggest that different modes of these interactions are possible. To distinguish between cobinding, tethering, and canonical binding, regions of HSF1 and ERα ChIP-seq peaks (and ChIA-PET reads) were analyzed whether each sequence contained a motif match of HSF1 (heat shock element [HSE]) and/or ERα (ERE) (Figure 5A; Supplementary file 4). Analyses of all peaks/reads existing after E2 treatment showed that most of them reflected canonical binding through the corresponding motif. We found 569 genes that could be independently co-regulated by both transcription factors. In addition, ERα and HSF1 may directly cooperate in the regulation of 275 genes: cobinding was found for 65 genes and possible tethering (i.e., the presence of both transcription factors in a given chromatin region containing only one motif) for 220 genes (Figure 5B and C). In the last group, it is also possible that ERα and HSF1 bound to chromatin in different sites can interact, leading to the formation of a chromatin loop. Moreover, various binding patterns were found in the regulatory region annotated to one gene (Figure 5D). This analysis showed that ERα and HSF1 can interact more frequently through tethering (or when each is bound to a different region) than cobinding.

Figure 5
Download asset Open asset
Classification of the ERα and HSF1-binding patterns in gene regulatory regions in estrogen-treated MCF7 cells.

(A) A workflow of the search for possible ERα and HSF1-binding patterns based on the presence (+) and absence (−) of estrogen-response element (ERE) and heat shock element (HSE) motifs within the binding/anchoring sites detected by ChIP-seq/ChIA-PET. (B) Graphical illustration of possible cooperation between ERα and HSF1 in the chromatin. (C) The number of genes potentially regulated by ERα and HSF1 via canonical binding, tethering, or cobinding. Peaks were annotated to the nearest gene transcription start site (several modes of regulation are possible for one gene). (D) Comparative analysis of 569, 220, and 65 genes potentially co-regulated by canonical ERα and HSF1 binding, tethering, and cobinding, respectively; examples of genes are shown in gray boxes.

Metastatic and nonmetastatic breast cancers differ in the level of HSF1 only in the ER+ group

Our in vitro analyses indicated that HSF1 could support the transcriptional action of ERα upon estrogen treatment. On the other hand, HSF1-regulated chaperones are necessary to keep estrogen receptors in an inactive state in the absence of ligands, which collectively indicated important functional crosstalk between both factors. Therefore, to further study the significance of the interaction between ERα and HSF1 in actual breast cancer, we utilized RNA-seq data deposited in TCGA database. The analysis revealed that the transcript level of HSF1 negatively correlated with the ESR1 transcript level, although this tendency was relatively weak (Figure 6A). Neither ESR1 nor HSF1 transcript levels had a significant prognostic value (Figure 6—figure supplement 1A). Therefore, out of all breast cancer cases, we selected four groups (numbered from I to IV) characterized by significantly different levels of ERα (mRNA and protein level) and HSF1 (mRNA) expression: ER−/HSF1low, ER−/HSF1high, ER+/HSF1low, and ER+/HSF1high (Figure 6B). These groups varied in molecular subtypes composition. In ER+ cancers (luminal A, luminal B, and normal-like), the HSF1low group was more homogenous (mostly luminal A) than the HSF1high group. In ER– cases (basal-like and HER2-enriched), the HSF1high group was more homogenous (mostly basal-like) (Figure 6C). Importantly, the exclusion of cases with moderate/intermediate expression of ESR1 or HSF1 enabled us to observe the effect of both transcription factors on the survival of breast cancer patients, although the expression of HSF1 alone had no significant effect on the survival in either ER– or ER+ group analyzed separately (Figure 6—figure supplement 1B). Nevertheless, the most divergent groups were ER+/HSF1low and ER−/HSF1high (better and worse prognosis, respectively; p=0.0044), which represented luminal A and basal-like enriched groups (Figure 6D, Figure 6—figure supplement 1B). The difference between ER+/HSF1low and ER−/HSF1high cancers was also clearly visible in the multidimensional scaling (MDS) plots where the cancer cases belonging to these groups were separated. MDS plotting generally separated ER+ cases from ER−/HSF1high cases, while ER−/HSF1low cases were scattered between them (Figure 6E). On the other hand, HSF1high and HSF1low cases were not separated, although they were slightly shifted against each other. When looking at molecular subtypes, it became apparent that ER−/HER2-positive cancers were separated from ER−/basal-like cancers and slightly overlapped with ER+ cancers. These analyses indicate collectively that HSF1 and ERα may affect survival and have stronger prognostic value if analyzed together but only when extreme expression values are taken into account.

Figure 6 with 2 supplements see all
Download asset Open asset
Relationship between ERα and HSF1 expression in breast cancer.

(A) Correlation of HSF1 and ESR1 transcript level in all TCGA breast cancers. Each dot represents one cancer case; log(CPM), log2-counts per million. (B) Cases with markedly different mRNA levels of ESR1 (additionally, protein level determined by immunohistochemistry [IHC] was considered) and HSF1 selected for further analyses (groups I–IV). (C) Characteristics of selected groups by the molecular subtypes of breast cancer. (D) Kaplan–Meier plots for all selected groups. (E) Multidimensional scaling (MDS) plots of selected cases with marked: ER and HSF1 statuses (left) or molecular subtypes (right). (F) Plots of HSF1 expression levels in groups defined by ER status and presence/absence of metastases (**p<0.001). Black dots represent mean values. (G) Metastatic and nonmetastatic cases in groups of patients defined by ER status and HSF1 expression level. *p<0.05. (H) The proportion of metastatic (red) to nonmetastatic (green) cases (and their number) in ER+, ER−, and NA groups with different levels of HSF1 expression (deciles). ER+/−, estrogen receptor-positive/negative; HSF1high, high HSF1 level, HSF1low, low HSF1 level; NA, not assigned to any group.

Since in vitro analyses showed an effect of HSF1 on E2-stimulated cell migration that may facilitate metastasis formation, we checked HSF1 levels in metastatic (defined as all cases with a nonzero number of positive lymph nodes or with distant metastases; 418 cases) and nonmetastatic (399 cases) breast cancers (data deposited in TCGA database). ER+ (defined by our criteria, Figure 6B) was the only group in which HSF1 expression level was higher in metastatic cases than in nonmetastatic ones (logFC = 0.32, p-value=0.0005) (Figure 6F). When groups of patients defined by ER status were analyzed for overrepresentation of metastatic tumors, we found that they might be more common among ER+ tumors (51.3% versus 39.6% in the ER− group; p-value=0.018, Fisher’s exact test) (Figure 6G, upper panel). Furthermore, only in the ER+ group a proportional increase in metastatic disease was observed with the increase in HSF1 expression (Figure 6H) and metastatic tumors were overrepresented in ER+/HSF1high (62% versus 44.3% in ER+/HSF1low) (p-value=0.059, Pearson’s chi-squared test, verified by chi-squared posthoc test) (Figure 6G, bottom). When all patients (split into groups by ER status from TCGA clinical data and HSF1 expression split by median value or intervals) were analyzed, among-groups differences were also present (p-value = 0.001) (Figure 6—figure supplement 2). These analyses suggest that the action of HSF1 and its effect on metastasis formation may differ in ER+ versus ER− breast cancers.

HSF1 increases the disparity of the transcriptome in ER-positive breast cancers

Furthermore, we analyzed global gene expression profiles in breast cancers with different ERα and HSF1 statuses. Differential expression tests between the above-selected groups of patients (Supplementary file 5) revealed that generally ERα had a much stronger influence on the transcriptome (i.e., ER+ versus ER−) than HSF1 (i.e., HSF1high versus HSF1low). Nevertheless, differences between ER+ and ER− cases were higher in the presence of high levels of HSF1, which implicates that HSF1 increases the disparity of the transcriptome of ER+ cancers. Also, the differences in the transcript levels between HSF1high and HSF1low cancers were higher in ER+ than ER− cases (Figure 7A). Remarkably, the most divergent were ER+/HSF1low and ER−/HSF1high cancers, which resembled the most significant differences in the survival probability (Figure 6D). Then, we looked at differences in numbers of differently expressed genes (DEGs) between patients’ groups. To eliminate the possible influence of the group size on DEGs, we repeated each test 10 times, randomly subsampling groups to an equal number of cases and averaging the number of DEGs. Furthermore, to check whether heterogeneity of selected groups regarding molecular subtypes could affect observed differences in gene expression profiles, only basal-like (ER−) and luminal A (ER+) cancers were included in these tests (Figure 7B). In general, these analyses also revealed that the number of genes differentiating ER+ and ER− cases was higher in HSF1high cancers, while the number of genes differentiating HSF1high and HSF1low cases was higher in ER+ cancers. The most divergent were again ER+/HSF1low and ER−/HSF1high cases while the most similar, ER−/HSF1low and ER−/HSF1high (Figure 7C). This tendency was maintained when groups with mixed molecular subtypes composition were analyzed as well as more homogenous cancer groups (i.e., only basal-like and luminal A). Furthermore, the prognostic value of both ESR1 and HSF1 was visible in such homogenous groups (Figure 6—figure supplement 1C), which may simply reflect the prognostic difference between the basal-like and luminal A (i.e., ER-negative and ER-positive) breast cancer subtypes. Also, HSF1high cases were dominant in basal-like cases, while HSF1low were dominant in luminal A cases. Further analyses showed that the level of HSF1 did not affect the survival of ER-positive luminal A cancers but may slightly worsen the prognosis of basal-like cancers (Figure 6—figure supplement 1C).

Figure 7 with 2 supplements see all
Download asset Open asset
HSF1 increases the disparity of the transcriptome of ER-positive breast cancer.

(A) Boxplots of fold changes (log fold change [logFC] absolute values) illustrating differences in gene expression between groups characterized in Figure 6; represented are the median, upper/lower quartiles, and the highest/lowest values (excluding outliers shown as dots). (B) Composition of ER+ and ER− groups with different levels of HSF1 reduced to one molecular subtype (luminal A and basal, respectively). (C) The number of differently expressed genes (y-axis) plotted cumulatively against the false discovery rate (FDR) value of differences (x-axis). Comparisons of ER+ and ER− cancer cases as well as HSF1high and HSF1low: all cases (upper graphs; for group indexes see panel A) and cases from pre-selected cancer subtypes (lower graphs; for group indexes, see panel B). (D) Gene set enrichment analyses showing differences between ER+ and ER− breast cancers with different HSF1 levels. Terms related to hormone signaling and metabolism and response to stimulus and protein processing in comparisons between groups that were selected in Figure 6B. Blue, a fraction of downregulated genes; red, a fraction of upregulated genes. (E) Differences in the expression of the E2-regulated gene set (as identified in MCF7 cells by RNA-seq; see Figure 2 and Figure 2—figure supplement 1) between breast cancers with different levels of ESR1 and HSF1 selected from TCGA database and qualified into four groups as shown in Figure 6B. Green boxes mark all possible comparisons between the ER+/HSF1high group to other groups. The black horizontal line separates genes up- and downregulated after E2 treatment in MCF7 cells.

Differences in gene expression profiles between pairwise compared groups of cancer were further illustrated on volcano plots that additionally separated upregulated and downregulated genes (Figure 7—figure supplement 1). Then we searched for the hypothetical influence of the HSF1 status on functions of ERα-related genes in actual cancer tissue. The gene set enrichment analysis identified terms related to estrogen response among the most significant ones associated with transcripts differentiating between ER+ and ER− cancers. It is noteworthy that terms related to spliceosomal complex assembly, especially the formation of a quadruple snRNP complex, were differentiating HSF1high and HSF1low cancers (Figure 7—figure supplement 2). The more detailed analysis focused on terms related to hormone signaling and metabolism showed differences between HSF1high and HSF1low cases when ER+ and ER− cancers were compared. These analyses indicate that HSF1 may enhance estrogen signaling. On the other hand, the analysis focused on terms related to response to stimulus and protein processing (i.e., functions presumed to be dependent on HSF1 action via the HSPs expression) revealed that most of them reached the statistical significance of differences between ER+/HSF1high and ER−/HSF1high cases (Figure 7D).

We additionally compared the expression of E2-regulated genes (the set identified in MCF7 cells by RNA-seq, i.e., 46 upregulated and 2 downregulated genes; Figure 2) in selected groups of breast cancers with different levels of ESR1 and HSF1. The analysis revealed the highest upregulation of PGR and LINC01016 genes in ER+ compared to ER− cancers (regardless of HSF1 status) (Figure 7E). It is noteworthy, however, that not all genes upregulated by E2 in MCF7 cells revealed an increased expression level in ER+ compared to ER− cancers. Especially, FOXC1 and LINC00511 were expressed at a higher level in ER− cancers. Moreover, regardless of ER status, cancers with high HSF1 levels revealed a higher expression of MYBL1 than cancers with low HSF1 levels. Furthermore, expression of a few genes systematically differentiated cancers with high levels of both factors (ER+/HSF1high) compared to cancers with the low level of at least one factor (including RAPGEFL1, AMZ1, KCNF1, HSPB8 upregulated, and CYP24A1, SIM1 downregulated in ER+/HSF1high cancers), which was consistent in all relevant comparisons (marked with green boxes in Figure 7E). Nevertheless, the observed features of gene expression profiles confirmed collectively that HSF1 affects the genomic action of ERα in breast cancer.

HSF1 functional knockout results in a better response to hydroxytamoxifen and palbociclib treatments in MCF7 cells

ER-positive breast cancers are frequently treated with tamoxifen, a selective estrogen receptor modulator. More recent therapeutic options include palbociclib, a selective inhibitor of the cyclin-dependent kinases CDK4 and CDK6, approved for women with advanced metastatic cancer. Thus, we studied the influence of HSF1 on the response of MCF7 and T47D cells to these drugs. Treatment of HSF1+ cells with 4-hydroxytamoxifen (4-OHT) resulted in slightly enhanced proliferation (Figure 8A). This may be a consequence of ERα activation (estimated by its phosphorylation at S118; Figure 8B) and induction of ERα-regulated genes (not shown) and is consistent with previous reports (Ali et al., 1993). Although HSF1 functional knockout by itself had different effects in both cell lines (MCF7 growth was inhibited, while T47D growth was enhanced after HSF1 knockout; see Figure 1D, Figure 1—figure supplement 1F and G, Figure 1—figure supplement 2B), treatment with 4-OHT did not result in increased proliferation, thus it gave better results than in HSF1+ cells (Figure 8A). 4-OHT slightly inhibited E2-stimulated cell proliferation, and the differences between HSF1+ and HSF1− cells reflected differences in response to E2. T47D cells were more resistant to palbociclib than MCF7 cells. The difference between HSF1+ and HSF1− was not significant in T47D cells while in MCF7 cells inhibitory concentration 50 (IC50) of palbociclib was more than twofold lower in HSF1− than in HSF1+ cells (Figure 8C). Palbociclib also inhibited E2-stimulated cell proliferation, yet only in MCF7 cells, it was slightly more effective in the absence of HSF1 (Figure 8A). These results show only some tendencies (the statistical significance depends on the tests used) but suggest that ER-positive breast tumors with low HSF1 expression may be more sensitive to treatment with 4-OHT and palbociclib than cases with high HSF1 levels.

Figure 8
Download asset Open asset
HSF1 functional knockout affects the response of cells to 4-hydroxytamoxifen (4-OHT) and palbociclib (Palbo).

(A) Growth of HSF1+ and HSF1− MCF7 and T47D cells (assessed using crystal violet staining). Cells were treated with DMSO (Ctr), 4-OHT (100 nM), Palbo (1 µM in MCF7 cells, 10 µM in T47D cells), and E2 for 6 days. Boxplots represent the median, upper and lower quartiles, maximum and minimum of absorbance ratio from three independent experiments (each in two technical replicates); **p<0.001, *p<0.05 (significance of differences versus the corresponding control – above the box, or between cell variants/treatments). (B) Western blot analysis of ERα level and its phosphorylated form (pS118) after E2 or/and (4-OHT) treatment in HSF1+ and HSF1− MCF7 cells. Ctr: DMSO; E2: 10 nM E2; 4-OHT: 100 nM 4-OHT; 4-OHT+ E2: 100 nM 4-OHT and 10 nM E2. All cells were incubated for 2 hr, E2 was added 1 hr before harvesting the cells. Actin (ACTB) was used as a protein loading control. Graphs below show the results of densitometric analyses (n = 3). *p<0.05 (significance of differences versus corresponding Ctr – above the bar, or between cell variants, versus the same treatment). (C) Viability of HSF1+ and HSF1− MCF7 and T47D cells treated with palbociclib and assessed by MTS. IC50 plots and values were generated with the Quest Graph IC50 Calculator.

Discussion

The precise mechanisms by which estrogens stimulate the proliferation of breast cancer cells are still unclear. We found that estrogen action may be supported by HSF1, a deficiency of which in ER-positive MCF7 breast cancer cells slows down the mitogenic effect of estrogen. This may be a consequence of a reduced level of ERα and transcriptional response to estrogen in these cells. In addition, analyses of the transcriptome of breast cancers from TCGA database showed the importance of HSF1 as evidenced by higher transcriptome disparity in ER-positive cases with a high expression of HSF1 rather than with low HSF1 levels. The effect of E2 and ERα on cell migration and metastasis also is unclear and published data are inconsistent. E2 was shown to suppress the invasion of ER-positive breast cancer cells (Padilla-Rodriguez et al., 2018) or to enhance breast cancer cell motility and invasion (Sanchez et al., 2010; Zheng et al., 2011; Ho et al., 2016; Vazquez Rodriguez et al., 2017). Correspondingly, ERα silencing or inhibition (by fulvestrant, a selective estrogen receptor degrader) was shown to enhance cell migration and invasion (Bouris et al., 2015; Gao et al., 2017) or to reduce motility (Bischoff et al., 2020). We showed that longer exposure to E2 induced cell scattering and increased mobility in ER-positive breast cancer cells and HSF1 deficiency could counteract these processes. It is noteworthy that responses to E2 and the effects of HSF1 are slightly different in MCF7 and T47D cells. Ameboid-like morphology and enhanced adhesion to collagens are induced by E2 in MCF7, while mesenchymal-like morphology is induced in T47D cells. Generally, T47D cells differ from MCF7 cells in response to estrogen, partially due to a lower level of ERα in T47D (Vydra et al., 2019). Moreover, T47D cells harbor a p53 missense mutation (L194F), which causes p53 stabilization (Lim et al., 2009). The mutant p53 exhibits gain-of-function activities in mediating cell survival, and this is likely the reason for the differences between T47D and MCF7 cells. Nevertheless, the data from TCGA showing a correlation between increasing levels of HSF1 and metastatic disease in ER-positive breast cancers support the observations from the in vitro model that HSF1 may affect migration.

The mechanism of supportive action of HSF1 in ER-positive cells was already proposed, by which upon E2 treatment HSF1 is phosphorylated via ERα/MAPK signaling, gains transcriptional competence, and activates several genes essential for breast cancer cell growth and/or ERα action (Vydra et al., 2019). Here, we found that HSF1 deficiency results in a weaker response to estrogen stimulus of many estrogen-induced genes. It is noteworthy that the reduced transcriptional response to estrogen could at least partially result from the enhanced binding of unliganded ERα to chromatin and higher basal expression of ERα-regulated genes. This suggests that HSF1-dependent mechanisms may amplify ERα action upon estrogen stimulation while inhibiting it in the absence of ligands. The proper action of ERs depends on HSF1-regulated chaperones, especially HSP90. As expected, the number of HSP90/ERα complexes decreased after ligand (E2) binding in cells with normal levels of HSF1. However, although HSP90 was downregulated in HSF1-deficient cells, more HSP90/ERα complexes were found both in untreated and estrogen-stimulated cells. Hence, increased activity of ERα in HSF1-deficient cells could not be explained by the reduced sequestration of unliganded ERα by HSP90. Accordingly, additional HSF1-dependent factors may influence the formation of these complexes. Nevertheless, because it is known that HSP90 inhibitors affected the ERα level (Fliss et al., 2000; Nonclercq et al., 2004; Fiskus et al., 2007; Wong and Chen, 2009), a decreased level of ERα observed in our experimental model may be a consequence of the decreased level of HSP90. Ligand-independent genomic actions of ERα are also regulated by growth factors that activate protein-kinase cascades, leading to phosphorylation and activation of nuclear ERs at EREs (Stellato et al., 2016). The involvement of HSF1 in the repression of estrogen-dependent transcription was reported in MCF7 cells treated with neuregulin (NRG1), the ligand for the HER2 (NEU/ERBB2) receptor tyrosine kinase (Khaleque et al., 2008). Interactions of HSF1 with the corepressor metastasis-associated protein 1 (MTA1) and several additional chromatin-modulating proteins were implicated in that process. Therefore, since the lack of HSF1 can alter the cellular context, it cannot be ruled out that HSF1 influences unliganded and liganded ERα by various mechanisms that have to be further investigated. Our observation from cell culture models that silencing or knockout of HSF1 has a different effect on ERα-regulated genes in the absence or presence of estrogen implicates that the consequences in real cancer may depend on the hormonal status of the patient, which is connected with age (pre-/postmenopausal) or use of contraceptive and hormone replacement therapies.

Transcriptional activation by ERα is a multistep process modulated by coactivators and corepressors. Cofactors interact with the receptor in a ligand-dependent manner and are often part of large multiprotein complexes that control transcription by recruiting components of the basal transcription machinery, regulating chromatin structure, and/or modifying histones (Welboren et al., 2009; Kovács et al., 2020; Pescatori et al., 2021). Liganded ERα may bind directly to DNA (to ERE), and indirectly via tethering to other transcription factors such as FOS/JUN (AP1), STATs, ATF2/JUN, SP1, and NFκB (Björnström and Sjöberg, 2005; Welboren et al., 2009; Heldring et al., 2011). It was established that direct ERE binding is required for most (75%) of the estrogen-dependent gene regulation and 90% of the hormone-dependent recruitment of ERα to genomic binding sites (Stender et al., 2010). Therefore, 10% of ERα binding occurs through tethering factors. Here, we found that HSF1 can potentially be an additional factor tethering liganded ERα to DNA. ERα has been shown to function via extensive chromatin looping to bring genes together for coordinated transcriptional regulation (Fullwood et al., 2009). Since ERα anchor sites were identified also in sites bound by HSF1 but not ERα, we propose that HSF1 may be a part of this ‘looping’ machinery. Other components in the same anchoring center are also possible. According to the data from ENCODE, the HSF1-binding sites may coincide with NR2F2, JUND, FOSL2, CEBPB, GATA3, MAX, HDAC2, etc. It is consistent with the finding that transcription factor binding in human cells occurs in dense clusters (Yan et al., 2013). In general, estrogen-induced HSF1 binding was weaker than ERα binding. However, PLA analyses indicated a large heterogeneity in a cell population regarding ERα and HSF1 interactions. The final transcriptional activity of ERα is modulated by interactions with various tethering factors, including HSF1. Therefore, we hypothesize that it can be modulated differently at the single-cell level by different cofactors and chromatin remodeling factors. Thus, the response measured on the whole-cell population is heterogeneous, while stochastic when a single cell is considered.

Some premises indicate that high levels of HSF1 may be associated with resistance of estrogen-dependent breast cancers to hormonal therapies based on antiestrogens. A significant association between high HSF1 expression and increased mortality among the ER-positive breast cancer patients receiving hormonal therapy was first noticed by Santagata et al., 2011, then confirmed by Gökmen-Polar and Badve, 2016. It was proposed that overexpression of HSF1 in ER-positive breast cancers was associated with a decreased dependency on the ERα-controlled transcriptional program for cancer growth (Silveira et al., 2021). However, this conclusion was based on experiments performed without estrogen stimulation. Our in vitro studies indicate that the influence of HSF1 on ERα action depends on the presence of the estrogen and HSF1 may repress the ERα-controlled transcriptional program only in the absence of the ligand. Nevertheless, we confirmed that HSF1-deficient cells responded better to 4-hydroxytamoxifen treatment than cells with normal HSF1 levels. Also, palbociclib, an inhibitor of CDK4/6, was more effective in these cells. Enhanced resistance to hormonal therapies could be mediated by HSF1-regulated genes. HSPs themselves can be prognostic factors in breast cancer and especially oncogenic properties of HSP90AA1 correlated with aggressive clinicopathological features and resistance to the treatment (Whitesell et al., 2014; Klimczak et al., 2019). Here, we have proposed a novel mechanism of the HSF1 action in ER-positive breast cancers, which is independent of typical HSF1-regulated genes. This mechanism assumes that HSF1 influences the transcriptional response to estrogen via the reorganization of chromatin structure in estrogen-responsive genes. This mode of HSF1 action may be important in all ERα-expressing cells. For example, ERα is a critical transcription factor that regulates epithelial cell proliferation and ductal morphogenesis during postnatal mammary gland development. It is noteworthy that HSF1 has been shown to promote mammary gland morphogenesis by protecting mammary epithelial cells from apoptosis and increasing their proliferative capacity (Xi et al., 2012).

Experimental data indicate that the level of HSF1 can be used to predict response to treatment, while HSF1 targeting may improve the efficacy of breast cancer treatment and prevent the development of metastases. High HSF1 nuclear levels (estimated by immunohistochemistry in patients with invasive breast cancer at diagnosis; in situ carcinomas and stage IV cancers were excluded from the outcome analysis) were previously associated with decreased survival specifically in ER-positive breast cancer patients (Santagata et al., 2011). However, in another study performed on samples from patients with ER-positive tumors, only a weak association was found between the HSF1 protein expression and poor prognosis (Gökmen-Polar and Badve, 2016). Nevertheless, both studies showed a significant correlation between HSF1 transcript levels and the survival in ER-positive breast cancer patients. In our analysis, using data from TCGA gene expression database, we did not observe such a correlation. ESR1 and HSF1 levels may be prognostic only if analyzed groups of patients are preselected regarding the high/low expression, and may reflect the differences between luminal and basal cancers. We found that in cancers with high expression of estrogen receptor (ER+), the HSF1low group consisted mainly of luminal A cases, which are known to have the best prognosis. Although high HSF1 levels slightly reduced the survival in ER+ cancer patients, they had a greater negative outcome on survival in ER-negative patients. In that group, HSF1high cases consisted mainly of the basal-like subtype, known to have a worse prognosis.

In conclusion, HSF1 and ERα cooperate in response to estrogen stimulation. The regulation is known to be mediated by HSF1-dependent chaperones, which are important for the proper ERα action (Echeverria and Picard, 2010). Additionally, however, estrogen via ERα and MAPK activates HSF1 (Vydra et al., 2019), which together with ERα forms new chromatin loops that enhance estrogen-stimulated transcription (Figure 9). This may be affected by other factors (acting differently in individual cells). Moreover, HSF1 may be involved in the repression of unliganded ERα. Furthermore, genes activated by ERα and HSF1 play an important role in regulating the growth and spread of estrogen-dependent tumors.

Figure 9
Download asset Open asset
Model of cooperation between ERα and HSF1 in response to estrogen (E2) stimulation.

Both ERα and HSF1 are kept in an inactive state by the complexes of HSP90, p23, and immunophilins (I). The binding of E2 to ERα is connected with the release of the chaperone complex and activation of ERα, leading to the phosphorylation of MEK1/2 followed by HSF1 activation. Oligomers of active transcription factors can bind to DNA and cooperate in the regulation of the transcription either directly or through chromatin reorganization. This may be influenced by other factors (differently in individual cells).

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Cell line (Homo sapiens)MCF7American Type Culture CollectionCat#:HTB-22; RRID:VCL_0031
Cell line (H. sapiens)T47DEuropean Collection of Authenticated Cell CulturesCat#:85102201; RRID:CVCL_0553
Transfected construct (human)Edit-R Human HSF1 crRNAsDharmacon, Horizon Discovery Group CompanyCat#:CM-012109-02-0002Transfected construct (human)5′GGTGTCCGGGTCGCTCACGA
Transfected construct (human)Edit-R Human HSF1 crRNAsDharmacon, Horizon Discovery Group CompanyCat#:CM-012109-05-0002Transfected construct (human)AAAGTGGTCCACATCGAGCA
Transfected construct (human)Edit-R Human HSF1 crRNAsDharmacon, Horizon Discovery Group CompanyCat#:CM-012109-03-0002Transfected construct (human)GTGGTCCACATCGAGCAGGG
Transfected construct (human)Edit-R tracrRNADharmacon, Horizon Discovery Group CompanyCat#:U-002005Transfected construct
AntibodyAnti-HSF1 (rabbit polyclonal)Enzo, Life Sciences, Famingdale, NYCat#:ADI-SPA-901; RRID:AB_1083465WB (1:4000),IF (1:300)PLA (1:300)ChIP (4 μg/sample)
AntibodyAnti-HSF1 (E-4) (mouse monoclonal)Santa Cruz Biotechnology, Inc, Inc, Dallas, TXCat#:sc-17757; RRID:AB_627753PLA (1:200)
AntibodyAnti-ERα (estrogen receptor α) (D8H8) (rabbit monoclonal)Cell Signaling Technology, Danvers, MACat#:8644; RRID:AB_2617128WB (1:1000)PLA (1:200)
AntibodyAnti-ERα (mouse monoclonal) ERalphaDiagenode, Liège, BelgiumCat#:C15100066; RRID:AB_2716575PLA (1:200)ChIP (4 μg/sample)IF (1:200)
AntibodyAnti-phosphoERα (S118) (16J4) (mouse monoclonal)Cell Signaling Technology, Danvers, MACat#:2511; RRID:AB_331289WB (1:2000)
AntibodyAnti-ACTB (AC-15)(HRP) (mouse monoclonal)Sigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:A3854; RRID:AB_262011WB (1:25,000)
AntibodyAnti-HSP90 (rabbit polyclonal)Enzo, Life Sciences, USA Famingdale, NYCat#:ADI-SPA-836; RRID:AB_10615944WB (1:2000)PLA (1:200)IF (1:200)
AntibodyAnti-HSP70 (HSPA1) (mouse monoclonal)Enzo, Life Sciences, Famingdale, NYCat#:ADI-SPA-810; RRID:AB_10616513WB (1:2000)
AntibodyAnti-HSP105 (rabbit polyclonal)BioVision, Milpitas, CACat#:3390-100; RRID:AB_2264190WB (1:600)
AntibodyAnti-HSPB8/HSP22 (rabbit polyclonal)Cell Signaling Technology, Danvers, MACat#:3059; RRID:AB_2248643WB (1:1000)
AntibodyAnti-TDAG51 (PHLDA1) (RN-E62) (mouse monoclonal)Santa Cruz Biotechnology, Inc, Inc, Dallas, TXCat#:sc-23866; RRID:AB_628117WB (1:1000)
AntibodyAnti-EGR3 (A-7) (mouse monoclonal)Santa Cruz Biotechnology, Inc, Inc, Dallas, TXCat#:sc-390967; RRID:AB_2894831WB (1:1000)
AntibodyAnti-HSC70 (HSPA8) (B-6) (mouse monoclonal)Santa Cruz Biotechnology, Inc, Inc, Dallas, TXCat#:sc-7298; RRID:AB_627761WB (1:5000)
AntibodyAnti-mouse IgG (HRP)Millipore, Billerica, MACat#:AP124P; RRID:AB_90456WB (1:5000)
AntibodyAnti-rabbit IgG (HRP)Millipore, Billerica, MACat#:AP132P; RRID:AB_90264WB (1:2000)
AntibodyAnti-rabbit IgG (Alexa Fluor 488)Abcam, Cambridge, Great BritainCat#:ab150077; RRID:AB_2630356IF (1:200)
AntibodyAnti-mouse IgG (Alexa Fluor 594)Abcam, Cambridge, Great BritainCat#:ab150116; RRID:AB_2650601IF (1:200)
Recombinant DNA reagentpLVX-shRNA1 vectorClontech/Takara Bio USA, Inc.Cat#:632177Lentivirus construct to express a small hairpin RNA (shRNA)
Recombinant DNA reagentpLVX-shHSF1This paperpLVX-shRNA1 vector encoding shRNA specific for HSF1
Recombinant DNA reagentEdit-R hCMV-PuroR-Cas9 Expression PlasmidDharmacon, Horizon Discovery Group CompanyCat#:U-005100-120Cas9 expression vector
Sequence-based reagentqPCR primersThis paperSee Supplementary files 6–8
Sequence-based reagentshRNAThis paperSee Materials and methods
Peptide, recombinant proteineSpCas9-GFP proteinSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:ECAS9GFPPR
Commercial assay or kitDuolink In Situ Red Kit Mouse/RabbitSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:DUO92101
Commercial assay or kitALDEFLUOR KitSTEMCELL TechnologiesCat#:01700
Commercial assay or kitThe iDeal ChIP-seq Kit for Transcription FactorsDiagenodeCat#:C01010055
Commercial assay or kitDirect-Zol RNA MiniPrep KitZymo ResearchCat#:R2052
Commercial assay or kitμMacs Streptavidin KitMiltenyi Biotec, Bergisch Gladbach, GermanyCat#:130-074-101
Commercial assay or kitCellTiter 96 AQueous One Solution AssayPromega; Madison, WICat#:G3580
Commercial assay or kitSuperSignal West Pico PLUS Chemiluminescent SubstrateThermo Fisher Scientific, Waltham, MACat#:34577
Commercial assay or kitQIAseq Ultralow Input Library KitQiagen, Venlo, NetherlandsCatt#:180492
Commercial assay or kitECM Cell Adhesion Array kitSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:ECM540
Commercial assay or kitPCR Master Mix SYBR GreenA&A Biotechnology, Gdynia, PolandCat#:2008-100A
Chemical compound, drug17 beta-estradiolSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:E4389
Chemical compound, drug4-HydroxytamoxifenSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:T176
Chemical compound, drugPalbociclib, hydrochloride saltLC Laboratories, Woburn, MACat#:P-7788
Chemical compound, drug4-ThiouridineCayman Chemical, Ann Arbor, MICat#:16373-100
Chemical compound, drugPuromycinSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:P8833
Chemical compound, drugPhalloidin-TRITCSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:P1951IF (1:800)
Software, algorithmAdobe Photoshop CS6AdobeVersion 13.0.1; RRID:SCR_014199
Software, algorithmImageJNIHRRID:SCR_003070
Software, algorithmSamtoolsdoi:10.1093/bioinformatics/btp352RRID:SCR_002105
Software, algorithmR softwareR Foundation for Statistical ComputingPackage v.3.6.2; RRID:SCR_001905
Software, algorithmDESeq2doi:10.1158/0008-5472.CAN-13-1070RRID:SCR_015687
Software, algorithmNOISeqdoi:10.1093/nar/gkv711Package v.3.12; RRID:SCR_003002
Software, algorithmFastQC softwarehttps://www.bioinformatics.babraham.ac.uk/projects/fastqcRRID:SCR_014583
Software, algorithmHisat2doi:10.1038/nmeth.3317Version 2.0.5; RRID:SCR_015530
Software, algorithmFeatureCountsdoi:10.1093/bioinformatics/btt656Version 1.6.5; RRID:SCR_012919
Software, algorithmChIPpeakAnnodoi:10.1186/1471-2105-11-237Version 3.24.2; RRID:SCR_012828
Software, algorithmdeepTools2doi:10.1093/nar/gkw257Version 3.5.0; SCR_016366
Software, algorithmBowtie2doi:10.1038/nmeth.1923Version 2.2.9; SCR_016368
Software, algorithmMEME Suitedoi:10.1093/nar/gkv416Version. 5.4.1; RRID:SCR_001783
Software, algorithmMACS softwaredoi:10.1038/nprot.2012.101Version 1.4.2; RRID:SCR_013291
Software, algorithmBedtools softwaredoi:10.1093/bioinformatics/btq033RRID:SCR_006646
Software, algorithmMedCalc Statistical SoftwareMedCalc Software Ltd, Ostend, BelgiumVersion 19.2.1; RRID:SCR_015044
Software, algorithmChIPseekerBioconductor packageVersion 1.26.2; RRID:SCR_021322
Software, algorithmTCGAbiolinks packagedoi:10.1093/nar/gkv1507Version 2.14; RRID:SCR_017683
Software, algorithmedgeR packagedoi:10.1093/bioinformatics/btp616Version 3.28.1; RRID:SCR_012802
Software, algorithmStatisticaTIBCO Software IncRRID:SCR_014213
Software, algorithmMSigDBdoi:10.1073/pnas.0506580102RRID:SCR_016863
OtherDeoxyribonuclease IWorthington Biochemical CorporationCat#:LS006333
OtherRNAClean XP beadsBeckman Coulter Life Science, Indianapolis, INCat#:A63987
OtherMTSEA-biotin-XXBiotium, Fremont, CACat#:90066
OtherDAPI stainInvitrogenCat#:D13061 µg/ml
OtherDharmaFECT DuoDharmacon, Horizon Discovery Group CompanyCat#:T-2010Transfection reagent
OtherViromer CRISPRLipocalyx GmbH, Halle (Saale), GermanyCat#:VCr-01LB-01Transfection reagent
OthercOmplete Protease Inhibitor CocktailSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:4693116001
OtherPhosSTOP (phosphatase inhibitor tablets)Sigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:4906837001
OtherCollagen ISigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:804592
OtherCollagen IVSigma-Aldrich, Merck KGaA, Darmstadt, GermanyCat#:C55333
OtherFibronectinCorning, NYCat#:354008

Cell lines and treatments

Request a detailed protocol

Human ERα-positive MCF7 and T47D breast cancer cell lines were purchased from the American Type Culture Collection (ATCC, Manassas, VA) and the European Collection of Authenticated Cell Cultures (ECACC, Porton Down, UK), respectively. Cells were cultured in DMEM/F12 medium (Merck KGaA, Darmstadt, Germany) supplemented with 10% FBS (EURx, Gdansk, Poland) and routinely tested for mycoplasma contamination. For heat shock, logarithmically growing cells were placed in a water bath at a temperature of 43°C for 1 hr. The cells were allowed to recover for the indicated time in a CO2 incubator at 37°C. For estrogen treatment, cells were seeded on plates and the next day the medium was replaced into a phenol-free medium supplemented with 5% or 10% dextran-activated charcoal-stripped FBS (PAN-Biotech GmbH, Aidenbach, Germany). 17β-estradiol (E2; #E4389, Merck KGaA) was added 48 hr later to a final concentration of 10 nM unless otherwise stated for the indicated time. For longer E2 treatments, the medium was changed every two days. In the case of treatments with palbociclib (hydrochloride salt, #P-7788, LC Laboratories, Woburn, MA) and 4-hydroxytamoxifen (#T176, Merck KGaA), an equal volume of DMSO was added as vehicle control. The growth media were not replaced either before or after treatments. Working solutions were prepared fresh before each experiment in a culture medium (without antibiotics).

HSF1 down-regulation using shRNA

Request a detailed protocol

The shRNA target sequence for human HSF1 (NM_005526.4) was selected using the RNAi Target Sequence Selector (Clontech, Mountain View, CA). The target sequences were shHSF1 - 5′GCA GGT TGT TCA TAG TCA GAA-3′ (1994–2013 in NM_005526.4), shHSF1.2–5′CCT GAA GAG TGA AGA CAT A (526–544), and shHSF1.3–5′ CAG TGA CCA CTT GGA TGC TAT (1306–1326). The negative control sequence was 5′-ATG TAG ATA GGC GTT ACG ACT. Sense and antisense oligonucleotides were annealed and inserted into the pLVX-shRNA vector (Clontech) at BamHI/EcoRI site. Infectious lentiviruses were generated by transfecting DNA into HEK293T cells and virus-containing supernatant was collected. Human MCF7 cells were transduced with lentiviruses following the manufacturer’s instructions and selected using a medium supplemented with 1 μg/ml puromycin (Life Technologies/Thermo Fisher Scientific, Waltham, MA, USA).

HSF1 functional knockout using the CRISPR/Cas9 editing system

Request a detailed protocol

To remove the human HSF1 gene, Edit-R Human HSF1 (3297) crRNA, Edit-R tracrRNA, and Edit-R hCMV-PuroR-Cas9 Expression Plasmid (Dharmacon, Lafayette, CO) were introduced into MCF7 cells using DharmaFECT Duo (6 µg/ml) (Dharmacon) according to producer’s instruction. Transfected cells were enriched by puromycin (2 µg/ml) selection for 4 days. Afterward, single clones were obtained by limiting dilution on a 96-well plate. The efficiency of the HSF1 knockout was monitored by western blot. Out of 81 tested clones, 2 individual clones with the HSF1 knockout (KO#1 and KO#2) and 6 pooled control clones (MIX) were chosen for the next experiments. Among individually tested HSF1-targeting crRNAs, only two were effective (target sequences: GTGGTCCACATCGAGCAGGG and AAAGTGGTCCACATCGAGCA, both in exon 3 on the plus strand). For validation experiments, a new model was created using DNA-free system: Edit-R Human HSF1 (3297) crRNAs (GGTGTCCGGGTCGCTCACGA in exon 1 on the minus strand and AAAGTGGTCCACATCGAGCA in exon 3 on the plus strand), Edit-R tracrRNA (Dharmacon), and eSpCas9-GFP protein (#ECAS9GFPPR, Merck KGaA) were introduced into MCF7 and T47D cells using Viromer CRISPR (Lipocalyx GmbH, Halle [Saale], Germany) according to the manual provided by the producer. Single clones were obtained by limiting dilution on a 96-well plate. The efficiency of the HSF1 knockout was monitored by western blot and confirmed by sequencing (Genomed, Warszaw, Poland). Five (T47D) or six (MCF7) individual unaffected clones (HSF1+) or with the HSF1 functional knockout (HSF1−) were pooled each time before analyses.

Protein extraction and western blotting

Request a detailed protocol

Whole-cell extracts were prepared using RIPA buffer supplemented with cOmplete protease inhibitors cocktail (Roche) and phosphatase inhibitors PhosSTOP (Roche, Indianapolis, IN). Proteins (20–30 μg) were separated on 10% SDS-PAGE gels and blotted to a 0.45 μm pore nitrocellulose filter (GE Healthcare, Europe GmbH, Freiburg, Germany) using Trans Blot Turbo system (Bio-Rad, Hercules, CA) for 10 min. Primary antibodies against HSF1 (1:4000, ADI-SPA-901), HSP90 (1:2000, ADI-SPA-836), and HSP70 (1:2000, ADI-SPA-810), all from Enzo Life Sciences (Farmingdale, NY), HSP105 (1:600, #3390-100, BioVision, Milpitas, CA), ERα (1:2000, #8644), phosphoERα (S118) (1:2000, #2511), HSPB8 (1:1000, #3059), all from Cell Signaling Technology (Danvers, MA), PHLDA1 (1:1000, #sc-23866), EGR3 (1:1000, #sc-390967), HSPA8/HSC70 (1:5000, #sc-7298), all from Santa Cruz Biotechnology (Dallas, TX), and ACTB (1:25,000, #A3854, Merck KGaA) were used. The primary antibody was detected by an appropriate secondary antibody conjugated with horseradish peroxidase (Thermo Fisher Scientific) and visualized by ECL kit (Thermo Fisher Scientific) or WesternBright Sirius kits (Advansta, Menlo Park, CA). Imaging was performed on x-ray film or in a G:BOX chemiluminescence imaging system (Syngene, Frederick, MD). The experiments were repeated in triplicate, and blots were subjected to densitometric analyses using ImageJ software to calculate relative protein expression after normalization with loading controls (statistical significance of differences was calculated using t-test).

Total and nascent RNA isolation, cDNA synthesis, and RT-qPCR

Request a detailed protocol

For nascent RNA labeling, 500 μM of 4-thiouridine (Cayman Chemical, Ann Arbor, MI) was added to control and E2-treated cells for the duration of the treatment (4 hr). Next, total RNA was isolated using the Direct-Zol RNA MiniPrep Kit (Zymo Research, Irvine, CA), digested with DNase I (Worthington Biochemical Corporation, Lakewood, NJ), and cleaned with RNAClean XP beads (Beckman Coulter Life Science, Indianapolis, IN). 5 µg of total RNA from each sample were taken for nascent RNA fraction isolation using methane thiosulfonate (MTS) chemistry according to Duffy and Simon, 2016. After the biotinylation step using MTSEA-biotin-XX (Biotium, Fremont, CA), s4U-RNA was cleaned with RNAClean XP beads and isolated using μMacs Streptavidin Kit (Miltenyi Biotec, Bergisch Gladbach, Germany) as described (Garibaldi et al., 2017). Total RNA (1 μg) and nascent RNA (isolated from 5 μg of total RNA) from each sample were converted into cDNA as described (Kus-Liskiewicz et al., 2013). Quantitative PCR was performed using a Bio-Rad C1000 Touch thermocycler connected to the head CFX-96. Each reaction was performed at least in triplicates using PCR Master Mix SYBRGreen (A&A Biotechnology, Gdynia, Poland). Expression levels were normalized against GAPDH, ACTB, HNRNPK, HPRT1, if not stated otherwise. The set of delta-Cq replicates (Cq values for each sample normalized against the geometric mean of four reference genes) for control and tested samples were used for statistical tests and estimation of the p-value. Shown are median, maximum, and minimum values of a fold change versus untreated control. The primers used in these assays are described in Supplementary file 6.

Clonogenic assay

Request a detailed protocol

Cells were plated onto 6-well dishes (1 × 103 cells per well) and cultured for 14 days. Afterward, cells were washed with the phosphate-buffered solution (PBS) and fixed with methanol. Colonies were stained with 0.2% crystal violet, washed, and air-dried. Colonies were counted manually.

Proliferation test

Request a detailed protocol

Cells (2 × 104 cells per well) were seeded and cultured in 12-well plates. At the indicated time, cells were washed with PBS, fixed in cold methanol, and rinsed with distilled water. Cells were stained with 0.1% crystal violet for 30 min, rinsed with distilled water extensively, and dried. Cell-associated dye was extracted with 1 ml of 10% acetic acid. Aliquots (200 μl) were transferred to a 96-well plate and the absorbance was measured at 595 nm (Synergy2 microtiter plate reader, BioTek Instruments, Winooski, VT). Grow curves are shown as the ratio of the absorbance on days 2, 4, and 6 against day 0 and were calculated from 3 to 6 independent experiments, each in 2–3 technical replicates.

Cell viability assay

Request a detailed protocol

The effect of drug treatment on cell viability was determined colorimetrically using the CellTiter 96 AQueous One Solution Cell Proliferation Assay (Promega, Madison, WI) according to the manufacturer’s protocol. Cells seeded into 96-well plates (4 × 103 MCF7 cells and 1 × 104 T47D cells per well) were incubated with palbociclib in concentrations ranging from 0 to 100 µM in DMSO for the next 72 hr (max. DMSO concentration <0.5%). The experiment was performed 3–5 times with three replicates for each concentration of the tested compound. IC50 values were determined using the Quest Graph IC50 Calculator (AAT Bioquest, Inc, 28 September 2021, https://www.aatbio.com/tools/ic50-calculator) with the option ‘Set minimum response to zero’.

F-actin staining

Request a detailed protocol

Cells were treated with E2 for 14 days, then 5 × 104 cells per well were plated on fibronectin-coated Nunc Lab-Tek II chambered coverglass (#155383, Nalge Nunc International, Rochester, NY) and allowed to grow for 24 hr. Cells were briefly washed with PBS, fixed for 10 min with 4% paraformaldehyde (PFA) solution in PBS, washed with PBS (3 × 5 min), treated with 0.1% Triton-X100 in PBS for 5 min, and washed again in PBS (3 × 5 min). F-actin was visualized by incubation with phalloidin–tetramethylrhodamine B isothiocyanate (1:800, #P1951, Merck KGaA) while nuclei were counterstained with DAPI. Images were taken using Carl Zeiss LSM 710 confocal microscope with ZEN navigation software. The experiments were performed in triplicate. Cells were counted in 10 randomly selected fields for each replicate.

Cell adhesion tests

Request a detailed protocol

Cells were cultured for 14 days with or without E2 (10 nM). The ECM Cell Adhesion Array kit (#ECM 540, Merck KGaA) was used according to the manufacturer’s instructions. Adhesion to collagen I and IV (#804592 and #C5533, Merck KGaA) was independently validated on collagen-coated (1 µg/cm2) 96-well plates. Cells (7 × 104 cells per well) were allowed to adhere for 30 min. The wells were washed three times with PBS, fixed with cold methanol, and stained with 0.1% crystal violet. The adsorbed dye was extracted with a 10% acetic acid solution for 5 min, and measurement was performed at 595 nm on the Synergy2 microtiter plate reader (BioTek Instruments).

Transwell migration test (Boyden chamber assay)

Request a detailed protocol

Transwell chambers (with 8 μm pore size membrane, Becton Dickinson) were coated with fibronectin (10 μg/ml, Becton Dickinson). Cells (8 × 104) were suspended in a HEPES-buffered serum-free medium containing 0.1% BSA, seeded in the top of the chambers, and placed in the wells containing medium supplemented with 10% FBS. After 24 hr, the inserts were washed with PBS, fixed with cold methanol, rinsed with distilled water, and stained with 0.1% crystal violet for 30 min. The cells on the upper surface of the inserts were gently removed with a cotton swab. Cells that migrated onto the lower surface were counted under a microscope in five random fields; all experiments were performed three times (three technical repetitions each).

ALDEFLUOR assay

Request a detailed protocol

The assay was performed using a kit from STEMCELL Technologies (Vancouver, Canada, #01700) according to the protocol. Cells (6 × 105) were harvested by trypsinization and resuspended in 1 ml ALDEFLUOR Buffer. After the addition of 5 μl of BODIPY-aminoacetaldehyd e (BAAA), the substrate for aldehyde dehydrogenase (ALDH), and a brief mixing, 500 μl of the cell suspension (3 × 105) was immediately transferred to another tube supplemented with 5 μl of diethylaminobenzaldehyde (DEAB), a specific inhibitor of ALDH, and pipetted to mix evenly. Tubes were incubated at 5% CO2, 37°C for 60 min. Cells were collected by centrifugation and resuspended in ALDEFLUOR Buffer. Analyses were performed using the BD Canto III cytometer (Becton Dickinson, Franklin Lakes, NJ).

Cell cycle distribution

Request a detailed protocol

Cells (3 × 105 per well) were plated onto 6-well plates. The next day medium was replaced and cells were grown for an additional 48 hr. Afterward, cells were harvested by trypsinization, rinsed with PBS, and fixed with ice-cold 70% ethanol at −20°C overnight. Cells were collected by centrifugation, resuspended in PBS containing RNase A (100 µg/ml), and stained with 100 µg/ml propidium iodide solution. DNA content was analyzed using flow cytometry to monitor the cell cycle changes.

Immunofluorescence (IF)

Request a detailed protocol

Cells were plated onto Nunc Lab-Tek II chambered coverglass (#155383, Nalge Nunc International, Rochester, NY) and fixed for 15 min with 4% PFA solution in PBS, washed, treated with 0.1% Triton-X100 in PBS for 5 min, and washed again in PBS (3 × 5 min). IF imaging was performed using primary antibodies: anti-HSP90 (1:200; ADI-SPA-836, Enzo Life Science), anti-HSF1 (1:300; ADI-SPA-901, Enzo Life Sciences), or anti-ESR1 (1:200; C15100066, Diagenode) and secondary Alexa Fluor (488 or 594) conjugated antibodies (Abcam). Finally, the DNA was stained with DAPI. Images were taken using Carl Zeiss LSM 710 confocal microscope with ZEN navigation software.

Proximity Ligation Assay

Request a detailed protocol

To detect the ERα/HSP90 and ERα/HSF1 interactions, the Duolink In Situ Proximity Ligation Assay (PLA) (Merck KGaA) was used according to the manufacturer’s protocol. Cells were plated onto Nunc Lab-Tek II chambered coverglass (#155383, Nalge Nunc International) 1 day before the experiment. Cells were fixed for 15 min with 4% PFA solution in PBS, washed in PBS, and treated with 0.1% Triton-X100 in PBS for 5 min. After washing, slides were incubated in Blocking Solution and immunolabeled (overnight, 4°C) with primary antibodies diluted in the Duolink Antibody Diluent: rabbit anti-HSP90 (1:200; #ADI-SPA-836, Enzo Life Science) and mouse anti-ERalpha (1:200; #C15100066, Diagenode, Liège, Belgium) or mouse anti-HSF1 (1:200; #sc-17757, Santa Cruz Biotechnology) and rabbit anti-ERα (1:200; #8644, Cell Signaling Technology), as well as rabbit anti-HSF1 (1:300; ADI-SPA-901, Enzo Life Sciences) and mouse anti-ERalpha; negative controls were proceeded without one primary antibody or both. Then the secondary antibodies with attached PLA probes were used. Signals of analyzed complexes were observed using Carl Zeiss LSM 710 confocal microscope with ZEN navigation software; red fluorescence signal indicated proximity (<40 nm) of proteins recognized by both antibodies (Fredriksson et al., 2002). Z-stacks images (12 slices; 5.5 μm) were taken at ×630 magnification. From each experimental condition, spots from 10 to 15 images were identified using Photoshop (Red Channel → Select → Color Range) and counted (Picture → Analysis → Record the measurements). Next, the mean number of spots per cell (nucleus, cytoplasm) in each image was calculated. Experiments were repeated three times.

Statistical analyses in in vitro experiments

Request a detailed protocol

Outliers were determined using the Grubbs, Tukey criterion, and QQ plot. For each dataset, the normality of distribution was assessed by the Shapiro–Wilk test and, depending on data distribution homogeneity of variances, was verified by the Levene test or Brown–Forsythe test. For analysis of differences between compared groups with normal distribution, the quality of mean values was verified by the ANOVA test with a pairwise comparison done with the HSD Tukey test or Games–Howell test and Tamhane test depending on the homogeneity of variance. In the case of non-Gaussian distribution, the Kruskal–Wallis ANOVA was applied for the verification of the hypothesis on the equality of medians with Conover–Iman’s test for pairwise comparisons. p=0.05 was selected as a statistical significance threshold.

Global gene expression profiling

Request a detailed protocol

Total RNA was isolated from all MCF7 cell variants (untreated, treated with 10 nM E2 for 4 hr, conditions based on Vydra et al., 2019) using the Direct-Zol RNA MiniPrep Kit (Zymo Research) and digested with DNase I (Worthington Biochemical Corporation). For each experimental point, RNAs from three biological replicates were first tested by RT-qPCR for the efficiency of treatments. They were sequenced separately for HSF1+ and HSF1− cell variants or pooled before sequencing for WT, SCR, shHSF1, MIX, KO#1, and KO#2 cell variants. cDNA libraries were sequenced by Illumina HighSeq 1500 (run type: paired-end, read length: 2 × 76 bp). Raw RNA-seq reads were aligned to the human genome hg38 in a Bash environment using hisat2 v 2.0.5. (Kim et al., 2015) with Ensembl genes transcriptome reference. Aligned files were processed using Samtools (v. 1.13) (Li et al., 2009). Furthermore, reads aligned in the coding regions of the genome were counted using FeatureCounts (v. 1.6.5) (Liao et al., 2014). Further data analyses were carried out using the R software package (v. 3.6.2; R Foundation for Statistical Computing; http://www.r-project.org). Read counts were normalized using DESeq2 (v. 1.32.0) (Lowe et al., 2014), then normalized expression values were subjected to differential analysis using NOISeq package (v. 3.12) (Tarazona et al., 2015) (E2 versus Ctr in all cell variants separately). To find common genes between samples, lists of differentiating genes were compared and Venn diagrams were performed (package VennDiagram v. 1.6.20 from CRAN). Heatmaps of normalized read counts or log2 fold changes (E2 versus Ctr) for genes shared between samples were generated (package pheatmap v. 1.0.12 from CRAN). The hierarchical clustering of genes was based on Euclidean distance. Colors are scaled per row. Gene set enrichment analysis was performed as follows: from the count matrices, we filtered out all the genes with less than 10 reads in each of the libraries. Then, we analyzed the gene-level effects of E2 stimulation of cells with normal/decreased HSF1 levels, performing the DESeq2 test for paired samples, with pairs defined by the cell variant (HSF1+ and HSF1−, and separate analysis of three cell variants with normal-HSF1 level: WT, SRC, MIX, and three cell variants with decreased HSF1 level: shHSF1, KO#1, and KO#2). Finally, we performed the gene set enrichment analysis in the same way as for TCGA data (see below for details) – for each test, genes were ranked according to their minimum significant difference (MSD), CERNO test from tmod package was used to find enriched terms, and tmodPanelPlot function was used to visualize the results. The raw RNA-seq data were deposited in the NCBI GEO database; accession nos. GSE159802 and GSE186004.

Chromatin immunoprecipitation and ChIP-qPCR

Request a detailed protocol

The ChIP assay was performed according to the protocol from the iDeal ChIP-seq Kit for Transcription Factors (Diagenode) as described in detail in Vydra et al., 2019. For each IP reaction, 30 µg of chromatin and 4 μl of mouse anti-ERalpha monoclonal antibody (C15100066, Diagenode) was used. For negative controls, chromatin samples were processed without antibody (mock-IP). Obtained DNA fragments were used for global profiling of chromatin-binding sites or gene-specific ChIP-qPCR analysis using specific primers covering the known EREs. The set of delta-Cq replicates (difference of Cq values for each ChIP-ed sample and corresponding input DNA) for control and test sample were used for ERα-binding calculation (as a percent of input DNA) and estimation of the p-values. ERE motifs in individual peaks were identified using MAST software from the MEME Suite package (v. 5.1.1) (Bailey et al., 2015). The sequences of used primers are presented in Supplementary file 7.

Global profiling of chromatin-binding sites

Request a detailed protocol

In each experimental point, four ChIP biological replicates (each from 30 μg of input chromatin) were collected and combined in one sample before DNA sequencing. Immunoprecipitated DNA fragments and input DNA were sequenced using the Illumina HiSeq 1500 system and QIAseq Ultralow Input Library Kit (run type: single read, read length: 1 × 65 bp). Raw sequencing reads were analyzed according to standards of ChIP-seq data analysis as described below. Quality control of reads was performed with FastQC software (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and low-quality sequences (average phred <30) were filtered out. Remained reads were aligned to the reference human genome sequence (hg19) using the Bowtie2.2.9 program (Langmead and Salzberg, 2012). Individual peaks (Ab-ChIPed samples versus input DNA) and differential peaks (17β-estradiol-treated versus untreated cells) were detected using MACS software (v. 1.4.2) (Feng et al., 2012), whereas the outcome was annotated with Homer package (v. 4.11) (Heinz et al., 2010). Peak intersections and their genomic coordinates were found using Bedtools software (Quinlan and Hall, 2010). The input DNA was used as a reference because no sequences were obtained using a mock-IP probe. The locations of identified ERα-binding sites were compared to genomic coordinates of E2-induced HSF1 peaks from our previous ChIP-seq analysis (NCBI GEO database; accession no. GSE137558). We defined ERα/HSF1-binding sites as ‘common’ if at least the center of one peak was within the corresponding peak. Dot plots showing peak size distribution were generated using MedCalc Statistical Software (v. 19.2.1; MedCalc Software Ltd, Ostend, Belgium; https://www.medcalc.org; 2020). Coverage of bam files was normalized using deepTools (v. 3.5.0; bamCoverage and bamCompare functions) (Ramírez et al., 2016), with scaling factors normalizing the coverage to 1 million reads. ChIP-seq heatmaps were prepared using peakHeatmap function from ChIPseeker Bioconductor package (v. 1.26.2), with margins of 3000 nucleotides upstream and downstream from the promoter. Venn diagrams peak overlap statistics and permutation tests were generated using the ChIPpeakAnno package (v. 3.24.2) (Zhu et al., 2010). The raw ChIP-seq data were deposited in the NCBI GEO database; accession no. GSE159724.

MEME-ChIP analyses

Request a detailed protocol

The consensus DNA sequences for ERα binding were identified in silico by Motif Analysis of Large Nucleotide Datasets (MEME-ChIP, v. 5.1.1) (Bailey et al., 2015) using a 150 bp region centered on the summit point and visualized by CentriMo (Local Motif Enrichment Analysis) (Bailey and Machanick, 2012).

Classification of canonical binding, cobinding, and tethered binding of HSF1 and ERα

Request a detailed protocol

Data from ERα-related chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) (Fullwood et al., 2009; GSE18046) were processed as follows: raw paired-end tags detected in the same genome localizations were combined to individual peaks using GenomicRanges 1.42.0R package (‘reduce’ function to merge overlapping ranges), obtaining genomic coordinates of the ERα binding/anchor sites, then identified peaks were annotated with Homer package (v. 4.11) (Heinz et al., 2010). The locations of detected anchors were compared to genomic coordinates of HSF1 and ERα peaks (versus input) from ChIP-seq analysis of E2-treated MCF7 cells (GSE137558 and GSE159724, respectively). We defined ERα_ChIA-PET, HSF1_ChIP-seq, and ERα_ChIP-seq binding sites as overlapping if the center of each peak was within the two corresponding peaks. Additionally, for the sequence of all peaks from each analysis, a search for ERE and HSE motifs was performed using MAST software from the MEME Suite package (v. 5.4.1) (Bailey et al., 2015) with position weight matrix (PWM) of ERα and HSF1 from the JASPAR database (Fornes et al., 2020). Based on the presence of HSF1/ERα in genomic locus and/or motif matching, three types of binding regions, including canonical binding, cobinding, and tethered binding regions, were identified within ChIP-seq/ChIA-PET peaks. A procedure is illustrated in Figure 5A (lists of all peaks with annotation and information about the presence of motifs are presented in Supplementary file 4).

Chromosome conformation capture assay (3C)

Request a detailed protocol

The procedure was carried out according to the protocol from Deng and Blobel, 2017. In brief, 1 × 107 cells per sample were trypsinized and fixed with 1% formaldehyde in 1× PBS. Crosslinking was quenched by 0.125 M glycine and cells were lysed (10 mM Tris pH 8.0, 10 mM NaCl, 0.2% NP-40, protease inhibitors). Cell nuclei were resuspended in HindIII RE buffer (10 mM Tris pH 8.0, 50 mM NaCl, 10 mM MgCl2, 100 μg/ml BSA) and incubated sequentially with 0.3% SDS (1.5 hr) and 1.8% Triton X-100 (1.5 hr) at 37°C with rotation. Chromatin was cleaved using 450U HindIII restriction enzyme (BioLabs, Ipswich, MA) at 37°C overnight and diluted 15-fold in ligation buffer (50 mM Tris pH 7.5, 10 mM MgCl2, 10 mM DTT, 1% Triton X-100, 100 μg/ml BSA). Ligation was carried out using 4000U T4-DNA ligase (EURx) at 16°C overnight, in the presence of 1 mM ATP. All samples were de-crosslinked (65°C, overnight with mixing), RNase A and Proteinase K treated, and DNA was isolated using standard Phenol/Chloroform/Isoamyl alcohol purification method. Precipitated DNA was dissolved in 10 mM Tris pH 8.0 and used as a template in PCR analyses. The primers used are listed in Supplementary file 8.

Analysis of human patient TCGA data (performed using R v. 3.6.2)

Data retrieval

Request a detailed protocol

Clinical and RNA-seq (HTSeq counts) data from TCGA breast cancer (BRCA) project were downloaded (1102 total samples) and prepared using TCGAbiolinks package (v. 2.14) (Colaprico et al., 2016). An additional file with clinical data containing ER receptor status, ‘nationwidechildrens.org_clinical_patient_brca.txt’(Anaya, 2021), was downloaded directly from the GDC repository (https://portalgdccancer.gov). Molecular subtype classification (according to Berger et al., 2018) was retrieved through TCGAbiolinks.

Cases selection

Request a detailed protocol

Counts were log CPM normalized with the cpm function from the edgeR package (v. 3.28.1) (Robinson et al., 2010). Then we selected four groups (numbered I–IV) of patients: ER-positive/negative with high/low HSF1 expression level using the following clinical and expression (log CPM) criteria: ER-positive if er_status_by_ihc: ‘Positive,’ er_status_ihc_Percent_Positive: ‘90–99%’ and expression level of ESR1 >6, ER-negative if er_status_by_ihc: ‘Negative’ and expression level of ESR1 <6, HSF1high (low) if the expression of HSF1 was above 67 (below 33) percentile across all TCGA_BRCA cases. We also excluded cases classified to HER2-enriched and basal-like subtypes from the ER-positive group, and luminal A (LumA), luminal B (LumB), and normal-like subtypes from the ER-negative group. In reduced groups (numbered IR–IVR), only the luminal A and basal-like cases were analyzed.

Survival analysis was performed using the survfit function from the survival package (v. 3.1–8) and plotted with the ggsurvplot function from the survminer package (v. 0.4.6) (Therneau and Grambsch, 2000).

MDS plots were used to visualize differences between patients. We performed MDS with MDSplot function from edgeR and plotted the results with ggplot2 (v. 3.3.2) (Wickham, 2016).

HSF1 expression and the occurrence of metastases

Request a detailed protocol

As metastatic we considered the cases fulfilling any of the following conditions in the clinical data: (1) containing the greater-than-0 value in the number_of_lymphnodes_positive_by_he column and/or (2) with the presence of metastases confirmed in any of the following columns: metastatic_site_at_diagnosis, metastatic_site_at_diagnosis_other, or distant_metastasis_present_ind2. Differential expression of HSF1 between metastatic and nonmetastatic tumors was assessed using the glmQLFTest function from edgeR. Statistical significance of the differences shown by contingency tables was assessed using fisher.test() function for 2 × 2 contingency tables and chisq.test() for the greater tables. Correction for the multiple testing was done using the Benjamini and Hochberg method (FDR = p.adjust(p, method = "fdr")).

Differential expression analysis was done with the edgeR package (Robinson et al., 2010). Lowly expressed genes were filtered out by filterByExpr function with default parameters, resulting in 24,696 genes kept for statistical analysis. Then we performed a quasi-likelihood F test for all groups’ combinations one-versus-one and two obvious two-versus-two cases: ER+ versus ER− and HSF1-high versus HSF1-low (using mean expression levels for joined groups, with the weight of 0.5). p-Values were corrected for multiple testing using the Benjamini and Hochberg method.

Comparison of the number of genes differentially expressed between groups

Request a detailed protocol

To compare the size of differences identified in each test, we plotted the cumulative distributions of false discovery rate (FDR). Each test was repeated 10 times with the groups randomly subsampled to equal size of 30 (or 28 in case of reduced groups) to avoid the p-values being affected by the group size inequality. Results were averaged and plotted with ggplot2 (Wickham, 2016).

Gene set enrichment analysis

Request a detailed protocol

For gene set enrichment analysis, we selected Hallmark, BioCarta, Reactome, and PID gene sets from MSigDB (v7.0) (Subramanian et al., 2005) and merged it with the list of pathways downloaded from KEGG. DESeq2 was used to calculate log-fold changes with its standard error (Love et al., 2014). Then all genes were ordered according to their MSD calculated as |logFC| – 2*logFC_standard_error and tested for enrichment using the CERNO test (Zyla et al., 2019) from the tmod package (v. 0.44) (Weiner, 2016). The most significant results (effect size >0.65, p-value<0.001 at least in one comparison) and results for gene sets related to the biological processes of interest were visualized with the tmodPanelPlot function.

Data availability

Sequencing data have been deposited in GEO under accession codes GSE159802, GSE159724, and GSE186004.

The following data sets were generated
    1. Vydra N
    2. Janus P
    3. Widlak W
    4. Stokowy T
    (2020) NCBI Gene Expression Omnibus
    ID GSE159802. Heat Shock Factor 1 (HSF1) supports the ESR1 action in breast cancer (RNA-seq).
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159802
    1. Vydra N
    2. Janus P
    3. Widlak W
    4. Stokowy T
    (2020) NCBI Gene Expression Omnibus
    ID GSE159724. Heat Shock Factor 1 (HSF1) supports the ESR1 action in breast cancer (ChIP-seq).
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159724
    1. Vydra N
    2. Janus P
    3. Widlak W
    4. Stokowy T
    (2021) NCBI Gene Expression Omnibus
    ID GSE186004. Heat Shock Factor 1 (HSF1) regulates the ESR1 action in breast cancer (RNA-seq).
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE186004
The following previously published data sets were used
    1. Vydra N
    2. Widlak W
    3. Stokowy T
    (2019) NCBI Gene Expression Omnibus
    ID GSE137558. Heat Shock Factor 1 (HSF1) acquires transcriptional competence under 17b-estradiol in ERa-positive breast cancer cells [ChIP-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE137558

References

    1. Ali S
    2. Metzger D
    3. Bornert JM
    4. Chambon P
    (1993)
    Modulation of transcriptional activation by ligand-dependent phosphorylation of the human oestrogen receptor A/B region
    The EMBO Journal 12:1153–1160.
  52. Book
    1. Ran R
    2. Lu A
    3. Xu H
    4. Tang Y
    5. Sharp FR
    (2007) Heat‐Shock Protein Regulation of Protein Folding, Protein Degradation, Protein Function, and Apoptosis
    In: Lajtha A, Chan PH, editors. Handbook of Neurochemistry and Molecular Neurobiology: Acute Ischemic Injury and Repair in the Nervous System. Boston, MA: Springer US. pp. 89–107.
    https://doi.org/10.1007/978-0-387-30383-3_6

Article and author information

Author details

  1. Natalia Vydra

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing - review and editing
    Contributed equally with
    Patryk Janus and Paweł Kus
    For correspondence
    natalia.vydra@io.gliwice.pl
    Competing interests
    No competing interests declared

  2. Patryk Janus

    1. Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    2. Department of Systems Biology and Engineering, Silesian University of Technology, Gliwice, Poland
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Supervision, Validation, Visualization, Writing – original draft, Writing - review and editing
    Contributed equally with
    Natalia Vydra and Paweł Kus
    Competing interests
    No competing interests declared

  3. Paweł Kus

    Department of Systems Biology and Engineering, Silesian University of Technology, Gliwice, Poland
    Contribution
    Formal analysis, Methodology, Visualization, Writing – original draft
    Contributed equally with
    Natalia Vydra and Patryk Janus
    Competing interests
    No competing interests declared

  4. Tomasz Stokowy

    Department of Clinical Science, University of Bergen, Bergen, Norway
    Contribution
    Data curation, Formal analysis
    Competing interests
    No competing interests declared

  5. Katarzyna Mrowiec

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared

  6. Agnieszka Toma-Jonik

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Formal analysis, Investigation, Methodology, Validation
    Competing interests
    No competing interests declared

  7. Aleksandra Krzywon

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4796-5478

  8. Alexander Jorge Cortez

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1284-2638

  9. Bartosz Wojtas

    Laboratory of Molecular Neurobiology, Nencki Institute of Experimental Biology, Polish Academy of Sciences, Warsaw, Poland
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  10. Bartłomiej Gielniewski

    Laboratory of Molecular Neurobiology, Nencki Institute of Experimental Biology, Polish Academy of Sciences, Warsaw, Poland
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  11. Roman Jaksik

    Department of Systems Biology and Engineering, Silesian University of Technology, Gliwice, Poland
    Contribution
    Formal analysis, Supervision
    Competing interests
    No competing interests declared

  12. Marek Kimmel

    1. Department of Systems Biology and Engineering, Silesian University of Technology, Gliwice, Poland
    2. Departments of Statistics and Bioengineering, Rice University, Houston, United States
    Contribution
    Funding acquisition, Supervision
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8161-890X

  13. Wieslawa Widlak

    Maria Sklodowska-Curie National Research Institute of Oncology, Gliwice Branch, Wybrzeże Armii Krajowej, Gliwice, Poland
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Visualization, Writing – original draft, Writing - review and editing
    For correspondence
    wieslawa.widlak@io.gliwice.pl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8440-9414

Funding

National Science Centre, Poland (2014/13/B/NZ7/02341)

  • Natalia Vydra

National Science Centre, Poland (2015/17/B/NZ3/03760)

  • Wieslawa Widlak

National Science Centre, Poland (2018/29/B/ST7/02550)

  • Marek Kimmel

European Social Fund (POWR.03.02.00-00-I029)

  • Paweł Kuś
  • Alexander Jorge Cortez

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Mrs. Krystyna Klyszcz and Urszula Bojko for technical assistance. This research was funded by the National Science Centre, Poland; grant numbers 2014/13/B/NZ7/02341 to NV (functional in vitro studies), 2015/17/B/NZ3/03760 to WW (genomic studies), and 2018/29/B/ST7/02550 to MK (analyses of TCGA data). PK and AJC were co-financed by the European Union through the European Social Fund (grant POWR.03.02.00-00-I029). The funding sources were not involved in study design, data collection and interpretation, or the decision to submit the work for publication.

Copyright

© 2021, Vydra et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,774
    views
  • 305
    downloads
  • 16
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Natalia Vydra
  2. Patryk Janus
  3. Paweł Kus
  4. Tomasz Stokowy
  5. Katarzyna Mrowiec
  6. Agnieszka Toma-Jonik
  7. Aleksandra Krzywon
  8. Alexander Jorge Cortez
  9. Bartosz Wojtas
  10. Bartłomiej Gielniewski
  11. Roman Jaksik
  12. Marek Kimmel
  13. Wieslawa Widlak
(2021)
Heat shock factor 1 (HSF1) cooperates with estrogen receptor α (ERα) in the regulation of estrogen action in breast cancer cells
eLife 10:e69843.
https://doi.org/10.7554/eLife.69843

Share this article

https://doi.org/10.7554/eLife.69843

Categories and tags

Research organism

Further reading

    1. Cancer Biology
    2. Genetics and Genomics

    Recurrent disruption of tumour suppressor genes in cancer by somatic mutations in cleavage and polyadenylation signals

    Yaroslav Kainov, Fursham Hamid, Eugene V Makeyev
    Research Article

    The expression of eukaryotic genes relies on the precise 3'-terminal cleavage and polyadenylation of newly synthesized pre-mRNA transcripts. Defects in these processes have been associated with various diseases, including cancer. While cancer-focused sequencing studies have identified numerous driver mutations in protein-coding sequences, noncoding drivers – particularly those affecting the cis-elements required for pre-mRNA cleavage and polyadenylation – have received less attention. Here, we systematically analysed somatic mutations affecting 3'UTR polyadenylation signals in human cancers using the Pan-Cancer Analysis of Whole Genomes (PCAWG) dataset. We found a striking enrichment of cancer-specific somatic mutations that disrupt strong and evolutionarily conserved cleavage and polyadenylation signals within tumour suppressor genes. Further bioinformatics and experimental analyses conducted as a part of our study suggest that these mutations have a profound capacity to downregulate the expression of tumour suppressor genes. Thus, this work uncovers a novel class of noncoding somatic mutations with significant potential to drive cancer progression.

    1. Cancer Biology
    2. Immunology and Inflammation

    Tissue-resident natural killer cells support survival in pancreatic cancer through promotion of cDC1-CD8 T activity

    Simei Go, Constantinos Demetriou ... Eric O Neill
    Research Article

    The immunosuppressive microenvironment in pancreatic ductal adenocarcinoma (PDAC) prevents tumor control and strategies to restore anti-cancer immunity (i.e. by increasing CD8 T-cell activity) have had limited success. Here, we demonstrate how inducing localized physical damage using ionizing radiation (IR) unmasks the benefit of immunotherapy by increasing tissue-resident natural killer (trNK) cells that support CD8 T activity. Our data confirms that targeting mouse orthotopic PDAC tumors with IR together with CCR5 inhibition and PD1 blockade reduces E-cadherin positive tumor cells by recruiting a hypoactive NKG2D-ve NK population, phenotypically reminiscent of trNK cells, that supports CD8 T-cell involvement. We show an equivalent population in human single-cell RNA sequencing (scRNA-seq) PDAC cohorts that represents immunomodulatory trNK cells that could similarly support CD8 T-cell levels in a cDC1-dependent manner. Importantly, a trNK signature associates with survival in PDAC and other solid malignancies revealing a potential beneficial role for trNK in improving adaptive anti-tumor responses and supporting CCR5 inhibitor (CCR5i)/αPD1 and IR-induced damage as a novel therapeutic approach.

    1. Cancer Biology

    Increased inflammatory signature in myeloid cells of non-small cell lung cancer patients with high clonal hematopoiesis burden

    Hyungtai Sim, Hyun Jung Park ... Murim Choi
    Research Article

    Clonal hematopoiesis of indeterminate potential (CHIP) allows estimation of clonal dynamics and documentation of somatic mutations in the hematopoietic system. Recent studies utilizing large cohorts of the general population and patients have revealed significant associations of CHIP burden with age and disease status, including in cancer and chronic diseases. An increasing number of cancer patients are treated with immune checkpoint inhibitors (ICIs), but the association of ICI response in non-small cell lung cancer (NSCLC) patients with CHIP burden remains to be determined. We collected blood samples from 100 metastatic NSCLC patients before and after ICI for high-depth sequencing of the CHIP panel and 63 samples for blood single-cell RNA sequencing. Whole exome sequencing was performed in an independent replication cohort of 180 patients. The impact of CHIP status on the immunotherapy response was not significant. However, metastatic lung cancer patients showed higher CHIP prevalence (44/100 for patients vs. 5/42 for controls; p = 0.01). In addition, lung squamous cell carcinoma (LUSC) patients showed increased burden of larger clones compared to lung adenocarcinoma (LUAD) patients (8/43 for LUSC vs. 2/50 for LUAD; p = 0.04). Furthermore, single-cell RNA-seq analysis of the matched patients showed significant enrichment of inflammatory pathways mediated by NF-κB in myeloid clusters of the severe CHIP group. Our findings suggest minimal involvement of CHIP mutation and clonal dynamics during immunotherapy but a possible role of CHIP as an indicator of immunologic response in NSCLC patients.