A prognostic glycolysis-related gene signature in osteosarcoma: implications for metabolic programming, immune microenvironment, and drug response

Introduction

Osteosarcoma (OS), an aggressive malignancy originating in the bone or cartilage that predominantly affects children and adolescents, is third leading cause of malignancy in this age group (Ritter & Bielack, 2010). The clinical manifestations include local pain, swelling, functional impairment, and pathological fractures, which are often accompanied by metastases to the lungs and other sites (Chen et al., 2021). The development of OS involves multifaceted genetic and environmental factors, including chromosomal abnormalities, growth factors, and cytokines; however, the precise mechanisms remain unclear. Currently, diagnosis of OS primarily relies on pathology and radiology, such as X-ray, CT, MRI, PET, and bone scans, enabling the delineation of tumor characteristics and metastatic spread, while histopathological assessment aids in subtype classification and molecular profiling (Bielack et al., 2009). However, diagnostic challenges persist, stemming from limitations in imaging sensitivity and specificity, sampling inadequacies, variability in pathological criteria, and the absence of reliable biomarkers. Treatment strategies for OS include surgical resection, chemotherapy, and radiotherapy; nevertheless, owing to its aggressive nature and therapeutic resistance, the prognosis remains dismal, with a 5-year survival rate of 60–70% (Zhu & Hou, 2020). Therefore, elucidating the pathogenesis of OS and identifying novel therapeutic targets are pivotal research endeavors and clinical imperatives.

Glycolysis is an anaerobic metabolic pathway in which glucose is broken down into lactate or ethanol, while producing a small amount of energy (Ganapathy-Kanniappan & Geschwind, 2013). Glycolysis typically only occurs under anaerobic conditions in normal cells, whereas in tumor cells, rapid glycolysis, known as the Warburg effect, occurs even in the presence of sufficient oxygen (Fernie, Carrari & Sweetlove, 2004). Increasing evidence suggests that glycolysis provides tumor cells with a rapid energy source while also generating large amounts of lactate, leading to the acidification of the tumor microenvironment (TME), inhibition of immune system responses, and promotion of tumor growth and metastasis (Li et al., 2023a). Recent studies have revealed that OS cells also exhibit the Warburg effect, relying on glycolysis to sustain rapid proliferation and invasion, even under aerobic conditions. The regulatory mechanisms and functional roles of glycolysis in OS cells remain incompletely understood; however, several studies suggested that glycolysis plays a significant role in the development of OS and represents a potential therapeutic target for this disease (Feng, Ou & Hao, 2022). Therapeutic strategies targeting glycolysis, including inhibiting the activity of glycolytic enzymes, reducing glucose uptake and utilization, and ameliorating TME acidification; these have been shown to exert several effects in vitro and in animal models; however, their clinical application requires further research and validation. Furthermore, research also revealed that glycolysis-related genes (GRGs) are closely associated with the occurrence, progression, and prognosis of OS, as evidenced by their expression levels and functional status within tumors. Therefore, GRGs can serve as crucial biomarkers for the diagnosis, prognosis prediction, and treatment of OS, suggesting their utility in personalized therapeutic strategies.

Recent studies have shown that the glycolytic activity of tumor cells can acidify the microenvironment by releasing metabolic products such as lactate, thereby inhibiting the function of immune cells and promoting immune evasion (Jiang et al., 2019). Therefore, understanding the impact of glycolysis-related genes on the immune microenvironment is crucial for understanding the progression mechanisms of osteosarcoma.

Here, patients with OS were stratified into two molecular subgroups according to GRGs, with subsequent exploration of variances in patient prognosis and TME across these subgroups. A risk signature composed of GRG was formulated based on the DEGs identified between the molecular subtypes. Furthermore, the association between the GRG risk signature and drug sensitivity, along with the distribution of immune cell types in the single-cell RNA sequencing (scRNA-seq) data, was examined relative to the four GRG risk signatures. Finally, reverse transcription quantitative PCR (RT-qPCR) was employed for the cellular-level validation of risk signatures. A simplified representation of the workflow is shown in Fig. 1.

Materials and methods

Results

Identification of glycolysis molecular subtypes in OS

The NFM algorithm was used to stratify patients with OS into subgroups based on 288 GRGs identified by univariate Cox analysis. The optimal clustering stability was achieved when K = 2 (Figs. 2A, 2B). Of the patients, 27 patients and 49 patients were clustered into Clusters 1 and 2, respectively. Heatmap visualization of the glycolytic gene expression revealed notable differences between the two clusters (Fig. 2C). Notably, patients in Cluster 2 exhibited significantly better overall survival than those in Cluster 1 (P = 0.027, Fig. 2D), indicating distinct molecular subtypes of OS characterized by glycolysis gene expression profiles. As shown in Figs. 2E, 2F, the GSEA have results shown that the genes in the Cluster 1 were mainly enriched in the signaling pathways, including “PPAR SIGNALING PATHWAY”, “NEUROACTIVE LIGAND RECEPTOR INTERACTION”, and “HYPERTROPHIC CARDIOMYOPATHY”, whereas the genes in the Cluster 2 were primarily related to “DILATED CARDIOMYOPATHY”, “GLYCEROLIPID METABOLISM”, “NEUROACTIVE LIGAND RECEPTOR INTERACTION”.

TME infiltration of molecular clusters in OS

As shown in Figs. 3A–3C, the CIBERSORT algorithm was used to explore the immune differences between the two molecular subtypes. Significant differences were observed in the abundance of various immune cell types between Clusters 1 and 2. A significant increase in the proportion of monocytes/macrophages was observed in the high glycolytic activity subtype, with a significant enrichment of immunosuppressive cells, including B cells, CD4_naive, CD4_T, CD8_naive, central_memory, cytotoxic, dendritic cells (DC), exhausted, gamma_delta, inducible regulatory T cells (iTregs), mucosal-associated invariant T cells (MAITs), monocytes, natural killer (NK) cells, NKT cells, naive T regulatory cells (nTregs), Th17 cells, and Th2 cells. Spearman’s correlation analysis was conducted to investigate the potential associations between various immune cell types and OS. Figure 3D illustrates the correlation results, revealing that the immune infiltration percentage of B cells was negatively associated with central_memory, gamma_delta, CD4_naive, exhausted, CD8_naive, Th17, Th2, and nTreg cell types, while exhibiting a positive association with neutrophils and monocytes. The immune infiltration percentage of monocytes was negatively correlated with gamma_delta, CD4_naive, central_memory, CD4_T, nTreg, Th2, Th17, CD8_naive, and exhausted cell types, whereas it was positively correlated with cytotoxic, MAITs, and B cell types. Notably, DC immune infiltration was consistently negatively correlated with exhausted, CD8_naive, Th17, Th2, and nTreg cell types. These results suggest that immune cell types in TME infiltration are a crosstalk link between the identified molecular clusters in OS.

DEGs and functional analyses

A total of 239 DEGs were detected, of which 238 were downregulated and one was upregulated in Cluster 2 compared to Cluster 1. The heatmap diagram shows the distribution of DEGs in Clusters 1 and 2 (Figs. 4A, 4B). Signaling pathway analysis suggested the screened out DEGs were closely associated with energy metabolism and glucometabolism, such as “glucose”, “6-phospofructo-2 kinase activity”, “AMP binding”, “NAD activity” (Fig. 4C). GO enrichment analysis showed that the DEGs were mainly enriched in the following biological processes (BPs): glycolysis, gluconeogenesis, and carbohydrate metabolism. The enriched cellular components (CCs) were the nuclear outer membrane, nuclear inner membrane, and extracellular exosome. Finally, the enriched molecular functions (MFs) were glucose, fructose-6-phosphate, and AMP binding (Fig. 4D). As shown in Fig. 5A, the PPI network of the DEGs was constructed using the STRING database, with a combined score of >0.4, comprising 215 nodes and 2018 interaction edges. Based on topological parameter calculations, GCK, GPI, PFKM, PGK1, INS, LDHA, G6PD, ENO1, PFKL, and GOT2 were screened out to a high degree. The DEGs in the PPI network were mainly enriched in GO terms related to the glycolysis metabolism, including “Carbohydrate metabolic process”, “Response to hypoxia”, “Glycolysis”, “Regulation of purine nucleotide metabolic process”, and “Glycosaminoglycan metabolism” (Fig. 5B). The PPI network of the modular analysis was performed using the MCODE algorithm (Fig. 4C), and the PPI network was divided into eight modules. As shown in Fig. 4C, and Table 2, the GO terms of the eight modules were mainly associated with the glycolysis metabolism, such as “Metabolism of carbohydrates”, “Regulation of Glucokinase-by-Glucokinase Regulatory Protein”, “Glycolysis/Gluconeogenesis”, “HIF-1 signaling pathway”, “Extracellular matrix organization”, “Phenylalanine metabolism”, and “Glycosaminoglycan biosynthesis”.

Development of risk signature of GRGs in OS and drug sensitivity prediction

Univariate Cox regression analysis of DEGs between molecular clusters was conducted to develop a GRGs risk signature for OS. Figures 6A, 6B and Table 3 demonstrate that four genes (CHPF, RRAGD, TPR, and VCAN) were significantly correlated with patient prognosis and identified as “risk signatures” (P < 0.01). The risk scores of all samples in the datasets using this four-gene risk signature were calculated as follows: (CHPF expression × 1.174) + (RRAGD expression × 0.852) + (TPR expression × 1.006) + (VCAN expression × 1.369). After ranking the risk scores and using the median score as the cutoff, the samples were classified into high- and low-risk groups. Survival analysis indicated lower survival rates in the high-risk group than in the low-risk group, with a p-value of 2.069 × 10⁻¹³ (Fig. 6C). Subsequently, the survival rates of patients in the two groups were compared. As shown in Fig. 6D, patients with low death and high survival rates were classified into the low-risk group. Drug sensitivity analysis using the “OncoPredict” package in R showed that patients in the low-risk group were more sensitive to AZD8055 and paclitaxel, while patients in the high-risk group were more sensitive to PD0325901, suggesting that these drugs respectively have the potential to improve outcomes in both low-risk and high-risk groups (Figs. 6E–6F). These results demonstrated that our established risk signature consisting of four GRGs could serve as a reliable predictor of OS in patients.

Table 2:

Modular analysis of PPI network.

NO.	Description	GO terms
1	MCODE 1 (Red)	Metabolism of carbohydrates
2	MCODE 2 (Blue)	Regulation of Glucokinase by Glucokinase Regulatory Protein
3	MCODE 3 (Green)	Glycolysis/Gluconeogenesis
4	MCODE 4 (Purple)	HIF-1 signaling pathway
5	MCODE 5 (Yellow)	Extracellular matrix organization
6	MCODE 6 (Orange)	Phenylalanine metabolism
7	MCODE 7 (Pink)	Histidine metabolism
8	MCODE 8 (Brown)	Glycosaminoglycan biosynthesis

DOI: 10.7717/peerj.19369/table-2

scRNA-seq analysis of characteristic GRGs

We investigated the expression patterns of the developed GRG risk signature within the TME of OS at the single-cell level using the TISCH database. Twenty-eight diverse cell subsets were identified and categorized into eight distinct subsets based on literature reviews. These subsets comprised CD4Tconv, CD8Tex, endothelial cells, fibroblasts, malignant cells, mono/macrophages, osteoblasts, and plasmocytes, as shown in Figs. 8A, 8B. Significant differences were observed in the distribution of GRGs among identified cell subsets. In Figs. 8C–8F, the expression levels of CHPF, RRAGD, and VCAN were uniformly distributed across diverse cell types, including CD4Tconv, CD8Tex, endothelial cells, fibroblasts, malignant cells, mono/macrophages, osteoblasts, and plasmocytes. Conversely, TPR expression exhibited high variability among the aforementioned immune cell types, providing evidence of substantially elevated expression of the developed GRG risk signature within specific OS cellular populations.

Table 3:

Multivariate cox regression analysis identified four GRGs signature.

Genes	Coefficient	HR	HR.95L	HR.95H	p value
CHPF	1.174	1.414	1.176	1.701	0.00023
RRAGD	0.852	1.379	1.086	1.752	0.00828
TPR	1.006	1.219	1.076	1.381	0.00180
VCAN	1.369	1.275	1.095	1.484	0.00173

DOI: 10.7717/peerj.19369/table-3

Discussion

OS, a primary malignant bone tumor that occurs globally, predominantly affects pediatric and young adult populations. Despite advancements in treatment, including surgery and comprehensive modalities, poor prognosis stems from the limitations of current therapies. However, most parameters used to construct prognostic models focus solely on the genome or transcriptome, neglecting the energy metabolism processes. Glycolysis, a critical energy metabolic pathway in cells, significantly influences the proliferation, invasion, and metastasis of tumor cells. This study established the GRGs risk signature as a crucial biomarker for prognosis and treating OS, facilitating the assessment of its biological behavior, prognosis prediction, and personalized therapeutic guidance for OS.

Initially, we applied the “NFM” algorithm to stratify OS patients into subgroups based on 288 glycolysis genes identified through univariable Cox analysis. Subsequently, we classified OS patients into two molecular subtypes based on the median gene expression in the training cohort (GSE16091 and GSE39058). Survival analysis revealed significant differences between the two molecular subtypes, with Cluster 1 demonstrating the poorest prognosis (P = 0.027). Moreover, notable differences were evident in TME infiltration, including B cells, cytotoxic, DC, exhausted, gamma_delta, iTregs, MAITs, monocytes, NK cells, and NKT cells. The glycolytic pathway serves as the primary energy source and is crucial for the growth and proliferation of OS cells. Mounting evidence suggests that glycolytic pathway activation in OS cells, coupled with suppressed tricarboxylic acid cycle (TCA) function, is indicative of the Warburg effect, thereby promoting tumorigenic activity and correlating with poor prognosis in patients with OS (Feng, Ou & Hao, 2022). Additionally, GO enrichment analysis of modules in the PPI network revealed significant terms such as “Metabolism of carbohydrates” and “Regulation of Glucokinase-by-Glucokinase Regulatory”, suggested that key metabolic genes of these modules’ were tightly associated with hypoxia.

Univariate Cox regression analysis of DEGs between the molecular clusters was conducted to develop the GRGs risk signature in OS, and four genes, including CHPF, RRAGD, TPR, and VCAN, were screened as relevant genes for the GRGs risk signature construction. As shown in the Figs. 5 and 6, the constructed GRGs risk signature had better predictive ability in the training and validation cohorts and could be utilized as an independent prognostic factor for patients with OS. CHPF functions as an enhancer protein that participates in glucose metabolism and glycolysis regulation, modulates hexokinase activity, and promotes the phosphorylation of glucose to glucose-6-phospate. Emerging evidence from both in vitro and in vivo studies supports the role of CHPF in promoting tumor growth and progression in various cancers, such as bladder, breast, and gastric cancers (Lin et al., 2021; Liao et al., 2021; Zhong et al., 2023), further underscoring its significance as a potential therapeutic target. Shen et al. (2022) further confirmed that CHPF promotes the development of OS by targeting SKP2 and activating the AKT signaling pathway, suggesting that it is a promising candidate target for treating OS. RRAGD, a member of the Rag family of small GTPases, functions as a guanine nucleotide-binding protein. Increasing evidence has shown that activation of mTORC1 by the RRAGD/RRAGA complex promotes cellular growth and proliferation by regulating protein synthesis, metabolism, and autophagy; through its role in mTORC1 signaling (Schlingmann et al., 2021), RRAGD coordinates cellular responses to nutrient availability, growth factors, and cellular stressors, thereby influencing various physiological processes, especially glucose uptake and metabolism in cancer cells (Lee et al., 2018). Lawrence et al. (2018) also found that the aberrant expression of activity of RRAGD in OS cells may contribute to the dysregulation of glycolytic metabolism, thereby promoting tumor growth and progression, suggesting that targeting RRAGD-mediated signaling pathways involved in glycolysis represents a promising therapeutic strategy for inhibiting OS progression. TPR, which encodes a tetratricopeptide repeat protein, regulates glycolytic metabolism by modulating the activity of hexokinase and phosphofructokinase, thereby influencing the rate of glycolysis and production of ATP and metabolic intermediates (Snow & Paschal, 2014). Dysregulation of TPR is associated with tumor progression and aggressiveness. TPR interacts with various oncogenic proteins, including transcription factors, and signaling molecules, to promote cell proliferation, survival, and metastasis (Kato et al., 2020). Additionally, TPR reportedly modulates cellular responses to hypoxia and nutrient deprivation. These are common features of the TME (Liu, Zhou & Tang, 2023), and aberrant expression or activity of TPR promotes glycolytic metabolism, providing cancer cells with energy and biosynthetic precursors necessary for their rapid proliferation. Furthermore, TPR may contribute to the resistance of OS cells to chemotherapy and radiotherapy, further exacerbating tumor progression and treatment challenges (Soman et al., 1991). Castillo et al. (2014) found that papillomavirus binding factor (PBF) interacts with the TPR and might in OS genesis by deregulating apoptotic mechanisms and controlling cellular transcription. VCAN has multifaceted roles in cellular processes, including modulating the expression and activity of proteins and enzymes in the glycolysis pathway; regulating glucose uptake, glycolytic flux, and the production of ATP; and influencing metabolic intermediates essential for cellular energy and biosynthesis. VCAN promotes cancer cell proliferation, survival, and invasion by modulating cell–matrix interactions, extracellular matrix remodeling, and signaling pathways involved in the tumor and metastasis of OS (Li et al., 2023b), contributing to the establishment of an immunosuppressive tumor microenvironment and enabling cells to evade immune surveillance (Zhou et al., 2023). Allen et al. (2021) developed a zebrafish model to study metastasis, isolating extravasated canine OS cells for RNA-seq analysis, and confirmed the dysregulation of KRAS signaling, immune pathways, and ECM organization, with VCAN upregulation implicated in OS metastasis. However, the specific molecular mechanisms of VCAN in relation to glycolytic metabolism and its biological role in OS still require further investigation. Our investigation further established that four-gene risk signature served as a reliable predictive indicator of prognosis in the validation cohort and by confirming cellular-level status via external dataset verification and RT-qPCR experiments.

Our study found that the increased lactate levels in the high glycolytic activity subtype may form an immunosuppressive microenvironment by inhibiting the function of CD8+ T cells and promoting the infiltration of Tregs. This suggests that targeting the glycolytic pathway may enhance the efficacy of immunotherapy. Previous studies have confirmed that lactate derived from glycolysis can induce the differentiation of monocytes into tumor-promoting M2 macrophages (Yin et al., 2024), which may be an important reason for the poor prognosis of this subtype.

The potential utility of the GRGs risk signature in predicting response to chemotherapy agents was further validated using the GDSC and CCLE databases. As shown in the Figs. 6E, 6F, patients in the low-risk group were more sensitive to AZD8055 and paclitaxel, whereas those in the high-risk group demonstrated increased sensitivity to PD0325901. AZD8055, a highly specific mTOR kinase inhibitor, exhibits promising antitumor activity by competitively binding to the ATP-binding cleft of mTOR (He et al., 2021). Preclinical studies have also indicated the potential antitumor activity of AZD8055, as it can inhibit the downstream signaling of both the mTORC1 and mTORC2 complexes, thereby regulating cellular growth, proliferation, motility, and survival in various cancer types, including bladder (Hu et al., 2023), colon (Chen et al., 2018), and hepatocellular cancers (Hu et al., 2014). Paclitaxel, a well-established antineoplastic agent, functions by binding and stabilizing microtubule proteins, disrupting tumor cell mitosis, and inhibiting cell proliferation and survival; as such, it is widely used in the treatment of various cancers (Dan, Raveendran & Baby, 2021). Evidence studies demonstrated that the paclitaxel exerts its therapeutic effects on OS through a combination of mechanisms including microtubule stabilization (Li et al., 1999), cell cycle arrest (Liu et al., 2010), apoptosis induction (Ho et al., 2021), anti-angiogenic effects (Bonzi et al., 2015), and modulation of cellular metabolism (Duan et al., 2014). PD0325901 acts as a potent suppressor of the MEK/ERK signaling pathway (El-Hoss et al., 2014), exerting anticancer effects through the inhibition of angiogenesis, enhancing the efficacy of PD-1 inhibitors, and potentially synergizing with other therapeutic agents (Luo et al., 2021). Additionally, scRNA-seq showed that OS cells gradually evolved into tumors characterized by high GRGs risk signature scores, and functional annotation of the GRGs confirmed tumorigenic pathways and inactive immunogenic pathways in OS cells with high GRGs scores. The TME is also a pivotal determinant influencing the proliferation, migration, and invasion of immune cells in OS. Although computational predictions align with preclinical evidence, future studies will validate precision treatment strategies guided by risk stratification in PDX models to promote translational applications.

Compared to previous studies, our research demonstrates several unique advantages. For instance, the study by Poudel and Koks focused on identifying differentially expressed genes in OS using FFPE and fresh tissue samples, highlighting genes related to ECM degradation and cell cycle regulation (Poudel & Koks, 2024). While this study provides a comprehensive transcriptomic analysis, our work specifically targets GRGs, which are central to cancer metabolism and have not been thoroughly explored in previous studies. Ho et al. also performed whole transcriptome analysis to identify differentially regulated networks in OS, with a focus on ECM-related pathways (Ho et al., 2017b). Our study complements these findings by emphasizing the role of glycolysis in OS progression and identifying a GRG signature with prognostic value. The analysis of repetitive DNA elements in OS by Ho et al. provides valuable insights into the genomic landscape of OS (Ho et al., 2017a). However, our study focuses on GRGs and their metabolic implications, offering a different perspective on the disease’s molecular mechanisms. Rothzerg et al. identified upregulated antisense long non-coding RNAs in OS, highlighting their potential role in tumorigenesis (Rothzerg et al., 2021). Our study, in contrast, focuses on GRGs and their prognostic significance, providing a complementary view of the molecular drivers of OS. Additionally, our study incorporates single-cell RNA sequencing (scRNA-seq) validation and drug sensitivity analysis using OncoPredict, which are not addressed in the aforementioned studies. These methods enhance the reliability and translational potential of our findings.

This study also presents a number of limitations that warrant discussion. First, the risk- score model was developed using retrospective data sourced from public databases, and additional prospective real-world evidence is required. Second, the construction of a risk score model relying solely on a single signature is unavoidable because numerous important prognostic genes in OS may have been overlooked. Third, the functional roles of TPR and VCAN in glycolytic reprogramming and immune evasion remain hypothetical, as direct experimental validation was not performed. Therefore, more diverse sample data from patients include independent datasets like TCGA or prospective cohorts at with varying risk levels are necessary to validate the robustness of this analysis. Moreover, experimental investigations are warranted to explore the relationship and mechanism between the risk score, immune response, and the predictive outcomes of drug sensitivity.

Conclusions

In conclusion, patients with OS were categorized into two molecular subtypes based on the GRGs expression matrix, with notable distinctions observed in the TME and immune cell infiltration between these subtypes. DEGs were screened out between the two molecular subtypes, and a GRGs risk signature of OS consisting of four genes (CHPF, RRAGD, TPR, and VCAN) was constructed. Furthermore, the robustness of this risk signature was validated in an independent cohort and corroborated at the cellular level using external dataset verification, scRNA-seq analysis, and RT-qPCR experiments. Collectively, these findings indicate that CHPF, RRAGD, TPR, and VCAN may be valuable for evaluating the prognosis and facilitating the clinical application of these biomarkers in the prognosis of OS. We propose a hypothesis that within the tumor microenvironment, glucose-related genes (GRGs) modulate the glycolytic process of tumor cells, thereby influencing the infiltration and function of immune cells, and consequently participating in immune evasion and drug response.

Supplemental Information