Dataset description
Single-cell transcriptome profiles of GBM were collected from two publicly available datasets. The scRNA-seq dataset by Abdelfattah et al. consisted of 40 tumor fragments from 16 GBM patients, including 11 newly diagnosed (termed as ndGBM in our research) and 5 recurrent (termed as rGBM in our research) cases. Additionally, we incorporated an scRNA-seq dataset from Ravi et al., which involved sorting lymphoid and myeloid cell populations (CD45+/CD3+) from eight patients with de novo GBM. Clinical details of the GBM patients are provided in Table S1. The clinical information was sourced from an open-access database and did not require an ethical declaration. The ethical approval was obtained for the original studies, with details available in the respective publications. To ensure data quality, cells with low read depths or high mitochondrial gene expression were excluded, and potential doublets were removed using the SoupX and DoubletFinder R packages as detailed in ‘Method’ section. In total, our analysis included 252,452 high-quality cells, with 211,376 cells from the Abdelfattah et al. dataset and 41,076 cells from the Ravi et al. dataset.
Main cell type annotation
After removal of the batch effect between samples, 252,452 single cells were clustered into 33 clusters with a resolution value of 0.8 in the FindClusters function. Cell types were annotated based on classic markers identified in previous studies and SingleR package. These cell types included B cells, characterized by the expression of CD19 and CD79A; T_NK cells, comprising both T and NK cells, identified by markers such as CD3D/G and NCAM1; myeloid cells, expressing LYZ, ITGAM, and CD14; and stromal cells including pericytes (expressing ACTA2 and PDGFRB), endothelial cells (expressing PECAM), and oligodendrocytes (expressing OLIG2 and MBP). Additionally, malignant glioma cells were characterized by the expression of SOX2, OLIG1, GFAP, and S100B (Fig. 1A, B, C). Cells from two datasets were integrated well (Fig. 1D). We quantified the relative tissue enrichment of major cell types in each sample and dataset (Fig. 1E). Consistent with prior studies, the proportions of cancer, stromal, and immune cells varied significantly across samples, potentially due to inherent differences in tumor phenotypes or the specific locations within the tumor from which biopsies were taken. Specifically, in the specimens from Ravi et al., nearly all identified cells were of immune origin, with only a minimal presence of stromal cells, except in the sample 'NC_P234_JAK' (Fig. 1E, Fig. S1A). All samples demonstrated a pronounced enrichment of myeloid cells, whereas the fractions of T_NK and B cells were relatively low. In contrast, the samples from the study by Ravi et al. predominantly contained glioma and myeloid cells, with other immune cells comprising approximately 10.6% of all cells analyzed. Notably, three fragments from ndGBM_01 (ndGBM_01_A, ndGBM_01_C, and ndGBM_01_D) and one from ndGBM_06 exhibited more than 80% glioma cells. Additionally, samples from recurrent GBM cases (rGBM_03, rGBM_04, and rGBM_05) contained higher proportions of glioma cells compared to their newly diagnosed counterparts (Fig. 1E). These results align with previous findings from flow cytometry and single-cell analyses (Abdelfattah et al., 2022; Fu et al., 2020; Müller et al., 2017), validating the accuracy of our cell type identification.
Inter- and intra-tumoral molecular heterogeneity of glioma cells
To explore the molecular heterogeneity of malignant glioma cells, we conducted de novo clustering of these cells into 34 distinct clusters (Fig. S1B). Consistent with expectations, glioma cells exhibited patient-to-patient heterogeneity, similar to many other cancer types (Fig. 1F). We carried out a copy number variation (CNV) analysis using inferCNV on 15 GBM fragments, each containing 2000 glioma cells. T_NK and Myeloid cells (randomly selected 2000 cells each cell type) were used for reference. This analysis revealed both shared and fragment-specific mutations across different regions within the same patient and among different patients (Fig. 1G), aligning with previously documented inter- and intra-tumoral genomic heterogeneity in GBM (Abdelfattah et al., 2022; Ravi et al., 2022). Notably, all fragments showed common alterations such as loss of chromosomes 10/10q, 13/13q, and gain of chromosome 7/7q, which are recurrent copy number changes in human GBMs (Verhaak et al., 2010). While major copy number changes were typically consistent across different fragments within each tumor, unique indels and mutational patterns also distinguished individual fragments in each patient (Fig. 1H). These findings corroborate the whole exome sequencing results reported by Abdelfattah et al (Abdelfattah et al., 2022).
Functional heterogeneityof TAM in GBM
Myeloid cells, including microglia and BMDM, constitute the predominant stromal compartment in the TME of GBM. To elucidate the cellular and molecular diversity within myeloid populations, we analyzed 114,351 myeloid cells through de novo clustering. This approach enabled the identification of three distinct myeloid cell subtypes characterized by unique gene expression profiles (Fig. 2A). TAM was noted for expressing well-known microglia markers such as P2RY12, TMEM119, BHLHE41, SORL1, SPRY1, and SRGAP28, alongside classical macrophage markers including C1AC, C1QA, and APOE. Within the BMDM subset, the cDC cell type marked by expression of traditional dendritic cell markers including CD1C, CLEC9A, CLEC10A, and LAMP3. Additionally, monocyte-derived cells were distinguished by high levels of S100A8, S100A9, and VCAN (Fig. 2B). All patients and fragments contributed to each myeloid cell types.
TAMs were further subjected to a second round of de novo clustering. In contrast to the traditional in vitro classifications of M1- or M2-like macrophages, TAMs in GBM displayed diverse behaviors in vivo (Fig. S2A), consistent with observations in other cancer types (Eisenbarth and Wang, 2023; Müller et al., 2017). To identify distinct TAM subtypes, we annotated these subtypes based on lineage markers (Fig. 2C, D). The TAM_MRC1 subtype was defined by high levels of MRC1, an established M2 polarization marker, as well as other M2-associated markers like CD163 and CD204/MSR1. Previous research has highlighted the role of MRC1, also known as ITGAM, in driving M2 polarization of macrophages potentially through the FAK/Akt1/β-catenin signaling pathway (Weng et al., 2019). On the other hand, the TAM_ISG15 subset exhibited increased levels of ISG15, along with markers linked to M1 polarization such as CD86, CXCL9, CXCL10, IRF1, and CD40. Interestingly, almost all identified clusters and subtypes co-expressed SPP1, C1QA, C1QB, and C1QC genes (Fig. 2D). This discovery contrasts with findings in colorectal cancer and breast cancer, where SPP1 expression typically shows mutual exclusivity with C1QC expression (Chen et al., 2021).
The functional phenotypes of TAMs are generally classified within an in vitro M1/M2 dual polarization framework. Subsequently, the AddModuleScore function from the Seurat package was employed to quantify the signatures of M1, M2, and pro-inflammatory gene sets (Supplementary Table S2). Notably, TAM_ISG15 macrophages displayed relatively higher M1 signatures, whereas TAM_MRC2 macrophages exhibited relatively higher levels of canonical M2 signatures, as expected. Additionally, TAM_NLRP3 macrophages demonstrated enhanced pro-inflammatory signatures (Fig. 2E). Collectively, these findings challenge the adequacy of the conventional in vitro polarization model and underscore the complexity of TAM phenotypes within the TME.
We subsequently analyzed the regulatory networks underlying each TAM subset using the SCENIC. We identified specific transcription factor (TF) regulons for each TAM subset (Fig. 2F, G). For instance, the TAM_ISG15 subset exhibited upregulation of positive regulators including ETV7, HESX1, IRF1, IRF2, IRF7, and IRF9, suggesting its potential role in enhancing immune responses. Similarly, the TAM_NLRP3 subset showed increased activity of ATF3, ARID5B, ECR1, and ECR3. Notably, in the TAM_MCR1 subset, the KLR8 and NRL regulons were significantly activated. These findings suggest potential candidates that could be driving the functional differences observed among TAM subsets, contributing to a deeper understanding of their heterogeneity and functions.
To investigate the association of TAM_MRC1 and TAM_ISG15 subtypes with specific signaling pathways, we performed GSEA analysis using the 'Hallmarks' gene set from the MsigDB. Our analysis indicated that the TAM_MRC1 subtype was significantly enriched in pathways such as HALLMARK_ANGIOGENESIS, HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_COAGULATION, and HALLMARK_GLYCOLYSIS. On the other hand, the TAM_ISG15 subtype exhibited significant enrichment in HALLMARK_INTERFERON_ALPHA_RESPONSE, HALLMARK_INTERFERON_GAMMA_RESPONSE, HALLMARK_KRAS_SIGNALING_DN, and HALLMARK_INFLAMMATORY_RESPONSE (Fig. 2H). It is well known that M2-like macrophages are recognized for secreting various factors that stimulate angiogenesis, including Vascular Endothelial Growth Factor (VEGF), Platelet-Derived Growth Factor (PDGF), and Transforming Growth Factor-Beta (TGF-β). These factors enhance endothelial cell proliferation and migration, crucial processes in angiogenesis (Xia et al., 2023). Epithelial-mesenchymal transition (EMT) is a key biological mechanism where epithelial cells transition into a mesenchymal state, involving loss of cell-cell adhesion and acquisition of migratory and invasive characteristics. This transformation is vital in cancer progression as it aids invasion and metastasis. M2 macrophages play a role in promoting EMT by releasing factors like TGF-β and IL-10, which induce EMT in nearby epithelial cells. This reciprocal interaction establishes a feedback loop where EMT supports M2 macrophage polarization, and conversely, M2 macrophages enhance EMT, contributing to a more aggressive tumor phenotype (Kim et al., 2021; Kim et al., 2023).
Detailed reclustering of T and NK cells
T and NK cells play critical roles in regulating tumorigenesis and cancer progression, but their characterization in GBM remains elusive. In our study, we classified the T and NK cells into five distinct subtypes based on specific cellular signature markers (Fig. 3A). These subtypes include T cells (n = 26,193), characterized by the expression of CD3D and CD3G; two NK subtypes, NK_CD56dim (n = 1,345) and NK_CD56bright (n = 993), identified by the relative expression levels of NCAM1 (CD56) and FCGR3A (CD16); NKT cells (n = 613), which exhibited markers of both T and NK cells; and ILC cells (n = 97), distinguished by the expression of KIT (Fig. 3B). Consistent with previous studies, we found significant variability in the proportions of these cell types across samples, with T cells being a substantially larger fraction than the other four cell types (Fig. 3C).
Detailed reclustering of CD4 + T and CD8 + T cells
To analyze T cells, we performed two rounds of unsupervised clustering. Initially, we distinguished between CD4 + T (CD4T) and CD8 + T (CD8T) cells based on the expression of CD4 and CD8A (Fig. 3D). Subsequently, the CD4 + T cells were further categorized into three subsets. The CD4T_FOXP3 cells, characterized by elevated expression of FOXP3, CTLA4, and TIGHT, were identified as regulatory T (Treg) cells. The CD4T_CCR7 subset, with high levels of CCR7 and SELL, was classified as Naive T cells. The CD4T_IL21 subset, known for significant IL21 expression, was recognized as helper T cells (Fig. 3D, E). The CD8 + T cells were also divided into three groups. The CD8T_LAG3 subset, distinguished by high expression of LAG3, PDCD1, and TIGHT, was identified as exhausted T (Trex) cells. The CD8T_GZMB subset, with robust expression of GZMB, GZMA, and PRF1, was classified as effector T (Teff) cells. The CD8T_GZMK subset, showing high levels of GZMK, was identified as effector memory (Tem) T cells (Fig. 3D, E). We utilized the AddModuleScore function to calculate regulatory, exhaustion, and cytotoxic scores for these T subsets, respectively. As expected, the Treg (CD4T_FOXP3) cells displayed significantly higher regulatory scores (Fig. 3F, G). The tumor-infiltrating CD8T_LAG3 cells exhibited a notably higher exhaustion score, whereas the effector CD8T_GZMB subset showed relatively higher cytotoxicity scores. Moreover, the CD8T_LAG3 (Tex) cells also demonstrated the second highest levels of cytotoxic effector signatures, indicating retained effector capabilities (Fig. 3H). Overall, the T cells displayed a dysfunctional state, emphasizing the complex immune responses within the GBM TME.
Detailed analysis of NK cells uncovering CD56dim_DNAJB1 dysfunctional state
Next, we aimed to delineate the specific characteristics of NK cells within the TME of GBM. Two distinct NK cell types were identified based on the expression levels of NCAM1 (CD56) and FCGR3A (CD16): CD56brightCD16lo and CD56dimCD16hi (Fig. 3A). CD56dimCD16hi NK cells primarily kill target cells through perforin and granzyme release, while CD56brightCD16lo NK cells have immunoregulatory functions and produce cytokines. In our study, CD56dimCD16hi NK cells showed increased expression of cytotoxic effector genes, while CD56brightCD16lo NK cells exhibited upregulation of specific cytokines, as expected.
Further analysis subdivided the CD56brightCD16lo NK cells into two subsets: NK_CD56bri_KLRK1 and NK_CD56bri_CXCR6 (Fig. 4A, B). The NK_CD56bri_KLRK1 subtype, distinguished by elevated levels of KLRK1, also known as NKG2D. KLRK1 is essential for enabling NK cells to activate, recognize, and eliminate tumor cells by interacting with stress-induced ligands on transformed cells (Wilton et al., 2019). This interaction is vital for triggering NK cells to attack and destroy tumor cells, a crucial mechanism for controlling tumor growth and promoting tumor regression (Terrén and Borrego, 2022; Wang et al., 2021). The NK_CD56bri_CXCR6 subtype is notable for its expression of cytokine genes such as CXCR6, CCL4, CCL4L2, CCL3, and CCL3L1. High expression of CXCR6 likely contributes to the increased infiltration of NK cells into the TME, facilitating their intense activation and potentially enhancing antitumor activity (Wang et al., 2021).
Similarly, the CD56dimCD16hi NK cells were classified into NK_CD56dim_CCL3 and NK_CD56dim_DNAJB1 subtypes (Fig. 4C, D). The NK_CD56dim_CCL3 subtype, characterized by high levels of CCL3, likely enhances antitumor immunity by promoting the recruitment of NK cells and other immune cells, such as DC and CD8 + T cells, thereby facilitating increased IFNγ production. This mechanism is crucial for activating and recruiting additional immune cells to the tumor site, potentially enhancing the efficacy of therapies such as anti-PD-1. The CD56dim_DNAJB1 subset was particularly noted for expressing stress response genes such as DNAJB1, HSPA1A, and HSP90AA1. This subset also showed elevated levels of NR4A1, identified as a key mediator of T cell dysfunction, aligning with characteristics of the recently described TaNK population (Tang et al., 2023).
We further analyzed the gene expression signature to investigate the functional differences between CD56dim_CCL3 and CD56dim_DNAJB1 subtypes. Notably, the CD56dim_DNAJB1 subset exhibited relatively higher cytotoxicity and stress scores (Fig. 4E). This suggests a unique role for the CD56dim_DNAJB1 subset in stress response and inflammation, potentially influencing the overall functionality of the NK cell population. To investigate the specific enriched signaling pathways, we performed GSEA analysis using the KEGG and 'Hallmarks' gene set from the MsigDB. Our analysis indicated that compared to CD56dim_CCL3, CD56dim_DNAJB1 showed significant upregulation in PD-L1 and PD-1 checkpoint pathway and HALLMARK_MTORC1_SIGNALING pathway (Fig. 4F). Recent studies have shown that activation of the PD-1 receptor hinders the cytotoxic functions of NK cells, leading to immunosenescence. The PD-1/PD-L1 checkpoint pathway interacts with different immunosuppressive cells, including Tregs, myeloid-derived suppressor cells, and M2 macrophages. Furthermore, the mTOR signaling pathway plays a crucial role in enhancing PD-L1 protein expression and facilitating immunosuppression (Salminen, 2024).
Complex intercellular communication networks in the GBM TME
To investigate cell-cell interactions and specific ligand-receptor pairs within the GBM TME, we performed CellphoneDB analyses, focusing on glioma alongside five TAM and four NK cell subtypes. Our analysis revealed significant communication between tumor and immune cells, particularly notable interactions between glioma cells and three TAM subsets: TAM_NLRP3, TAM_MRC1, and TAM_ISG15. Our results demonstrated that glioma cells exhibit elevated levels of amyloid precursor protein (APP). Additionally, high expression of CD74 was observed across five TAM subtypes, establishing a specific APP-CD74 ligand-receptor interaction. Research in various cancer models has highlighted that CD74 + TAMs can undermine immune surveillance, diminish the effectiveness of immune responses, and foster an immunosuppressive TME. Furthermore, the analysis indicated high expression of Amphiregulin (AREG), a member of the epidermal growth factor family, in four NK cell subsets. Epidermal Growth Factor Receptor (EGFR), which serves as the receptor for AREG as well as for EREG and HBEGF, was prominently expressed in glioma cells. Previous researches have shown that interactions mediated by EGFR and AREG may augment the invasiveness and metastatic potential of tumors (Paez et al., 2004; Singh et al., 2016). We also evaluated several inhibitory receptors on the surface of NK cells, including NKG2A, KLRB1 (CD161), CLEC2D, TIGIT, and PDCD1 (PD-1), suggesting the impaired cytotoxic activity of NK cells in the TME. This comprehensive analysis underscores the complex interplay of immune evasion mechanisms at play within the GBM TME, revealing potential targets for therapeutic intervention.