Jiatong Li
Yang Bai2†
He Zhang
Ting Chen
Guanning Shang- 1Department of Orthopedics, Shengjing Hospital of China Medical University, Shenyang, China
- 2Department of Nursing, Shengjing Hospital of China Medical University, Shenyang, China
Introduction: Osteosarcoma is a common type of bone cancer characterized by a poor prognosis due to its metastatic nature. The tumor microenvironment (TME) plays a critical role in tumor metastasis and therapy response. Therefore, our study aims to explore the metastatic mechanism of osteosarcoma, potentially opening new avenues for cancer treatment.
Methods: In this study, we collected data from the GSE152048, GSE14359, and GSE49003 datasets. Differentially expressed genes (DEGs) were identified in osteosarcoma cases with primary and metastatic features using R software and the limma package. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to investigate metastasis-related genes. A protein–protein interaction (PPI) network was established using the STRING database to further analyze these metastasis-associated genes. The abundances of different cell types with a mixed cell population were estimated using the CIBERSORT approach. The scRNA-seq data were analyzed by the Seurat package in R software, and intercellular communications were elucidated using the CellChat R package.
Results: In this study, 92 DEGs related to metastasis were identified, including 41 upregulated and 51 downregulated genes in both the GSE14359 and GSE49003 datasets. Metastasis-associated pathways were identified, including those involving the cyclin-dependent protein kinase holoenzyme complex, transferase complex, transferring phosphorus-containing groups, SCF ubiquitin ligase complex, and the serine/threonine protein kinase complex. KEGG and PPI network analyses revealed 15 hub genes, including Skp2, KIF20A, CCNF, TROAP, PHB, CKS1B, MCM3, CCNA2, TRIP13, CENPM, Hsp90AB1, JUN, CKS2, TK1, and KIF4A. Skp2 has been known as an E3 ubiquitin ligase involved in osteosarcoma progression. The proportion of CD8+ T cells was found to be higher in metastatic osteosarcoma tissues, and high expression of PHB was associated with a favorable prognosis in osteosarcoma patients. Additionally, 23 cell clusters were classified into eight cell types, including chondrocytes, MSC, T cells, monocytes, tissue stem cells, neurons, endothelial cells, and macrophages. The 15 hub genes were expressed across various cell types, and interactions between different cell types were observed.
Conclusion: Our study reveals the intricate communication between tumor microenvironment components and tumor metastasis in osteosarcoma.
1 Introduction
Osteosarcoma is a common type of bone cancer that often occurs in the long bones of the legs and arms and is more prevalent among children and young adults. Symptoms of osteosarcoma include pain, swelling, and a noticeable lump in the affected bones (1). Treatment options include surgical excision, chemotherapy, radiation therapy, and targeted therapy (2, 3). Early detection and intervention are crucial for achieving better outcomes, which makes exploring the underlying mechanism of osteosarcoma vital (4). Although the exact cause of osteosarcoma is not fully understood, certain genetic and epigenetic factors may increase the risk of its development (5–7). It is essential to uncover the mechanisms of osteosarcoma development and reoccurrence in order to develop effective therapy for patients (8, 9).
Tumor metastasis is a key factor influencing therapeutic outcomes and prognosis in cancer patients (10, 11). Osteosarcoma is particularly characterized by its highly metastatic feature. Due to this tendency for metastasis, patients with osteosarcoma often experience poor treatment responses and worse prognosis. The 5-year survival rate for osteosarcoma is about 70%, but it drops to 20%–30% for patients with metastatic disease (12). The lungs are the most common site for osteosarcoma metastasis, although it can also spread to other bones or organs via the bloodstream or lymphatic system (13). Early detection of metastasis is critical for improving the chances of successful management in osteosarcoma patients.
The tumor microenvironment (TME) has gained significant attention because it plays a crucial role in tumor progression, metastasis, and treatment outcomes (14, 15). The TME refers to the surrounding microenvironment in which tumors exist, including both cellular and noncellular components such as noncancerous cells, blood vessels, and extracellular matrix (ECM) components. Noncancerous cells include immune cells, fibroblasts, and other cell types. Blood vessel formation supplies oxygen and nutrients to tumor cells and provides a pathway for tumor metastasis (16). The TME of osteosarcoma includes osteoblasts, osteocytes, osteoclasts, pericytes, endothelial cells, fibroblasts, mesenchymal stem cells, and various immune cells, such as lymphoid and myeloid cells, as well as the ECM (17, 18).
A better understanding of the TME in osteosarcoma, including the interactions between cancer cells and noncancerous cells and the mechanisms of metastasis, is crucial for developing effective tumor treatments (19). Tumor cells interact with their surrounding TME in complex ways, influencing their biological behavior and response to therapies. Immunotherapy often targets the TME to enhance the immune response against cancer by focusing on specific components within the TME. Immunotherapy includes treatments such as checkpoint inhibitors, cytokines, cancer vaccines, and chimeric antigen receptor (CAR) T-cell therapy (20). CAR T-cell therapy has demonstrated promising results in combating osteosarcoma (21). Targeting the TME has opened new avenues for cancer treatment. In this study, we used single-cell RNA sequencing (scRNA-seq) to reveal the communications between tumor microenvironment components and tumor metastasis in osteosarcoma.
2 Materials and methods
2.1 Data collection
The GSE152048 dataset from the GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152048) database was downloaded. The scRNA-seq data were obtained from the GSE152048 dataset. We selected two primary osteosarcoma samples (BC2, BC22) and two metastasis osteosarcoma samples (BC10, BC17) for further analysis. Additionally, the GSE14359 and GSE49003 datasets were downloaded from the GEO database. The GSE14359 dataset includes 10 primary osteosarcoma samples and 8 metastatic osteosarcoma samples, while the GSE49003 dataset includes 6 primary osteosarcoma samples and 6 metastatic osteosarcoma samples.
2.2 Differential and functional analysis
Differentially expressed genes (DEGs) were explored in osteosarcoma cases with primary and metastatic features using R software and the limma package. A cutoff of |log2 (fold change)| > 0.5 and p < 0.05 was used to determine the DEGs. Metastasis-related genes were dissected when DEGs were shared between the GSE14359 and GSE49003 datasets. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to examine metastasis-related genes. Several R packages, including “clusterProfiler”, “org.Hs.eg.db”, and “enrichplot” were used for functional analysis. A p-value < 0.05 was considered a significant enrichment.
2.3 Protein–protein interaction network construction
A protein–protein interaction network was established using the STRING database to analyze the metastasis-associated genes. The network was reconstructed and visualized using Cytoscape software. The cytohubba plugin was employed to identify hub genes.
2.4 Tumor microenvironment analysis
Tumor purity was calculated for each tumor sample. An estimate algorithm was used to detect the scores for immune cells and stromal cells (22). The estimate R package (https://bioinformatics.mdanderson.org/estimate/rpackage.html) was utilized to determine the stromal score, immune score, ESTIMATE score, and tumor purity of different clusters in osteosarcoma patients from the GSE14359 and GSE49003 datasets.
2.5 Tumor immune infiltration analysis
The abundance of different cell types was explored and estimated in a mixed-cell population using the CIBERSORT approach, developed by the Newman group (23). The CIBERSORT approach was applied to estimate the differences in 22 immune cell types among osteosarcoma patients in the GSE14359 and GSE49003 datasets. This analysis was performed using the R package “e1071” (Version: 1.7-3) as a precondition.
2.6 scRNA-seq data processing and analysis
The scRNA-seq data were further analyzed using the Seurat package (version 3.1.5; http://satijalab.org/seurat/) in R software (version 3.6.1) for each sample (24). In each sample, the Seurat object containing gene expression data was represented using the Read10× () function. Three quality control criteria were applied to exclude low-quality cells: (1) genes observed in fewer than 10 cells; (2) cells with fewer than 200 total observed genes; and (3) cells with more than 5% of genes expressed in mitochondria.
Gene expression data were normalized by converting values to the natural logarithm after multiplying the gene fraction by 10,000. The top 2,000 highly variable genes were selected, scaled, and analyzed using principal component analysis (PCA). Significant principal component (PC) values for cell clustering were determined using Dimheatmap and JackStrawPlot. Subsequently, t-distributed stochastic neighbor embedding (t-SNE) was performed to identify cell classification, and the maker genes of each cluster were screened by the FindAllMarkers function, with cutoff values of |log2 fold change(FC)| > 1, the cell population ratio > 0.25, and an adjusted p-value < 0.05. The SingleR (version: 1.0.6) package [16] was used for cell cluster annotation, drawing on data from the celldex package HumanPrimaryCellAtlasData. The monocle (version: 2.14.0) package (25) was used to reconstruct cell differentiation trajectories. Dimension reduction analysis was conducted using the reduceDimension function with reduction_method = “DDRTree” and max_components = 2. Characteristic genes of different cell states were filtered using the criteria of |log2FC| > 1 and adjusted p-value < 0.05 for downstream analysis.
2.7 Cell communication analysis using CellChat
Intercellular communications were clarified through ligand–receptor interactions and analyzed using the R package CellChat (26). Firstly, all ligand–receptor pairs were required to be expressed in the 8 cell subpopulations, and a ligand–receptor subnetwork was developed and saved in the human ligand–receptor pair database CellChatDB. Secondly, a biologically significant cell–cell communication network was constructed, mediated by ligand–receptor interactions. Thirdly, interacting ligand–receptor pairs related to the TNF, ANGPTL, CXCL, and IL families were further explored to examine the relationships among cell types.
3 Results
3.1 Identification of metastasis-associated genes
To identify metastasis-associated genes in osteosarcoma, we selected two osteosarcoma datasets in the GEO database: GSE14359 and GSE49003. These datasets include numerous primary and metastatic osteosarcoma tissues. In the GSE14359 dataset, 1,195 genes were upregulated and 1,435 genes were downregulated, while in the GSE49003 dataset, 561 genes were upregulated and 580 genes were downregulated (Figures 1A, B). Among these, 92 DEGs were shared between both datasets, with 41 DEGs upregulated and 51 DEGs downregulated (Figure 2A).
Figure 1. Identification of metastasis-related genes. (A) Heatmaps depicting metastasis-related genes in the GSE14359 and GSE49003 datasets. (B) Volcano plot illustrating differentially expressed genes in primary and metastatic osteosarcoma tissues.
Figure 2. Identification of metastasis-related pathways. (A) In total, 41 DEGs were upregulated and 51 DEGs were downregulated in both the GSE14359 and GSE49003 datasets. (B) GO enrichment analysis of the 92 DEGs. (C) KEGG pathway analysis of the 92 DEGs. (D) PPI network of the DEGs generated using the STRING database.
3.2 Identification of metastasis-associated pathways
GO and KECG enrichment analyses were performed on the 92 DEGs (Figures 2B, C). The GO analysis revealed that the regulation of cyclin-dependent protein kinase activity, growth plate cartilage chondrocyte differentiation, chondrocyte development involved in endochondral bone morphogenesis, and chondrocyte differentiation were enriched in the 92 DEGs (Figure 2B). In terms of cellular components, the cyclin-dependent protein kinase holoenzyme complex, transferase complex, transferring phosphorus-containing groups, SCF ubiquitin ligase complex, serine/threonine protein kinase complex, and G-protein beta/gamma-subunit complex were found to be enriched (Figure 2B). Regarding molecular function, cyclin-dependent protein serine/threonine kinase regulator activity, protein kinase regulator activity, kinase regulator activity, electron transfer activity, and S100 protein binding were involved (Figure 2B).
3.3 KEGG analysis
KEGG pathway analysis showed that the 92 DEGs were involved in regulating the cGMP-PKG and NOD-like receptor signaling pathways (Figure 2C). Furthermore, the protein–protein interaction (PPI) network of the DEGs was generated using the STRING database (Figure 2D). We found the top 15 hub genes by cytoHubba plug-in of the Cytoscape software, including Skp2, KIF20A, CCNF, TROAP, PHB, CKS1B, MCM3, CCNA2, TRIP13, CENPM, Hsp90AB1, JUN, CKS2, TK1, and KIF4A (Figure 2D). Among these, CCNA2 and KIF20A emerged as the top two candidates within the PPI network (Figure 2D). Additionally, we observed that the proportion of CD8+ T cells was higher in metastatic tissues compared to primary osteosarcoma tissues (Figure 3A). Notably, high expression of PHB was associated with a favorable prognosis in osteosarcoma patients, as shown by Kaplan–Meier survival analysis (Figure 3B).
Figure 3. CD8+ T cells are associated with metastasis. (A) Difference in the fraction of 22 immune cells across different tumor types. The proportion of CD8+ T cells was higher in metastatic tissues. (B) Kaplan–Meier survival analysis showing that high expression of PHB is associated with a favorable prognosis in osteosarcoma patients.
3.4 scRNA-seq data analysis
Next, we analyzed a total of 24,375 single-cell samples. First, we used Pearson’s correlation to measure the relationship between sequencing depth and mitochondrial gene expression (Figure 4A). Sequencing depth had a positive correlation with the number of genes, with a coefficient of 0.91 (Figure 4B). We also demonstrated the number of genes and sequencing depth in 24,375 cells from four osteosarcoma patients (Figures 4C–E). A volcano plot highlighted the highly variable genes across the cells, including MYL1, ACTC1, MYLPF, TNNC2, TNNT2, MYH3, MYOG, and PLA2G2A (Figure 4F). The principal component analysis (PCA) was performed based on these highly variable genes (Figure 4G). The determination of significant components was guided by the JackStraw function, with 20 principal components selected (Figure 4H). We then applied the t-SNE algorithm for nonlinear dimension reduction, ultimately identifying 23 clusters within the single-cell samples (Figure 5A). These 23 cell clusters were marked by a total of 4,164 genes, with the top three marker genes of each cluster shown in a heatmap (Figure 5B). Furthermore, the 23 cell clusters were classified into eight cell types: chondrocytes, MSC, T cells, monocytes, tissue stem cells, neurons, endothelial cells, and macrophages (Figure 5C).
Figure 4. scRNA-seq data analysis. (A) Pearson correlation between sequencing depth and mitochondrial gene expression. (B) The sequencing depth has a positive correlation with the number of genes, with a coefficient of 0.91. (C–E) The number of genes and sequencing depth across 24,375 cells are shown. (F) Volcano plot data show highly variable genes across cells, including MYL1, ACTC1, MYLPF, TNNC2, TNNT2, MYH3, MYOG, and PLA2G2A. (G) Principal component analysis was performed based on these highly variable genes. (H) Determination of principal component value using the JackStraw function.
Figure 5. Cell clusters and cell types are identified. (A) The t-SNE algorithm was applied to reduce the nonlinear dimension, identifying 23 clusters within the single-cell samples. (B) Heatmap data showing the top three marker genes for each cluster. (C) The 23 cell clusters were classified into eight cell types, including chondrocytes, MSC, T cells, monocytes, tissue stem cells, neurons, endothelial cells, and macrophages.
3.5 Analysis of top 15 hub genes
We further examined the expression and function of the top 15 hub genes across the eight cell types in osteosarcoma (Figures 6A, B). We found that TK1, CCNF, KIF4A, and KIF20A had higher expression levels in neurons, MSC, and chondrocytes (Figure 6B). Skp2, CCNA2, TRIP13, CENPM, and TROAP were highly expressed in neurons. JUN was more highly expressed in T cells and endothelial cells, while its expression was lower in neurons. Hsp90AB1 showed higher expression in chondrocytes and T cells but lower expression levels in tissue stem cells, monocytes, and MSC (Figures 6B, C). These 15 hub genes were expressed across various cell types (Figure 7), suggesting that they play different roles in osteosarcoma development (Figure 7).
Figure 6. Expression analysis of the top 15 hub genes. (A) Distribution of 23 clusters in primary and metastatic tissues. (B) Bubble plot showing the expression levels of the top 15 hub genes across the 23 cell clusters. (C) tSNE maps illustrating the expression of the top 15 genes within the 23 cell clusters.
Figure 7. Across various cell types, 15 hub genes are expressed.
3.6 Cell–cell communications
Finally, we analyzed cell–cell communication between different cell types using the CellChat tool. We found that there was an interaction between different cell types (Figure 8A). Each cell type interacted differently with other cell types (Figure 8B). Using netVisual bubble plots, we observed significant interactions between certain cell groups and others (Figures 9A, B). We also compared the outgoing or incoming signaling associated with each cell population (Figures 10A, B). The incoming communication of target cells and outgoing communication patterns of secreting cells highlighted the correspondence between cell groups and signaling pathways (Figure 10C). Additionally, we discovered the contribution of different cell populations to all signaling pathways (Figure 10D).
Figure 8. Comparison of cell–cell interactions among various cell types. (A) CellChat analysis of cell–cell communication across various cell types, showing an overview of the different interactions in osteosarcoma. Interaction number (left side) and interaction strength (right side) are illustrated. (B) Cell–cell communication analysis across eight cell types.
Figure 9. Interactions between various cell groups. (A) NetVisual bubble plot showing interactions between specific cell groups. (B) Heatmaps depicting interactions among eight cell types.
Figure 10. Outgoing and incoming signaling in each cell population. (A) Comparison of outgoing and incoming signaling associated with each cell population. (B) Signaling dots are illustrated. (C) Incoming and outgoing communication patterns across various cell types. (D) Contribution of cell populations to signaling pathways, with outgoing communication patterns of secreting cells (top) and incoming communication patterns of target cells (bottom) illustrated.
4 Discussion
The TME plays a crucial role in the progression and metastasis of osteosarcoma (27). Comprehensive bioinformatics analyses have been utilized to explore the mechanisms underlying osteosarcoma development, prognosis, immune microenvironment, and drug sensitivity. For example, one bioinformatics study identified the oncogenic function of FoxM1 and its association with TME, prognosis, and drug resistance in osteosarcoma (28). Another study identified 28 metastasis-related genes in osteosarcoma, including 16 upregulated and 12 downregulated genes. GO analysis suggested that these genes were enriched in Wnt and kinase pathways, bone morphogenesis, and development. KEGG pathway analysis showed that these genes were involved in the Wnt, JAK-STAT pathway, and cytokine–cytotoxic receptor interactions and were also associated with immunity in osteosarcoma patients (29). In this study, we identified 41 upregulated and 51 downregulated DEGs in both the GSE14359 and GSE49003 databases.
ScRNA-seq is a powerful tool for evaluating intratumoral heterogeneity in various cancers, including osteosarcoma (30–32). For instance, the ScRNA transcriptome has revealed the heterogeneity and immunosuppressive roles of regulatory T cells (Tregs) in osteosarcoma, with patients displaying high Tregs risk showing sensitivity to axitinib, sunitinib, and sorafenib. This finding suggests that Treg heterogeneity could open a new window for osteosarcoma treatment (33). Li et al. used ScRNA-seq data to identify natural killer (NK) cell marker genes, finding that low-risk osteosarcoma patients had high levels of infiltrating immune cells, particularly naïve CD4 and CD8 T cells, in an immunosuppressive microenvironment (34). Yi et al. reported that a proangiogenic tumor-associated macrophages (TAMs) gene risk signature could serve as a prognostic biomarker and therapeutic target in osteosarcoma (35). He et al. analyzed scRNA-seq data from lung metastasis samples and identified intratumoral heterogeneity in osteosarcoma lung metastasis, highlighting the abundance of T cells, particularly CD8+ T cells, which showed low immune checkpoints (36). In our study, we found that CD8+ T-cell proportions were higher in metastatic tissues than in primary osteosarcoma tissues.
RNA-sequencing of individual cells from lung metastatic tissues, primary tissues, and recurrent osteosarcoma tissues revealed 11 cell clusters in osteosarcoma, with lower osteoclast infiltration in chondroblastic, recurrent, and lung metastatic tissues, suggesting intratumoral heterogeneity and an immunosuppressive microenvironment in advanced osteosarcoma (37). Another study using single-cell transcriptomics observed the complexity of the TME in osteosarcoma patients without chemotherapy, identifying nine major cell types and highlighting T-cell depletion as a key feature (38). Moreover, one study integrating scRNA-seq, microarray, and bulk RNA-seq data identified tumor stem cell-related genes associated with survival, TME, and immune infiltration status in osteosarcoma, such as CKLF, DKK1, and Myc (39). Based on the scRNA-seq data, 11 clustered cell types were identified, and cell–cell interactions were found among various cell types in recurrent, metastatic, and primary osteosarcomas. NF-kB pathway was critically involved in controlling the TME in osteosarcoma (40). Herein, we identified eight cell types, including chondrocytes, MSC, T cells, monocytes, tissue stem cells, neurons, endothelial cells, and macrophages in osteosarcoma using scRNA-seq.
Based on scRNA-seq and bulk RNA-seq, Qin et al. revealed ATG16L1 as a potential immune signature and prognostic factor in osteosarcoma (41). Similarly, several survival-related genes, including MUC1, COL13A1, KAZALD1, and JAG2, were identified as prognostic markers in osteosarcoma by scRNA-seq and TARGET analysis (42). Osteoclast differentiation-related genes such as LOXL1, SERPINE2, FAM207A, TPM1, ST3GAL4, and S100A13 were also linked to prognosis (43). Additionally, genes like COL6A1, COL6A3, and MIF were associated with lung metastasis, indicating that a highly invasive subcluster in osteosarcoma correlated with a poor prognosis (44). Our study reported that high expression of PHB was associated with a favorable prognosis in osteosarcoma patients.
Using scRNA-seq, a chemoresistant risk-scoring model was established, demonstrating that scRNA-seq could be a valuable tool for evaluating chemoresistance in osteosarcoma (45). Genes such as EGFL7 and VEGF were differentially expressed in osteosarcoma (46), and NR4A1 was found to play a key role in osteosarcoma pathophysiology (47). Single-cell transcriptomics have also been employed to explore the functions of cancer-associated fibroblasts in the TME of recurrent osteosarcoma (48). A five-gene panel, including BAMBI, TMCC2, NOX4, DKK1, and CBS, was revealed as a prognostic model via scRNA-seq and bulk RNA-seq analysis (49). Similar studies integrating bulk RNA-seq and scRNA-seq have been conducted to explore the TME and metabolic pattern in osteosarcoma (50). Moreover, scRNA-seq has been used to characterize TME phenotypes in pleural effusion and tumor tissues from advanced-stage osteosarcoma patients (51). Seven genes, including ANXA1, FKBP11, SP7, TPM1, FDPS, IFITM5, and SQLE, were identified as prognostic biomarkers in mesenchymal stem cells using scRNA-seq (52). LPAR5 was identified as a potential indicator for TME remodeling and a therapeutic target for osteosarcoma based on scRNA-seq data (53). Additionally, five key genes, including CD4, RUNX2, OMD, COL9A3, and JUN, were found to be associated with osteosarcoma progression (54). GNG4 (55) and C1Q+ tumor-associated macrophages (56) were reported to predict osteosarcoma prognosis using scRNA-seq. Our study identified 15 hub genes expressed across various cell types, including TK1, CCNF, KIF4A, KIF20A, Skp2, CCNA2, TRIP13, CENPM, and TROAP.
Among these hub genes, Skp2 has been reported as an E3 ubiquitin ligase involved in osteosarcoma progression. Studies have shown that suppression of Skp2 attenuates cell proliferation and invasion in osteosarcoma, while overexpression of Skp2 enhances cell growth and motility (57, 58). Additionally, Skp2 has been implicated in promoting epithelial-mesenchymal transition (EMT) in osteosarcoma cells (59) and facilitating stemness and tumor progression through interaction with p27 (60). Targeting Skp2 inhibition has demonstrated anticancer effects, leading to survival benefits in osteosarcoma (61). It has been known that retinoblastoma (RB) and p53 pathways are critically involved in aging, senescence, and tumorigenesis (62, 63). Skp2 knockout has been shown to trigger immune infiltration and improve prognosis in Rb1/p53-deficient mice with osteosarcoma (64).
Recently, Lu et al. used integrative analysis of single-cell and bulk transcriptome data to discover neutrophil-related genes in osteosarcoma, such as C3AR1 and FCER1G (65). Zheng et al. employed single-cell transcriptomics to uncover chemotherapy-mediated remodeling of the TME in osteosarcoma (66). Yi et al. identified a pro-protein synthesis osteosarcoma subtype using scRNA-se1, which could be useful for predicting prognosis and treatment (67). Additionally, Truong et al. discovered a single-cell differentiation landscape in osteosarcoma, which may aid in developing targeted therapy (68).
5 Conclusion
In this study, we identified metastasis-associated genes and pathways in osteosarcoma. Additionally, we verified the presence of 8 cell types, including chondrocytes, MSC, T cells, monocytes, tissue stem cells, neurons, endothelial cells, and macrophages, in osteosarcoma using scRNA-seq. Furthermore, 15 hub genes were observed in various cell types. Our study provides valuable insights into the metastatic mechanisms of osteosarcoma. However, it is important to acknowledge the limitations of this study. The findings from scRNA-seq have not been validated through cellular experiments and mouse studies. Additionally, metastasis-associated genes and pathways in osteosarcoma have not been clarified through clinical studies. To validate our findings, in vitro experiments and in vivo studies are required. It is also necessary to examine the expression of the 15 hub genes in metastatic tissues of osteosarcoma patients. By unraveling the complex interactions between osteosarcoma cells and the surrounding microenvironment at the single-cell level, we can identify novel therapeutic targets to inhibit metastasis and improve treatment outcomes. Moreover, scRNA-seq could be instrumental in uncovering mechanisms of resistance to current therapies, leading to the development of combination therapies that overcome resistance. Future research should explore how the TME evolves during tumor progression and how certain microenvironmental niches contribute to osteosarcoma metastasis.
Data availability statement
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.
Author contributions
JL: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Writing – original draft. YB: Formal analysis, Investigation, Methodology, Resources, Software, Writing – original draft. HZ: Data curation, Formal analysis, Investigation, Methodology, Software, Writing – original draft. TC: Conceptualization, Investigation, Project administration, Supervision, Validation, Visualization, Writing – review & editing. GS: Conceptualization, Project administration, Supervision, Validation, Visualization, Writing – review & editing.
Funding
The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.
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.
References
1. Beird HC, Bielack SS, Flanagan AM, Gill J, Heymann D, Janeway KA, et al. Osteosarcoma. Nat Rev Dis Primers. (2022) 8:77. doi: 10.1038/s41572-022-00409-y
2. Gill J, Gorlick R. Advancing therapy for osteosarcoma. Nat Rev Clin Oncol. (2021) 18:609–24. doi: 10.1038/s41571-021-00519-8
3. Yan P, Wang J, Yue B, Wang X. Unraveling molecular aberrations and pioneering therapeutic strategies in osteosarcoma. Biochim Biophys Acta Rev Cancer. (2024) 1879:189171. doi: 10.1016/j.bbcan.2024.189171
4. Twenhafel L, Moreno D, Punt T, Kinney M, Ryznar R. Epigenetic changes associated with osteosarcoma: A comprehensive review. Cells. (2023) 12(12):1595. doi: 10.3390/cells12121595
5. Zhang Y, Xu Y, Bao Y, Luo Y, Qiu G, He M, et al. N6-methyladenosine (m6A) modification in osteosarcoma: expression, function and interaction with noncoding RNAs - an updated review. Epigenetics. (2023) 18:2260213. doi: 10.1080/15592294.2023.2260213
6. Almansa-Gomez S, Prieto-Ruiz F, Cansado J, Madrid M. Autophagy modulation as a potential therapeutic strategy in osteosarcoma: current insights and future perspectives. Int J Mol Sci. (2023) 24(18):13827. doi: 10.3390/ijms241813827
7. Qin S, Wang Y, Ma C, Lv Q. Competitive endogenous network of circRNA, lncRNA, and miRNA in osteosarcoma chemoresistance. Eur J Med Res. (2023) 28:354. doi: 10.1186/s40001-023-01309-x
8. Wang S, Ren Q, Li G, Zhao X, Zhao X, Zhang Z. The targeted therapies for osteosarcoma via six major pathways. Curr Mol Pharmacol. (2024) 17:e210823220109. doi: 10.2174/1874467217666230821142839
9. Ji Z, Shen J, Lan Y, Yi Q, Liu H. Targeting signaling pathways in osteosarcoma: Mechanisms and clinical studies. MedComm (2020). (2023) 4:e308. doi: 10.1002/mco2.308
10. Entenberg D, Oktay MH, Condeelis JS. Intravital imaging to study cancer progression and metastasis. Nat Rev Cancer. (2023) 23:25–42. doi: 10.1038/s41568-022-00527-5
11. Li YQ, Sun FZ, Li CX, Mo HN, Zhou YT, Lv D, et al. RARRES2 regulates lipid metabolic reprogramming to mediate the development of brain metastasis in triple negative breast cancer. Mil Med Res. (2023) 10:34. doi: 10.1186/s40779-023-00470-y
12. Sheng G, Gao Y, Yang Y, Wu H. Osteosarcoma and metastasis. Front Oncol. (2021) 11:780264. doi: 10.3389/fonc.2021.780264
13. Du X, Wei H, Zhang B, Wang B, Li Z, Pang LK, et al. Molecular mechanisms of osteosarcoma metastasis and possible treatment opportunities. Front Oncol. (2023) 13:1117867. doi: 10.3389/fonc.2023.1117867
14. Hou B, Chen T, Zhang H, Li J, Wang P, Shang G. The E3 ubiquitin ligases regulate PD-1/PD-L1 protein levels in tumor microenvironment to improve immunotherapy. Front Immunol. (2023) 14:1123244. doi: 10.3389/fimmu.2023.1123244
15. Walsh LA, Quail DF. Decoding the tumor microenvironment with spatial technologies. Nat Immunol. (2023) 24:1982–93. doi: 10.1038/s41590-023-01678-9
16. Tian B, Du X, Zheng S, Zhang Y. The role of tumor microenvironment in regulating the plasticity of osteosarcoma cells. Int J Mol Sci. (2022) 23(24):16155. doi: 10.3390/ijms232416155
17. Baron M, Drohat P, Crawford B, Hornicek FJ, Best TM, Kouroupis D. Mesenchymal stem/stromal cells: immunomodulatory and bone regeneration potential after tumor excision in osteosarcoma patients. Bioeng (Basel). (2023) 10(10):1187. doi: 10.3390/bioengineering10101187
18. Zeng J, Peng Y, Wang D, Ayesha K, Chen S. The interaction between osteosarcoma and other cells in the bone microenvironment: From mechanism to clinical applications. Front Cell Dev Biol. (2023) 11:1123065. doi: 10.3389/fcell.2023.1123065
19. Jin J, Cong J, Lei S, Zhang Q, Zhong X, Su Y, et al. Cracking the code: Deciphering the role of the tumor microenvironment in osteosarcoma metastasis. Int Immunopharmacol. (2023) 121:110422. doi: 10.1016/j.intimp.2023.110422
20. Xiong D, Zhang L, Sun ZJ. Targeting the epigenome to reinvigorate T cells for cancer immunotherapy. Mil Med Res. (2023) 10:59. doi: 10.1186/s40779-023-00496-2
21. Cheng S, Wang H, Kang X, Zhang H. Immunotherapy innovations in the fight against osteosarcoma: emerging strategies and promising progress. Pharmaceutics. (2024) 16(2):251. doi: 10.3390/pharmaceutics16020251
22. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612
23. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. (2019) 37:773–82. doi: 10.1038/s41587-019-0114-2
24. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. (2018) 36:411–20. doi: 10.1038/nbt.4096
25. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402
26. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
27. Nirala BK, Yamamichi T, Petrescu DI, Shafin TN, Yustein JT. Decoding the impact of tumor microenvironment in osteosarcoma progression and metastasis. Cancers (Basel). (2023) 15(20):5108. doi: 10.3390/cancers15205108
28. Shi S, Wang Q, Du X. Comprehensive bioinformatics analysis reveals the oncogenic role of FoxM1 and its impact on prognosis, immune microenvironment, and drug sensitivity in osteosarcoma. J Appl Genet. (2023) 64:779–96. doi: 10.1007/s13353-023-00785-5
29. Qin S, Li L, Liu D. Metastasis-related gene signature associates with immunity and predicts prognosis accurately in patients with osteosarcoma. Aging (Albany NY). (2023) 15:7219–36. doi: 10.18632/aging.204902
30. Thomas DD, Lacinski RA, Lindsey BA. Single-cell RNA-seq reveals intratumoral heterogeneity in osteosarcoma patients: A review. J Bone Oncol. (2023) 39:100475. doi: 10.1016/j.jbo.2023.100475
31. Feng DC, Zhu WZ, Wang J, Li DX, Shi X, Xiong Q, et al. The implications of single-cell RNA-seq analysis in prostate cancer: unraveling tumor heterogeneity, therapeutic implications and pathways towards personalized therapy. Mil Med Res. (2024) 11:21. doi: 10.1186/s40779-024-00526-7
32. Tang L, Huang ZP, Mei H, Hu Y. Insights gained from single-cell analysis of chimeric antigen receptor T-cell immunotherapy in cancer. Mil Med Res. (2023) 10:52. doi: 10.1186/s40779-023-00486-4
33. Cheng D, Zhang Z, Mi Z, Tao W, Liu D, Fu J, et al. Deciphering the heterogeneity and immunosuppressive function of regulatory T cells in osteosarcoma using single-cell RNA transcriptome. Comput Biol Med. (2023) 165:107417. doi: 10.1016/j.compbiomed.2023.107417
34. Li Q, Huang X, Zhao Y. Prediction of prognosis and immunotherapy response with a novel natural killer cell marker genes signature in osteosarcoma. Cancer Biother Radiopharm. (2023). doi: 10.1089/cbr.2023.0103
35. Yi C, Li Z, Zhao Q, Gong D, Zhao S, Chen Z, et al. Single-cell RNA sequencing pro-angiogenic macrophage profiles reveal novel prognostic biomarkers and therapeutic targets for osteosarcoma. Biochem Genet. (2023) 62(2):1325–46. doi: 10.1007/s10528-023-10483-w
36. He M, Jiang X, Miao J, Feng W, Xie T, Liao S, et al. A new insight of immunosuppressive microenvironment in osteosarcoma lung metastasis. Exp Biol Med (Maywood). (2023) 248:1056–73. doi: 10.1177/15353702231171900
37. 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
38. Liu Y, Feng W, Dai Y, Bao M, Yuan Z, He M, et al. Single-cell transcriptomics reveals the complexity of the tumor microenvironment of treatment-naive osteosarcoma. Front Oncol. (2021) 11:709210. doi: 10.3389/fonc.2021.709210
39. Xu A, Qian C, Lin J, Yu W, Jin J, Liu B, et al. Cell differentiation trajectory-associated molecular classification of osteosarcoma. Genes (Basel). (2021) 12(11):1685. doi: 10.3390/genes12111685
40. Wu R, Dou X, Li H, Sun Z, Li H, Shen Y, et al. Identification of cell subpopulations and interactive signaling pathways from a single-cell RNA sequencing dataset in osteosarcoma: A comprehensive bioinformatics analysis. Front Oncol. (2022) 12:853979. doi: 10.3389/fonc.2022.853979
41. Qin Z, Luo K, Liu Y, Liao S, He J, He M, et al. ATG16L1 is a potential prognostic biomarker and immune signature for osteosarcoma: A study based on bulk RNA and single-cell RNA-sequencing. Int J Gen Med. (2022) 15:1033–45. doi: 10.2147/IJGM.S341879
42. Feleke M, Feng W, Rothzerg E, Song D, Wei Q, Koks S, et al. Single-cell RNA-seq identification of four differentially expressed survival-related genes by a TARGET: Osteosarcoma database analysis. Exp Biol Med (Maywood). (2022) 247:921–30. doi: 10.1177/15353702221080131
43. Shao H, Ge M, Zhang J, Zhao T, Zhang S. Osteoclasts differential-related prognostic biomarker for osteosarcoma based on single cell, bulk cell and gene expression datasets. BMC Cancer. (2022) 22:288. doi: 10.1186/s12885-022-09380-z
44. Guo J, Tang H, Huang P, Guo J, Shi Y, Yuan C, et al. Single-cell profiling of tumor microenvironment heterogeneity in osteosarcoma identifies a highly invasive subcluster for predicting prognosis. Front Oncol. (2022) 12:732862. doi: 10.3389/fonc.2022.732862
45. Zeng Z, Li W, Zhang D, Zhang C, Jiang X, Guo R, et al. Development of a chemoresistant risk scoring model for prechemotherapy osteosarcoma using single-cell sequencing. Front Oncol. (2022) 12:893282. doi: 10.3389/fonc.2022.893282
46. Feleke M, Feng W, Song D, Li H, Rothzerg E, Wei Q, et al. Single-cell RNA sequencing reveals differential expression of EGFL7 and VEGF in giant-cell tumor of bone and osteosarcoma. Exp Biol Med (Maywood). (2022) 247:1214–27. doi: 10.1177/15353702221088238
47. Liu W, Hao Y, Tian X, Jiang J, Qiu Q. The role of NR4A1 in the pathophysiology of osteosarcoma: A comprehensive bioinformatics analysis of the single-cell RNA sequencing dataset. Front Oncol. (2022) 12:879288. doi: 10.3389/fonc.2022.879288
48. Huang X, Wang L, Guo H, Zhang W, Shao Z. Single-cell transcriptomics reveals the regulative roles of cancer associated fibroblasts in tumor immune microenvironment of recurrent osteosarcoma. Theranostics. (2022) 12:5877–87. doi: 10.7150/thno.73714
49. Yang J, Zhang J, Na S, Wang Z, Li H, Su Y, et al. Integration of single-cell RNA sequencing and bulk RNA sequencing to reveal an immunogenic cell death-related 5-gene panel as a prognostic model for osteosarcoma. Front Immunol. (2022) 13:994034. doi: 10.3389/fimmu.2022.994034
50. Huang R, Wang X, Yin X, Zhou Y, Sun J, Yin Z, et al. Combining bulk RNA-sequencing and single-cell RNA-sequencing data to reveal the immune microenvironment and metabolic pattern of osteosarcoma. Front Genet. (2022) 13:976990. doi: 10.3389/fgene.2022.976990
51. Zhang Z, Ji W, Huang J, Zhang Y, Zhou Y, Zhang J, et al. Characterization of the tumour microenvironment phenotypes in Malignant tissues and pleural effusion from advanced osteoblastic osteosarcoma patients. Clin Transl Med. (2022) 12:e1072. doi: 10.1002/ctm2.1072
52. Jiang H, Du H, Liu Y, Tian X, Xia J, Yang S. Identification of novel prognostic biomarkers for osteosarcoma: a bioinformatics analysis of differentially expressed genes in the mesenchymal stem cells from single-cell sequencing data set. Transl Cancer Res. (2022) 11:3841–52. doi: 10.21037/tcr-22-2370
53. He Y, Zhou H, Huang X, Qu Y, Wang Y, Pei W, et al. Infiltration of LPAR5(+) macrophages in osteosarcoma tumor microenvironment predicts better outcomes. Front Immunol. (2022) 13:909932. doi: 10.3389/fimmu.2022.909932
54. Wang Y, Qin D, Gao Y, Zhang Y, Liu Y, Huang L. Identification of therapeutic targets for osteosarcoma by integrating single-cell RNA sequencing and network pharmacology. Front Pharmacol. (2022) 13:1098800. doi: 10.3389/fphar.2022.1098800
55. Jiang X, Tang F, Zhang J, He M, Xie T, Tang H, et al. High GNG4 predicts adverse prognosis for osteosarcoma: Bioinformatics prediction and experimental verification. Front Oncol. (2023) 13:991483. doi: 10.3389/fonc.2023.991483
56. Tu J, Wang D, Zheng X, Liu B. Single-cell RNA datasets and bulk RNA datasets analysis demonstrated C1Q+ tumor-associated macrophage as a major and antitumor immune cell population in osteosarcoma. Front Immunol. (2023) 14:911368. doi: 10.3389/fimmu.2023.911368
57. Ding L, Li R, Han X, Zhou Y, Zhang H, Cui Y, et al. Inhibition of Skp2 suppresses the proliferation and invasion of osteosarcoma cells. Oncol Rep. (2017) 38:933–40. doi: 10.3892/or.2017.5713
58. Ding L, Li R, Sun R, Zhou Y, Zhou Y, Han X, et al. S-phase kinase-associated protein 2 promotes cell growth and motility in osteosarcoma cells. Cell Cycle. (2017) 16:1547–55. doi: 10.1080/15384101.2017.1346760
59. Ding L, Wang C, Cui Y, Han X, Zhou Y, Bai J, et al. S-phase kinase-associated protein 2 is involved in epithelial-mesenchymal transition in methotrexate-resistant osteosarcoma cells. Int J Oncol. (2018) 52:1841–52. doi: 10.3892/ijo.2018.4345
60. Wang J, Aldahamsheh O, Ferrena A, Borjihan H, Singla A, Yaguare S, et al. The interaction of SKP2 with p27 enhances the progression and stemness of osteosarcoma. Ann N Y Acad Sci. (2021) 1490:90–104. doi: 10.1111/nyas.14578
61. Wang J, Ferrena A, Zhang R, Singh S, Viscarret V, Al-Harden W, et al. Targeted inhibition of SCF(SKP2) confers anti-tumor activities resulting in a survival benefit in osteosarcoma. Oncogene. (2024) 43:962–75. doi: 10.1038/s41388-024-02942-4
62. Huang Y, Che X, Wang PW, Qu X. p53/MDM2 signaling pathway in aging, senescence and tumorigenesis. Semin Cancer Biol. (2024) 101:44–57. doi: 10.1016/j.semcancer.2024.05.001
63. Pareek A, Kumar D, Pareek A, Gupta MM, Jeandet P, Ratan Y, et al. Retinoblastoma: An update on genetic origin, classification, conventional to next-generation treatment strategies. Heliyon. (2024) 10:e32844. doi: 10.1016/j.heliyon.2024.e32844
64. Ferrena A, Wang J, Zhang R, Karadal-Ferrena B, Al-Hardan W, Singh S, et al. SKP2 knockout in rb1/p53-deficient mouse models of osteosarcoma induces immune infiltration and drives a transcriptional program with a favorable prognosis. Mol Cancer Ther. (2024) 23:223–34. doi: 10.1158/1535-7163.MCT-23-0173
65. Lu J, Rui J, Xu XY, Shen JK. Exploring the role of neutrophil-related genes in osteosarcoma via an integrative analysis of single-cell and bulk transcriptome. Biomedicines. (2024) 12(7):1513. doi: 10.3390/biomedicines12071513
66. Zheng X, Wu W, Zhao Z, Zhang X, Yu S. Single-cell transcriptomic insights into chemotherapy-induced remodeling of the osteosarcoma tumor microenvironment. J Cancer Res Clin Oncol. (2024) 150:356. doi: 10.1007/s00432-024-05787-2
67. Yi C, Liu J, Zhao S, Gong D, Xu B, Li A, et al. Identification of a pro-protein synthesis osteosarcoma subtype for predicting prognosis and treatment. Sci Rep. (2024) 14:16475. doi: 10.1038/s41598-024-67547-z
Keywords: osteosarcoma, immunotherapy, therapy, RNA sequencing, TME
Citation: Li J, Bai Y, Zhang H, Chen T and Shang G (2024) Single-cell RNA sequencing reveals the communications between tumor microenvironment components and tumor metastasis in osteosarcoma. Front. Immunol. 15:1445555. doi: 10.3389/fimmu.2024.1445555
Received: 07 June 2024; Accepted: 22 August 2024;
Published: 11 September 2024.
Edited by:
Xiangpeng Dai, Jilin University, ChinaReviewed by:
Wangpan Shi, University of California, San Diego, United StatesDong Ren, UC Irvine Medical Center, United States
Rong Li, Xinjiang Medical University, China
Copyright © 2024 Li, Bai, Zhang, Chen and Shang. 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: Ting Chen, chenting88613@163.com; Guanning Shang, gnshang@cmu.edu.cn
†These authors have contributed equally to this work