Unfolded Protein Response–Related Signature Associates With the Immune Microenvironment and Prognostic Prediction in Osteosarcoma

Background: Osteosarcoma is a highly malignant bone tumor commonly occurring in adolescents with a poor 5-year survival rate. The unfolded protein response (UPR) can alleviate the accumulation of misfolded proteins to maintain homeostasis under endoplasmic reticulum stress. The UPR is linked to the occurrence, progression, and drug resistance of tumors. However, the function of UPR-related genes (UPRRGs) in disease progression and prognosis of osteosarcoma remains unclear.

Methods: The mRNA expression profiling and corresponding clinical features of osteosarcoma were acquired from TARGET and GEO databases. Consensus clustering was conducted to confirm different UPRRG subtypes. Subsequently, we evaluated the prognosis and immune status of the different subtypes. Functional analysis of GO, GSEA, and GSVA was used to reveal the molecular mechanism between the subtypes. Finally, four genes (STC2, PREB, TSPYL2, and ATP6V0D1) were screened to construct and validate a risk signature to predict the prognosis of patients with osteosarcoma.

Result: We identified two subtypes according to the UPRRG expression patterns. The subgroup with higher immune scores, lower tumor purity, and active immune status was linked to a better prognosis. Meanwhile, functional enrichment revealed that immune-related signaling pathways varied markedly in the two subtypes, suggesting that the UPR might influence the prognosis of osteosarcoma via influencing the immune microenvironment. Moreover, prognostic signature and nomogram models were developed based on UPRRGs, and the results showed that our model has an excellent performance in predicting the prognosis of osteosarcoma. qPCR analysis was also conducted to verify the expression levels of the four genes.

Conclusion: We revealed the crucial contribution of UPRRGs in the immune microenvironment and prognostic prediction of osteosarcoma patients and provided new insights for targeted therapy and prognostic assessment of the disease.

Introduction

Osteosarcoma is one of the most frequent primary malignant bone tumors in adolescents, often occurring in fast-growing long bones (Tang et al., 2008; Arndt et al., 2012). With significant advances in surgical options, neoadjuvant chemotherapy, and diagnostic imaging, the overall survival of patients with osteosarcoma has been significantly prolonged (Li et al., 2015; Wakamatsu et al., 2019; Hao et al., 2021). Nevertheless, the outcome of metastatic and recurrent patients remains unsatisfactory due to its intensely aggressive feature, with a 5-year mortality ratio of below 20% (Ferrari et al., 2005). In addition, high tumor heterogeneity, leading to chemoresistance in some patients, remains challenging clinically (Xiao et al., 2015). Therefore, there is a need for the discovery of novel targets for targeted treatment and better prognostic markers for osteosarcoma.

The unfolded protein response (UPR) is a highly conserved homeostatic response in all eukaryotic cells and can help cells mitigate the accumulation of misfolded proteins in the endoplasmic reticulum (ER) (Hetz et al., 2020). Under ER stress conditions, such as hypoxia, nutritional deprivation, acidosis, and inflammatory stimuli, the UPR can activate three sensor proteins (IRE1α, PERK, and ATF6) to initiate specific signal transduction cascades which regulate the rate of protein production for maintaining ER homeostasis (Vanacker et al., 2017). Lately, extensive studies have indicated that the UPR acts a dual function in tumor occurrence and development. In earlier stages of tumorigenesis, the UPR is capable of exerting antitumor effects to hinder tumor transformation, while in the established tumors, tumor cells can induce chronic UPR to relieve ER stress-induced apoptosis and develop drug resistance to maintain tumor survival (Ma and Hendershot, 2004; Yang et al., 2019). Aberrant activation of the UPR was found in a wide range of tumors, including bladder cancer, cutaneous melanoma, and liver cancer (Houessinon et al., 2016; Wan et al., 2020; Zhu et al., 2021). Apart from a direct impact on tumor biology, the UPR also has the ability to remodel the tumor immune microenvironment (TIME) to regulate the crosstalk between immune cells, which serves an essential function in immune surveillance and immune escape (Zanetti et al., 2022). Mahadevan et al. demonstrated that the UPR can modulate the phenotype of dendritic cells and CD8⁺ T cells to facilitate tumor growth. Notably, previous studies have confirmed the UPR to be relevant to the growth, prognosis, and drug resistance of osteosarcoma (Mahadevan et al., 2012; Yan et al., 2015; Shi et al., 2022). However, the role of the UPR-related gene (UPRRG) sets on prognostic prediction and immune infiltration in patients with osteosarcoma remains unclear.

Recently, following the advancement of the high-throughput sequencing technology for tumor genomics, probing new molecular patterns through bioinformatics approaches offers new insights for tumor treatment and prognosis evaluation (Qian et al., 2021). In the present study, we performed a comprehensive analysis of the prognosis and immune microenvironment of different molecular subtypes based on the expression of UPRRGs and revealed potential functional and signaling pathways. Moreover, we developed and validated a novel prognostic signature for predicting the prognosis of osteosarcoma patients via integrating risk scores and clinical features.

Methods and Materials

Data Collection

The mRNA expression profiles and corresponding clinical features of 85 osteosarcoma patients were acquired from the TARGET database as a training cohort (https://ocg.cancer.gov/programs/target). UPRRGs were extracted from hallmark gene sets from the Molecular Signature Database (Zhang et al., 2021). In addition, mRNA expression profiles and relevant clinical information of 53 osteosarcoma patients from GSE21257 were obtained from the GEO database to serve as an external validation cohort (https://www.ncbi.nlm.nih.gov/geo/). The clinical information for osteosarcoma patients is presented in Supplementary Table S1.

Classification of Molecular Subtypes

To assess the biological functions of UPRRGs in OS, we first identified 15 prognosis-related UPRRGs for osteosarcoma based on a univariate Cox regression analysis. Subsequently, consensus clustering analysis was conducted based on the expression matrix of these 15 genes with the R package “ConsensusClusterPlus”. For the analysis, 80% of the total samples of the target dataset are included in each iteration and reiterated 1,000 times to ensure cluster stability. The average consistency value and clinical significance within the clustering group were used together with the optimal number of clusters. Principal component analysis (PCA) was performed to confirm the gene expression patterns in the different clusters. KM curve and log-rank tests were employed to evaluate the prognosis of different clusters. The heat map was used to show the gene expression and clinical characteristics of the different clusters.

TIME Evaluation

The tumor microenvironment scores for individual samples in the different subtypes were assessed by the ESTIMATE algorithm (Yoshihara et al., 2013). The TIMER algorithm was conducted to comprehensively evaluate the abundance of immune infiltrating cells in each sample (Li et al., 2020). Moreover, the infiltration abundance of 28 immune cell types for an individual sample was tested by the single sample gene set enrichment analysis (ssGSEA) algorithm. The expression of immune checkpoint (ICP) genes in different subtypes was also evaluated.

Identification of DEG and Enrichment Analysis

Differentially expressed genes (DEGs) in different subtypes were analyzed via the R package “limma”, and log2 (Foldchange) > 1.5 and FDR<0.05 were considered as the threshold. Gene Ontology (GO) terms of DEGs were enriched by using the R package “clusterProfiler” and visualized via Metascape (https://metascape.org) (Zhou et al., 2019). Gene set variation analysis (GSVA) was conducted to verify functional pathway variations between the different subtypes according to GO: the biological process (BP). Also, Gene Set Enrichment Analysis (GSEA) was employed to analyze hallmark gene sets for different subtypes (Subramanian et al., 2005). |NES| > 1, NOM p-val < 0.05, and FDR<0.25 were taken as the threshold.

Establishment and Validation of the Prognostic Signature

The previously obtained 15 prognosis-related UPRRGs were screened via least absolute shrinkage and selection operator (LASSO) Cox regression based on the R package “glmnet”, and the minimum lamba is taken as the optimal value. Then, a multivariate Cox analysis was conducted to further optimize and establish the prognostic signature. The risk score in the training and validation cohorts was calculated with the following formula: Risk score = ∑ⁱ_n (Coef_i * X_i). In the training cohort, all patients were classified into low- and high-risk groups by the median risk score, and the overall survival (OS) between the two groups was investigated using the KM curve and log-rank tests. In addition, receiver operating characteristic (ROC) curves were applied to measure the effectiveness of prognostic models. Furthermore, the abovementioned formula was also used in patients from GSE23257 to generate a risk model to validate the prognostic benefit. Finally, we integrated different clinical characteristics including gender, age, metastatic status, and disease site to assess if the risk score was an independent prognostic element in osteosarcoma patients based on univariate and multivariate Cox regression.

Construction and Calibration of the Nomogram

A nomogram model was established to forecast the 3-year and 5-year survival probability of osteosarcoma patients via integrating risk scores and clinical profiles. The C-index, ROC curves, and calibration plots were employed to survey the predictive performance of the constructed nomogram in both cohorts.

Patient Sample Collection and Cell Culture

A total of six osteosarcoma patient tissues and corresponding adjacent normal tissues were obtained from the operation and stored in liquid nitrogen banks. The study was approved by the institutional review board of Xijing Hospital, the Fourth Military Medical University, and all patients provided signed informed consent. Moreover, human osteosarcoma cell lines HOS and MG-63 were procured from Procell Life Science & Technology Co. Ltd and cultured in DMEM medium (Gibco, Shanghai, China) containing 10% fetal bovine serum (FBS; Gibco, Shanghai, China) and 1% (v/v) penicillin/streptomycin (Gibco, Shanghai, China). The osteoblast cell line hFOB 1.19 was provided by Dr. Jianping Bai of Xijing Hospital and cultured in DMEM/F12 medium with 10% FBS, 0.3 mg/ml G418 (Procell, Wuhan, China), and 1% (v/v) penicillin/streptomycin. At 37 C, 5% CO₂ environment, all cell lines were incubated.

Quantitative Real-Time PCR (qRT-PCR)

The TRIzol method was utilized to extract and purify RNA from tissues and cells. Then, a cDNA synthesis kit (Takara, China) was applied to reverse transcribe the RNA. The TB Green Premux Ex TaqⅡ (Tli RNaseH Plus) was used for qRT-PCR with the Bio-Rad CFX96 Real-Time PCR system (Bio-Rad, USA). The internal control was GADPH. The primer sequences of the candidate genes are shown in Table 1, and the analysis was conducted three times for all genes.

TABLE 1

TABLE 1. Primer sequences of the candidate genes.

Statistics

R software (version 4.0.5), SPSS 21.0 software, and GraphPad Prism 8 were carried out for all statistical analyses. The t-test was applied for two groups. One-way ANOVA was applied to three groups. p < 0.05 was taken as statistically significant. ∗p < 0.05; ∗∗p < 0.01; and ∗∗∗p < 0.001.

Results

Identification of UPRRG Molecular Subtypes via Consensus Analysis

A total of 113 UPRRGs were retrieved from hallmark gene sets, out of which we identified 15 prognosis-related UPRRGs for osteosarcoma, according to univariate COX analysis (Supplementary Table S2, 3). According to the expression profile of these genes, consensus clustering analysis was used to ascertain subgroups of osteosarcoma patients in the training database. k = 2 is considered the best category number of clusters, depending on the average consistency value and clinical significance within the clustering group (Figure 1A and Supplementary Figure S1). PCA analysis revealed a relatively apparent distinction between the two subtypes (Figure 1B). The heat map illustrated the gene expression profile and clinical characteristics of the two subtypes (Figure 1C). Moreover, we noticed that patients in cluster 2 experienced a dismal prognosis to that in cluster 1 (p<0.001, Figure 1D). Previous studies have proven that the UPR can coordinate the crosstalk between immune cells and tumor cells in the TIME to exert immunosurveillance and immunosuppressive functions to influence tumor prognosis (Vanacker et al., 2017). Thus, we then evaluated the differences in the TIME across different subtypes. The ESTIMATE algorithm indicated that cluster 1 had higher immune score (p = 8.4e-5), stromal score (p = 8.1e-7), ESTIMATE score (p = 9.6e-7), and lower tumor purity (p = 1.1e-6), as compared to cluster 2 (Figure 2A). In addition, the TIMER algorithm discovered that the abundance of most immune infiltrating cells in cluster1 was significantly increased than in cluster 2, including dendritic cells, neutrophil cells, CD4⁺ T cells, neutrophil cells, and CD8⁺ T cells, while B cells showed the opposite result (Figure 2B). As shown in Figure 2C, ssGSEA analysis found that the abundance of 20 immune cell types was significantly increased in cluster 1 compared to cluster 2. In addition, we also observed that CD274, LAG3, HAVCR2, and PDCD1 were expressed at an elevated level in cluster 1 than in cluster 2 (Figure 2D). Our findings suggested that the prognosis of different subtypes may be affected by the TIME.

FIGURE 1

FIGURE 1. Identification of molecular subtypes of UPRRGs by consensus clustering. (A). Clustering heat map at k = 2. (B) PCA plot between the two subtypes. (C) Heat map of the UPR-related gene expression and clinical features in the two subtypes. (D) Survival curves for the two subgroups.

FIGURE 2

FIGURE 2. Tumor microenvironment in the two subtypes. (A) Stromal score, immune score, ESTIMATE score, and tumor purity based on ESTIMATE algorithm. (B) Six immune cell abundance assessments by the TIMER algorithm. (C) Twenty-nine immune cell abundance assessments by the ssGSEA algorithm. (D) Differences in immune checkpoints between the two subtypes. *p < 0.05; **p < 0.01; ***p < 0.001.

Functions and Pathway Annotations of DEGs for UPRRG Subtypes

To reveal the potential mechanisms regulating the TIME between different subtypes, we performed DEGs analysis on the two clusters. The results indicated that 121 genes were upregulated and 157 genes were downregulated in cluster 2 to cluster 1 (Figure 3A). GO analysis implied that these DEGs were primarily associated with inflammatory response, leukocyte activation, and extracellular matrix (Figure 3B). These findings implied that UPR could influence the TIME and prognosis of osteosarcoma via modulation of immune-associated pathways. We then used GSVA and GSEA analysis to explore the functional differences in the two clusters. GSEA analysis demonstrated that coagulation, inflammatory response, and IL6/JAK/STAT3 signaling were markedly upregulated in cluster 1 (Figure 3C). GSVA analysis showed that positive regulation of calcium ion import, regulation of the apoptotic process involved in the development, and some immune-related pathways were significantly upregulated in cluster 1 to cluster 2 (Figure 3D). Therefore, we speculated that the UPR plays an essential role in regulating the immune function, thus contributing to the prognosis of osteosarcoma.

FIGURE 3

FIGURE 3. Differentially expressed genes and functional enrichment analyses. (A) Volcano plot showing the DEGs between the two subgroups. (B) Bubble plot exhibited the functional enrichment of DEGs through GO analysis. (C) GSEA shows the hallmark gene sets in the two subgroups. (D) Heat map depicted the results of the GSVA analysis.

Construction and Validation of the UPRRG Risk Signature

To further elucidate the prognostic predictive role of UPRRGs in osteosarcoma patients, LASSO regression analysis was applied to screen for potential genes, and eight genes were identified by the minimal lambo value (Figures 4A,B). Finally, based on the results of previous screening, four UPRRGs (STC2, PREB, TSPYL2, and ATP6V0D1) were identified to establish a risk model for osteosarcoma via the multivariate Cox regression analysis (Figure 4C and Supplementary Table S4). The following formula was applied to generate the risk score for each sample: Risk score = -1.523 × ATP6V0D1 + 0.903× PREB +0.586 × STC2 -0.760 × TSPYL2. All patients in the training cohort were classified into high- and low-risk groups by the median risk scores, and patients in the high-risk group showed shorter survival intervals than those in the low-risk group (Figure 4D). KM curves demonstrated a poor prognosis in the high-risk group to the low-risk group (Figure 4E). Moreover, time-dependent ROC analysis found the AUC values of 1, 3, and 5 years were 0.84, 0.87, and 0.83, respectively, which suggested our risk signature showed excellent predictive performance (Figure 4F). Notably, we also observed that the low-risk group saw higher immune scores (p = 3.5e-4), stromal scores (p = 1.1e-4), ESTIMATE scores (p = 3.3e-5), and lower tumor purity (p = 2.9e-5) relative to the high-risk group, which suggested that the TIME might be strongly linked to prognosis in different risk groups (Figure 4G). Then, the relevance between the risk signature and clinical features was also evaluated, and the results revealed metastatic patients had significantly higher risk scores than non-metastatic patients (p = 0.03), while no differences were found in any other clinical characteristics (Figure 5A). When patients were reclassified for metastatic status, there was a significantly improved prognosis for patients in the low-risk group over those in the high-risk group (Figures 5B,C). These findings supported that the UPRRG risk signature might be strongly correlated with the metastasis status in patients with osteosarcoma. Furthermore, multivariate Cox analysis demonstrated that the risk score and metastatic status were independent prognostic factors for osteosarcoma patients, which meant that the UPRRGs risk model was applicable to predict survival in osteosarcoma (Figures 5D,E).

FIGURE 4

FIGURE 4. Construction of the risk signature based on UPRRGs in the training cohort. (A) Eight optimal UPRRGs filtered by LASSO analysis. (B) Distribution of risk scores and patient status in the two risk groups; (C). Heat map showing the expressions of four candidate genes. (D) Survival curves for the two risk groups. (E) Time-dependent ROC curve of the risk model. (F) Tumor microenvironment analysis in the two risk groups through the ESTIMATE algorithm. *p < 0.05; **p < 0.01; ***p < 0.001.

FIGURE 5

FIGURE 5. Correlation of the risk signature with clinical features in the training cohort. (A) Differences in risk scores among osteosarcoma patients by age, gender, and metastatic status. (B–C) Survival curves for patients with osteosarcoma regrouped by metastatic status. (D–E) Univariate and multivariate cox regression analyses for integrating risk characteristics and clinical features.

Moreover, we further tested the prognostic performance of the UPRRG risk signature in a validation cohort. As shown in Figure 6A, the patients in the validation cohort were clearly separated into different risk groups via the abovementioned formula. The heat map demonstrated these four genes’ expression profiles in subgroups (Figure 6B). The KM curve likewise showed that the low-risk group had a better prognosis (p = 0.04) (Figure 6C). ROC curves suggested the risk signature had better prediction accuracy at 1 and 3 years (Figure 6D). Similarly, the ESTIMATE algorithm obtained results consistent with the training cohort, which further confirmed the role of the TIME in the UPRRG risk signature (Figure 6E).

FIGURE 6

FIGURE 6. Validation of the constructed risk signature in the verification cohort. (A) Distribution of risk scores and patient status in different risk groups. (B) Heat map displayed the expressions of four candidate genes in the verification cohort. (C) Survival curves of the two risk groups. (D) Time-dependent ROC curve in the verification cohort. (E) Tumor microenvironment analysis by the ESTIMATE algorithm. *p < 0.05; **p < 0.01; ***p < 0.001.

Construction and Validation of the Nomogram Prediction Model

A nomogram model was constructed to improve the accuracy of predicting the prognosis of osteosarcoma patients at 3 and 5 years by integrating risk scores and clinical characteristics (Figure 7A). Then, we validated the predictive efficacy of the nomogram model in the two cohorts. The C-index values for the training and validation cohorts were 0.88 and 0.87, respectively. The ROC curves revealed AUC values of 0.93 and 0.90 for the training cohort at 3 and 5 years, respectively, and for the validation cohort, ROC curves also exhibited excellent prediction accuracy (Figures 7B,C). Moreover, the calibration curves for the training and validation cohort showed that the nomogram model has a strong predictive capacity for the prognosis of osteosarcoma patients at 3 and 5 years (Figures 7D,E). Collectively, the aforementioned results pointed out that the UPRRG nomogram model has high predictive accuracy and can be applied to predict the prognosis of osteosarcoma patients.

FIGURE 7

FIGURE 7. Construction and evaluation of the nomogram. (A). Nomogram for predicting the prognosis of patients with osteosarcoma. (B) Calibration for 3-and 5-year OS in the training cohort. (C) Calibration for 3-and 5-year OS in the verification cohort. (D) ROC analysis for 3-and 5-year OS in the training cohort. (E) ROC analysis for 3-and 5-year OS in the verification cohort.

Verification of Candidate Genes by qRT-PCR Analysis

To certify the expression levels of these four candidate genes, we performed a qPCR analysis in patients' tissue and cell lines. The results revealed that the expression levels STC2 and PREB were elevated clearly in osteosarcoma tissues than in normal tissues, whereas the ATP6V0D1 expression level was significantly downregulated in tumor tissues compared with normal tissues (Figure 8A). The results of the cell further identified that the expression levels of PREB and STC2 were higher in osteosarcoma cells than in normal osteoblasts, while TSPYL2 and ATP6V0D1 exhibited the opposite results (Figure 8B).

FIGURE 8

FIGURE 8. Candidate gene validation. (A) Expression levels of candidate genes in tumor and normal tissues from osteosarcoma patients. (B) Expression levels of candidate genes in different cell lines. *p < 0.05; **p < 0.01; ***p < 0.001.

Discussion

The unfolded protein response is an adaptive signaling pathway that regulates the endoplasmic reticulum homeostasis (Hetz et al., 2020). The UPR has been extensively studied in tumorigenesis, and its abnormal activation is involved in various stages of tumorigenesis and progression (Ma and Hendershot, 2004). Chaiyawat et al. revealed that plenty of UPR-related proteins were highly expressed in osteosarcoma patients and closely associated with chemoresistance, suggesting that targeting UPR pathways might be promising for the treatment of osteosarcoma (Chaiyawat et al., 2019). Ji et al. found that PERK could induce autophagy in osteosarcoma by inhibiting the mTORC1 pathway to counteract ER stress-induced cell apoptosis (Ji et al., 2015). Yan et al. demonstrated that the UPR could inhibit cisplatin-mediated apoptosis by triggering the NF-κB pathway, contributing to drug resistance in osteosarcoma (Yan et al., 2015). Shi et al. identified aberrant activation of UPR-related pathways in osteosarcoma and built a prognostic model for differentially expressed genes (Shi et al., 2022). Although previous studies have shown the UPR was implicated in the development, progression, and treatment of osteosarcoma, the function of UPRRGs on immune infiltration and prognostic prediction in patients with osteosarcoma is not clear.

In this research, we determined two molecule subtypes of osteosarcoma based on prognosis-related UPRRGs and presented the relevance of different subtypes to clinical profiles. Our study found that different subtypes had different survival intervals and TIME. The TIME consists of various components, including immune and stromal cells, which together contribute to tumorigenesis, progression, and prognosis (Zhao et al., 2018). The UPR is a critical factor in regulating the balance of immune dynamics in the tumor microenvironment which can directly affect innate and adaptive immune responses to participate in tumor progression (Martins et al., 2016). Our results indicated that cluster 1 has a better prognosis with higher immune and stromal scores compared to cluster 2. Previous studies have demonstrated that a higher immune score and stromal score were connected to a better prognosis in osteosarcoma, which is in agreement with our findings (Qian et al., 2021). Further analysis indicated that the abundance of multiple immune infiltrating cells was distinctly greater in cluster 1 than in cluster 2, according to the TIMER, indicating that the active immune status of the UPRRG subtypes might be closely better associated with osteosarcoma prognosis. ssGSEA analysis found that the abundance of a variety of tumor-infiltrating lymphocytes was expressed at obviously higher levels in cluster 1 than in cluster 2. Previous studies have demonstrated that the UPR acts as a critical mediator of tumor immunity as an appropriate UPR can induce immune cells to eliminate tumor cells, whereas the sustained UPR could induce immune cell apoptosis to enhance tumor proliferation and invasion (Vanacker et al., 2017). Hence, we believed that the elevated abundance of multiple TILs in cluster 1 might be due to appropriate UPR, which induced activation of TILs to deliver anti-tumor immunity. Nevertheless, the sustained UPR led to massive TIL exhaustion, causing a poor prognosis in patients with osteosarcoma in Cluster 2. Moreover, we also noted that as a consequence of a large number of immune infiltration cells in cluster 1, the expressions of several ICPs were correspondingly elevated in cluster 1 compared to Cluster 2, implying that patients in cluster 1 might be more sensitive to ICP inhibitors. These findings suggested that the identification of UPRRG subtypes may provide a new clinical strategy for prognostic evaluation and targeted therapy of osteosarcoma.

To reveal the molecular mechanisms underlying regulating TIME between different UPRRG subtypes, we carried functional enrichment analysis of the DEGs between the two subtypes. GO analysis suggested that the DEGs were mainly participating in immune-related pathways including inflammatory response and leukocyte activation. In addition, GSEA analysis also confirmed significant enrichment of some immune-related pathways in cluster 1, including coagulation, inflammatory response, and IL6/JAK/STAT3 signaling. The UPR is proven to not only affect the growth and survival of tumor cells but also play an essential role in remodeling the TIME (Zanetti et al., 2022). Batista et al. identified that the UPR participates in the macrophage polarization in the TIME by activating the IRE1α/XBP1 axis, resulting in upregulation of IL-6, IL-23, arginase 1, CD86, and PD-L1 that lead to local immune dysregulation (Batista et al., 2020). Zanetti et al. showed that tumor cells regulate immune cell phenotypes through UPR activation of dendritic cells and T cells to promote tumor proliferation (Mahadevan et al., 2012). In addition, Medel et al. also demonstrated that dendritic cells could enhance CD8⁺ T-cell–specific responses through activation of the IRE1α/XBP1 axis to exert an anti-tumor effect (Medel et al., 2018). Taken together, these findings revealed that the UPR can modulate immune cell infiltration through various immune signaling pathways to improve the prognosis of patients with osteosarcoma. Moreover, GSVA analysis revealed that aside from immune-related pathways, apoptosis and calcium homeostasis-related pathways were also upregulated in cluster 1. These results suggested that the UPR may exert antitumor effects through multiple pathways between UPRRG subtypes. Overall, our data showed that patients with osteosarcoma can activate the UPR in the presence of ER stress to modulate tumorigenesis and progression through multiple pathways and that targeting UPR might be a promising treatment strategy for osteosarcoma.

Next, to assess the role of UPRRGs in predicting the prognosis of osteosarcoma, we constructed a prognostic signature to predict the survival of osteosarcoma patients via four genes (STC2, PREB, TSPYL2, and ATP6V0D1). Our results found that high expressions of STC2 and PREB were linked to high risk (risk factors), whereas high expressions of TSPYL2 and ATP6V0D1 were linked to low risk (protective factors). STC2 encodes a glycoprotein that performs an essential function in the development and invasion of multiple tumors (Li S. et al., 2021). Previous studies have demonstrated endoplasmic reticulum stress could activate PERK-ATF4 to induce the upregulated expression of STC2 to inhibit cell apoptosis (Ito et al., 2004). Chen et al. found that STC2 could promote tumor proliferation by activating the AKT-ERK pathway, and increased STC2 was strongly correlated with poor prognosis in colorectal tumors (Chen et al., 2016). PREB can encode transcription factors that bind and activate the basal prolactin promoter activity. Murao et al. showed that PREB could act as a transcription factor for TNF-αand IL-1β to regulate the expression of monocyte chemoattractant protein-1, suggesting that PREB plays an active role in immune responses (Murao et al., 2009). In addition, PREB is a member of a eukaryotic family of WD-repeat proteins involved in many biological activities, including vesicle trafficking, RNA processing, and signal transduction (Taylor Clelland et al., 2000). TSPYL2, a member of the TSPY-L nucleosome assembly protein-1 superfamily, can exert anti-tumor effects by inhibiting the cell cycle and regulating DNA damage. Previous studies have found that TSPYL2 maintains the G1 checkpoint function by inducing p21 transcription to modulate DNA damage and inhibit cellular growth (Tao et al., 2011). Liu et al. indicated that TSPYL2 could inhibit SIRT1-mediated FOXO3 deacetylation to reduce gefitinib resistance and inhibit DNA damage, implying that TSPYL2 is a promising therapeutic target (Liu et al., 2022). ATP6V0D1 is an encoded protein involved in vacuolar ATPase formation, which has a crucial role in the modulation of the acidic microenvironment (Lu et al., 2021). Numerous research studies have found that dysregulation of V-ATPase is related to tumor growth and invasion (Whitton et al., 2018). Targeting V-ATPase could upregulate ER stress-related markers and inhibit mTOR signaling to exert anticancer effects (Kitazawa et al., 2017).

Survival analysis and ROC curve presented that the UPRRG risk signature demonstrated satisfactory predictive accuracy in both cohorts and that the low immune status in the high-risk group was strongly connected with poor prognosis. Moreover, our results identified that risk score was an independent prognostic factor for osteosarcoma patients. Subsequently, we performed a nomogram model to better predict the prognosis of osteosarcoma patients via risk score and clinical characteristics. Previous studies have reported several nomogram models to predict the prognosis of osteosarcoma (Qian et al., 2021). Li et al. built a nomogram model by autophagy-related genes with AUC values of 0.735 and 0.726 at 3 and 5 years, respectively (Li J. et al., 2021). Wen et al. developed a 3-gene nomogram model with 3-year and 5-year AUC values of 0.853 and 0.818, respectively (Wen et al., 2020). Wu et al. established a hypoxic nomogram model with an AUC value of only 0.73 (Wu et al., 2021). In our study, the constructed nomogram model had 3-year and 5-year AUC values of 0.93 and 0.90, respectively, which were consistently superior to other models. We presented that the UPRRG nomogram model could better predict the prognosis of patients with osteosarcoma than existing models. Collectively, these findings offered new options for personalized treatment and prognostic prediction of osteosarcoma.

Despite the many strengths of the current study, there are notable limitations. First, the UPRRG risk signature was constructed based on the TARGET and GEO databases, which may be biased due to the limited number of patients. Second, this study lacked some clinical information relevant to the prognosis of osteosarcoma, such as tumor pathological grade, which constrained clinical variables that can be incorporated into the nomogram model. A larger, multicenter, prospective clinical cohort will be needed to further evaluate the clinical merit of our findings in the future.

Conclusion

Our study identified two molecular subtypes through consensus clustering based on prognosis-related UPRRGs, and the two subtypes exhibited different survival times and immune statuses. Functional analysis revealed that the UPRRG subtypes might influence the progression and prognosis of osteosarcoma patients through immune-related pathways. Moreover, a novel prognostic model based on UPRRGs was constructed and validated to better predict the prognosis of patients with osteosarcoma. We elucidated the important function of UPRRGs in the development and prognosis of osteosarcoma, shedding light on new insights for targeted therapy and clinical decision-making in patients with osteosarcoma.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Xi-Jing Hospital, The Fourth Military Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

HF and ZZ have designed the study. ZZ and XL have written the manuscript, and DC, JD, ZM, YS, and LW have analyzed the data. All authors have read and approved it for publication.

Funding

This work was supported by the National Natural Science Foundation of China (No. 31971272).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.911346/full#supplementary-material

References

Arndt, C. A. S., Rose, P. S., Folpe, A. L., and Laack, N. N. (2012). Common Musculoskeletal Tumors of Childhood and Adolescence. Mayo Clin. Proc. 87 (5), 475–487. doi:10.1016/j.mayocp.2012.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Batista, A., Rodvold, J. J., Xian, S., Searles, S. C., Lew, A., Iwawaki, T., et al. (2020). IRE1α Regulates Macrophage Polarization, PD-L1 Expression, and Tumor Survival. PLoS Biol. 18 (6), e3000687. doi:10.1371/journal.pbio.3000687

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaiyawat, P., Sungngam, P., Teeyakasem, P., Sirikaew, N., Klangjorhor, J., Settakorn, J., et al. (2019). Protein Profiling of Osteosarcoma Tissue and Soft Callus Unveils Activation of the Unfolded Protein Response Pathway. Int. J. Oncol. 54 (5), 1704–1718. doi:10.3892/ijo.2019.4737

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Zeng, X., He, Y., Wang, X., Liang, Z., Liu, J., et al. (2016). STC2 Promotes the Epithelial-Mesenchymal Transition of Colorectal Cancer Cells through AKT-ERK Signaling Pathways. Oncotarget 7 (44), 71400–71416. doi:10.18632/oncotarget.12147

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrari, S., Smeland, S., Mercuri, M., Bertoni, F., Longhi, A., Ruggieri, P., et al. (2005). Neoadjuvant Chemotherapy with High-Dose Ifosfamide, High-Dose Methotrexate, Cisplatin, and Doxorubicin for Patients with Localized Osteosarcoma of the Extremity: a Joint Study by the Italian and Scandinavian Sarcoma Groups. Jco 23 (34), 8845–8852. doi:10.1200/jco.2004.00.5785

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, Y., An, R., Xue, Y., Li, F., Wang, H., Zheng, J., et al. (2021). Prognostic Value of Tumoral and Peritumoral Magnetic Resonance Parameters in Osteosarcoma Patients for Monitoring Chemotherapy Response. Eur. Radiol. 31 (5), 3518–3529. doi:10.1007/s00330-020-07338-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Hetz, C., Zhang, K., and Kaufman, R. J. (2020). Mechanisms, Regulation and Functions of the Unfolded Protein Response. Nat. Rev. Mol. Cell. Biol. 21 (8), 421–438. doi:10.1038/s41580-020-0250-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Houessinon, A., Gicquel, A., Bochereau, F., Louandre, C., Nyga, R., Godin, C., et al. (2016). Alpha-fetoprotein Is a Biomarker of Unfolded Protein Response and Altered Proteostasis in Hepatocellular Carcinoma Cells Exposed to Sorafenib. Cancer Lett. 370 (2), 242–249. doi:10.1016/j.canlet.2015.10.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, D., Walker, J. R., Thompson, C. S., Moroz, I., Lin, W., Veselits, M. L., et al. (2004). Characterization of Stanniocalcin 2, a Novel Target of the Mammalian Unfolded Protein Response with Cytoprotective Properties. Mol. Cell. Biol. 24 (21), 9456–9469. doi:10.1128/mcb.24.21.9456-9469.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, G.-r., Yu, N.-c., Xue, X., and Li, Z.-g. (2015). PERK-mediated Autophagy in Osteosarcoma Cells Resists ER Stress-Induced Cell Apoptosis. Int. J. Biol. Sci. 11 (7), 803–812. doi:10.7150/ijbs.11100

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitazawa, S., Nishizawa, S., Nakagawa, H., Funata, M., Nishimura, K., Soga, T., et al. (2017). Cancer with Low Cathepsin D Levels Is Susceptible to Vacuolar (H+ )-ATPase Inhibition. Cancer Sci. 108 (6), 1185–1193. doi:10.1111/cas.13240

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Guo, Z., Wang, Z., Fan, H., and Fu, J. (2015). Does Microwave Ablation of the Tumor Edge Allow for Joint-Sparing Surgery in Patients with Osteosarcoma of the Proximal Tibia? Clin. Orthop. Relat. Res. 473 (10), 3204–3211. doi:10.1007/s11999-015-4447-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Tang, X., Du, Y., Dong, J., Zhao, Z., Hu, H., et al. (2021a). Establishment of an Autophagy-Related Clinical Prognosis Model for Predicting the Overall Survival of Osteosarcoma. BioMed Res. Int. 2021, 1–17. doi:10.1155/2021/5428425

CrossRef Full Text | Google Scholar

Li, S., Huang, Q., Li, D., Lv, L., Li, Y., and Wu, Z. (2021b). The Significance of Stanniocalcin 2 in Malignancies and Mechanisms. Bioengineered 12 (1), 7276–7285. doi:10.1080/21655979.2021.1977551

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fu, J., Zeng, Z., Cohen, D., Li, J., Chen, Q., et al. (2020). TIMER2.0 for Analysis of Tumor-Infiltrating Immune Cells. Nucleic Acids Res. 48 (W1), W509–w514. doi:10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Li, C., Yu, C., Chen, Z., Zhao, C., and Ye, L. (2022). TSPYL2 Reduced Gefitinib Resistance and DNA Damage Repair via Suppressing SIRT1-Mediated FOXO3 Deacetylation. Future Med. Chem. 14 (6), 407–419. doi:10.4155/fmc-2021-0136

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, J., Ma, J., Hao, Z., and Li, W. (2021). HPS6 Regulates the Biogenesis of Weibel-Palade Body in Endothelial Cells through Trafficking V-ATPase to its Limiting Membrane. Front. Cell. Dev. Biol. 9, 743124. doi:10.3389/fcell.2021.743124

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., and Hendershot, L. M. (2004). The Role of the Unfolded Protein Response in Tumour Development: Friend or Foe? Nat. Rev. Cancer 4 (12), 966–977. doi:10.1038/nrc1505

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahadevan, N. R., Anufreichik, V., Rodvold, J. J., Chiu, K. T., Sepulveda, H., and Zanetti, M. (2012). Cell-Extrinsic Effects of Tumor ER Stress Imprint Myeloid Dendritic Cells and Impair CD8+ T Cell Priming. PLoS One 7 (12), e51845. doi:10.1371/journal.pone.0051845

PubMed Abstract | CrossRef Full Text | Google Scholar

Martins, A. S., Alves, I., Helguero, L., Domingues, M. R., and Neves, B. M. (2016). The Unfolded Protein Response in Homeostasis and Modulation of Mammalian Immune Cells. Int. Rev. Immunol. 35 (6), 457–476. doi:10.3109/08830185.2015.1110151

PubMed Abstract | CrossRef Full Text | Google Scholar

Medel, B., Costoya, C., Fernandez, D., Pereda, C., Lladser, A., Sauma, D., et al. (2018). IRE1α Activation in Bone Marrow-Derived Dendritic Cells Modulates Innate Recognition of Melanoma Cells and Favors CD8+ T Cell Priming. Front. Immunol. 9, 3050. doi:10.3389/fimmu.2018.03050

PubMed Abstract | CrossRef Full Text | Google Scholar

Murao, K., Imachi, H., Yu, X., Muraoka, T., Hosami, N., Dobashi, H., et al. (2009). The Transcriptional Factor PREB Mediates MCP-1 Transcription Induced by Cytokines in Human Vascular Endothelial Cells. Atherosclerosis 207 (1), 45–50. doi:10.1016/j.atherosclerosis.2009.03.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, H., Lei, T., Hu, Y., and Lei, P. (2021). Expression of Lipid-Metabolism Genes Is Correlated with Immune Microenvironment and Predicts Prognosis in Osteosarcoma. Front. Cell. Dev. Biol. 9, 673827. doi:10.3389/fcell.2021.673827

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, C., Zhao, F., Zhang, T., Xu, D., Hao, Z., Cui, F., et al. (2022). A Novel Prognostic Signature in Osteosarcoma Characterised from the Perspective of Unfolded Protein Response. Clin. Transl. Med 12 (3), e750. doi:10.1002/ctm2.750

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. U.S.A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, N., Song, W.-X., Luo, J., Haydon, R. C., and He, T.-C. (2008). Osteosarcoma Development and Stem Cell Differentiation. Clin. Orthop. Relat. Res. 466 (9), 2114–2130. doi:10.1007/s11999-008-0335-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, K. P., Fong, S. W., Lu, Z., Ching, Y. P., Chan, K. W., and Chan, S. Y. (2011). TSPYL2 Is Important for G1 Checkpoint Maintenance upon DNA Damage. PLoS One 6 (6), e21602. doi:10.1371/journal.pone.0021602

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor Clelland, C. L., Craciun, L., Bancroft, C., and Lufkin, T. (2000). Mapping and Developmental Expression Analysis of the WD-Repeat Gene Preb. Genomics 63 (3), 391–399. doi:10.1006/geno.1999.6089

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanacker, H., Vetters, J., Moudombi, L., Caux, C., Janssens, S., and Michallet, M.-C. (2017). Emerging Role of the Unfolded Protein Response in Tumor Immunosurveillance. Trends Cancer 3 (7), 491–505. doi:10.1016/j.trecan.2017.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wakamatsu, T., Kakunaga, S., Takenaka, S., Outani, H., Hamada, K., Imura, Y., et al. (2019). Prognostic Implication of Adjuvant/neoadjuvant Chemotherapy Consisting of Doxorubicin and Ifosfamide in Patients with Extraskeletal Osteosarcoma. Int. J. Clin. Oncol. 24 (10), 1311–1319. doi:10.1007/s10147-019-01475-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, Q., Jin, L., and Wang, Z. (2020). Comprehensive Analysis of Cancer Hallmarks in Cutaneous Melanoma and Identification of a Novel Unfolded Protein Response as a Prognostic Signature. aging 12 (20), 20684–20701. doi:10.18632/aging.103974

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, C., Wang, H., Wang, H., Mo, H., Zhong, W., Tang, J., et al. (2020). A Three-Gene Signature Based on Tumour Microenvironment Predicts Overall Survival of Osteosarcoma in Adolescents and Young Adults. Aging 13 (1), 619–645. doi:10.18632/aging.202170

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitton, B., Okamoto, H., Packham, G., and Crabb, S. J. (2018). Vacuolar ATPase as a Potential Therapeutic Target and Mediator of Treatment Resistance in Cancer. Cancer Med. 7 (8), 3800–3811. doi:10.1002/cam4.1594

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Xu, J., Jin, M., Jiang, X., Li, J., Li, X., et al. (2021). Development and Verification of a Hypoxic Gene Signature for Predicting Prognosis, Immune Microenvironment, and Chemosensitivity for Osteosarcoma. Front. Mol. Biosci. 8, 705148. doi:10.3389/fmolb.2021.705148

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, X., Wang, W., Zhang, H., Gao, P., Fan, B., Huang, C., et al. (2015). Individualized Chemotherapy for Osteosarcoma and Identification of Gene Mutations in Osteosarcoma. Tumor Biol. 36 (4), 2427–2435. doi:10.1007/s13277-014-2853-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, M., Ni, J., Song, D., Ding, M., and Huang, J. (2015). Activation of Unfolded Protein Response Protects Osteosarcoma Cells from Cisplatin-Induced Apoptosis through NF-Κb Pathway. Int. J. Clin. Exp. Pathol. 8 (9), 10204–10215.

PubMed Abstract | Google Scholar

Yang, H., Liang, S.-Q., Xu, D., Yang, Z., Marti, T. M., Gao, Y., et al. (2019). HSP90/AXL/eIF4E-regulated Unfolded Protein Response as an Acquired Vulnerability in Drug-Resistant KRAS-Mutant Lung Cancer. Oncogenesis 8 (9), 45. doi:10.1038/s41389-019-0158-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanetti, M., Xian, S., Dosset, M., and Carter, H. (2022). The Unfolded Protein Response at the Tumor-Immune Interface. Front. Immunol. 13, 823157. doi:10.3389/fimmu.2022.823157

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, F., Feng, D., Wang, X., Gu, Y., Shen, Z., Yang, Y., et al. (2021). An Unfolded Protein Response Related Signature Could Robustly Predict Survival Outcomes and Closely Correlate with Response to Immunotherapy and Chemotherapy in Bladder Cancer. Front. Mol. Biosci. 8, 780329. doi:10.3389/fmolb.2021.780329

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Shao, Q., Zhu, H., Xu, H., Long, W., Yu, B., et al. (2018). Resveratrol Ameliorates Lewis Lung Carcinoma‐bearing Mice Development, Decreases Granulocytic Myeloid‐derived Suppressor Cell Accumulation and Impairs its Suppressive Ability. Cancer Sci. 109 (9), 2677–2686. doi:10.1111/cas.13720

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, K., Xiaoqiang, L., Deng, W., Wang, G., and Fu, B. (2021). Identification of a Novel Signature Based on Unfolded Protein Response-Related Gene for Predicting Prognosis in Bladder Cancer. Hum. Genomics 15 (1), 73. doi:10.1186/s40246-021-00372-x

PubMed Abstract | CrossRef Full Text | Google Scholar