Prognostic and Functional Significant of Heat Shock Proteins (HSPs) in Breast Cancer Unveiled by Multi-Omics Approaches

Simple Summary In this study, we investigated the expression pattern and prognostic significance of the heat shock proteins (HSPs) family members in breast cancer (BC) by using several bioinformatics tools and proteomics investigations. Our results demonstrated that, collectively, HSPs were deregulated in BC, acting as both oncogene and onco-suppressor genes. In particular, two different HSP-clusters were significantly associated with a poor or good prognosis. Interestingly, the HSPs deregulation impacted gene expression and miRNAs regulation that, in turn, affected important biological pathways involved in cell cycle, DNA replication, and receptors-mediated signaling. Finally, the proteomic identification of several HSPs members and isoforms revealed much more complexity of HSPs roles in BC and showed that their expression is quite variable among patients. In conclusion, we elaborated two panels of HSPs that could be further explored as potential biomarkers for BC progression and prognosis. Abstract Heat shock proteins (HSPs) are a well-characterized molecular chaperones protein family, classified into six major families, according to their molecular size. A wide range of tumors have been shown to express atypical levels of one or more HSPs, suggesting that they could be used as biomarkers. However, the collective role and the possible coordination of HSP members, as well as the prognostic significance and the functional implications of their deregulated expression in breast cancer (BC) are poorly investigated. Here, we used a systematic multi-omics approach to assess the HSPs expression, the prognostic value, and the underlying mechanisms of tumorigenesis in BC. By using data mining, we showed that several HSPs were deregulated in BC and significantly correlated with a poor or good prognosis. Functional network analysis of HSPs co-expressed genes and miRNAs highlighted their regulatory effects on several biological pathways involved in cancer progression. In particular, these pathways concerned cell cycle and DNA replication for the HSPs co-expressed genes, and miRNAs up-regulated in poor prognosis and Epithelial to Mesenchymal Transition (ETM), as well as receptors-mediated signaling for the HSPs co-expressed genes up-regulated in good prognosis. Furthermore, the proteomic expression of HSPs in a large sample-set of breast cancer tissues revealed much more complexity in their roles in BC and showed that their expression is quite variable among patients and confined into different cellular compartments. In conclusion, integrative analysis of multi-omics data revealed the distinct impact of several HSPs members in BC progression and indicate that collectively they could be useful as biomarkers and therapeutic targets for BC management.


Introduction
Breast cancer (BC) represents the most common type and the leading cause of death of cancer among females and accounts for~16% of all cancers [1]. The etiology of BC is quite complex, involving several genetic and epigenetic changes [2]. Consequently, breast cancer is a heterogeneous disease with several subtypes of different cellular compositions, molecular alterations, as well as clinical behavior [3,4]. Several factors, such as histological grade, type and size of the tumor, lymph node metastasis, estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER), are generally considered as prognostic factors, but are insufficient to provide useful information for clinical management [5,6]. Despite the achieved improvements, the prognosis of BC patients is still a poor predictor and the identification of more reliable biomarkers must be explored.
Heat shock proteins (HSPs) are one of the largest groups of molecular chaperones that assist the correct folding of partially folded or denatured proteins and prevent the formation of potential aggregates in the cells, promoting their proteasomal degradation [7]. HSPs were first discovered as stress-inducible proteins against physical (temperature elevation) or chemical (increase or decrease in pH, salinity, or oxygen concentration) stressors, performing a wide range of functional activities: Modulation of their synthesis, regulation of kinases activation, participation in signal transduction pathways and rRNA processing [8]. Under physiological conditions, HSPs perform chaperone functions on a broad array of client proteins, including transcription factors, nuclear hormone receptors, viral proteins and signaling mediators. Depending on their molecular weight and main functions, six major subfamilies of HSPs have been reported: HSP110, HSP90, HSP70, HSP60 (chaperonins), HSP40 (DNAJ-class proteins), and HSP20 (small heat shock proteins) [9]. Within each subfamily, some members are constitutively expressed, inducible regulated, and/or targeted to different compartments [7,10,11]. Genes encoding HSPs are also transcriptionally regulated by a variety of physiologic processes not typically associated with cell stress, including cell cycle, cell proliferation, and differentiation [12].
Aberrant expression of HSPs has been reported in a wide range of human tumors, including breast, endometrial, ovarian, colon, lung, and prostate [13,14]. In cancer cells, the HSPs network is extensively remodeled and could participate in the functional metabolism of tumor cells, protect them from harmful factors, provide an immunogenic context and allow tumor cells to tolerate genetic alterations, which would otherwise be fatal [15][16][17]. The expression of HSPs has been associated with tumor cell proliferation and differentiation, as well as with resistance to apoptosis and poor prognosis [18]. These observations, while intriguing, have not resolved whether the association between cancer evolution and HSPs is causal or correlative.
Here, we aimed to more precisely estimate the relationship between HSP expression and prognosis in breast cancer, and highlight the functional implications of their deregulation by using a multi-omics integrative analysis.
An in silico analysis on differential expression of HSP family members in breast and other tumors was firstly performed, and the differentially expressed members (50 upregulated and 26 down-regulated) were checked for genetic alterations and promoter methylation status. The prognostic value of HSP members and association with clinical variables were also assessed. Next, to highlight the important biological pathways deregulated as a consequence of HSPs expression, the lists of HSPs co-expressed genes and miRNAs extracted by different cross-platforms were used to perform functional analyses. Finally, proteomics investigations on a large sample set of breast cancer tissues showed additional complexity of the HSPs network in BC, with different protein isoforms identified for each HSP gene, probably performing additional functions and/or interactions. Our findings identified two distinct HSP-clusters with opposite functions in breast cancer: The first one with 19 members, was significantly associated with a poor prognosis and affected both cell cycle and DNA replication. The second HSP-cluster included 10 members associated with good prognosis, presumably affecting receptor signaling and Epithelial to Mesenchymal Transition (ETM). Interestingly, the HSPs-network involved also several miRNAs and protein isoforms that, in turn, could add new perspectives on HSPs contribution to BC. To the best of our knowledge, this is the first attempt exploring the potential mechanisms of the HSPs-mRNAs-miRNAs-proteins network, underlining that HSP-based regulation of BC progression is a multi-level process. We propose that the identified HSP-clusters can be used to stratify patients according to their predicted clinical outcome.

Expression Analysis of HSP Members Using UALCAN and Oncomine
The mRNA expression levels of each HSP member in tumors and their normal tissue counterparts were analyzed using UALCAN (https://ualcan.path.uab.edu/, accessed date 26 February 2021) and Oncomine database (https://www.oncomine.org/resource/ login.html, accessed date 26 February 2021). UALCAN is a comprehensive, user-friendly, and interactive web resource for analyzing cancer OMICS data (TCGA, MET500, and CPTAC) [19]. For each query, the database provides graphs and plots depicting the expression profile of protein and miRNAs between normal and cancer tissues and evaluates epigenetic regulation of gene expression by promoter methylation. The expression level of HSP members was normalized as transcript per million reads, and a p-value of no more than 0.01 calculated through the Student's t-test was considered to be significant. Oncomine, an online microarray database, can analyze the mRNA expression differences between tumor and normal tissues in 20 different human cancers. The thresholds were set as follows: p-value: 0.01; fold change: 2; gene rank: 10%; analysis type: cancer vs. normal analysis; data type: mRNA.

Genetic Alterations and Epigenetic Regulation of HSP Members Using cBioPortal and UALCAN
The frequency of HSPs alterations (amplification, deep deletion, and missense mutations) in BC patients was assessed using the OncoPrint tool of cBioPortal (http://www. cbioportal.org, accessed date 26 February 2021). cBioportal is an interactive open-source platform, that provides large-scale cancer genomics data sets database [20]. Promoter methylation status was analyzed using UALCAN.

Evaluation of the Prognostic Value of HSP Members Using KM Plotter Database and Human Protein Atlas
The correlation between the mRNA expression levels of HSP members and the survival probability of BC patients was analyzed using the Kaplan-Meier Plotter database (http://kmplot.com/analysis/, accessed date 26 February 2021) [21]. For each HSP gene, Distant Metastasis Free Survival (DMFS), Overall Survival (OS), and Relapse Free Survival were assessed. The Affimetrix ID probes were entered into the database by using the multigene classifier. Patients were divided into high and low expression groups by using the best cutoff of mRNA expression. The database provided statistical plots for individual genes. The correlation between HSPs protein expression levels via immunohistochemistry analysis and overall survival of BC patients was analyzed using the Human Protein Atlas (https://www.proteinatlas.org/, accessed date 26 February 2021) [22]. For each HSP protein, the immunohistochemical staining was evaluated and scored based on staining intensity (negative, weak, moderate, or strong) and the fraction of stained cells (<25%, 25~75%, >75%). The staining quantity of each protein via IHC was determined as the percentage of stained cells in 10 high-power fields. All annotation data and immunohistochemistry images analyzed in the present investigation and all anti-body validation data are publicly available at https://www.proteinatlas.org/, accessed date 26 February 2021.

Relationship between HSP Members and Clinical-Pathological Parameters Using Bc-GenExMiner and GOBO
The correlation between the mRNA expression of HSP members and different clinicalpathological parameters, such as ER/PR/HER receptor status and lymph node (N) involvement, were evaluated using the Breast cancer Gene-Expression Miner v4.5 database (bc-GenExMiner v4.5) [23], an on-line statistical mining tool (http://bcgenex.ico.unicancer.fr, accessed date 26 February 2021) of published annotated BC transcriptomic data (DNA microarrays [n = 10,716] and RNA-seq [n = 4712]). bcGenExMiner offers the possibility to explore gene expression, prognostic and correlation analyses, providing various kinds of plots. The association between HSP members and grading (G1/G2/G3) was performed by Gene expression-based Outcome for Breast cancer Online (GOBO database) [24]. GOBO (http://co.bmc.lu.se/gobo, accessed date 26 February 2021) enables a rapid assessment of gene expression levels, the identification of co-expressed genes and association with the outcome for single genes, gene sets, or gene signatures in an 1881-sample breast cancer data set, generated on Affymetrix U133A microarrays.

Analysis of HSPs Co-Expressed Genes and Pathways Enrichment Analysis
The co-expression profiles of the HSPs associated with poor and good prognosis were retrieved using three different online platforms, namely UALCAN (https://ualc an.path.uab.edu/, accessed date 26 February 2021), GOBO (http://co.bmc.lu.se/gobo/ coexpressed_genes.pl, accessed date 26 February 2021) and bcGenExMiner (http://bcgenex. centregauducheau.fr/BC-GEM/GEM-Requete.php?mode=5, accessed date 26 February 2021). Results were statistically analyzed using Pearson's correlation coefficient (p ≥ 0.4 and p ≤ −0.4). Positive and negative HSPs co-expressed genes were properly sorted to obtain two lists of up-regulated genes in poor prognosis (positive co-expressed genes for HSPs correlated with poor prognosis and negative co-expressed genes for HSPs correlated with good prognosis) and up-regulated in good prognosis (positive co-expressed genes for HSPs correlated with good prognosis and negative co-expressed genes for HSPs correlated with poor prognosis). Each gene list was submitted into the FunRich database [25], (http://www.funrich.org/, accessed date 26 February 2021) a stand-alone software tool used mainly for functional enrichment and interaction network analysis of genes and miRNA. Enriched biological pathways were ranked by p-value and the top six significant pathways were exhibited as bar charts. A p-value < 0.05 was regarded as statistically significant.

Interactome Construction Using STRING Database
STRING (https://string-db.org, accessed date 26 February 2021) is a database of known and predicted protein-protein interactions. The interactions include direct (physical) and indirect (functional) associations derived from computational prediction, knowledge transfer between organism, and interactions aggregated from other (primary) databases analysis [26].

HSP-Mediated Regulatory Network of miRNAs-mRNAs
The list of HSPs co-expressed miRNAs was retrieved using bcGenExMiner. Pearson's correlation coefficient was considered significant when p ≥ 0.4 and p ≤ −0.4. Positive and negative HSPs co-expressed miRNAs were properly sorted to obtain two lists of upregulated miRNAs in a poor and in good prognosis. Within each list, only the miRNAs recurring at least three times were selected for miRNAs-mRNAs network construction using miRNet (www.mirnet.ca, accessed date 26 February 2021), an easy-to-use web-based tool designed for creation, customization, visual exploration, and functional interpretation of miRNA-target interaction networks [27]. miRNet was also used to predict miRNAs acting on HSP-cluster associated with poor and good prognosis.

Identification of HSP Proteins and Isoforms in Breast Cancer Tissues
Two-dimensional (2D)-IPG electrophoresis was performed on 63 protein extracts from breast cancer tissues obtained following surgical interventions during the years 2004-2010 at the "La Maddalena" Hospital of Palermo, as previously described [28][29][30][31][32][33]. Briefly, the surgical samples homogenated overnight at 4 • C with RIPA buffer, were centrifugated, and the supernatants were dialyzed against ultrapure distilled water, lyophilized and resuspended in ISOT buffer. Aliquots of 45 µg of total proteins were rehydrated in rehydration buffer containing 8 M urea, 2% CHAPS, 10 mM DTE, and 0.5% carrier ampholytes (Resolyte 3.5-10). The electrophoretic separation was performed on 18 cm long strips with a pH range of 3-10. The strips were then equilibrated in a solution containing 50 mM Tris-HCl pH 6.8, 6 M urea, 0.5% SDS, 30% Glycerol, 130 mM DTE and 135 mM Iodoacetamide and then separated on 9-16% linear gradient polyacrylamide gels (SDS-PAGE), with a constant current of 20 mA/gel, as previously described [34][35][36]. The gels were silver stained and analyzed with the dedicated Image-Master 2D Platinum software. Proteins of interest were excised from the gel and the identity was assigned by peptide mass fingerprinting using the Voyager DE MALDI-TOF mass spectrometer as described [37][38][39][40][41][42]. The expression level of protein spots was calculated as the volume of the spots (i.e., integration of optical density over the spot area), relative to the sum of the volume of all spots on each gel (% Vol). Measurements of relative expression levels of individual protein spots were normalized in each proteomic map for actin content (N % V), as previously reported [43].

Results
We systematically analyzed by data mining 95 heat shock genes (and the protein members that they encode), grouped into 6 sub-families (Table 1)

Gene Expression Analysis of HSP Family Members between Normal and Cancer Tissues
The expression profiles of HSP family members in breast cancer were determined using UALCAN database. Among the 95 HSP members, 76 were significantly deregulated in cancer tissues when compared to the corresponding normal tissues. In particular, 50 HSPs, were found up-regulated ( Figure 1A, Table 2) while 26 HSPs, were found downregulated in breast cancer than normal tissues ( Figure 1B, Table 2). Interestingly, the number of transcripts, per million, was spanning between 0.07 for DNAJA5G and 3500 for HSPB1 and CRYAB, suggesting that some members are expressed at very low levels (HSPA1L, HSPB9, DNAJB7, DNAJ5B, DNAJC6, DNAJC27, and DNAJC28) while others are expressed at high levels (HSPA1A, HSPA8, HSPB6, HSP90AA1, HSP90AB1, and HSP90B1). Next, Oncomine meta-analysis was employed to assess the comprehensive expression of HSP members across 20 cancer types ( Figure S1). The database compared the HSPs gene expression levels in cancer versus normal samples across a wide variety of datasets in different cancer types. A total of 1460 analyses showed a significant statistical difference for mRNA expression in cancer versus normal samples (the selected thresholds for the analyses were: p-value: 1 × 10 −4 ; fold change: 2; gene rank: top 10%; data type: mRNA; sample type: clinical specimens). Overlapping results were recorded between UALCAN and Oncomine platforms except for HSPA1A, HSPA1L, HSPA2, DNAJA4, DNAJB2, DNAJC4, DNAJC11, DNAJC21 and GAK found up-regulated in breast cancer but collectively downregulated in the other tumors, and DNAJC8, DNAJC24, SACS and HSPB6 down-regulated in breast cancer and collectively up-regulated in other tumors. The obtained results clearly indicated that the HSP members are involved in carcinogenesis, acting both as oncogenes or suppressor genes, and suggested that some of them could exert specific roles in different cancer types. Moreover, a high agreement was obtained using different bioinformatics platforms, supporting the robustness of the results.

Genetic and Epigenetic Alterations of HSP Members
It is well-recognized that the expression levels of selected genes could depend on genetic (amplifications or deletions) and epigenetic (promotor methylation status) alterations. To verify whether the deregulated expression of HSP members in breast cancer could be caused by one or a combination of these factors, the genetic mutations of HSP members, in a cohort of breast cancer patients using cBioPortal web, were analyzed. The mutation rate was higher for the HSP members up-regulated in BC (Figure 2A), spanning from 0.4% to 12%. Among these, TRAP1, HSPA6, CCT2, CCT3, DNAJA3, DNAJC5, DNAJC5B and DNAJC19 showed a percentage of alterations higher than 3%. A lower mutation rate was recorded for HSP members down-regulated in BC ( Figure 2B), spanning from 0.2% to 4%. Among these, only DNAJB13 showed a percentage of alteration higher than 3%. Concerning the type of HSPs alteration, amplification was the major one, followed by deletion. A very low percentage of missense mutations was detected in all the analyzed samples. Promoter methylation status was analyzed from TCGA data through UALCAN database. Among the differentially expressed HSPs, 26 over 50 up-regulated ( Figure 3A) and 17 over 26 down-regulated ( Figure 3B), showed statistical differences in promoter methylation status between normal and breast cancer samples. Different beta value cut-off indicates hyper-methylation [beta value: 0.7-0.5] or hypo-methylation [beta-value: 0.3-0.25]. Interestingly, based on beta values, 23/26 HSPs were found with the hypo-methylated promoter and 9/17 HSPs were validated with the hyper-methylated promoter. According to expression pattern analysis among the HSPs up-regulated in BC, 15/26 showed lower promoter methylation while among the HSPs down-regulated in BC, 10/17 showed higher promoter methylation in tumor samples compared with normal samples. These data suggested that the deregulated HSPs gene expression associated with BC may be partly due to genetic and epigenetic alterations.

HSPs Expression and Clinical Outcome
To analyze the prognostic value of the HSP members, Kaplan Meier-plotter database restricted to BC was searched. For each HSP gene (using the Affimetrix ID probe listed in Table 1), the survival outcomes were evaluated as Relapse Free Survival (RSF), Distant Metastasis Free Survival (DMFS), and Overall Survival (OS), (Table S1). Patients were divided into two groups (high and low) based on the expression levels of individual HSPs by selecting the best cut-off. In this case, the software computed all possible cut-off values between the lower and the upper quartiles and the best performing threshold was used as the cut-off value. HSPs expression levels were considered as significantly associated with prognosis when all survival outcomes (RSF, DMSF, and OS) were coherently significant (Table 3, Figure S2). The increased expression levels of 21 HSPs were significantly associated with a worse prognosis, whereas high transcriptional levels of 13 HSPs favored a good prognosis. Among these, only HSPA1B, DNAJB1, DNAJB8, DNAJB12, and DNAJC16 were not differentially expressed in BC, and thus, not further considered. Interestingly, when the survival analyses were re-evaluated, including the panel of worse prognosis-associated HSPs, better prognostic values than the individual ones were obtained ( Figure 4A). Similar results ( Figure 4B) were obtained including the panel of good prognosis-associated HSPs. These results revealed that collectively, the selected clusters of 19 HSPs and 10 HSPs could be represent robust biomarkers for poor and good prognosis in breast cancer, than individual members.  Table 3. List of HSP members significantly associated with clinical outcome, including Relapse Free Survival (RFS), Overall Survival (OS), and Distant Metastasis Free Survival (DMFS), evaluated by KM Plotter database. The Log Rank value was marked in red when higher expression of selected HSP was associated with poor prognosis and in the black when higher expression of selected HSP was associated with a good prognosis. The probe expression value cut-off for each analysis is also reported.

Relationship between HSP Members and Clinical-Pathological Parameters
The relationship between the expression levels of HSP-clusters differentially associated with prognosis and current clinical-pathological parameters was assessed using bcGeneXMiner and GOBO databases (Table 4). Clinical parameters included immunocytochemical expression of estrogen receptor (ER), progesterone receptor (ER), human epidermal growth receptor 2 (HER), lymph node (N) metastases, and grading (G1-G2-G3). The HSP-cluster associated with poor prognosis showed significant correlations with ER-/PR-, HER+, N+, and G3 tumors, which clinically identify more aggressive tumors. On the contrary, for the HSP-cluster associated with a good prognosis, an inverse trend of correlation with ER+/PR+, HER-N-and G1 tumors was found. Accordingly, it is well known that HER-, N-and G1 tumors are less aggressive. The obtained results confirm that the identified clusters can group patients with homogeneous clinical characteristics and underline their robustness in predicting prognosis in subgroups of patients. Table 4. Relationship between HSPs mRNA expression in BC and clinical-pathological parameters, including estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth receptor 2 (HER) status, lymph node (N) involvement, and grading (G1-G2-G3). Data were retrieved for each HSP using bcGeneXMiner and the grading was assessed using GOBO database. p < 0.05 was considered significant (p < 0.05 *, p < 0.01 **, p < 0.001 ***).

HSPs-Gene Co-Expression Networks
To unveil the mechanisms by which HSPs deregulation affect breast cancer development and influence prognosis, we identified the predicted genes that are transcriptionally co-expressed with individual HSPs via computational analysis. The co-expressed genes for the HSP-cluster that are significantly associated with a poor and good prognosis were queried to GOBO, UALCAN, and bc-GenExMiner database, respectively, using a Pearson correlation score of ≥ 0.4. The lists of the derived positive and negative HSPs co-expressed genes (Table S2), were appropriately sorted to obtain up-regulated genes in poor prognosis and up-regulated genes in good prognosis. Functional enrichments in the HSP-regulated networks (biological pathways) were highlighted by using the Fun Rich tool. Although each HSP was significantly associated with different genes in different databases, they were implicated in similar biological functions. Overall, biological processes were included in the mitotic cell cycle, DNA replication and cell cycle checkpoint for the co-expressed HSP genes up-regulated in poor prognosis, independently from the used database ( Figure 5A, Figure S3A), suggesting that these HSPs-regulated biological processes might play a significant role in the initiation and progression of breast cancer. Less robust results were obtained for the co-expressed HSP genes up-regulated in good prognosis ( Figure 5B, Figure S3B): the biological processes in which these genes are involved were Epithelial to Mesenchymal Transition (ETM) for the gene list derived from GOBO and bc-GenExMiner, and receptor signaling network (PAR-1, c-MET, VEGFR1, VEGFR2, HER2) for the gene list derived from UALCAN. Interestingly, when the list of common genes from the three bioinformatics platforms ( Figure 5C,D) was used to find functional protein-protein interactions by using the STRING analysis tool ( Figure 5E,F), a significant network, involved in cell cycle and DNA replication, was obtained for the co-expressed HSP genes up-regulated in poor prognosis.
Again, the protein-protein interaction network derived from the co-expressed HSP genes up-regulated in good prognosis, was less robust and affected the Estrogen-dependent gene expression. Collectively, the pathway analysis of HSP co-expressed genes showed a specific signature for poor prognosis, affecting cell cycle and DNA replication.  . (A,B) The top six biological pathways in which the HSPs co-expressed genes up-regulated in poor and good prognosis, derived from GOBO, were significantly involved. (C,D) Venn diagrams of overlapping HSPs co-expressed genes up-regulated in poor and good prognosis, derived from GOBO, UALCAN, and bc-GenExMiner. (E,F) Protein-protein interaction network of the common HSPs co-expressed genes derived from the three cross-platforms GOBO, UALCAN and bc-GenExMiner using STRING database. The network was performed with medium confidence (combined score > 0.4) and disconnected nodes were deleted.

Regulatory Network Analysis of HSPs-miRNA-mRNA
To further explore the molecular mechanisms responsible for HSPs effects on BC, we analyzed the networks of HSPs co-expressed miRNAs generated by the bc-GenExMiner database. Positive and negative HSPs co-expressed miRNAs were properly sorted to obtain two lists of up-regulated miRNAs in poor and in good prognosis. Only the miRNAs repeated at least three times were selected for miRNA-mRNA network construction using miRNet. A total of 14 miRNAs were clustered within the miRNAs network up-regulated in poor prognosis (Table 5). Interestingly, the target genes were involved in the cell cycle ( Figure S4A). A total of 20 miRNAs were clustered within the miRNAs network up-regulated in good prognosis ( Table 5). The target genes were implicated in the negative regulation of DNA-dependent transcription ( Figure S4B). Interestingly, among the HSPs coexpressed miRNAs up-regulated in poor prognosis 5 miRNAs (hsa-mir-320b-2, hsa-mir-545, hsa-mir-651, hsa-mir-570, hsa-mir-643) were significantly up-regulated in breast cancer than normal tissues, as verified by using UALCAN, while among the HSPs associated miRNAs up-regulated in good prognosis, 10 miRNAs (hsa-mir-143, hsa-mir-196b, hsa-mir-145, hsamir-150, hsa-mir-1228, hsa-mir-1910, hsa-mir-1247, hsa-mir-411, hsa-mir-370, hsa-mir-770) were significantly down-regulated in BC compared to normal tissues. As the mRNAs-miRNAs network was created considering the HSPs clusters associated with prognosis, a greater interconnection mRNAs-miRNAs was recorded for the cluster associated with poor prognosis (Figure S4C), compared to the network obtained for HSPs-cluster associated with good prognosis ( Figure S4D). As expected, collectively they were involved in protein folding. Several miRNAs, selected as key components of the networks, showed significant differences between breast cancer than normal tissues. In particular, 12 of them were up-regulated in breast cancer tissues, while 6 of them were down-regulated in breast cancer. These results, for the first time, indicated that the HSPs deregulation in BC affected epigenetic networks involving also several microRNAs, pointing that the complex network triggered by HSPs deregulation provides a distinctive molecular portrait of each tumor. Table 5. List of the hubs miRNAs involved in the HSPs-miRNAs-mRNAs network. The red bolded miRNAs were significantly up-regulated in breast cancer than normal tissues; the black bolded miRNAs were significantly down-regulated in breast cancer than normal tissues. A p-value < 0.01 was regarded as significant.

Proteomics Expression of HSP Members and Prognostic Significance
Complementing the transcriptomics results with proteomics ones allows studying globally the HSPs-network within cancer cells and is also necessary to validate the existence of different protein isoforms as well as the exact localization of proteins, tightly linked to their function. The prognostic value of HSPs protein expression, evaluated by immunohistochemistry (IHC) staining, was retrieved using The Human Protein Atlas (HPA) database ( Figure 6A). High expression levels of 5 HSP members (HSPA9, HSP90AA1, TCP1, CCT4 and CCT6A) were significant associated with shorter overall survival, while a high expression level of HSPA2 was associated with a good prognosis. Moreover, taking advantage of our previous proteomic data, we performed a comprehensive screening of protein spots, in different proteomic maps, and collectively the following HSP protein members were identified: HYOU1 (HYOU1), HSP90AA1 (HS90A), HSP90AB (HS90B), HSP90B1 (ENPL), HSPA1A (HS71A), HSPA4 (HSP74, 2 isoforms), HSPA5 (BIP), HSPA8 (HSP7C, 4 isoforms), HSPA9 (GRP75), HSPD1 (CH60, 4 isoforms), TCP1 (TCPA), CCT2 (TCPB), CCT3 (TCPG), CCT5 (TCPE), CCT6A (TCPZ), HSPB1 (HSP27, 7 isoforms), HSPE1 (CH10, 2 isoforms). Figure 6B shows a representative proteomic map of a breast surgical tissue silver-stained, with cropped windows where the 17 identified HSPs proteins and 14 isoforms are marked with labels corresponding to the gene name. The protein identity was assessed by Maldi-Tof (Table S3). Different isoforms of the same protein were labeled by alphabetic letters starting from the more acidic one. Several HSP proteins and isoforms showed high variability of expression between analyzed patients ( Figure 6C), suggesting that different proteins and isoforms should be used for patient stratification. Finally, the immunohistochemical analysis using the HPA database gives us the possibility to define HSPs protein localization in different compartments at a single-cell and subcellular level and also provided important spatial information in the context of neighboring cells. Figure 6D reports immunohistochemical staining of BC tumors showing weak, moderate and strong expression of selected HSPs, based on the major number of identified isoforms in proteomics maps. Interestingly, antibody immunoreactivity was observed in several cellular compartments, such as cytoplasm, plasma membrane, and nucleus. Collectively, our proteomic data on HSPs expression and prognostic value are in line with genomic data.
Although, proteomic investigations showed a more complex scenario in which different HSP-isoforms can operate within the cancer cells, creating new functional arrangements, difficult to predict based on gene expression. In the hope for the clinical use of HSP-clusters, both the gene and protein expression levels should be taken into consideration, for a more accurate prognosis.

Discussion
HSPs are a ubiquitous and highly conserved protein class, essential for cell protection from a wide range of harmful conditions, including heat shock, oxidative stress, inflammation, and proteotoxic stress (a specific form of deleterious stress causing increased levels of misfolded proteins). However, when the cells reach the limit of stress tolerance, they activate apoptosis or autophagy. HSPs play critical roles in inhibiting pro-apoptotic molecules through the modulation of several signaling cascades such as JNK, AKT, and NF-κB [44]. Several HSP members are closely associated with the onset and progression of several human cancer types [45][46][47][48][49], but the expression pattern and potential roles of HSPs in breast cancer, as well as their underlying regulatory mechanisms, remain poorly investigated. Only two recent studies investigated the prognostic significance of HSP members in BC, exclusively limited to the transcriptomic analysis [50,51]. Here, for the first time, we performed an integrated multi-omics analysis to assess the prognostic value and the functional implications of HSPs deregulation in BC. Using UALCAN database, we found that collectively HSP family members were deregulated in BC. In particular, 50 HSPs members were significantly up-regulated in cancer compared to normal tissues, while 26 HSPs members were down-regulated. These results indicated that HSPs could act both as oncogenes or tumor suppressor genes. Moreover, the Oncomine analysis of HSPs deregulation in different cancer types showed that HSPs exert both, pro-and anti-tumorigenic actions depending on the tumor type.
HSP members are coded by distinct chromosomal loci that are not clustered. The molecular basis of HSPs deregulation in BC has not been well-studied and may have multiple molecular etiologies. For example, associated with genetic alterations and/or epigenetic modifications. Collectively, we showed that the mutation rate (amplification followed by deletion) was higher in HSP members up-regulated in BC, while a lower mutation rate was detected in HSP members down-regulated in BC. Besides, HSPs promoter regions were found to be commonly hypo-methylated in HSP members up-regulated in BC, while a general tendency to hyper-methylated promoters was detected in HSP members down-regulated in BC. In alignment with our results, the CpG island located in the promoter regions of the DNAJC10 gene was found to be frequently hypo-methylated in breast cell lines [52]. We speculate that genetic alterations and epigenetic deregulations could in part explain the deregulation of HSPs in BC. Other mechanisms responsible for HSPs deregulation are: (1) The heat shock factor (HSF) family that ensures prompt transcriptional activation of HSP members under stress and equally precipitous switch-off after recovery [53]; (2) the genetic changes associated with tumor progression, producing over-expressed onco-proteins, mostly mutated and/or conformationally altered proteins, may elicit an HSPs response; (3) the aneuploidy (defined as abnormal chromosome number), and its associated abnormalities [54], such as chromosomal instability (CIN), could impair HSPs expression and function. Aneuploidy alters the relative dosage of genes on the affected chromosomes, leading to cellular stress response at least partially due to impaired protein folding capacity. This, in turn, affects the cellular protein quality control pathways that preserve homeostasis and induces chronic proteotoxic stress (broadly referring to the overburdening of cellular systems that maintain proper protein folding and homeostasis) [55,56]. Consequently, the aneuploid cells would exhibit an over-production of proteins, relative to the chaperone systems needed to fold nascent polypeptides or the degradation systems that remove misfolded or damaged proteins. This is of particular importance regarding the maintenance of the correct stoichiometry of macromolecular complexes, whose subunits may be encoded by genes on different chromosomes. Further studies are needed to more comprehensively explore the detailed molecular mechanisms of altered HSPs expression in BC.
The prognostic analysis, performed on HSPs gene expression levels and clinical outcome, identified two HSP-clusters, associated with a poor and good prognosis, respectively, in line with previous reports [50,51]. Analytically, when the HSP subfamilies were considered, the chaperonins subfamily (TCP1, CCT2, CCT3, CCT4, CCT6A, CCT7, and CCT8) was collectively associated with poor prognosis. Accordingly, the expression levels of different CCT members are reported as up-regulated in various cancers [57]. In particular, it was estimated that CCT members facilitate the folding of approximately 10% of newly synthesized proteins [58], and in cancer cells they could potentially fold more proteins, with substrates, including oncogenic proteins and mediators of oncogenesis [59]. A more heterogeneous involvement with prognosis was detected for HSP40 members, with 7 members significantly associated with poor prognosis (DNAJA1, DNAJB1, DNAJB11, DNAJC2, DNAJC5, DNAJC9, DNAJC17) and 10 members associated with good prognosis (DNAJB8, DNAJB12, DNAJB14, DNAJC4, DNAJC5G, DNAJC10, DNAJC12, DNAJC14, DNAJB16, and DNAJC19). HSP40 members acting as co-chaperones regulate the major functions of other HSPs (HSP90, HSP70, and HSP60) and could be able to switch the canonical functions of HSPs towards oncogenic roles. HSP40 members work as tumor suppressors and inhibit tumor growth [60]. Other members have been implicated in cancer development and metastasis, such as DNAJA1 (in glioblastoma and prostate), DNAJB11 (in ovarian tumors), and DNAJC9 (cervical) [61].
Interestingly, the HSP-cluster associated with poor prognosis showed significant correlations with ER-/PR-/HER+, lymph node metastases, and higher grading, identifying more aggressive tumors. On the contrary, the HSPs cluster associated with a good prognosis showed a significant correlation with ER+/PR+, HER-, lymph node-negative, and lower grading tumors. Collectively, these results confirm the robustness of the identified HSP clusters in predicting prognosis in clinical subgroups of patients.
To explore HSPs-related pathways altered in BC, genes co-expressed with HSPs were analyzed. We found that HSPs co-expressed genes associated with poor prognosis were enriched with cell cycle, DNA replication and cell cycle checkpoint, commonly oncogenic processes associated with cell proliferation. Therefore, we hypothesize that HSP-networks, affecting cell cycle, and DNA replication, may provide a robust measure of proliferation rate and more aggressive clinical course in BC. Less robust results were obtained for the HSPs co-expressed genes up-regulated in a good prognosis. The pathways enrichment analysis of these genes showed that they are involved in ETM, receptors signaling (PAR-1, c-MET, VEGFR1, VEGFR2, HER2) and Estrogen-dependent gene expression and thus further functional investigations are needed to validate these findings.
As proteins are considered the molecular effectors, the characterization of HSPs protein expression profiles represents the best chance of better understanding their molecular context-dependent functions. Proteomic data on HSPs expression and prognostic value are in line with the transcriptomic results, although proteomic investigations showed a more complex scenario in which different HSP-proteoforms (probably arising from alternatively spliced RNA transcripts and post-translational modifications (PTM)) can operate within the cancer cells, creating new functional arrangements, difficult to predict on gene expression. Based on IHC analyses, only 5 HSP members (HSPA9, HSP90AA1, TCP1, CCT4 and CCT6A) were significant associated with shorter overall survival, while a higher expression level of HSPA2 was associated with a good prognosis. This stringency in the prognostic value of HSP proteins could partly depend on the expression of multiple proteoforms, probably localized in different subcellular compartments, affecting a variety of intracellular signaling pathways. It is well-known that the HSPs activity is regulated by PTM such as phosphorylation, acetylation, methylation, ubiquitination, glycosylation, S-nitrosylation [71]. Since the addition and reciprocal removal of chemical groups can be triggered very rapidly, PTMs provide an efficient switch to precisely regulate the HSPs functions and activities. For instance, the HSP27 (HSPB1) function is modulated by phosphorylation at different serine residues, catalyzed by the mitogen-activated protein kinase MAPK2 and 3 [72]. Moreover, in clinical cancer tissues, various phosphorylation patterns of HSP27 have been found to associate with the aggressiveness of tumor phenotype [73]. Consequently, it has been proposed that HSP27 phosphorylation isoforms could also represent useful tumor markers, while a more comprehensive analysis of tissue samples will be required. Besides, each HSP is properly localized in a specific cellular compartment(s). For instance, HSP70 and HSP90 can be localized in the cytosol and the nucleus, grp78 (HSPA5) in the endoplasmic reticulum. Whereas, HSP60 (HSPD1) is found in mitochondria [74]. Some HSPs can also be found on the cell surface, such as HSP60 and HSP70, especially in lipid rafts (plasma membrane subdomains, containing high levels of cholesterol and glycosphingolipids) [75]. More recently, HSPs are considered key players in the intercellular cross-talk [76]. They may be secreted through Extracellular Vesicles (EVs), e.g., exosomes, which derive from endosomes and multivesicular bodies. HSP27, HSP60, HSC70, HSP70, and HSP90 were found in the extracellular environment where their role is believed to be immunogenic and may induce either a pro-or anti-inflammatory response. For example, EV containing HSP70 on their surface activates macrophages and natural killer cells. HSPs can be released from cells in free, soluble form, possibly via Golgi. In this form, HSPs circulate via blood throughout the organism and act in an endocrine fashion: [77].

Conclusions
Altogether, based on bioinformatics analysis and experimental validation, we provided a reliable integrated analysis of the HSPs expression pattern, prognostic value, and potential regulatory mechanisms in breast cancer. HSPs represent a multitasking protein family exerting multi-level functions within BC cells. HSP network is dynamic within cancer cells whose biochemistry and function depend on cellular and environment milieu. In this intricate scenario, two different HSP-clusters here identified might serve as a potential panel for prognosis evaluation. The usefulness of obtained data should be considered prospectively because, in the era of personalized therapy, the molecular characterization of the two HSP-clusters should be given the possibility to stratify patients with different prognoses. Moreover, several clinical trials are currently investigating the use of HSP-inhibitors or HSP-based vaccines for cancer therapy and immunotherapy. We believe that, rather than a single HSP member, multiple HSPs should be targeted so that they have a better chance of reverting the clinical phenotype. In fact, the simultaneous deregulation of several HSPs in cancer cells could be orchestrating new interactive circuits that define tumor progression.
However, further validations are needed to judge their comprehensive prognostic impact in BC.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study when enrolled. Data Availability Statement: Not applicable.