Integrated bioinformatics and experiments reveal the roles and driving forces for HSF1 in colorectal cancer

ABSTRACT Heat shock factor 1 (HSF1) has watershed significance in different tumors. However, the roles and driving forces for HSF1 in colorectal cancer (CRC) are poorly understood. Our study integrally analyzed the roles and driving forces for HSF1 in CRC by bioinformatics and experiments. The expression and prognostic characteristics of HSF1 were analyzed via UALCAN, GEPIA2, TISIDB, Prognoscan and HPA databases. Then, we analyzed the correlation between HSF1 expression and immune features via TIMER2 database. Subsequently, we explored the driving forces for HSF1 abnormal expression in CRC by bioinformatics and experiments. Our results showed that HSF1 was overexpressed and correlated with poor prognosis in CRC. And the expression of HSF1 was significantly correlated with multiple immune cell infiltration and was negatively correlated with immunomodulators such as programmed cell death 1 ligand 1(PD-L1). Along with many driver genes in particular TP53, super-enhancer, miRNA and DNA methylation were all responsible for HSF1 overexpression in CRC. Moreover, we demonstrated that β-catenin could promote the translation process of HSF1 mRNA by interacting with HuR, which could directly bind to the coding sequence (CDS) region of HSF1 mRNA. Collectively, HSF1 may be useful as a diagnostic and prognostic biomarker for CRC. HSF1 was closely correlated with immune features. Genetic and epigenetic alterations contributed to HSF1 overexpression in CRC. More importantly, we demonstrated that HSF1 may be regulated at the level of mRNA translation by β-catenin-induced HuR activity.


Introduction
Based on the latest data published by the International Agency for Research on Cancer (IARC), about 1.93 million new colorectal carcinoma (CRC) cases and 0.93 million deaths occurred in 2020 worldwide, accounting for about 10% of all new cancer cases and 9.4% of cancer deaths [1]. Furthermore, the morbidity and mortality of CRC is still rising in many countries with poor prognosis. Hence, the exploration of new biomarkers is essential for the early diagnosis and treatment of CRC.
With the development of high-throughput sequencing technology, huge amounts of databases such as TCGA (the cancer genome atlas) and GEO (gene expression omnibus) have emerged [2,3]. It is beneficial for us to analyze the genes of interest by bioinformatics, providing novel potential biomarkers. Heat shock factor 1 (HSF1) is a master transcription factor for modulation of the heat-shock response and the proteotoxic stress response. HSF1 has a carcinogenic role in regulating proliferation, survival, invasion and metastasis of cancer cells [4,5]. However, the roles of HSF1 in colorectal cancer are poorly understood. The latest research showed that HSF1 could induce protective autophagy for oxaliplatin resistance [6]. Our previous study reported that HSF1 promoted colorectal carcinogenesis by stimulating glutamine metabolism [7]. New exploration for the roles and driving forces of HSF1 in CRC is indispensable for targeting this molecular.
This study was dedicated to analyze the clinical characteristics for HSF1 in CRC by bioinformatics. More importantly, we were the first study to integrally reveal the driving forces for HSF1 overexpression in colorectal cancer by bioinformatics and experiments, providing theoretical basis for targeting HSF1.

Expression analysis
The UALCAN [8] and Oncomine database [9] were used to investigate the different expression of HSF1 in colorectal cancer and normal tissues. Concomitantly, HSF1 mRNA expression was also compared at different CRC stages using the UALCAN, TISIDB [10] and GEPIA2 database [11]. P value cutoff is 0.05.

Prognosis analysis
The relationship between the HSF1 expression and the prognostic significance including overall survival (OS), disease-free survival (DFS) and diseasespecific survival (DSS) in colorectal cancer was validated by dedicated prognostic database -PrognoScan [12]. Meanwhile, The Human Protein Atlas (HPA) database [13] (http://www. proteinatlas.org/) was used to further explore the prognostic value. P value cutoff is 0.05.

Functional enrichment analysis
The top 50 genes related with HSF1 in colon cancer and rectal cancer were speculated separately by UALCAN database. HSF1-related genes were input into DAVID [14] database for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis. GO analysis contains biological process (BP), molecular function (MF) and cellular component (CC). The related results were visualized in Bubble plots. P value cutoff is 0.05.

Immune feature analysis
The relationship between the HSF1 expression and the immune cell infiltration was analyzed by TIMER2 database [15]. P value cutoff is 0.05. The correlation between the expression of HSF1 and immunomodulators including PD-L1 (CD274), CTLA4 and PD-L2 (PDCD1LG2) was analyzed by TIMER2 and StarBase [16] database. P value cutoff is 0.05.

Mutations analysis
The somatic mutation ratio of HSF1 was elucidated using TIMER2 database and COSMIC database [17]. Then, the online tool-ACLBI (https:// www.aclbi.com/) was used to analyze the somatic mutation in 536 colorectal cancer samples from the TCGA database and explore the mutation characteristics of the top 20 high frequency mutated genes as well as HSF1. The data of cBioPortal database [18] showed the mutation types of HSF1 in CRC. The correlation between HSF1 mRNA expression and copy number variation was analyzed by cBioPortal database in colorectal cancer. P value cutoff is 0.05.

Driver gene and promoter methylation analysis
The correlation between HSF1 expression and driver gene was illustrated through the TCGA portal database [19]. MethPrimer [20] was used to evaluate the CpG Island in HSF1 promoter. Subsequently, the association between the HSF1 expression and DNA methylation status in colorectal cancer was evaluated using cBioPortal database. Ultimately, DiseaseMeth version 2.0 [21] and UALCAN were used to estimate methylation levels of HSF1 between the CRC and corresponding adjacent tissues. P < 0.05 was regarded as statistically significant.

Human tissue samples
Our study was approved by The Institute Research Medical Ethics Committee of Affiliated Hospital of Weifang Medical University. All patients signed informed consent forms. The samples of CRC patients were obtained from the Affiliated Hospital of Weifang Medical University from 1 January 2021 to 30 June 2021. None of the patients received other treatment before surgery. The samples were obtained immediately after surgery and stored at −80°C.

SiRNA transfection
For this, 1 × 10 5 cells was plated in 6-well plates and

Western blot
After siRNA transfection, the cells in the 6-well plates were lysed by RIPA buffer. Cell lysates were separated and transferred to PVDF membranes. The PVDF membranes were incubated with primary antibodies at 4°C overnight. Then, the membranes were incubated with the secondary antibody conjugated with HRP (Abcam, ab6721) at room temperature for 1 h. Finally, the membranes were visualized with chemiluminescence (Solarbio, PE0010).

Co-immunoprecipitation (CO-IP)
For this, 5 × 10 6 cells in 10 cm dish were lysed by NP40 buffer (Solarbio, N8031) containing 1× cocktail (Roche, 05892970001). Five hundred micrograms of protein was incubated with 5 µg RIgG or HSF1 antibody for overnight with rotating at 4°C. Then, 50 µl agarose beads was added to the mixture for 3 h with rotating at 4°C. The mixture was subjected to Western blot for detection.

RNA immunoprecipitation (RIP)
For this, 1 × 10 7 cells were lysed using RNA-Binding Protein Immunoprecipitation Kit (Millipore, No. 17-700) according to the manufacturer's instructions. The lysates were immunoprecipitated with 5 µg RIgG or HuR antibody for overnight with rotating at 4°C. Total RNA was extracted and was performed with reverse transcription and real-time quantitative PCR.

Results
Our study integrally analyzed the roles and driving forces for HSF1 in CRC by bioinformatics and experiments. We explored the expression and prognostic characteristics of HSF1 in CRC and performed the functional enrichment analysis as well as immune feature analysis. Subsequently, we explored the driving forces for HSF1 abnormal expression in CRC such as genetic alterations, epigenetic alterations and translation efficiency.

The expression and the prognostic characteristics of HSF1 in CRC
We first analyzed the expression of HSF1 mRNA in CRC by bioinformatics. HSF1 was significantly overexpressed in carcinoma than normal tissue by UALCAN and Oncomine database (Figure 1a-b). Moreover, HSF1 expression level was correlated with tumor stage and nodal metastasis status in colon adenocarcinoma (COAD) by GEPIA2, TISIDB, UALCAN database (Figure 1c-f). These results demonstrated that aberrant HSF1 expression may occur from early stages of COAD. Subsequently, we investigated the clinical significance of HSF1 expression in CRC patients. Based on the Kaplan-Meier survival curves, PrognoScan and HPA database analysis illustrated that high expression of HSF1 was significantly associated with poor disease-free survival (DFS) and overall survival (OS) in the CRC (Figure 1g and h). In summary, the aforementioned data indicated that HSF1 may be useful as a diagnostic and prognostic biomarker for CRC.

Functional enrichment analysis
To further explore the roles of HSF1 in CRC, we screened out HSF1-related genes involved in COAD and rectum adenocarcinoma (READ) respectively using UALCAN database. We discovered 596 common members and visualized the top 50 gene networks through the GeneMANIA [22] database (Figure 2a and b). Among them, ZC3H3, SCRIB, PUF60 and SHARPIN showed the highest correlation with HSF1 expression in COAD ( Figure 2c). While in READ, HSF1 expression was highly correlated with SCRIB, MAF1, ZC3H3 and SHARPIN ( Figure 2d). Next, we performed the functional enrichment analysis of all relevant genes via the DAVID database. Figure 2e showed the top 10 most significantly enriched GO and KEGG pathways. In biological process (BP) term, transcription, DNA-templated and rRNA processing were mostly enriched. Nucleus, cytoplasm and nucleoplasm were significantly enriched in cellular components (CC). The molecular function (MF) mainly contained the protein binding, DNA and RNA binding. As for KEGG pathway, endocytosis, thyroid hormone and Hippo signaling pathway were significantly enriched.

Immune feature analysis
With the popularity of immunotherapy, we wondered whether HSF1 was correlated with immune features. First, we analyzed the infiltration level of immune cells via TIMER2 database. The results showed that HSF1 expression was negatively correlated with CD4 + T cell, B cell, Treg cell, NK cell and myeloid dendritic cell in COAD (Figure 3a). And HSF1 expression was positively correlated with cancer associated fibroblast (CAF), myeloid derived suppressor cells (MDSC) and macrophage in COAD. While in READ, HSF1 expression was negatively correlated with CD4 + T cell and positively correlated with cancer associated fibroblast (CAF), CD8 + T cell and NK cell (Figure 3b). More importantly, the expression of HSF1 was negatively correlated with the expression of immunomodulators such as PD-L1(CD274), CTLA4 and PD-L2(PDCD1LG2) (Figure 3c-e). Collectively, these results revealed that HSF1 may be involved in tumor immune regulation.

Mutation status of HSF1 in CRC
Now that HSF1 has a critical role in CRC, we are keen to explore the driving forces for HSF1 to target it. In tumor cells, copy number variation (CNV) and mutation counts are important factors affecting the expression of related genes [23,24]. First, we illustrated the mutation features of HSF1 in tumors exploiting the TIMER2 database. The results demonstrated that the HSF1 mutation rate was low in all tumors including COAD-about 2.4% ( Figure S1A). The top 20 high-frequency mutated genes in colorectal carcinoma were shown by online tool-ACLBI ( Figure S1B) . In 536 colorectal cancer samples, the overall mutation rate was about 97.95% (525 out of 536) and the mutant frequency of APC was up to 75%, while the mutation rate of HSF1 was only 1%. The data of cBioPortal database showed the mutation types of HSF1 in CRC ( Figure S2A). We summarized the CNV rate and gene expression rate of HSF1 in CRC through Catalog of Somatic Mutation in Cancer (COSMIC) database analysis. As shown in Figure S1C and S1D, the CNV rate of HSF1 was just 4.32%, while the gene expression rate was 29.67%. These data demonstrated that HSF1 mutation is one of the mechanisms leading to the upregulation of HSF1 expression among CRC patients. However, CNV does not fully explain the overexpression of HSF1 in CRC. Therefore, we need to further investigate other mechanisms that mediate HSF1 upregulation in CRC.

The effect of genetic and epigenetic alterations on HSF1 overexpression
To further investigate the reasons for the high expression of HSF1 in CRC, we analyzed the driver genes that are the critical nodes of regulatory networks and signaling pathways [25]. By TCGA portal database, we found HSF1 expression was correlated with multiple driver genes including APC, TP53, KRAS and PIK3CA (Figure 4a). Noteworthy, among them, TP53 showed the highest correlation with HSF1 expression. Consistent with this result, HSF1 expression was significantly associated with mutant TP53 in the UALCAN database ( Figure S2B). Besides genetic changes, tumor is also attributed to epigenetic alterations including histone modification, non-coding RNAs and DNA methylation. Histone modification such as acetylation, methylation, phosphorylation, etc. is critical epigenetic mechanism. In particular, high levels of histone H3 lysine 27 acetylation (H3K27ac) lead to the formation of superenhancer, which is responsible for the control of key genes identity [26]. We observed that high H3K27ac signal exist in the transcription initiation site (TSS) of HSF1 by UCSC database (Figure 4b).
Additionally, the expression of HSF1 is positively with the expression of BRD4, the master reader that binds to acetylated histones and regulates gene transcription (Figure 4c and d).
Pharmacological inhibition or genetic deficiency of BRD4 could suppressed the expression of HSF1 (Figure 4e and f). All of these results implied the possible role of super-enhancer in regulating HSF1 expression.
Non-coding RNAs, especially microRNA (miRNA), could regulate gene expression by inducing mRNA degradation or inhibiting mRNA translation [27]. Subsequently, we screened out the potential miRNAs for HSF1 by TargetScan and miRWalk database [28,29]. There were about 33 common members (Figure 5a). Among them, only has-miR-378a-5p and has-miR-874-3p was downregulated and negatively correlated with the expression of HSF1 in CRC via StarBase database (Figure 5b-d). DNA methylation is closely related to gene expression [30]. Intriguingly, the MethPrimer website showed that there was a CpG Island in the HSF1 promoter region (Figure 5e). We further explored the relationship between HSF1 and DNA methylation in CRC by DiseaseMeth database. The analysis elucidated that HSF1 showed lower DNA methylation levels in CRC compared with normal tissues (figure 5f). A similar trend was observed in UALCAN database ( Figure S2C). In addition, cBioPortal database revealed an obviously negative correlation between HSF1 expression and DNA methylation level in CRC ( Figure S2D). Taken together, DNA methylation may provide a cause for HSF1 overexpression.

HSF1 is regulated at the level of mRNA translation
Many key controls of gene expression occur at the level of mRNA translation ensuring the synthesis of cellular proteins under specific conditions such as heat shock and starvation [31]. We explored the translatome information of HSF1 mRNA by POSTAR3 database [32]. We found that there were multiple open reading frames (ORFs) for HSF1 mRNA among different tumor cells (Figure 6a). Next, we analyzed the most common transcript-ENST00000528838 and discovered that the translatome signal of this transcript in HCT116 was significantly increased (Figure 6b). Meanwhile, the translation efficiency and density of this transcript was higher in many tumor cells including HCT116 (Figure 6c and d). The above data demonstrated that HSF1 may be regulated at the level of mRNA translation in CRC.

HuR co-expressed with HSF1 in CRC
RNA binding proteins (RBPs) are commonly considered to combine RNA via globular RNA binding domains (RBDs) and alter the fate or function of the target RNA [33]. They are involved in many aspects of posttranscriptional regulation such as RNA processing, RNA transport, localization, degradation and translation efficiency [34]. Now that the translation level of HSF1 mRNA in colorectal cancer is high, we wanted to explore the certain RNA binding protein. To validate this hypothesis, we first screened out common RBPs using the RBPDB online database primarily focusing on collecting experimentally validated RBPs and RBDs [35]. The filter criteria were set to human, number of experiments (≥20) and Human Homologs (≥10). Finally, we selected 2 most common RBPs-ELAVL1 (HuR) and HNRNPA1 (Figure 7a). However, there is no significant correlation between the expression of HNRNPA1 and HSF1 in CRC ( Figure S3A-B). While we found that the expression of HuR was positively correlated with that of HSF1 in CRC (Figure 7b-c). Furthermore, the expression levels of ELAVL1 were significantly higher in CRC than their corresponding normal tissues (Figure 7d and Figure  S3C). Prognostic analyses from the PrognoScan database confirmed that patients with higher expression of ELAVL1 mRNA had significantly shorter overall survival (OS) and disease-specific survival (DSS) (Figure 7e). Subsequently, the relative protein expression of HSF1 and HuR in eight pairs of CRC tissues and their adjacent normal tissues were analyzed by Western blot assay. The results demonstrated that the expression of HSF1 and HuR was both significantly elevated and positively correlated in most of CRC patients ( Figure 7f and Table 1). Table 2 shows the correlation between HSF1 expression and clinicopathological parameters in CRC patients (n = 8). Accordingly, HuR may be a potential RBP responsible for HSF1 overexpression.

HuR promoted the translation of HSF1 by binding to its CDS region
HuR is critical for regulating the translation process of oncogenic proteins [36]. To further verify whether HuR contributed to the translation of HSF1 mRNA, we first compared the protein expression of HSF1 and HuR in CRC cell lines. We found that both of them were upregulated in CRC cell lines compared with human normal intestinal epithelial cell CCD841 (Figure 8a). As expected, the protein level of HSF1 was downregulated after knockdown of HuR but not the mRNA level (Figure 8b-c). Normally, HuR interacted with one or several U or AU-rich elements (AREs) in the 3'-UTR (untranslated region) [37]. Interestingly, we observed that exogenous HSF1 with CDS region but not the 3'-UTR was also reduced after knockdown of HuR, demonstrating that HuR may promote the translation of HSF1 by binding to its CDS region (Figure 8d). Indeed, we found two fragments containing GUUUG sequences in the CDS region of HSF1 mRNA (Figure 8e). Moreover, we identified that HuR could directly bind to the two fragments of HSF1 CDS region as well as 3'-UTR by RNA immunoprecipitation (RIP) experiment (Figure 8f-g and Figure S4A). Collectively, these results implied  that the translation of HSF1 was not only entirely dependent on the 3'-UTR but HuR could also play its role by the CDS region.

β-Catenin co-operated with HuR to promote the translation of HSF1
Our previous study has reported that β-catenin (CTNNB1) can positively modulate HSF1 translation [38]. Therefore, we wonder whether there was a correlation between β-catenin and HuR. Similarly, exogenous HSF1 with CDS region but not the 3'-UTR was also reduced after knockdown of β-catenin ( Figure 9a). Additionally, as shown in Figure 9b-c and Figure S4B, there was a positive correlation between CTNNB1 and HuR in CRC. Nevertheless, knockdown of β-catenin did not significantly alter the expression of HuR (Figure 9d). Notably, we discovered an interaction between β-catenin and HuR via co-immunoprecipitation (CO-IP) experiment ( Figure 9e). Likewise, we found that the interaction between them scored highly by the HitPredict database [39] ( Figure S4C). More importantly, the binding level between HuR and HSF1 mRNA was decreased after silencing β-catenin (figure 9f). Consequently, we  hypothesized that HuR was able to stimulate the translation of HSF1 mRNA under the activity of βcatenin.

Discussion
HSF1 is well known to be an important regulator of the heat shock reaction. The role of HSF1 in cancer has received widespread attention in recent years. Different from the heat shock response, HSF1 could drive a transcriptional program to endorse the malignant phenotype of cancer cells [4]. Our previous studies had verified that HSF1 is highly expressed in colorectal cancer patients. In present study, we further explore the roles and driving forces for HSF1 overexpression in CRC by bioinformatics and experiments.
We confirmed that HSF1 may be useful as a diagnostic and prognostic biomarker for CRC. Meanwhile, we performed the functional enrichment analysis, revealing the underlying molecular mechanism. Notably, our results implied that HSF1 was closely correlated with immune features in CRC. There was a discrepancy about the immune cell infiltration in COAD and READ. Extraordinarily, a significant positive correlation was both observed in COAD and READ between HSF1 expression and the abundance of cancerassociated fibroblast (CAF). Latest study reported that HSF1 is also essential for extracellular matrix remodeling in CRC by CAF [40], highlighting HSF1 vital roles in regulating tumor microenvironment. Interestingly, HSF1 could induce PD-L1 expression and enhance tumor growth in breast cancer [41]. However, our findings demonstrated that in CRC HSF1 was negatively correlated with the expression of immunomodulators such as PD-L1, CTLA4 and PD-L2. This result is consistent with previous findings that HSF1 expression was negatively correlated with tumor mutation burden (TMB) and microsatellite instability (MSI) in CRC [42]. Accordingly, HSF1 could be a predictive biomarker for immunotherapy in CRC, which is beneficial for selecting suitable patients to receive immunotherapy.
The present study about the driving forces for HSF1 overexpression in CRC is still limited. Tumorigenesis is a multistep process that involves genetic and epigenetic alterations, contributing to aberrant proteins and subsequent abnormal functions [43,44]. In this study, our data indicated that mutation of HSF1 may be partly responsible for the upregulation of HSF1 in CRC patients. While overexpressed HSF1 was subject to multiple driver genes in CRC, especially TP53 which showed the highest correlation. This result is in accordance with Isermann's study that mutp53 could unleash HSF1 function by eliminating the repressive WTp53-HSF1 axis [45]. Subsequently, we explored the effect of epigenetic alterations on HSF1 expression including histone modification, non-coding RNAs and DNA methylation. Unexpectedly, we found a super-enhancer in the TSS of HSF1 mRNA. Different from typical enhance, superenhancer is a group of active-enhancer clusters which is more densely occupied by transcription factors and other histone regulators. Therefore, super-enhancer is indispensable for regulating the function of important genes and cell fate [46]. Targeting this super-enhancer may impose restrictions on the function of HSF1 in CRC. DNA methylation is one of the predominant epigenetic mechanisms affecting gene expression, which could regulate gene transcription levels depending on the methylation status of target genes [47,48]. Multiple databases showed a significant negative correlation between promoter methylation levels and HSF1 expression in CRC, providing an additional reason for HSF1 overexpression.
Protein expression is a complex process that is controlled by multiple steps, involving the transcription rate, mRNA degradation, translation regulation and protein degradation [49]. Bioinformatic analysis revealed a high translation efficiency of HSF1 mRNA in HCT116 cells. Our recent research has demonstrated that β-catenin inhibit miR-455-3p to increase m6A modification of HSF1 mRNA and facilitate its translation in CRC [38]. But there was still one problem unsolved that the expression of exogenous HSF1 without the 3'-UTR was also diminished after knockdown of β-catenin. Notably, in this study, we further found that this process was dependent on HuR. HuR binds transcripts primarily at U-or AU-rich RNA stretches which are usually located in the 3'-UTR of the target mRNA. Our results confirmed that besides 3 -UTR, HuR could also directly bind to the CDS region of HSF1 mRNA to promote its translation, providing new regulation mode for HuR. Previous studies showed that the effect of HuR was dependent on HSF1 expression [50]. These data supported a possibility that HSF1/HuR may constitute a feed-forward loop to coordinate their functions in CRC. Additionally, we identified a physical interaction between β-catenin and HuR, which has been reported previously [51,52]. More importantly, we revealed that knockdown of β-catenin attenuated the binding levels between HuR and HSF1 without influencing the expression of HuR, suggesting that HuR was able to stimulate the translation of HSF1 mRNA under the activity of βcatenin.
Collectively, our results shed light on the roles and driving forces for HSF1 overexpression in CRC. However, several issues still exist. First, the variation of some results could be due to differences in sample sizes, which is a limitation of online databases. Additionally, the specific role of β-catenin-HuR-HSF1 axis in CRC progression needs to be further confirmed by in vivo experiments.

Conclusions
HSF1 may be useful as a diagnostic and prognostic biomarker for CRC. HSF1 was closely correlated with immune features. Genetic and epigenetic alterations contributed to HSF1 overexpression in CRC. More importantly, we demonstrated that HSF1 may be regulated at the level of mRNA translation by β-catenin-induced HuR activity.

Highlights
• HSF1 may serve as a diagnostic and prognostic biomarker for CRC. • HSF1 was closely correlated with tumor immunity in CRC. • HSF1 overexpression in CRC may be attributed to genetic and epigenetic alterations. • HSF1 may be regulated at the level of mRNA translation by β-catenin-induced HuR activity.