Comprehensive analysis of hypoxia-related genes for prognosis value, immune status, and therapy in osteosarcoma patients

Osteosarcoma is a common malignant bone tumor in children and adolescents. The overall survival of osteosarcoma patients is remarkably poor. Herein, we sought to establish a reliable risk prognostic model to predict the prognosis of osteosarcoma patients. Patients ’ RNA expression and corresponding clinical data were downloaded from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) and Gene Expression Omnibus databases. A consensus clustering was conducted to uncover novel molecular subgroups based on 200 hypoxia-linked genes. A hypoxia-risk models were established by Cox regression analysis coupled with LASSO regression. Functional enrichment analysis, including Gene Ontology annotation and KEGG pathway analysis, were conducted to determine the associated mechanisms. Moreover, we explored relationships between the risk scores and age, gender, tumor microenvironment, and drug sensitivity by correlation analysis. We identified two molecular subgroups with significantly different survival rates and developed a risk model based on 12 genes. Survival analysis indicated that the high-risk osteosarcoma patients likely have a poor prognosis. The area under the curve (AUC) value showed the validity of our risk scoring model, and the nomogram indicates the model’s reliability. High-risk patients had lower Tfh cell infiltration and a lower stromal score. We determined the abnormal expression of three prognostic genes in osteosarcoma cells. Sunitinib can promote osteosarcoma cell apoptosis with down-regulation of KCNJ3 expression. In summary, the constructed hypoxia-related risk score model can assist clinicians during clinical practice for osteosarcoma prognosis management. Immune and drug sensitivity analysis can provide essential insights into subsequent mechanisms. KCNJ3 may be a valuable prognostic marker for osteosarcoma development.

Introduction

Osteosarcoma is widely acknowledged as the most common primary bone tumor, especially in adolescents (.8–1.1 per 100,000 per year at age 15–19 years) (Whelan et al., 2012; Valery et al., 2015). Osteosarcoma tends to occur in the metaphysis of long bones, with a high metastasis and recurrence rate. It is well-established that the lungs represent the most common metastatic site, posing a severe threat to the health of adolescents (Casali et al., 2018). Surgery, radiotherapy (RT), and multimodal chemotherapy (ChT) remain the mainstay of treatment for osteosarcoma patients (Casali et al., 2018). Despite recent advances in medicine, the prognosis of osteosarcoma is still poor. It has been reported that the aggregate event-free survival rate of patients with recurrent, unresectable osteosarcoma was only 12% at 4 months (Lagmay et al., 2016). In addition, due to genetic heterogeneity, patients with the same clinical manifestations and pathological types often have different prognoses after receiving the same treatment scheme (Wu et al., 2020). Therefore, a reliable new prognostic model is pivotal to ensuring targeted therapy and improving osteosarcoma treatment efficacy, which aligns with the concept of precision medicine.

Over the past years, experts have attached great importance to hypoxia, a feature of the tumor immune microenvironment. The rapid proliferation of cancer cells, desmoplastic fibrotic stroma, and non-functional angiogenesis lead to increased oxygen consumption and decreased oxygen supply (Gola et al., 2020; Tao et al., 2021a). Ample evidence has confirmed that tumor hypoxia has an indispensable effect on apoptosis, cell proliferation, vascularization, immune responses, metabolism, genomic instability, and metastasis (Wigerup et al., 2016). Moreover, hypoxia promotes dihydropyrimidine dehydrogenase expression in macrophages, resulting in chemotherapy resistance in colorectal cancer (Malier et al., 2021). It is generally believed that bone is prone to hypoxia. An increasing body of studies has shown that hypoxia in the tumor microenvironment is a driver for resistance to chemotherapy, tumor growth, cell survival, and metastasis in numerous solid tumors, such as osteosarcoma (Prudowsky and Yustein, 2020; Jiang et al., 2021). Studies have shown that hypoxia can induce HIF-1α, NUSAP1, and NDUFA4L2 expression and promotes the survival, epithelial-mesenchymal transition progression, migration, and invasion of osteosarcoma cells (Xu et al., 2020; Zhang et al., 2021). Hence, we should concentrate on differences in hypoxia genes to examine the biological mechanisms of osteosarcoma and reclassify tumor subtypes.

Hypoxia plays an indispensable role in tumor development and the anti-tumor process. However, the relationship between hypoxia and osteosarcoma prognosis remains further unexplored. We conducted a comprehensive analysis of hypoxia-related genes in osteosarcoma to understand the role of hypoxia in tumor genesis and development. Herein, we retrieved genes associated with hypoxia from the MSigDB database (Zhou et al., 2021). Consensus clustering analysis is a popular method for identifying new tumor molecular subtypes (Zuccato et al., 2022). Importantly, we used hypoxia gene expression as the input of consensus clustering to place a unique molecular subtype of osteosarcoma. We established a model for estimating patient prognostic risk based on this subtype. Furthermore, we analyzed the immune microenvironment and chemosensitivity of osteosarcoma. Finally, we further studied the expression and role of some predictive genes in osteosarcoma in vitro.

Materials and methods

Data collection

84 osteosarcoma samples, containing RNA-Seq data and corresponding clinical information, were downloaded from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database (https://ocg.cancer.gov/programs/target). The clinical characteristics of the samples are shown in Supplementary Table S1. The osteosarcoma dataset GSE21257 was obtained from the GEO database (https://www.ncbi.nlm.gov/geo) as a validation cohort. The hypoxia-linked genes (shown in Supplementary Table S2) were retrieved from the hallmark gene sets in the Molecular Signature Database3 (MSigDB).

Cluster analysis based on hypoxia-related genes

The R package “ConsensusClusterPlus” was performed to cluster the samples based on mRNA expression levels of 200 hypoxia-related genes. The lowest intergroup correlation and the highest intragroup correlation were shown when clustering variable (k) = 2. R packages “survival” and “survminer” (v3.6.1) were used to analyze the overall survival time between the two subtypes. The clinical characteristics and related genes between the two clusters were displayed in a heatmap using the R package.

Development and validation of a hypoxia subtype-related gene prognostic model

23 differentially expressed genes (DEGs) between the hypoxia-related subtypes were identified for downstream analyses using the screening criteria |log2FC| ≥ 2 and FDR < .05. We uncovered genes related to individual survival and prognosis of osteosarcoma through Cox regression analysis in the TARGET cohort. A p-value of .05 was set as the threshold. LASSO (Least absolute shrinkage and selector operation) Cox regression was conducted to uncover significant prognostic genes for prognostic model using the R package “glmnet”. Ultimately, we retained the 12 genes coupled with their coefficients, and the penalty parameter (λ) was determined via the minimum criteria. The risk score formula was calculated as follows:

$RiskScore={\displaystyle \sum }_{i=1}^{n}Coe{f}_i\hspace{0.17em}\ast \hspace{0.17em}{x}_i$

According to the median risk score, osteosarcoma patients from the TARGET cohort were stratified into low- and high-risk groups. Based on the 12-gene signature, we performed the PCA (principal component analysis) and t-SNE (t-distributed stochastic neighbor embedding). The Kaplan-Meier analysis was conducted to compare survival possibility and overall survival time between the high- and low-risk groups. The R package “timeROC” was used to draw the time-based ROC curve. The area under the curve (AUC) was calculated to assess the sensitivity and specificity of the risk score system. Additionally, we underwent external validation for the gene signature model in the GSE21257.

Independent prognostic analysis and a nomogram

We applied the univariate and multivariable Cox regression analysis (“survival” R package) to assess the risk score combined with clinical information (age, gender and metastatic status) in the TARGET cohort. Combined with prognostic signature, a nomogram was presented to predict 1, 3, and 5-year overall Survival (OS) of osteosarcoma patients.

Tumor immune microenvironment and functional enrichment analysis

We employed the single sample gene set enrichment analysis (ssGSEA) in low- and high-risk groups to explore the infiltrating scores of immune cells and immune pathways. Benjamini-Hochberg (BH) correction method was used to calculate the adjusted p-value. We utilized the estimate algorithm to calculate the infiltration levels of immune and stromal cells. We applied Spearman correlation analysis to analyze the relationships between risk score and immune and stromal cells. Screening of DEGs between the low- and high-risk groups was done using the criteria |log2FC|≥1 and FDR < .05. The package “clusterprofiler” was used for Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

Drug sensitivity analysis

We downloaded the NCI-60 human cancer cell lines from the CellMiner database (https://discover. nci.nih.gov/cellminer) (Luna et al., 2021). Then Pearson correlation analysis was conducted to determine the correlation between prognostic genes and drug sensitivity.

Cell lines and cultures

We purchased one human osteoblast cell line (hFOB1.19) and two human osteosarcoma cell lines (U20S and 143B) from the National Collection of Authenticated Cell Cultures (Shanghai, China). Dulbecco’s modified Eagle’s medium (DMEM, Gibco) contains 1% penicillin/streptomycin (Thermo Fisher Scientific) and 10% fetal bovine serum (FBS). We cultured hFOB1.19 cells in DMEM at 34°C with 5% CO₂ and U20S and 143B cells at 37°C with 5% CO₂.

Treatment and cell viability analysis

Osteosarcoma cells were treated with a range of concentrations of sunitinib (0, 10, 20, and 30 μM). According to operational guidelines, we used the Cell Counting Kit-8 (CCK-8) assay (MedChemExpress; NJ, the United States) to evaluate cell viability. Briefly, the second-generation osteosarcoma cells were grown on 96 well plates (5 × 10³ cells/well) and exposed to different concentrations of sunitinib. Osteosarcoma cells were incubated in DMEM at 37°C with 5% CO₂ for 24 h. Next, we rinsed cells with PBS and added 10 μL CCK8 to the medium for another 2 h. Finally, a microwell plate (Thermo Fisher) reader detected the absorbance at 450 nm.

TUNEL staining

Apoptotic DNA fragmentations of osteosarcoma cells were detected with One-step TUNEL in Situ Apoptosis Kit (E-CK-A321, Elabscience, Wuhan, China). Osteosarcoma cells were treated with sunitinib (0, 10, 20, and 30 μM) for 24 h. Following the manufacturer’s guidelines, cells were fixed in 4% paraformaldehyde for 30 min at room temperature and then incubated with 3% H2O2 and .2% Triton X-100 for 10 min. Next, we washed with PBS three times and stained the nuclei with DAPI. Finally, we used an Olympus fluorescence microscope (Tokyo, Japan) to observe apoptotic cells (TUNEL-positive cells).

Western blotting

We extracted total cellular protein from osteosarcoma cells using RIPA lysis buffer with 1 mM PMSF and determined their concentrations with the BCA Protein Assay Kit (Beyotime, Shanghai, China). Next, proteins were separated by PAGE and transferred to PVDF membranes. After blocking with skim milk (5%, w/v) for 2 h at room temperature, membranes were incubated with the following primary antibodies (1: 1,000) overnight at 4°C: RASGRP2 (ABclonal, A15381), KCNJ3 (Abcam, ab129182), ACTG2 (Abcam, ab231802), CASP3 (Proteintech, 19677-1-AP), Bcl-2 (Proteintech, 26593-1-AP), Bax (Proteintech, 50599-2-Ig), and GAPDH (ABclonal, AC001). Membranes were then incubated with HRP-labeled IgG secondary antibody (1:2000, Beyotime, A02080) for 2 h at 25°C. Protein bands on the membrane were visualized using the ECL Plus kit (Meilunbio). Finally, the band intensity was quantified via Image Lab 6.1 software (BioRad, Hercules, CA, USA).

Statistical analysis

R software (Version:3.6.1) and GraphPad Prism (Version:7.00) were conducted to perform all statistical analyses. Pearson chi-square test, t-test, the Mann-Whitney test, and analysis of variance (ANOVA) were used to calculate and compare the two groups. a p-value < .05 was statistically significant.

Results

Tumor classification based on the hypoxia-related genes

In the TARGET dataset, consensus clustering was conducted for all osteosarcoma patients to study the relationship between the expression of 200 hypoxia-related genes and osteosarcoma subtypes. After setting the clustering variable (k) value of 2, the patients could be stratified into two groups (Figure 1A). We found a significant difference in the overall survival rate between the two subtypes (p = .029, Figure 1B), and the cluster 2 had a significantly better overall survival than cluster 1. The relationship between the expression of these genes and the clinical characteristics, including gender, age, metastasis status, and primary tumor site, is displayed in a heatmap. However, no significant differences in clinical characteristics were observed between the two clusters (Figure 1C).

FIGURE 1

FIGURE 1. Tumour classification based on the hypoxia-related genes. (A) 84 osteosarcoma patients were grouped into two clusters according to the consensus clustering matrix (k = 2). (B) Kaplan–Meier overall survival curves for the two clusters. (C) A heatmap (blue: low expression level; red: high expression level) for the connections between clinicopathologic features and the clusters.

Construction of a prognostic gene model in the TARGET cohort

DEGs were identified between the two clusters using the screening criteria |log2FC| ≥ 2 and FDR < .05. Univariate Cox regression was used for the initial determination of survival-linked genes. 23 genes that met the criteria of p < .05 were retained for subsequent analyses. Of these, 4 genes (FAP, GFPT2, POSTN, ACTG2) were protective genes with an HR < 1, while the remaining 19 genes were associated with increased risk (HRs > 1) (Figure 2A). LASSO Cox regression was used to create a 12-gene signature based on the optimum λ value (Figures 2B, C). Computation of the risk score was done using the formula: risk score= (−.175* CYFIP2 exp.) + (−.503* RASGRP2 exp.) + (−.320* DKK1exp.) + (−.645* DLX2 exp.) + (−.130* GFPT2exp.) + (−.670* KCNJ3 exp.) + (−.019* ACTG2exp.) + (−.204* CHMP4C exp.) + (−.544* KLK1 exp.) + (−.588* NRXN1 exp.) + (−.305* ABCA4 exp.) + (−.734* CORT exp.). According to the median risk score, the 84 osteosarcoma patients were stratified into high-and low-risk subgroups (Figure 2D). PCA and the t-SNE analysis exhibited that patients with high or low risks were well separated into two classes (Figures 2E, F). In contrast with the high-risk group, the low-risk group experienced fewer deaths and exhibited longer survival (Figure 2G, left side of dotted line). The Kaplan-Meier curve showed that overall survival time and possibility were significantly lower in the high-risk group (p < .001, Figure 2H). Then ROC analysis was applied to determine the sensitivity and specificity of the prognostic model. We found that the value of the area under the ROC curve (AUC) was .749 for 1-year, .786 for 3-year, and .797 for 5-year survival prediction (Figure 2I).

FIGURE 2

FIGURE 2. Construction of risk signature in the TARGET cohort. (A) Univariate cox regression analysis of overall survival for each hypoxia-related gene, and 23 genes with p < .05. (B) LASSO regression of the 12 overall survival-related genes. (C) Cross-validation for tuning the parameter selection in the LASSO regression. (D) Distribution of patients based on the risk score. (E) PCA plot for osteosarcoma based on the risk score. (F) The t-SNE analysis based on the risk score. (G) The survival status for each patient (left side of the dotted line: low-risk population; right side of the dotted line: high-risk population). (H) Kaplan–Meier curves for the overall survival of patients between the high- and low-risk groups. (I) ROC curves demonstrated the predictive efficiency of the risk score.

External validation of the risk score

53 osteosarcoma patients from a GEO cohort (GSE21257) were extracted as the external validation set. The gene expression data were normalized with the R function “Scale”. Based on the median risk score of the TARGET model, 28 patients were regarded as the high-risk group, while the other 25 were at low risk (Figure 3A). High-risk patients had higher mortality and shorter overall survival time (Figure 3B on the right of the dotted line). The PCA plot and the t-SNE results showed that the two groups were separated (Figures 3C, D). Besides, there were significant differences in survival time between the high-OS and low-OS risk groups (p = .012, Figure 3E). ROC curve analysis of the GEO dataset showed that our model had an excellent prediction power, with AUC values for survival at 1, 2, and 3 years of .837, .678, and .689, respectively (Figure 3F).

FIGURE 3

FIGURE 3. Validation of the risk model in the GEO cohort. (A) Distribution of patients in the GEO cohort based on the median risk score of the TARGET cohort. (B) The survival status for each patient (left side of the dotted line: low-risk population; right side of the dotted line: high-risk population). (C) The PCA plot for osteosarcoma. (D) The t-SNE analysis based on the risk score. (E) Kaplan–Meier curves for comparison of the overall survival between low- and high-risk groups. (F) Time-dependent ROC curves for osteosarcoma.

Independent prognostic value of the risk model

We employed univariate and multivariable Cox regression to assess independent prognostic factors of the gene-based risk score and clinical characteristics. The univariate Cox regression data illustrated that the risk score and metastatic status were independent predictors of poor survival in the TARGET cohort (p < .001, Figure 4A). After correcting for confounding factors, the results of multivariate cox regression indicated that risk score (HR = 4.665, p < .001) and metastatic status (HR = 3.596, p < .001) were still independent predictors of the prognosis of osteosarcoma patients (Figure 4B). Moreover, we provided a heatmap to show gene expression and clinical characteristics between low- and high-risk groups (Figure 4C). Two genes (GFPT2 and ACTG2) were downregulated, while the remaining genes were upregulated in the high-risk subgroup.

FIGURE 4

FIGURE 4. Independence detection of the constructed risk prediction model. (A) Univariate analysis for the TARGET cohort. (B) Multivariate analysis for the TARGET cohort. (C) A heatmap (blue: low expression; red: high expression) for the connections between clinicopathologic features and the risk groups.

Nomogram

According to the prognostic model and clinical factors (age, gender, and metastatic status), we developed a risk estimation nomogram in the TARGET cohort (Figure 5A). 1-, 3-, and 5-year calibration curves showed that the nomogram consistently predicted the survival rate (Figure 5B).

FIGURE 5

FIGURE 5. Construction of the predictive model. (A) A prognostic model to predict overall survival in the TARGET cohort. (B) Calibration curves of the OS nomogram model in the TARGET set.

Immune status and tumor microenvironment

We calculated these two clusters’ enrichment scores of diverse immune cell sub-populations by ssGSEA. Significant differences in aDCs cells and Tfh cells were found between the high and low-risk groups (Figure 6A). Meanwhile, we indicated that the immune score of APC_co_stimulation was decreased in high-risk group, compared to the low-risk group (Figure 6B, p < .05). Moreover, correlation analysis illustrated that the risk score was positively correlated to stromal score (p < .001, Figure 6D), while no significant relationship was observed between the immune score and risk score (p > .05, Figure 6C).

FIGURE 6

FIGURE 6. Immune status between different risk groups and the association between risk score and tumor microenvironment. (A) Comparison of the enrichment scores of 16 types of immune cells between low- (blue box) and high-risk (red box) group in the TARGET cohort. (B) Comparison of the enrichment scores of 13 types of immune functions between low- (blue box) and high-risk (red box) group in the TARGET cohort. (C) The relationship between risk score and immune score. (D) The relationship between risk score and stromal score. (*p < .05).

Functional analyses of the risk model

To further assess differences in gene functions and cascades between the sub-clusters classified via the risk model, the R package “limma” was used to screen DEGs with an FDR < .05 along with |log2FC |≥1. We identified 73 DEGs between low-OS and high-OS risk groups in the Target cohort. Of these, 50 genes were upregulated, and 23 genes were downregulated in the high-risk group (Supplementary Table S3). These DEGs underwent GO and KEGG analyses. We found that biological processes linked to hypoxia (i.e., organization of extracellular matrix, organization of extracellular structure, and organization of external encapsulating structure) were significantly enriched. (Figures 7A, B).

FIGURE 7

FIGURE 7. Functional analysis based on the DEGs between the two-risk groups in the TARGET cohort. (A) Bubble graph for GO enrichment (the bigger bubble means the more genes enriched, and the increasing depth of red means the differences were more obvious; q-value: the adjusted p-value). (B) Barplot graph for KEGG pathways (the longer bar means the more genes enriched, and the increasing depth of red means the differences were more obvious).

Sensitivity to chemotherapy

We further investigated the sensitivity of 12 predictive hypoxia-related genes to chemotherapeutic drugs. At the same time, we downloaded data from the NCI-60 panel of human cancer cell lines. Results for the top 16 correlation analysis were shown based on the p-value (Figure 8). KCNJ3 is sensitive to multiple chemotherapeutic agents, including LOXO-101, NMS-E628, and sunitinib (all p < .001). Results showed that RASGRP2 is sensitive to nelarabine, hydroxyurea, asparaginase, cyclophosphamide, fludarabine, pipobroman, chlorambucil, and cladribine (all p < .001). Besides, ACTG2 is susceptible to zoledronate (p < .001). Please see Supplementary Table S4 for details.

FIGURE 8

FIGURE 8. Scatter plot of relationship between prognostic gene expression and drug sensitivity. The top 16 correlation analyses are shown based on the p-value. The horizontal axis represents the gene expression; The vertical axis represents changes in gene expression after administration.

Detection of three predictive hypoxia-related markers in osteosarcoma

We investigate the expression levels of three hypoxia-related genes in osteosarcoma cells using western blotting analysis. Results showed higher expression levels of RASGRP2 and KCNJ3 in two osteosarcoma cell groups (U20S and 143B) compared to the osteoblast cell group (hFOB). In contrast, ACTG2 was down-regulated in osteosarcoma groups (Figures 9A, B). Such results validated the abnormal expression of three prognostic genes in osteosarcoma.

FIGURE 9

FIGURE 9. The expression levels of three hypoxia-related genes between osteosarcoma cell lines and osteoblasts. (A) Western blotting of the expressions of RASGRP2, KCNJ3, and ACTG2 in hFOB, U20S, and 143B groups. GAPDH serves as an internal standard. The gels have been run under the same experimental conditions. (B) A histogram of the OD values of RASGRP2, KCNJ3, and ACTG2 in each group (n = 3 per group). The obtained data are represented as mean ± SE. Significance: ***p-value < .001, vs. hFOB group.

Sunitinib down-regulated KCNJ3 expression and promoted apoptosis in osteosarcoma cells

To evaluate the therapeutic effect of sunitinib on osteosarcoma, we treated two types of osteosarcoma cells (U20S and 143B) with various doses of sunitinib (10, 20, or 30 μM). Sunitinib demonstrated a cytotoxic effect on osteosarcoma cells in a dose-dependent manner (Figure 10A; Figure 11A). Western blotting was used to determine the level of KCNJ3 in U20S and 143B osteosarcoma cells with various doses of sunitinib. Results indicated the expression of KCNJ3 was down-regulated after sunitinib treatment (Figures 10B, C; Figures 11B, C). Western blot results demonstrated that the expression levels of Bax and CASP3 were enhanced with the increasing sunitinib concentration, while the expression of Bcl-2 decreased (Figures 10D, E); (Figures 11D, E). Additionally, TUNEL was used to assess sunitinib treatment for osteosarcoma cells. As expected, the number of TUNEL-positive cells was upregulated with the increase of sunitinib concentration (Figures 10F, G).

FIGURE 10

FIGURE 10. Sunitinib can decrease KCNJ3 expression and enhance apoptosis in U20S osteosarcoma cells. (A) Evaluation of U20S osteosarcoma cell viability using CCK-8 assay after exposure to various concentrations of sunitinib for 24 h. (B,C) The expression level of KCNJ3 protein of osteosarcoma cells in control and sunitinib treatment groups. (D,E) The protein expression levels of Bax, Bcl-2, and caspase3 in osteosarcoma cells in the sham (0μM) and sunitinib treatment groups. (F,G) TUNEL staining was used to detect osteosarcoma cell apoptosis after sunitinib treatment (bar: 50 μm; nuclei: blue; positive cells: green). All experiments were repeated in triplicates (n = 3). The obtained data are represented as mean ± SE. Significance: ^@ p-value < .05, ^@@ p-value < .01, ^@@@ p-value < .001, vs. sham (0 μM) group. ^# p-value < .05, ^## p-value < .01, ^### p-value < .001.

FIGURE 11

FIGURE 11. Sunitinib promotes apoptosis in 143B osteosarcoma cells with decreased expression of KCNJ3. (A) 143B osteosarcoma cell viability evaluation via CCK-8 assay. (B,C) The expression level of KCNJ3 protein. (D,E) The protein expression levels of Bcl-2, Bax, and caspase3. All experiments were repeated in triplicates (n = 3). The obtained data are represented as mean ± SE. Significance: ^@@@ p-value < .001, vs. sham (0 μM) group. ^# p-value < .05, ^## p-value < .01, ^### p-value < .001.

Discussion

Osteosarcoma is the most common primary tumor with early metastasis and rapid progression, accounting for the poor survival rates (Simpson et al., 2017). Despite remarkable medical advancements achieved in recent years, little progress has been made in improving osteosarcoma patients’ prognosis, considering that the current clinical approach is still based on conventional methods and is often ineffective (Ouyang et al., 2020). Precision medicine is widely acknowledged as a medical model that adjusts disease prevention and treatment methods according to the individual differences of each patient. In recent years, with the rapid development of tumor molecular biology and genomics, the understanding of tumor molecular phenotype has been improving. Accordingly, the development of targeted therapy based on specific molecular phenotypes has become the preferred treatment method for advanced cancer, accounting for the increase in the momentum of precision medicine for clinical tumor treatment (Konig et al., 2017). High throughput sequencing is the cornerstone of precision medicine. A comprehensive understanding of tumor molecular phenotype is significant for rapid and accurate diagnosis, treatment efficacy, and prognosis prediction (Heather and Chain, 2016). An increasing body of evidence shows that hypoxia is indispensable in promoting survival, epithelial-mesenchymal transition (EMT) progression, invasion of osteosarcoma cells, and drug resistance (Prudowsky and Yustein, 2020; Fu et al., 2021; Li et al., 2022a). Therefore, in this study, patients with osteosarcoma were classified based on the expression of hypoxia-related genes.

Based on these 200 hypoxia-related genes, osteosarcoma patients were divided into C1 and C2 types to determine the differentially expressed prognostic genes. To further evaluate the prognostic significance of these genes, univariate Cox regression and LASSO Cox regression analysis were applied to establish a 12-gene risk signature. The good performance of the signature was validated in an external dataset. CYFIP2 is well-recognized as a p53-driven pro-apoptotic protein with low expression levels in gastric cancer (Jackson et al., 2007; Jiao et al., 2017). Similar to the conclusions of He Y et al. (He et al., 2021), we believe that RASGRP2 is a risk factor for the prognosis of osteosarcoma patients. RASGRP2 can inhibit apoptosis by activating Rap1 to down-regulate the production of TNF-induced ROS (Takino et al., 2019). Another study demonstrated that DKK1 upregulation could accelerate the deterioration of bone microstructure related to the occurrence of femoral head necrosis and osteosarcoma (Chen et al., 2021). DLX2, a homeobox transcription factor, plays an indispensable role in tumor progression and metastasis in TGF-β-exposed cancer cells (Park et al., 2021). Overexpression of DLX2 has been documented to promote osteoblast differentiation (Zeng et al., 2020). Moreover, it has been reported that GFPT2 is involved in the developing colorectal cancer (Liu et al., 2020). Besides, the potassium channel gene KCNJ3 has been upregulated in non-small cell lung cancer, pancreatic cancer, and breast cancer. KCNJ3, a gene encoding G-protein activated inwardly rectifying K (+) channel (GIRK1), is related to lymph node metastasis and prognosis in breast cancer patients (Kammerer et al., 2016; Rezania et al., 2016). In addition, the differential expression of ACTG2 may play an essential role in the development of osteosarcoma (Lauvrak et al., 2013). Moreover, CHMP4C promotes the viability and motility of cervical cancer cells by regulating epithelial-mesenchymal transition (Lin et al., 2020). KLK1, a member of the Kalinin gene family, can increase extracellular matrix degradation and enhance tumor cells’ survival, proliferation, and invasion (Avgeris et al., 2012). Interestingly, the cell adhesion molecule neurexin-1 (NRXN1) is considered a new potential target for tumor cells (Li et al., 2022b; Feng et al., 2022). ABCA4 has been extensively studied in retinal diseases (Scortecci et al., 2021), with few studies on tumors warranting further research. We performed immunoblot in vitro for RASGRP2, KCNJ3, and ACTG2 genes to further verify our results’ reliability. The expression levels of RASGRP2 and KCNJ3 were upregulated in osteosarcoma cells, while the expression of ACTG2 was decreased.

Over the years, few studies have focused on immunotherapy and immune cell infiltration for osteosarcoma. Accordingly, exploring the immune cell infiltration and immune status in osteosarcoma is crucial to better understanding the molecular mechanism. By comparing the immune cell infiltration in the low-and high-risk groups, we found that the infiltration level of aDCs in the high-risk group was significantly higher. In comparison, the infiltration level of Tfh cells was considerably lower. Consistently, Tao Z et al. found that the expression of Tfh in thymoma tissues decreased while the expression of aDCs increased (Tao et al., 2021b). Indeed, the primary function of Tfh cells is to participate in information transmission during the process of B-cell differentiation and assist in B-cell activation. Tfh cells are abundant in inflammatory infiltrates of breast cancer, and the presence of Tfh cells can improve survival and reduce immunosuppression (Gu-Trantien et al., 2013). Besides, regulating Tfh cell infiltration may provide a therapeutic way for colorectal cancer treatment (Overacre-Delgoffe et al., 2021). Accordingly, we believe that high infiltration levels of Tfh cells in tumor tissues may benefit osteosarcoma patients’ prognosis. Overall, these results indicate that the tumor immune microenvironment plays an important role in tumor therapy.

We further conducted GO and KEGG pathway enrichment analysis to understand the possible biological functions of these prognostic genes. Our results illustrated that these genes were remarkably abundant in the extracellular matrix and extracellular structures, indicating a significant correlation between hypoxia-related genes and osteosarcoma prognosis. Subsequently, we evaluated gene therapy targets’ sensitivity to drugs for future treatment. Based on data analysis of sixty different cell lines, the elevated expressions of 12 prognostic genes influence sensitivity or resistance to chemotherapy drugs approved by the Food and Drug Administration. For example, with the increased KCNJ3 expression, cancer cells were sensitive to LOXO-101, NMS-E628, and sunitinib. CORT up-regulation is a useful biomarker for cancer cells sensitive to ifosfamide. At present, clinicians generally use ifosfamide in osteosarcoma chemotherapy. Therefore, these findings may provide a new perspective for precision treatment for osteosarcoma.

We found that KCNJ3 expression in tumor cells significantly decreased after administering sunitinib through drug sensitivity analysis. KCNJ3 expression was reduced in osteosarcoma cells (U20S and 143B) in a dose-dependent manner of sunitinib. Meanwhile, sunitinib significantly promoted osteosarcoma cell apoptosis. One study reported that sunitinib could reduce PD-L1 expression and remodel the immune system, thus inhibiting the migration and invasion of osteosarcoma cells (Duan et al., 2020). Another study showed that sunitinib can activate immune cells for sarcoma treatment (Ocadlikova et al., 2021). Local implantation of sunitinib and chlorin e6 can promote 143B osteosarcoma cell apoptosis (Yu et al., 2020). We believed sunitinib might provide better treatment for osteosarcoma with elevated KCNJ3 expression.

Nevertheless, there are some limitations to our study. Firstly, our risk score model was constructed and validated based on public databases. Therefore, we need a multicenter extensive sample survey to evaluate the clinical application of our model. Moreover, the mechanism of these 12 predictive genes isn’t precise in osteosarcoma. We will conduct comprehensive functional experiments and multi-omics analysis in future research.

Conclusion

In conclusion, we uncovered a novel molecular subgroup via consensus clustering in the present study. Subsequently, we established a prognostic model and nomogram for osteosarcoma patients to predict the prognosis based on 12 hypoxia-related DEGs. The prognostic signature reflected the different risk groups’ immune characteristics and chemotherapeutic drug sensitivity. Our model may guide the clinical application and patients’ prognostic management. Of note, KCNJ3 may be a reliable prognostic biomarker for clinical application.

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.

Author contributions

HY contributed to conception and design of the study. TH, ZW, ZZ, JL, and CX designed, performed the experiments, conducted the data analysis and wrote the manuscript. HY revised the manuscript. All authors contributed to read and approved the submitted version manuscript version.

Funding

This study was supported by grants from the Training Program for High-level Innovative Health Talents in Zhejiang Province.

Acknowledgments

We like to acknowledge the TARGET and the GEO (GSE21257) network for providing data.

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/fphar.2022.1088732/full#supplementary-material

References

Avgeris, M., Mavridis, K., and Scorilas, A. (2012). Kallikrein-related peptidases in prostate, breast, and ovarian cancers: From pathobiology to clinical relevance. Biol. Chem. 393 (5), 301–317. doi:10.1515/hsz-2011-0260

PubMed Abstract | CrossRef Full Text | Google Scholar

Casali, P. G., Bielack, S., Abecassis, N., Aro, H. T., Bauer, S., Biagini, R., et al. (2018). Bone sarcomas: ESMO-PaedCan-EURACAN clinical practice guidelines for diagnosis, treatment and follow-up. Ann. Oncol. 29 (4), iv79–iv95. doi:10.1093/annonc/mdy310

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Yang, C., Yang, Y., Liang, Q., Xie, K., Liu, J., et al. (2021). Targeting DKK1 prevents development of alcohol-induced osteonecrosis of the femoral head in rats. Am. J. Transl. Res. 13 (4), 2320–2330.

PubMed Abstract | Google Scholar

Duan, X. L., Guo, J. P., Li, F., Xiu, C., and Wang, H. (2020). Sunitinib inhibits PD-L1 expression in osteosarcoma by targeting STAT3 and remodels the immune system in tumor-bearing mice. Future Oncol. 16 (24), 1815–1824. doi:10.2217/fon-2019-0725

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, P., Li, Z., Li, Y., Zhang, Y., and Miao, X. (2022). Characterization of different subtypes of immune cell infiltration in glioblastoma to aid immunotherapy. Front. Immunol. 13, 799509. doi:10.3389/fimmu.2022.799509

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Bao, Q., Liu, Z., He, G., Wen, J., Liu, Q., et al. (2021). Development and validation of a hypoxia-associated prognostic signature related to osteosarcoma metastasis and immune infiltration. Front. Cell Dev. Biol. 9, 633607. doi:10.3389/fcell.2021.633607

PubMed Abstract | CrossRef Full Text | Google Scholar

Gola, C., Iussich, S., Noury, S., Martano, M., Gattino, F., Morello, E., et al. (2020). Clinical significance and in vitro cellular regulation of hypoxia mimicry on HIF-1α and downstream genes in canine appendicular osteosarcoma. Vet. J. 264, 105538. doi:10.1016/j.tvjl.2020.105538

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu-Trantien, C., Loi, S., Garaud, S., Equeter, C., Libin, M., de Wind, A., et al. (2013). CD4⁺ follicular helper T cell infiltration predicts breast cancer survival. J. Clin. Invest 123 (7), 2873–2892. doi:10.1172/JCI67428

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Y., Zhou, H., Wang, W., Xu, H., and Cheng, H. (2021). Construction of a circRNA-miRNA-mRNA regulatory network reveals potential mechanism and treatment options for osteosarcoma. Front. Genet. 12, 632359. doi:10.3389/fgene.2021.632359

PubMed Abstract | CrossRef Full Text | Google Scholar

Heather, J. M., and Chain, B. (2016). The sequence of sequencers: The history of sequencing DNA. Genomics 107 (1), 1–8. doi:10.1016/j.ygeno.2015.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackson, R. S., Cho, Y. J., Stein, S., and Liang, P. (2007). CYFIP2, a direct p53 target, is leptomycin-B sensitive. Cell Cycle 6 (1), 95–103. doi:10.4161/cc.6.1.3665

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, F., Miao, X. L., Zhang, X. T., Yan, F., Mao, Y., Wu, C. Y., et al. (2021). A hypoxia gene-based signature to predict the survival and affect the tumor immune microenvironment of osteosarcoma in children. J. Immunol. Res. 2021, 5523832. doi:10.1155/2021/5523832

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, S., Li, N., Cai, S., Guo, H., and Wen, Y. (2017). Inhibition of CYFIP2 promotes gastric cancer cell proliferation and chemoresistance to 5-fluorouracil through activation of the Akt signaling pathway. Oncol. Lett. 13 (4), 2133–2140. doi:10.3892/ol.2017.5743

PubMed Abstract | CrossRef Full Text | Google Scholar

Kammerer, S., Jahn, S. W., Winter, E., Eidenhammer, S., Rezania, S., Regitnig, P., et al. (2016). Critical evaluation of KCNJ3 gene product detection in human breast cancer: mRNA in situ hybridisation is superior to immunohistochemistry. J. Clin. Pathol. 69 (12), 1116–1121. doi:10.1136/jclinpath-2016-203798

PubMed Abstract | CrossRef Full Text | Google Scholar

Konig, I. R., Fuchs, O., Hansen, G., von Mutius, E., and Kopp, M. V. (2017). What is precision medicine? Eur. Respir. J. 50 (4), 1700391. doi:10.1183/13993003.00391-2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Lagmay, J. P., Krailo, M. D., Dang, H., Kim, A., Hawkins, D. S., Beaty, O., et al. (2016). Outcome of patients with recurrent osteosarcoma enrolled in seven phase II trials through children's cancer group, pediatric oncology group, and children's oncology group: Learning from the past to move forward. J. Clin. Oncol. 34 (25), 3031–3038. doi:10.1200/JCO.2015.65.5381

PubMed Abstract | CrossRef Full Text | Google Scholar

Lauvrak, S. U., Munthe, E., Kresse, S. H., Stratford, E. W., Namlos, H. M., Meza-Zepeda, L. A., et al. (2013). Functional characterisation of osteosarcoma cell lines and identification of mRNAs and miRNAs associated with aggressive cancer phenotypes. Br. J. Cancer 109 (8), 2228–2236. doi:10.1038/bjc.2013.549

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Mao, Q., Zhang, Z., Wang, Y., Chen, D., Chen, Z., et al. (2022). Identification of a five-gene panel to assess prognosis for gastric cancer. Biomed. Res. Int. 2022, 5593619. doi:10.1155/2022/5593619

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z. F., Meng, D. D., Liu, Y. Y., Bi, F. G., Tian, K., Xu, J. Z., et al. (2022). Hypoxia inducible factor-3α promotes osteosarcoma progression by activating KDM3A-mediated demethylation of SOX9. Chem. Biol. Interact. 351, 109759. doi:10.1016/j.cbi.2021.109759

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, S. L., Wang, M., Cao, Q. Q., and Li, Q. (2020). Chromatin modified protein 4C (CHMP4C) facilitates the malignant development of cervical cancer cells. FEBS Open Bio 10 (7), 1295–1303. doi:10.1002/2211-5463.12880

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Pan, Y., Ren, X., Zeng, Z., Sun, J., Zhou, K., et al. (2020). GFPT2 promotes metastasis and forms a positive feedback loop with p65 in colorectal cancer. Am. J. Cancer Res. 10 (8), 2510–2522.

PubMed Abstract | Google Scholar

Luna, A., Elloumi, F., Varma, S., Wang, Y., Rajapakse, V. N., Aladjem, M. I., et al. (2021). CellMiner Cross-Database (CellMinerCDB) version 1.2: Exploration of patient-derived cancer cell line pharmacogenomics. Nucleic Acids Res. 49 (D1), D1083–D1093. doi:10.1093/nar/gkaa968

PubMed Abstract | CrossRef Full Text | Google Scholar

Malier, M., Gharzeddine, K., Laverriere, M. H., Marsili, S., Thomas, F., Decaens, T., et al. (2021). Hypoxia drives dihydropyrimidine dehydrogenase expression in macrophages and confers chemoresistance in colorectal cancer. Cancer Res. 81 (23), 5963–5976. doi:10.1158/0008-5472.CAN-21-1572

PubMed Abstract | CrossRef Full Text | Google Scholar

Ocadlikova, D., Lecciso, M., Broto, J. M., Scotlandi, K., Cavo, M., Curti, A., et al. (2021). Sunitinib exerts in vitro immunomodulatory activity on sarcomas via dendritic cells and synergizes with PD-1 blockade. Front. Immunol. 12, 577766. doi:10.3389/fimmu.2021.577766

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, Z., Li, G., Zhu, H., Wang, J., Qi, T., Qu, Q., et al. (2020). Construction of a five-super-enhancer-associated-genes prognostic model for osteosarcoma patients. Front. Cell Dev. Biol. 8, 598660. doi:10.3389/fcell.2020.598660

PubMed Abstract | CrossRef Full Text | Google Scholar

Overacre-Delgoffe, A. E., Bumgarner, H. J., Cillo, A. R., Burr, A. H. P., Tometich, J. T., Bhattacharjee, A., et al. (2021). Microbiota-specific T follicular helper cells drive tertiary lymphoid structures and anti-tumor immunity against colorectal cancer. Immunity 54 (12), 2812–2824.e4. doi:10.1016/j.immuni.2021.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, H. R., Choi, Y. J., Kim, J. Y., Kim, I. G., and Jung, U. (2021). Repeated irradiation with gamma-ray induces cancer stemness through TGF-beta-DLX2 signaling in the A549 human lung cancer cell line. Int. J. Mol. Sci. 22 (8), 4284. doi:10.3390/ijms22084284

PubMed Abstract | CrossRef Full Text | Google Scholar

Prudowsky, Z. D., and Yustein, J. T. (2020). Recent insights into therapy resistance in osteosarcoma. Cancers (Basel) 13 (1), 83. doi:10.3390/cancers13010083

PubMed Abstract | CrossRef Full Text | Google Scholar

Rezania, S., Kammerer, S., Li, C., Steinecker-Frohnwieser, B., Gorischek, A., DeVaney, T. T., et al. (2016). Overexpression of KCNJ3 gene splice variants affects vital parameters of the malignant breast cancer cell line MCF-7 in an opposing manner. BMC Cancer 16, 628. doi:10.1186/s12885-016-2664-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Scortecci, J. F., Molday, L. L., Curtis, S. B., Garces, F. A., Panwar, P., Van Petegem, F., et al. (2021). Cryo-EM structures of the ABCA4 importer reveal mechanisms underlying substrate binding and Stargardt disease. Nat. Commun. 12 (1), 5902. doi:10.1038/s41467-021-26161-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Simpson, S., Dunning, M. D., de Brot, S., Grau-Roma, L., Mongan, N. P., and Rutland, C. S. (2017). Comparative review of human and canine osteosarcoma: Morphology, epidemiology, prognosis, treatment and genetics. Acta Vet. Scand. 59 (1), 71. doi:10.1186/s13028-017-0341-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Takino, J. I., Sato, T., Nagamine, K., and Hori, T. (2019). The inhibition of Bax activation-induced apoptosis by RasGRP2 via R-Ras-PI3K-Akt signaling pathway in the endothelial cells. Sci. Rep. 9 (1), 16717. doi:10.1038/s41598-019-53419-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, J., Yang, G., Zhou, W., Qiu, J., Chen, G., Luo, W., et al. (2021). Targeting hypoxic tumor microenvironment in pancreatic cancer. J. Hematol. Oncol. 14 (1), 14. doi:10.1186/s13045-020-01030-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, Z., Lu, C., Gao, S., Zhang, P., Chen, Y., Wang, Y., et al. (2021). Two types of immune infiltrating cells and six hub genes can predict the occurrence of myasthenia gravis in patients with thymoma. Bioengineered 12 (1), 5004–5016. doi:10.1080/21655979.2021.1958634

PubMed Abstract | CrossRef Full Text | Google Scholar

Valery, P. C., Laversanne, M., and Bray, F. (2015). Bone cancer incidence by morphological subtype: A global assessment. Cancer Causes Control 26 (8), 1127–1139. doi:10.1007/s10552-015-0607-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Whelan, J., McTiernan, A., Cooper, N., Wong, Y. K., Francis, M., Vernon, S., et al. (2012). Incidence and survival of malignant bone sarcomas in England 1979-2007. Int. J. Cancer 131 (4), E508–E517. doi:10.1002/ijc.26426

PubMed Abstract | CrossRef Full Text | Google Scholar

Wigerup, C., Pahlman, S., and Bexell, D. (2016). Therapeutic targeting of hypoxia and hypoxia-inducible factors in cancer. Pharmacol. Ther. 164, 152–169. doi:10.1016/j.pharmthera.2016.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z. L., Deng, Y. J., Zhang, G. Z., Ren, E. H., Yuan, W. H., and Xie, Q. Q. (2020). Development of a novel immune-related genes prognostic signature for osteosarcoma. Sci. Rep. 10 (1), 18402. doi:10.1038/s41598-020-75573-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W. N., Yang, R. Z., Zheng, H. L., Jiang, L. S., and Jiang, S. D. (2020). NDUFA4L2 regulated by HIF-1α promotes metastasis and epithelial-mesenchymal transition of osteosarcoma cells through inhibiting ROS production. Front. Cell Dev. Biol. 8, 515051. doi:10.3389/fcell.2020.515051

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Z., Xiao, Z., Shuai, X., and Tian, J. (2020). Local delivery of sunitinib and Ce6 via redox-responsive zwitterionic hydrogels effectively prevents osteosarcoma recurrence. J. Mater Chem. B 8 (30), 6418–6428. doi:10.1039/d0tb00970a

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, X., Wang, Y., Dong, Q., Ma, M. X., and Liu, X. D. (2020). DLX2 activates Wnt1 transcription and mediates Wnt/β-catenin signal to promote osteogenic differentiation of hBMSCs. Gene 744, 144564. doi:10.1016/j.gene.2020.144564

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Song, J., Xin, X., Sun, D., Huang, H., Chen, Y., et al. (2021). Hypoxia stimulates the migration and invasion of osteosarcoma via up-regulating the NUSAP1 expression. Open Med. (Wars) 16 (1), 1083–1089. doi:10.1515/med-2020-0180

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, C., Zhang, H., and Lu, L. (2021). Identification and validation of hypoxia-related lncRNA signature as a prognostic model for hepatocellular carcinoma. Front. Genet. 12, 744113. doi:10.3389/fgene.2021.744113

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuccato, J. A., Patil, V., Mansouri, S., Liu, J. C., Nassiri, F., Mamatjan, Y., et al. (2022). DNA methylation-based prognostic subtypes of chordoma tumors in tissue and plasma. Neuro Oncol. 24 (3), 442–454. doi:10.1093/neuonc/noab235

PubMed Abstract | CrossRef Full Text | Google Scholar