Yizhuo Wang1†
Jiyuan Zhu
Quan Qi- 1The Second Department of Orthopedic Surgery, First Affiliated Hospital of Harbin Medical University, Harbin, China
- 2The Fourth Department of Medical Oncology, Harbin Medical University Cancer Hospital, Harbin, China
- 3Department of pathology, First Affiliated Hospital of Harbin Medical University, Harbin, China
Background: Immunotherapy has shown considerable promise in cancer treatment, yet only a minority of osteosarcoma patients derive benefits from this approach. Hypoxia and lactate metabolism are two predominant characteristics of the tumor microenvironment. These features are crucial for molding the immune landscape and thus have the potential to act as predictive indicators for immunotherapy response.
Methods: Prognostic modeled genes were identified through univariate and multivariate Cox regression as well as LASSO regression analyses. The tumor microenvironment was evaluated using ESTIMATE, CIBERSORT, and ImmuCellAI analyses. Tide prediction and expression of immune checkpoints, MHC molecules, chemokines, interleukins, interferons, receptors, and other cytokines were utilized to estimate immunotherapy efficacy. Single-cell analysis was performed to demonstrate the expression of modeled genes among various immune cell types. Experimental validation was carried out to verify the expression and functions of SFXN4 and SQOR.
Results: A potent signature was constructed with 8 genes related to hypoxia and lactate metabolism, including MAFF, COL5A2, FAM162A, SQOR, UQCRB, SFXN4, PFKFB2 and COX6A2. A nomogram incorporating risk scores and other clinical features demonstrated excellent predictive capacity. Osteosarcoma patients with high-risk scores exhibited poor prognosis and more “cold” tumor characteristics. According to the ESTIMATE algorithm, these patients displayed lower immune, stromal, and ESTIMATE scores, partially attributed to inadequate infiltration of key immunocytes. The Ciborsort analysis similarly indicated that high-risk individuals had diminished infiltration of critical anti-tumor immune cells such as Cytotoxic T cells, CD4+ T cells, and NK cells. The low expression levels of certain immune checkpoints, MHC molecules, chemokines, interleukins, interferons, receptors, and other cytokines in high-risk cases suggested their unsatisfactory responses to immune treatment. Tide prediction further demonstrated that fewer individuals classified as high risk may exhibit sensitivity to immune checkpoint inhibitor therapy. Notably, SFXN4 was found to be highly expressed in osteosarcoma tissues and cells; it promoted the growth, migration, and invasion of osteosarcoma cells, while SQOR had the opposite effect.
Conclusion: Our research has developed a robust hypoxia- and lactate metabolism-related gene signature, providing a solid theoretical foundation for prognosis prediction, classification of “cold” and “hot” tumors, accessing immunotherapy response, and directing personalized treatment for osteosarcoma.
1 Introduction
Osteosarcoma is a highly aggressive bone sarcoma of mesenchymal origin, predominantly occurring in young individuals and often leading to mortality as a result of lung metastases (1–3). The conventional therapeutic methods routinely applied in osteosarcoma, including neoadjuvant chemotherapy, surgical resection, and postoperative adjuvant chemotherapy, have remained largely unchanged for many years (4–7). Compared to patients with local diseases whose survival rate at five years reached 65-70% following these treatments, those with metastatic disease had a significantly lower survival rate of only 20% (8–10). Hence, it is urgent to identify reliable tools for accurate initial diagnosis, introduce novel therapeutic strategies for osteosarcoma, and establish methods for ongoing monitoring of treatment efficacy (11).
Hypoxia signifies a reduction in the tissue partial pressure of oxygen [ptiO (2)], which is a critical characteristic of the tumor microenvironment (TME) commonly observed in multiple solid tumors (12). Any factor triggering an imbalance between the high proliferative capacity of cancer cells and oxygen availability, such as overlong diffusion distances and structural abnormalities in vascular vessels, can contribute to hypoxia (13). Notably, tumors can adapt to hypoxic conditions to promote growth, aggressiveness, angiogenesis, metastatic potential, and resistance to treatment through transcription factors, especially hypoxia-inducible factor 1α (HIF1α) (14). Given the impact of hypoxia on the malignant biological function of osteosarcoma, identifying key signature genes within the hypoxia gene set could facilitate the development of therapies targeting hypoxia-induced tumorigenesis in osteosarcoma.
Lactate is more than a mere metabolite and also acts as a signaling molecule whose intracellular and extracellular concentrations rise because of the essential effect of aerobic glycolysis in energy production (15, 16). Hypoxia invariably triggers lactate production, which subsequently accumulates in the TME, in turn facilitating hypoxic response (17). In addition, Zhang et al. has identified lactylation as a significant epigenetic modification modulating the transcription of numerous oncogenes and tumor suppressor genes, thereby influencing tumor progression (18). Several studies have highlighted the impact of lactic acid on tumor cell progression, metastasis, epithelial-mesenchymal transition, angiogenesis, and immune evasion (19–21). Considering the multifaceted roles of lactate, a comprehensive understanding of the mechanisms by which lactate metabolism-related genes operate in osteosarcoma will supply valuable insights into tumor reprogramming and therapeutic strategies.
Metabolic disorders in cancers are marked by hypoxia and subsequent elevated levels of lactate (22). This interplay contributes to the creation of a restrictive TME dampening the efficacy of immune effector cells, and ultimately weakens the immune system’s ability to combat tumors (23). Hypoxia plays a vital role in promoting tumor immune evasion by upregulating certain immunosuppressive immune checkpoints, cytokines, growth factors, and interleukins to accelerate tumor progression (24). Furthermore, hypoxia hinders immune response by reducing infiltration and function of CD8+ T cells, impairing maturation and activities of dendritic cells coupled with natural killer cells, evoking M2-like polarization of tumor-associated macrophages, and boosting regulatory T cells and myeloid-derived suppressor cells activities (25, 26). Simultaneously, lactate augments antitumor immunity by driving M2 polarization, enhancing stemness of CD8+ T cells, and promoting PD-1 expression in regulatory T cells (27–29).
In this study, our objective was to comprehensively elucidate the heterogeneity of the TME and precisely predict the effects of immunotherapy in osteosarcoma. In terms of the impact of hypoxia and lactate metabolism on tumor immunity, we identified 8 genes related to hypoxia and lactate metabolism as potential predictive markers and developed a hypoxia- and lactate metabolism-related gene signature that held promise for clinical application. Osteosarcoma patients stratified based on this prognostic model exhibited distinct survival outcomes, TME characteristics, immune cell infiltration, expression levels of immune checkpoints, immune response, and sensitivity to drugs.
2 Materials and methods
2.1 Data source and preprocessing
Transcriptional expression data and corresponding clinical information of osteosarcoma samples were sourced from the TARGET (https://ocg.cancer.gov/programs/target) database and GEO (http://www.ncbi.nlm.nih.gov/geo/) database. Transcriptome data of para-cancerous tissues were retrieved from the GTEx (https://xenabrowser.net/datapages/) database. Hypoxia-related genes (HRGs) and lactate metabolism-related genes (LMRGs) were obtained from The Molecular Signatures Database (MSigDB, https://www.gsea.msigdb.org/gsea/msigdb/index.jsp). By exploring the keywords “hypoxia” and “lactate”, we found 5 priority LMRG datasets and 7 major HRG datasets in the MSigDB database. After filtering out duplicate genes in the identified LMRG and HRG datasets, a total of 284 LMRGs and 493 HRGs were selected for the forthcoming study (Supplementary Table S1).
2.2 Identification of DEGs, DELMRGs and DEHRGs
The R package “limma” was utilized to screen differentially expressed genes (DEGs) meeting the criteria of absolute fold change greater than 2 and p values less than 0.05. Differentially expressed hypoxia-related genes (DEHRGs) were determined by intersecting DEGs with HRGs, while differentially expressed lactate metabolism-related genes (DELMRGs) were recognized by intersecting DEGs with LMRGs.
2.3 Development of a hypoxia- and lactate metabolism-related gene signature
Univariate Cox regression analysis was conducted based on the R package “survival” to pick out meaningful candidate genes associated with the overall survival time of patients among DELMRGs and DEHRGs. Genes were deemed statistically significant if their p values were less than 0.05. To mitigate overfitting and refine gene selection, the least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed. Finally, a total of 8 genes were chosen for inclusion in the gene signature based on multivariate Cox regression analysis.
2.4 Derivation of the prognosis classifier
The following equation was used to determine the risk score: risk score = (0.6049967×UQCRB) + (-0.4744143×SQOR) + (0.4064687×COL5A2) + (0.5114636×FAM162A) + (0.8507392×SFXN4) + (-1.187998×PFKFB2) + (0.4085068×MAFF) + (0.2279495×COX6A2). Subsequently, the risk score of each osteosarcoma sample in the three datasets were calculated. In accordance with the median risk score of each database, the osteosarcoma samples were then stratified into two groups: low-risk and high-risk.
2.5 Survival analysis and measurement of predictive ability
The R package “survminer” was utilized for conducting Kaplan-Meier survival analysis to compare survival between the two risk groups using the optimal cutoff value of gene expression. The prediction performance of this model was appraised by employing the “timeROC” package to generate receiver operating characteristic (ROC) curves and compute the corresponding area under the curve (AUC) values.
2.6 Construction of the prognostic nomogram
A nomogram was organized to predict the 1-, 3-, and 5-year survival rates of osteosarcoma patients. The predictive accuracy of the nomogram was rated using ROC curves, C-index, and calibration curves. The net benefit of the nomogram and other clinical characteristics was assessed by Decision curve analysis (DCA).
2.7 Functional annotation enrichment
The R package “clusterProfiler” was utilized for conducting GO functional enrichment analysis, while the DAVID (https://david.ncifcrf.gov/summary.jsp) online database was searched for KEGG functional enrichment analysis. The top 10 results from the GO functional enrichment analysis in terms of 3 distinct aspects were selected based on smaller p values. Similarly, the top 10 results from the KEGG analysis were also chosen. Bar charts and bubble charts were visualized using the “ggplot2” package.
2.8 Evaluation of tumor microenvironment pathways
The 29 pathways related to the tumor microenvironment (TME) were initially identified in a previous study (30). Subsequently, the R package “clusterProfiler” was used to calculate single-sample gene set enrichment analysis (ssGSEA) scores to quantify the pathway activity. Then, we identified statistically significant associations among the various TME pathways. The “ggcor” package was utilized to unmask the correlation between risk scores and the TME.
2.9 Immune status analyses
The immune cell abundance identifier (ImmuCellAI) was employed to evaluate the infiltration levels of 24 immune cell types in the TARGET database and GSE21257 dataset. Additionally, the CIBERSORT algorithm was conducted to appraise the infiltration abundance of 22 immune cell types. The correlation analysis was for investigating the relationship between the expression of the 8 HLMRGs and cell infiltration inferred from ImmuCellAI and CIBERSORT algorithm. Moreover, four key indicators of the tumor microenvironment components were evaluated through the “estimate” R package.
2.10 Immunotherapy and molecular therapy response prediction
Chemokines, interleukins, interferons, receptors, and other cytokines were sourced from a published article (31). The TIDE online tool was employed to forecast potential clinical response to immunotherapy in osteosarcoma patients. A higher TIDE prediction score typically indicates increased immune evasion by the tumor and reduced likelihood of positive response to immune checkpoint inhibitors (ICIs). Genes associated with T-cell-inflamed gene expression profile (GEP), Th1/IFNγ gene signature, and cytolytic activity were obtained from a prior study and their scores were computed using the “ssGSEA” algorithm to predict response to immune checkpoint therapy (30).
2.11 Drug sensitivity analysis
The diversity in medication sensitivity was investigated by conducting a comparison of the IC50 values of more than 100 chemotherapeutic agents through the R package “pRRophetic”. Setting a significance threshold of p smaller than 0.005, we identified several efficacious chemotherapy agents in both risk groups.
2.12 Single-cell analysis
The scRNA sequence matrix was imported using the R package “Seurat.” To acquire high-quality scRNA-seq data, specific filtering conditions were applied for every cell: only cells with over 500 UMI counts were included; cells expressing fewer than 300 genes and equal to or more than 5000 genes were excluded; cells with over 10% mitochondrial gene expression were also eliminated. Normalization of the scRNA-seq data was carried out by applying the “NormalizeData” function, with the normalization method set to “LogNormalize”. The top 3,000 highly variable genes were recognized through the “FindVariableFeatures” function. Subsequently, principal component analysis (PCA) was carried out to lessen the dimensionality of the raw data based on the top 3,000 genes using the “RunPCA” function. Principal components with statistical significance were discerned using “JackStraw” function, and the top 14 meaningful principal components were selected for cell clustering based on variance explained. Cell clustering analysis was conducted using the “FindNeighbors” and “FindClusters” functions, followed by uniform manifold approximation and projection (UMAP) analysis through the “RunUMAP” function. Batch correction was implemented by applying the “Harmony” package. The “FindAllMarkers” function was carried out to identify marker genes in each cluster. The cell annotation was from a previously published paper (32).
2.13 Cell culture
Human osteosarcoma cell lines (HOS, 143B) and human normal osteoblast cell lines (hFOB1.19) were purchased from Procell Life Science and Technology Co., Ltd. (Wuhan, China). Cells were cultured with 10% fetal bovine serum (PAN, German) and 1% penicillin-streptomycin solution (SEVEN, China) in a humidified environment at 5% CO2 and 37°C.
2.14 Transfections
Tumor cells were cultured in a six-well plate at a density of 3×105 cells per well in appropriate culture medium and incubated overnight with agitation. Upon reaching a cell density of 60 to 70%, the cells were washed with PBS and supplemented with serum-free culture medium, followed by the addition of the plasmid-lipo2000 mixture to each well. The six-well plates were then incubated in the incubator for 6 hours before being switched to serum-containing culture medium. The sequence of the siRNA purchased from General Biol (China) is as follows:GAUCAAAGCUAGAGUGACUTT (si-SQOR-1); ATCTTTACCTTCCCAAATACTCC (si-SQOR-2); GGACAAAGGCUGUUAGAGATT (si-SFXN4-1); CUGACUGGCCCUUGGAUUATT (si-SFXN4-2).
For virus transfection, 143B cells were seeded in 24-well plates at a density of 1 × 105 cells per well and incubated for 18 to 24 hours. The original medium was replaced with 2 ml of fresh medium containing 6 µg/ml polybrene, along with an appropriate quantity of virus (Genechem, Shanghai, China). After four hours, an additional 2 ml of fresh medium was added to dilute the polybrene. Forty-eight hours post-transfection, the cells were screened using new medium supplemented with puromycin. Following another forty-eight hours, the medium was again replaced with fresh medium before harvesting the cells.
2.15 CCK8 assay
100 μL of suspension containing 5000 transfected cells was dispensed into each well of a 96-well plate along with 10 μL of CCK8 solution (Beyotime, China). The plate was then placed in a cell culture incubator for 1 hour, following which absorbance readings were taken at 450 nm and recorded.
2.16 Cell colony assay
3600 osteosarcoma cells were equally distributed into six-well plates and incubated at 37°C with 5% CO2 for 14 days with regular medium changes. Post-incubation, the cells were fixed and stained using 4% paraformaldehyde and 0.1% crystal violet for 20 minutes each, after which images were captured and data documented.
2.17 Scratch assay
Upon the cell density reaching 80 to 90 percent on a six-well plate, a scratch was created using a pipette tip perpendicular to the plate surface. Following a wash with PBS to remove unattached cells, microscopic images were captured at 0 and 48 hours.
2.18 Transwell assay
The transfected HOS and 143B cell lines were cultured with the serum-free DMEM and serum-free 1640, respectively, in a Transwell upper chamber. Corresponding culture medium containing 10% FBS was added to the lower chamber. The cells were incubated at 37°C with 5% CO2 for 48 hours, fixed with formaldehyde, stained with crystal violet, and visualized under a microscope for analysis.
2.19 qRT-PCR
The total RNA was extracted from HOS and 143B as well as hFOB1.19 using TRIzol reagent (Invitrogen, USA). Following the manufacturer’s instructions, the PrimeScriptTM RT reagent kit (Takara, Japan) was used to create the cDNA by reverse transcription. The relative expression of SQOR and SFXN4 was measured according to the 2-ΔΔCt algorithm by normalizing the samples to GAPDH utilizing a LightCycler 480 Fluorescence Quantification System (Roche, Basel, Switzerland).
2.20 Western blotting
The Radioimmunoprecipitation (RIPA) lysis buffer (Beyotime, China) was used to isolate the proteins and then transferred the proteins to onto PVDF membranes (Millipore, Billerica, MA, USA) by electrophoretic separation. The membranes were closed for one hour at room temperature with 5% nonfat milk in TBST before being incubated at 4° with a particular primary antibody. After overnight incubation the bands were visualized using the ECL kit. Specific antibodies included: SFXN4 (YT4298, Immunoway Biotechnology Company, China), SQOR (21112-1-AP, Proteintech Group Inc., Wuhan, China).
2.21 Tissue specimens and immunohistochemistry
The expression of two genes (SQOR and SFXN4) significantly affecting prognosis was confirmed through immunohistochemical experiments. Four osteosarcoma patients from the First Affiliated Hospital of Harbin Medical University voluntarily donated tumor tissues and normal bone tissues. All of them signed informed consent forms before tissue donation. Prior to surgery, none of them received radiation therapy or chemotherapy. The hospital ethics committee gave approval for this research. Slices of tissues were fixed, dehydrated, and paraffin-embedded. Then, the tissues were dewaxed and hydrated through rinses with xylene, anhydrous ethanol, ethanol, and distilled water after baking at 60°C. For 12 minutes, tissues were submerged in 0.3% H2O2 to deactivate endogenous peroxidase. Sections were sealed with serum and then treated with the appropriate antibodies for an overnight period at 4°C. Next, all sections were exposed to secondary antibodies for 30 minutes at 37°C. Two pathologists employed 3,3-diaminobenzidine (DAB) reagent (Boster, Wuhan, China) to stain every tissue piece. The antibodies utilized in the aforementioned SQOR (1:200) and SFXN4 (1:200) experiments were both from. Microscopic images were captured and analyzed using cellSns software.
2.22 In vivo experiment
The BALB/c nude mice were procured from Huachuang Sino (Jiangsu, China) and housed at the Animal Experiment Center of the First Hospital of Harbin Medical University. A total of ten four-week-old nude mice were randomly assigned to two groups and subcutaneously injected bilaterally with 1×106 cells each of 143B wild-type, SFXN4 knockdown, and SQOR knockdown cells. Tumor volume was monitored and measured continuously for a period of 15 days post-implantation. Following this observation period, the mice were euthanized. All procedures involving the nude mice adhered strictly to the regulations set forth by the Ethics Committee of the First Hospital of Harbin Medical University (Batch Number: IRB-AF/SC-0402.0).
2.23 Statistical analysis
The Wilcoxon signed-rank test was applied to recognize differences between two groups. Spearman’s rank correlation coefficient was employed to assess the correlation between risk scores and pathway enrichment scores calculated by the “ssGSEA” algorithm. A p value less than 0.05 was considered as statistically significant. * denotes p < 0.05, ** represents p < 0.01, *** means p < 0.001 and **** indicates p < 0.0001.
3 Results
3.1 Construction of the prognostic HLMRGS and validation of its predictive ability
The research process was illustrated in Figure 1. The differential analysis of transcriptome data between osteosarcoma patients in the TARGET database and normal specimens in the GTEx database revealed 2925 upregulated genes and 3361 downregulated genes (Figure 2A). A Venn diagram depicted 234 differentially expressed hypoxia-related genes (DEHRGs) and 129 differentially expressed lactate metabolism-related genes (DELMRGs), with 6 shared genes between them (Figure 2B). Univariate Cox regression analysis identified 51 genes out of the DEHRGs and DELMRGs that were associated with prognosis, including 20 DEHRGs and 31 DELMRGs (Supplementary Table S2). Subsequent LASSO Cox regression analysis was carried out to reduce overfitting and filtered 8 genes (Figures 2C, D). Multivariate Cox regression analysis confirmed that these 8 identified genes, including 2 protective elements (SQOR and PFKFB2) and 6 hazardous factors (MAFF, COL5A2, FAM162A, UQCRB, SFXN4, and COX6A2), could independently predict the overall survival of osteosarcoma samples (Figure 2E). Eventually, an optimal prognostic gene signature named hypoxia- and lactate metabolism-related gene signature (HLMRGS) was generated with those 8 determined genes. The correlation of 8 genes was portrayed as follows (Figure 2F). Among these genes, MAFF, COL5A2, FAM162A and SQOR were DEHRGs, while UQCRB, SFXN4, PFKFB2 and COX6A2 were DELMRGs. According to the Kaplan-Meier (KM) survival analysis of 8 modeled genes, elevated expression levels of COL5A2, COX6A2, UQCRB, SFXN4, FAM162A and MAFF sharply shortened the survival time of osteosarcoma patients (Supplementary Figures S1A–H). Conversely, high expression levels of PFKFB2 and SQOR were associated with prolonged survival. The risk scores derived from the expression levels of these 8 modeled genes effectively categorized samples into low- and high-risk groups. The findings indicated that high-risk patients experienced worse clinical outcomes (Figures 2G, H). Receiver operating characteristic (ROC) analysis reflecting the sensitivity and specificity of gene signatures showed that the area under curve (AUC) values of the HLMRGS were 0.916 at 1 year, 0.932 at 3 years and 0.950 at 5 years (Figure 2I). After the establishment of HLMRGS within the TARGET database, its prognostic significance and predictive capacity were verified externally in the GSE39055 and GSE21257 datasets, thereby affirming both its repeatability and reliability (Supplementary Figures S2A–F). Altogether, the collective findings unequivocally validated the prognostic value of the HLMRGS and its predictive potential.
Figure 1. Flow chart of this study.
Figure 2. Development of the HLMRGS and verification of its predictive ability in the TARGET database. (A) Volcano chart of 6286 DEGs between osteosarcoma and normal samples. (B) Venn diagram depicting the junction of HRGs, LRGs, and DEGs. (C, D) LASSO penalized Cox regression, and (E) multivariate Cox regression analyses of the DEHRGs and DELMRGs. (F) The association network diagram and types of 8 HLMRGs. (G) The risk scores, clinical outcomes, and expression of 8 modeled genes of osteosarcoma samples divided by median risk score. (H) KM survival analysis for the two risk cohorts. (I) The ROC curve of HLMRGS at 1, 3, 5 years.
3.2 Stratified analysis by different clinical characteristics of osteosarcoma patients
To further elucidate the universal applicability of the HLMRGS, we conducted survival outcome predictions for osteosarcoma patients stratified by various clinical characteristics using the HLMRGS. Osteosarcoma samples from the TARGET database were grouped into diverse subcategories bottomed on age (≥18 vs. <18 years), gender (male vs. female), metastatic stage (metastatic vs. nonmetastatic), primary tumor location (leg vs. arm/pelvis), and specific tumor site (lower limb vs. upper limb vs. pelvis/sacrum/ilium) (Supplementary Table S3). KM survival analysis was performed across these subgroups, illustrating that individuals with high-risk scores exhibited significantly poorer outcomes amidst most subgroups, including age, gender, and metastatic stage (Supplementary Figures S3A–K).
3.3 Clinical value of the HLMRGS in osteosarcoma
After conducting univariate and multivariate Cox regression analyses, a nomogram was developed to validate the clinical significance of the HLMRGS. The univariate Cox regression analysis indicated that metastasis, primary tumor site, specific tumor site, and risk scores were relevant to overall survival time of osteosarcoma specimens from the TARGET database (Figure 3A). In addition, the multivariate Cox proportional hazard models revealed that metastasis, specific tumor site, and risk scores were independent prognostic factors for osteosarcoma patients (Figure 3B). Then, a nomogram model was structured based on these three independent prognostic indicators to predict the 1-, 3-, and 5-year survival probabilities of osteosarcoma samples (Figure 3C). The calibration plots depicted a satisfactory agreement between predicted and actual outcomes with a high C-index of 0.915 (Figure 3D). The AUC values of the nomogram were a whopping 0.916 at 1 year, 0.932 at 3 years and 0.950 at 5 years (Figure 3E). The decision curve analysis (DCA) curve suggested that the nomogram provided the highest net benefit, strongly confirming its clinical suitability in osteosarcoma patients (Supplementary Figure S3L). In particular, the HLMRGS displayed superior accuracy in comparison to previously published gene signatures related with hypoxia or lactate metabolism, as evidenced by its highest AUC value and C-index (Figure 3F). Overall, these results highlighted the applicability of the HLMRGS in clinical practice for osteosarcoma patients.
Figure 3. Construction of the clinical nomogram and validation of its predictive power in the TARGET database. (A, B) The forest diagram of univariate and multivariate Cox regression analyses comprising risk scores and distinctive clinical features. (C) The nomogram based on three factors with statistical significance, involving metastasis, specific tumor site and risk scores. Total scores were summed up by points of each variable. The survival probabilities in three different time periods are indicated on the axis at the bottom. (D) Calibration curves showing a tremendous fit between the predictive outcomes and actual outcomes. (E) ROC curves of the nomogram depicting the prognostic performance of HLMRGS. (F) Comparison of AUC values and C-index between the HLMRGS and other gene signatures related with hypoxia or lactate metabolism.
3.4 Functional enrichment analysis based on the HLMRGS
To explore the functional enrichment disparities between the two cohorts delineated by the HLMRGS, we conducted GO and KEGG enrichment analyses based on the differentially expressed genes (DEGs) between the low- and high-risk groups in the TARGET database. A total of 84 DEGs were discerned between these two groups, with 47 genes highly expressed and 37 genes downregulated (Supplementary Figure S4A). KEGG pathway analysis revealed evident enrichment of DEGs in pyrimidine metabolism, drug metabolism, EGFR tyrosine kinase inhibitor resistance and TGF-β signaling pathways (Supplementary Figure S4B). Additionally, GO pathway analysis exhibited predominant enrichment in the regulation of cellular response to growth factor stimulus, immunoglobulin complex, and signaling receptor activator activity (Supplementary Figure S4C).
3.5 Heterogeneity of the TME between the two risk groups
To unravel the heterogeneity within the TME, we did single-sample gene set enrichment analysis (ssGSEA) to estimate the activity of 29 TME pathways encompassing immune, matrix, and biological characteristics. Osteosarcoma patients were stratified into two cohorts according to their risk scores and distinct clinical features to compare the activity of 29 TME pathways. The two subgroups separated by only five factors namely, sex, age, metastasis, tumor site, and risk scores, demonstrated differences in TME pathway activities (Figure 4A and Supplementary Figure S5A). Notably, the risk scores-based group showed the most significant disparities in the activity of a total of 17 TME pathways compared to other groups. In addition, only 11 TME pathways linked with immunity exhibited obvious differences when categorized by risk scores (Figure 4B and Supplementary Figure S5B). These pathways reflected elevated activity levels in the low-risk group, as observed in both the TARGET database and the GSE21257 dataset. The pathways identified included MHC-II, coactivation molecules, checkpoint molecules, effector cells, NK cells, T cells, cancer-associated fibroblasts, myeloid cell traffic, macrophage and DC traffic, effector cell traffic and Th2 signature. More importantly, the risk scores displayed the strongest correlation with all 29 TME pathways, with a high degree of interrelation among the majority of TME gene signatures (Figure 4C and Supplementary Figure S5C).
Figure 4. Heterogeneity of the tumor microenvironment between the two risk groups in the TARGET database. (A) The differentially expressed activity of 29 TME pathways between the two subgroups using “limma” package. (B) Heatmap and box line plot depicting the activity discrepancies of 29 TME pathways between the two cohorts. (C) Multidimensional correlation plot demonstrating the relationship of 29 TME pathways, and risk score as well as other clinical features. ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001.
Since the risk scores were primarily associated with immunity, we aimed to quantify the infiltration of diverse immune cell types to further elucidate the heterogeneity of the TME by using three different methodologies: ImmuCellAI, Cibersort and ESTIMATE (Figure 5A and Supplementary Figure S6A). Among the 8 key genes comprising the HLMRG, the protective factors SQOR and PFKFB2 exhibited positive correlations with the abundance of most infiltrating immune cells, while the risk factors FAM162A and COL5A2 showed negative associations with immune cell infiltration in both the TARGET and GSE21257 cohorts (Figure 5A and Supplementary Figure S6B). In accordance with the ImmuCellAI database, low-risk patients had a noticeably enriched population of some immune-activated cell types in both the TARGET and GSE21257 cohorts, such as CD4+ T cells, CD8+ T cells, cytotoxic T cells, dendritic cells and natural killer cells, manifesting the latent capacity of humoral immunity and cellular immunity in low-risk patients (Figure 5A and Supplementary Figure S6C). Regrettably, the distribution differences of most immune cells between the two risk groups obtained by the CIBERSORT algorithm were not statistically significant in both the TARGET and GSE21257 cohorts (Figure 5A and Supplementary Figure S6D). ESTIMATE analysis revealed that the low-risk group displayed increased stromal, immune, and ESTIMATE scores in the TARGET database (Figure 5B). The results were consistent in the GSE21257 cohort, with the exception of stromal scores, where statistical significance was not observed (Supplementary Figure S6E). Furthermore, the risk scores were inversely correlated with stromal, immune, and ESTIMATE scores in both the TARGET and GSE21257 cohorts (Figure 5C and Supplementary Figure S6F). In brief, the osteosarcoma patients with low-risk scores exhibited more “hot” tumor features. These results suggested that the two groups classified by the HLMRGS displayed quite different infiltration of certain key immune cells, which may be one of the explanations for the distinction in survival outcomes between them.
Figure 5. Divergences in the immune landscape between the low- and high-risk groups in the TARGET dataset. (A) Heatmap and box line plot showing the variances in immune infiltration with the ImmuCellAI and CIBERSORT algorithms. Bubble diagram depicting the correlation between the modeled genes expression and ImmuCellAI, and CIBERSORT scores. (B) Comparison of stromal, immune, and ESTIMATE scores between the two groups. (C) Correlation analysis of risk scores with stromal, immune, and ESTIMATE scores. ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001.
3.6 Prediction of response to immunotherapy and chemotherapy
We further investigated the differential expression of immune checkpoints, MHC molecules, and immunomodulators in low- and high-risk groups to predict the efficacy of immunotherapy in osteosarcoma samples. As essential prognostic markers of immunotherapy response, some immune checkpoints were considerably upregulated in the low-risk patients in the TARGET and GSE21257 cohorts, including BTN3A1, LAG3, CD209, TNFSF4, BTN2A2, CD70, TNFSF14, HAVCR2, CD96, CD200R1, LGALS9, TIGIT, PDCD1LG2, and SIRPA (Figure 6A and Supplementary Figure S7A). Additionally, a remarkably negative correlation was observed between the risk scores and the expression levels of these immune checkpoints (Figure 6B). Furthermore, a majority of MHC molecules, such as HLA-DQA1, HLA-DOA, HLA-DMA, HLA-DRA, HLA-E, HLA-DPB1, HLA-DPA1, and HLA-DMB, were expressed at considerably high levels in the low-risk group (Figure 6C and Supplementary Figure S7B). Meanwhile, the levels of most chemokines, interleukins, interferons, receptors, and other cytokines were also increased in the low-risk group, implying the enhanced sensitivity to immunomodulators and better response to immunotherapy compared to those of high-risk group (Figure 6D). The IFNG scores, Merck18 scores and T-cell dysfunction scores were found to be higher in the low-risk group of the TARGET and GSE21257 cohorts, while MDSC scores were lower in this group according to data from the TIDE website (Figure 7A). Remarkably, TAM M2 scores were increased in the high-risk group of the GSE21257 cohort (Supplementary Figure S7C). Subsequent analysis of immunotherapy response from the TIDE website exhibited a higher percentage of responders in the low-risk group (54.8% vs. 47.6%), suggesting that high-risk patients may have a tendency to escape immune surveillance and exhibit a less favorable response to immune checkpoint inhibitors (ICIs) (Figure 7B). Additionally, the T-cell-inflamed gene expression profile (GEP) score, Th1/IFN gene signature score, cytolytic activity score, and tumor infiltrating lymphocyte score were significantly elevated in the low-risk samples (Figures 7C–F and Supplementary Figures S7D–G) and displayed a negative correlation with risk scores (Figures 7G–J). The assessment of chemotherapy response through the computation of half-maximal inhibitory concentration (IC50) values for more than 100 chemotherapeutic agents indicated that high-risk patients were more likely to benefit from pictilisib, BAY.61.3606, doramapimod, CGP.082996, Elesclomol, and GW.441756, while the low-risk group was predicted to exhibit greater sensitivity to ponatinib, serdemetan, TCS JNK 6o, MG.132, dactolisib, and WH.4.023 (Figures 8A–L). These findings disclosed that the efficacy of immunotherapy and chemotherapy differed in the HLMRGS-based groups, and low-risk patients tended to have improved immunotherapy outcomes.
Figure 6. Evaluation of immunotherapy response of the low- and high-risk groups in the TARGET database. (A) Box plot showing the expression variations in immune checkpoints between the two cohorts. (B) Multidimensional correlation plot illustrating the relationship of risk scores and the gene expression of immune checkpoints reaching statistically significant differences in both the TARGET and GSE21257 cohorts. (C) Expression divergences of MHC molecules of the two subgroups. (D) Heatmap illustrating the differences in the gene expression of chemokines, interleukins, interferons, receptors, and other cytokines between the two risk groups. ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001, ∗∗∗∗ p < 0.0001.
Figure 7. Differences in TIDE prediction and pathway activities of immunotherapy-related pathways between the low- and high-risk groups in the TARGET database. (A) TIDE prediction of the immunotherapy effectiveness of the two cohorts. (B) The bar graph depicting the proportion of patients sensitive and resistant to immunotherapy. (C–F) Violin plots presenting the distinctions of cytolytic activity score, tumor infiltrating lymphocyte score, Th1/IFN gene signature score, and T-cell-inflamed gene expression profile (GEP) score. (G–J) Correlation analysis between the risk score and the four scores mentioned above. ∗ p < 0.05, ∗∗ p < 0.01, ∗∗∗ p < 0.001.
Figure 8. Discrepancies in treatment response to chemotherapy drugs of osteosarcoma patients between the low- and high-risk groups in the TARGET dataset. (A–L) Box plots showing the IC50 values for clinical drugs with statistical significance amidst the two risk groups. ** p < 0.01, *** p < 0.001, **** p < 0.0001.
3.7 Single-cell analysis of the HLMRGs in TME-associated cells
Considering the tight association of the risk scores with TME, we examined the expression of 8 HLMRGs in various TME-associated cells at the single-cell level using the GSE162454 dataset. Following the implementation of initial quality control procedures, a gene expression profile was generated for a total of 30,660 cells encompassing 33,538 genes obtained from the primary tumors of 6 osteosarcoma samples (Figures 9A, B). The Harmony package was employed to mitigate batch effects (Figure 9C). Through unbiased clustering of the cells, 19 main clusters were identified in parallel, based on uniform manifold approximation and projection (UMAP) analyses of their gene profiles (Figure 9D). Subsequent cell annotation was performed using markers specific to different cell types as described in a previously published study, identifying 11 distinct cell populations denoted as B cells, chondroblastic cells, endothelial cells, myeloid cells, NK cells, osteoblastic cells, malignant cells, T cells, pericytes, proliferation cells, and fibroblast cells (Figure 9E). In general, UQCRB was found to be the most highly expressed gene, while COX6A2 exhibited the lowest expression across all TME-related cell subtypes (Figure 9F). Remarkably, SQOR displayed significant variation in gene expression between malignant cells and osteoclasts (Figures 9F, G). The tumor suppressor genes PFKFB2 showed quite low expression in malignant cells, which was consistent with its role as a protective factor in osteosarcoma as predicted by bioinformatics.
Figure 9. Single-cell analysis of 8 model genes in the GSE162454 dataset. (A–C) Quality control and de-batching effect of the single-cell sequencing raw data. (D) The UMAP plot of 30,660 cells colored by the 19 clusters. (E) The UMAP plot of 30,660 cells colored by 11 different cell subtypes. (F, G) The expression of 8 modeled genes across various cell subtypes.
3.8 Experimental validation of the expression and biological functions of SFXN4 and SQOR
We selected two genes that have not been reported in osteosarcoma for further experimental validation. These genes were chosen based on their classification as protective and hazardous factors, respectively (Figure 2E). SQOR was picked out from two protective elements in the multivariate Cox regression on account of its visibly low p value. Additionally, the decision to focus on SFXN4 was supported by its obviously high HR value compared to other modeled genes. The RT-qPCR, western blot and immunohistochemistry (IHC) analyses confirmed that SFXN4 displayed elevated expression at both RNA and protein levels in osteosarcoma tissues and cells, whereas SQOR exhibited heightened abundance in normal tissues and cells (Figures 10A–C). Transfection of short interfering RNA fragments into HOS and 143B cell lines resulted in decreased expression levels of SFXN4 and SQOR (Figures 10D, E). Knockdown of SFXN4 led to decreased cell proliferation, migration, and invasion, while knockdown of SQOR had the opposite effect (Figures 10F, G, 11A–D). Then, we performed in vivo experiments to further investigate the functions of SFXN4 and SQOR in osteosarcoma. As shown in Figures 12A–C, knockdown of SQOR enhanced tumor growth, whereas knockdown of SFXN4 yielded an opposing effect (Figures 12D–F).
Figure 10. Expression and biological functions on growth of SFXN4 and SQOR in osteosarcoma. (A) SFXN4 and SQOR expression in para-cancerous and osteosarcoma tissues by IHC. (B, C) The expression of SFXN4 and SQOR at RNA and protein levels in normal osteoblast cell and osteosarcoma cell lines through qRT-PCR and western blotting assays. (D, E) Effect of knocking down SFXN4 and SQOR at RNA and protein levels in osteosarcoma cell lines by qRT-PCR and western blotting assays. (F, G) CCK8 and colony formation assays demonstrating the proliferation capacity ability after knocking down SFXN4 and SQOR in osteosarcoma cell lines. ∗ p < 0.05, ∗∗ p < 0.01.
Figure 11. Effect of SFXN4 and SQOR on invasion and metastasis of osteosarcoma. (A, B) Wound healing assay revealing the impact of SFXN4 and SQOR knockdown on osteosarcoma cell invasiveness. (C, D) The invasive and migratory abilities of HOS and 143B cell lines were assessed by Transwell assays after SFXN4 and SQOR knockdown. ∗ p < 0.05, ∗∗ p < 0.01.
Figure 12. Knockdown of SFXN4 and SQOR affected osteosarcoma growth in vivo. 143B cells were transfected with lentivirus expressing either shNC, shSFXN4, or shSQOR. (A) Knockdown of SFXN4 significantly inhibited tumor growth. (B, C) Volume changes and weight of shNC and shSFXN4 tumors. (D) Knockdown of SQOR promoted tumor growth. (E, F) Volume changes and weight of shNC and shSQOR tumors. **p < 0.01.
4 Discussion
Hypoxia and lactate metabolism are two widely recognized characteristics of the tumor microenvironment (TME) that profoundly affect tumor progression, immune infiltration, and treatment response (33–38). This study has introduced a novel gene signature related to hypoxia and lactate metabolism (HLMRGS), representing the first identification of such a signature in osteosarcoma. Moreover, the research also underscored the clinical relevance of the HLMRGS concentrated on its latent force for predicting prognosis and immunotherapy response.
Previous studies primarily stressed on gene signatures related to hypoxia or lactate metabolism in osteosarcoma (39, 40). However, a gene signature that integrates both hypoxia and lactate metabolism for prognostic evaluation and treatment guidance in osteosarcoma has not yet been reported. The area under the curve (AUC) values of the HLMRGS surpassed those of existing models, indicating its superior predictive performance. Besides, compared to other models related to hypoxia or lactate metabolism, our AUC values stood out as the highest, showcasing exceptional predictive capability (41–47). This suggested that the HLMRGS was well suitable for osteosarcoma patients and served as a valuable tool for predicting clinical outcomes.
Immune cells are the cellular foundation of immunotherapy (48). The abundance and variety of tumor-infiltrating immune cells are intricately linked to therapeutic efficacy prediction and clinical outcomes (49). The immune cell composition of individuals with the same cancer may vary significantly inside the tumor microenvironment (TME), emphasizing the importance of characterizing immune infiltrates and their functional status for enhancing response rates (50). In our study, we observed that lower risk scores were closely correlated with more infiltration of crucial antitumor cells, including NK cells, CD4+ T cells, CD8+ T cells, cytotoxic T cells, and dendritic cells (DCs) within the TME. The cytotoxic activity of CD8+ T cells is a potent effector mechanism in tumor destruction (51), while CD4+ T cells are pivotal in initiating and coordinating innate and antigen-specific immune response (52). NK cells can recruit DCs to tumors to strengthen CD8+ T-cell response and improve the efficacy of immunotherapy by cooperating with T cells based on their complementary functions in tumor immunity (53). Besides, the low-risk group exhibited increased immune scores according to the ESTIMATE algorithm, indicating a more favorable TME enriched with abundant immune cells. It has been reported that osteosarcoma can be classified as an immune “hot” tumor or “cold” tumor or somewhere in between (54). Generally, “hot” tumors demonstrate higher response rates to immunotherapy (55). In our research, the tumors of low-risk patients were more likely to be classified as “hot” tumors more sensitive to immunotherapy due to more infiltration of crucial anti-tumor cells.
As a prevailing immunotherapy, immune checkpoint inhibitors (ICIs) treatment takes into account the genetic background of tumors by utilizing biomarker-based patient selection, thereby hopefully increasing the proportion of osteosarcoma patients who benefit from the ICIs treatment (5). In our study, we observed that the subset classified as “hot” exhibited an adaptive immune resistance mechanism characterized by the upregulation of T cell inhibitory immune-checkpoint proteins, specifically TIGIT (32). Besides, PDCD1LG2 blockade may become potential immunotherapeutic intervention targets enhancing the cytotoxic effects against osteosarcoma (56). Overall, low-risk patients were more likely to be categorized as individuals with “hot” tumors and harbored a better prognosis after ICIs treatment thanks to more infiltration of tumor-killing cells with reversible dysfunction.
Additionally, low-risk patients exhibited increased expression levels of most HLA genes, interleukin-2 (IL-2), and heightened activity in the M2-like tumor-associated macrophages (TAMs) pathway, as well as an elevated gene expression profile (GEP) score. The majority of “cold” tumors displayed decreased expression levels of HLA antigen-presenting genes, leading to the absence of TCR productive clonality and subsequent immune evasion (57). The combination of IL-2 and chemotherapy demonstrate a three-year survival rate exceeding 40% among patients (58). TAMs may create an immunosuppressive environment to participate in the malignant progression of osteosarcoma (59). As a novel gene signature biomarker for patient's response to ICIs, a higher GEP score implies worse treatment effects (60). These results were in line with the findings of the TIDE prediction analysis, indicating that low-risk individuals may derive greater benefits from immunotherapy interventions.
Most modeled genes included in the HLMRGS have been demonstrated to be connected to tumor microenvironment or immunotherapy efficiency across various cancers in prior studies. Here, our bioinformatic correlation analysis suggested that SQOR was positively correlated with the infiltration abundance of nearly all immune cells, with particularly strong correlations observed for NK T cell and Gamma delta T cell infiltration. NK T cells are a specialized subpopulation of T cells that expressed both T cell receptor and NK cell receptor (61). Beyond their direct cytotoxic effects on cancer cells, they also modulate the immune response by secreting a large number of cytokines and chemokines (e.g., IL-4, IFN-V, etc.) (62). Gamma delta T cells are recognized as the most tumor-killing immune cell; they not only induce apoptosis in tumor cells but also activate NK cells through 4-1BBL expression (63). Hence, we hypothesized that SQOR may activate the interaction network between NK T cells and Gamma belta T cells, thereby increasing the abundance of immune cell infiltration, altering the tumor microenvironment, and enhancing the patient’s response to immunotherapy. SFXN4 has been shown to be associated with multimodal immune cell infiltration in hepatocellular carcinoma (64). Unfortunately, our analysis showed that this correlation did not hold true in osteosarcoma. Nevertheless, the results of our in vivo and in vitro experiments supported a critical role of SFXN4 in tumor progression. However, the molecular mechanism of SFXN4 in the malignant progression of osteosarcoma needs to be further investigated. COL5A2 was found to enhance the infiltration abundance of tumor-associated macrophages while concurrently diminishing the population of CD8 T cells in prostate cancer (65). FAM162A has been incorporated into the prediction model of lung squamous cell carcinoma, and its expression has been taken into account for calculating the m6A score to predict the immunotherapy response of patients (66). UQCRB can ensure the normal function of mitochondrial complex III and rectify hypoglycemia and lactic acidosis in the tumor microenvironment of colorectal cancer (67). COX6A2 was reported to be incorporated into the oxidative phosphorylation gene model of osteosarcoma for predicting the prognosis of patients and the infiltration of immune cells (68). Last but not the least, we conducted in vivo experiments to investigate the effects of SQOR and SFXN4 on the growth of osteosarcoma. The results demonstrated that SQOR exhibited a tumor-suppressive function, whereas SFXN4 facilitated tumor growth. These findings aligned with our bioinformatics predictions, indicating that both SQOR and SFXN4 may serve as potential drug targets for future clinical treatment of osteosarcoma. This opened new avenues for targeted therapy aimed at osteosarcoma patients. However, further research is required to determine whether SQOR and SFXN4 influence the growth of osteosarcoma through the tumor microenvironment and tumor immunity.
Our study had several limitations. The signature was developed and validated through bioinformatics analyses using publicly available databases. Although we have assessed its reliability in external validation sets, a large-scale clinical trial is necessary for further validation. Moreover, additional research on the efficacy of immunotherapy and drug sensitivity assays is needed for clinical evaluation in future studies. Finally, the specific molecular mechanisms by which SQOR and SFXN4 affect the tumor microenvironment and immunotherapeutic response need to be further investigated.
5 Conclusions
In summary, this study developed and validated a powerful hypoxia- and lactate metabolism-related gene signature that effectively predicted patient’s prognosis and response to immunotherapy. In vivo and in vitro experiments demonstrated that the expression of SQOR could suppress the growth of osteosarcoma, whereas the expression of SFXN4 promoted it. Additionally, in vitro experiments also verified that SQOR inhibited the metastasis and invasiveness of osteosarcoma cells, while SFXN4 had the opposite effect. These findings will contribute to understanding the tumor microenvironment in osteosarcoma and guiding the selection of immune checkpoint inhibitors in osteosarcoma treatment strategies, thus improving individualized therapy and enlarging the population benefiting from immunotherapy.
Data availability statement
The original contributions presented in the study were included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Ethics statement
The studies involving humans were approved by Ethics Committee of First Hospital of Harbin Medical University (Batch Number: IRB-AF/SC-0402.0). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. The animal study was approved by The Ethics Committee of the First Hospital of Harbin Medical University (Batch Number: IRB-AF/SC-0402.0). The study was conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
YW: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft. XW: Data curation, Investigation, Resources, Validation, Writing – review & editing. YL: Formal analysis, Methodology, Validation, Writing – review & editing. JX: Software, Supervision, Writing – review & editing. JZ: Methodology, Resources, Software, Writing – review & editing. YZ: Supervision, Writing – review & editing. QQ: Data curation, Formal analysis, Investigation, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was supported by the Heilongjiang Disability Welfare Foundation (2020HX034) and Baiqiu’en Charity Foundation (2020IIT071).
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/fimmu.2024.1467052/full#supplementary-material.
Abbreviations
AUC, Area under curve; CAF, Cancer associated fibroblast; CI, Confidence interval; DAB, Diaminobenzidine; DCA, Decision curve analysis; DC, Dendritic cells; DEHRGs, Differentially expressed hypoxia-related genes; DELMRGs, Differentially expressed lactate metabolism-related genes; EGFR, Epidermal growth factor receptor; GC, Germinal center; GEO, Gene expression omnibus; GEP, Gene expression profile; GO, Gene ontology; GTEx, Genotype-tissue expression; HLMRGS, Hypoxia and lactate metabolism-related gene signature; HR, Hazard ratio; IC50, Half-maximal inhibitory concentration; ICIs, Immune checkpoint inhibitors; IFN, Interferon; IFNG, Interferon gamma; ImmuCellAI, Immune cell abundance identifier; KM, Kaplan-Meier; KEGG, Kyoto encyclopedia of genes and genomes; LASSO, Least absolute shrinkage and selection operator; MAC, Monocyte–derived macrophages; MALT, Mucosa-associated lymphoid tissue; MDSC, Myeloid-derived suppressor cell; MHC, Major histocompatibility complex; mo-TAM, Tumor-associated macrophages; MSI, Microsatellite instability; MSigDB, Molecular signatures database; NK, Natural killer; OS, Overall survival; PCA, Principal component analysis; ROC, Receiver operating characteristic; ssGSEA, Single sample gene set enrichment analysis; TARGET, Therapeutically applicable research to generate effective treatments; Th1, T helper 1; Th2, T helper 2; TIDE, Tumor immune dysfunction and exclusion; TISCH, Tumor immune single cell hub; TME, Tumor microenvironment; Treg, Regulatory T cells.
References
1. Caudill JS, Arndt CA. Diagnosis and management of bone Malignancy in adolescence. Adolesc Med State Art Rev. (2007) 18:62–78, ix.
2. Durfee RA, Mohammed M, Luu HH. Review of osteosarcoma and current management. Rheumatol Ther. (2016) 3:221–43. doi: 10.1007/s40744-016-0046-y
3. Mirabello L, Troisi RJ, Savage SA. Osteosarcoma incidence and survival rates from 1973 to 2004: data from the Surveillance, Epidemiology, and End Results Program. Cancer. (2009) 115:1531–43. doi: 10.1002/cncr.24121
4. Carrle D, Bielack SS. Current strategies of chemotherapy in osteosarcoma. Int Orthop. (2006) 30:445–51. doi: 10.1007/s00264-006-0192-x
5. Gill J, Gorlick R. Advancing therapy for osteosarcoma. Nat Rev Clin Oncol. (2021) 18:609–24. doi: 10.1038/s41571-021-00519-8
6. Misaghi A, Goldin A, Awad M, and Kulidjian AA. Osteosarcoma: a comprehensive review. SICOT J. (2018) 4:12. doi: 10.1051/sicotj/2017028
7. Tsukamoto S, Errani C, Angelini A, and Mavrogenis AF. Current treatment considerations for osteosarcoma metastatic at presentation. Orthopedics. (2020) 43:e345–58. doi: 10.3928/01477447-20200721-05
8. Allison DC, Carney SC, Ahlmann ER, Hendifar A, Chawla S, Fedenko A, et al. A meta-analysis of osteosarcoma outcomes in the modern medical era. Sarcoma. (2012) 2012:704872. doi: 10.1155/2012/704872
9. Luetke A, Meyers PA, Lewis I, and Juergens H. Osteosarcoma treatment - where do we stand? A state of the art review. Cancer Treat Rev. (2014) 40:523–32. doi: 10.1016/j.ctrv.2013.11.006
10. Saraf AJ, Fenger JM, Roberts RD. Osteosarcoma: accelerating progress makes for a hopeful future. Front Oncol. (2018) 8:4. doi: 10.3389/fonc.2018.00004
11. Botter SM, Neri D, Fuchs B. Recent advances in osteosarcoma. Curr Opin Pharmacol. (2014) 16:15–23. doi: 10.1016/j.coph.2014.02.002
12. Jing X, Yang F, Shao C, Wei K, Xie M, Shen H, et al. Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol Cancer. (2019) 18:157. doi: 10.1186/s12943-019-1089-9
13. Korbecki J, Kojder K, Kapczuk P, Kupnicka P, Gawrońska-Szklarz B, Gutowska I, et al. The effect of hypoxia on the expression of CXC chemokines and CXC chemokine receptors-A review of literature. Int J Mol Sci. (2021) 22(2):843. doi: 10.3390/ijms22020843
14. Rankin EB, Nam JM, Giaccia AJ. Hypoxia: signaling the metastatic cascade. Trends Cancer. (2016) 2:295–304. doi: 10.1016/j.trecan.2016.05.006
15. Liberti MV, Locasale JW. The warburg effect: how does it benefit cancer cells? Trends Biochem Sci. (2016) 41:211–8. doi: 10.1016/j.tibs.2015.12.001
16. San-Millan I, Brooks GA. Reexamining cancer metabolism: lactate production for carcinogenesis could be the purpose and explanation of the Warburg Effect. Carcinogenesis. (2017) 38:119–33. doi: 10.1093/carcin/bgw127
17. Lee DC, Tang Z, Huang H, Zhou G, Cui C, Weng Y, et al. A lactate-induced response to hypoxia. Cell. (2015) 161:595–609. doi: 10.1016/j.cell.2015.03.011
18. Zhang D, Tang Z, Huang H, Zhou G, Cui C, Weng Y, et al. Metabolic regulation of gene expression by histone lactylation. Nature. (2019) 574:575–80. doi: 10.1038/s41586-019-1678-1
19. Brand A, Singer K, Koehl GE, Kolitzus M, Schoenhammer G, Thiel A, et al. LDHA-associated lactic acid production blunts tumor immunosurveillance by T and NK cells. Cell Metab. (2016) 24:657–71. doi: 10.1016/j.cmet.2016.08.011
20. Brown TP, Ganapathy V. Lactate/GPR81 signaling and proton motive force in cancer: Role in angiogenesis, immune escape, nutrition, and Warburg phenomenon. Pharmacol Ther. (2020) 206:107451. doi: 10.1016/j.pharmthera.2019.107451
21. Vegran F, Boidot R, Michiels C, Sonveaux P, Feron O. Lactate influx through the endothelial cell monocarboxylate transporter MCT1 supports an NF-kappaB/IL-8 pathway that drives tumor angiogenesis. Cancer Res. (2011) 71:2550–60. doi: 10.1158/0008-5472
22. Yang X, Zhang W, Zhu W. Profiling of immune responses by lactate modulation in cervical cancer reveals key features driving clinical outcome. Heliyon. (2023) 9:e14896. doi: 10.1016/j.heliyon.2023.e14896
23. Riera-Domingo C, Audigé A, Granja S, Cheng WC, Ho PC, Baltazar F, et al. Immunity, hypoxia, and metabolism-the menage a trois of cancer: implications for immunotherapy. Physiol Rev. (2020) 100:1–102. doi: 10.1152/physrev.00018
24. Wu Q, You L, Nepovimova E, Heger Z, Wu W, Kuca K, et al. Hypoxia-inducible factors: master regulators of hypoxic tumor immune escape. J Hematol Oncol. (2022) 15:77. doi: 10.1186/s13045-022-01292-6
25. Boutilier AJ, Elsawa SF. Macrophage polarization states in the tumor microenvironment. Int J Mol Sci. (2021) 22(13):6995. doi: 10.3390/ijms22136995
26. Mortezaee K, Majidpoor J. The impact of hypoxia on immune state in cancer. Life Sci. (2021) 286:120057. doi: 10.1016/j.lfs.2021.120057
27. Feng Q, Liu Z, Yu X, Huang T, Chen J, Wang J, et al. Lactate increases stemness of CD8 + T cells to augment anti-tumor immunity. Nat Commun. (2022) 13:4981. doi: 10.1038/s41467-022-32521-8
28. Kumagai S, Koyama S, Itahashi K, Tanegashima T, Lin YT, Togashi Y, et al. Lactic acid promotes PD-1 expression in regulatory T cells in highly glycolytic tumor microenvironments. Cancer Cell. (2022) 40:201–218 e9. doi: 10.1016/j.ccell.2022.01.001
29. Murray PJ. Macrophage polarization. Annu Rev Physiol. (2017) 79:541–66. doi: 10.1146/annurev-physiol-022516-034339
30. Hu D, Zhang Z, Zhang Y, Huang K, Li X. Identification of immune related molecular subtypes and prognosis model for predicting prognosis, drug resistance in cervical squamous cell carcinoma. Front Genet. (2023) 14:1137995. doi: 10.3389/fgene.2023.1137995
31. Zhao C, Xiong K, Adam A, Ji Z, Li X. Necroptosis identifies novel molecular phenotypes and influences tumor immune microenvironment of lung adenocarcinoma. Front Immunol. (2022) 13:934494. doi: 10.3389/fimmu.2022.934494
32. Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, et al. Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma. Nat Commun. (2020) 11:6322. doi: 10.1038/s41467-020-20059-6
33. Hou Y, Guo H, Cao C, Li X, Hu B, Zhu P, et al. Single-cell triple omics sequencing reveals genetic, epigenetic, and transcriptomic heterogeneity in hepatocellular carcinomas. Cell Res. (2016) 26:304–19. doi: 10.1038/cr.2016.23
34. Li Q, Wang R, Yang Z, Li W, Yang J, Wang Z, et al. Molecular profiling of human non-small cell lung cancer by single-cell RNA-seq. Genome Med. (2022) 14:87. doi: 10.1186/s13073-022-01089-9
35. Li W, Tanikawa T, Kryczek I, Xia H, Li G, Wu K, et al. Aerobic glycolysis controls myeloid-derived suppressor cells and tumor immunity via a specific CEBPB isoform in triple-negative breast cancer. Cell Metab. (2018) 28:87–103 e6. doi: 10.1016/j.cmet.2018.04.022
36. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. (2011) 27:1739–40. doi: 10.1093/bioinformatics/btr260
37. Wu F, Fan J, He Y, Xiong A, Yu J, Li Y, et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun. (2021) 12:2540. doi: 10.1038/s41467-021-22801-0
38. Zito Marino F, Bianco R, Accardo M, Ronchi A, Cozzolino I, Morgillo F, et al. Molecular heterogeneity in lung cancer: from mechanisms of origin to clinical implications. Int J Med Sci. (2019) 16:981–9. doi: 10.7150/ijms.34739
39. Wu Z, Han T, Su H, Xuan J, Wang X. Comprehensive analysis of fatty acid and lactate metabolism-related genes for prognosis value, immune infiltration, and therapy in osteosarcoma patients. Front Oncol. (2022) 12:934080. doi: 10.3389/fonc.2022.934080
40. Zhang W, Lyu P, Andreev D, Jia Y, Zhang F, Bozec A. Hypoxia-immune-related microenvironment prognostic signature for osteosarcoma. Front Cell Dev Biol. (2022) 10:974851. doi: 10.3389/fcell.2022.974851
41. Chen D, Huang H, Zang L, Gao W, Zhu H, Yu X, et al. Development and verification of the hypoxia- and immune-associated prognostic signature for pancreatic ductal adenocarcinoma. Front Immunol. (2021) 12:728062. doi: 10.3389/fimmu.2021.728062
42. Huang L, Zeng X, Liang W, Chen J, Zhong C, Cai W, et al. Dissecting the role of lactate metabolism LncRNAs in the progression and immune microenvironment of osteosarcoma. Transl Oncol. (2023) 36:101753. doi: 10.1016/j.tranon.2023.101753
43. Li J, Qiao H, Wu F, Sun S, Feng C, Li C, et al. A novel hypoxia- and lactate metabolism-related signature to predict prognosis and immunotherapy responses for breast cancer by integrating machine learning and bioinformatic analyses. Front Immunol. (2022) 13:998140. doi: 10.3389/fimmu.2022.998140
44. Li Y, Mo H, Wu S, Liu X, Tu K. A novel lactate metabolism-related gene signature for predicting clinical outcome and tumor microenvironment in hepatocellular carcinoma. Front Cell Dev Biol. (2021) 9:801959. doi: 10.3389/fcell.2021.801959
45. Shi R, Bao X, Unger K, Sun J, Lu S, Manapov F, et al. Identification and validation of hypoxia-derived gene signatures to predict clinical outcomes and therapeutic responses in stage I lung adenocarcinoma patients. Theranostics. (2021) 11:5061–76. doi: 10.7150/thno.56202
46. Zhang B, Tang B, Gao J, Li J, Kong L, Qin L. A hypoxia-related signature for clinically predicting diagnosis, prognosis and immune microenvironment of hepatocellular carcinoma patients. J Transl Med. (2020) 18:342. doi: 10.1186/s12967-020-02492-9
47. Zhong J, Shen X, Zhou J, Yu H, Wang B, Sun J, et al. Development and validation of a combined hypoxia and ferroptosis prognostic signature for breast cancer. Front Oncol. (2023) 13:1077342. doi: 10.3389/fonc.2023.1077342
48. Rui R, Zhou L, He S. Cancer immunotherapies: advances and bottlenecks. Front Immunol. (2023) 14:1212476. doi: 10.3389/fimmu.2023.1212476
49. Duan Q, Zhang H, Zheng J, Zhang L. Turning Cold into Hot: Firing up the Tumor Microenvironment. Trends Cancer. (2020) 6:605–18. doi: 10.1016/j.trecan.2020.02.022
50. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. (2013) 39:782–95. doi: 10.1016/j.immuni.2013.10.003
51. Baxevanis CN, Perez SA, Papamichail M. Cancer immunotherapy. Crit Rev Clin Lab Sci. (2009) 46:167–89. doi: 10.1080/10408360902937809
52. Speiser DE, Chijioke O, Schaeuble K, Münz C. CD4(+) T cells in cancer. Nat Cancer. (2023) 4:317–29. doi: 10.1038/s43018-023-00521-2
53. Kyrysyuk O, Wucherpfennig KW. Designing cancer immunotherapies that engage T cells and NK cells. Annu Rev Immunol. (2023) 41:17–38. doi: 10.1146/annurev-immunol-101921-044122
54. Yang H, Zhao L, Zhang Y, Li FF. A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma. Cancer Med. (2021) 10:5696–711. doi: 10.1002/cam4.4117
55. Zemek RM, De Jong E, Chin WL, Schuster IS, Fear VS, Casey TH, et al. Sensitization to immune checkpoint blockade through activation of a STAT1/NK axis in the tumor microenvironment. Sci Transl Med. (2019) 11(501):eaav7816. doi: 10.1126/scitranslmed.aav7816
56. McEachron TA, Triche TJ, Sorenson L, Parham DM, Carpten JD. Profiling targetable immune checkpoints in osteosarcoma. Oncoimmunology. (2018) 7:e1475873. doi: 10.1080/2162402x.2018.1475873
57. Wu CC, Beird HC, Andrew Livingston J, Advani S, Mitra A, Cao S, et al. Immuno-genomic landscape of osteosarcoma. Nat Commun. (2020) 11:1008. doi: 10.1038/s41467-020-14646-w
58. Meazza C, Cefalo G, Massimino M, Daolio P, Pastorino U, Scanagatta P, et al. Primary metastatic osteosarcoma: results of a prospective study in children given chemotherapy and interleukin-2. Med Oncol. (2017) 34:191. doi: 10.1007/s12032-017-1052-9
59. Huang Q, Liang X, Ren T, Huang Y, Zhang H, Yu Y, et al. The role of tumor-associated macrophages in osteosarcoma progression - therapeutic implications. Cell Oncol (Dordr). (2021) 44:525–39. doi: 10.1007/s13402-021-00598-w
60. Wang Y, Tong Z, Zhang W, Zhang W, Buzdin A, Mu X, et al. FDA-approved and emerging next generation predictive biomarkers for immune checkpoint inhibitors in cancer patients. Front Oncol. (2021) 11:683419. doi: 10.3389/fonc.2021.683419
61. Crouse J, Xu HC, Lang PA, Oxenius A. NK cells regulating T cell responses: mechanisms and outcome. Trends Immunol. (2015) 36:49–58. doi: 10.1016/j.it.2014.11.001
62. Liu X, Li L, Si F, Huang L, Zhao Y, Zhang C, et al. NK and NKT cells have distinct properties and functions in cancer. Oncogene. (2021) 40:4521–37. doi: 10.1038/s41388-021-01880-9
63. Mamedov MR, Vedova S, Freimer JW, Sahu AD, Ramesh A, Arce MM, et al. CRISPR screens decode cancer cell pathways that trigger gammadelta T cell detection. Nature. (2023) 621:188–95. doi: 10.1038/s41586-023-06482-x
64. Du Z, Zhang Z, Han X, Xie H, Yan W, Tian D, et al. Comprehensive analysis of sideroflexin 4 in hepatocellular carcinoma by bioinformatics and experiments. Int J Med Sci. (2023) 20:1300–15. doi: 10.7150/ijms.86990
65. Ren X, Chen X, Fang K, Zhang X, Wei X, Zhang T, et al. COL5A2 promotes proliferation and invasion in prostate cancer and is one of seven gleason-related genes that predict recurrence-free survival. Front Oncol. (2021) 11:583083. doi: 10.3389/fonc.2021.583083
66. Li P, Xiong P, Li X, Zhang X, Chen X, Zhang W, et al. Tumor microenvironment characteristics and prognostic role of m(6)A modification in lung squamous cell carcinoma. Heliyon. (2024) 10:e26851. doi: 10.1016/j.heliyon.2024.e26851
67. Kim HC, Chang J, Lee HS, Kwon HJ. Mitochondrial UQCRB as a new molecular prognostic biomarker of human colorectal cancer. Exp Mol Med. (2017) 49:e391. doi: 10.1038/emm.2017.152
Keywords: hypoxia, lactate metabolism, osteosarcoma, gene signature, immune checkpoint inhibitors, immunotherapy, tumor microenvironment, single-cell sequencing
Citation: Wang Y, Wang X, Liu Y, Xu J, Zhu J, Zheng Y and Qi Q (2024) A novel hypoxia- and lactate metabolism-related prognostic signature to characterize the immune landscape and predict immunotherapy response in osteosarcoma. Front. Immunol. 15:1467052. doi: 10.3389/fimmu.2024.1467052
Received: 19 July 2024; Accepted: 10 October 2024;
Published: 06 November 2024.
Edited by:
Fernando Torres Andón, Institute of Biomedical Research of A Coruña (INIBIC), SpainReviewed by:
Yuqin Tang, Henan Provincial People’s Hospital, ChinaTianyun Qiao, Fourth Military Medical University, China
Copyright © 2024 Wang, Wang, Liu, Xu, Zhu, Zheng and Qi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Quan Qi, qiquan@hrbmu.edu.cn
†These authors have contributed equally to this work and share first authorship