Prognostic and Immunomodulatory Roles of BNIP3 in Osteosarcoma Revealed by Integrated Single‐Cell and Bulk Transcriptomic Profiling

TL;DR

This study identifies BNIP3 as a harmful gene in osteosarcoma, linked to poor survival and immune evasion, suggesting it could be a new therapeutic target.

Contribution

The study reveals BNIP3 as a novel prognostic and immunomodulatory gene in osteosarcoma using integrated single-cell and bulk RNA sequencing.

Findings

BNIP3 is the only significant prognostic gene associated with poor survival in osteosarcoma patients.

High BNIP3 expression correlates with reduced immune cell infiltration and immune-related pathways.

BNIP3-high patients show reduced sensitivity to potential therapeutic agents.

Abstract

Osteosarcoma is a highly aggressive bone tumor with a complex tumor microenvironment (TME) that contributes to its progression and therapeutic resistance. In this study, we integrated single‐cell RNA sequencing (scRNA‐seq) and bulk RNA‐seq datasets to characterize the TME and identify key prognostic genes in osteosarcoma. Using scRNA‐seq data from 16 osteosarcoma samples, we defined eight major cell types within the TME and performed functional enrichment analyses. Through weighted gene co‐expression network analysis (WGCNA) focused on an inflammation‐related gene signature, we identified the yellow module as the most correlated with inflammation. By intersecting tumor‐upregulated genes with WGCNA‐derived genes, we identified BNIP3 as the only significant prognostic gene associated with poor survival in both the TARGET and GSE21257 cohorts. Functional annotation revealed that high BNIP3 expression is negatively correlated with immune‐related pathways and immune cell infiltration, including T cells, B cells, NK cells, and neutrophils. Additionally, BNIP3‐high patients exhibited a reduced sensitivity to several potential therapeutic agents. Our findings highlight BNIP3 as a hazardous gene in osteosarcoma, with important roles in immune evasion and prognosis, suggesting its potential as a therapeutic target.

Figures

Functional annotation of microenvironment cells in osteosarcoma at the scRNA‐seq level. Heat map shows the pathway enrichment results of microenvironment cells based on GOBP, KEGG, and WikiPathway terms.

Drug prediction of BNIP3 in osteosarcoma. The drug sensitivity of nine drugs in BNIP3‐based groups in the TARGET cohort.

Functional annotation of BNIP3 in osteosarcoma. GSEA of 12 GO terms related to BNIP3 in the TARGET cohort.

Pan‐cancer survival analysis of BNIP3 across 33 tumor types.

Full text

1. Introduction

Osteosarcoma, the most prevalent primary malignant bone tumor in children and adolescents, presents a significant clinical challenge due to its aggressive nature, high metastatic potential, and unsatisfactory long‐term survival rates despite advances in multimodal therapies, including surgery and chemotherapy [1]. Its pronounced heterogeneity, both at the genetic and cellular levels, contributes to therapeutic resistance and disease relapse, underscoring the urgent need for a deeper molecular understanding to identify novel prognostic markers and therapeutic targets [2].

The tumor microenvironment (TME) is now recognized as a critical orchestrator of osteosarcoma progression, metastasis, and treatment response [3]. It is a complex ecosystem comprising tumor cells, immune cells (e.g., T cells, macrophages, and myeloid‐derived suppressor cells), stromal cells (including fibroblasts and endothelial cells), and noncellular components. Dynamic interactions within the TME can foster an immunosuppressive milieu, promote angiogenesis, and facilitate tumor cell invasion [4]. Among the hallmarks of cancer, inflammation plays a particularly pivotal role in shaping the TME, influencing immune cell recruitment, activation, and function, thereby affecting patient prognosis [5]. Immunotherapies, such as immune checkpoint inhibitors, have shown limited efficacy in osteosarcoma, which is often characterized as an immunologically “cold” tumor with a low immune cell infiltration and a suppressive TME. This highlights the critical need to identify key drivers of immune evasion. However, a precise, high‐resolution dissection of the inflammatory landscape and its key regulatory genes within the osteosarcoma TME remains incomplete.

The advent of single‐cell RNA sequencing (scRNA‐seq) technology has revolutionized our ability to deconvolute this complexity, enabling the characterization of cellular diversity, identification of rare cell populations, and inference of cell–cell communication at an unprecedented resolution [6]. Concurrently, bulk transcriptomic analyses of large patient cohorts provide robust statistical power for linking molecular features to clinical outcomes [7]. Integrating these complementary approaches offers a powerful strategy to bridge cellular phenotypes with patient‐level prognosis and to discover genes that are not only differentially expressed but also functionally central within disease‐relevant networks.

In this study, we leveraged integrated scRNA‐seq data from 16 osteosarcoma samples to comprehensively map the cellular composition and functional states of the TME. We combined this with bulk RNA‐seq data from the TARGET and GEO databases to perform a systems biology analysis. Specifically, we applied weighted gene co‐expression network analysis (WGCNA) to define gene modules co‐regulated with an inflammation‐related gene signature. By intersecting module genes with those upregulated in tumor cells, we aim to pinpoint key drivers situated at the nexus of tumor cell intrinsic properties and the inflammatory microenvironment. Our analysis identified BNIP3 (BCL2 interacting Protein 3), a gene involved in autophagy, apoptosis, and mitochondrial function, as a top candidate [8]. We subsequently validated its prognostic significance and conducted in‐depth explorations of its association with the immune landscape and drug sensitivity in osteosarcoma. This work establishes BNIP3 as a critical immunomodulatory prognostic factor, providing new insights into the pathogenesis of osteosarcoma and highlighting a potential target for therapeutic intervention.

2. Materials and Methods

2.1. Data Collection

The scRNA‐seq datasets, GSE152048 [9] and GSE162454 [10], of osteosarcoma were collected from the GEO database. The bulk RNA sequencing dataset of osteosarcoma was collected from the TARGET database and the GEO database (GSE21257) [11], respectively. The bulk RNA sequencing dataset of pan‐cancer data was collected from the TCGA database.

2.2. Data Analysis

The “inflammation”‐related gene list was obtained from the HALLMARK database [12]. The R package Seurat [13] was implemented to process the scRNA‐seq dataset. We systematically addressed batch effects using the Harmony algorithm. For quality control, cells were filtered to retain those with 200 < nFeature_RNA < 6000 and mitochondrial content of < 20%. The R package InferCNV was used to define tumor cells. Uniform manifold approximation and projection (UMAP) was used to present the distribution of cells. We used manual cell annotation based on canonical markers. The pathway enrichment analysis based on GOBP, KEGG, and WikiPathway terms of microenvironment cells was conducted using the R package SCP. The R package WGCNA was used to identify the “inflammation”‐related genes [14]. The connection between BNIP3 and stromal score, immune score, ESTIMATE score, and tumor purity using the ESTIMATE algorithm [15, 16]. The correlation between BNIP3 and immune cells was predicted using Porpimol′s study [17]. The drug prediction of BNIP3 was performed using the R package pRRophetic [18] based on the pharmacogenomic data from the GDSC database.

2.3. Statistical Analysis

All statistical analyses were performed using R and relevant bioinformatics packages. All p values were two‐sided, and statistical significance was set at p < 0.05 unless otherwise specified.

3. Results

3.1. Definition and Functional Annotation of Microenvironment Cells in Osteosarcoma at the scRNA‐seq Level

The integration of 16 osteosarcoma samples is shown in Figure 1a. The immune cells and nonimmune cells in the microenvironment of osteosarcoma were defined (Figure 1b). Subsequently, eight cell types, including tumor cells, T/NK cells, myeloid cells, fibroblasts, endothelial cells, osteoclasts, pericytes, and B cells, in the microenvironment of osteosarcoma were defined (Figure 1c). The pathway enrichment results based on GOBP, KEGG, and WikiPathway terms indicated that extracellular matrix organization, ECM–receptor interaction, and ossification were related to tumor cells (Figure 2). The enrichment of the “gastric cancer network” likely reflects the shared oncogenic signaling pathways common across cancers rather than tissue‐specific biology.

Figure 1. Definition of microenvironment cells in osteosarcoma at the scRNA‐seq level. (a) UMAP shows the integration of 16 osteosarcoma samples. (b) UMAP shows the major cell categories in the microenvironment of osteosarcoma. (c) UMAP shows the minor cell categories in the microenvironment of osteosarcoma.(a)(b)(c)

Functional annotation of microenvironment cells in osteosarcoma at the scRNA‐seq level. Heat map shows the pathway enrichment results of microenvironment cells based on GOBP, KEGG, and WikiPathway terms.

3.2. WGCNA on “Inflammation” Signature in Osteosarcoma

The optimal scale‐free topology model was determined at the soft threshold of four (Figure 3a). The dendrogram of the hierarchy of gene modules linked to the “inflammation” signature is shown in Figure 3b. The “inflammation” signature was found to correlate most with the yellow module (R = 0.73, P = 4e − 15) (Figure 3c). A significant association between gene significance and module membership in the yellow module is shown in Figure 3d.

Figure 3WGCNA on “inflammation” signature in osteosarcoma. (a) The association between scale‐free topology model fit, mean connectivity, and soft threshold. (b) Dendrogram shows the hierarchy of gene modules linked to the “inflammation” signature. (c) The correlation between “inflammation” signature and gene modules. (d) The association between gene significance and module membership in the yellow module.(a)(b)(c)(d)

3.3. Identification of BNIP3 as a Hazardous Gene in Osteosarcoma

Tumor‐upregulated genes were defined by comparing the tumor cell cluster with all other nonmalignant cell clusters in the integrated scRNA‐seq data. Two hundred ninety‐eight upregulated genes in tumor cells compared with other microenvironment cells were identified. Four hundred ninety‐eight yellow module‐derived genes were extracted. Eight intersected genes between genes upregulated in tumor cells and WGCNA‐derived genes were finally identified (Figure 4a). Univariate Cox regression analysis confirmed BNIP3 to be the only prognostic intersected gene (p < 0.05) (Figure 4b). Accordingly, osteosarcoma patients with high BNIP3 expression were associated with worse survival outcomes in the TARGET cohort with a p < 0.0001 (Figure 4c) and the GSE21257 cohort with a p = 0.031 (Figure 4d).

Figure 4. Identification of BNIP3 as a hazardous gene in osteosarcoma. (a) Venn plot shows the intersected genes between genes upregulated in tumor cells and WGCNA‐derived genes. (b) Univariate Cox regression analysis on eight intersected genes. (c) Survival curves of BNIP3‐based groups in the TARGET cohort. (d) Survival curves of BNIP3‐based groups in the GSE21257 cohort.(a)(b)(c)(d)

3.4. Drug Prediction of BNIP3 in Osteosarcoma

The drug sensitivity of nine drugs, including Navitoclax (1011), GSK1904529A (1093), cyclophosphamide (1512), Pevonedistat (1529), AGI‐6780 (1634), TAF1_5496 (1732), JAK_8517 (1739), OF‐1 (1853), and sepantronium bromide (1941), was significantly lower in osteosarcoma patients with high BNIP3 expression in the TARGET cohort (Figure 5). These results represent computational predictions that require experimental validation.

Drug prediction of BNIP3 in osteosarcoma. The drug sensitivity of nine drugs in BNIP3‐based groups in the TARGET cohort.

3.5. Functional Annotation of BNIP3 in Osteosarcoma

Twelve GO‐based pathways, including T cell activity, B cell activity, NK cell activity, and immune response, were significantly negatively related to BNIP3 in the TARGET cohort (Figure 6).

Functional annotation of BNIP3 in osteosarcoma. GSEA of 12 GO terms related to BNIP3 in the TARGET cohort.

3.6. Immune Characteristics of BNIP3 in Osteosarcoma

BNIP3 exhibited a significantly negative correlation with stromal score, immune score, ESTIMATE score, and tumor purity (Figure 7a). Additionally, immune cells, including T cells, B cells, NK cells, and neutrophils, were all negatively correlated with BNIP3 expression (Figure 7b).

Figure 7. Immune characteristics of BNIP3 in osteosarcoma. (a) The correlation between BNIP3 and microenvironment scores. (b) Heat map shows the correlation between BNIP3 and immune cells.(a)(b)

3.7. Pan‐Cancer Survival Analysis of BNIP3 Across 33 Tumor Types

To investigate whether the prognostic value of BNIP3 observed in osteosarcoma is consistent across other malignancies, we extended our analysis to a pan‐cancer cohort using the TCGA database (Figure 8). The forest plot revealed that BNIP3 expression significantly correlates with overall survival in multiple cancer types. Consistent with our findings in the bone microenvironment (SARC), BNIP3 was identified as a significant risk factor in other major solid tumors, including breast invasive carcinoma (BRCA) and skin cutaneous melanoma (SKCM). This suggests that the adverse prognostic role of BNIP3 is not limited to osteosarcoma but may share underlying mechanisms with these aggressive epithelial and melanocytic tumors.

Pan‐cancer survival analysis of BNIP3 across 33 tumor types.

4. Discussion

Our study provides a comprehensive transcriptomic analysis of the osteosarcoma TME at single‐cell resolution, identifying BNIP3 as a novel prognostic and immunomodulatory gene. The integration of scRNA‐seq data allowed us to define eight major cell types within the TME, including tumor cells, immune subsets, and stromal components, highlighting the cellular complexity of osteosarcoma. Through WGCNA, we identified an inflammation‐related gene signature associated with the yellow module, which was enriched for genes involved in immune and stromal regulation. BNIP3 emerged as a key intersection gene between tumor‐upregulated genes and the inflammation‐associated module, and it was the only gene significantly associated with poor survival in two independent cohorts. This underscores its potential as a robust prognostic biomarker in osteosarcoma [19].

The exclusive significance of BNIP3 among the intersected candidates suggests it occupies a critical functional nexus‐linking tumor cell biology and inflammatory signaling. Functionally, BNIP3 expression was negatively correlated with multiple immune‐related pathways and immune cell infiltration, suggesting its role in creating an immunosuppressive TME. The negative association with stromal, immune, and ESTIMATE scores further supports the idea that BNIP3 may contribute to a “cold” tumor immune phenotype, which is often linked to resistance to immunotherapy [20]. We hypothesize that BNIP3‐mediated autophagy and mitophagy in osteosarcoma cells may facilitate immune evasion through several potential mechanisms: (1) by reducing mitochondrial reactive oxygen species (ROS) and associated immunogenic cell death; (2) by limiting the release of damage‐associated molecular patterns (DAMPs) and mitochondrial DNA that can stimulate antitumor immunity; and (3) by promoting metabolic adaptations that alter the immunomodulatory secretome of tumor cells, collectively impairing immune cell recruitment and function. These findings align with previous studies implicating BNIP3 in autophagy and apoptosis, but extend its role to immune modulation in osteosarcoma.

Moreover, drug sensitivity analysis revealed that high BNIP3 expression correlates with reduced sensitivity to several agents, including Navitoclax and cyclophosphamide, indicating potential mechanisms of chemoresistance. For Navitoclax, BNIP3′s interaction with BCL‐2 family proteins and its prosurvival mitophagic function may provide an escape pathway. For cyclophosphamide, which induces oxidative stress, BNIP3‐mediated clearance of damaged mitochondria could enhance tumor cell survival. This highlights the need for targeted therapeutic strategies that could overcome BNIP3‐mediated resistance [21, 22].

It is important to note that BNIP3′s role in cancer is context‐dependent, acting as either a prosurvival or prodeath factor in different malignancies. In osteosarcoma, our consistent findings across cohorts suggest its prosurvival, adaptive functions likely dominate, promoting tumor resilience and a suppressive TME. The observed association between high BNIP3 and poor prognosis in osteosarcoma appears to be consistent with its hazardous role in some other cancers like glioblastoma, though it can be favorable in contexts like colorectal cancer, underscoring the importance of the tissue‐specific microenvironment.

Limitations of this study include its retrospective nature and reliance on publicly available datasets. Furthermore, mRNA expression may not always correlate perfectly with protein activity, and the drug sensitivity predictions are based on pan‐cancer cell line data, which may not fully capture osteosarcoma‐specific biology. Functional validation through in vitro and in vivo experiments is required to confirm the mechanistic role of BNIP3 in osteosarcoma progression and immune evasion. Additionally, the clinical translation of BNIP3 as a biomarker or therapeutic target warrants further investigation in prospective cohorts. Future studies should also evaluate whether BNIP3 expression is an independent prognostic factor when accounting for clinical variables like metastasis.

5. Conclusion

In summary, our integrative multiomics analysis identifies BNIP3 as a critical prognostic biomarker and immunomodulatory gene in osteosarcoma. High BNIP3 expression is associated with poor survival, reduced immune infiltration, and altered drug sensitivity, highlighting its potential as a therapeutic target. Although no specific BNIP3 inhibitors are currently in clinical development, our findings highlight the need for future research to explore strategies targeting BNIP3 or its downstream pathways to overcome therapy resistance. These findings provide new insights into the inflammatory landscape of osteosarcoma and pave the way for future studies aimed at targeting BNIP3 to improve clinical outcomes.

Author Contributions

Kangjun Yu and Li Li designed the study. Jiamiao Li and Wei Wang conducted analyses and drafted the paper.

Funding

No funding was received for this manuscript.

Disclosure

All authors critically reviewed and approved the final manuscript.

Consent

The authors have nothing to report.

Conflicts of Interest

The authors declare no conflicts of interest.