Single-cell transcriptomics reveals immunosuppressive microenvironment and highlights tumor-promoting macrophage cells in Glioblastoma

In our study, the scRNA-seq datasets were collected from two published studies with open-assess, so the ethical declaration is not applicable. However, ethical approvals for the original studies were obtained, and details can be found in the respective publications. For the datasets from Abdelfattah et al. [2], human tumor tissues were obtained under Institutional Review Board (IRB)-approved protocols (Pro00014547) at Houston Methodist Hospital, Houston, Texas and MD Anderson Cancer Center (PA 19-0661) following national guidelines. All patients signed informed consent during clinical visits before surgery and sample collection. For the datasets from Ravi et al. [25], the local ethics committee of the University of Freiburg approved the data evaluation, imaging procedures, and experimental design (protocol 100020/09 and 472/15_160880). These studies were approved by an institutional review board.

Single-cell transcriptome profiles of GBM were collected from two publicly available datasets. The scRNA-seq dataset from the study by Abdelfattah et al. [2] contains 44 tumor fragments from 18 glioma patients, including two low-grade gliomas (grades II), 11 ndGBM (grade IV) and 5 rGBM (grade IV). To minimize the patients’ variance, we downloaded the raw fastq sequences from 40 high-grade GBMs datasets from GSE182109↗. We also collected a scRNA-seq dataset from Ravi et al. [25], where lymphoid and myeloid populations (CD45 + /CD3+) were sorted from eight patients diagnosed with de novo glioblastoma. We downloaded the matrix Seurat object from https://osf.io/4q32e/https://doi.org/10.17605/OSF.IO/4Q32E↗.

The raw Illumina sequencing reads from 40 samples sequenced by Abdelfattah et al. were aligned to GRCh38 (human) using Cell Ranger (V7.2.0) software with default parameters. Subsequently, genes were quantified as UMI counts using Cell Ranger. Downstream analysis was performed on filtered feature counts generated by Cell Ranger, and low-quality single cells containing < 200 expressed genes or > 20% mitochondrial transcripts were removed. Additionally, genes expressed in fewer than three single cells were removed.

SoupX (V1.6.0) [26] was used to estimate and remove cell-free mRNA contamination in droplet-based sc-RNAseq data. The raw_feature_bc_matrix and filtered_feature_bc_matrix resulted results from Cell Ranger were used. Additionally, we identified potential single-cell doublets using DoubletFinder (V2.0.3) (https://github.com/chris-mcginnis-ucsf/DoubletFinder↗), with an expectation of a 7.5% doublet rate assuming Poisson Statistics.

Then scRNA-seq data from different samples were merged into Seurat (V4.2.0) [27] R package, then normalized using the NormalizeData function and scaled on 2000 most variable features using the FindVariableFeatures and ScaleData functions. Computed principal components were batch corrected for sample variations using the Harmony R package (V1.0). We used batch-corrected PCs as input for Louvainbased graphing and chose resolution parameters with 0.5, 0.8 and 1.0. FindAllMarkers function was used to identify cluster-specific marker genes and visualization with dot and feature plots. The genes expressed in each cluster and SingleR (V1.4.1) [28] were examined to identify the cell types.

We used clusterProfiler [29] package to identify and visualize enriched pathways in our subsets. Differentially expressed genes (DE genes) between cell types were identified with FindMarkers or FindAllMarkers function with default parameters. Subsequently, we employed the enrichGO and enrichKEGG functions to identify enriched pathways for Gene Ontology (GO) enrichment analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, respectively.

Copy-number variations (CNVs) were estimated by the inferCNV software package (https://github.com/broadinstitute/inferCNV↗). The software aligns genes to their chromosomal location and applies a moving average to the relative expression values, with a sliding window of 100 genes within each chromosome. Malignant glioma cells from some patients were used as observation datasets, and T_NK and myeloid cells (randomly selected 2000 cells for each cell type) were used as normal reference set.

To evaluate the functional differences among NK cell subsets, we created an extensive compilation of cytotoxicity, inflammatory, and stress-related gene sets and activating and inhibitory receptors from previous research. The gene sets mentioned in Zhang et al.‘s [30] article (refer to S2 Table for more information) were included. We computed the scores for each gene set using the AddModuleScore function in Seurat.

To investigate potential interactions among various cell types within the GBM TME, a cell-cell communication analysis was conducted using CellPhoneDB (V3.0) [31]. This analysis included cells annotated as glioma, five TAM subtypes and four NK subtypes. CellPhoneDB facilitated the assessment of cell-cell interactions by analyzing enriched receptor-ligand interactions, which depend on the expression of a specific receptor by one cell type and its corresponding ligand by another. Subsequently, we selected interactions that were biologically significant.

To identify key transcription factors (TFs) across various cell types, a cis-regulatory analysis was conducted using pySCENIC (V0.11.2) [32,33]. The analysis focused on five subtypes of TAM. TFs were identified utilizing GENIE3, grouped into modules known as regulons, and analyzed via RcisTarget employing default gene-motif rankings. Subsequently, regulon activity for each cell in the dataset was assessed using AUCell.

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 study) cases. Additionally, we incorporated a 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 S1 Table. The clinical information was sourced from an open-access database and did not require an ethical declaration. 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.

After the 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 the 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. 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–1C). 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, S1A Fig). 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 [2,11,12], validating the accuracy of our cell type identification.

To explore the molecular heterogeneity of malignant glioma cells, we conducted de novo clustering of these cells into 34 distinct clusters (S1B Fig). Consistent with expectations, glioma cells exhibited patient-to-patient heterogeneity, similar to many other cancer types (Fig 1F). We carried out a CNV analysis using inferCNV on 15 GBM fragments, each containing 2000 glioma cells. T_NK and Myeloid cells (randomly selected 2000 cells for each cell type) were used for reference. This analysis revealed 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 [2,25]. 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 [34]. 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 exome sequencing results of Abdelfattah et al. [2].

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 is marked by the 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 type.

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 (S2A–S2C Fig), consistent with observations in other cancer types [12,35]. To identify distinct TAM subtypes, we annotated these subtypes based on lineage markers (Fig 2C, 2D). The TAM_MRC1 subtype was defined by high levels of MRC1, an established M2 polarization marker, and 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 [36]. 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 [8].

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 (S2 Table). 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). 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 TF regulons for each TAM subset (Fig 2F, 2G). 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 drive 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 [37]. 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 [38,39].

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_CD56^dim (n = 1,345) and NK_CD56^bright (n = 993), identified by the relative expression levels of NCAM1 (CD56) and FCGR3A (CD16). NKT cells (n = 613) exhibited markers of both T and NK cells. 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).

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 FOXP3, CTLA4, and TIGHT expression, 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, 3E). 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, 3E). 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, 3G). 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.

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): CD56^brightCD16^lo and CD56^dimCD16^hi (Fig 3A). CD56^dimCD16^hi NK cells primarily kill target cells through perforin and granzyme release, while CD56^brightCD16^lo NK cells have immunoregulatory functions and produce cytokines. In our study, CD56^dimCD16^hi NK cells showed increased expression of cytotoxic effector genes, while CD56^brightCD16^lo NK cells exhibited upregulation of specific cytokines, as expected.

Further analysis subdivided the CD56^brightCD16^lo NK cells into two subsets: NK_CD56^bri_KLRK1 and NK_CD56bri_CXCR6 (Fig 4A, 4B). The NK_CD56^bri_KLRK1 subtype is 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 [40]. 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 [41,42]. The NK_CD56^bri_CXCR6 subtype is notable for expressing 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 [42].

Similarly, the CD56^dimCD16^hi NK cells were classified into NK_CD56^dim_CCL3 and NK_CD56^dim_DNAJB1 subtypes (Fig 4C, 4D). The NK_CD56^dim_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 Interferon-gamma(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 CD56^dim_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 [43].

We further analyzed the gene expression signature to investigate the functional differences between CD56^dim_CCL3 and CD56d^im_DNAJB1 subtypes. Notably, the CD56dim_DNAJB1 subset exhibited relatively higher cytotoxicity and stress scores (Fig 4E). This suggests a unique role for the CD56^dim_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 CD56^dim_CCL3, CD56^dim_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 is crucial in enhancing PD-L1 protein expression and facilitating immunosuppression [44].

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 (Fig 5A, 5B). Our results demonstrated that glioma cells exhibit elevated levels of amyloid precursor protein (APP) (Fig 5C). 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 (Fig 5C). 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 (Fig 5D). Previous research have shown that interactions mediated by EGFR and AREG may augment the invasiveness and metastatic potential of tumors [45,46]. We also evaluated several inhibitory receptors on the surface of NK cells, including NKG2A, KLRB1 (CD161), CLEC2D, TIGIT, and PDCD1, suggesting the impaired cytotoxic activity of NK cells in the TME. This comprehensive analysis underscores the complex interplay of immune evasion mechanisms within the GBM TME, revealing potential targets for therapeutic intervention.