RNA expression data and clinical information of osteosarcoma patients were downloaded from TARGET (https://ocg.cancer.gov/programs/target) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases. The following were the inclusion criteria: (a)osteosarcoma samples with gene expression matrix; (b)samples with clinical information such as age, gender, survival time, survival status, and whether metastasis occurred; (c)samples with expression values for more than half of the genes. Based on the above criteria, 4 eligible osteosarcoma cohorts (GSE21257, GSE39055, GSE16102, and TARGET-OS)were collected for further analysis. The batch effect of non-biotechnical bias was corrected using the “ComBat” package. The cohorts GSE21257, GSE39055 and GSE16102 were merged into the meta-cohort. 331 cancer stem cells-associated genes (CSCRGs) were obtained from the molecular marker database, of which 228 genes were expressed in the TARGET-OS cohort and meta-cohort. In the TARGET-OS cohort, we selected 25 CSCRGs for further studies using univariate COX analysis (P<0.05) ( Supplementary Tables 1 and 2 ).

The consensus clustering was performed using the “ConsensusClusterPlus” package based on the expression matrix of the 25 CSCRGs (39). Differentially expressed genes (DEGs) between three clusters were analyzed using ‘limma’ package with the cutoff criteria of P < 0.05.

The Gene Ontology (GO) enrichment was performed using the “clusterprofiler” package. The geneset “h.all.v7.5.1.symbols” was obtained from MSigDB. Mariathasan et al. had constructed a geneset in which genes related to certain biological processes were stored (40). Based on the above gene sets, we performed Gene set variation analysis (GSVA) by using the “GSVA” and “limma” packages to show alterations in signaling pathways among three clusters (41). Adjusted P value less than 0.05 was defined as statistically significant. The ESTIMATE algorithm was used to calculate the immune score, stromal score and tumor purity.

Cohorts of 23 immune infiltrating cells and 13 immune-related functions were obtained (42). And the score calculated with ssGSEA was used to express the relative abundance of different immunocytes and immune-related functions in every case (43). The abundance of six immunocytes were analyzed using the TIMER algorithm.

In order to quantify the stemness of individual tumors, we constructed a score system called CSC score, which was constructed in the following steps. We selected overlapping DEGs found in distinct CSC clusters and performed prognostic analyses of individual genes by using univariate Cox regression. Extraction of genes with remarkable prognosis was used to construct the CSC score by principal component analysis (PCA). Principal component1 and Principal component2 were used in the construction of the CSC score. We also defined CSC score using a similar approach from previous studies (44, 45). CSC score =ΣPC1_i + PC2_i, i is the expression of genes associated with the CSC phenotype. We stratified the tumors into CSC score low and high subgroups using the surv-cutpoint function in the ‘survival’ package.

We used the ‘pRRophetic’ package to assess the sensitivity to different CSC clusters to small molecule drugs. In addition, the CellMiner database was utilized to assess the association between CSCRGs and different drugs (46).

Based on one-class logistic regression (OCLR) algorithm, the stemness index model trained from the Progenitor Cell Biology Consortium database was used to calculate tumor stemness (47, 48). The stemness index can be used to measure how similar tumor cells are to stem cells, with stemness index being a value between 0 (lowest) and 1 (highest). The closer the stemness index is to 1, the stronger the stem cell properties. We calculated transcriptome feature scores for the cohorts using the same Spearman correlation.

A total of 10 OS tissues were collected, all from Peking University People’s Hospital. The samples were examined by three experienced pathologists. All patients provided informed consent, and the study protocol was approved by the Ethics Commttee of Peking University People’ s Hospital (2019PHB198-01). Immunohistochemical examination was performed with MEF2C antibody (10056-1-AP, proteintech).

Human osteosarcoma cell line (143B) was purchased from the American Type Culture Collection (ATCC). 143B cells was cultured in DMEM (Gibco) containing 10% fetal bovine serum (FBS, Gibco) and 1% penicillin and streptomycin (Gibco). 143B cells were cultured in a humidified incubator with 5% CO2 at 37°C. Si-MEF2C (Suzhou, China, sequences: 5’GACAAGGAAUGGGAGGAUA3’) and GP-transfect-Mate (Suzhou, China) were used for transfection.

143B cells were inoculated at a density of 1000 cells/well in six-well ultra-low attachment plates (Corning). And the cells were cultured in DMEM/F12 medium (Gibco) containing N2 medium (Invitrogen), human EGF (20ng/ml, PeproTech) and human bFGF (20ng/ml, PeproTech) for 14 days. Spheres were observed in size and the number of spheres formed was calculated.

50 µl of vitronectin (PeproTech) or fibronectin (Biocoat) was added to each well of a 96-well plate and incubated overnight at 4°C. Unbound proteins were washed with PBS and closed with PBS containing 2% BSA for 2 hours at 37°C. 143B cells were inoculated at a density of 10000 cells/well in 96-well plate (Corning). And the cells were cultured in DMEM (Gibco) for 1 hour. Unbound cells were washed with PBS, fixed with 4% paraformaldehyde and stained with 0.1% crystal violet. Finally, the number of adherent cells was observed under the microscope.

Statistical analyses were performed via R (version 4.1.2), and survival analyses were performed using the Kaplan-Meier method. The Student’s t test was used for normally distributed variables and the Wilcoxon rank sum test was used for non-normally distributed variables. The Kruskal-Wallis test and one-way ANOVA were used for the non-parametric and parametric methods, respectively (49). Correlations coefficients between the expression of CSCRGs and the TME infiltrating immunocytes were calculated by Spearman analysis. We used univariate Cox regression analysis to compute the hazard ratio (HR) of CSCRGs and CSC-related signature genes. To verify whether the constructed risk scores can be used as an independent prognostic factor independently of other clinical traits. Univariate and multivariate COX analyses of patients’ age, gender, presence of metastasis, and CSC score were performed. The predictive performance of the CSC score was assessed by using receiver operating characteristic (ROC) curves, and the area under the curve (AUC) was calculated by the ‘timeROC’ package. P-value was bilateral and P < 0.05 was defined as statistically significant difference.

In this study, we finally selected 25 cancer stem cells related genes (CSCRGs) by univariate COX analysis in TARGET-OS cohort and later investigated the role of these genes in osteosarcoma. Metascape analysis and GO enrichment analysis were performed on 25 CSCRGs, which were seen to be enriched in multiple stemness-related regulatory pathways, and the results were shown in Figure 1A and Figure S1A . In addition, Spearman analysis was used to assess the relevance between 25 CSCRGs ( Figure S1B ). COX regression analysis was used to analyze the relationship of 25 CSCRGs with patient prognosis in osteosarcoma. The forestplot showed that FOLR1, SEMA3B, SEMA4G, MYC, OVOL1 and MEF2C were considered as risk factors ( Figures S1C, D ). The above analyses showed 25 CSCRGs played important roles in the development and progression of osteosarcoma.

Data of 94 osteosarcoma patients in the TARGET-OS cohort were used for analysis. The network showed the interactions of the 25 CSCRGs and their prognostic significance for osteosarcoma patient ( Figure 1B ). The results indicated that these 25 CSCRGs were mutually regulated and played crucial roles in the development of osteosarcoma. We utilized consensus clustering algorithm to stratify samples into distinct CSC clusters based on the expression of the 25 CSCRGs. Consequently, we identified three distinct clusters, including 28 samples in cluster A, 39 samples in cluster B and 27 samples in cluster C ( Figures S2A, B , Supplementary Tables 3, 4). Prognostic analysis for these three clusters indicated that CSC cluster A showed an excellent survival advantage, while CSC cluster C had the worst prognosis in TARGET-OS cohort ( Figures 1C, D ). We also noted that there were remarkable differences in expression levels of the 25 CSCRGs between distinct clusters ( Figures S2C, D ). ANXA6, ENG, EVI2B, SEMA3E, and SPI1 were significantly elevated in CSC cluster A, CDK6, ETV6, MYC, PHACTR4, PTN, SEMA3A, and SOX4 were evidently increased in CSC cluster B, and SEMA4G was evidently increased in CSC cluster C ( Figures S2C, D ).

To explore the biological behaviors underlying these different CSC clusters, we performed GSVA enrichment analysis ( Figures 1E, F ). The results showed that CSC cluster A was markedly abundant in immune activation-related processes, such as complement, inflammatory and interferon gamma response, and the TNFA pathway, KARS pathway and apoptosis pathway were also enriched in CSC cluster A. Based on the above findings, we presumed that the better prognosis of CSC cluster A might be related to its high immune infiltration. Furthermore, we quantified the stromal score, immune score and tumor purity for the three clusters using ESTIMATE algorithm. The analysis showed that CSC cluster A had the highest immune score and stromal score, followed by CSC cluster B and C ( Figures 2A, B ). Conversely, CSC cluster B and C had higher tumor purity compared to CSC cluster A, suggesting that tumors in CSC cluster A were surrounded by more non-tumor components (immunocytes and stromal cells) ( Figure 2C ). In addition, a heat map was built by ssGSEA to visualize and compare the abundance of 23 immunocytes under different clusters ( Figure 2D ). The great majority of immunocytes such as anti-tumor lymphocyte cell subpopulations and NK T cells were mainly enriched in the CSC cluster A. We further described the immune infiltration profile using TIMER2.0 and observed consistent results ( Figures 2E–H ). The results showed that neutrophils, macrophages and myeloid dendritic cells were mainly enriched in the CSC cluster A, while CD8+ T cells were enriched in CSC cluster A and C. We also calculated 13 immune-related function scores using ssGSEA, and the results suggested that majority of immune-related functions were enriched in CSC cluster A ( Figure 2I ). Curiously, CSC cluster C also had a higher infiltration of immune cells but did not show the same survival advantage. DCs have been shown in previous studies to be responsible for antigen presentation and initial T cell activation, bridging innate and adaptive immunity, and their activation is dependent on high expression levels of MHC molecules, co-stimulatory factors and adhesion factors. Therefore, we compared the expression of MHC molecules, co-stimulatory molecules and adhesion molecules in the three CSC subtypes and found that most molecules including CD40, CD80, CD86, HLA-C, HLA-DMC, HLA-DMB, HLA-DPB1, HLA-DRA, HLA-E, and ICAM1 were significantly elevated in CSC cluster A ( Figure S2E ). Thus, we hypothesize that although CSC cluster C had a higher immune infiltration, its antigen-presenting ability and ability to activate DCs were weaker compared to CSC cluster A. Therefore, CSC cluster C had a poorer survival prognosis compared to CSC cluster A. The results from GSVA analysis demonstrated that Pan-F-TBRS, antigen processing machinery, immune checkpoint, and CD8 T effector were enriched in CSC cluster A further corroborating our hypothesi ( Figure S2F ). Taking into account that PD-L1 is a proven biomarker to predict immunotherapy response (37), we identified a significant upregulation of PD-L1 expression levels in CSC cluster A ( Figure S2G ). Based on these findings, we identified three subtypes with distinct immune infiltration characteristics.

In addition, Spearman analysis was used to assess the specific relationship between the 25 CSCRGs and immunocyte infiltration or immune-related function ( Figures 3A, B ). High expression of EVI2B, ENG and SPI1 was markedly associated with enhanced immunocyte infiltration and immune-related function, whereas MEF2C, PHACTR4, MYC, PTN, SEMA3A and SEMA4G high expression showed a negative correlation with the immunocyte infiltration and immune-related function level. Among the 25 CSCRGs, we focused on MEF2C, which was found to be markedly negatively correlated with substantial immune infiltration and immune-related function levels and the expression high of MEF2C was significantly negatively correlated with patient prognosis ( Figure 3C ). We firstly compared the overall level of immune cell infiltration in patients with high and low MEF2C expression. Patients with low MEF2C expression had higher immune scores, indicating that TME immune cell infiltration was markedly increased in patients with low MEF2C expression ( Figure 3D ). We then compared the differences in 23 immunocytes between the two subgroups with low and high MEF2C expression ( Figure 3E ). We found that high MEF2C expression was remarkably negatively associated with the levels of infiltration of multiple immune cells, including regulatory T cells, T follicular helper cells, type 1 T helper cells, macrophages, MDSCs, natural killer cells, activated DCs and CD8+ T cells ( Figure 3E ). DCs are responsible for antigen presentation and initial T cell activation, bridging the gap between innate and adaptive immunity [52]. According to these findings, we hypothesize that MEF2C may inhibit the cytotoxic T lymphocytes and activated DCs, thereby hindering the intratumoral anti-tumor immune response. Based on the above findings, we speculated that the expression of MEF2C could influence the prognosis of patients by affecting the infiltration of multiple immune cells.

To further confirm that the typing based on 25 CSCRGs was also applicable to other datasets, we performed validation with the GSE21257 cohort. Similar to the clustering results of the TARGET-OS cohort, unsupervised clustering also revealed three completely different clusters of the 25 CSCRGs in the GSE21257 cohort ( Supplementary Table 5 , Figures S3A, B ). The mRNA expression of 25 CSCRGs was significantly different in the three cluster ( Figure 4A ). Prognostic analysis indicated that the survival of CSC cluster A was better than that in CSC cluster B and C ( Figure S3C ). We also quantified the stromal score, immune score and tumor purity for the three clusters using ESTIMATE algorithm. The results showed that CSC cluster A had the highest immune score compared to CSC clusters B and C ( Figure 4B ), and CSC clusters A and B had a higher stromal score compared to CSC cluster C ( Figure 4C ). And consistent with the previous result, CSC cluster A had lower tumor purity ( Figure 4D ). We then analyzed the differences in immune infiltration between the three clusters, and we found the vast majority of immunocytes such as anti-tumor lymphocyte cell subpopulations and NK T cells were largely enriched in the CSC cluster A ( Figure 4E ). And consistent with previous results, most immune-related functions were enriched in CSC cluster A ( Figure 4F ). Similarly, the results of gsva analysis remained consistent with previous ones, showing enrichment of CD8 T-effects, immune checkpoints and Pan-F-TBRS in the CSC cluster A ( Figure S3D ).

Previously, tumor samples were divided into three subtypes associated with CSC based on 25 CSCRGs, and to further explore the genetic alterations in these phenotypes, we determined 104 CSC-related DEGs using limma package in the TARGET-OS cohort ( Figure 5A , Supplementary Table 6 ). We performed GO enrichment analysis on these DEGs, which were seen to be enriched in multiple stemness-related and immune regulation pathways ( Figure 5B ). Furthermore, the above analysis supported that DEGs were closely associated with tumor immunity and cancer stemness, and thus might be considered as CSCRGs. We performed consensus clustering analysis on the 104 CSC phenotype-related genes obtained to further validate this regulatory mechanism and obtained three gene clusters ( Figures S4A, B ). We called these clusters as CSC gene cluster A-C. These results suggested that three distinct CSC subtypes indeed existed in osteosarcoma ( Figure 5C ). Among the three CSC gene clusters, 25 CSCRGs were found to be significantly differentially expressed, which was consistent with the expected results ( Figure S4C ). Survival analysis further showed prognostic differences between three CSC gene clusters. CSC gene cluster B was shown to be related to better prognosis, while CSC gene cluster A was proven to be related to worse prognosis ( Figures 5D, E ). To explore the roles of the 104 genes in immune infiltration, we examined the differences in 23 TME immune cells in the three clusters and showed that the most immunocytes increased in CSC gene cluster B ( Figure S4D ). Similarly, the result suggested that most immune-related functions were enriched in CSC gene cluster B ( Figure S4E ). Meanwhile, we found that PDL1 expression was significantly upregulated in CSC gene cluster B compared with CSC gene cluster A ( Figure S4F ). The above results once again suggested that CSC-related genes played non-negligible regulatory roles in the formation of different TME landscapes.

While previous studies have found important roles in prognosis and regulation of immune infiltration for CSC-related genes, these findings were only applicable to assess patient populations and could not be used to evaluate individual patients. Taking into account individual differences, we constructed a score system to quantitate CSC-related subtypes in single osteosarcoma patients based on the discovery of CSC-related signature genes, which we called CSC score. We constructed an alluvial diagram to illustrate the workflow of CSC score construction ( Figure 6A ). The result showed that CSC gene cluster A was related to higher CSC scores, and CSC gene cluster B was associated with lower CSC scores ( Figure S5A ). Notably, consistent with the expected results, CSC cluster B and C showed a higher CSC score than CSC cluster A ( Figure S5B ). We examined the relationship between CSC scores and certain biometric scores by Spearman analysis. The result showed that CSC score was negatively related to immune activation-related processes ( Figure 6B ). CSC score was also significantly negatively related to immune score ( Figure S5C ). Similarly, patients in the low CSC score subgroup had a higher degree of immune infiltration and were more enriched for immune-related functions compared to patients in the high CSC score subgroup ( Figures S5D, E ). The above analyses clearly indicated that low CSC score was remarkably correlated with immune infiltration. Based on the above findings, we concluded that the CSC score could better assess the CSC-related subtypes of individual tumors and further assess the tumor immune infiltration characteristics. Furthermore, we attempted to determine the value of CSC score in predicting patient prognosis. Patients were separated into low and high subgroups with a cutoff value of -0.3269, and patients with low CSC score had a better prognosis ( Figures 6C,D ). ROC curve analysis result verified the predictive advantage of the CSC score system ( Figure S6A ). We analyzed multiple clinical traits of patients using multivariate Cox regression and found that the CSC score system could potentially serve as an independent prognostic factor of osteosarcoma ( Figure S6B ). Futhermore, PD-L1 expression level was significantly higher in the group with low CSC score ( Figure S6C ). The constructed CSC scoring system was validated by meta-cohort, and patients with low CSC score indicated better prognosis ( Figure 6E , Supplementary Table 7 ). To further validate the reliability of the CSC score system, we used a cohort (GSE21257) from the meta cohort to explore the association between CSC score and patient prognosis. Consistent with the results above, patients with low CSC score showed a significant survival advantage relative to patients with high score ( Figure 6F ). Similarly, ROC curve analysis result verified the predictive advantage of the CSC score system ( Figure S6D ). The above results strongly indicated that CSC scores could represent the CSC pattern of osteosarcoma patients and predict the prognosis of osteosarcoma patients.

CSC is associated with a variety of clinical traits such as tumor metastasis. We compared the differences of CSC score between different clinical subgroups in the TARGET cohort. Accordingly, we found that tumor metastasis was significantly related to higher CSC score, implying that these patients were characterized by poor immune infiltration, with a poorer clinical outcome ( Figure 7A ). And there were no differences in scores between age and gender subgroups ( Figures S6E, F ). In particular, we found that low CSC score subgroup showed a significant survival advantage both in patients who had developed tumor metastases and in those who did not ( Figures 7B,C ). Additional results showed that the predictive ability of the CSC score was not interfered by whether the tumor had metastasized, and the low CSC score group consistently showed a significant survival advantage in both patients with and without metastasis ( Figures 7B, C ). These results suggest that the CSC score can also be used to assess certain clinical symptoms such as tumor metastasis in osteosarcoma patients. Considering the strong correlation between CSC score and immune response in the above results, we investigated whether CSC score could predict treatment response to anti-PDL1 in patients in an independent immunotherapy cohort. We performed a systematic search that resulted in the inclusion of an immunotherapy cohort: uroepithelial carcinoma intervening with atezolizumab (IMvigor210 cohort). The result indicated that patients with low CSC score scores gained a survival advantage ( Figure 7D ). And there was a significant therapeutic advantage of anti-PDL1 treatment in patients with low CSC scores ( Figure 7E ). The above results suggested that the CSC score might be used to assess the therapeutic response and clinical prognosis of immunotherapy.

Using the OCLR algorithm, the stemness index was computed for individual samples based on the gene expression profile of the patient. We compared the mRNAsi between different CSC clusters and showed that mRNAsi was significantly elevated in CSC cluster C ( Figure 7F ).Similarly, mRNAsi was significantly elevated in CSC gene cluster A ( Figure 7G ). We also found that patients with high MEF2C expression possessed higher mRNAsi ( Figure 7H ). Considering that mRNAsi represents the stemness of individual patients, we analyzed the correlation between CSC score and mRNAsi. The results revealed a significant positive relationship between CSC score and mRNAsi in the TARGET cohort ( Figure 7I ). However, the result showed a weak correlation between CSC score and mRNAsi (R=0.22). Therefore, unlike the previous result that the CSC score of CSC cluster A was significantly lower than cluster B, the mRNAsi of CSC cluster A was not significantly different from cluster B. Similarly, there was a significant positive correlation between CSC score and mRNAsi in both the meta cohort and the GSE21257 cohort ( Figures 7J, K ). The above results again demonstrated that the CSC score was a reliable scoring system that could represent the stemness of individual patients.

We performed a drug sensitivity analysis and identified 74 small molecule drugs that may be used in the treatment of osteosarcoma ( Supplementary Table 8 ). Sensitivity of 12 drugs in different CSC clusters was depicted in Figures 8A-L . The results showed that CSC cluster A was sensitive to BIRB.0796 and OSI.906, while CSC clutser B was sensitive to Methotrexate, Sorafenib, Sunitinib and Nilotinib, and CSC clutser C was more sensitive to Cisplatin, Doxorubicin, Imatinib, CHIR.99021, Pazopanib and Vinorelbine were more sensitive ( Figures 8A-L ). We then evaluated the relationship between 25 CSCRGs expressions and drug sensitivity using the CellMiner database ( Figure S7 , Supplementary Table 9 ). The above findings indicated that exploring CSC subtypes in osteosarcoma patients could be used to guide the clinical use of drugs.

Given that MEF2C was found to affect prognosis, immune infiltration and stemness in osteosarcoma patients in our previous study, we verified the expression of MEF2C in clinical tissues of osteosarcoma patients by immunohistochemical experiments. The results showed that MEF2C was significantly expressed in clinical tissues of osteosarcoma patients ( Figure 9A ).

We used a sphere formation assay to examine the effect of MEF2C on tumor cell stemness ( Figures 9B, C ). The results showed that the si-MEF2C cells formed fewer sphere ( Figure 9B ) and the size of sphere was smaller ( Figure 9C ) compared to the control group. These results can tentatively demonstrate the effect of MEF2C on the maintenance of stemness of CSC in osteosarcoma. Our previous analysis showed that MEF2C affected immune infiltration in the TME. We also verified that MEF2C affects the adhesion of tumor cells to the extracellular matrix by cell adhesion assays ( Figures 9D, E ).The results showed that the adhesion of si-MEF2C cells to the extracellular matrix such as vitronectin ( Figure 9D ) and fibronectin ( Figure 9E ) was reduced compared with the control group.