New insights into the role of mitophagy related gene affecting the metastasis of osteosarcoma through scRNA-seq and CRISPR-Cas9 genome editing

Abstract

Background

Osteosarcoma (OSA), the most common primary bone malignancy, poses significant challenges due to its aggressive nature and propensity for metastasis, especially in adolescents. Mitophagy analysis can help identify new therapeutic targets and combined treatment strategies.

Methods

This study integrates single-cell sequencing (scRNA-seq) data and bulk-seq to identify mitophagy-related genes (MRGs) associated with the progression of OSA metastasis and analyze their clinical significance. scRNA-seq data elucidates the relationship between mitophagy and OSA metastasis, employing “CellChat” R package to explore intercellular communications and report on hundreds of ligand-receptor interactions. Subsequently, the combination of bulk-seq and CRISPR-Cas9 gene editing identifies mitophagy-related biomarker associated with metastatic prognosis. Finally, validation of the relationship between mitophagy and OSA metastasis is achieved through cellular biology experiments and animal studies.

Results

The distinct mitophagy activity of various mitochondria manifests in diverse spatial localization, cellular developmental trajectories, and intercellular interactions. OSA tissue exhibits notable heterogeneity in mitophagy within osteoblastic OSA cells. However, high mitophagy activity correlates consistently with high metastatic potential. Subsequently, we identified three critical genes associated with mitophagy in OSA, namely RPS27A, TOMM20 and UBB. According to the aforementioned queue of genes, we have constructed a mitophagy_score (MIP_score). We observed that it consistently predicts patient prognosis in both internal and external datasets, demonstrating strong robustness and stability. Furthermore, we have found that MIP_score can also guide chemotherapy, with varying sensitivities to chemotherapeutic agents based on different MIP_score. It is noteworthy that, through the integration of CRISPR-Cas9 genome-wide screening and validation via cellular and animal experiments, we have identified RPS27A as a potential novel biomarker for OSA.

Conclusions

Our comprehensive analysis elucidated the profile of mitophagy throughout the OSA metastasis process, forming the basis for a mitophagy-related prognostic model that addresses clinical outcomes and drug sensitivity following OSA metastasis. Additionally, an online interactive platform was established to assist clinicians in decision-making (https://mip-score.shinyapps.io/labtan/). These findings lay the groundwork for developing targeted therapies aimed at improving the prognosis of OSA patients.

Introduction

Osteosarcoma (OSA) is the most common primary bone malignancy, with incidence peaks occurring during adolescence [1, 2]. It is believed to develop at some point during the differentiation process of mesenchymal stem cells (MSCs) into pre-osteoblasts or from osteoblast precursors [3]. There are three major subtypes of conventional OSA that reflect the predominant type of matrix within the tumor, including osteoblastic (76–80%), chondroblast (10–13%) and fibroblastic (10%) [4]. OSA often invades and metastasizes to other tissues (especially the lungs), with up to 15–20% of patients presenting with metastatic lesions at the initial diagnosis [5]. Studies have shown that the five-year survival rate of patients with carcinoma in situ who undergo surgery and chemotherapy is about 65%, compared with only 25% for patients with metastatic cancer [6, 7]. Extensive research indicates a tight correlation between the metastasis of OSA and molecular mechanisms involving complex karyotypes, genomic instability, and intricate protein signaling pathway interactions [8]. Henceforth, exploring the molecular mechanisms underpinning the metastasis of OSA and identifying novel biomarkers for diagnosis and treatment embodies an innovative therapeutic strategy.

Mitophagy is a selective autophagic process crucial for maintaining cellular homeostasis by eliminating damaged mitochondria [9, 10]. A diversity of mitophagy mechanisms has been elucidated, with distinct stimuli inducing mitophagy through varying pathways. Broadly, mitophagy can be categorized into PRKN (parkin RBR E3 ubiquitin protein ligase)-dependent and PRKN-independent pathways [11]. In cancer, similar to non-selective autophagy, it is widely acknowledged that mitophagy plays a dual role: it suppresses tumorigenesis but can also facilitate tumor metastasis, drug resistance, and the maintenance of cancer stem cell properties [12,13,14]. In Parkin-deficient mice, the growth and metastasis of melanoma are inhibited, suggesting that Parkin-dependent mitophagy plays a tumor-promoting role [15]. BNIP3-dependent mitophagy can facilitate cell invasion and resistance to apoptosis by remodeling the cytoskeleton, indicating that mitophagy can promote tumor metastasis [16]. Consequently, due to its involvement in tumor adaptability, the targeting of the mitophagy pathway may present itself as a novel and highly efficacious strategy for addressing metastatic OSA. Nevertheless, pertinent research in this area remains scarce.

Traditional bulk-seq, which rely on bulk cell populations, often lack the necessary resolution to discern specific cell types and fail to capture the intricate intratumoral heterogeneity present in OSA [17]. In contrast, single-cell RNA sequencing (scRNA-seq) represents a cutting-edge approach, enabling the detection of cellular heterogeneity at the single-cell level and facilitating the exploration of intercellular communication [17, 18]. The Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) system, an adaptive defense mechanism in numerous bacteria and archaea, has been repurposed as a powerful tool for precise genome editing [19]. Prior studies have demonstrated the efficacy of CRISPR-Cas9 in identifying genes that are critical for the proliferation and survival of cancer cells, thereby constructing cancer dependency map (DepMap) that unveils new methodologies for pinpointing genes with potential as predictive biomarkers and therapeutic targets in oncology [20,21,22]. Consequently, integrating scRNA-seq, CRISPR-based genetic screens, and machine learning approaches to identify pivotal biomarkers of mitophagy in OSA presents a groundbreaking avenue for advancing our understanding and treatment of this malignancy.

In this investigation, we have amalgamated diverse omics analyses, encompassing scRNA-seq, bulk-seq (from both tissues and cell lines), and CRISPR-Cas9, to comprehensively explicate the heterogeneity of OSA under varying mitophagy states and its interplay with the tumor microenvironment. Moreover, we have devised an RNA sequencing-based mitophagy scoring framework termed MIP_score. By merging this scoring framework with CRISPR-Cas9, we have successfully pinpointed RPS27A as an indispensable gene significantly associated with OSA metastasis.

Materials and methods

Data collection and preprocessing

To assess the prognostic impact of MRGs in OSA patients, we downloaded the scRNA-seq sequencing dataset (GSE152048) from the GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152048). Furthermore, we acquired RNA-seq data and clinical follow-up details, including age, sex, primary tumor location, metastatic state at diagnosis, survival duration, and status, for 88 OSA patients from TCGA. The immunohistochemical (IHC) results for glioma and melanoma are sourced from the HPA database (https://www.proteinatlas.org/). The 29 MRGs are provided in Additional File 1: Table S1.

scRNA-seq analysis

To analyze the scRNA-seq dataset, we loaded it into R using the Seurat R package (version 4.0.3) and utilized the functions Read10X and CreateSeuratObject to generate a Seurat object [23]. We filtered out single cells from subsequent analysis based on the following criteria: (1) cells containing less than 300 expressed genes and (2) cells with a mitochondrial rate greater than 10%. The data from the expression table were normalized and regressed according to the UMI value of each sample to obtain the scaled data. To plot high-dimensional data, we utilized two dimensions reduction algorithms. First, we used the harmony function to integrate the cells and remove batch effects. Then, we performed UMAP to reduce dimensionality and visualize the cell clusters. Dimension reduction is not only useful for visualizing high-dimensional data but also for reducing random noise in the data, allowing for better distinction of any biological signals in subsequent analyses and avoiding the curse of dimensionality. To identify marker genes, we used the “FindMarkers” function with the Wilcox rank-sum test algorithm. We retained marker genes that met the following criteria: logFC greater than 0.25, P-value less than 0.05, and min.pct greater than 0.1. To calculate the gene set activity, we compared the expression of mitophagy-related genes against a set of control genes with expression levels implemented in Seurat via the “AddModuleScore” function, and divided the cells into three groups (high, middle, low) [24].

Cell communication analysis

Intercellular communication analysis was integrated using the CellChat R package with default parameters [25]. The analysis focused on the communication between OSA stromal cells and osteoblastic cells, classified based on varying levels of mitophagy-related gene expression. The resulting cell-cell communication network was visualized by quantifying communication pathways, assessing information flow, and identifying specific pathways within the relevant cell types. Subsequently, specific communication pathways of interest were selected for further visualization.

Pseudotime trajectory analysis

We employed the Monocle algorithm (v2.18.0) to assess the dynamic states of osteoblastic cells through pseudotime trajectory analysis [26, 27]. Monocle utilized an unsupervised algorithm to integrate single-cell whole-transcriptome profiles, enabling the identification of a developmental trajectory that represents the progression of individual cells during differentiation. The “reduce Dimension” approach was utilized to calculate the CellDataSet object, and the “DDRTree” function was applied to reduce dimensionality. Subsequently, we utilized the “orderCells” function to extract the two most informative features. We then used these features as coordinate axes to visualize the trajectory. Additionally, we performed branched expression analysis modeling (BEAM) to identify genes with branch-dependent expression patterns in order to elucidate the mechanisms underlying fate decisions.

Construction risk scoring model based on TCGA and GEO data

We performed multivariate Cox analysis of DEGs via using the R package “survival” to identify MRG with prognostic value. We constructed a risk model based on TCGA-train cohort via the expression of MRG and corresponding coefficient was defined as the MIP_score = ∑in (coefficient × expression). Based on the MIP_score we utilized clipping values for receiver operating characteristic (ROC) curves to define the threshold values and divided the patients into two different groups using the median MIP_score. We created the survival curves by Kaplan–Meier (KM) analysis and used the “survival” R package by the log rank test to accessed the accuracy of the prediction, and constructed ROC curves for the risk scores using the “timeROC” R package. TCGA-test cohort and GSE21257 data continued to repeat the above steps to further confirm the prognostic accuracy of this risk score model. Additionally, to enhance prediction accuracy, we integrated gender, age, metastasis status, and risk scores of OSA patients from the TCGA database into the model. The nomogram was constructed to intuitively predict OSA patient survival, enabling clinical application of the risk-scoring model. Moreover, the nomogram’s stability was rigorously validated.

CRISPR-Cas9 assay was performed using CERES analysis

The DepMap database utilizes CERES, a computational method that combines gene copy number and sgRNA dropout data, to evaluate cellular dependency on genes, particularly in cancer cells. Through CRISPR-Cas9 screening, which allows for high-throughput gene knockouts, researchers identify essential genes for tumor proliferation and survival, providing potential therapeutic targets. The CERES scores reflect the influence of essential and nonessential genes across different cancer cell lines. Essential genes, especially those with negative CERES scores, indicate a higher likelihood of being critical for specific cancer cell line survival. In our study, we focused on 989 cell lines encompassing 32 types of cancer, analyzing the CERES scores of the RPS27A, TOMM20, and UBB genes within these cell lines. We utilized a threshold of < -1 to identify genes that exhibit significant dependency in cancer cells [28].

Drug sensitivity evaluation

The prophetic R package is utilized to predict patient responses to chemotherapy drugs [29]. By leveraging gene expression and drug sensitivity data, it calculates each patient’s drug response within high-risk and low-risk groups in the MIP_score. A lower IC50 score indicates higher drug sensitivity and, consequently, a more favorable therapeutic outcome.

Pan-cancer prognosis

LASSO algorithm and Cox proportional hazard regression for the MRG implemented a model to assess the association of each gene expression of mitophagy signature with the DFI, DSS, OS, and PFI. We subsequently obtained each patient mitophagy score based on the Cox proportional hazard regression model of MRG. To explore the association between MRG and malignant characteristics, we utilized the z-score algorithm to quantify the tumor’s involvement in RPS27A promotion, mitophagy, apoptosis, and hypoxia. This algorithm integrates gene expression of specific features to reflect pathway activity.

Cell culture

In accordance with previous studies, apart from the normal human osteoblast cell line hFOB1.19, the human OSA cell lines (143B, SJSA-1, MNNG, MG63, and Saos-2) were all cultured in DMEM (Gibco, USA) supplemented with 10% fetal bovine serum (BI, Israel) and 1% penicillin G-streptomycin (Gibco, USA) [30]. The hFOB1.19 cells were maintained in DMEM/F-12 (Gibco, USA) enriched with 10% fetal bovine serum and 0.3 mg/ml G418 (Invitrogen, Carlsbad, CA, USA).

Cell transfection

Cell transfection was performed using an established protocol detailed in the literature [31]. For shRNA transfection, EndoFectin™ Max (GeneCopoeia, Guangzhou, China) was utilized according to the manufacturer’s guidelines. The transfection reagent was added to Opti-MEM (Gibco, Carlsbad, California, United States) containing shRNA and allowed to incubate for 5 min. The mixture was then left to stand for 15 min. After 8 h, the culture medium was replaced, and the cells were harvested 24 h after shRNA transfection.

Real-time reverse transcription polymerase chain reaction (RT-PCR)

Using the RNAiso Plus kit (TaKaRa, Japan), RNA was extracted and purified employing the phenol-chloroform method. Subsequently, the Evo M-MLV RT Kit (Accurate Biotechnology, China) was used to reverse transcribe the RNA into cDNA. Finally, RT-PCR was conducted with the SYBR Green Premix Pro Taq HS qPCR kit (Accurate Biology, China) to measure RNA expression levels. The primers utilized in this study are listed in Additional File 1: Table S2.

Cell proliferation assay

Cell viability was assessed as per established protocols. We employed a CCK8 assay (CCK-8, GLPBIO, USA) for this purpose. Briefly, approximately 2 × 10³ OSA cells were seeded into individual wells of 96-well plates. Following the designated treatment, 10 µl of CCK8 reagent was added to each well. After a 2-hour incubation period at 37℃, absorbance readings were taken using a microplate reader at 450 nm.

Immunofluorescence

The mitochondrial transmembrane potential changes in cells were assessed using the JC-1 mitochondrial membrane potential detection kit (Biyuntian, China). After washing the cells with PBS, JC-1 was incubated with the cells for 20 min, and the fluorescent signals were observed under the AX NIS-Elements 5.4 (Nikon, Japan).

Fix, permeabilize, and encapsulate the cultivated OSA cells. After incubating with the primary antibody overnight at 4 °C, incubate with the secondary antibody for 1 h at 37 °C, protecting from light. Employ 0.1 g/ml of DAPI (Beyotime) to stain the cell nuclei and 25 nM Mito-Tracker Red CMXRos (MTRC, Beyotime, China) to label the mitochondria. Following rinsing, cellular imaging was conducted utilizing a confocal microscope.

Flow cytometry

Flow cytometry was employed to assess apoptosis and cell cycle progression. Logarithmically growing cells were trypsinized without EDTA and subjected to three washes with PBS. Apoptosis was quantified according to the manufacturer’s instructions, utilizing the Annexin V-APC Apoptosis Detection Kit (BestBio, China). Each experiment was repeated at least three times.

Transwell assays

As mentioned in our previous research [30], we utilized Transwell migration chambers or Matrigel-coated Transwell invasion chambers (BD Biosciences, USA). The transfected cells were seeded into the upper chambers containing 200 µl of serum-free DMEM, while the lower chambers were filled with 600 µl of DMEM supplemented with 10% FBS. Cells that had migrated or invaded the lower surface were photographed using an inverted light microscope.

Xenograft tumors

Five-week-old male BALB/c athymic nude mice were obtained from the Experimental Animal Center of Southern Medical University in Guangzhou, China, which were housed in designated facilities. Subcutaneous inoculations of 143B cells (5 × 10⁶ cells) suspended in 100µL of PBS were administered to each group. The cells were transfected with either Con-shRNA or RPS27A-shRNA. Weekly assessments of tumor progression and volume were conducted, culminating in live imaging of the mice at the study’s conclusion.

Intravenous injections of 143B cells containing Con-shRNA or RPS27A-shRNA (2 × 10⁶ cells) were administered into the tail veins of the subjects to establish a model of pulmonary tumor metastasis. Photographic documentation of the tumors was undertaken prior to the conclusion of the experimental procedures. All methods pertaining to animals have received formal approval (ID: LAEC-2022–043) from the Ethics Committee of Zhujiang Hospital, affiliated with Southern Medical University, and are conducted in accordance with the ethical guidelines established therein.

Immunohistochemistry (IHC)

As previously delineated, the determination of target protein expression in tissue samples derived from xenograft models was conducted via IHC analysis, as elucidated in reference. Incubate the tissue with primary antibodies against PCNA (1:300), E-cadherin (1:200), and Vimentin (1:300) simultaneously. Utilize an optical microscope (Olympus, Japan) for the observation and imaging of tissues.

Statistical analysis

To assess inter-group differences, we employed t-tests or one-way analysis of variance (ANOVA). Survival curves were calculated using the Kaplan-Meier method, and survival information was compared through log-rank tests. Statistical evaluations were conducted using R version 4.2.2 and SPSS version 25.0. Statistical significance was considered when the p-value was less than 0.05, and the results were presented as mean values and corresponding standard deviations (SD).

Result

Heterogeneity in OSA metastasis revealed by scRNA-seq

Figure 1 provides an overview of the study’s experimental design and data analysis pipeline. To investigate the cellular heterogeneity resulting from the expression of mitophagy-related genes, which are associated with metastasis and poor prognosis in OSA, we conducted scRNA-seq experiments on OSA tissues. The scRNA-seq data can be accessed openly via GEO. The table detailing the clinical characteristics of OSA patients reveals that samples BC2, BC3, BC5, BC6, BC16, BC21, and BC22 were diagnosed as primary cancer, while samples BC10 and BC17 exhibited metastatic cancer. For each donor, bone cells are subjected to stringent quality control protocols, encompassing tissues from seven primary OSA lesions and two from pulmonary metastatic OSA lesions (Fig. 2A). To mitigate potential batch effects from the nine donors, we employed PCA and UMAP algorithms within Seurat. The resultant plot in Fig. 2B demonstrates a significant separation within the merged data when categorized by dataset in the first two principal components. This suggests that cells are primarily clustered by cell type, which accounts for the majority of variance observed in the data.

Fig. 1

The overview of this study. Schematic diagram of the experimental design and analysis

Fig. 2

Single-cell transcriptome analysis of metastatic OSA. (A) UMAP plot illustrates cells from metastatic and primary OSA following stringent quality control measures. (B) UMPA plot of 12 cell clusters (B cells, chondroblastic cells, endothelial cells, fibroblasts, myeloid cells, myoblasts, NKT cells, Osteoblastic cells, Osteoblast Proliferating Cells, Osteoclasts, Pericytes). (C) Expression of diagnostic maker gene in each cell cluster via bubble charts. (D) Proportional distribution of cell types in metastatic and primary OSA. (E) UMAP plots illustrating mitophagy pathway activity distribution across donors

Using a set of marker genes identified from prior investigations, we performed comprehensive cell type annotation and characterized 11 distinct categories, including osteoclasts (ACP5, CTSK, and MMP9), B cells (MS4A1 and JCHAIN), chondroblastic cells (ACAN, SOX9, and PTH1R), endothelial cells (PECAM1 and VWF), fibroblasts (DCN and COL1A1), myeloid cells (FCGR3A, CD74, and CD14), myoblasts (MYL1 and MYLPF), NKT cells (NKG7 and GNLY), osteoblastic cells (RUNX2, CDH11, and IBSP), osteoblast proliferating cells (MKI67, TOP2A, and PCNA), and pericytes (RGS5 and ACTA2) (Fig. 2F). To validate the identities of the OSA clusters, we analyzed a list of known subtype markers. Dot size indicates the cell percentage expressing the gene, while the color scale represents the average normalized expression level within each cluster. (Fig. 2C-F).

Furthermore, our analysis revealed variations in the proportion of cell types between primary and metastatic OSA tissue. Specifically, the proportion of osteoblastic cells and NKT cells was lower in primary OSA compared to metastasis, while the proportion of osteoclasts and myeloid cells was higher in primary OSA tissue (Fig. 2D). These findings emphasize the demand for further research on the role of these cell types in OSA progression and metastasis, as they suggest potential differences in the cellular composition at different sites.

Dynamic modulations of mitophagy pathway activity during metastasis of OSA

Osteoblastic and chondroblastic OSA are identified as the predominant subtypes of conventional OSA in clinical practice. Significant heterogeneity in mitophagy pathway activity is present among the examined cell types at varying OSA stages. Chondroblastic and osteoblastic OSA cells, in particular, exhibit the obvious contrast in mitophagy pathway activity, with the latter demonstrating more pronounced activity in correlation with malignancy (Fig. 3A). Consequently, this study focused predominantly on osteoblastic cells. Their transcriptomic data were reanalyzed and depicted via UMAP visualization (Fig. 3B). The mitophagy pathway activity in osteoblastic cells’ clusters was quantified using the “AddModuleScore” function and is illustrated in Fig. 3C. osteoblastic_cluster 3 showed the highest mitophagy activity, whereas osteoblastic_cluster 11 had the least. Notably, osteoblastic_cluster 2, osteoblastic_cluster 3, and osteoblastic_cluster 4 showed the most substantial differences. Specifically, osteoblastic_cluster 3 demonstrated lower expression in primary tumors but higher expression in metastasis, whereas osteoblastic_cluster 2 and osteoblastic_cluster 4 exhibited the opposite pattern (Fig. 3D). Additionally, a distinct distribution in the cell types, categorized by different mitophagy activity expressions, was evident between the metastatic and primary groups (Fig. 3E). These findings suggest a potential correlation between the mitophagy and OSA, indicating that higher mitophagy is associated with an increased propensity for OSA metastasis.

Fig. 3

Enrichment analysis of mitophagy in OSA. (A) The violin plot displays the expression of mitophagy pathway score in subtype cells of the metastasis and primary groups. (B) The identities of 15 clusters of osteoblastic cells visualized by UMAP. (C) Percentage of cell population in 15 osteoblastic cell subtypes present in each of the primary and metastasis tissue. (D) The violin plot illustrates the expression of mitophagy pathway score in osteoblastic cell subtypes. (E) UMAP plots depict the mitophagy pathway score and distribution of subtype cells in the metastatic and primary groups, with darker colors indicating higher mitophagy pathway score

Mitophagy-dependent modulates intercellular communication in osteoblastic OSA cells

We conducted cell-communication analysis between different types using the “CellChat” R package to reveal the influence of mitophagy on the microenvironment of the OSA. Osteoblastic cells were classified by mitophagy activity using quartiles: the top 25% as high, the bottom 25% as low, and the middle 50% as medium. We then aggregated the number of interactions and their respective weights to construct the intercellular communication network. Figure 4A and B illustrate significant disparities in the intensity of incoming and outgoing signals and the number of their interactions, indicating that the activity of mitophagy in distinct osteoblasts OSA cell can regulate biological functions through intercellular communication. Notably, osteoblastic cells with low mitophagy exhibited higher strength in both incoming and outgoing signals compared to those with high mitophagy (Fig. 4C-D). Drawing upon the ligand-receptor (LR) paradigm, we delve into how mitochondrial autophagy achieves intercellular communication by modulating intercellular LR interactions. Osteoblastic OSA cells with low mitophagy activity enhanced cellular communication with pericytes, NKT cells, myeloid cells, and fibroblasts via SPP1-(ITGAV + ITGB1), SPP1-(ITGAV + ITGB5), and SPP1-CD44 interactions; they also demonstrate a more active cell communication capacity with fibroblasts through POSTN-(ITGAV + ITGB5) interactions. Notably, osteoblastic OSA cells with high mitophagy activity display stronger intercellular communication with endothelial cells through VEGFA-VEGFR1, VEGFA-VEGFR1R2, VEGFA-VEGFR2, and VEGFB-VEGFR1 connections (Fig. 4E). These findings highlight how distinct mitophagy directly impacts cell communication between osteoblastic OSA cells and other cells within the tumor microenvironment, thus influencing tumor metastasis.

Fig. 4

The illustration depicts the intercellular communication network between OSA cells exhibiting varying levels of mitophagy activity and other co-localized cell types. (A-B) Interactions between mitophagy pathway activity subtypes of osteoblastic cells and other cell types, depicted in a diagram. The line thickness represents the strength of their interaction weight/intensity. (C) Point plots illustrating outgoing and incoming interaction signal pathways across cell types. (D) The heatmap illustrates the intensity of intercellular interactions among 13 distinct cell types. (E) Upregulated receptor-ligand interaction networks between osteoblastic cells of different mitophagy pathway activity types and other cell types. The horizontal axis represents ligand-receptor cell pairs, while the vertical axis represents different signaling pairs

To deeper elucidate the mechanisms of OSA metastasis, we analyzed signal distribution in certain subgroups and highlighted significant signaling molecules—PTN, PERIOSTIN, MIF, VISFATIN, VEGF—that play roles in apoptosis, cytokine activity, extracellular matrix and receptor binding, growth factor activity, and enzyme inhibition. Our research specifically indicates that myoblasts are the most active in the PTN signaling pathway, followed by OSA cells. The osteoblastic OSA cells, with varying levels of mitophagy, function as mediator in distinct manners (Fig. 5A-B). Notably, osteoblastic OSA cells with low mitophagy in the PERIOSTIN signaling pathway interact more robustly with tumor-associated fibroblasts (Fig. 5C-D). Within the MK signaling pathway, osteoblastic OSA cells exhibiting reduced mitophagy are more likely to act as receptors (Fig. 5E-F). Conversely, osteoblastic OSA cells with high mitophagy exhibit stronger interactions with myeloid cells, NKT cells, and B cells in the MIF pathway, leaning towards a sender role (Fig. 5G-H). Compared to osteoblastic OSA cells with low mitophagy, those with high mitophagy demonstrate more prominent cell communication with endothelial cells in the VISFATIN and VEGF pathways (Fig. 5I-L). Therefore, osteoblastic OSA cells with different levels of mitophagy activity can influence OSA metastasis by modulating the tumor microenvironment.

Fig. 5

The heatmap displays the relative likelihood of osteoblastic OSA cells, categorized by mitophagy types, assuming four distinct roles (sender, receiver, mediator, and influencer) within the signaling pathways. Color intensity corresponds to the magnitude of cellular impact. The Circos diagram depicts interrelationships between osteoblastic cells of different mitophagy types and other cell types across various signaling pathways

Pseudotime analysis revealed modulation of osteoblastic OSA cells in response to mitophagy alterations

Subsequently, pseudotime analysis was employed to assess whether mitophagy influenced the developmental trajectory of osteoblastic OSA cells. Figure 6A illustrates osteoblastic OSA cells labeled with distinct colors. The pseudotime analysis reveals a temporal progression of cell differentiation, indicated by a gradient of blue shades, with darker tones representing earlier differentiation stages. The result in Fig. 6B reveal that over time, osteoblastic OSA cells differentiate from right to left, with the darker shades of blue indicating the earlier differentiated cell state 1, while cell states 2 and 3 represent later stages of osteoblastic cells differentiation. Figure 6C illustrates the distribution of primary OSA cells and metastasis OSA cells during the differentiation process. Figure 6D showcases the distribution of different subgroups of osteoblastic along the cellular differentiation trajectory. By integrating the differences in mitophagy activity, we discovered that osteoblastic_cluster 3, characterized by the highest metabolic activity, is predominantly located in the later stages of the cellular trajectory, corresponding to metastasis OSA tissue. On the other hand, the osteoblastic_cluster 11, exhibiting the lowest metabolic activity, is primarily distributed in the early stages of the cellular trajectory, indicating primary OSA tissue. These findings strongly support the close association between mitophagy activity and the metastasis of OSA.

Fig. 6

Pseudotime trajectory analysis of clusters of the osteoblastic OSA cells. (A) Pseudotime trajectory differentiation plot of osteoblastic OSA cells. (B) Pseudotime trajectory illustrating distinct differentiation states of osteoblastic OSA cells. (C) Osteoblastic OSA cell differentiation in metastatic and primary OSA. (D) Pseudotime trajectory analysis plot depicting differentiation based on distinct clusters of osteoblastic OSA cells

Construction of the MIP_score model

Given the critical interplay between mitophagy and OSA, there is an urgent need to identify precise biomarkers. Utilizing multivariable Cox regression analysis, we have pinpointed three pivotal mitophagy-associated prognostic genes: RPS27A, TOMM20, and UBB. Every patient acquired a MIP_score according to the expression of survival-associated MRGs and their coefficients (Fig. 8A-B). The MIP_score was calculated using the following equation: MIP_score = (RPS27A × 1.50103539587517 + TOMM20 × 1.10434214530872 − UBB×1.19608412496501). Based on the median MIP_score from the TCGA-train dataset, patients were stratified into low-risk and high-risk cohorts. Kaplan-Meier survival analysis revealed that OSA patients in the low-risk group exhibited significantly improved overall survival (Fig. 7A). Time-dependent ROC analysis further underscored the prognostic robustness of the MIP_score, with AUC values of 0.974, 0.848, and 0.808 at 1, 3, and 5 years, respectively (Fig. 7B). Figure 7C depicts the distribution of risk scores, survival status, and the expression heatmap of three candidate genes in OSA patients. Moreover, the TCGA-OSA dataset was divided into a training set and an internal validation set in a 1:1 ratio, with GSE21257 employed as an external validation cohort. Our findings demonstrate that the MIP_score accurately predicts OSA prognosis across both internal and external datasets(Fig. 7D–I).

Fig. 7

The prognostic model of mitophagy-associated genes was rigorously evaluated and validated in both the training and validation cohorts. (A–C) The distribution of Kaplan-Meier survival curves, ROC curves, and risk scores in the training cohort. (D–I) The distribution of Kaplan-Meier survival curves, ROC curves, and risk scores in the validation cohort

Fig. 8

Establishing a clinical nomogram for MRGs. (A-B) Regression coefficients of 3 genes obtained in Cox regression. (C) Developing a prognostic nomogram to predict the 3-year and 5-year survival rates of OSA patients. (D) Conducting decision curve analysis to compare nomogram, MIP_score, and various clinical features. (E-F) Assessing the ROC curves of 3-year and 5-year nomogram for the predictive model. (G-I) Differences in pathway activities between the high- and low-risk groups scored by GSVA and GSEA

Furthermore, to facilitate clinicians in swiftly understanding the classification of OSA patients, we have developed a web server (https://mip-score.shinyapps.io/labtan/) (Figure S1). This server allows users to input the expression levels of RPS27A, UBB, and TOMM20, automatically generating survival graphs and predicting survival time.

Evaluation of the three hub genes involved in MRGs

Based on the results of the multivariable Cox analysis, we constructed a prognostic nomogram model within the TCGA cohort. The findings confirm the independence of the MIP_score as a significant prognostic factor (Fig. 8C). Decision curve analysis (DCA) demonstrated the notable clinical utility and potential of this MIP_score (Fig. 8D). Furthermore, ROC analysis established the accuracy of the MIP_score in predicting clinical outcomes. The AUC values for the MIP_score at 3 years and 5 years were 0.738 and 0.755, respectively, indicating excellent prognostic performance (Fig. 8E-F). These findings underscore the significance of investigating three hub MRGs (RPS27A, TOMM20, and UBB) and their MIP_score.

Subsequently, we employed GSVA and GSEA to identify the pathways through which MIP_score modulates the malignant behavior of OSA by affecting. Our findings indicate that the low-risk group is associated with cellular apoptosis, oxidative damage, inflammatory response, JAK_STAT signaling pathway, cytokines, and PD1 tumor immunotherapy-related processes, while the high-risk group exhibits enrichment in processes related to urea metabolism, TCA metabolism, steroid metabolism, cell adhesion, and cell differentiation (Fig. 8G-I).

To delve deeper into the immunological landscape mediated by MIP_score, we applied multiple deconvolution techniques including TIMER, MCPCOUNTER, QUANTISEQ, CIBERSORT-ABS, XCELL, EPIC, and CIBERSORT. We have observed that the MIP_score is associated with reduced infiltration of most immune cells, including macrophages, CD8 + T cells, CD4 + T cells, monocytes, NK cells, mast cells, and neutrophils. These findings also substantiate that the MIP_score serves as an independent risk factor, exerting its influence on immune infiltration to modulate the biological changes in OSA (Fig. 9A).

Fig. 9

(A) The abundance of immune cell infiltration between the high- and low-risk groups in the OSA cohorts. (B-L) The MIP_score constructed using three mitophagy genes was compared with previous studies, including the C-index, ROC curves, and survival curves

To further assess the performance of the MIP_score in predicting overall survival, we conducted a comparative analysis with other established prognostic models. Our research findings indicate that the MIP_score consistently outperforms other prognostic features across all datasets. Consequently, our ROC and C-index analyses affirm the enduring precision and stability of MIP_score in prognosticating clinical outcomes in OSA, thereby emphasizing its significant potential for clinical implementation (Fig. 9B-L).

Chemotherapy response prediction and mitophagy related genes signature

To explore the potential relationship between our constructed MIP_score and drug sensitivity, we employed the pRRophetic package to predict the IC50 for various chemotherapeutic agents among OSA patient risk groups (TCGA and GSE21257). This analysis showed that drugs such as Bexarotene, JNK.Inhibitor, KIN001.135, GSK269962A, Parthenolide, Shikonin, and Vorinostat exhibited significantly higher efficacy in low-score patients, indicating increased drug resistance in the high-score group due to the need for higher drug concentrations to achieve similar inhibitory effects (Figure S2). Furthermore, IC50 analysis demonstrated that patients with a high-mitophagy score were less responsive to standard anti-OSA chemotherapy compared to those with a low-mitophagy score, suggesting a greater potential clinical benefit for the latter. These findings suggest a potential correlation between our model genes and drug sensitivity, offering invaluable insights for personalized therapeutic strategies in OSA.

Enhanced mitophagy in OSA cell lines with higher metastatic potential

To further investigate the role of mitophagy in the metastasis of OSA, we selected the highly metastatic cell line 143B and the low-metastatic cell line Mg63 to evaluate mitophagy-related markers. The JC-1 assay revealed a decrease in mitochondrial membrane potential in the high-metastatic 143B cells (Fig. 10A). Furthermore, the co-staining of lysosomes and mitochondria indicated a higher level of mitophagy in 143B than in Mg63 (Fig. 10B). These findings affirm a positive correlation between increased mitophagy and enhanced metastatic potential in OSA.

Fig. 10

The relationship between the metastatic ability of cell lines and mitophagy. A The mitochondrial membrane potential was assessed using JC-1 to examine the differences between cell lines with high and low metastatic capabilities. B Lysosomal labeling with Lyso-Tracker followed by co-staining with Mito-Tracker was conducted to visualize alterations in mitochondrial and lysosomal structures within OSA cells

Inhibiting RPS27A expression suppresses metastasis by reducing mitophagy in OSA cells

Based on the results obtained from CRISPR-Cas9, we identified three mitophagy related genes (RPS27A, TOMM20, and UBB) as essential for the survival of nine OSA cell lines. The results revealed that RPS27A is the most significant among these genes, exerting a notable impact on eight cell lines, including OS252, SJSA1, 143B, SAOS2, U2OS, G-292, and HSOS1. Consequently, RPS27A was selected as the focal point for subsequent experiments (Fig. 11A-C). We utilized RT-PCR analysis to verify the differential mRNA expression of RPS27A between normal osteoblasts and OSA cell lines. As illustrated in Fig. 11D, the mRNA expression levels of RPS27A in OSA cell lines were significantly higher than in normal osteoblasts, with the highest expression observed in the 143B cell line. To examine the impact of RPS27A on OSA, we employed RPS27A-specific shRNA to downregulate RPS27A expression in 143B cells. As shown in Fig. 11E, RT-PCR analysis indicated that all shRPS27A constructs effectively downregulated RPS27A expression. We conducted CCK8 cell proliferation assays, Transwell invasion assays, and flow cytometry to evaluate the impact of RPS27A silencing on OSA cell proliferation, invasiveness, metastasis, and apoptosis (Fig. 11F-K). The results indicated that, compared to the corresponding controls, transfection with shRPS27A led to decreased proliferation, reduced invasion and migration abilities, and increased apoptosis of OSA cells.

Fig. 11

The role of RPS27A in OSA cells is investigated. A-C The CERES scores for UBB, TOMM20, and RPS27A across nine OSA cell lines. D Quantitative analysis of RPS27A expression in five OSA cell lines via RT-PCR. E Assessment of RPS27A expression following sh-RPS27A treatment using RT-PCR. F CCK-8 assays demonstrate that RPS27A knockdown suppresses the proliferation of OSA cells. G-I Transwell assays illustrate the effect of RPS27A knockdown on the metastatic capacity of OSA cells. J-K Flow cytometry-based apoptosis assays evaluate the impact of RPS27A knockdown on OSA cell apoptosis

Additionally, we investigated the interplay between RPS27A and mitophagy in OSA cells. Autophagosome formation, assessed via endogenous LC3B staining, revealed that the shRPS27A significantly downregulated LC3B expression (Fig. 12A). Confocal microscopy was utilized to investigate the alterations in mitophagy. Immunofluorescence (IF) co-localization of mitochondria and lysosomes confirmed that the reduction in RPS27A expression attenuated mitophagy (Fig. 12B). These results highlight the suppressive role of RPS27A-specific shRNA in mitophagy within OSA cells.

Fig. 12

Reduction of RPS27A is implicated in the regulation of mitophagy within OSA cells. A Immunofluorescence staining was utilized to assess differential expression of LC3B subsequent to RPS27A knockdown. B Lyso-Tracker (green) and Mito-Tracker (red) were co-stained to observe the changes of mitophagy in OSA cells

The influence of RPS27A on OSA was further corroborated using a xenograft mouse model. 143B cells were transduced with lentiviruses encoding LV-RPS27A or a control vector (LV-vector). Puromycin selection was employed to isolate lentivirus-infected 143B cells in vitro. These selected cells were subsequently injected subcutaneously into nude mice to establish xenografts. Our findings revealed that RPS27A downregulation significantly inhibited OSA growth in these mice (Fig. 13A-D). Immunohistochemical (IHC) analysis of the harvested xenografts was performed to evaluate the expression of PCNA, Vimentin, and E-cadherin. Pathological examination indicated that, compared to the control group, shRPS27A-treated mice displayed increased in E-cadherin-positive cells and a marked reduction in active PCNA and Vimentin-positive cells (Fig. 13E). To further elucidate the role of RPS27A in OSA cell metastasis in vivo, a mouse lung metastasis model was established. Bioluminescent imaging revealed a significant reduction in lung metastasis following RPS27A knockdown in OSA cells (Fig. 13F-G). These results indicate that RPS27A plays a crucial role in regulating the malignant behavior of OSA cells by influencing mitophagy.

Fig. 13

In vivo experimental validation of RPS27A. A-B Employing an in vivo imaging system to evaluate the effects of RPS27A on cellular proliferation. C-D Investigating the impact of RPS27A on tumor growth and mass in nude mice models. E Assessing the expression levels of PCNA, E-cadherin, and Vimentin proteins across various groups via IHC. F-G Utilizing an in vivo imaging system to analyze the influence of RPS27A on cellular metastasis

Prognostic role of MIP_score and RPS27A in a pan-cancer cohort

Given the prognostic significance of MIP_score in osteosarcoma, we believe it possesses broader potential applications in pan-cancer research. Prognostic analyses of MIP_score across 32 cancer types reveal a positive correlation between MIP_score and poor prognosis in LGG (lower grade gliomas), UVM (uveal melanomas), and LUSC (lung squamous cell carcinomas) (Figure S3). Moreover, CRISPR-Cas9 knockout results from 989 cell lines across these cancer types indicate that the loss of RPS27A affects the survival of 96.7% of cell lines (956 lines), while UBB and TOMM20 impact 9.1% (89 lines) and 3.79% (34 lines) of cell lines, respectively (Table S3). Consequently, RPS27A emerges as a biomarker worthy of further exploration in pan-cancer studies. CRISPR-cas9 results in LGG and UVM also suggest that RPS27A could serve as a biotherapeutic target (Figure S4 A-D). Immunohistochemical results from gliomas and melanomas reveal significantly elevated RPS27A expression in more aggressive gliomas and highly metastatic melanomas (Figure S4 E-F). Additionally, pan-cancer analyses demonstrate a positive correlation between RPS27A expression and mitochondrial autophagic activity, apoptosis, and hypoxia (Figure S5). Collectively, these findings support the use of the MIP_score as a prognostic indicator for both gliomas and melanomas, while also suggesting that RPS27A may serve as a valuable prognostic target. This underscores the importance of investigating mitochondrial autophagy in the study of these tumors.

Discussion

The limited understanding of OSA pathogenesis presents a challenge in enhancing clinical outcomes, notably due to the tendency of OSA to metastasize to the lungs, significantly influencing patient mortality and posing challenges for clinicians [32]. Metastasis and chemoresistance are principal causes of mortality in OSA patients, yet their underlying mechanisms are not fully elucidated [33]. Mitophagy, a selective form of autophagy, plays a pivotal role in cellular homeostasis by targeting and degrading damaged or dysfunctional mitochondria [11, 12, 34]. This process has been implicated in the metastatic progression of various malignancies, including breast cancer, colorectal cancer, and gastric cancer [35,36,37]. However, the relationship between mitophagy and metastasis in OSA remains inadequately explored. Consequently, it is imperative to investigate the molecular mechanisms underlying mitophagy in OSA metastasis and to identify reliable biomarkers for its effective monitoring and potential therapeutic targeting.

Through scRNA-seq analysis, we discerned marked heterogeneity in osteoblastic OSA cells within both primary and metastatic OSAs. Mitophagy exhibited the most significant variation in osteoblastic OSA, prompting us to focus on these cells. To investigate the regulation heterogeneity of mitophagy in OSA cells, we utilized the harmony integration algorithm to categorize the cells into distinct subtypes. Pseudotime analysis further suggested that elevated mitophagy correlates with a propensity for metastasis. These findings underscore the potential of inhibiting OSA metastasis by regulating mitophagy, identifying viable therapeutic targets, and designing corresponding pharmacological interventions.

During tumor metastasis, intercellular communication is vital for coordinating cellular behavior. Through mitophagy pathway scoring, we identified that OSA cells classified as “mitophagy-high” and “mitophagy-low” can dynamically modulate their interaction with the tumor microenvironment. Our analysis of cellular communication revealed that mitophagy pathway activity is intricately associated with endothelial cells, particularly through the VEGF-related pathways, where heightened mitophagy correlates with increased activity in VEGF signaling and receptor-ligand interactions. Moreover, we discovered that the mitophagy pathway is more active in the MIF pathway, facilitating tighter communication with NKT cells, B cells, and endothelial cells. MIF is a multifunctional cytokine that activates various signaling pathways, including ERK, MAPK, and Akt, promoting cancer cell proliferation and metastasis [38, 39]. VISFATIN primarily affects the production of NAD + phosphoribosyltransferase (NAMPT), enhancing tumor metastasis and angiogenesis through PI3K/Akt, MAPK, and c-Abl/STAT3 pathways in macrophages and endothelial cells [40]. Our findings indicate that OSA cells with high mitophagy activity exert a more pronounced influence on endothelial cells, a phenomenon rarely reported. These results collectively highlight the pivotal role of mitophagy in orchestrating the interactions between OSA cells and the tumor microenvironment, providing insights into potential therapeutic strategies to remodel immune responses and improve clinical outcomes in patients with high mitophagy activity.

Based on the aforementioned research, the relationship between mitotophagy and OSA metastasis has been delineated. We must now endeavor to elucidate precisely which mitophagy genes can impact the clinical treatment and prognosis of patients. By utilizing data from TCGA and GSE21257 for training and validation, we identified three pivotal genes among the mitophagy related genes—RPS27A, TOMM20, and UBB—and established a mitophagy risk model based on calculated coefficients. RPS27A, an established proto-oncogene, significantly influences the cell cycle, apoptosis, and proliferation, correlating with a poor prognosis in diverse cancer types. Li et al. depicted that the suppression of RPS27A enhances RPL11’s binding to MDM2, which inhibits MDM2-mediated ubiquitination and degradation of p53, pinpointing RPS27A as a contributor to tumorigenesis [41]. Furthermore, it has been reported that RPS27A is overexpressed in renal cancer, breast cancer, and colon cancer, indicating its potential as a clinical diagnostic and therapeutic target [42]. TOMM20, a receptor involved in the recognition and translocation of mitochondrial proteins from the cytosol into the mitochondria, is strongly associated with the malignancy level across various cancer types [43]. Recent research has demonstrated that high-grade human chondrosarcoma tumors exhibit elevated expression levels of TOMM20 in comparison to low-grade tumors [44]. This is attributed to TOMM20’s role in promoting proliferation, migration, apoptosis resistance, and chemotherapy resistance, thereby supporting our study. UBB, a gene belonging to the ubiquitin family, plays a role in virtually all cellular processes and is acknowledged for its heightened expression levels in numerous cancer types [45]. Experimental evidence suggests UBB’s involvement in the pathogenesis and progression of various disorders, including gynecological cancers, advanced renal cell carcinoma, and gastric cancer [45,46,47]. However, Wang J et al.’s research downregulation of UBB potentiates renal cell carcinoma cell proliferation, tumor burden, and angiogenesis [46]; thus, UBB could serve as a potential cancer inhibitor. The performance was rigorously evaluated against both external and in-house datasets, utilizing retrospectively corrected samples to ensure the robustness of the MIP_score across a spectrum of experimental conditions. More critically, utilizing the MIP_score enables us to pinpoint variations that might influence the chemotherapeutic response in OSA. Of particular interest is the potential of JNK inhibitors, tyrosine kinase inhibitors (Axitinib), p38 MAPK inhibitors (BIRB.0796), topoisomerase inhibitors (Camptothecin), ROCK inhibitors (GSK269962A), and Aurora kinase inhibitors (ZM.447439) to identify patient sensitivity through MIP_score analysis, thereby enhancing its clinical translational potential and implementation. Consequently, the identification of mitophagy genes as prospective therapeutic targets at both the transcriptomic and proteomic levels could bear significant clinical implications, facilitating the stratification of OSA patients using MIP_score for optimized treatment strategies and improved prognostic outcomes.

Additionally, employing CRISPR-Cas9 technology, we utilized the DepMap database to conduct whole-genome functional inactivation screenings across various cell lines. Our findings indicated that the mitophagy marker RPS27A significantly affects patient prognosis and is intricately involved in the initiation and progression of OSA, thus presenting a potential biomarker for the diagnosis and inhibition of OSA metastasis. Furthermore, our experimental results demonstrated that reduced expression of RPS27A can influence OSA metastasis by modulating mitophagy, which is consistent with the above conclusion.

While our research offers significant insights into the nexus between mitophagy and OSA metastasis, certain limitations must be acknowledged. Primarily, our analysis is based on retrospective data, highlighting the imperative for future studies to ascertain the clinical relevance of our findings. Additionally, the intricate and multifaceted mechanisms underlying OSA metastasis necessitate more extensive foundational and clinical research to elucidate the role of mitophagy-related genes in this process. Lastly, despite our preliminary use of CRISPR-Cas9 to discover and verify the role of RPS27A in OSA metastasis, there is a lack of phase III randomized controlled trials to confirm this gene’s role in decision-making for patients exempt from adjuvant chemotherapy (e.g., using Oncotype DX). Thus, high-quality, adequately followed, large-sample, multicenter randomized controlled trials are needed to substantiate our results.

Conclusion

During the course of conducting a series of experiments utilizing sc-RNA seq, CRISPR-Cas9 methodologies, and comprehensive omics data analysis, we meticulously delineated the genetic characteristics associated with mitophagy. Our investigations have unveiled the pivotal role of mitophagy in OSA metastasis prognosis and treatment efficacy. The development of the MIP_score holds promising potential for personalized therapeutic strategies in the future, including prognostic assessments, identification of therapeutic targets, and elucidation of chemotherapy resistance mechanisms, thus propelling forward the field of precision medicine in OSA clinical management.