- Open Access
Abnormal signal pathways and tumor heterogeneity in osteosarcoma
Journal of Translational Medicine volume 21, Article number: 99 (2023)
Osteosarcoma (OS) is the most frequent and aggressive primary malignant sarcoma among adolescents and chemotherapy has not substantially progressed for decades. New insights into OS development and therapeutic strategies are urgently needed.
We analyzed integrated single-cell transcriptomes, bulk RNA-seq, and microarray data from Gene Expression Omnibus (GEO) datasets. We also used Weighted Gene Co-expression Network Analysis (WGCNA), Gene set enrichment analysis (GSEA), and Gene set variation analysis (GSVA), along with Simple ClinVar and Enrichr web servers.
The findings of integrated single-cell analysis showed that OS arises from imperfect osteogenesis during development. Novel abnormalities comprised deficient TGFβ and P53 signal pathways, and cell cycle pathway activation, and a potentially new driver mutation in the interferon induced transmembrane protein 5 (IFITM5) that might function as a pathogenic factor in OS. Osteosarcoma is characterized by oncocyte heterogeneity, especially in immunogenic and adipocyte-like subtypes that respectively promote and hamper OS treatment. Etoposide is a promising chemotherapeutic that provides palliation by affecting the subtype of OS and correcting the abnormal pathways.
Various abnormal signal pathways play indispensable roles in OS development. We explored the heterogeneity and underlying mechanisms of OS and generated findings that will assist with OS assessment and selecting optimal therapies.
Osteosarcoma (OS) is an extremely aggressive sarcoma for which chemotherapeutic strategies have not substantially progressed for decades [1, 2]. The mechanism of OS awaits explanation and new molecular targets are urgently needed to improve the prognosis of patients. The causes of OS remain unclear in terms of oncogenetic mechanisms, but it might be due to genetic abnormalities, changes in tumor suppressor pathways (such as p53 and Rb), telomerase, and alternative telomere lengthening (ALT) inside bone cells . Osteosarcoma driver mutations have been found in TP53, RB, PTCH1, MYC, NOTCH1, BRCA2, APC, and PRKAR1A genes. Clinical evaluations of targeted agents have yielded disappointing results [4, 5]. However, some new driver genes are currently under development as potential therapies for OS.
Interrelated biological signaling pathways in OS comprise WNT/βcatenin, Hedgehog, mTOR, and RANKL/NF-κB . Achievements in the underlying molecular biology of OS have contributed to transformative advances in understanding this type of malignancy. Osteosarcoma originates from bone mesenchymal stem cells (MSCs) with osteoblastic lineage commitment. However, understanding the molecular mechanisms of OS development remains insufficient [3, 7]. For example, the TGFβ signal pathway exerts dual effects in terms of tumor prevention and carcinogenesis, depending on the timing of tumor development [8,9,10].
Osteosarcoma is significantly heterogeneous at the genomic, transcriptomic, and epigenetic levels. Malignant OS cells have stemness properties that are closely involved in chemotherapy resistance, relapse, or refractory and metastatic processes [11,12,13]. Current approaches to intra-tumor heterogeneity mainly depend on cellular genomic or transcriptional sequencing. However, this lacks high resolution in authenticating the complex cellular subtypes and intra-tumoral heterogeneity in OS. Single-cell RNA-sequencing (scRNA-seq) is a promising tool that can analyze the composition of heterogeneous cell populations and lineage hierarchies [14, 15].
Here we mapped single-cell (sc) RNA sequences of human embryonic long bones, normal bone tissue, and OS at single-cell high resolution to address the development of OS and oncocyte heterogeneity. The results of our integrated bioinformatic analysis identified and validated the process of osteogenesis imperfecta, potential new driver genes, abnormal signaling pathways, transcriptomic heterogeneity, and targeted therapeutic agents for OS.
Single cell RNA sequencing (scRNA-seq) data processing
We analyzed scRNA-seq data using R statistical software (v. 4.1.2, The R Foundation for Statistical Computing, Vienna, Austria), the Seurat v4.0 toolkit and MPLAB Harmony v1.0. Public scRNA-seq data were acquired from the GSE143753, GSE196678, and GSE162454 datasets in the Gene Expression Omnibus (GEO). Low-quality cells (< 3; features, < 200; > 10% mitochondrial genes) were filtered, the data were initially normalized, then batch effects were removed using LogNormalize (features = 3000) and the Harmony function (max.iter.harmony = 20). Dimensions were reduced using Uniform Manifold Approximation and Projection (UMAP). We found marker genes for cell groups in datasets using the COSine similarity-based marker gene identification (COSG) package . We quantified marker genes using Z-scores and the results are shown in a heatmap. We calculated copy number variations to identify malignant oncocytes using the inferCNV package.
Cell types were identified in each subpopulation based on their known lineage markers: mesenchymal cells, vimentin (VIM); immune cells, cluster of differentiation (CD)45; myeloid cells: CD83, CD14, CD68; lymphoid cells, CD3 and B-cell maturation antigen (BCMA); endothelial cells, von Willebrand factor, (VWF) and carcinoembryonic antigen-related cell adhesion molecule-1 (CEACAM1); myogenic cells, myogenic differentiation 1 (MYOD1) and myogenin (MYOG); bone-related cells, collagen type 1 alpha 1 (COL1A1) and alpha 2 (COL1A2) chains, lumican (LUM), SRY-Box Transcription Factor 9 (SOX9), alkaline phosphatase (ALP) and RUNX Family Transcription Factor 2 (RUNX2)).
Pseudotime and trajectory analysis
Single-cell pseudotime and trajectory analyses were constructed using the Monocle3 toolkit (v. 2.14.0). Evolutionary processes were organized into potentially discontinuous trajectories by the learn_graph function. Pseudotime was defined using the order_cells function with a selected node representing development. Genes that were differentially expressed over these trajectories were then identified [17, 18].
Weighted gene co-expression network analysis (WGCNA)
We also searched related public expression profiles in 12 MSCs, 3 osteoblasts, and 84 OS samples (GSE33383). A co-expression network was constructed to identify co-expressed modules using the WGCNA package in R. The expression matrix was restricted to only the top 25% of expressed genes according to variance analyses. Relationships between gene sets (modules) were explored using a hierarchical cluster dendrogram. A clustering tree based on the eigengenes of modules calculated the dissimilarity of the module eigengenes (MEs). Associations between gene sets (modules) and clinical features were assessed using Pearson correlations.
Clinical datasets and resources
We downloaded 53 samples including expression profiles and clinical outcomes of patients with OS from GSE19743 for further analysis. Except for the expression profile, the clinical information mainly included living status, survival duration, Huvos grading scores, and metastatic status. Patients were assigned to high- (n = 27) and low- (n = 26) risk groups based on the median numbers of samples. We applied Cox regression and Kaplan–Meier curves to estimate overall survival.
Microarray data processing
We used the gene set with the expression profile GSE84500 to obtain significantly differentially expressed genes (DEGs). A data matrix was downloaded to annotate the probe into gene symbol sets. Significance was analyzed using Limma, or the Deseq2 package. The most significant changes in gene expression were identified using false discovery rates (FDR < 0.01) and absolute fold change (FC > 1).
Gene set enrichment analysis (GSEA)
We assessed gene set enrichment using GSEA . Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway gene sets were curated from the Molecular Signatures Database (Human MSigDB database v2022.1.Hs). Immunogenic and adipocyte-like subtype molecular signatures were created by identifying significant marker genes in scRNA-seq. Normalized enrichment scores (NES) > 1 and p < 0.05 reliably filtered significant pathways.
Gene set variation analysis (GSVA)
We established the molecular signatures of OS oncocyte subtypes depending on the significant marker genes, which were acquired by single cell analysis . We also calculated final scores for further analyses using GSVA (v. 1.47.0).
Databases and datasets
Simple ClinVar (https://simple-clinvar.broadinstitute.org/) is a web-server application that summarize variant, gene, gene- and disease-wise statistics based on the entire ClinVar database in a dynamic and user-friendly web-interface . Enrichr (https://maayanlab.cloud/Enrichr/) is a robust web-server that contains many types datasets, among which we used the MGI Mammalian Phenotype, Jensen DISEASES, and BioPlanet pathway databases, as well as the PanglaoDB Augmented, and DrugMatrix datasets [22, 23]. Statistical significance was set at < 0.05.
Osteosarcoma arises from osteogenesis imperfecta during development
We integrated data from 8-week-old long bones (GSE143753), human normal bone tissue (GSE196678) and OS (GSE162454) to determine whether a developmental disorder is associated with OS. After quality control and batch correction, we obtained 60,553 cells for downstream analysis. After dimensionality reduction, UMAP-based cell clustering identified myeloid, lymphoid, endothelial, myogenic, and bone-related cell clusters (Fig. 1A). Based on the lineage markers (Fig. 1B), then split them depending on their tissue origins. We found that bone mesenchymal stem cells (BMCs), normal osteocytes, and OS oncocytes mapped together, indicating that oncocytes might be derived from BMSCs due to a developmental disorder (Fig. 1C). We compared transcriptional profiles to identify DEGs between normal osteocytes and OS oncocytes during development. We identified 349 DEGs (log2FC > 1; P < 0.01), of which 143 were elevated in OS oncocytes (Fig. 1D; Additional file 2: Table S1). We analyzed the upregulated DEGs using the MGI Mammalian Phenotype and Jensen DISEASES databases to gain insight into their functional relevance. We found enrichment of decreased compact bone thickness, decreased trabecular bone volume, decreased length of long bones, and osteogenesis imperfecta (Fig. 1E and F; black arrows). This further supported the notion of a developmental disorder in OS. We analyzed trajectories to gain insight into the origin and development of bone-related cells (Fig. 2A and B). The results revealed partial BMSCs at the start of the developmental trajectory. This indicated that BMSCs can develop into OS.
Interferon induced transmembrane protein 5 (IFITM5) is a potential new driver gene in OS development
Among the functional ontologies, only IFITM5 played a crucial role in the process of OS oncogenesis, including osteogenesis imperfecta, and decreased bone thickness, volume, and length (Fig. 2C). Furthermore, IFITM5 and Metallothionein 2A (MT2A) were the most significantly upregulated and downregulated genes, respectively, among the DEGs between normal osteocytes and OS oncocytes (Fig. 1D). The trajectory analyses also offered a clearer view of gene dynamics along a single path. We found that IFITM5 significantly increased from BMSCs to OS, whereas MT2A did not (Fig. 2D). The pseudotime-dependent IFITM5 gene changed as cells progressed along the oncocyte trajectory. Therefore knowing the order in which IFITM5 changes can lead to the generation of new models of OS development (Fig. 2E and F). We confirmed the IFITM5 variant using Simple ClinVar. We summarized various mutation categories of IFITM5 and found that a missense mutation was the most prevalent type (Fig. 2G), and that three variations of ITITM5, including the point mutation (c.-14C > T), the missense c.119C > T (p.S40L), and c.143A > G (p.N48S), can lead to pathogenic states according to the Simple ClinVar dataset and current literature (Fig. 2H). These findings together indicated that IFITM is a potential new driver gene in OS development that warrants further investigation.
Osteosarcoma was attributed to complicated regulation of abnormal signal pathways during development
We further investigated OS development by analyzing the expression profiles of 12 MSCs, 3 osteoblasts, and 84 OS samples (GSE33383). We established correlation networks and identified co-expression modules using WGCNA. Nine co-expression modules with 4,637 genes were obtained using the Dynamic Tree Cut algorithm (Fig. 3A). The clustering tree mainly showed two branches based on the module eigengenes (Fig. 3B). Correlations between modules and samples of MSCs, osteoblasts, and OS are shown in (Fig. 3C). The heat map shows that MEgreen and MEblue modules correlated negatively with MSCs and osteoblasts, and positively with OS. The MEgery module correlated positively with MSCs and osteoblasts, but negatively with OS. We analyzed BioPlanet pathways using the co-expressed genes in the MEgreen, MEblue and MEgery modules to gain insight into the signal pathway relevance (Additional file 3: Table S2). The three most significant pathways in each module, namely, Cell cycle, Antigen processing and presentation, and TGFβ regulation of extracellular matrix were enriched (Fig. 3D‒F). The module correlation coefficients for MSCs were consistent with those of osteoblasts, in contrast to OS (Fig. 3G and H). That is, OS was attributed to constitutive activation of the Cell cycle and Antigen processing and presentation pathways, and deficient TGFβ regulation of extracellular matrix signals. The atypical signal pathways also had complex connections to regulate bone development (ConsensusPathDB; Fig. 3I).
Transcriptomic heterogeneity of oncocytes in OS
We analyzed the transcriptomic heterogeneity of OS oncocytes by re-clustering bone-related cells in OS (n = 8758 cells) and found six distinct transcriptomic sub-clusters (Fig. 4A). The positive expression of COL1A1, COL1A2, ALP and the absence of CD45 revealed that the six sub-clusters were bone-related oncocytes (Additional file 1: Fig. S1 B). The copy number variation (CNV) was estimated using inferCNV (with myogenic cells as a reference) to distinguish oncocytes and normal cells (Fig. 4B). The CNV from scRNA-seq revealed numerous alterations in all sub-clusters compared with the reference cells (Fig. 4C). We determined the top 100 differentially expressed transcripts among the six sub-populations and defined their cell types in the Enrichr datasets (Additional file 4: Table S3). A heat map shows the hub genes in each sub-cluster (Fig. 4D). Six sub-populations were annotated and visualized by UMAP as neuron-like, immunogenic, fibroblastic, chondroblastic, adipocyte-like, and osteoblastic oncocytes (Fig. 4A). The top 100 gene signatures of each sub-population were assessed using GSVA and the oncocyte subtype abundance was estimated from 53 human transcriptomic profiles (GSE21257). Kaplan–Meier estimates associated higher immunogenic and adipocyte-like scores with longer progression-free survival (p = 0.0202), whereas and unfavorable prognoses, respectively (p = 0.0065). Other subtypes had no significant prognostic value (Fig. 4E). IFITM5 is a potential new driver gene in OS development in our study. Correlation analysis was performed to further identify the relationship between tumor subtypes and IFITM5. We found that the expression of IFITM5 has a strongly significant positive correlation with adipocyte-like and chondroblastic subtypes and a significant negative correlation with immunogenic and fibroblastic subtypes. Therefore, the high expression of IFITM5 was associated with the shorter overall survival subtype which is in line with the previous study (Additional file 1: Fig. S1A) .
Deficient TGFβ signaling leads to adipocyte-like subtype of OS during development
The adipocyte-like OS subtype was associated with a poor prognosis. We profiled gene expression in human MSCs (hMSCs) in the presence of adipogenic/osteogenic factors (GSE84500) to assess the development of adipocyte-like subtype of OS. A comparison of transcriptional profiles between panels of hMSC_BMP2 + IBMX and hMSC_BMP2 + IBMX + TGFβ. Principal Component Analysis (PCA) initially revealed correlations (Fig. 5A). We identified 469 upregulated transcripts (log2FC > 1, P < 0.01), of which 240 were elevated in the hMSC_BMP2 + IBMX group (Additional file 5: Table S4) and named Adipocytes in PanglaoDB Augmented Datasets and 229 that were abundantly expressed in the hMSC_BMP2 + IBMX + TGFβ group (Additional file 6: Table S5) were labeled Chondrocytes, Osteoblasts, and Fibroblasts in the same datasets (Fig. 5B). Meanwhile, BioPlanet pathways confirmed activated TGFβ regulation of extracellular matrix in the hMSC_BMP2 + IBMX + TGFβ group. Therefore, we judged that BMP2 and IBMX cause hMSCs to differentiate into adipocytes, TGFβ inhibits this process and redirects these cells to differentiate into chondrocytes, osteoblasts, and fibroblasts. The TGFβ signal was defective during OS development (Fig. 3C and F), and adipocytes could not differentiate into normal bone cells without the TGFβ signal. This stalled adipocyte development and led to the adipocyte-like subtype of OS (Fig. 5D). In clinical translation, patients with metastasis are more likely to have high-risk scores (Fig. 5E). We applied the Huvos grading system to judge the effectiveness of neoadjuvant chemotherapy on OS; higher scores reflect larger necrotic areas in the clinical data (GSE21257). However, the risk score of the adipocyte-like subtype did not differ among grades, indicating that the adipocyte-like subtype is resistant to neoadjuvant chemotherapy (Fig. 5F).
Etoposide exerts effective palliation by affecting the subtype of OS and correcting abnormal pathways
Patients with OS and higher-, than lower-risk scores for the immunogenic subtype survived longer (Fig. 4E) and metastasis was less likely to develop (Fig. 6A). Furthermore, the Huvos system associated the immunogenic subtype with higher grades and sensitivity to neoadjuvant chemotherapy (Fig. 6B). We applied the immunogenic subtype gene signature to the DrugMatrix dataset that contains comprehensive information about gene expression in rats under treatment with various drugs. Transcriptional changes caused by thioguanine and etoposide were enriched (Fig. 6C). An extensive literature review revealed that etoposide is a feasible candidate as salvage therapy for relapsed and metastatic OS. A high-throughput transcriptomic screen of OS cells (U-2 OS) exposed to seven concentrations of etoposide was evaluated by GSEA (GSE200845) to determine the underlying mechanism of etoposide. The GSEA enrichment plot associated DEGs with the activation of immunogenic subtype signatures (p = 0.000, NES = 1,80; Fig. 6D left). The heat map also shows that the top 10 genes with immunogenic subtype signatures gradually increased as the dose increased (Fig. 6D right). Meanwhile, the enrichment plot and a heat map showed that the adipocyte-like subtype signatures correlated negatively with the dose-dependent profile of etoposide treatment (p = 0.015; NES = − 1.49; Fig. 6E). We also used GSVA to determine a dose-dependent tendency between immunogenic and adipocyte-like subtype signatures. The results were the same as those of the GSEA analysis (p < 0.01, Fig. 6F). We identified the high-risk genes checkpoint kinase 2 (CHEK2) and retinoblastoma transcriptional corepressor 1 (RB1), which represented 46.15% and 15.39% of the OS driver genes respectively, and play crucial roles in the cell cycle pathway, which we confirmed was activated (Figs. 6G; 3D). The mutation rate of the important tumor suppressor gene tumor protein 53 (TP53) was also 15.39% in OS (Fig. 6G). Etoposide can inhibit the cell cycle pathway, and activate the P53 signaling pathway (Fig. 6H and I). The immunogenic OS subtype originated in MSCs, because it correlated closely with Antigen processing and presentation as described in Fig. 3E (Fig. 6J). Therefore, etoposide has promise as a palliative therapeutic by promoting immunogenic subtypes, inhibiting adipocyte-like subtypes, and correcting the abnormal cell cycle and p53 signaling pathways.
Osteosarcoma is the most frequent primary malignant sarcoma in adolescents. Given the treatment dilemma, new insights into the OS development and therapeutic strategies are urgently needed [25, 26]. Here, we discovered using integrated single-cell analysis, that OS arises from osteogenesis imperfecta during development. We identified deficient TGFβ and P53 signal pathways, an activated cell cycle pathway, and a potentially novel driver gene IFITM5 mutant as potential pathogenic factors of OS. We characterized OS heterogeneity and found that the immunogenic and adipocyte-like subtypes respectively exerted beneficial and detrimental effects on OS. Etoposide is a promising chemotherapeutic drug that achieves palliation by affecting the OS subtypes and correcting aberrant pathways.
Osteosarcoma is derived from MSCs and is characterized by osteogenesis imperfecta [27, 28]. The IFITM5 gene encodes an osteoblast-specific membrane protein that is an established positive regulatory factor during bone mineralization [29, 30]. Single base variants (c.119 C > T), c.143A > G and (c.-14 C > T) in the coding region of IFITM5 have been identified in patients with osteogenesis imperfecta type V [31, 32]. Furthermore, IFITM5 is overexpressed in abnormal bone hyperplasia in rat primary osteoblasts, UMR106 rat osteosarcoma cells, human primary osteoblasts and Saos-2 human osteosarcoma cells [33,34,35]. There are already proven cases published about osteosarcoma occurring in osteogenesis imperfecta due to IFITM5 mutation [36, 37]. The skeletal disorders caused by alterations in the IFITM5 gene, including c.-14C > T, c.119C > T (p.S40L), c.143A > G (p.N48S), while more case series are required to establish detailed genotype-phenotypes for these alterations in the IFITM5 gene [32, 38]. Here, we determined that IFITM is a potential new driver gene in OS development depending on pseudotime using trajectory analysis and the ClinVar database.
The various abnormal signal pathways that prompt MSCs towards abnormal bone growth also play indispensable roles in OS development [6, 25]. The dysregulated cell cycle (CHEK2, RB1 mutant) and p53 (TP53 deficiency) signaling pathways often result in pathogenesis and aberrant OS growth. The p53 signal pathway acts as a tumor suppressor and exerts a crucial role in safeguarding our body from developing OS . There are many studies reported that human osteosarcomas can have a deletion, mutations, and/or rearrangements of the p53 gene, which may cause loss of normal constraints on cellular growth, cell cycle, senescence and metabolism [40,41,42].P53 also regulates osteogenic, chondrogenic, myogenic, adipogenic differentiation of MSCs [43, 44]. The present findings are also in line with the fact that Rb and p53 knockout in bone-marrow-derived MSCs (BM-MSCs) results in OS-like tumors that can be serially transplanted [45, 46]. The TGF-β signal favors normal bone formation in the mesenchymal osteoblastic lineage [47, 48]. However, TGF-β seems to mainly have a pro-tumoral effect on OS [49,50,51,52]. Several clinical strategies to block TGF-β signaling pathways, such as neutralizing antibodies (GC-1008), ligands or receptors inhibitors (AP12009), or chemical compounds (SB-431542, or LY2157299), have not been successful against OS [10, 53,54,55,56]. In contrast to previous findings, we have the new perspective that a TGF-β deficiency is a major cause of OS and leads to a detrimental adipocyte-like subtype. This is an alternative experimental approach to treating patients with OS. Briefly, the TGF-β signal pathway in MSCs guards bone formation in the early stage and then compromises their osteogenic differentiation in the late stage, therefore it dictates the conditions for bone development [47, 57]. In this context, TGF-β signal pathway mutations and/or alterations mutations to TGF-β cascade components have been associated with several bone disorders and many carcinomas [58, 59].Several high-throughput genomic and transcriptomic studies have delineated the intra-tumoral heterogeneity of OS [11, 60,61,62]. Zhang et.al identified seven OS tumor cell clusters with three differentiation branches by using single-cell RNA sequencing of conventional OS and cancellous bone (CB) samples, which have different prognoses and possible drug sensitivities . Jiang et.al classified OS into four subtypes according to the genomic, epigenomic, and transcriptomic data, while our classification depended on the single cell sequencing. We used different methods to identify the heterogeneity of OS, but partially with similar molecular features and clinical prognosis, such as Immune activated (S-IA) and our immunogenic subtype both have immunological traits and better clinical prognosis. Their Immune suppressed (S-IS) subtype has activated adipogenesis and fatty acid metabolism-related pathways and encodes the fatty acid scavenger receptor CD36 which is in line with our adipocyte-like subtype . We explored the heterogeneity and the underlying mechanism of OS that could help to provide better, customized therapy. We identified six oncocyte subtypes of OS, of which four did not affect prognosis (Additional file 1: Fig. S1C and D). The immunogenic and adipocyte-like subtypes were respectively associated with a better and worse prognosis. To our knowledge, this is a novel molecular signature classification of OS, which will facilitate OS assessment, and choosing the optimal therapy. For example, etoposide is associated with the activation and inhibition, respectively, of the immunogenic and adipocyte-like subtype signatures, and combining it with other strategies has shown promising anti-tumour activity in clinical trials [65,66,67,68].
However, it should be noted that the study we examined is limited to bioinformatics analysis, and a significant number of challenges still need to be overcome for a successful in vitro and vivo study. Large-scale samples including both single-cell sequencing and clinical outcomes should be taken into account in the future. In addition, Some different signaling networks in osteosarcoma, including RANKL/RANK, Wnt, Notch, PI3K/Akt/mTOR, and mechanotransduction pathways, contribute to osteosarcoma progression and metastasis and tumor heterogeneity, which certainly looks worthy of further investigation.
We attributed the occurrence of OS to an IFITM5 mutant, deficient TGFβ signaling, continuous activation of the cell cycle signal pathway, and P53 signal inhibition during development. Among six OS subtypes, the immunogenic subtype was beneficial, whereas the adipocyte-like subtype was detrimental to health. Etoposide could be a useful palliative because it altered the OS subtype, inhibited the cell cycle pathway, and improved the p53 signal pathway.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information files.
Alternative telomere lengthening
Antigen presenting cell
B-cell maturation antigen (endothelial cells)
Bone marrow mesenchymal stem cells
Bone morphogenetic protein 2
Breast cancer gene 2
Cluster of differentiation
Carcinoembryonic antigen-related cell adhesion molecule-1
Checkpoint kinase 2
Copy number variation
Collagen type 1 alpha 1
Collagen type 1 alpha 2
COSine similarity-based marker Gene identification
Gene Expression Omnibus
Gene set enrichment analysis
Gene set variation analysis
Human mesenchymal stem cells
Interferon induced transmembrane protein 5
Mammalian target of rapamycin
Myogenic differentiation 1
Nuclear factor kappa B
Neurogenic locus notch homolog protein 1
Protein Kinase CAMP-Dependent Type I Regulatory Subunit Alpha
Protein patched homolog 1
Receptor activator of nuclear factor kappa beta
Retinoblastoma transcriptional corepressor 1
RUNX Family Transcription Factor 2
SRY-Box Transcription Factor 9
Transforming growth factor beta
Transcriptional repressor 53
Uniform Manifold Approximation and Projection
Von Willebrand factor
Weighted Gene Co-expression Network Analysis
Wingless-related integration site
Xie L, Ji T, Guo W. Anti-angiogenesis target therapy for advanced osteosarcoma. Oncol Rep. 2017;38:625–36.
Yu D, Zhang S, Feng A, Xu D, Zhu Q, Mao Y, Zhao Y, Lv Y, Han C, Liu R. Methotrexate, doxorubicin, and cisplatinum regimen is still the preferred option for osteosarcoma chemotherapy: a meta-analysis and clinical observation. Medicine. 2019;98:e15582.
Gorlick R, Anderson P, Andrulis I, Arndt C, Beardsley GP, Bernstein M, Bridge J, Cheung N-K, Dome JS, Ebb D. Biology of childhood osteogenic sarcoma and potential targets for therapeutic development: meeting summary. Clin Cancer Res. 2003;9:5442–53.
Walkley CR, Qudsi R, Sankaran VG, Perry JA, Gostissa M, Roth SI, Rodda SJ, Snay E, Dunning P, Fahey FH. Conditional mouse osteosarcoma, dependent on p53 loss and potentiated by loss of Rb, mimics the human disease. Genes Dev. 2008;22:1662–76.
Czarnecka AM, Synoradzki K, Firlej W, Bartnik E, Sobczuk P, Fiedorowicz M, Grieb P, Rutkowski P. Molecular biology of osteosarcoma. Cancers. 2020;12:2130.
Saraf AJ, Fenger JM, Roberts RD. Osteosarcoma: accelerating progress makes for a hopeful future. Front Oncol. 2018;8:4.
Cortini M, Avnet S, Baldini N. Mesenchymal stroma: role in osteosarcoma progression. Cancer Lett. 2017;405:90–9.
Roberts AB, Wakefield LM. The two faces of transforming growth factor β in carcinogenesis. Proc Natl Acad Sci. 2003;100:8621–3.
Principe DR, Doll JA, Bauer J, Jung B, Munshi HG, Bartholin L, Pasche B, Lee C, Grippo PJ. TGF-β: duality of function between tumor prevention and carcinogenesis. JNCI J Natl Cancer Inst. 2014;106:djt369.
Colak S, Ten Dijke P. Targeting TGF-β signaling in cancer. Trends Cancer. 2017;3:56–71.
Schiavone K, Garnier D, Heymann M-F, Heymann D. The heterogeneity of osteosarcoma: the role played by cancer stem cells. Stem Cells Heterogeneity Cancer. 2019;1139:187–200.
Mutsaers AJ, Walkley CR. Cells of origin in osteosarcoma: mesenchymal stem cells or osteoblast committed cells? Bone. 2014;62:56–63.
Mohseny AB, Szuhai K, Romeo S, Buddingh EP, Briaire-de Bruijn I, de Jong D, van Pel M, Cleton-Jansen AM, Hogendoorn PC. Osteosarcoma originates from mesenchymal stem cells in consequence of aneuploidization and genomic loss of Cdkn2. J Pathol. 2009;219:294–305.
Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, Rotem A, Rodman C, Lian C, Murphy G. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352:189–96.
Mereu E, Lafzi A, Moutinho C, Ziegenhain C, McCarthy DJ, Álvarez-Varela A, Batlle E, Grün D, Lau JK, Boutet SC. Benchmarking single-cell RNA-sequencing protocols for cell atlas projects. Nat Biotechnol. 2020;38:747–55.
Dai M, Pei X, Wang X-J. Accurate and fast cell marker gene identification with COSG. Briefings Bioinform. 2022;23:bbab579.
Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn JL. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32:381–6.
Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, Trapnell C. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14:979–82.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50.
Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J Biomed Mater Res, Part A. 2021;109:1004–14.
Pérez-Palma E, Gramm M, Nürnberg P, May P, Lal D. Simple ClinVar: an interactive web server to explore and retrieve gene and disease variants aggregated in ClinVar database. Nucleic Acids Res. 2019;47:W99–105.
Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44:W90-97.
Xie Z, Bailey A, Kuleshov MV, Clarke DJB, Evangelista JE, Jenkins SL, Lachmann A, Wojciechowicz ML, Kropiwnicki E, Jagodnik KM, et al. Gene set knowledge discovery with Enrichr. Curr Protoc. 2021;1: e90.
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.
Rathore R, Van Tine BA. Pathogenesis and current treatment of osteosarcoma: perspectives for future therapies. J Clin Med. 2021;10:1182.
Raymond AK, Jaffe N. Osteosarcoma multidisciplinary approach to the management from the pathologist’s perspective. Pediatr Adolesc Osteosarcoma. 2009;152:63–84.
Broadhead ML, Sivaji S, Balogh Z, Choong PF. Osteosarcoma: from molecular biology to mesenchymal stem cells. Osteosarcoma Biol Behav Mech 2017.
Abarrategi A, Tornin J, Martinez-Cruzado L, Hamilton A, Martinez-Campos E, Rodrigo JP, González M, Baldini N, Garcia-Castro J, Rodriguez R. Osteosarcoma: cells-of-origin, cancer stem cells, and targeted therapies. Stem Cells Int. 2016;2016:3631764.
Moffatt P, Salois P, Gaumond MH, St-Amant N, Godin É, Lanctôt C. Engineered viruses to select genes encoding secreted and membrane-bound proteins in mammalian cells. Nucleic Acids Res. 2002;30:4285–94.
Hanagata N, Takemura T, Monkawa A, Ikoma T, Tanaka J. Phenotype and gene expression pattern of osteoblast-like cells cultured on polystyrene and hydroxyapatite with pre-adsorbed type-I collagen. J Biomed Mater Res Part A. 2007;83:362–71.
Hanagata N, Li X, Morita H, Takemura T, Li J, Minowa T. Characterization of the osteoblast-specific transmembrane protein IFITM5 and analysis of IFITM5-deficient mice. J Bone Miner Metab. 2011;29:279–90.
Pachajoa H, Giraldo-Ocampo S. A patient with bone fragility, multiple fractures, osteosarcoma, and the variant c. 143A> G in the IFITM5 gene: a case report. Orthop Res Rev. 2022;14:453–8.
Liu BY, Lu YQ, Han F, Wang Y, Mo XK, Han JX. Effects of the overexpression of IFITM5 and IFITM5 c.-14C> T mutation on human osteosarcoma cells. Oncol Lett. 2017;13:111–8.
Moffatt P, Gaumond MH, Salois P, Sellin K, Bessette MC, Godin É, de Oliveira PT, Atkins GJ, Nanci A, Thomas G. Bril: a novel bone-specific modulator of mineralization. J Bone Miner Res. 2008;23:1497–508.
Gorski JP, Wang A, Lovitch D, Law D, Powell K, Midura RJ. Extracellular bone acidic glycoprotein-75 defines condensed mesenchyme regions to be mineralized and localizes with bone sialoprotein during intramembranous bone formation. J Biol Chem. 2004;279:25455–63.
Mathew SE, Santhanam M, Madhuri V. Interferon-induced transmembrane protein 5 mutation causing type-V osteogenesis imperfecta: a case report. JBJS Case Connector. 2015;5: e15.
Maiya S, Grimer R, Ramaswamy R, Deshmukh N. Osteosarcoma occurring in osteogenesis imperfecta tarda. Int Orthop. 2002;26:126–8.
Zhytnik L, Maasalu K, Duy BH, Pashenko A, Khmyzov S, Reimann E, Prans E, Kõks S, Märtson A. IFITM5 pathogenic variant causes osteogenesis imperfecta V with various phenotype severity in Ukrainian and Vietnamese patients. Hum Genomics. 2019;13:1–11.
Velletri T, Xie N, Wang Y, Huang Y, Yang Q, Chen X, Chen Q, Shou P, Gan Y, Cao G. P53 functional abnormality in mesenchymal stem cells promotes osteosarcoma development. Cell Death Dis. 2016;7:e2015–e2015.
Fu H-L, Shao L, Wang Q, Jia T, Li M, Yang D-P. A systematic review of p53 as a biomarker of survival in patients with osteosarcoma. Tumor Biology. 2013;34:3817–21.
Gokgoz N, Wunder JS, Mousses S, Eskandarian S, Bell RS, Andrulis IL. Comparison of p53 mutations in patients with localized osteosarcoma and metastatic osteosarcoma. Cancer. 2001;92:2181–9.
Miller CW, Aslo A, Tsay C, Slamon D, Ishizaki K, Toguchida J, Yamamuro T, Lampkin B, Koeffler HP. Frequency and structure of p53 rearrangements in human osteosarcoma. Can Res. 1990;50:7950–4.
Molchadsky A, Shats I, Goldfinger N, Pevsner-Fischer M, Olson M, Rinon A, Tzahor E, Lozano G, Zipori D, Sarig R. p53 plays a role in mesenchymal differentiation programs, in a cell fate dependent manner. PLoS ONE. 2008;3: e3707.
Molchadsky A, Rivlin N, Brosh R, Rotter V, Sarig R. p53 is balancing development, differentiation and de-differentiation to assure cancer prevention. Carcinogenesis. 2010;31:1501–8.
Cheng L, Wang C, Jing J. Cell cycle kinases in osteosarcoma: potential for therapeutic intervention. Curr Pharm Des. 2016;22:4830–4.
Rubio R, Gutierrez-Aranda I, Sáez-Castillo A, Labarga A, Rosu-Myles M, González-García S, Toribio ML, Menendez P, Rodriguez R. The differentiation stage of p53-Rb-deficient bone marrow mesenchymal stem cells imposes the phenotype of in vivo sarcoma development. Oncogene. 2013;32:4970–80.
Janssens K, Ten Dijke P, Janssens S, Van Hul W. Transforming growth factor-β1 to the bone. Endocr Rev. 2005;26:743–74.
Juárez P, Guise TA. TGF-β in cancer and bone: implications for treatment of bone metastases. Bone. 2011;48:23–9.
Lamora A, Talbot J, Mullard M, Brounais-Le Royer B, Redini F, Verrecchia F. TGF-β signaling in bone remodeling and osteosarcoma progression. J Clin Med. 2016;5:96.
Lamora A, Talbot J, Bougras G, Amiaud J, Leduc M, Chesneau J, Taurelle J, Stresing V, Le Deley MC, Heymann MF. Overexpression of smad7 blocks primary tumor growth and lung metastasis development in osteosarcoma. Clin Cancer Res. 2014;20:5097–112.
Yang R-S, Wu C-T, Lin K-H, Hong R-L, Liu T-K, Lin K-S. Relation between histological intensity of transforming growth factor-β isoforms in human osteosarcoma and the rate of lung metastasis. Tohoku J Exp Med. 1998;184:133–42.
Xu S, Yang S, Sun G, Huang W, Zhang Y. Transforming growth factor-beta polymorphisms and serum level in the development of osteosarcoma. DNA Cell Biol. 2014;33:802–6.
Mohammad KS, Javelaud D, Fournier PG, Niewolna M, McKenna CR, Peng XH, Duong V, Dunn LK, Mauviel A, Guise TA. TGF-β-RI kinase inhibitor SD-208 reduces the development and progression of melanoma bone metastasesTGF-β inhibition to treat melanoma bone metastases. Can Res. 2011;71:175–84.
Hau P, Jachimczak P, Schlingensiepen R, Schulmeyer F, Jauch T, Steinbrecher A, Brawanski A, Proescholdt M, Schlaier J, Buchroithner J. Inhibition of TGF-β 2 with ap 12009 in recurrent malignant gliomas: from preclinical to phase I/II studies. Oligonucleotides. 2007;17:201–12.
Morris JC, Tan AR, Olencki TE, Shapiro GI, Dezube BJ, Reiss M, Hsu FJ, Berzofsky JA, Lawrence DP. Phase I study of GC1008 (fresolimumab): a human anti-transforming growth factor-beta (TGFβ) monoclonal antibody in patients with advanced malignant melanoma or renal cell carcinoma. PLoS ONE. 2014;9: e90353.
Rodon Ahnert J, Baselga J, Calvo E, Seoane J, Brana I, Sicart E, Gueorguieva I, Cleverly A, Lahn M, Pillay S. First human dose (FHD) study of the oral transforming growth factor-beta receptor I kinase inhibitor LY2157299 in patients with treatment-refractory malignant glioma. J Clin Oncol. 2011;29:3011–3011.
Maeda S, Hayashi M, Komiya S, Imamura T, Miyazono K. Endogenous TGF-β signaling suppresses maturation of osteoblastic mesenchymal cells. EMBO J. 2004;23:552–63.
Alliston T, Piek E, Derynck R. TGF-beta family signaling in skeletal development, maintenance, and disease. Cold Spring Harbor Monogr Ser. 2008;50:667.
Levy L, Hill CS. Alterations in components of the TGF-β superfamily signaling pathways in human cancer. Cytokine Growth Factor Rev. 2006;17:41–58.
Martin JW, Squire JA, Zielenska M. The genetics of osteosarcoma. Sarcoma. 2012;2012:627254.
Chen X, Bahrami A, Pappo A, Easton J, Dalton J, Hedlund E, Ellison D, Shurtleff S, Wu G, Wei L. Recurrent somatic structural variations contribute to tumorigenesis in pediatric osteosarcoma. Cell Rep. 2014;7:104–12.
Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, Wang Y, Zhang Z, Yuan T, Ding X. Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma. Nat Commun. 2020;11:1–17.
Adamopoulos C, Gargalionis AN, Basdra EK, Papavassiliou AG. Deciphering signaling networks in osteosarcoma pathobiology. Exp Biol Med. 2016;241:1296–305.
Jiang Y, Wang J, Sun M, Zuo D, Wang H, Shen J, Jiang W, Mu H, Ma X, Yin F. Multi-omics analysis identifies osteosarcoma subtypes with distinct prognosis indicating stratified treatment. Nat Commun. 2022;13:1–17.
Gaspar N, Venkatramani R, Hecker-Nolting S, Melcon SG, Locatelli F, Bautista F, Longhi A, Lervat C, Entz-Werle N, Casanova M. Lenvatinib with etoposide plus ifosfamide in patients with refractory or relapsed osteosarcoma (ITCC-050): A multicentre, open-label, multicohort, phase 1/2 study. Lancet Oncol. 2021;22:1312–21.
Akazawa R, Umeda K, Saida S, Kato I, Hiramatsu H, Sakamoto A, Arakawa Y, Sumiyoshi S, Okamoto T, Moritake H. Temozolomide and etoposide combination for the treatment of relapsed osteosarcoma. Jpn J Clin Oncol. 2020;50:948–52.
Gaspar N, Occean B-V, Pacquement H, Bompas E, Bouvier C, Brisse HJ, Castex M-P, Cheurfa N, Corradini N, Delaye J. Results of methotrexate-etoposide-ifosfamide based regimen (M-EI) in osteosarcoma patients included in the French OS2006/sarcome-09 study. Eur J Cancer. 2018;88:57–66.
Perret A, Dômont J, Chamseddine AN, Dumont SN, Verret B, Briand S, Court C, Lazure T, Adam J, Ngo C. Efficacy and safety of oral metronomic etoposide in adult patients with metastatic osteosarcoma. Cancer Med. 2021;10:230–6.
We are very grateful to the participating universities and China Scholarship Council.
This work was supported by the National Natural Science Foundation of China (81702667), Natural Science Foundation of Shandong (ZR202111050005), China Postdoctoral Science Foundation (8206300728). This work also funded by China Scholarship Council.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Fig.S1.
A The correlation analysis between IFITM5 expression and each tumor subtype. B Markers used to define subpopulations. C Box plots show the correlation between metastasis and other subtype scores. D Huvos grades of neuron-like, fibroblastic, chondroblastic, and osteoblastic subtypes.
Additional file 2: Table S1.
Different genes between normal osteocytes and OS oncocytes.
Additional file 3: Table S2.
Signatures of each module in WGCNA analysis.
Additional file 4: Table S3.
Signatures of each OS subtype.
Additional file 5: Table S4.
Significantly different genes of the hMSC_BMP2+IBMX group.
Additional file 6: Table S5.
Significantly different genes of the hMSC_BMP2+IBMX+TGFβ group.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Sun, Y., Zhang, C., Fang, Q. et al. Abnormal signal pathways and tumor heterogeneity in osteosarcoma. J Transl Med 21, 99 (2023). https://doi.org/10.1186/s12967-023-03961-7
- Abnormal signal pathway
- Oncocyte heterogeneity