IDH wild-type glioma has been observed with aggressive phenotypes and poor prognosis compared to the IDH mutant glioma (Agnihotri et al., 2014; Picca et al., 2018). In order to identify the underlying mechanisms that drive the development and progression of IDH wild-type glioma, here we performed scRNA-seq analysis of one surgically resected of glioma tumor without IDH1 and IDH2 mutations. First, we showed that this tumor belongs to a grade III anaplastic astrocytoma based on its histopathological features (Figure 1A). Besides, the immunohistochemical staining results revealed that this tumor shows strong expression of the GFAP, a classic marker of astrocytoma. Moreover, this tumor also exhibited positive staining for the IDH1 (Figure 1A).

By use of the 10x Genomics scRNA-seq platform, a total of 5,942 high quality cells were identified in this tumor with an average of 2,265 genes detected per cell. Overall, five distinct cell clusters could be distinguished (Figure 1B). We observed that most cells are GFAP⁺ astrocytoma cells, as well as relatively low proportions of PTRPC⁺ immune cells and MOG⁺ oligodendrocytes in this tumor (Figure 1C). Moreover, we found that tumor-associated macrophages (TAMs) predominate over T cells and B cells in this tumor (Figure 1D). However, we did not observe other cell types in this tumor, such as the tissue-resident astrocytes, vascular endothelial cells and neurons. In addition, we observed that the expression levels of IDH1, IDH2, EGFR and TP53 are relatively low in tumor cells, whereas the expression levels of CDK4 and PDGFRA are relatively high (Figure 1E).

In this study, we identified a total of 1,765 infiltrating immune cells in this tumor, which could be further divided into five subpopulations, including three TAM, one T cell and one B cell subpopulations (Figure 2A and Supplementary Figure 1A). Notably, TAMs were observed to be predominant in the TME of this tumor, accounting for over 79% of immune cells. Next, we found that T cells are characterized by high expression of genes related to NKT cells, such as NKG7, IL7R and MATK (Figure 2B and Supplementary Figure 1B). Moreover, we found that B cells show high expression of IGKC and MZB1 genes, but low expression of CD79A and CD79B (Figure 2C). In addition, we observed high expression of universal markers of macrophage (C1QA, C1QB, APOE and TREM2) in all three TAM subpopulations (Figure 2D). Although the three clusters of TAMs have been defined by unsupervised clustering analysis, all the TAMs could be defined as CD45^high CX3CR1⁺ CD11b^low F4/80^low CCR2⁻ macrophages (Supplementary Figure 1C). Interestingly, we found that a subset of TAMs also express genes encoding the markers of glioma cells, such as GFAP, DLL3, OLIG1 and SOX4 (Figure 2E). Consistently, immune cells expressing markers of tumor cells have also been reported in prostate cancer (Chen et al., 2021).

Finally, through comparing gene expression profiling of each TAM subpopulation, we observed that some genes, encoding many cytokines and chemokines, are significantly upregulated in the TAM1 and TAM3 subpopulations, indicating that these TAMs are in an activated state (Figure 2F). For example, we observed high expression of inflammatory molecules (e.g., CCL2, CCL3, CCL4, and IL1B) and M2-like polarization markers (e.g., CLEC7A, SGK1, ATF3, CXCL2, CXCL8, and VEGFA) in the TAM1 and TAM3 cells (Figures 2G,H). Collectively, these results indicated that TAMs in this glioma do not fit into a typical M1 or M2 phenotype, instead they represent a mixed M1/M2-like polarized phenotype.

In order to dissect the intra-tumoral heterogeneity of cancer cells in this glioma, we further re-clustered all the cancer cells, which also could be divided into five subpopulations (Figure 3A). Previously, high-grade glioma has been identified to be composed of mixed cancer cell populations, including astrocyte-like (AC-like), neural progenitor-like (NPC-like), oligodendrocyte progenitor-like (OPC-like) and mesenchymal-like (MES-like) cells (Phillips et al., 2006; Neftel et al., 2019). Consistently, we found that cells in cluster 0 and three display oligodendroglia lineage markers, such as OLIG1, BCAN, PLP1 and PLLP, which refers to the OPC-like cells (Figures 3B,C). Notably, these cells also show high expression of NPC-like cell signature, including SOX4 and DLL3 (Figure 3B), reflecting the potential of NPC-like toward OPC-like cell differentiation. Here we named these cell clusters as OPC1 and OPC2. Next, we observed that cells in cluster two show high expression of cell cycle-related genes, including CDK1, CCNB1, CCNB2, and UBE2C, which refers to the cycling cell subpopulation (Figure 3E). Interestingly, we observed that these cycling cells also show features of the OPC- and NPC-like cells (Figure 3B).

Next, we found that astrocytic cell markers are widely expressed in all cancer cells, indicating that these cancer cells show features of AC-like cells (Figure 3D). In addition, we observed that cells in cluster one and four exhibit high expression levels of mesenchymal-related genes (e.g., CD44 and VIM), referring to the MES-like cells (Figure 3F). Notably, we found that cells in cluster four shows an increased expression of genes related to the hypoxia (e.g., HILPDA), tumor angiogenesis (e.g., VEGF and NRN1) and glycolytic (e.g., ENO2) (Figure 3F). Therefore, we defined these two MES-like subpopulations as hypoxia-independent (MES1, cluster 1) and hypoxia-dependent (MES2, cluster 4) cells. Together, our findings suggested that this high-grade glioma tumor contains a mixed population of four major cell states.

In order to trace the dynamic evolution of glioma cells during tumor progression, we performed trajectory analysis of the cancer cells. Overall, we identified four branches (represented as “B1,” “B2,” “B3,” and “B4”) and four states of cancer cells (Figure 4A). We found that the trajectory’s root (B1) is composed by majority of the cycling cells (cluster 2) and some OPC- and NPC-like cells. Moreover, we observed that the B3 and B4 branches are composed of most cells from the MES-like cells (cluster one and 4), indicating that cells in these branches are derived from cells in B1 branch (Figure 4B). Therefore, we could postulate that glioma cells evolve from B1 through branch point 1, which further differentiate into two branches, the B2 (OPC- and NPC-like cells) and the B4 (MES-like cells). Therefore, these results indicated the malignant progression of glioma cells, spanning from the cycling, OPC/NPC-like, eventually to MES-like cells. Furthermore, we displayed the expression levels of representative genes in distinct cell states along the trajectory tree (Figure 4C), and also the expression of representative markers of cycling, OPC/NPC-like and MES-like cells over pseudo-time (Figure 4D). In addition, the differentially expressed genes in cycling, OPC/NPC-like and MES-like cells were plotted (Figure 4E). Together, these findings indicated that a dynamic evolution of cancer cells drives the intra-tumoral heterogeneity of glioma.

In order to identify mechanisms that drive the evolution of glioma cells, we investigated the cell-cell communication mediated by ligand-receptor interactions across all cell components. Accordingly, a total of 512 pairs of interactions were identified in this glioma (Supplementary Table 1). We observed that all cell subpopulations are connected with each other (Figure 5A and Supplementary Figure 2A). First, we focused on exploring the relationship between tumor and stromal cells. Many interactions between cancer cells and other cell types in TME were observed, especially the TAM and T cells (Figure 5B). For example, we found that all cancer cell subpopulations express macrophage migration inhibitory factor (MIF), which could bind to CD74 receptor on the TAM and T cells (Figure 5B and Supplementary Figure 2B, C). Besides, we found increased expression of angiogenic CXC chemokines in the TAM1 macrophages, including CXCL2 and CXCL8 (Supplementary Figure 2C, D). Notably, MIF has been reported to be overexpressed in several cancers, and the interaction between MIF and CD74 has also been shown to promote cell proliferation and tumor growth (Xu et al., 2008; McClelland et al., 2009; Cheng et al., 2011; Zheng et al., 2012).

Furthermore, we investigated the roles of stromal cells, especially the TAMs, in regulating cancer cells. We observed multiple ligands and receptors interactions related to tumor progression. For example, we observed the interaction between growth arrest-specific 6 (GAS6) expressed by TAMs and its cognate receptor TYRO3 expressed on the cycling tumor cells (Figure 5C and Supplementary Figure 2E), which was reported to positively correlate with tumor growth (Wu et al., 2018). Also, we observed OSM/OSMR and OSM/LIFR interaction pairs between the TAMs and MES-like cells (Figure 5C and Supplementary Figure 2F), which was consistent with previous findings in glioma (Hara et al., 2021). Notably, we found that SPP1/CD44 interaction pair stands out among all interaction pairs, which primarily mediates the crosstalk between the TAMs and MES-like cells and shows the highest score (Figures 5C–E). OPN encoded by the SPP1 gene has been shown to play critical roles in tumor progression by interacting with different receptors, such as integrin and CD44 receptors (Rangaswami et al., 2006). Here we found that the SPP1 gene mainly expresses in the macrophages rather than tumor cells, whereas the CD44 primarily expresses in the MES-like cells (Figure 5F), suggesting that TAM-secreted OPN may regulate the mesenchymal phenotype of glioma by interacting with CD44 on tumor cells. Furthermore, we confirmed the strong expression of CD44, OPN and VIMENTIN in this glioma by immunohistochemical analysis (Figure 5G). Collectively, these findings indicated that the macrophage-mediated SPP1/CD44 interaction might contribute to the induction of MES-like glioma cells.

Finally, we investigated the clinical significance of increased expression of SPP1 and CD44 in glioma. We found that high SPP1 expression correlate with a shorter overall survival time of glioma patients (Figure 6A). Similarly, high CD44 expression was observed to correlate with decreased overall survival time of glioma patients (Figure 6B). Moreover, we showed that there is a strongly positive correlation between SPP1 and CD68 expression levels in glioma samples (Figure 6C), Similarly, positive correlation between VIM and CD68 expression levels in glioma samples was observed (Figure 6D). Importantly, we also observed a significant correlation between CD44 and VIM expression in glioma samples (Supplementary Figure 3). These results suggested that high SPP1 and VIM expression also associates with an increased infiltration of macrophages in glioma. Consistently, we found that high CD68 and VIM expression correlate with poor prognosis of glioma patients (Figures 6E,F). Taken together, these findings suggested that MES-like cells enriched tumors show increased infiltration of macrophages and a worse prognosis.

In addition, we analyzed the infiltration of immune cells in human glioma tissues. Accordingly, we found that the infiltration of B cells is significantly lower in the tumors with high VIM expression than those with low VIM expression (Figure 6G). Moreover, we observed that the infiltration of CD4 naïve T cells and activated NK cells significantly decreases, whereas Treg cells and resting NK cells significantly increases (Figure 6G ). Importantly, we found that M2 macrophages are more abundant in the tumors with high VIM expression (Figure 6G). Together, these findings indicated that glioma tumors showing MES-like features are characterized by an immunosuppression microenvironment.

The glioma specimen was collected from a patient who received surgery in Shanghai General Hospital. The patient did not receive any treatment before the surgery and provided informed consent. This study was performed under the protocol (2019SQ175) approved by the Institutional Review Board at Shanghai General Hospital and in accordance with ethical guidelines (Declaration of Helsinki). The fresh tumor specimen was processed for scRNA-seq analysis in 4 h after resection. The histology of this tumor was confirmed by qualified pathologists.

Fresh tissue was transported in cold PBS on ice. Tissue was washed three times with PBS to remove blood, then cut into tiny pieces (∼1–2 mm³) and enzymatically dissociated by 1 mg/ml Type IV collagenase for 45 min at 37°C. The dissociation process was carried out in a shaking incubator with. Undigested tissue masses were removed using 70-μM strainers (Corning), and red blood cells (RBC) were removed by use of RBC lysis buffer (Biolegend). The dissociated cells were washed twice with PBS containing 0.04% bovine serum albumin (BSA, Sigma-Aldrich). Finally, cell viability was examined by Trypan blue (Invitrogen) staining.

In this study, 7,000 cells were loaded to generate the library. The library was constructed using 10X Genomics Chromium Single Cell 3′ kit (V2 chemistry) following the protocol provided by the manufacturer. The library was sequenced on the Illumina HiSeq platform using paired end 150 bp mode.

The raw sequencing reads were mapped to the human reference genome (GRCh38) using the CellRanger (version 4.0.0) with default parameters. The Gene-Barcode matrices were generated by calculating unique molecular identifiers (UMIs) and filtering non-cell-related barcodes, which contain the barcoded cells and gene expression counts. The output was then imported into Seurat (version 4.0.3) R toolkit for quality control (QC) and downstream function analysis (Zhou and Jin, 2020). Unless otherwise specified, all functions were run with default parameters. In particular, low quality cells (< 200 genes per cell, < 3 cells per gene and > 10% mitochondrial genes) were excluded. We used the LogNormalize function to normalize the data, which normalized the feature expression level through the total expression of each cell. We applied the FindVariableFeatures function to count highly variable features and returned the selected number of features (2000 by default). Principal component analysis (PCA) was used to reduce dimension, followed by clustering analysis in PCA space using a graph-based clustering method. t-SNE was used to visualize the identified clusters. Finally, the Seurat functions, including DotPlot, VlnPlot, FeaturePlot and Heatmap, were used to visualize the gene expression, and the FindAllMarkers function was employed to identify markers for a specific cluster against all remaining clusters.

The glioma tissue was fixed in 10% neutral buffered formalin (NBF) and embedded in paraffin. 5 µM thick sections were prepared and subjected to hematoxylin and eosin (H&E) and immunohistochemical (IHC) staining using standard protocols. The following primary antibodies were used: anti-GFAP (Abcam, ab68428), anti-IDH1 (Abcam, ab256557), anti-Vimentin (Cell Signaling Technology, 5,741), anti-CD44 (Abcam, ab157107), and anti-OPN (Abcam, ab214050). Peroxidase AffiniPure Goat anti-rabbit lgG secondary antibody (Jackson ImmunoResearch, 111-035-003) was used. Images were acquired using Leica DM6 B microscope.

Pseudotime trajectory analysis was performed using Monocle (version 2.20.0). Briefly, a group of ordering genes with differential expression among dictinct tumor cell clusters were selected through differentialGeneTest Function and filtered by q value < 0.01 and p value < 0.01. Then cells were projected onto the space and arranged into a trajectory with branch points. Cells in the same segment of the trajectory were grouped into the same state.

The CellChat (version 1.1.0) R package was used to infer, analyze and visualize cell-cell communication among cellular components in the glioma ecosystem. The Seurat normalized expression matrix was used for ligand-receptor interaction analysis. A known list of ligand-receptor pairs was obtained from CellChatDB, a database of literature-supported ligands and receptors interactions in both mouse and human. We first identified overexpressed ligands or receptors among cell types. The gene expression data were then projected onto the protein-protein interaction (PPI) network. If the ligands or receptors overexpressed, the interactions between the overexpressed ligands and receptors were recognized. The communication probability was inferred by calculating the communication probability of all ligands and receptors interactions associated with each signal pathway. The circle plot, bubble diagram or heatmap was used to visualize cellular interaction of specific signal pathway or ligand-receptor.

The cBioportal database (http://www.cbioportal.org/public-portal/) was used to access gene expression and survival information of glioma patients (Merged low-grade gliomas and GBM, TCGA, mRNA Expression (RNA seq V2 RSEM)) (McLendon et al., 2008; Brennan et al., 2013; Ceccarelli et al., 2016). First, the median value of the SPP1, CD44, CD68 and VIM genes in all tumors were calculated. Based on the gene expression value of each gene in individual tumor, tumors were allocated into high expression group and low expression group. Kaplan-Meier estimator was used to estimate the overall survival time of patients in each group. Statistical significances between the high and low expression group were calculated by Log-Rank (Mantel-Cox) test using GraphPad Prism 5.

In order to determine the infiltration of immune cells in human glioma tissues, the bulk RNA-seq data of human gliomas from the TCGA database was used. The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm was employed to infer the proportion of infiltrated immune cells as described previously (Newman et al., 2015). The quality of tumor samples was first examined and only samples with p value <0.05 were selected for next analysis. The scores of infiltrated immune cells in each sample were visualized by boxplot.