Exploring the relationship between metabolism and immune microenvironment in osteosarcoma based on metabolic pathways

Abstract

Background: Metabolic remodeling and changes in tumor immune microenvironment (TIME) in osteosarcoma are important factors affecting prognosis and treatment. However, the relationship between metabolism and TIME needs to be further explored.

Methods: RNA-Seq data and clinical information of 84 patients with osteosarcoma from the TARGET database and an independent cohort from the GEO database were included in this study. The activity of seven metabolic super-pathways and immune infiltration levels were inferred in osteosarcoma patients. Metabolism-related genes (MRGs) were identified and different metabolic clusters and MRG-related gene clusters were identified using unsupervised clustering. Then the TIME differences between the different clusters were compared. In addition, an MRGs-based risk model was constructed and the role of a key risk gene, ST3GAL4, in osteosarcoma cells was explored using molecular biological experiments.

Results: This study revealed four key metabolic pathways in osteosarcoma, with vitamin and cofactor metabolism being the most relevant to prognosis and to TIME. Two metabolic pathway-related clusters (C1 and C2) were identified, with some differences in immune activating cell infiltration between the two clusters, and C2 was more likely to respond to two chemotherapeutic agents than C1. Three MRG-related gene clusters (GC1-3) were also identified, with significant differences in prognosis among the three clusters. GC2 and GC3 had higher immune cell infiltration than GC1. GC3 is most likely to respond to immune checkpoint blockade and to three commonly used clinical drugs. A metabolism-related risk model was developed and validated. The risk model has strong prognostic predictive power and the low-risk group has a higher level of immune infiltration than the high-risk group. Knockdown of ST3GAL4 significantly inhibited proliferation, migration, invasion and glycolysis of osteosarcoma cells and inhibited the M2 polarization of macrophages.

Conclusion: The metabolism of vitamins and cofactors is an important prognostic regulator of TIME in osteosarcoma, MRG-related gene clusters can well reflect changes in osteosarcoma TIME and predict chemotherapy and immunotherapy response. The metabolism-related risk model may serve as a useful prognostic predictor. ST3GAL4 plays a critical role in the progression, glycolysis, and TIME of osteosarcoma cells.

Keywords: Metabolism; Osteosarcoma; Prognosis; ST3GAL4; Treatment response; Tumor immune microenvironment; Vitamin and cofactor.

Fig. 1

Identification of key metabolic pathways in osteosarcoma. A Kaplan–Meier curves depict the overall survival difference between pathway activity-high and pathway activity-low groups in the TARGET cohort. Red representing the pathway activity-high group and blue representing the pathway activity-low group. B Correlations between four key metabolic pathways and immune infiltration scores, only correlations that are significant are shown. C Correlations between four key metabolic pathways and abundance of 28 immune cells, only correlations that are significant are shown. D Correlations between four key metabolic pathways and expression of immune checkpoint genes. E Correlations between four key metabolic pathways and known core biological pathway scores. Correlation coefficients are calculated by Spearman’s correlation analysis, with red representing negative correlations and blue representing positive correlations. Blank represents a correlation P-value > 0.05

Fig. 2

PPI network and hub genes in the vitamin & cofactor metabolic pathway. A PPI network of vitamin & cofactor metabolic pathway genes according to the STRING database. B The top 10 hub genes of vitamin & cofactor metabolic pathway genes. C Correlations among eight matched hub genes in the TARGET cohort. Red representing negative correlations and blue representing positive correlations. Blank represents a correlation P-value > 0.05. D Univariate Cox regression analysis of overall survival for eight hub genes. E The two gene modules identified from the PPI network

Fig. 3

Metabolic pathway-related clusters and the relationship between clusters and TIME in osteosarcoma. A Consensus heatmap based on four key metabolic pathways in the TARGET cohort. B The heatmap of four key metabolic pathways between C1 and C2. C Differences of four key metabolic pathways between C1 and C2. D, E Kaplan–Meier curves depict the OS (D) and RFS (E) difference between C1 and C2. Red representing the C1 patients and blue representing the C2 patients. F Differences of immune checkpoint genes expression between C1 and C2. G The heatmap of 28 immune cells between the two clusters and the correlations of the clusters and clinical parameters. H Differences of the abundance of 28 immune cells between C1 and C2. I Differences of ImmuneScore, StromalScore and tumor purity between C1 and C2. J Differences of core biological pathway activity between C1 and C2. *P < 0.05, **P < 0.01, ****P < 0.0001

Fig. 4

Enrichment analysis of MRGs and identification of MRG-related gene clusters. A The top eight enriched terms in GO enrichment analysis of MRGs. B The KEGG pathway analysis networks of MRGs. C Consensus heatmap based on MRGs in the TARGET cohort. D The heatmap of four key metabolic pathways among GC1-3. E Differences of four key metabolic pathways among GC1-3. F, G Kaplan–Meier curves depict the OS (F) and RFS (G) difference among GC1-3. Red representing the GC1 patients, blue representing the GC2 patients, and yellow representing the GC3 patients. *P < 0.05, **P < 0.01

Fig. 5

The relationship between MRG-related gene clusters and TIME in osteosarcoma. A Differences of ImmuneScore, StromalScore and tumor purity among GC1-3. B The heatmap of 28 immune cells among the three gene clusters and the correlations of the gene clusters and clinical parameters. C Differences of the abundance of 28 immune cells among GC1-3. D Differences of immune checkpoint genes expression among GC1-3. E Differences of core biological pathway activity among GC1-3. F GSEA enrichment plot based on the HALLMARK gene set showing the relatively significantly enriched pathways in GC3 patients. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Fig. 6

Construction of a metabolism-related risk model and the relationship between the risk model and TIME. A Coefficients for the 17 genes in the risk model. B Kaplan–Meier curve depicts the OS difference between high-risk and low-risk groups in the TARGET cohort. Red representing the high-risk group and blue representing the low-risk group. C Kaplan–Meier curve depicts the OS difference between high-risk and low-risk groups in the GSE16091 cohort. Red representing the high-risk group and blue representing the low-risk group. D ROC curve analysis of the risk score for OS in the TARGET cohort. E Kaplan–Meier curve depicts the RFS difference between high-risk and low-risk groups in the TARGET cohort. Red representing the high-risk group and blue representing the low-risk group. F ROC curve analysis of the risk score for RFS in the TARGET cohort. G Differences of ImmuneScore, StromalScore and tumor purity between high and low risk groups. H Differences of the abundance of 28 immune cells between high and low risk groups. I Differences of immune checkpoint genes expression between high and low risk groups. J Differences of core biological pathway activity between high and low risk groups. K Alluvial diagram of metabolism clusters, MRG-related gene clusters, and risk levels. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Fig. 7

Immunotherapy response and drug sensitivity analysis. A Differences of TIDE-related scores between C1 and C2. B Differences of TIDE-related scores among GC1-3. C Correlations of the risk score with TIDE-related scores. D Rate of predicted clinical response to ICB immunotherapy in different clusters and risk levels. E Differences in IC50 values of cisplatin, cyclophosphamide, gemcitabine and sorafenib between C1 and C2. F Differences in IC50 values of cisplatin, cyclophosphamide, gemcitabine and sorafenib among GC1-3. G Correlations of the risk score with IC50 values of cisplatin, cyclophosphamide, gemcitabine and sorafenib. *P < 0.05, **P < 0.01, ****P < 0.0001

Fig. 8

ST3GAL4 is highly expressed in malignant cells and is closely associated with the TIME of osteosarcoma. A The dot plot shows the expression of 41 characteristic genes in 11 cell clusters. The size of the dots indicates the proportion of cells expressing a specific marker, and the color indicates the average expression level of the markers. MSC, mesenchymal stem cell; TIL, tumor-infiltrating lymphocyte. B The t-SNE plot of the 11 main cell types in osteosarcoma. C Feature plot for ST3GAL4. The color legend shows the normalized expression levels of the genes. D Violin plot showing the normalized expression levels of ST3GAL4 across the 11 cell types. E Expression of ST3GAL4 among different cell types in the GSE162454 cohort. The color indicates the average expression level of ST3GAL4. F Kaplan–Meier curves depict the OS and RFS difference between high-ST3GAL4 and low-ST3GAL4 groups in the TARGET cohort. Red representing the high-ST3GAL4 group and blue representing the low-ST3GAL4 group. G Differences of the expression of ST3GAL4 between C1 and C2. H Differences of the expression of ST3GAL4 among GC1-3. I Differences of immune checkpoint genes expression between high-ST3GAL4 and low-ST3GAL4 groups. *P < 0.05. J Differences of ImmuneScore, StromalScore and tumor purity between high-ST3GAL4 and low-ST3GAL4 groups. K Correlations of the expression of ST3GAL4 with MDSC score and TAM-M2 score. L Rate of predicted clinical response to ICB immunotherapy in high-ST3GAL4 and low-ST3GAL4 groups. M Correlations between the expression of ST3GAL4 and known core biological pathway scores. Correlation coefficients are calculated by Spearman’s correlation analysis, with red representing negative correlations and blue representing positive correlations. Blank represents a correlation P-value > 0.05

Fig. 9

Protein expression of ST3GAL4 in osteosarcoma tissues and its effects on proliferation, invasion and migration of osteosarcoma cells. A IHC staining images of ST3GAL4 in osteosarcoma tissues (#6, n = 14) and corresponding normal tissues (#2, n = 5). The IHC scores indicated that the protein expression of ST3GAL4 was higher in tumor tissues. B Kaplan–Meier curve depicts the RFS difference between high and low ST3GAL4 protein group in the Xiangya cohort. Red representing the high ST3GAL4 protein group and blue representing the low ST3GAL4 protein group. C Folded line plots showing the effect of ST3GAL4 knockdown and overexpression on the proliferation of MG-63 and U2OS cells. The blue line represents the control group and the yellow line represents the knockdown/overexpression group. D Transwell chamber experiments showing the effect of ST3GAL4 knockdown and overexpression on the invasion of MG-63 and U2OS cells. Scale bar: 100 μm. E Scratch assays showing the effect of ST3GAL4 knockdown and overexpression on the migration of MG-63 and U2OS cells. F Colony formation assays showing the effect of ST3GAL4 knockdown and overexpression on the ability of colony formation of MG-63 and U2OS cells. Data are represented as mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001

Fig. 10

ST3GAL4 regulates the glycolysis of tumor cells and the M2 polarization of macrophages in osteosarcoma. A, B Seahorse assays indicated that the knock down of ST3GAL4 inhibited glycolysis in osteosarcoma cells. Left, representative curve; Right, quantification of basal ECAR and maxi ECAR. ECAR, extracellular acidification rate. C RT-qPCR analysis is shown for PD-L1 and M2 marker CD206 in macrophages. D Flow cytometry analysis is shown for expression of CD206 in macrophages cultured with si-NC or si-ST3GAL4 tumor cells. Shown are representative plots and quantification of the percentage of CD206 positive cells in total macrophages