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

Decision letter

  1. Maureen E Murphy
    Reviewing Editor; The Wistar Institute, United States
  2. Anna Akhmanova
    Senior Editor; Utrecht University, Netherlands
  3. Mary Ann Allen
    Reviewer; Howard Hughes Medical Institute, University of Colorado, United States
  4. Sean W Fanning
    Reviewer; University of Chicago, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

Thank you for sending your article entitled "Heat Shock Factor 1 (HSF1) as a new tethering factor for ESR1 supporting its action in breast cancer" for peer review at eLife. Your article is being evaluated by 3 peer reviewers, and the evaluation is being overseen by a Reviewing Editor and Maureen Murphy as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Mary Allen (Reviewer #2); Sean Fanning (Reviewer #3).

Heat Shock Factor 1 (HSF1) as a new tethering factor for ESR1 supporting its action in breast cancer.

The authors present an interesting genomics approach to understanding the role of heat shock factor-1 (HSF1) in breast cancer cells. Namely, they show that HSF-1 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 HSF-1 and ERα tether within the genome to enhance the transcription of a subset of genes associated with disease progression. They show relevance to the breast cancer patient through comparing their data to publicly available data. However there are concerns about some of the conclusions reached in the manuscript, and how some of the statistical analyses have been performed.

Major Concerns:

1. The RNA-seq comparisons were done by counting the number of genes that were different. They note that the basal level of several of the genes targeted by ESR1 is higher in HSF1 deficient cells. If that is true then my expectation is fold-change will be lower after ligand addition. Mathematically, fold change is a ratio, and increasing the denominator decreases the ratio value. It is unclear if they removed the ER target genes with higher basal levels in HSF1 deficient cells that the higher induction they note would still be present. To clarify my confusion, I would like to know if the difference in the presence of ligand due to the difference before the ligand was added or due to what happens after the ligand is added. Also, Figure 4 seems to imply that the level of ESR1 is lower in the HSF1- cells. Is that a factor in the RNA-seq response to E2?

2. The authors show that HSF1 deficient cells have a higher number of ESR1 peaks with more tags in the absence of E2. However, a spike-in was not used, so the change seen in the ChIP-seq could be a technical artifact. For example, formaldehyde crosslinking varies from plate to plate regardless of the sample type. Sequencing depth of the ChIP-seq samples was not discussed and unless the depth is comparable, graphs should be on normalized tags, not tags.

3. All figures that compare chip regions should be randomizing one of the two data sets and asking what is the expectation for overlap. For instance, they talk about how often HSF1 bound sites are also ESR1 bound. Knowing if that amount of interaction is expected by random chance is only possible if they shuffle the positions of one of the datasets and ask for overlap rate. (Or using an equivalent statical test like bedtools jaccard or reldist).

4. When comparing the HSF1 ChIP-seq the ESR1 ChIP-seq to the ChIA-pet their intersections, and overlap statics should be used to determine if the interactions in 5D are just random overlaps or are seen more than one would expect by chance. (See comment above on shuffling and bedtools tests.)

5. The relationship is not clear between their results regarding the sequencing of HSF1 deficient cells and the human tumor data. The molecular data implies HSF1 regulates ESR1 in wild-type cells. The authors should add to the discussion potential reasons why are ESR1+/HSF1low and ER-/HSFhigh are the most divergent groups.

6. Title: In contrast to the classical DNA-binding mode, where ER binds directly to EREs, tethering factors such as AP-1 family transcription factors (PMIDs: 11162939, 21964465) and Sp1 (PMID: 16651265), have been shown to recruit ER to their target genes, which lack canonical EREs. The proposed model suggests that HSF1 acts as a cooperative factor for ER, where both transcription factors bind DNA independently (Figure 8). The authors also stated that "co-binding of both factors in the same DNA region is not critical in the regulation of the ESR1 transcriptional activity" (pg 12, ln 286-287). Thus, the mechanism of HSF1 action does not meet the basic requirements of a tethering factor as suggested in the title.

7. Analysis of data: A strength of this manuscript is the use of RNA-seq to examine the effect of HSF1 loss on the estrogen-regulated transcriptional program in MCF7 cells. However, the data presentation is unnecessarily complicated (Figure 1A). Moreover, data presented as fold change E2/Ctr (Figure 2B, 2C and 2E) can be misleading and should be avoided in this case, because HSF1 KO appears to have gene-specific effects on estrogen-deprived cells.

8. The MGAT3/FOS/CYP24A1 expression profile, where HSF1 loss reduces gene expression with and without estrogen stimulation (Figure 2F), supports a role of HSF1 as an ER collaborator or cooperative factor. In contrast, the predominant expression profile, which is exhibited by GREB1/PGR/KDM4B, where HSF1 loss enhances gene expression in estrogen-deprived cells (Figure 2F), supports an inhibitory role for HSF1 under estrogen-deprived conditions, which undermines the role of HSF1 as a cooperative tethering factor for ER at these estrogen-induced genes. Together, these findings are interesting, but do not show conclusively that the efficacy or potency of E2 is reduced in HSF1-deficient cells.

9. The authors examined the effect of HSF1 KO on the ER-HSP90 interaction using PLA. The authors show conclusively that HSF1 loss enhanced the ER-HSP90 interaction (Figure 4). This finding is novel and compelling, but evidence that HSP90 contributes to the HSF1 KO phenotypes was not presented, which is a major weakness in the paper.

10. The authors used ChIP-seq to compare ER and HSF1 binding sites across the MCF7 cell genome. ER and HSF1 binding sites show < 3% overlap in estrogen-deprived cells, with about 4-fold more overlapping binding sites observed upon E2-stimulation (Figure 5A, 5B), which indicates that both factors generally bind DNA independently, with some cooperative binding where ER recruitment was increased. Comparing ChIP-seq recruitment data with ChIA-PET interaction data, suggests interaction between some distinct overlapping and non-overlapping HSF1 and ER binding sites (Figure 5D). It is not clear exactly how prevalent these long-range interactions are across the genome. Using PLA, the authors showed that E2 enhanced the ER-HSF1 interaction (Figure 5F), suggesting that ER can drive chromatin organization that enables these long interactions or chromatin looping. However, chromatin conformation capture (3C) showed that some of these long-range interactions occur in HSF1 KO cells (Figure 5E), suggesting that ER does not require HSF1 to drive chromatin looping. Can the authors clarify these findings into a cogent model?

11. In breast cancer, ER+ status is an independent and mostly favorable prognostic marker (PMIDs: 6488142, 28601929). High HSF1 expression is associated with mortality in ER+ breast cancer (PMID: 22042860, 27713164). Prognostic impact of HSF1 expression on ER- cases was not obvious (PMID: 27713164). Here, the authors analyzed the TCGA database for prognostic value of ER and HSF1 mRNA and suggest that ER- cases with high HSF1 mRNA levels have the poorest outcome (Figure 6D), suggesting an ER-independent role for HSF1 in breast cancer. However, this conclusion is probably a biased reflection of the prognostic difference between basal and luminal breast cancer subtypes (Figure 6C) (PMID: 26693050). Moreover, how this finding supports a HSF1 role as a cooperative or tethering factor for ER is not clear, and the prognostic value of HSF1 mRNA levels in ER+ breast cancer was not obvious from the data presented.

12. With regard to their conclusion that HSF1 increases the diversity of the transcriptome in ER-positive breast cancers, the rationale for this analysis is unclear, and the idea that HSF1 affects ER-target genes has already been demonstrated in the KO model using RNA-seq (Figure 2).

13. Stylistic changes are needed: Much of the text in the figures are too small to read. It was hard to tell the differences between the A, B, C, D groups and the a, b, c, d groups in figure 7.

14. In figure 5, they authors should be subsampling the "all" groups to the same number as the smaller groups to see if the distributions are the same.

15. Figures 5A and 5B are trying to show information for the HSF1 and ESR1 overlapping peaks. However, for that data to be useful we need to understand non-overlapping peak data. Also, is the figure showing tags from ESR1 ChIP-seq or HSF1 ChIP-seq?

16. Important concern: The authors should be using statistical tests to see if overlap in data sets are above random chance.

17. The authors generated HSF-1 knockouts for MCF7 and T47D breast cancer cells. However, the vast majority of their work was performed with only the T47D cells. These cell lines show different dependencies on ERα and antiestrogen differently affect receptor stability between them. More work should have been done outside of one knockdown and proliferation experiment to show that HSF-1 is important in both cell lines. It would give credence to its role in breast cancer cells beyond MCF7 cells.

18. The authors used an aldefluor assay to show that HSF-1 affects stem-progenitor populations. They did not show any additional data for this line of inquiry. If this phenotype is affected by HSF-1 they should examine EMT signatures along with cellular invasion, and migration.

19. Therapeutic approaches to targeting ERα were not discussed in this manuscript. However, the role of HSF-1 in the therapeutic response to antiestrogens would be of great importance to the field. Showing the impact of HSF-1 knockdown to ERα transcriptional activity and cellular proliferation in the presence of clinically important hormone therapies fulvestrant and 4-hydroxytamoxifen as well as CDK4/6 inhibitors like palbociclib would greatly enhance the impact of this paper.

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

Author response

Major Concerns:

1. The RNA-seq comparisons were done by counting the number of genes that were different. They note that the basal level of several of the genes targeted by ESR1 is higher in HSF1 deficient cells. If that is true then my expectation is fold-change will be lower after ligand addition. Mathematically, fold change is a ratio, and increasing the denominator decreases the ratio value. It is unclear if they removed the ER target genes with higher basal levels in HSF1 deficient cells that the higher induction they note would still be present. To clarify my confusion, I would like to know if the difference in the presence of ligand due to the difference before the ligand was added or due to what happens after the ligand is added. Also, Figure 4 seems to imply that the level of ESR1 is lower in the HSF1- cells. Is that a factor in the RNA-seq response to E2?

We are aware that the lower gene activation (fold-change) by estrogen in HSF1-deficient cells could partially result from the higher baseline expression level of some genes; data presented in old Figure 2F (revised Figure 2—figure supplement 1G). To elucidate that, we performed additional RNAseq analyses on the new HSF1-deficient cell model, which revealed that there are together 209 genes up-regulated by E2 at least two-fold in either HSF1+ or HSF1− cells. Higher basal expression in HSF1− cells (versus HSF1+) was observed in the case of 79 genes. In the case of 59 genes (~75%) smaller fold change after E2 treatment was observed in HSF1− cells. Lower basal expression in untreated HSF1− cells was observed in the case of 125 genes. In this set, fold change after E2 treatment was smaller in HSF1− cells in the case of 82 genes (~63%). Stronger activation by E2 (measured as fold change) in HSF1− cells was observed in the case of 67 out of 209 genes (32%), which included 25%/37% of genes with higher/lower basal expression in HSF1− cells. This suggests that the fraction of genes with a lower basal expression may be more responsive to E2 stimulation than the fraction with a higher basal expression. When we focused only on the E2-stimulated genes identified in all RNA-seq analyzes, i.e. 46 genes, 37 (~80%) responded less effectively in HSF1− cells. Significantly higher/lower basal expression in HSF1− cells (versus HSF1+) was observed in the case of 10/9 such genes. Hence, higher basal expression in HSF1-deficient cells was not the major cause for the weaker response to E2. Nevertheless, the new RNA-seq analyses generally confirmed observations from the previous model (differences are in details), that the lack of HSF1 generally results in the inhibition of transcriptional response to E2, yet some genes respond to E2 stronger in HSF1-deficient cells. Thus, the response is gene-specific.

We agree that the altered response to E2 in HSF1-deficient cells may be a consequence of the changed level of ERα (ESR1). Among two HSF1 KO clones, KO#1 had a higher level of ERα than KO#2 (shown in the revised Figure 3—figure supplement 2C). KO#1 cells responded to E2 better than KO#2 cells, although weaker than WT and control MIX cells (revised Figure 2—figure supplement 1A).

The lower level of ERα can be a consequence of the reduced HSP90 level, which is dependent on HSF1. This is in line with previous works that investigated the effects of HSP90 inhibitors on ERα (Fliss et al., 2000) (Nonclercq et al., 2004) (Fiskus et al., 2007) (Wong and Chen, 2009). However, it has been postulated recently that overexpression of HSF1 is associated with decreased levels of ERα. The authors suggested that high levels of HSF1 triggered the degradation of ERα through the proteasome (Silveira et al., 2021). Hence, the relationship between HSF1 and ERα levels remains unclear.

Importantly, the level of ERα is variable and depends on the culture conditions. The reduced level of ERα in HSF1-deficient MCF7 cells is better visible in phenol red-free medium (Author response image 1A) , which is recommended for all experiments with E2 (to avoid activation of ERα by phenol red). The level of ERα also decreases after longer E2 treatment in phenol red-free medium (Author response image 1B) .

Author response image 1
Download asset Open asset
ERα level (assessed by western blot) is decreased in HSF1− cells (new model), especially in media without phenol red (A) and in wild-type MCF7 cells after treatment with E2 (B).

Cells were seeded on plates and the next day the medium was replaced into a phenol-free medium supplemented with 10% dextranactivated charcoal-stripped FBS. 48 hours later cells were collected (A) or E2 was added to a final concentration of 10 nM for the indicated time (B).

2. The authors show that HSF1 deficient cells have a higher number of ESR1 peaks with more tags in the absence of E2. However, a spike-in was not used, so the change seen in the ChIP-seq could be a technical artifact. For example, formaldehyde crosslinking varies from plate to plate regardless of the sample type. Sequencing depth of the ChIP-seq samples was not discussed and unless the depth is comparable, graphs should be on normalized tags, not tags.

For experiments in which sequence depth differs between input and treatment samples, the MACS software linearly scales the total control tag count to be the same as the total ChIP tag count. All ERα ChIP-seq samples were normalized together and have a comparable number of reads. In the revised version, we additionally normalized ChIP-seq data by scaling factor using bamCoverage tool and obtained the same results. An enhanced binding of ERα in untreated HSF1-deficient cells and weaker enrichment after E2 stimulation were observed in many binding sites (examples of normalized peaks are shown in the revised Figure 3D). Additionally, the results were validated and confirmed by ChIP-qPCR (at least in regions well above the background level) on another HSF1 deficient cell model (Figure 3E). Nevertheless, since the binding of ERα in E2-untreated cells is relatively weak, we cannot say with certainty that it will not be influenced by other factors. Thus, we have toned down our conclusions:

“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.”

instead of:

“Therefore, we validated ChIP-seq results and confirmed that in strongly-responsive ESR1 binding sites deficiency of HSF1 correlated with enhanced binding of unliganded ESR1 and weaker enrichment of ESR1 binding upon estrogen stimulation.”

3. All figures that compare chip regions should be randomizing one of the two data sets and asking what is the expectation for overlap. For instance, they talk about how often HSF1 bound sites are also ESR1 bound. Knowing if that amount of interaction is expected by random chance is only possible if they shuffle the positions of one of the datasets and ask for overlap rate. (Or using an equivalent statical test like bedtools jaccard or reldist).

4. When comparing the HSF1 ChIP-seq the ESR1 ChIP-seq to the ChIA-pet their intersections, and overlap statics should be used to determine if the interactions in 5D are just random overlaps or are seen more than one would expect by chance. (See comment above on shuffling and bedtools tests.)

Points 3-4. We extended the analysis of ChIP-seq data using ChIPpeakAnno and bedtools Jaccard as suggested by the Reviewer. Using ChIPpeakAnno we prepared a Venn diagram that better shows the overlap between ERα and HSF1 binding regions after E2 treatment (revised Figure 4B). Using the peakPermTest function we performed a permutation test, which revealed that there is a significant overlap between two sets of peaks (P-value: 0.0099). The revised manuscript was supplemented with this information. Also, the final Jaccard statistic reflecting the similarity of the two sets, as well as the number of intersections indicated that the overlap of two sets of peaks is not accidental (p <0.01, n_intersections = 213, with the degree of similarity at the level of 0.0203251).

5. The relationship is not clear between their results regarding the sequencing of HSF1 deficient cells and the human tumor data. The molecular data implies HSF1 regulates ESR1 in wild-type cells. The authors should add to the discussion potential reasons why are ESR1+/HSF1low and ER-/HSFhigh are the most divergent groups.

The molecular data implies that HSF1 regulates ERα and may support cell growth and migration. Using cell culture models, we have shown that HSF1 silencing or knockout may have different effects on ERα regulated genes depending on the presence of estrogen. This implicates that the hormonal status connected with age (pre/post-menopausal), or use of contraceptive and hormone replacement therapies may determine the effect of HSF1 and ERα on breast cancer prognosis. Generally, the difference between ER+/HSF1low and ER-/HSF1high groups correlates with the differences between breast cancer subtypes, especially when groups reduced to basal-like and luminal A cancers were analyzed. Basal-like (ER-negative) cases are known to have a worse prognosis, while luminal A (ER-positive) cases have the best prognosis. It is noteworthy, however, that HSF1high cases dominated in basal-like cases, while HSF1low dominated in luminal A cases. Thus, the level of HSF1 may be a part of the molecular signature of these subtypes. Extended analysis of the survival probability of TCGA breast cancer patients indicated that the level of HSF1 did not affect the survival of ER-positive luminal A cancers and worsened the prognosis of basal-like cancers. However, these groups are relatively small, and no statistical significance is observed (revised Figure 6—figure supplement 1C). The reason why HSF1 is connected with a worse prognosis in ER-negative breast cancers is speculative and requires further investigation. Gene ontology analysis revealed differences in terms related to spliceosomal complex assembly, especially the formation of a quadruple snRNP complex, between HSF1high and HSF1low cancers (revised Figure 7—figure supplement 2). Thus, addressing the hypothetical role of HSF1 in alternative splicing could open a new direction in this area. These points were added to the revised results and discussion.

6. Title: In contrast to the classical DNA-binding mode, where ER binds directly to EREs, tethering factors such as AP-1 family transcription factors (PMIDs: 11162939, 21964465) and Sp1 (PMID: 16651265), have been shown to recruit ER to their target genes, which lack canonical EREs. The proposed model suggests that HSF1 acts as a cooperative factor for ER, where both transcription factors bind DNA independently (Figure 8). The authors also stated that "co-binding of both factors in the same DNA region is not critical in the regulation of the ESR1 transcriptional activity" (pg 12, ln 286-287). Thus, the mechanism of HSF1 action does not meet the basic requirements of a tethering factor as suggested in the title.

It was established that direct ERE binding is required for most of (75%) estrogen-dependent gene regulation and 90% of hormone-dependent recruitment of ERα to genomic binding sites (Stender et al., 2010). Therefore, only 10% of hormone-dependent ERα binding putatively happens through different tethering factors. Taking into account that for example FOS/JUN (AP1), STATs, ATF2/JUN, SP1, NFκB can also be found in such tethering sites, 2.5% of all ERα binding sites shared with HSF1 seems to be rational. However, to distinguish between canonical binding, cobinding, and tethered binding regions of HSF1 and ERα within ChIP-seq peaks (and ChIAPet reads), we analyzed whether the peak contained a motif match of HSF1 (HSE, heat shock element) and/or ERα (ERE, estrogenresponsive element). We applied MAST software from the MEME Suite package to scan DNA sequences corresponding to ChIP-seq peaks with the PWM (position weight matrix) of HSF1/ERα to determine the direct binding positions of the transcription factor within ChIP-seq peaks. Canonical binding was defined when the transcription factor was bound directly to DNA within its ChIP-seq peaks and its binding motif was present in the peak. Cobinding was defined when both factors were simultaneously bound to DNA within the region and both motifs were present. If ERα bound indirectly to DNA and its ChIP-seq peaks contained only HSF1 binding motifs (and vice versa), tethered binding could be expected. However, in this situation it is also very likely that both factors interact with each other when bound to DNA in different regions, leading to looping of the chromatin. Furthermore, such interactions can be mediated by other co-factors. Given that further research is needed to prove tethering, in the revised version of the manuscript we only suggest that it is possible. Thus, the title of the manuscript was changed.

We have already proposed in the first version of the manuscript that HSF1 may only be a part of the “looping” machinery and other components in the same anchoring center are also possible. Currently, new information was added to the revised discussion:

“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).”

7. Analysis of data: A strength of this manuscript is the use of RNA-seq to examine the effect of HSF1 loss on the estrogen-regulated transcriptional program in MCF7 cells. However, the data presentation is unnecessarily complicated (Figure 1A). Moreover, data presented as fold change E2/Ctr (Figure 2B, 2C and 2E) can be misleading and should be avoided in this case, because HSF1 KO appears to have gene-specific effects on estrogen-deprived cells.

We assume that the remark concerns old Figure 2A (Figure 1A was a western blot). Indeed, this figure mainly showed differences between models and could be confusing. Figures2B, C, and F showed normalized read counts (row z-score in B, C), Figure 2E showed fold changes, while Figure 2F was intended to show the normalized expression levels before and after E2 treatment. In the revised manuscript, we performed additional RNA-seq analyses of cells that were initially used for validation yet finally appeared as the most relevant/representative model. Therefore, for the simplicity of the RNA-seq data presentation, we focused in the main manuscript on the results from this new experimental model. These new data are presented in the revised Figure 2, while data from other models are presented as supplementary data (revised Figure 2—figure supplement 1). All experimental models are included to illustrate the common and specific effects for each one (each generates different side effects). In the revised Figure 2A, we showed differences in the response to estrogen treatment in HSF1+ and HSF1− cells (the Venn diagram). Then, we selected only these genes, which were responding to E2 in all our experimental models similarly, and showed the changes in the expression on the heatmap with hierarchical clustering and combined diagram with fold change and normalized expression levels (other panels of the revised Figure 2). Hence, the presentation of data was simplified in the revised manuscript as much as possible.

8. The MGAT3/FOS/CYP24A1 expression profile, where HSF1 loss reduces gene expression with and without estrogen stimulation (Figure 2F), supports a role of HSF1 as an ER collaborator or cooperative factor. In contrast, the predominant expression profile, which is exhibited by GREB1/PGR/KDM4B, where HSF1 loss enhances gene expression in estrogen-deprived cells (Figure 2F), supports an inhibitory role for HSF1 under estrogen-deprived conditions, which undermines the role of HSF1 as a cooperative tethering factor for ER at these estrogen-induced genes. Together, these findings are interesting, but do not show conclusively that the efficacy or potency of E2 is reduced in HSF1-deficient cells.

We propose that the regulation of ERα pathway by HSF1 takes place at two levels: through HSF1dependent chaperones, which are important for the proper ERα action, and at the level of chromatin organization, which may be influenced by other factors (differently in individual cells), thus the final result (effect on transcription of a specific gene in an individual cell) could be difficult to predict. Analyzes of the presence of HSE and ERE motifs in HSF1 and ERα ChIP-seq peaks and ChIA-PET reads showed that both cobinding and tethering (or interactions of both factors directly binding to DNA in distinct regions) are possible. However, other gene-specific cofactors may interact with both transcription factors in unstimulated and estrogen-stimulated cells and they may be modified differently. In the revised manuscript, we additionally compared the changes after E2 treatment in HSF1+ and HSF1− MCF7 cells created using the DNA-free CRISPR/Cas9 method. Analyzes showed again that transcriptional response to E2 was lower in HSF1− cells. This can not be simply explained by the lower level of ERα in HSF1− cells and approximately 27% of genes similarly regulated in both cell variants responded stronger (fold change) in HSF1− cells. Importantly, the response is genespecific and different modes of regulation are observed.

9. The authors examined the effect of HSF1 KO on the ER-HSP90 interaction using PLA. The authors show conclusively that HSF1 loss enhanced the ER-HSP90 interaction (Figure 4). This finding is novel and compelling, but evidence that HSP90 contributes to the HSF1 KO phenotypes was not presented, which is a major weakness in the paper.

We decided to show the above mentioned results, although they do not fit the simplest model that could be expected in the case of HSF1-deficient cells. Hence, we propose that the activity of ERα in these cells may also be influenced by other HSF1-dependent factors (yet unidentified), different then HSP90. This topic requires further research, which, in our opinion, was beyond the scope of this manuscript. HSP90 comprises 1–2% of the total cellular protein in unstressed mammalian cells and several hundred of HSP90 client substrates are estimated (Sanchez, 2012). Therefore, inhibition of HSP90 chaperone function results in the repression of cell growth and apoptosis in cancer cells, which is mediated by different HSP90 client proteins. Also, the ERα level was affected by HSP90 inhibitors (Fliss et al., 2000); (Nonclercq et al., 2004); (Fiskus et al., 2007); (Wong and Chen, 2009). Further experiments that could demonstrate how HSP90 contributes to HSF1-deficient phenotypes should be carefully planned to distinguish between effects mediated by ERα and other client proteins. Therefore, at the moment we cannot propose any meaningful experiments that could be completed within the time given for the revision. Hence, we decided to move Figure 4 to the supplement (Figure 3—figure supplement 2; also to provide a space for additional new Figures).

10. The authors used ChIP-seq to compare ER and HSF1 binding sites across the MCF7 cell genome. ER and HSF1 binding sites show < 3% overlap in estrogen-deprived cells, with about 4-fold more overlapping binding sites observed upon E2-stimulation (Figure 5A, 5B), which indicates that both factors generally bind DNA independently, with some cooperative binding where ER recruitment was increased. Comparing ChIP-seq recruitment data with ChIA-PET interaction data, suggests interaction between some distinct overlapping and non-overlapping HSF1 and ER binding sites (Figure 5D). It is not clear exactly how prevalent these long-range interactions are across the genome. Using PLA, the authors showed that E2 enhanced the ER-HSF1 interaction (Figure 5F), suggesting that ER can drive chromatin organization that enables these long interactions or chromatin looping. However, chromatin conformation capture (3C) showed that some of these long-range interactions occur in HSF1 KO cells (Figure 5E), suggesting that ER does not require HSF1 to drive chromatin looping. Can the authors clarify these findings into a cogent model?

In our opinion, HSF1 modulates the estrogenic response and the action of ERα, but it is not a critical factor (see the response to point 6 for detail). HSF1 works on two levels: indirectly, through the HSF1-regulated genes (e.g. HSP90), and directly through the influence on the organization of chromatin, which is mainly emphasized in our work. Some long-range interactions occur in HSF1deficient cells since they can be mediated by other proteins (3C is showing only connections between genomic loci regardless of which proteins are involved in them). 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). To clarify these findings we initially planned to conduct additional experiments – ChIP-3C using anti-HSF1 and anti-ERα Abs in MCF7 and T47D cells, both HSF1-proficient and HSF1deficient. However, within the time given for the revision, these new experiments were not concluded due to several technical problems. Nevertheless, by re-analyzing the presence of HSE and ERE in HSF1 and the revision ChIP-seq peaks or ChIA-PET reads, we provided a model of possible interactions shown in the new Figure 5.

11. In breast cancer, ER+ status is an independent and mostly favorable prognostic marker (PMIDs: 6488142, 28601929). High HSF1 expression is associated with mortality in ER+ breast cancer (PMID: 22042860, 27713164). Prognostic impact of HSF1 expression on ER- cases was not obvious (PMID: 27713164). Here, the authors analyzed the TCGA database for prognostic value of ER and HSF1 mRNA and suggest that ER- cases with high HSF1 mRNA levels have the poorest outcome (Figure 6D), suggesting an ER-independent role for HSF1 in breast cancer. However, this conclusion is probably a biased reflection of the prognostic difference between basal and luminal breast cancer subtypes (Figure 6C) (PMID: 26693050). Moreover, how this finding supports a HSF1 role as a cooperative or tethering factor for ER is not clear, and the prognostic value of HSF1 mRNA levels in ER+ breast cancer was not obvious from the data presented.

Inconsistent data regarding the prognostic role of HSF1 are present in the available literature and databases, which are discussed below. Moreover, the presentation of own results was improved in the revised manuscript, as abbreviated below.

Our analyses (shown in revised Figure 6—figure supplement 1) and analyses of TCGA data available in the Human Protein Atlas (https://www.proteinatlas.org/ENSG00000185122-HSF1/pathology/breast+cancer https://www.proteinatlas.org/ENSG00000091831-ESR1/pathology/breast+cancer) indicate that neither HSF1 nor ERα are prognostic in breast cancer when there is no pre-selection of patients. A similar analysis performed using Kaplan-Meier Plotter, Pan-cancer RNA-seq, confirmed these results (https://kmplot.com/analysis/index.php?p=service&cancer=pancancer_rnaseq). On the other hand, however, the better prognosis of breast cancer patients with high ESR1 levels (logrank p < 1E-16) or low HSF1 levels (p = 0.00089) is visible in Kaplan-Meier Plotter, Breast Cancer analyses based on Affymetrix gene chip (in the case of HSF1, only when Auto select best cutoff is used to split patients) (https://kmplot.com/analysis/index.php?p=service&cancer=breast#).

Importantly, however, HSF1 and ERα may have prognostic value when TCGA dataset analyses are restricted to selected subtypes. For example, high HSF1 level is connected with better prognosis in stage 1 (logrank p = 0.021), white race (p = 0.049), and Asian race (p = 0.016), while with worse prognosis in stage 3 (p = 0.05), stage 4 (p = 0.0043), grade 1 (p = 0.021), grade 2 (p = 0.013), and African American race (p = 0.0075). High ESR1 level is connected with with better prognosis in stage 3 (logrank p = 0.015), stage 4 (p = 0.022), grade 1 (p = 0.02), grade 2 (p = 1.3e-05), grade 3 (p = 0.042), African American race (p = 0.0098), and Asian race (p = 1e-04), while with worse prognosis in stage 1 (p = 0.013), and white race (p = 0.034). In our selected groups caucasians (white race) predominates. On the other hand, caucasians are more likely to be diagnosed with luminal A, while African Americans are more likely to be diagnosed with basal-like.

To delineate the clinical and prognostic significance of HSF1 in breast cancer two analyses have been performed so far, by (Santagata et al., 2011) and (Gökmen-Polar and Badve, 2016)(mentioned by the reviewer and discussed in our manuscript). In the paper by Santagata et al., the association of nuclear HSF1 status (high vs low, detected by IHC) with clinical parameters and survival outcomes in 1,841 breast cancer patients (diagnosed between 1976 and 1996) were investigated. They revealed a significant correlation between high HSF1 levels and mortality of ERpositive patients (1416 cases) but not ER-negative ones ( triple-negative 267 cases, HER+ 193 cases). However, in the other 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ökmenPolar and Badve, 2016). Studies by Sanatagata et al., and Gokmen-Polar and Badve showed a significant correlation between HSF1 transcript levels and the survival in ER-positive breast cancer patients. However, in the paper of Satagata et al., data from (van de Vijver et al., 2002) were utilized. In that study, patients were preselected to have a good or bad prognosis to better evaluate the predictive power of the prognosis profile: patients (295 cases) with primary breast carcinomas as having a gene-expression signature associated with either a poor prognosis or a good prognosis were chosen. All patients had stage I or II breast cancer and were younger than 53 years old; 151 had lymph-node-negative disease, and 144 had lymph-node-positive disease. Thus, they are not representative of all breast cancers. Gokmen-Polar and Badve analyzed publicly available Affymetrix U133A gene expression data, thus a bit earlier than those available in TCGA. Our results based on the TCGA cohort are different from those previously published. At the moment, we can only speculate on the cause. Generally, different sets of patients were compared. We should keep in mind that early screening tests and new treatments have improved the survival rate of breast cancer in women (Iwase et al., 2021) and this may be one of the reasons for the difference between these cohorts. Especially HER2-positive cancers have a much better prognosis nowadays.

We are aware that in our selected ER+/− groups, differences in survival may simply reflect the prognostic difference between breast cancer subtypes. In our analysis, we excluded cases with moderate HSF1 and ERα levels and compared only the most divergent groups to better see the possible effect of HSF1 and ERα on the transcriptome. Such selection enabled us to reveal a statistically significant difference in survival between ER+ and ER− patients as well as HSF1+ and HSF1− (the effect was not observed when all TCGA cases were analyzed). However, the effect of HSF1 was poorly visible in each ER group (but slightly separated ER− cases, regardless of the size of the analyzed groups, see Figure 6—figure supplement 1A, B, C). Significantly, luminal A cases were dominant in our HSF1low/ER+ group. Increased expression of HSF1 was connected with the enrichment of the ER+ group with luminal B cases, known to have a slightly worse prognosis than luminal A. Analyzes of the survival probability in the reduced (selected) groups (defined in Figure 7B) containing only basal-like or luminal A cancers showed that high levels of HSF1 may worsen the prognosis of ER-negative patients but not of luminal A (revised Figure 6—figure supplement 1C). Nonetheless, experimental data (e.g. HSF1 impact on cell migration or response to anti-cancer treatments) supports the idea that HSF1 may influence the probability of survival.

We also analyzed pre/postmenopausal breast cancers for survival probability, obtaining different results depending on ER+ and ER− inclusion criteria. However, we do not have clear conclusions from these analyzes because there is a lack of important information about the patient's condition that may have an impact on survival (among others obesity, contraceptive or hormone replacement therapies, etc.). Considering that the HSF1-dependent survival rate may be multivariable, in the revised manuscript we focused on other aspects of HSF1 action in breast cancer that were documented in new set of data. Thus, the subtitle “The combined expression level of ESR1 and HSF1 can be used to predict the survival of breast cancer patients” was changed to: “Metastatic and nonmetastatic breast cancers differ in the level of HSF1 only in the ER+ group”. The abstract also was changed accordingly.

12. With regard to their conclusion that HSF1 increases the diversity of the transcriptome in ER-positive breast cancers, the rationale for this analysis is unclear, and the idea that HSF1 affects ER-target genes has already been demonstrated in the KO model using RNA-seq (Figure 2).

This is true, but using patient/cancer-derived data we wanted to support the results obtained in the cellular model.

13. Stylistic changes are needed: Much of the text in the figures are too small to read. It was hard to tell the differences between the A, B, C, D groups and the a, b, c, d groups in figure 7.

The initial idea was about maintaining connections between initial and reduced patient groups (patients from group “A” were selected to group “a”, etc.) and avoiding typing the full description of the groups each time, which appeared confusing. Hence, the group description was changed in the revised figure: initial groups are numbered from I to IV and reduced groups are marked with superscript R (IR to IVR). In the revision, figures are uploaded separately (not embedded in the text), thus we believe that their quality and readability are better.

14. In figure 5, they authors should be subsampling the "all" groups to the same number as the smaller groups to see if the distributions are the same.

15. Figures 5A and 5B are trying to show information for the HSF1 and ESR1 overlapping peaks. However, for that data to be useful we need to understand non-overlapping peak data. Also, is the figure showing tags from ESR1 ChIP-seq or HSF1 ChIP-seq?

16. Important concern: The authors should be using statistical tests to see if overlap in data sets are above random chance.

Points 14-16. In the revised manuscript we tested the distributions of ERα and HSF1 ChIP-seq peaks using Jaccard bedtool as well as peakPermTest function from the ChIPpeakAnno tool, which revealed that there is a significant overlap between two sets of peaks. The revised manuscript was supplemented with this information (see the response to Main Concerns 3 and 4 for detail). Also, we included the Venn diagram in the relevant figure (revised Figure 4B) showing better overlapping and nonoverlapping peak data.

17. The authors generated HSF-1 knockouts for MCF7 and T47D breast cancer cells. However, the vast majority of their work was performed with only the T47D cells. These cell lines show different dependencies on ERα and antiestrogen differently affect receptor stability between them. More work should have been done outside of one knockdown and proliferation experiment to show that HSF-1 is important in both cell lines. It would give credence to its role in breast cancer cells beyond MCF7 cells.

18. The authors used an aldefluor assay to show that HSF-1 affects stem-progenitor populations. They did not show any additional data for this line of inquiry. If this phenotype is affected by HSF-1 they should examine EMT signatures along with cellular invasion, and migration.

Points 17-18. In the revised manuscript, we showed additional data documenting the effect of HSF1 on the phenotype of both MCF7 and T47D cells (revised Figure 1, Figure 1—figure supplement 1, and Figure 1—figure supplement 2). In addition to the proliferation and clonogenic assay, we showed results from estrogen-stimulated migration in the Boyden chamber, cell adhesion, and morphology. For simplicity, we decided to show in the main manuscript only the results from the HSF1-deficient MCF7 model obtained with the DNA-free CRISPR/Cas9 system, which appeared the most relevant and representative (revised Figure 1). The same functional tests performed on other models are shown in the revised supplement (revised Figure 1—figure supplement 1 and Figure 1—figure supplement 2).

HSF1 deficiency in MCF7 cells may result in a slight reduction in the proliferation rate under standard conditions. Unlike MCF7 cells, HSF1-deficient T47D cells grew slightly faster than HSF1proficient cells. T47D cells harbor a p53 missense mutation (L194F). Although p53 mutations were found in approximately 5% of all breast cancer cases, L194F mutation is only present in 0.38% of cases (TCGA data), so we believe that the results obtained on this line, while valuable, are less representative. L194F mutation causes p53 stabilization (Lim et al., 2009) and generally T47D cells differ from MCF7 cells in response to estrogen (lower level of ERα in T47D; (Vydra et al., 2019)). Also, the mutant p53 exhibits gain-of-function activities in mediating cell survival of T47D cells. This is likely the reason for the differences in proliferation between cells.

Both cell lines (T47D and MCF7) also respond differently to prolonged estrogen treatment, especially in acquired cell morphology: amoeboid-like morphology was dominant among the scattered MCF7 cells, while mesenchymal-like morphology was dominant in T47D cells. HSF1 deficiency resulted in inhibition of estrogen-mediated cell migration (as assessed by Boyden assay), but this effect was more pronounced in MCF7 cells. However, the differences in EMT markers were not documented in these cell lines. Results of in vitro analysis that showed an effect of HSF1 on E2stimulated cell migration prompted us to analyze the TCGA data for HSF1 expression and the presence of metastases at diagnosis. We found that an elevated HSF1 expression correlated with the occurrence of metastatic disease only in ER-positive breast cancer patients (revised Figure 6F-H and Figure 6—figure supplement 2), which suggests that HSF1 in cooperation with ERα may facilitate metastasis formation.

19. Therapeutic approaches to targeting ERα were not discussed in this manuscript. However, the role of HSF-1 in the therapeutic response to antiestrogens would be of great importance to the field. Showing the impact of HSF-1 knockdown to ERα transcriptional activity and cellular proliferation in the presence of clinically important hormone therapies fulvestrant and 4-hydroxytamoxifen as well as CDK4/6 inhibitors like palbociclib would greatly enhance the impact of this paper.

Thank you for this comment, which prompted us to perform additional experiments related to the role of HSF1. 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 (Gokmen-Polar and Badve, 2016). Recently, it has been shown that overexpression of HSF1 induces resistance to antiestrogens (tamoxifen, fulvestrant) in MCF7 cells (Silveira et al., 2021). Our new results (included in the revised manuscript) also showed that the response to hydroxytamoxifen is better in HSF1− than in HSF1+ cells, although the HSF1− knockout itself has a different effect on the growth of MCF7 and T47D cells (revised Figure 8). Palbociclib is a relatively new therapeutic option and the effect of HSF1 on its effectiveness has not yet been investigated. We found that HSF1 functional knockout was connected with higher sensitivity to palbociclib but this effect was better visible in MCF7 than T47D cells (which generally were more resistant to palbociclib than MCF7 cells). These results 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.

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

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.

Senior Editor

  1. Anna Akhmanova, Utrecht University, Netherlands

Reviewing Editor

  1. Maureen E Murphy, The Wistar Institute, United States

Reviewers

  1. Mary Ann Allen, Howard Hughes Medical Institute, University of Colorado, United States
  2. Sean W Fanning, University of Chicago, United States

Publication history

  1. Received: April 28, 2021
  2. Preprint posted: May 7, 2021 (view preprint)
  3. Accepted: November 15, 2021
  4. Accepted Manuscript published: November 16, 2021 (version 1)
  5. Version of Record published: December 24, 2021 (version 2)

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,338
    Page views
  • 261
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

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

Categories and tags

Research organism

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology

    Comprehensive characterization of tumor microenvironment in colorectal cancer via molecular analysis

    Xiangkun Wu, Hong Yan ... Li Liang
    Research Article

    Colorectal cancer (CRC) remains a challenging and deadly disease with high tumor microenvironment (TME) heterogeneity. Using an integrative multi-omics analysis and artificial intelligence-enabled spatial analysis of whole-slide images, we performed a comprehensive characterization of TME in colorectal cancer (CCCRC). CRC samples were classified into four CCCRC subtypes with distinct TME features, namely, C1 as the proliferative subtype with low immunogenicity; C2 as the immunosuppressed subtype with the terminally exhausted immune characteristics; C3 as the immune-excluded subtype with the distinct upregulation of stromal components and a lack of T cell infiltration in the tumor core; and C4 as the immunomodulatory subtype with the remarkable upregulation of anti-tumor immune components. The four CCCRC subtypes had distinct histopathologic and molecular characteristics, therapeutic efficacy, and prognosis. We found that the C1 subtype may be suitable for chemotherapy and cetuximab, the C2 subtype may benefit from a combination of chemotherapy and bevacizumab, the C3 subtype has increased sensitivity to the WNT pathway inhibitor WIKI4, and the C4 subtype is a potential candidate for immune checkpoint blockade treatment. Importantly, we established a simple gene classifier for accurate identification of each CCCRC subtype. Collectively our integrative analysis ultimately established a holistic framework to thoroughly dissect the TME of CRC, and the CCCRC classification system with high biological interpretability may contribute to biomarker discovery and future clinical trial design.

    1. Cancer Biology

    Defining function of wild-type and three patient specific TP53 mutations in a zebrafish model of embryonal rhabdomyosarcoma

    Jiangfei Chen, Kunal Baxi ... Myron S Ignatius
    Research Article

    In embryonal rhabdomyosarcoma (ERMS) and generally in sarcomas, the role of wild-type and loss or gain-of-function TP53 mutations remains largely undefined. Eliminating mutant or restoring wild-type p53 is challenging; nevertheless, understanding p53 variant effects on tumorigenesis remains central to realizing better treatment outcomes. In ERMS, >70% of patients retain wild-type TP53, yet mutations when present are associated with worse prognosis. Employing a kRASG12D-driven ERMS tumor model and tp53 null (tp53-/-) zebrafish, we define wild-type and patient-specific TP53 mutant effects on tumorigenesis. We demonstrate that tp53 is a major suppressor of tumorigenesis, where tp53 loss expands tumor initiation from <35% to >97% of animals. Characterizing three patient-specific alleles reveals that TP53C176F partially retains wild-type p53 apoptotic activity that can be exploited, whereas TP53P153D and TP53Y220C encode two structurally related proteins with gain-of-function effects that predispose to head musculature ERMS. TP53P153D unexpectedly also predisposes to hedgehog expressing medulloblastomas in the kRASG12D-driven ERMS-model.

    1. Cancer Biology

    Proteomic characteristics reveal the signatures and the risks of T1 colorectal cancer metastasis to lymph nodes

    Aojia Zhuang, Aobo Zhuang ... Chen Ding
    Research Article Updated

    The presence of lymph node metastasis (LNM) affects treatment strategy decisions in T1NxM0 colorectal cancer (CRC), but the currently used clinicopathological-based risk stratification cannot predict LNM accurately. In this study, we detected proteins in formalin-fixed paraffin-embedded (FFPE) tumor samples from 143 LNM-negative and 78 LNM-positive patients with T1 CRC and revealed changes in molecular and biological pathways by label-free liquid chromatography tandem mass spectrometry (LC-MS/MS) and established classifiers for predicting LNM in T1 CRC. An effective 55-proteins prediction model was built by machine learning and validated in a training cohort (N=132) and two validation cohorts (VC1, N=42; VC2, N=47), achieved an impressive AUC of 1.00 in the training cohort, 0.96 in VC1 and 0.93 in VC2, respectively. We further built a simplified classifier with nine proteins, and achieved an AUC of 0.824. The simplified classifier was performed excellently in two external validation cohorts. The expression patterns of 13 proteins were confirmed by immunohistochemistry, and the IHC score of five proteins was used to build an IHC predict model with an AUC of 0.825. RHOT2 silence significantly enhanced migration and invasion of colon cancer cells. Our study explored the mechanism of metastasis in T1 CRC and can be used to facilitate the individualized prediction of LNM in patients with T1 CRC, which may provide a guidance for clinical practice in T1 CRC.