Introduction

Thyroid carcinoma (THCA) is the most prevalent endocrine cancer, and its prevalence is rising globally [1]. Although the majority of patients have a good prognosis, a few recurrences do occur in THCA patients [2]. Therefore, there is an urgent need to identify promising predictors for prognosis and find a novel target for therapies. Biomarkers have recently been discovered to be important in the diagnosis, progression, and prognosis prediction of thyroid cancer, according to several studies [3, 4].

The fifth member of the heat shock protein family A (HSPA5) gene encodes an endoplasmic reticulum (ER) chaperone protein [binding immunoglobulin protein (BiP)] from the heat shock protein 70 (Hsp70) family [5]. BiP is a protein that is mostly present in the ER and may be found in several different organs. The ability of BiP to modify its function via post-translational modifications is one of its most essential functions. These modifications play critical roles in maintaining the steady state of the protein and allowing cells to respond to environmental changes [5]. Given that BiP is a critical regulator in the endoplasmic reticulum, variations in its expression, activation, or inhibition have been associated with cancer [5]. According to recent findings in head and neck cancer, HSPA5 is known to negatively regulate lysosomal activity by ubiquitination of mitochondrial ubiquitin ligase 1 (MUL1) [6]. Moreover, HSPA5 may play a role in granulin-epithelin precursor (GEP)-mediated cancer development in hepatocellular carcinoma, according to the study by Yip et al. [7]. Xia et al. concluded that HSPA5 is tightly linked to lung cancer development and bad prognosis [8]. Miura et al. suggested that increased expression of HSPA5/BiP is a new predictor of good prognosis in individuals who have advanced thymic cancer [9]. A study by Langer et al. indicated that in early tumour stages, elevated expression of HSPA5 could be responsible for limiting local tumour development [10]. Nevertheless, the functional role of HSPA5 in THCA is unknown.

Here, we examine the relationship between HSPA5 expression and clinicopathological features of THCA patients, evaluate the prognosis merits of HSPA5 in THCA, and develop a prognosis prediction model for THCA. We also look at the difference in HSPA5 expression between THCA tumour tissues and normal tissues. The possible roles of HSPA5 were investigated using protein-protein interaction (PPI) network building and functional enrichment analysis for the high and low expression groups of HSPA5. A further investigation into the probable role of HSPA5 in THCA was carried out using an immune infiltration analysis method. The purpose of this research is to investigate the expression of HSPA5 in THCA, as well as the potential diagnostic and therapeutic benefits of HSPA5 in THCA.

Material and methods

Data acquisition and bioinformatics analysis

We obtained data for 502 samples from the Cancer Genome Atlas (TCGA) THCA project. The 502 samples consisted of 58 para-carcinoma tissues and 502 tumour tissues. The data comprised clinicopathological features and RNA-sequencing data in level 3 HTSeq-FPKM format. The clinicopathological features include age, gender, race, TNM stage, pathologic stage, histological type, residual tumour, extrathyroidal extension, primary neoplasm focus type, neoplasm location, thyroid gland disorder history, overall survival (OS) event, and progression-free interval (PFI) event. Patients without corresponding clinical information were not included in the summary of the clinicopathological features (Tab. 1). RNAseq data in FPKM (fragments per kilobase per million) format were converted to TPM (transcripts per million reads) format for analysis. The gene expression profiles of GSE33630 and GSE58545 from the Gene Expression Omnibus (GEO) collection were utilized to verify the difference of HSPA5 expression between normal and malignant tissues in THCA patients.

Data analysis and visualization were performed using R (version 3.6.3). An alpha level of 0.05 was chosen as the threshold for statistical significance.

Differential expression of HSPA5

We used rms package, survival package, and basic package to compare the HSPA5 expression level between malignant tissues and normal tissues and assess the relationship between the expression of HSPA5 and clinicopathological features. Statistical methods included Shapiro-Wilk normality test, Levene’s test, t-test, Dunn’s test, Kruskal–Wallis test and Wilcoxon rank sum test, chi-square test, Fisher’s exact test, and logistic regression.

The relationship between HSPA5 expression and prognosis

We first performed ROC analysis (using the pROC package for analysis and the ggplot2 package) and Kaplan-Meier survival analysis (using the survminer package and the survival package) to evaluate the diagnostic and prognostic values of HSPA5 in THCA patients [11]. Subsequently, we performed univariate and multivariate regression analysis (using the survivor package) to explore the correlation between HSPA5 expression and PFI in THCA patients [11]. An alpha level of 0.05 was chosen as the threshold for statistical significance.

Construction of prognosis prediction model

Univariate and multivariate regression analyses (survivor package) were used to construct nomograms for prediction of PFI in THCA patients, and calibration plots were used for evaluation of the model accuracy [11]. The parameters in calibration visualization were as follows: the number of samples in each group calculated repeatedly (80); the number of repeat calculation (200); method (boot).

PPI network and enrichment analysis

To investigate HSPA5-related proteins, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was utilized to construct a network of PPI. Subsequently, we looked into the possible route and related proteins of HSPA5 using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG).

Relationship between HSPA5 expression and immune cell infiltration

Single-sample Gene Set Enrichment Analysis (ssGSEA; GSVA package) was performed to analyse the potential pathway related to HSPA5 in THCA [12]. This investigation comprised a total of 24 immune cells [13]. Statistical significance is indicated by the following symbols: ns — p ≥ 0.05; * — p < 0.05; ** — p < 0.01; *** — p < 0.001.

Results

Clinicopathological features

The data of 502 patients were acquired from TCGA. The clinicopathological features included age, gender, race, TNM stage, pathologic stage, histological type, residual tumour, extrathyroidal extension, primary neoplasm focus type, neoplasm location, thyroid gland disorder history, overall survival (OS) event, and PFI event. The details are presented in Table 1.

The relationship between HSPA5 expression and clinicopathological features of THCA

It was shown in Figure 1A that HSPA5 expression was reduced in THCA tissues (p < 0.001). HSPA5 expression was significantly reduced in THCA tissues (p < 0.001), as shown by a comparison of 58 THCA samples to their matched normal thyroid tissues in Figure 1B. Figure 1C shows the statistical difference in HSPA5 expression between different pathological types of thyroid cancer tissues and reveals that HSPA5 expression levels in normal thyroid tissues are significantly higher than those in different histological types of tumour tissues. Figure 1DE reveals that HSPA5 expression was validated in GSE60502 and GSE62232 (p < 0.001). As shown in Figure 1F, the AUC of HSPA5 was 0.783. Clinicopathological characteristics and HSPA5 expression levels are shown to be associated in Figure 2. HSPA5 was less expressed in T3 and T4 stage (p < 0.001), N1 stage (p < 0.001), pathologic III and IV stage (p < 0.001), unifocal type (p = 0.016), extrathyroidal extension (p < 0.001), R1 and R2 residual tumour (p = 0.03), and tumour recurrence in PFI event (p = 0.002). Also, the logistic regression results in Table II show that T stage (p < 0.001), N stage (p < 0.001), pathologic stage (p < 0.001), residual tumour (p = 0.037), extrathyroidal extension (p < 0.001), and primary neoplasm focus type (p < 0.01) have statistically significant differences.

Figure 1. Expression of HSPA5 in thyroid carcinoma (THCA) and normal or matched normal thyroid tissues. A. THCA tissues and normal thyroid tissues; B. THCA tissues and matched normal thyroid tissues; C. Different pathological types of THCA and normal thyroid tissues; D. THCA tissues and normal thyroid tissues in GSE58545; E. THCA tissues and normal thyroid tissues in GSE33630; F. Receiver operating characteristic (ROC) curve

Figure 2. HSPA5 expression in thyroid carcinoma (THCA) was significantly associated with clinical characteristics

The relationship between HSPA5 expression and recurrence of THCA patients

Figure 3A shows that among THCA patients reduced expression of HSPA5 was related to short PFI (p = 0.002). Furthermore, the results of subgroup analysis suggested that low expression of HSPA5 in white race (p = 0.011), female (p = 0.003), N1 stage (p = 0.001), pathologic stage I (p = 0.02), residual tumour R0 (p = 0.018), histological classical type (p = 0.027), and multifocal type (p = 0.013) were also correlated with short PFI of THCA patients (Fig. 3). Univariate Cox analysis shows that low HSPA5 levels (p = 0.042) are related to high T stage (p = 0.001), M stage (p < 0.001), pathologic stage (p < 0.001), and extrathyroidal extension (p = 0.023) in PFI events (Fig. 4A). Multivariate Cox analysis suggested that only M stage was independently associated with the PFI of THCA patients (p = 0.031; Fig. 4B).

Figure 3. Low HSPA5 expression in subgroups of clinical characteristics was associated with poor progression-free interval (PFI) in patients with thyroid carcinoma (THCA). A. All patients; B. White race; C. Female; D. N1 stage; E. Pathologic stage I; F. Residual tumour R0; G. Histological classical type; H. Multifocal type

Figure 4. Association of progression-free interval (PFI) and clinical characteristics in The Cancer Genome Atlas (TCGA) thyroid carcinoma (THCA) patients (Cox regression). A. Univariate analysis; B. Multivariate analysis

Nomograms for prediction of PFI in THCA patients

To increase the predictive accuracy of HSPA5 and to create an easy-to-use clinical model for predicting the chance of recurrence in THCA patients, we combined HSPA5 with significant clinicopathological features from univariate and multivariate regression analyses. The larger the total of the scores for each component in the nomogram, as shown in Figure 5A, the worse the PFI for the patient at 1, 3, and 5 years. A THCA patient with low HSPA5 expression and T3&T4 stage, for example, would receive a total of 90 points (27.5 points for low HSPA5 expression and 62.5 points for T3&T4 stage), with predicted 1-year, 3-year, and 5-year PFI rates of 95% to 90%, 90% to 80%, and 80%, respectively. In addition, we constructed calibration plots to evaluate the accuracy of the prediction model. The calibration plots (Fig. 5B) show good agreement between the nomogram’s predictions and the actual results. For the PFI nomogram, the concordance indices and bootstrapped 95% confidence intervals were 0.731 (0.690-0.771).

Figure 5. Nomogram to predict recurrence of the disease in thyroid carcinoma (THCA) patients. A. Prediction model of progression-free interval (PFI) rates in THCA patients; B. Accuracy verification of prognostic prediction model

PPI network and enrichment analysis

In this study, we used the STRING tool to build a network of HSPA5 and its related gene. The top 10 genes associated with HSPA5 included HSP90B1, EIF2AK3, ERN1, ATF6, MANF, CANX, DNAJC10, CALR, A2M, and DNAJB11 (Fig. 6). Their correlation scores were greater than or equal to 0.98. Furthermore, GO/KEGG enrichment analysis revealed that these genes were present in 229 pathways (217 GO and 12 KEGG). Among the 229 pathways, we included the top 3 low P-value datasets in GO cellular component (CC), GO molecular function (MF), GO biological process (BP), and KEGG, respectively. Our results show that the included enriched pathways were almost associated with endoplasmic reticulum function. The details are presented in Figure 7.

Figure 6. Protein-protein interaction comprehensive analysis of HSPA5. A. A network of HSPA5 and its potential co-expression genes was analysed by the STRING tool. The results were visualized by Cytoscape; B. Detailed information of HSPA5-related genes was listed

Figure 7. Enrichment analysis of HSPA5. A. The bubble plot showed the signalling pathway of HSPA5; B. A network between the potential co-expression genes of HSPA5 and its related pathway

The relationship between HSPA5 expression and immune infiltration

Except for NK CD56 bright cells, pDC, TFH, and Th17 cells, there is a significant connection between HSPA5 expression and infiltration levels of the rest of the immune cells (p < 0.05). The details are described in Figure 8.

Figure 8. The relationship between HSPA5 expression thyroid carcinoma (THCA) and immune cells. A. Differences in immune cell enrichment scores between HSPA5 high and low expression groups (Wilcoxon rank sum test); B. Correlation between HSPA5 expression in THCA and immune cells

Discussion

This study relies on bioinformatics analyses to explore HSPA5 expression in THCA and its potential therapeutic and prognostic merits for new directions of the therapeutic target. We found that variations in HSPA5 expression level were associated with the prognosis of THCA patients. Decreased HSPA5 expression indicated a poor prognosis in THCA. Furthermore, our findings revealed that the immune infiltration level of multiple immune cells were linked to the expression level of HSPA5. As a result, our findings contribute to a better understanding of the possible function of HSPA5 in immunotherapy and its role as a biomarker in THCA.

In our research, we built the PPI network to confirm the related gene and potential pathway of HSPA5. We discovered that the top 10 genes related to HSPA5 were all enriched in the pathways related to endoplasmic reticulum function (Fig. 7), such as protein processing in endoplasmic reticulum, response to endoplasmic reticulum stress, response to topologically incorrect protein, and response to unfolded protein. As an endoplasmic reticulum partner molecule, HSPA5 plays a crucial role in cell viability. [14–16]. Many previous investigations found that HSPA5 plays a major part in carcinoma progression and metastasis [17–22]. Looking back at the studies on HSPA5, we found that elevated expression of HSPA5 was identified in multiple distinct cancer types when compared to normal tissue and was correlated with the progression and poor prognosis in hepatocellular carcinoma, oesophageal squamous cell carcinoma, endometrial carcinoma, gastric adenocarcinoma, and breast cancer [23–27]. Interestingly, some studies have suggested an opposite role for HSPA5 in cancers. For instance, Miura et al. found that increased HSPA5 level is a predictor of a good prognosis in patients with advanced thymic carcinoma [9]. Wei et al. discovered that low expression of HSPA5 enhances cell migration in hepatocellular carcinoma cells and demonstrated that HSPA5 silencing enhances cancer metastasis through upregulation of vimentin in tumour cells [28]. A high HSPA5 level is related to a good prognosis in some cancers, most likely because cell-surface BIP, which is expressed by the HSPA5 gene, is a receptor for extrinsic tumour suppressor factors and plays an essential role in apoptosis [29, 30]. Moreover, BIP on the cell surface can operate as a pro-apoptotic regulator that collaborates with other tumour-suppressive proteins, whereas intracellular BiP possesses cell-protective functions [31].

Cancer immunotherapy has recently emerged as a popular research topic [32–35]. Understanding immuno-infiltrating cells can aid in exploring new immunotherapies. Nevertheless, few studies have been done on the link between HSPA5 and immunocytes in THCA. Therefore, in this work we investigated the relationship between HSPA5 and immunocyte infiltration. Our results show that HSPA5 expression was associated with infiltration of aDC, B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, iDC, macrophages, mast cells, neutrophils, NK CD56dim cells, T cells, T helper cells, Tcm, Tem, Th1 cells, Th2 cells, TReg, NK cells, and Tgd. These associations could point to possible ways by which HSPA5 inhibits the function of aDC, B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, iDC, macrophages, mast cells, neutrophils, NK CD56dim cells, T cells, T helper cells, Tcm, Tem, Th1 cells, Th2 cells, and TReg, and stimulates the activity of Tgd and NK cells.

Our research has some limitations. First, the number of normal samples from the TCGA database was relatively insufficient. Second, we have not yet thoroughly explored the direct mechanism of HSPA5 participation in THCA. These findings should be confirmed in future research. As a result, in the future, we are going to create cell and animal models to test our claims further.

Conclusions

To summarize, low HSPA5 expression level is related to poorer clinicopathological features, shorter PFI and higher immune infiltration level of multiple immune cells. These findings indicate that HSPA5 can be employed as a biomarker to predict the prognosis of patients with THCA.

Author’s contributions

Conception and design: WD. Collection and collation of data: WD, DD, and HH. Data analysis and interpretation: WD and HH. Manuscript writing and revisions: WD, DD, and HH. Final approval of manuscript: All authors. Accountable of all aspects of work: All authors.

Funding

This study was supported by the Science and Technology Program of Traditional Chinese Medicine in Zhejiang Province (2021ZB208) and Zhejiang Medical and Health Science and Technology Plan Project (2021KY927).

Acknowledgements

The authors would like to thank all the researchers who provided the data for this study.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the TCGA and GEO databases.

Ethical approval and consent to participate

All data in this research were from the TCGA and GEO public databases, so ethical approval and informed consent were exempted.