Establishing and Validating an Aging-Related Prognostic Signature in Osteosarcoma

Abstract

Aging is an inevitable process that biological changes accumulate with time and results in increased susceptibility to different tumors. But currently, aging-related genes (ARGs) in osteosarcoma were not clear. We investigated the potential prognostic role of ARGs and established an ARG-based prognostic signature for osteosarcoma. The transcriptome data and corresponding clinicopathological information of patients with osteosarcoma were obtained from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Molecular subtypes were generated based on prognosis-related ARGs obtained from univariate Cox analysis. With ARGs, a risk signature was built by univariate, least absolute shrinkage and selection operator (LASSO), and multivariate Cox regression analyses. Differences in clinicopathological features, immune infiltration, immune checkpoints, responsiveness to immunotherapy and chemotherapy, and biological pathways were assessed according to molecular subtypes and the risk signature. Based on risk signature and clinicopathological variables, a nomogram was established and validated. Three molecular subtypes with distinct clinical outcomes were classified based on 36 prognostic ARGs for osteosarcoma. A nine-ARG-based signature in the TCGA cohort, including BMP8A, CORT, SLC17A9, VEGFA, GAL, SSX1, RASGRP2, SDC3, and EVI2B, has been created and developed and could well perform patient stratification into the high- and low-risk groups. There were significant differences in clinicopathological features, immune checkpoints and infiltration, responsiveness to immunotherapy and chemotherapy, cancer stem cell, and biological pathways among the molecular subtypes. The risk signature and metastatic status were identified as independent prognostic factors for osteosarcoma. A nomogram combining ARG-based risk signature and metastatic status was established, showing great prediction accuracy and clinical benefit for osteosarcoma OS. We characterized three ARG-based molecular subtypes with distinct characteristics and built an ARG-based risk signature for osteosarcoma prognosis, which could facilitate prognosis prediction and making personalized treatment in osteosarcoma.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Osteosarcoma is the most common sarcoma mainly occurring in teenagers and young adults around the whole world [1]. It is highly aggressive, and the annual incidence of osteosarcoma is approximately 4.4 per million [2]. Though effective advancement in the treatment and prevention of osteosarcoma has been made over the past decades, overall survival (OS) of osteosarcoma patients is far from satisfactory [3]. Osteosarcoma patients are prone to metastases, and patients with metastatic osteosarcoma have a lower survival rate [4]. More than half of osteosarcoma patients die from metastasis [5]. Meanwhile, osteosarcoma patients with a poor response to chemotherapy usually have an unfavorable prognosis [6]. The application of a risk-adapted strategy to evaluate the therapy response and prognosis ahead of therapy could improve personalized treatment. Therefore, identifying novel prognostic markers for improving the OS of osteosarcoma patients is particularly important.

Aging is characterized by gradual functional deterioration over time and results in increasing susceptibility to a variety of diseases, including cancers and cardiovascular, neurodegenerative, metabolic, and neoplastic diseases [7–10]. Cellular senescence is closely associated with aging [11, 12], and accumulating evidence suggested that senescence cells have a highly complex effect on tumors and can be both beneficial and detrimental. It can not only irreversibly arrest cell growth and inhibit cancer development [13–15] but also have the opposite effect and promote tumor malignancy via the secretion of senescence-associated secretory phenotype factors by the paracrine pathway [16, 17]. In the generation and regulation of cell aging, aging-related genes (ARGs) play a crucial role. It was confirmed that ARGs can not only inhibit tumors by regulating the cellular senescence of tumor cells but also potentially stimulate the initiation, metastasis, and development of tumors [13, 14, 18–20]. It has attracted great attention to identify key ARG characteristics and induce senescence of tumor cells [13]. The diagnostic or prognostic value of the ARG-based signature as biomarkers in malignancy has been widely studied [21, 22]. But the underlying mechanism and prognostic value of ARGs in osteosarcoma remain unclear, and a satisfactory ARG-based prognostic signature for osteosarcoma patients has not been reported.

In this study, by taking advantage of the TCGA database, we established molecular subtypes for osteosarcoma on the basis of prognosis-related ARGs and constructed an ARG-based risk signature for OS prediction of osteosarcoma patients. We compared the difference in clinicopathological features, immune infiltration, immune checkpoints, immunotherapy and chemotherapy response, and biological pathways based on molecular subtypes and the risk signature. The prognostic significance of the ARG-based prognostic signature was validated in GEO cohorts, and a nomogram composed of some clinicopathological factors and the risk signature was generated for providing an accurate prediction of osteosarcoma prognosis.

2. Materials and Methods

2.1. Data Collection and Preparation

Gene expression files and matched clinicopathological data of osteosarcoma patients were retrieved from TCGA and GEO. We removed cases without complete clinical information, follow-up of shorter than 30 days, or cases without status information. 85 patients from the TCGA cohort were retained as a training set, and 86 samples from the GSE21257 (53 cases) and GSE16091 (33 cases) cohorts were utilized for validation. The median value was taken as the gene expression value when a gene ID corresponded to multiple probes in the GEO cohorts or when multiple gene symbols existed in the TCGA cohort. From the Human Aging Genomic Resources (https://genomics.senescence.info/), a total of 307 human ARGs were obtained (Supplementary Table 1).

2.2. Consensus Clustering

The univariate Cox regression analysis filtered prognosis-correlated ARGs, and a heatmap was used to display the correlation among these prognosis-related ARGs. Consensus clustering analysis was conducted to generate ARG-related molecular subtypes using the “ConsensusClusterPlus” R package [23]. To compare the prognosis among the clusters, Kaplan-Meier (K-M) analysis was conducted.

2.3. Developing and Validating the Prognostic ARG Signature

The differentially expressed ARGs (deARGs) among the molecular subtypes were acquired by “limma” package under the screening conditions as the $\text{false}\text{discovery}\text{rate}\text{FDR}<0.05$ and $\log 2\text{fold}\text{change}\text{FC}>\log 22$ [24]. Then, prognosis-related ARGs were identified from the deARGs by univariate Cox regression analysis in the TCGA cohort. LASSO Cox regression analysis was applied to shrink the ARG number. Finally, to establish the prognostic signature, multivariate Cox regression analysis was then implemented. The formula for a risk signature was as follows: $$\begin{array}{}\text{(1)}& \text{Risk}\text{Score}=\sum _{k=0}^n{\beta }_i\times {\text{Exp}}_i,\end{array}$$

where ${\beta }_i$ represents the coefficient and $\text{Ex}{\text{p}}_i$ represents the normalized expression level of a gene. Osteosarcoma cases in the TCGA cohort and GEO cohort were divided into the high- and low-risk groups, according to the optimal cut-off point value of the risk score determined by maximally selected rank statistics using the “maxstat” R package (https://cran.r-project.org/web/packages/maxstat/index.html); here, a significant difference in prognosis between the two groups was detected. OS differences between the risk groups were compared by K-M analysis. The accuracy of the risk signature prediction was estimated by time-dependent ROC analysis in two sets (the validation and training sets).

2.4. Immune Infiltration, Chemotherapeutic Sensitivity, and Immunotherapy Response Predictions

The number of immune cells in the TCGA cohort was calculated using the xCell algorithm (https://xcell.ucsf.edu/) [25] that conducted cell type enrichment analysis for 64 immune cells based on gene expression data. The “estimate” R package was applied to estimate and extrapolate immune and stromal cell fraction in tumor samples [26]. In different groups, immune checkpoint expressions were compared. According to clusters and risks, to predict the clinical response to immune checkpoint inhibitors, the Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) algorithm was applied [27]. For predicting the chemosensitivity of osteosarcoma to several common anticancer drugs (methotrexate, paclitaxel, cisplatin, and doxorubicin), the “pRRophetic” R package was employed for determining half-maximal inhibitory concentration (IC₅₀) values from osteosarcoma gene expression levels [28].

2.5. Gene Set Enrichment Analysis (GSEA)

To analyze the potential differences of pathways on the basis of molecular subtypes and ARG-related risk groups, GSEA was performed using hallmark gene sets as reference. The “GSVA” R package was used to carry out GSEA and pathways with the $\text{FDR}<0.05$ being significantly enriched [29].

2.6. Establishment of a Predictive Nomogram

The decision tree model was applied to classify subgroups based on age, gender, metastatic status, and risk score by using the “rpart” R package (https://cran.r-project.org/web/packages/rpart/index.html). The independent prognostic factors of OS for osteosarcoma were identified by multivariate Cox regression analysis. By the “rms” R package (https://cran.r-project.org/web/packages/rms/index.html), we developed a nomogram integrating independent prognostic clinicopathological factors and the risk signature in the TCGA cohort. To evaluate the prediction accuracy between the actual observations and the predicted 1-, 3-, and 5-year OS probabilities, calibration curves were utilized. Time-dependent ROC curves assessed the nomogram discriminate ability. Decision curve analysis (DCA) tested the clinical applicability of the nomogram using “rmda” R package [30].

2.7. mRNA Expression-Based Stemness Indices (mRNAsi)

Based on the one-class logistic regression (OCLR) algorithm, the stemness index model trained from the Progenitor Cell Biology Consortium database was used to calculate tumor stemness. The stemness index can be used to measure how similar tumor cells are to stem cells, with the stemness index being a value between 0 (lowest) and 1 (highest). The closer the stemness index is to 1, the stronger the stem cell properties. We calculated transcriptome feature scores for the cohorts using the same Spearman correlation.

2.8. Statistical Analysis

The R software (v3.6.3) performed all of the statistical studies. The correlation matrices were conducted using the Pearson or Spearman correlation. The Wilcoxon test was conducted for two-group comparisons. Survival differences were compared using K-M curves with a log-rank test. A $P$ value < 0.05 was considered statistically significant.

3. Results

3.1. Molecular Subtypes of ARG in Osteosarcoma

We firstly identified 36 ARGs which were associated with osteosarcoma prognosis (Figure 1(a)). According to the expression profiles of the 36 ARGs, TCGA osteosarcoma patients were classified into C1, C2, and C3 with distinct outcomes (Figures 1(b)–1(d)). Those in the C3 subtype showed a longer survival than C1 and C2 subtypes ($P<0.0001$, Figure 1(e)). Additionally, the expression level of the 36 prognosis-related ARGs is illustrated in Figure 1(f). Twelve genes belong to the “risk” group and are overall highly expressed in the C1 subtype, while the rest 24 genes were commonly upregulated in the C3 subtype and regarded as the “protective” genes. Moreover, as shown in Figure 1(g), osteosarcoma patients in the C1 subtype had higher tumor metastasis and mortality rate, followed by C2 and C3. Those data showed that 3 subtypes had clinical significance and may provide value for clinical diagnosis.

[figure(s) omitted; refer to PDF]

3.2. Differences in Immune Infiltration, Immunotherapy and Chemotherapy Response, and Biological Pathways among Molecular Subtypes

The relative abundance of 64 immune cells in the TCGA cohort was assessed. Our data demonstrated that the estimated proportion of 40.625% (26/64) of immune cells were significantly different among the three subtypes (Figure 2(a)). The C3 subtype showed a significantly higher immune score than the C1 and C2 subtypes, indicating the highest immune infiltration of the C3 subtype (Figure 2(b)). Meanwhile, there were significant differences of the expression in the immune checkpoint-related gene among different subtypes, as shown in Figure 2(c). Subsequently, the TIDE algorithm was applied to assess the potential response to immunotherapy in the three molecular subtypes. As shown in Figure 2(d), we found that the C2 subtype has a higher TIDE score than C2 and C3 subtypes, indicating that osteosarcoma patients in the C1 and C3 clusters would be more reactive to immunotherapy compared with those in the C2 cluster. The C2 subtype was characterized by a significantly higher T cell exclusion and cancer-associated fibroblast (CAF) score than C1 and C3 subtypes, while the C3 subtype presented significantly higher myeloid-derived suppressor cells (MDSCs) and tumor-associated macrophages (TAMs), as well as a lower T cell dysfunction score. In addition, osteosarcoma patients in the C1 subtype showed a significantly lower IC₅₀ to cisplatin than those in the C2 and C3 subtypes, whereas the C2 subtype showed a significantly lower IC₅₀ to doxorubicin and paclitaxel than C1 and C3 subtypes, indicating that osteosarcoma patients in C1 were more sensitive to cisplatin, while patients in the C2 subtype were more sensitive to doxorubicin and paclitaxel (Figure 2(e)). Moreover, C1 had high mRNAsi (Supplementary Figure 1A).

[figure(s) omitted; refer to PDF]

3.3. Biological Pathways among Different Molecular Clusters

GSEA was performed to reveal the potential biological functions of the genes in the three molecular clusters. The results showed that the enriched pathways in these three subtypes are closely related to immunity. As shown in Figure 3(a), there are 31, 24, and 24 pathways enriched in the C1, C2, and C3 subtypes, respectively. In the C1 subtype, the enriched terms mainly consisted of suppressed immune-related pathways and activated cell cycle-related pathways. The immune-related pathways enriched in the C2 subtype were commonly suppressed, while the HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_GLYCOLYSIS, HALLMARK_HYPOXIA, and HALLMARK_UV_RESPONSE_DN were activated. Most of the enriched immune-related pathways in the C3 subtype were activated. Additionally, the expression of HIPPO, MYC, and RAS was also significantly different among the three clusters (Figure 3(b)), which manifested as a higher expression level in the C2 subtype than the C1 and C3 subtypes.

[figure(s) omitted; refer to PDF]

3.4. Establishment and Evaluation of a Prognostic ARG-Based Risk Signature

495 differentially expressed ARGs were obtained by intersecting C1 vs. C2, C1 vs. C3, and C2 vs. C3, using the “limma” package. A total of 104 ARGs ($P<0.01$) which contribute to osteosarcoma prognosis were screened from the 495 ARGs using univariate Cox regression analysis (Figure 4(a)). According to the results of the LASSO analysis, nineteen prognostic ARGs were further screened out based on the optimal lambda value (0.0594) (Figures 4(b) and 4(c)). Ultimately, we identified 9 ARGs in the risk score by multivariate Cox regression analysis (Figure 4(d)). $\text{Risk}\text{score}=+0.506^{\ast }\text{CORT}-0.714^{\ast }\text{RASGRP}2-0.822^{\ast }\text{SDC}3+0.728^{\ast }\text{BMP}8\text{A}+0.271^{\ast }\text{GAL}+0.408^{\ast }\text{SLC}17\text{A}9-1.{04}^{\ast }\text{EVI}2\text{B}+0.{37}^{\ast }\text{VEGFA}-0.365^{\ast }\text{SSX}1$. Based on the risk score, osteosarcoma patients were successfully separated into the high- and low-risk groups (Figure 5(a)). Our results showed that patients with high risk exhibited shorter OS ($P<0.0001$; Figure 5(c)) than those with low risk. The predictive efficacy of the risk signature was verified by the ROC curve (Figure 5(b)), with the AUC value for 1-, 3-, and 5-year OS being 0.87, 0.92, and 0.9, respectively. This was further validated in the GSE21257 and GSE16091 cohorts, as shown in Figures 5(d)–5(g). In addition, we determine the mRNAsi between the high group and the low group and found that there was no significance between the high group and the low group in the TCGA, GSE21257, and GSE16091 cohorts (Supplementary Figure 1B-D).

[figure(s) omitted; refer to PDF]

3.5. Correlation of Clinicopathological Characteristics with Risk Signature

We analyzed the relationship between clinicopathological features and the ARG-based risk model. The results showed that the risk score was significantly different in groups classified by age, metastasis, status, and molecular clusters (Figure 6(a)). Meanwhile, the high-risk group was characterized by a significantly higher metastasis, mortality rate, and proportion of C1 subtype than the low-risk group (Figure 6(b)). There was no significant difference between the high- and low-risk groups with regard to age and gender. Stratification analysis was conducted according to age, gender, and metastasis, and K-M survival analyses revealed that high-risk patients had significantly unfavorable survival outcomes compared with low-risk patients (Figure 6(c)). It suggested that the prognostic value of risk signature was applicable to other clinical features.

[figure(s) omitted; refer to PDF]

3.6. Immune Infiltration and Biological Pathways between Low- and High-Risk Groups

As shown in Figure 7(a), eleven of the 64 tumor immune cell types showed significantly different abundance between the high- and low-risk groups. In the low-risk group, the multinucleated variant endothelial cells, macrophages M1, macrophages, lymphatic endothelial cells, lentivirus-induced dendritic cells (iDC), hematopoietic stem cells (HSC), fibroblasts, endothelial cells, and chondrocytes had higher infiltration levels than the high-risk group, while plasma cells and melanocytes in the high-risk group were significantly higher compared with the low-risk group in the TCGA cohort. The correlations between risk score and immune cell infiltration were heatmapped in Figure 7(b). It showed that the risk score was significantly negatively correlated with the infiltration level of M1 macrophages, lymphatic endothelial cells, and HSC. Additionally, compared with the high-risk group, patients in the low-risk group showed a higher stromal score and immune score, as well as a higher ESTIMATE score (Figure 7(c)). Biological pathways with a correlation greater than 0.3 are illustrated in Figure 7(d), and the majority of biological pathways were negatively correlated with the risk score. These findings indicated that the low group had higher immune infiltration.

[figure(s) omitted; refer to PDF]

3.7. Risk Signature Predicts Chemotherapy and Immunotherapy Response

We assessed the practicability of the risk signature in guiding systemic therapies in TCGA cohort. The expression of 11 immune checkpoints was significantly different between the high- and low-risk groups, including CD200R1, CD274, CD48, HAVCR2, HHLA2, LAG3, LAIR1, LGALS9, TMIGD2, TNFSF14, and TNFSF9 (Figure 8(a)). As shown in Figure 8(b), patients in the high-risk group had a higher MDSC score than those in the low-risk group. Nonetheless, there was no significant difference between the two risk groups with regard to the CAF score, TAM.M2 score, T cell exclusion score, T cell dysfunction score, and TIDE score. The risk score showed significant positive correlations with the MDSC score ($r=0.45,P=1.6e-5$), TAM.M2 score ($r=0.24,P=0.025$), and T cell dysfunction score ($r=0.44,P=2.3e-5$), as shown in Figure 8(c). We also found that the estimated IC₅₀ values of doxorubicin were significantly higher in the high-risk group compared with that in the low-risk group ($P=0.0072$). However, there was no significant difference between the risk score and the IC₅₀ of methotrexate, paclitaxel, and cisplatin (Figure 8(d)). These findings indicated that a high risk score was related to elevated doxorubicin sensitivity.

[figure(s) omitted; refer to PDF]

3.8. Construction and Validation of a Nomogram

We constructed a decision tree based on age, gender, metastasis, and risk score of osteosarcoma patients in the TCGA cohort (Figure 9(a)), and the results showed that the patients can be stratified into four distinct groups (lowest, low, mediate, and high) using a decision tree on only risk score, gender, and metastatic status. K-M survival analysis showed that there was a significant difference in OS among the four groups (Figure 9(b)). Patients in the “lowest” and “low” subgroups belong to the ARG-based low-risk group (Figure 9(c)). In addition, C1 and C2 subtypes occupy more than the C3 subtype in the “highest” group (Figure 9(d)). Multivariate Cox regression analysis demonstrated that the risk score was the most significant independent prognostic factor of osteosarcoma ($\text{HR}=9.68$, 95% CI: 4.73-19.80, $P=5.41e-10$), followed by metastatic status ($\text{HR}=2.88$, 95% CI: 1.29-6.47, $P=0.0102$) (Figures 9(e) and 9(f)). Therefore, a nomogram is then generated using the risk score and metastatic status to predict the OS of osteosarcoma patients (Figure 9(g)). The calibration plot demonstrated that the nomogram can effectively forecast the actual survival outcomes (Figure 9(h)). Moreover, the DCA curve and time-ROC analysis demonstrated that the nomogram and risk signature had better prognostic capacity than other clinicopathological features, as shown in Figures 9(i) and 9(j).

[figure(s) omitted; refer to PDF]

4. Discussion

Osteosarcoma is the most common bone sarcoma with high heterogeneity and has various subtypes based on morphological and molecular characteristics [31]. Herein, we established three ARG-based patterns with distinct diversity for osteosarcoma patients. The cellular constitutions of stromal and immune cells in the tumor microenvironment (TME) are involved in osteosarcoma progression, chemotherapy resistance, and immunosuppressive activities [32]. In this study, the C3 subtype had a distinctly higher stromal and immune score than C1 and C2, implying the potential differences among the clusters with regard to the progression and response to chemotherapy and immunotherapy. Previous evidence suggested that patients with higher stromal or immune scores had a favorable OS of osteosarcoma [33]. Consistently, the C3 subtype exhibited prolonged survival time and lower mortality and metastasis rates. Although advances in chemotherapy improved the prognosis of osteosarcoma patients [34], drug resistance would result in worse clinical outcomes [35]. Here, patients in the C2 cluster showed more responsiveness to paclitaxel and doxorubicin, while osteosarcoma patients in the C1 subtype were higher reactive to cisplatin. These data suggested that these patients were more likely to benefit from these chemotherapy drugs.

T cell immune checkpoint molecules serve as promising immunotherapeutic targets for tumor, and their inhibitors have dramatically changed the therapeutic landscape of osteosarcoma patients with metastasis or recurrence [36]. Nevertheless, few parts of osteosarcoma patients could benefit from immunotherapies [37]. C1 and C3 subtypes presented increased infiltration levels of many immune cells and higher TIDE scores. Meanwhile, our GSEA results revealed that immune-related pathways were significantly inhibited in the C1 subtype and activated in the C3 subtype. These findings suggested that patients in C1 and C2 subtypes were more likely to respond to immunotherapies. Furthermore, our data found significant differences among the molecular subtypes with regard to M2-TAMs, MDSCs, CAFs, and cytotoxic T lymphocytes (CTLs), which contribute to tumor growth, invasion, and metastasis and mediate the communication between malignant cells, immune cells, and stromal cells [38].

With the development of omic technology and public databases, a variety of risk signatures for prognosis prediction in osteosarcoma has been established on the basis of clinicopathological factors and omic characteristics. Accumulating evidence has elucidated the potential capability of radiomics [39], DNA methylation [40], and immune-related genes [41] for predicting the prognosis of osteosarcoma. Recently, ARG-based prognostic models for various cancers have been established and displayed satisfactory performance for prognosis prediction [42, 43]. At present, the prognostic role of the ARGs in osteosarcoma is still unclear, and an ARG-based risk model has not been developed. Here, we constructed and verified an ARG-based risk signature, which possessed a good performance in the prediction of clinical outcomes of osteosarcoma patients. ARGs can not only inhibit tumor growth but also promote tumor invasion and metastasis [13]. In our ARG-based risk model, osteosarcoma patients with a high risk score were characterized by higher metastasis and poorer prognosis, indicating that ARGs in osteosarcoma were closely associated with cancer progression. Further, we constructed a predictive nomogram integrating the risk score and metastatic status for osteosarcoma, which presented good calibration and discrimination and had adequate ability to predict survival outcomes in patients with osteosarcoma.

Accumulating studies have reported the associations between the ARGs in our risk signature and osteosarcoma. Bone morphogenetic protein 8a (BMP8A), a number of the bone morphogenetic protein ligand, encodes a secreted ligand of the TGF-β superfamily of proteins. It has been proved to promote survival and drug resistance in clear cell renal cell carcinoma (ccRCC) [44], but studies about its effect on osteosarcoma are rare. Cortistatin (CORT) encodes a neuropeptide that is structurally similar to somatostatin. A recent report showed that CORT could inhibit the proliferation of the human thyroid carcinoma cell line, indicating a possible inhibitory role of CORT in cancer development [45]. Solute carrier family 17 member 9 (SLC17A9) was involved in the progression of colorectal cancer (CRC) and breast cancer and was recognized as a potential biomarker for outcome prediction of CRC and BC patients [46, 47]. Vascular endothelial growth factor A (VEGFA) is a potent angiogenic factor for blood vessel formation, and research has proved that VEGFA participates in the angiogenesis and progression of osteosarcoma [48]. Galectin-1 (GAL-1) functions in tissue development, cell proliferation, and immunoregulation. The expression of galectin-1 was capable to be used to differentiate small-cell osteosarcoma from Ewing sarcoma [49]. SSX family member 1 (SSX1) has been associated with stem cell migration, suggesting a potential biologically important role in the metastatic phenotype [50]. RAS guanyl releasing protein 2 (RASGRP2) encodes a brain-enriched nucleotide exchanged factor, and research demonstrated that abnormal expression of RASGRP2 in lung adenocarcinoma correlated with the infiltration level of immune cells [51]. Syndecan 3 (SDC3) encodes a protein that belongs to the syndecan proteoglycan family. Upregulation of SDC3 contributes to perineural invasion and poor outcomes in pancreatic cancer [52]. A recent study has confirmed that the ecotropic viral integration site 2B (EVI2B) could be used as a new prognostic biomarker for metastatic melanoma [53]. However, the functions and underlying mechanisms of these ARGs have not been exhaustively investigated in osteosarcoma.

Cellular senescence and immune infiltration in the TME were proved to contribute to the response of immunotherapy [54], but the correlations between immune infiltration and senescence in osteosarcoma remain poorly understood. Our data revealed that osteosarcoma patients with low risk scores showed an apparent increase in immune cell infiltration, which manifests as a significantly elevated infiltration of multinucleated variant endothelial cells, macrophages M1, macrophages, lymphatic endothelial cells, iDC, HSC, fibroblasts, endothelial cells, and chondrocytes. Moreover, we found that the risk score was negatively correlated with the immune-related pathways. A higher expression level of ten immune checkpoints was observed in the low-risk group, indicating that the patients with lower risk scores had a higher potential to benefit from immunotherapy. However, osteosarcoma patients in the low-risk group had a lower TNFSF9 expression, which was associated with immune response and cell growth in osteosarcoma [55]. Furthermore, evidence suggested that chemosensitivity was regulated by ARGs in osteosarcoma patients [56]. Our data revealed that osteosarcoma patients in the high-risk group were more reactive to doxorubicin. A recent study showed that BMP8A can diminish chemotherapy sensitivity in ccRCC by promoting Nrf2 phosphorylation and activating TRIM24 [44], implying its potential role in drug resistance in osteosarcoma.

Nevertheless, several limitations in our study should be acknowledged. First, the molecular subtypes and ARG-based risk signature were generated using retrospective data from public databases. Therefore, it should be validated in more prospective and multicenter osteosarcoma cohorts in the future. Second, we only investigated the potential prognostic value of the ARG risk signature, so further studies are required to explore the underlying mechanisms of the signature in the development of osteosarcoma.

5. Conclusions

In summary, this study identified molecular subtypes based on ARGs and developed an ARG-based survival prognostic model for osteosarcoma. The difference in immune landscape, biological functions, drug sensitivity, and immunotherapy response according to molecular subtypes and risk groups were analyzed. A nomogram combining the novel ARGs-based risk model and metastatic status was constructed. It may serve as a clinical tool for making personalized therapeutic treatments and forecasting prognoses for osteosarcoma patients.

Authors’ Contributions

Yibo Ma and Shuo Zheng equally contributed.

References

[1] B. A. Lindsey, J. E. Markel, E. S. Kleinerman, "Osteosarcoma overview," Rheumatology and Therapy, vol. 4 no. 1, pp. 25-43, DOI: 10.1007/s40744-016-0050-2, 2017.

[2] S. Avnet, T. Chano, A. Massa, G. Bonuccelli, S. Lemma, L. Falzetti, G. Grisendi, M. Dominici, N. Baldini, "Acid microenvironment promotes cell survival of human bone sarcoma through the activation of cIAP proteins and NF- κ B pathway," American Journal of Cancer Research, vol. 9 no. 6, pp. 1127-1144, 2019.

[3] N. Rainusso, L. L. Wang, J. T. Yustein, "The adolescent and young adult with cancer: state of the art -- bone tumors," Current Oncology Reports, vol. 15 no. 4, pp. 296-307, DOI: 10.1007/s11912-013-0321-9, 2013.

[4] D. C. Dean, S. Shen, F. J. Hornicek, Z. Duan, "From genomics to metabolomics: emerging metastatic biomarkers in osteosarcoma," Cancer Metastasis Reviews, vol. 37 no. 4, pp. 719-731, DOI: 10.1007/s10555-018-9763-8, 2018.

[5] R. L. Siegel, K. D. Miller, A. Jemal, "Cancer statistics, 2019," CA: a Cancer Journal for Clinicians, vol. 69 no. 1,DOI: 10.3322/caac.21551, 2019.

[6] H. Yuan, Y. Gao, "MicroRNA-1908 is upregulated in human osteosarcoma and regulates cell proliferation and migration by repressing PTEN expression," Oncology Reports, vol. 34 no. 5, pp. 2706-2714, DOI: 10.3892/or.2015.4242, 2015.

[7] T. Shavlakadze, M. Morris, J. Fang, S. X. Wang, J. Zhu, W. Zhou, H. W. Tse, R. Mondragon-Gonzalez, G. Roma, D. J. Glass, "Age-related gene expression signature in rats demonstrate early, late, and linear transcriptional changes from multiple tissues," Cell Reports, vol. 28 no. 12, pp. 3263-3273.e3, DOI: 10.1016/j.celrep.2019.08.043, 2019.

[8] M. Armanios, R. de Cabo, J. Mannick, L. Partridge, J. van Deursen, S. Villeda, "Translational strategies in aging and age-related disease," Nature Medicine, vol. 21 no. 12, pp. 1395-1399, DOI: 10.1038/nm.4004, 2015.

[9] B. A. Benayoun, E. A. Pollina, P. P. Singh, S. Mahmoudi, I. Harel, K. M. Casey, B. W. Dulken, A. Kundaje, A. Brunet, "Remodeling of epigenome and transcriptome landscapes with aging in mice reveals widespread induction of inflammatory responses," Genome Research, vol. 29 no. 4, pp. 697-709, DOI: 10.1101/gr.240093.118, 2019.

[10] K. Jr Smetana, L. Lacina, P. Szabo, B. Dvořánková, P. Brož, A. Šedo, "Ageing as an important risk factor for cancer," Anticancer Research, vol. 36 no. 10, pp. 5009-5018, DOI: 10.21873/anticanres.11069, 2016.

[11] D. Yin, K. Chen, "The essential mechanisms of aging: irreparable damage accumulation of biochemical side-reactions," Experimental Gerontology, vol. 40 no. 6, pp. 455-465, DOI: 10.1016/j.exger.2005.03.012, 2005.

[12] C. López-Otín, M. A. Blasco, L. Partridge, M. Serrano, G. Kroemer, "The hallmarks of aging," Cell, vol. 153 no. 6, pp. 1194-1217, DOI: 10.1016/j.cell.2013.05.039, 2013.

[13] A. Calcinotto, J. Kohli, E. Zagato, L. Pellegrini, M. Demaria, A. Alimonti, "Cellular senescence: aging, cancer, and injury," Physiological Reviews, vol. 99 no. 2, pp. 1047-1078, DOI: 10.1152/physrev.00020.2018, 2019.

[14] L. Mosteiro, C. Pantoja, N. Alcazar, R. M. Marión, D. Chondronasiou, M. Rovira, P. J. Fernandez-Marcos, M. Muñoz-Martin, C. Blanco-Aparicio, J. Pastor, G. Gómez-López, A. De Martino, M. A. Blasco, M. Abad, M. Serrano, "Tissue damage and senescence provide critical signals for cellular reprogramming in vivo," Science, vol. 354 no. 6315,DOI: 10.1126/science.aaf4445, 2016.

[15] S. HeN, E. Sharpless, "Senescence in health and disease," Cell, vol. 169 no. 6, pp. 1000-1011, DOI: 10.1016/j.cell.2017.05.015, 2017.

[16] B. Bottazzi, E. Riboli, A. Mantovani, "Aging, inflammation and cancer," Seminars in Immunology, vol. 40, pp. 74-82, DOI: 10.1016/j.smim.2018.10.011, 2018.

[17] J. P. Coppé, C. K. Patil, F. Rodier, Y. Sun, D. P. Muñoz, J. Goldstein, P. S. Nelson, P. Y. Desprez, J. Campisi, "Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p 53 tumor suppressor," PLoS Biology, vol. 6 no. 12, pp. 2853-2868, DOI: 10.1371/journal.pbio.0060301, 2008.

[18] S. Lee, C. A. Schmitt, "The dynamic nature of senescence in cancer," Nature Cell Biology, vol. 21 no. 1, pp. 94-101, DOI: 10.1038/s41556-018-0249-2, 2019.

[19] S. C. Johnson, P. S. Rabinovitch, M. Kaeberlein, "mTOR is a key modulator of ageing and age-related disease," Nature, vol. 493 no. 7432, pp. 338-345, DOI: 10.1038/nature11861, 2013.

[20] L. Galluzzi, I. Vitale, S. A. Aaronson, J. M. Abrams, D. Adam, P. Agostinis, E. S. Alnemri, L. Altucci, I. Amelio, D. W. Andrews, M. Annicchiarico-Petruzzelli, A. V. Antonov, E. Arama, E. H. Baehrecke, N. A. Barlev, N. G. Bazan, F. Bernassola, M. J. M. Bertrand, K. Bianchi, M. V. Blagosklonny, K. Blomgren, C. Borner, P. Boya, C. Brenner, M. Campanella, E. Candi, D. Carmona-Gutierrez, F. Cecconi, F. K. M. Chan, N. S. Chandel, E. H. Cheng, J. E. Chipuk, J. A. Cidlowski, A. Ciechanover, G. M. Cohen, M. Conrad, J. R. Cubillos-Ruiz, P. E. Czabotar, V. D’Angiolella, T. M. Dawson, V. L. Dawson, V. de Laurenzi, R. de Maria, K. M. Debatin, R. J. DeBerardinis, M. Deshmukh, N. di Daniele, F. di Virgilio, V. M. Dixit, S. J. Dixon, C. S. Duckett, B. D. Dynlacht, W. S. el-Deiry, J. W. Elrod, G. M. Fimia, S. Fulda, A. J. García-Sáez, A. D. Garg, C. Garrido, E. Gavathiotis, P. Golstein, E. Gottlieb, D. R. Green, L. A. Greene, H. Gronemeyer, A. Gross, G. Hajnoczky, J. M. Hardwick, I. S. Harris, M. O. Hengartner, C. Hetz, H. Ichijo, M. Jäättelä, B. Joseph, P. J. Jost, P. P. Juin, W. J. Kaiser, M. Karin, T. Kaufmann, O. Kepp, A. Kimchi, R. N. Kitsis, D. J. Klionsky, R. A. Knight, S. Kumar, S. W. Lee, J. J. Lemasters, B. Levine, A. Linkermann, S. A. Lipton, R. A. Lockshin, C. López-Otín, S. W. Lowe, T. Luedde, E. Lugli, M. MacFarlane, F. Madeo, M. Malewicz, W. Malorni, G. Manic, J. C. Marine, S. J. Martin, J. C. Martinou, J. P. Medema, P. Mehlen, P. Meier, S. Melino, E. A. Miao, J. D. Molkentin, U. M. Moll, C. Muñoz-Pinedo, S. Nagata, G. Nuñez, A. Oberst, M. Oren, M. Overholtzer, M. Pagano, T. Panaretakis, M. Pasparakis, J. M. Penninger, D. M. Pereira, S. Pervaiz, M. E. Peter, M. Piacentini, P. Pinton, J. H. M. Prehn, H. Puthalakath, G. A. Rabinovich, M. Rehm, R. Rizzuto, C. M. P. Rodrigues, D. C. Rubinsztein, T. Rudel, K. M. Ryan, E. Sayan, L. Scorrano, F. Shao, Y. Shi, J. Silke, H. U. Simon, A. Sistigu, B. R. Stockwell, A. Strasser, G. Szabadkai, S. W. G. Tait, D. Tang, N. Tavernarakis, A. Thorburn, Y. Tsujimoto, B. Turk, T. vanden Berghe, P. Vandenabeele, M. G. Vander Heiden, A. Villunger, H. W. Virgin, K. H. Vousden, D. Vucic, E. F. Wagner, H. Walczak, D. Wallach, Y. Wang, J. A. Wells, W. Wood, J. Yuan, Z. Zakeri, B. Zhivotovsky, L. Zitvogel, G. Melino, G. Kroemer, "Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018," Cell Death and Differentiation, vol. 25 no. 3, pp. 486-541, DOI: 10.1038/s41418-017-0012-4, 2018.

[21] N. Zeng, C. Guo, Y. Wang, L. Li, X. Chen, S. Gao, F. Jiang, B. Cao, "Discovery of lipid metabolism-related genes for predicting tumor immune microenvironment status and prognosis in prostate cancer," Journal of Oncology, vol. 2022,DOI: 10.1155/2022/8227806, 2022.

[22] L. Zheng, Y. Yang, X. Cui, "Establishing and validating an aging-related prognostic four-gene signature in colon adenocarcinoma," BioMed Research International, vol. 2021,DOI: 10.1155/2021/4682589, 2021.

[23] M. D. Wilkerson, D. N. Hayes, "ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking," Bioinformatics, vol. 26 no. 12, pp. 1572-1573, DOI: 10.1093/bioinformatics/btq170, 2010.

[24] M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, G. K. Smyth, "limma powers differential expression analyses for RNA-sequencing and microarray studies," Nucleic Acids Research, vol. 43 no. 7,DOI: 10.1093/nar/gkv007, 2015.

[25] D. Aran, Z. HuA, J. Butte, "xCell: digitally portraying the tissue cellular heterogeneity landscape," Genome Biology, vol. 18 no. 1,DOI: 10.1186/s13059-017-1349-1, 2017.

[26] K. Yoshihara, M. Shahmoradgoli, E. Martínez, R. Vegesna, H. Kim, W. Torres-Garcia, V. Treviño, H. Shen, P. W. Laird, D. A. Levine, S. L. Carter, G. Getz, K. Stemke-Hale, G. B. Mills, R. G. Verhaak, "Inferring tumour purity and stromal and immune cell admixture from expression data," Nature Communications, vol. 4 no. 1,DOI: 10.1038/ncomms3612, 2013.

[27] P. Jiang, S. Gu, D. Pan, J. Fu, A. Sahu, X. Hu, Z. Li, N. Traugh, X. Bu, B. Li, J. Liu, G. J. Freeman, M. A. Brown, K. W. Wucherpfennig, X. S. Liu, "Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response," Nature Medicine, vol. 24 no. 10, pp. 1550-1558, DOI: 10.1038/s41591-018-0136-1, 2018.

[28] P. Geeleher, N. Cox, R. S. Huang, "pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels," PLoS One, vol. 9 no. 9,DOI: 10.1371/journal.pone.0107468, 2014.

[29] S. Hänzelmann, R. Castelo, J. Guinney, "GSVA: gene set variation analysis for microarray and RNA-seq data," BMC Bioinformatics, vol. 14 no. 1,DOI: 10.1186/1471-2105-14-7, 2013.

[30] X. Yang, H. Sun, Y. Song, L. Yang, H. Liu, "Diagnostic and prognostic values of upregulated SPC25 in patients with hepatocellular carcinoma," Peer J., vol. 8,DOI: 10.7717/peerj.9535, 2020.

[31] K. Schiavone, D. Garnier, M. F. Heymann, D. Heymann, "The heterogeneity of osteosarcoma: the role played by cancer stem cells," Advances in Experimental Medicine and Biology, vol. 1139, pp. 187-200, DOI: 10.1007/978-3-030-14366-4_11, 2019.

[32] S. Smeland, S. S. Bielack, J. Whelan, M. Bernstein, P. Hogendoorn, M. D. Krailo, R. Gorlick, K. A. Janeway, F. C. Ingleby, J. Anninga, I. Antal, C. Arndt, K. L. B. Brown, T. Butterfass-Bahloul, G. Calaminus, M. Capra, C. Dhooge, M. Eriksson, A. M. Flanagan, G. Friedel, M. C. Gebhardt, H. Gelderblom, R. Goldsby, H. E. Grier, R. Grimer, D. S. Hawkins, S. Hecker-Nolting, K. Sundby Hall, M. S. Isakoff, G. Jovic, T. Kühne, L. Kager, T. von Kalle, E. Kabickova, S. Lang, C. C. Lau, P. J. Leavey, S. L. Lessnick, L. Mascarenhas, R. Mayer-Steinacker, P. A. Meyers, R. Nagarajan, R. L. Randall, P. Reichardt, M. Renard, C. Rechnitzer, C. L. Schwartz, S. Strauss, L. Teot, B. Timmermann, M. R. Sydes, N. Marina, "Survival and prognosis with osteosarcoma: outcomes in more than 2000 patients in the EURAMOS-1 (European and American Osteosarcoma Study) cohort," European Journal of Cancer, vol. 109, pp. 36-50, DOI: 10.1016/j.ejca.2018.11.027, 2019.

[33] P. M. Alves, J. A. A. de Arruda, D. A. C. Arantes, S. F. S. Costa, L. L. Souza, H. A. R. Pontes, F. P. Fonseca, R. A. Mesquita, C. F. W. Nonaka, E. F. Mendonça, A. C. Batista, "Evaluation of tumor-infiltrating lymphocytes in osteosarcomas of the jaws: a multicenter study," Virchows Archiv, vol. 474 no. 2, pp. 201-207, DOI: 10.1007/s00428-018-2499-6, 2019.

[34] C. Zhang, J. H. Zheng, Z. H. Lin, H. Y. Lv, Z. M. Ye, Y. P. Chen, X. Y. Zhang, "Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of osteosarcoma," Aging, vol. 12 no. 4, pp. 3486-3501, DOI: 10.18632/aging.102824, 2020.

[35] X. Xiao, W. Wang, Y. Li, D. Yang, X. Li, C. Shen, Y. Liu, X. Ke, S. Guo, Z. Guo, "HSP90AA1-mediated autophagy promotes drug resistance in osteosarcoma," Journal of Experimental & Clinical Cancer Research, vol. 37 no. 1,DOI: 10.1186/s13046-018-0880-6, 2018.

[36] C. Chen, L. Xie, T. Ren, Y. Huang, J. Xu, W. Guo, "Immunotherapy for osteosarcoma: fundamental mechanism, rationale, and recent breakthroughs," Cancer Letters, vol. 500,DOI: 10.1016/j.canlet.2020.12.024, 2021.

[37] W. Yu, Y. Wang, J. Zhu, L. Jin, B. Liu, K. Xia, J. Wang, J. Gao, C. Liang, H. Tao, "Autophagy inhibitor enhance ZnPc/BSA nanoparticle induced photodynamic therapy by suppressing PD-L1 expression in osteosarcoma immunotherapy," Biomaterials, vol. 192, pp. 128-139, DOI: 10.1016/j.biomaterials.2018.11.019, 2019.

[38] A. Ferguson, K. Chen, "Analysis of transcriptional profiling of immune cells at the single-cell level," Methods in Molecular Biology, vol. 2111, pp. 47-57, DOI: 10.1007/978-1-0716-0266-9_4, 2020.

[39] S. Zhao, Y. Su, J. Duan, Q. Qiu, X. Ge, A. Wang, Y. Yin, "Radiomics signature extracted from diffusion-weighted magnetic resonance imaging predicts outcomes in osteosarcoma," Journal of Bone Oncology, vol. 19,DOI: 10.1016/j.jbo.2019.100263, 2019.

[40] X. Zhang, Y. Zheng, G. Li, C. Yu, T. Ji, S. Miao, "Identifying four DNA methylation gene sites signature for predicting prognosis of osteosarcoma," Translational Cancer Research, vol. 9 no. 11, pp. 7299-7309, DOI: 10.21037/tcr-20-3204, 2020.

[41] W. Liu, X. Xie, Y. Qi, J. Wu, "Exploration of immune-related gene expression in osteosarcoma and association with outcomes," JAMA Network Open, vol. 4 no. 8,DOI: 10.1001/jamanetworkopen.2021.19132, 2021.

[42] Q. Xu, Y. Chen, "An aging-related gene signature-based model for risk stratification and prognosis prediction in lung adenocarcinoma," Frontiers in Cell and Development Biology, vol. 9,DOI: 10.3389/fcell.2021.685379, 2021.

[43] T. Yue, S. Chen, J. Zhu, S. Guo, Z. Huang, P. Wang, S. Zuo, Y. Liu, "The aging-related risk signature in colorectal cancer," Aging, vol. 13 no. 5, pp. 7330-7349, DOI: 10.18632/aging.202589, 2021.

[44] Y. P. Yu, L. C. Cai, X. Y. Wang, S. Y. Cheng, D. M. Zhang, W. G. Jian, T. D. Wang, J. K. Yang, K. B. Yang, C. Zhang, "BMP8A promotes survival and drug resistance via Nrf 2/TRIM24 signaling pathway in clear cell renal cell carcinoma," Cancer Science, vol. 111 no. 5, pp. 1555-1566, DOI: 10.1111/cas.14376, 2020.

[45] P. Cassoni, G. Muccioli, T. Marrocco, M. Volante, E. Allia, E. Ghigo, R. Deghenghi, M. Papotti, "Cortistatin-14 inhibits cell proliferation of human thyroid carcinoma cell lines of both follicular and parafollicular origin," Journal of Endocrinological Investigation, vol. 25 no. 4, pp. 362-368, DOI: 10.1007/BF03344019, 2002.

[46] L. Yang, Z. Chen, W. Xiong, H. Ren, E. Zhai, K. Xu, H. Yang, Z. Zhang, L. Ding, Y. He, X. Song, J. Liu, "High expression of SLC17A9 correlates with poor prognosis in colorectal cancer," Human Pathology, vol. 84, pp. 62-70, DOI: 10.1016/j.humpath.2018.09.002, 2019.

[47] J. Li, T. Su, L. Yang, L. Deng, C. Zhang, Y. He, "High SLC17A9 expression correlates with poor survival in gastric carcinoma," Future Oncology, vol. 15 no. 36, pp. 4155-4166, DOI: 10.2217/fon-2019-0283, 2019.

[48] L. Zhang, Z. Lv, J. Xu, C. Chen, Q. Ge, P. Li, D. Wei, Z. Wu, X. Sun, "MicroRNA-134 inhibits osteosarcoma angiogenesis and proliferation by targeting the VEGFA/VEGFR1 pathway," The FEBS Journal, vol. 285 no. 7, pp. 1359-1371, DOI: 10.1111/febs.14416, 2018.

[49] I. Machado, G. J. A. López, S. Navarro, E. Mayordomo, K. Scotlandi, P. Picci, A. Llombart-Bosch, "Galectin-1 (GAL-1) expression is a useful tool to differentiate between small cell osteosarcoma and Ewing sarcoma," Virchows Archiv, vol. 462 no. 6, pp. 665-671, DOI: 10.1007/s00428-013-1423-3, 2013.

[50] G. Cronwright, K. Le Blanc, C. Götherström, P. Darcy, M. Ehnman, B. Brodin, "Cancer/testis antigen expression in human mesenchymal stem cells: down-regulation of SSX impairs cell migration and matrix metalloproteinase 2 expression," Cancer Research, vol. 65 no. 6, pp. 2207-2215, DOI: 10.1158/0008-5472.CAN-04-1882, 2005.

[51] S. Zhao, X. Jin, S. Xu, "Expression of RASGRP2 in lung adenocarcinoma and its effect on immune microenvironment," Zhongguo Fei Ai Za Zhi, vol. 24 no. 6, pp. 404-411, DOI: 10.3779/j.issn.1009-3419.2021.105.01, 2021.

[52] N. A. Afratis, D. Nikitovic, H. A. Multhaupt, A. D. Theocharis, J. R. Couchman, N. K. Karamanos, "Syndecans - key regulators of cell signaling and biological functions," The FEBS Journal, vol. 284 no. 1, pp. 27-41, DOI: 10.1111/febs.13940, 2017.

[53] S. Yonekura, K. Ueda, "EVI2B is a new prognostic biomarker in metastatic melanoma with IFNgamma associated immune infiltration," Cancers, vol. 13 no. 16,DOI: 10.3390/cancers13164110, 2021.

[54] W. Lin, X. Wang, Z. Xu, Z. Wang, T. Liu, Z. Cao, X. Feng, Y. Gao, J. He, "Identification and validation of cellular senescence patterns to predict clinical outcomes and immunotherapeutic responses in lung adenocarcinoma," Cancer Cell International, vol. 21 no. 1,DOI: 10.1186/s12935-021-02358-0, 2021.

[55] Y. Tsuda, C. Tanikawa, T. Miyamoto, M. Hirata, V. Yodsurang, Y. Z. Zhang, S. Imoto, R. Yamaguchi, S. Miyano, H. Takayanagi, H. Kawano, H. Nakagawa, S. Tanaka, K. Matsuda, "Identification of a p53 target, CD137L , that mediates growth suppression and immune response of osteosarcoma cells," Scientific Reports, vol. 7 no. 1,DOI: 10.1038/s41598-017-11208-x, 2017.

[56] N. Oka, A. Komuro, H. Amano, S. Dash, M. Honda, K. Ota, S. Nishimura, T. Ueda, M. Akagi, H. Okada, "Ascorbate sensitizes human osteosarcoma cells to the cytostatic effects of cisplatin," Pharmacology Research & Perspectives, vol. 8 no. 4,DOI: 10.1002/prp2.632, 2020.