Introduction

Glioblastoma (GBM) is a fatal disease without effective long-term treatment options. The current standard of care consists of tumor resection followed by adjuvant chemoradiotherapy resulting in a median overall survival of only 14 months [1]. One of the hallmarks in GBM progression is the high rate of neovascularization. The GBM-induced aberrant vessels not only nourish glioma cells, but also provide a specialized niche for tumor-associated stromal and immune cells such as monocyte-derived macrophages (MdMs), yolk sac-derived microglia (MG) (together termed glioma-associated macrophages/microglia, GAMs), and peripheral adaptive immune cells. This immune tumor microenvironment (iTME) paradoxically acts in an immunosuppressive manner and promotes tumor progression [2]. For example, clinical trials of systemic T cell checkpoint blockade showed only disappointing results [3, 4], which was attributed in part to the immunosuppressive components of the GBM iTME. The origin of GAMs, infiltration of peripherally derived macrophages across the blood-brain-barrier (BBB) or recruitment of tissue-resident MG to the tumor site, as well as their contribution to gliomagenesis are studied intensively [2, 5–7]. Hence, major efforts have been undertaken to describe the GBM iTME on a single cell level [5, 6, 8], or dissect the composition and changes upon disease stages, recurrence and immunotherapy specifically within the GAM compartment [9–11]. However, human data on the composition of the iTME in different tumor regions (contrast enhancing tumor center versus peripheral infiltration zone) remain scarce [12, 13].

To study the region-dependent cellular diversity within individual GBMs, we performed single-cell RNA sequencing (scRNA-seq) on patient-matched biopsies from the tumor center and the peripheral infiltration zone of five primary GBM patients. Additionally, peripheral blood mononuclear cells (PBMC) of the same patients were included to explore the transcriptional changes occurring during tumor infiltration of circulating immune cells.

Our analysis revealed a regionally distinct transcription profile of MG and MdMs and an impaired activation signature in the tumor-peripheral cytotoxic-cell compartment. Comparing tumor-infiltrating CD8⁺ T cells with PBMC-derived, identified CX3CR1^high and CX3CR1^int CD8⁺ T cells with effector and memory phenotype, respectively, enriched in blood but absent in the iTME. Tumor CD8⁺ T cells displayed features of tissue-resident memory T cells and were characterized by an exhaustion phenotype. This work provides a regionally-resolved map of transcriptional states in glioma-associated cell types complemented by patient-matched PBMCs, dissecting the composition and molecular diversity of the iTME in GBM.

Results

scRNA-seq analysis of paired tumor center, periphery and PBMC samples

Fresh, neurosurgically resected tissue from five primary, treatment naïve GBM patients were harvested (Figure 1A, Supplementary File 1). According to the 2021 WHO Classification of Tumors of the Central Nervous System [14], in which the term glioblastoma designates only IDH-wildtype grade 4 tumors, we will hence use the term grade 4 glioma, as we included as well IDH-mutant grade 4 tumors (Supplementary File 1). The tumor center was defined as contrast enhancing, whereas the tumor periphery was defined as T2 hyperintense by magnetic resonance imaging (MRI)-guided, navigated surgical resection (Figure 1B). Increased cellular density of the center vs. periphery samples was confirmed by nuclear DAPI staining on matched histological micrographs of the resected tissue specimens used for scRNA-seq (Figure 1C). As outlined in Figure 1A, we separately processed patient tumor and blood samples and enriched them for immune cells by fluorescence-activated cell sorting (FACS) (Figure 1 – figure supplement 1A and 1B). The three samples per patient (center, periphery and PBMC) were loaded on different wells of a 10x Genomics Chromium system for a targeted recovery of 10,000 cells. Due to technical issues cells from the center sample of patient BTB 609 could not be collected.

In total we analyzed 45,466 cells that passed initial quality control and filtering, comprising 8,254 cells from tumor center, 5,954 cells from tumor periphery and 31,258 PBMCs, with 6,354 to 10,957 cells per patient (Supplementary File 2; Figure 1 – figure supplement 1C-1E, Supplementary File 3). All cells were projected onto a two dimensions t-distributed stochastic neighbor embedding (tSNE) [15]. As we observed a good overlap of cells across patients for most of the dataset (see Methods) (Figure 1 – figure supplement 2B – 2F), we chose not to perform any correction for patient-specific effects. Using hierarchical clustering, the cells were partitioned into clusters (Figure 2 – figure supplement 1A, 1B) which were then annotated into nine distinct cell types for the immune subset, including two transcriptionally distinct MG subsets (MG_1 and MG_2) and four cell types for the CD45 negative subset (Figure 2A; Figure 2 – figure supplement 1C, 1D; Supplementary File 3).

Figure 1

Figure 2

In more detail, our annotation strategy made use of the relative expression patterns of known marker genes (Figure 2B), and of cluster-specific genes (Figure 2 – figure supplement 2). Additionally, unbiased cell-type annotation using whole-transcriptome comparisons to reference bulk and single-cell datasets was performed with the SingleR package. We used as reference (i) a public bulk RNA-seq dataset of sorted immune cell types from human PBMC samples [16], which helped in identifying the major immune cell lineages: B cells, CD4⁺ and CD8⁺ T cells, dendritic cells (DCs), monocytes and NK cells (Figure 2 – figure supplement 1E). (ii) A bulk RNA-seq dataset of sorted immune cell types from the TME of human gliomas [5] was helpful to separate GBM-associated CD4⁺ and CD8⁺ T cells, as well as microglia and MdMs (Figure 2 – figure supplement 1F). (iii) A 10X genomics scRNA-seq dataset of the innate immune TME of seven newly diagnosed GBM patients [9] confirmed the annotation of the innate immune subset (Figure 2 – figure supplement 1I), which was also supported by a microglia and a macrophage signature scores defined using gene lists obtained from [7] (Figure 2 – figure supplement 1J, 1K). (v) Finally a Smartseq2 scRNA-seq dataset of IDH1^wt tumors [17] was used to characterize the CD45^neg population (Figure 2 – figure supplement 1N).

To perform a differential expression analysis between tumor sites, we stratified the analysis by cell type to decrease the influence of differential abundance patterns (analyzed separately, see below). A principal component analysis (PCA) on the pseudo-bulk aggregated samples (see Methods) confirmed that the major source of variation was due to differences between lymphoid, MG/myeloid cells, and CD45^- cells (PCs 1 and 2; Figure 2C). There was no clear association to the IDH1 status of patients on these or deeper components. However, our study was neither designed nor powered to find regional signatures within the iTME depending on IDH1 status, but rather to identify common transcriptional differences within the iTME between tumor center, periphery and PBMC of IDH1^wt and IDH1^mut grade 4 glioma.

Differential abundance analysis between tumor center and periphery did not reveal significant changes, but some populations displayed suggestive differences: the abundance of both CD4⁺ T cell and monocytic clusters increased in periphery compared to center, while the abundance of MdM clusters decreased (Figure 2D). This has been observed by others as well [12, 18]. Testing the abundance at a finer scale with miloR provided us interesting insights. Indeed, within the innate immune population (Figure 2E) we confirmed the decreased abundance of MdMs in the tumor periphery and noticed an increased abundance of the MG_1 subset of MG in the tumor periphery. (Figure 2F, 2G).

MG display regionally resolved transcriptional profiles that differ from those of DCs and MdMs

Differential expression analysis between MG from tumor center and periphery revealed downregulation of inflammatory-related genes in the peripheral MG (FCGBP and CCL20), downregulation of genes associated with canonical interferon (IFN) responses (IFI6, IFI27, STAT1, ISG15), and cell proliferation (STMN1) as well as downregulation of scavenger receptor gene CD163 (Figure 3A and Supplementary File 4). The latter was shown to have a potential role in the phagocytic response of MG to beta-amyloid depositions identified in single-nucleus RNA-seq of postmortem human brain[19].

Figure 3

Among many upregulated genes whose function is not yet known, we found Inhibitor of DNA-Binding 1, also known as Inhibitor of Differentiation 1 (ID1) to be increased in the peripheral MG (Figure 3A). ID1 is well described in GBM progression, treatment resistance and glioma stem cell biology [20]. However, new evidence has emerged, linking ID1 to suppression of the anti-tumor immune response in the myeloid compartment and promoting tumor progression [21].

To further explore the underlying biological processes differing between MG in the two compartments, we conducted a gene set enrichment analysis (GSEA) on the results of the differential expression analysis using Gene Ontology (GO) database (Biological Processes). This revealed overall a significant downregulation of GO categories involved in antigen processing and presentation via MHC-I and MHC-II in the peripheral MG relative to the center MG, as well as downregulation of amino acid metabolism and TNF-α signaling pathway (Figure 3B), which supported the conclusion of an impaired activation state of peripheral MG.

When comparing the transcriptional profiles of the other innate immune phagocytic populations, DCs and MdMs, we observed for both cell types upregulation of pro-inflammatory related gene sets in the tumor periphery compared to tumor center (Figure 3C, 3D, Figure 3 - figure supplement 1A, 1B and Supplementary File 4). Within DCs, the induced pro-inflammatory phenotype was also observed when comparing DCs from tumor periphery to PBMC-derived ones (Figure 3 -figure supplement 1C, 1D and Supplementary File 6). Taken together, DCs and MdMs seem to have a proinflammatory phenotype in the glioma periphery as opposed to MG. However, they are less abundant there, at least as far as MdMs are concerned (Figure 2E-2G).

The iTME of grade 4 glioma harbors two transcriptionally distinct MG subpopulations

Unsupervised hierarchical sub-clustering of the MG population revealed two transcriptionally distinct iTME MG subsets, which we termed MG_1 and MG_2, respectively (Figure 3E, Figure 3 – figure supplement 1E and Supplementary File 5). The MG_1 cluster was highly enriched for the activator protein-1 (AP-1) family of transcription factors including FOS, FOSB, JUN, JUNB, MAF and MAFB (Figure 3F and Supplementary File 5), which convey a surveilling phenotype to adult MG, but are also involved in numerous processes including cell growth, differentiation and immune activation [22]. Specifically, FOSB gene products were implicated in the excitotoxic MG activation by regulating complement C5a receptor expression [23]. Yet, concomitant upregulation of anti-inflammatory Krüppel-like factor 2 (KLF2) [24] and Dual Specificity Protein Phosphatase 1 (DUSP1), an inhibitor of innate inflammation by negatively regulating the mitogen-activated protein kinase (MAPK) pathway [25], together with increased expression of anti-proliferative genes like RHOB, BTG1 and BTG2 paint a more complex picture of these cells. Particularly, BTG1 was identified as an activation-induced apoptotic sensitizer in MG after exposure to inflammatory stimuli [26], serving as an autoregulatory mechanism and possibly hinting towards an exhausted state in these MG_1 cells. GSEA for differences between MG_1 and MG_2 clusters using the MSigDB Hallmark collection of major biological pathways [27] revealed downregulation of many MG effector functions in the MG_1 population including (1) inflammation (“Complement”, “Allograft Rejection”, “Reactive Oxygen Species Pathway”), (2) immune cell activation (“IFN-α Response”, “IFN-ψ Response”, “IL6 JAK STAT3 Signaling”) and (3) immunometabolism (“Fatty Acid Metabolism”, Oxidative Phosphorylation”, “Glycolysis”) (Figure 3G). As we examined the expression of previously described reactivity markers of MG including C1QA, FCGR1A, CD14, HLA-DRA, TREM2 and Ferritin (FTH1) [28–31] (Figure 3H and Figure 3 – figure supplement 1F), and established MG homeostatic genes like CX3CR1, HEXB and SPI1 (PU.1) (Supplementary File 5), we noted a reduced expression of these genes in the MG_1 cluster compared to MG_2 cells, while the anti-inflammatory transcription factors NR4A1 [32] NR4A2 [33] were highly upregulated (Figure 3 – figure supplement 1G). Additionally, when comparing the distribution of MG subsets between sites among total MG, we observed an increased abundance of MG_1 cells in the tumor periphery, which might at least partially explain the previously observed impaired activation state of total peripheral MG (Figure 3 – figure supplement 1H). To rule out this explanation we stratified the differential expression analysis between tumor periphery and center by MG subtype and observed similar genes downregulated in the tumor periphery in each subtype (e.g. CCL20, FCGBP) (Supplementary File 5), arguing that the non-reactive phenotype is a common feature for all peripheral MG and is not only driven by a differential abundance pattern of one MG subpopulation relative to the other.

The tumor peripheral cytotoxic cell compartment exhibits an impaired activation signature

Next, we investigated the regional differences in the lymphoid compartment composed of CD4⁺ and CD8⁺ T cells and natural killer (NK) cells. We observed only very few significant changes in the transcriptomic profiles of CD4⁺ T cells between tumor center and periphery (Figure 4 – figure supplement 1A and Supplementary File 4). Yet, comparing peripheral CD8⁺ T cells with CD8⁺ T cells from tumor center revealed 110 differentially expressed genes (43 genes upregulated and 67 genes downregulated) (Figure 4A and Supplementary File 4). Many downregulated genes in the peripheral CD8⁺ T cells associated with canonical IFN responses (IFI6, IFI27, MX1, STAT1, EPSTI1 PARP9, ISG15) [34] cell proliferation (STMN1, CENPF, HELLS, NUSAP1 and DNPH1) and T cell co-stimulation (CD28, TMIGD2 (CD28H), TNFRSF4 (OX40), CD27 and TNFRSF18 (GITR)) (Figure 4A). Contrary to our expectations, we saw upregulation of CTLA4 in the center CD8⁺ T cells which acts as a negative costimulatory molecule. However, unlike other costimulatory receptors, such as CD27 and CD28, CTLA-4 is not constitutively expressed on T lymphocytes [35]. but only induced following T cell activation, along with positive costimulatory molecules such as OX40 and GITR. In addition, upregulation of CTLA-4 requires entry into the cell cycle [35]. In line with that, we detected an upregulation of proliferative genes in center CD8⁺ T cells. In summary, CTLA-4 induction in center CD8⁺ T cells rather suggested T cell activation than exhaustion, especially since other inhibitory receptors like PDCD1 (PD-1), LAG3 and HAVCR2 (TIM-3) were not differentially expressed between sites. Moreover, we did not observe differential expression of genes involved in CD8⁺ T cell effector functions like cytotoxicity (e.g., GZMK, GZMB, KLRG1, PRF1) or cytokines (e.g., CCL5, XCL1, XCL2, IL10). Yet, we noted upregulation of inhibitory genes (TGFB1 and FCRL6 [36]) in the peripheral CD8⁺ T cells, suggesting that a pool of activated, proliferating and IFN-responsive CD8⁺ T cells is present in the tumor center, however, absent in the tumor periphery.

Figure 4

Similar trends were observed for the peripheral NK cell population with peripherally reduced IFN response (MX1 and IFI44L), and proliferative genes (STMN1, HELLS, CENPF, PTTG1 and DNPH1), downregulated stimulatory receptors (TMIGD2 (CD28H) and TNFRSF18 (GITR)), and reduced NF-κB signaling (NFKB1 and RELB) (Figure 4B and Supplementary File 4). Although, we observed upregulation of key genes associated with NK cell effector function in the periphery (e.g., FCGR3A (CD16), FGFBP2, ITGB2, GZMH and KIR2DS4), increased expression of inhibitory receptors like LILRB1 and KLRG1, the latter especially in co-expression with chemokine receptor CX3CR1, identified the peripheral NK cells rather to be terminally differentiated with impaired cytotoxic capabilities [37]. This was in line with the observed abrogated cytokine activity profile in the peripheral NK cells with reduced expression of key factors like XCL1, XCL2, LTB and CKLF. In summary, our data revealed an impaired activation signature in the peripheral cytotoxic cell compartment.

CX3CR1 labels a specific CD8⁺ T cell population in the circulation of grade 4 glioma patients

Next, we investigated the relationships between circulating CD8⁺ T cells and those from the tumor milieu and, more specifically, the peripheral, infiltration zone characterized by an abrogated CD8⁺ T-cell IFN response and activation signature. Strikingly, there were large transcriptomic differences between PBMC and periphery CD8⁺ T cells (Figure 5A), with 1,417 differentially expressed genes (864 genes upregulated in the tumor periphery and 553 genes upregulated in PBMC) (Figure 5B, Supplementary File 6).

Figure 5

Interestingly, one of the key genes upregulated in PBMC CD8⁺ T cells was the chemokine receptor CX3CR1 (Figure 5B). Flow cytometry of an additional matched glioma grade 4 patient cohort confirmed an increased abundance of CX3CR1⁺ CD8⁺ T cells in PBMC compared to almost absent CX3CR1⁺ CD8⁺ T cells in tumor periphery (Figure 5C, Figure 5 – figure supplement 1A, Supplementary File 1). Unsupervised hierarchical sub-clustering of CD8⁺ T cells revealed that CX3CR1 expression labelled a specific subpopulation among PBMC CD8⁺ T cells (Figure 5D, 5E).

Recently, expression of CX3CR1 was demonstrated to distinguish memory CD8⁺ T cells with cytotoxic effector function in healthy donors and different inflammation related conditions [38–41]. Further characterization of classical central memory (T_cm) and effector memory (T_em) populations by varying surface expression levels of CX3CR1 identified a novel CX3CR1^int subpopulation, termed peripheral memory (T_pm). T_pm cells underwent frequent homeostatic divisions, re-acquired CD62L, homed to lymph nodes, and predominantly surveyed peripheral tissues compared to T_cm and T_em [39]. In our dataset, the circulating CX3CR1⁺ CD8⁺ T cells indeed displayed a core signature of memory CD8⁺ T cells with effector function, comprising expression of LFA-1 (IGAL-ITGB2), TBX2 (Tbet), SELL (CD62L), GZMB (Granzyme B) and PRF1 (Perforin 1) (Figure 5F), separating them from circulating CX3CR1^- CD28^high, IL7R^high and CD27^high naïve CD8⁺ T cells (Figure 5G and Figure 5 – figure supplement 1B). Flow cytometric analysis confirmed T_eff to be CX3CR1^high, with negligible expression levels in the naïve CD8⁺ T cells, whereas the identified memory CD8⁺ T cells (T_em and T_pm) were CX3CR1^int (Figure 5H, 5I). Collectively, surface expression analysis of CX3CR1 identified a subset of CX3CR1^high T_eff and CX3CR1^int memory (T_em, T_pm) CD8⁺ T cells in the circulation of grade 4 glioma patients with potentially elevated tissue surveilling properties in the case of T_pm, which are, however, largely absent in the tumor microenvironment.

CD8⁺ T cells in the tumor periphery share features with tissue-resident memory T cells (T_rm)

We next examined the differing transcriptional and surface-specific features between tumor infiltrating and circulating CD8⁺ T cells. Surface staining for CD45RA and CD45RO, discriminating naïve/effector from memory T cells, attributed a predominant CD45RO⁺ memory phenotype to the tumor infiltrating CD8⁺ T cells (Figure 6A, 6B). Interrogation of the transcriptomic profile of these cells revealed a key marker expression signature consistent with tissue-resident memory T cells (T_rm): Expression of cellular adhesion molecules (integrins) ITGA1 (CD49a) and ITGAE (CD103), tissue retention marker CD69, chemokine receptors implicated in tissue-homing CXCR3, CXCR6 and CCR5 [42] and transcription factors, ZNF683 (Hobit) and PRDM1 (Blimp1) as well as reduced expression of TBX21 (Tbet) and EOMES [43], strongly suggested a T_rm phenotype for these cells (Figure 6C, 6D and Figure 6 – figure supplement 1A). Co-expression analysis of paired PBMC and tumor periphery samples using flow cytometry showed that CD69⁺ CD103^- and CD69⁺ CD103⁺ cells are the dominant CD8⁺ T cell populations in the tumor periphery (Figure 6E, 6F). Combined, these data strongly suggest a T_rm phenotype for the CD8⁺ T cells in the tumor periphery.

Figure 6

Previous reports of T_rm populating the brain in the aftermath of central or peripheral infections concluded that brain T_rm cells surveil the brain tissue and mediate protection by rapid activation and enhanced cytokine production [42]. Indeed, CD8⁺ T cells in the tumor periphery showed increased expression of genes belonging to costimulatory pathways, including ICOS, TNFRSF4 (OX40) and TNFRSF9 (4-1BB) (Figure 6 – figure supplement 1B, Supplementary File 6), albeit accompanied by high levels of inhibitory receptors PDCD1 (PD-1), LAG3, HAVCR2 (TIM-3) and CTLA4 (Figure 6G). Moreover, expression of genes coding for cytotoxic molecules, including Granzyme B and Perforin 1 were decreased in the peripheral CD8⁺ T cells, suggesting a compromised killing capacity of these cells. And lastly, CD8⁺ T cells in the tumor periphery exhibited a transcription factor profile of exhausted T cells with high expression of NR4A1, MAF and IRF4 (Figure 6G and Figure 6 – figure supplement 1C), which were implicated in T cell dysfunction and exhaustion [44, 45]. Collectively, these data indicate that CD8⁺ T cells in the glioma periphery share features with T_rm cells. However, inhibitory receptor expression, functional molecules and transcriptional signature ascribe an exhausted phenotype to these cells.

Noteworthy, we observed high upregulation of similar genes in the comparison tumor periphery vs. PBMC for CD4⁺ T cells as for CD8⁺ T cells (Figure 6 – figure supplement 1D). These included transcription factor family NRA41-3, identified as key mediator of T cell dysfunction [45], Dual Specificity Protein Phosphatase 2/4 (DUSP2, DUSP4) described as negative regulators of mitogen-activated protein (MAP) kinase superfamily and associated with impaired T cell effector activity [46] and T cell senescence [47], and transcription factor CREM which was implicated in IL-2 suppression [48]. These genes could potentially identify pan T cell dysfunction markers within the GBM iTME [49].

Interrogation of cell-cell interactions revealed critical role of SPP1-mediated crosstalk between MG and lymphocytes in the tumor periphery

We finally investigated cell-cell interactions based on ligand-receptor expression levels using the CellChat platform [50]. Considering that MG and lymphocytes displayed an impaired activation signature in the tumor periphery, we focused our analysis on the tumor-peripheral crosstalk between these cells (Figure 7A). This revealed SPP1 (Osteopontin) as a leading potential cell-cell interaction mediator between MG and lymphocytes (Figure 7A and 7B). MG SPP1-mediated signaling was as well among the most significant interactions, when investigating cell-cell communication across all cell types and both sites (Figure 7 – figure supplement 1A, 1B). Further, we found that SPP1 is mainly expressed by MG rather than glioma cells, contrary to previous reports [51] (Figure 7C, Figure 7 – figure supplement 1C, 1D). MG SPP1 conveys different interactions, depending on the recipient cell binding receptor expression profile. The predicted interactions of MG SPP1 with NK cells could be mediated via the integrin complex ITGA4-ITGB1 (CD49d-CD29) (Figure 7A, 7C), whereas CD4⁺ and CD8⁺ T cells exhibited strong interactions with MG via SPP1/CD44 interaction (Figure 7A, 7B). The SPP1/CD44 axis was recently described to suppress T cell activation and proliferation [52]. The identified MG SPP1-mediated interactions might represent potential targets to modulate MG-lymphocyte crosstalk in the tumor periphery.

Figure 7

Discussion

In this study, we combined single-cell RNA sequencing and flow cytometry to interrogate the regional leukocyte activation signature in patient-matched biopsies from contrast-enhancing tumor center, infiltrative peripheral rim, and blood PBMCs of grade 4 glioma patients. Our analyses revealed a distinct, regionally dependent transcriptional profile for most of the investigated cell populations. While peripheral MG and cytotoxic cells predominantly displayed an impaired activation signature, MdMs showed pro-inflammatory traits in the tumor periphery, however, were less abundant there compared to the tumor center, which was reported by others as well [12, 13]. Supplemented with transcriptional analysis of paired PBMC samples, we provide an in-depth characterization of the three main immunological compartments of grade 4 glioma.

Previous studies focused on the description of the TME of grade 4 glioma, which also considered regional differences, yet they focused primarily at neoplastic cells rather than the immune compartment [12]. Others investigated the differences in the iTME composition between primary and metastatic brain tumors [5, 6]. Interestingly, the two latter ones reported differences in the iTME composition between IDH1^wt and IDH1^mut glioma. Of note, both authors included low-grade and even pre-treated recurrent glioma patients into the IDH1^mut group, representing a quite heterogenous patient cohort. In this study, we aimed at providing a representative selection of primary, treatment-naïve grade 4 glioma patients including IDH1^wt and IDH1^mut to identify common transcriptional differences within the iTME between tumor center, periphery and PBMC of grade 4 glioma. How far these regional differences vary between IDH1^wt and IDH1^mut grade 4 glioma merits further investigation.

We identified a transcriptionally distinct MG subcluster, MG_1, which was enriched in the infiltrative tumor periphery and displayed an anti-inflammatory/non-reactive phenotype. A similar MG subpopulation expressing a comparable gene signature was recently described to be enriched in Alzheimer’s disease patients [53].

Additionally, the peripheral cytotoxic cell compartment exhibited an impaired activation state, including a downregulated IFN response signature in CD8⁺ T cells. Induction of an IFN response state was described as a consequence of T cell receptor-mediated IFN-ψ production, likely serving as an autocrine response and inducing the proliferative program [34]. Hence, the reduced autocrine IFN-responsive state in the tumor peripheral CD8⁺ T cells, together with downregulated proliferative and co-stimulatory genes emphasized their impaired activation in the peripheral infiltration zone. Recently, we showed that the response to immunotherapy in GBM is indeed region dependent. For this, we cultured GBM explants in perfusion bioreactors and treated with anti-CD47, anti-PD1, or their combination which induced an IFN-γ response only in the tumor center, but not periphery [54]. Adding experimental support to the here described impaired activation signature in the tumor periphery.

By exploring the transcriptional trajectory of CD8⁺ T cells from the blood circulation into the immunosuppressive TME of the tumor periphery, we uncovered CX3CR1^high and CX3CR1^int effector and memory CD8⁺ T cells, respectively, to be highly enriched in the PBMC, but absent in the iTME. Recently, adoptive transfer studies of CX3CR1⁺ CD8⁺ T cells in a melanoma mouse model significantly suppressed tumor growth [41]. Others identified increased frequencies of CX3CR1⁺ CD8⁺ T cells in non-small cell lung and melanoma patients who responded to anti-PD-1 therapy, where these cells exhibited migratory capabilities into the tumor site followed by potent tumor rejection [40, 41]. Thus, the authors proposed T cell CX3CR1 expression as a predictor of response to ICI therapy. Even though CX3CR1⁺ CD8⁺ T cells might not be specific to GBM, as they are found in healthy donor and different inflammation related conditions [38–41], the absence of these effector and potentially ICI therapy responsive CD8⁺ T cells in the glioma TME could serve as an additional explanation for the disappointing outcome of clinical trials using ICI in glioma patients.

The observed T_rm exhaustion phenotype of the glioma residing CD8⁺ T cells was recently reported as well for tumor-infiltrating PD-1^high CD8⁺ T cells in hepatocellular carcinoma [44]. Whether these glioma-associated CD8⁺ T cells really possess tumor-specificity requires further study. Particularly in the light of a recent study by Smolders and colleagues who reported a consistent brain-resident CD8⁺ T cell population in a miscellaneous autopsy cohort of patients with neurological disorders excluding brain malignancies (Alzheimer’s disease, Parkinson’s disease, dementia, depression, multiple sclerosis), as well as patients with no known brain disease. These brain-resident CD8⁺ T cells displayed a remarkably consistent T_rm phenotype [55]. The authors further showed high expression of inhibitory receptors CTLA-4 and PD1 on the brain-resident CD8⁺ T_rm cells, which is in line with the core phenotypic signature of T_rm cells from other tissues [56, 57]. Yet, the brain CD8⁺ T_rm cells showed a preserved inflammatory potential with substantial production of IFN-ψ and TNF-α upon ex vivo stimulation. They concluded that extensive immune activation with release of highly neurotoxic lytic enzymes, such as perforin and granzyme B, harmfully impacts the brain parenchyma and should be tightly controlled, whilst maintaining the capability to elicit a fast inflammatory response when a neurotropic virus threatens the CNS [55]. Therefore, inhibitory receptors like PD-1 and CTLA-4 on brain CD8⁺ T_rm cells may support CNS homeostasis by preventing uncontrolled T cell reactivity, and the availability of the receptor ligands may determine their inhibitory effect. While this may represent a well-balanced equilibrium under healthy conditions, the tumor setting leads to its disruption with upregulation of inhibitory ligands like PD-L1 on glioma cells and CD86 on GAMs, leading to the dysfunctional state seen in the glioma-residing CD8⁺ T_rm cells.

Another study comprehensively showed, that peripheral infections generate antigenspecific CD8⁺ T_rm cells in the brain, mediating protection against CNS infections [42]. These data could implicate that the glioma-associated CD8⁺ T cells are devoid of tumor-specific reactivity, but rather represent a pre-existing T cell population generated after peripheral infections, which acquired a dysfunctional state upon glioma formation. To test this hypothesis, further characterization of these cells is required, including analysis of T cell receptor clonality and tumor-specificity by patient-matched T cell/glioma-sphere co-culture assays.

Lastly, our cell-cell interaction analysis revealed signaling pathways between peripheral MG and lymphocytes potentially inducing the observed impaired activation signature. The predicted interaction between MG SPP1 and NK cells integrin complex ITGA4-ITGB1 (CD49d-CD29), might mediate NK cell adhesion and migration [58]. This may facilitate interaction of inhibitory NK receptors KLRB1 and CD94/NKG2A with MG C-type lectin-related ligands and HLA-E, respectively, which could explain the observed impaired activation state of peripheral NK cells. Moreover, SPP1/CD44 interaction in T lymphocytes was described to suppress cell activation and proliferation [52]. In a comprehensive approach where transcriptional states of human MG were mapped during aging and disease, SPP1 was found to be differentially expressed in aging-microglia. This was associated to a doubling of the abundance of SPP1⁺ GAMs in GBM samples compared to MG from age-matched controls [59]. These observations support a possible role of MG SPP1 in glioma progression. More recently, myeloid-derived osteopontin (encoded by SPP1) was shown to trigger a chronic activation of NFAT2 in tumor-reactive CD8⁺ T cells, leading to T cell dysfunction and exhaustion in an experimental mouse model of GBM [60]. This adds an additional layer of mechanistic evidence for the SPP1/Osteopontin-mediated signaling axis leading to CD8⁺ T cell dysfunction in GBM.

Limitations of our study include the limited patient number, thereby our study was neither designed nor powered to explore differences in neoplastic cells, given the high inter- and intra-patient variability in glioma cells [12]. Importantly, our dataset establishes a starting point for further interrogation and provides a first analysis of the transcriptional landscape of the major immune populations in grade 4 glioma within three important regional compartments. Further, we confirmed the observed phenotype of CD8⁺ T cells in the blood and tumor periphery by flow cytometry in a cohort of ten additional patients, addressing possible generalization concerns. Together, we provide a regionally resolved map of leukocyte activation in the TME and blood circulation from grade 4 glioma patients, helping the research community to uncover novel therapeutic strategies to combat this fatal disease.

Methods

Ethics statement

Human adult GBM tissue samples were obtained at the Neurosurgical Clinic of the University Hospital of Basel, Switzerland in accordance with the Swiss Human Research Act and institutional ethics commission (EKNZ 02019-02358). All patients gave written informed consent for tumor biopsy collection and signed a declaration permitting the use of their biopsy specimens in scientific research, including storage in our brain tumor biobank (Req-2019-00553). All patient identifying information was removed and tissue was coded for identification.

Glioma tissue dissociation

Resected glioma tissue samples were immediately placed on ice and transferred to the laboratory for single cell dissociation within 2-3 h after resection. Human brain tissue was manually minced using razor blades and enzymatically dissociated at 37°C for 30 minutes with 1 mg/ml collagenase-4 (#LS004188, Worthington Biochemical Corporation, USA) and 250 U/ml DNAse1 (#10104159001, Roche, Switzerland) in a buffer containing Hank’s Balanced Salt Solution (HBSS) with Ca²⁺/Mg²⁺, 1% MEM non-essential amino acids (Gibco, USA), 1 mM sodium pyruvate (Gibco), 44 mM sodium bi-carbonate (Gibco), 25 mM HEPES (Gibco), 1% GlutaMAX (Gibco) and 1% antibiotic-antimycotic (Sigma-Aldrich, USA). Cells were filtered and separated from dead cells, debris and myelin by a 0.9 M sucrose (#84100, Sigma Aldrich) density gradient centrifugation. Upon ACK-lysis for removal of erythrocytes (#A1049201, Gibco) the now generated single-cell suspension (SCS) was washed, counted and frozen in Bambanker (#BB01, Nippon Genetics, Germany) in liquid nitrogen until use.

PBMCs (Peripheral blood mononuclear cells) preparation

Patient blood samples were directly placed on ice and transferred to the laboratory for PBMC isolation. Blood samples were centrifuged to separate buffy coat from plasma and erythrocytes, followed by standard density gradient centrifugation protocol (#17144002, Ficoll-Paque PLUS, Cytiva, USA) to isolate PBMCs. PBMCs were frozen in Bambanker (#BB01, Nippon Genetics, Germany) in liquid nitrogen until use.

FACS sorting for single cell RNA sequencing (scRNA-seq)

Cryopreserved tumor digests from glioma samples (center and periphery), as well as autologous PBMCs were thawed and washed with excess ice-cold 1xPBS and spun down at 350xg for 5 min. Subsequently, the cells were stained with Live/Dead (APC-Cy7 (Near IR), # L34976, Thermo Fischer) and a cocktail of fluorescently-conjugated antibodies CD11b (FITC, clone M1/70, #101206, BioLegend) and CD45 (FITC, clone 2D1, #368508, BioLegend), and large debris were removed with a 40-μm strainer. All samples were acquired on the BD FACS ARIA Fusion III (Becton Dickinson GmbH, Germany). For single-cell RNA-seq experiments, live and single gated cells were sorted into non-immune cell (CD45^-CD11b^-) and immune cell (CD45⁺CD11b⁺) populations. Both populations were directly sorted into Eppendorf tubes with 1xPBS supplemented with 1% BSA for single cell RNA sequencing.

Single cell RNA sequencing (scRNA-seq) – Library preparation and sequencing

Single-cell RNA-seq was performed using Chromium Single Cell 3’ GEM, Library & Gel Bead Kit v3 (#CG000183, 10x Genomics, Pleasanton, CA, USA) following the manufacturer’s protocol. Briefly, non-immune cells and immune cells were mixed at a defined ratio of 1:4. Roughly 8000-10000 cells per sample, diluted at a density of 100–800 cells/μL in PBS plus 1% BSA determined by Cellometer Auto 2000 Cell Viability Counter (Nexelom Bioscience, Lawrence, MA USA), and were loaded onto the chip. The quality and concentration of both cDNA and libraries were assessed using an Agilent BioAnalyzer with High Sensitivity kit (#5067-4626, Agilent, Santa Clara, CA USA) and Qubit Fluorometer with dsDNA HS assay kit (#Q33230, Thermo Fischer Scientific, Waltham, MA USA) according to the manufacturer’s recommendation. For sequencing, samples were mixed in equimolar fashion and sequenced on an Illumina HiSeq 4000 with a targeted read depth of 50,000 reads/cell and sequencing parameters were set for Read1 (28 cycles), Index1 (8 cycles), and Read2 (91 cycles).

Single cell RNA sequencing (scRNA-seq) – Computational analysis

The dataset was analyzed by the Bioinformatics Core Facility, Department of Biomedicine, University of Basel. Read quality was controlled with the FastQC tool (version 0.11.5). Sequencing files were processed using the Salmon Alevin tool (v 1.3.0) [61] to perform quality control, sample demultiplexing, cell barcode processing, pseudo-alignment of cDNA reads to the human Gencode v35 reference and counting of UMIs. Parameters --keepCBFraction 1 and --maxNumBarcodes 100000 were used.

Processing of the UMI counts matrix was performed using the Bioconductor packages DropletUtils (version 1.8.0) [62, 63], scran (version 1.16.0) [64, 65] and scater (version 1.16.2) [66], following mostly the steps illustrated in the OSCA book (http://bioconductor.org/books/release/OSCA/) [65, 67]. Filtering for high-quality cells was done based on library size (at least 2,000 UMI counts per cell), the number of detected genes (at least 700 genes detected) and the percentage of reads mapping to mitochondrial genes (larger than 0% and lower than 15%), based on the distribution observed across cells. Low-abundance genes with average counts per cell lower than 0.006 were filtered out. The presence of doublet cells was investigated with the scDblFinder package (version 1.2.0), and suspicious cells were filtered out (score > 0.6). After quality filtering, the resulting dataset consisted of UMI counts for 15,523 genes and 45,466 cells, ranging from 803 to 9,121 per sample.

UMI counts were normalized with size factors estimated from pools of cells created with the scran package quickCluster() function [64, 68]. To distinguish between genuine biological variability and technical noise we modeled the variance of the log-expression across genes using a Poisson-based mean-variance trend. The scran package denoisePCA() function was used to denoise log-expression data by removing principal components corresponding to technical noise. Consistent to the findElbowPoint() function from the PCAtools Bioconductor package, this led us to retain the top 5 principal components (PCs), explaining 48.3% of the total variance in the dataset, for the clustering and dimensionality reduction steps.

Since the excluded deeper PCs (PC8 and 10) were associated to patient-specific effects (which in our experimental setup are also confounded with batch effects), this choice dispensed us from performing an additional correction for patient-specific effects. A quantitative test of overlap of cells across patients was made using the CellMixS package [69] (version 1.12.0) developed to quantify the effectiveness of batch correction methods. We used the cell-specific mixing (CMS) score, which highlighted a very good overlap across cells from different patients in the lymphoid compartment and for the monocytes. The myeloid compartment displayed a slightly elevated patient-specific structure, while it was most pronounced for the CD45 negative subset.

A t-stochastic neighbor embedding (t-SNE) was built with a perplexity of 50 using the top most variable genes (141 genes with estimated biological variance > 0.3, excluding genes with highest proportion of reads in the ambient RNA pool estimated from empty droplets), and the denoised principal components as input. Clustering of cells was performed with hierarchical clustering on the Euclidean distances between cells (with Ward’s criterion to minimize the total variance within each cluster [70]; package cluster version 2.1.0). The number of clusters used for following analyses was identified by applying a dynamic tree cut (package dynamicTreeCut, version 1.63-1) [71], resulting in 10 with argument deepSplit set to 1, or 22 clusters with argument deepSplit set to 2. This clustering was validated with an alternative clustering approach using a graph-based approach and Louvain algorithm for community detection (using the Find-Neighbors() and FindClusters() functions from the Seurat package [72] (version 4.3.0), with a resolution of 0.7) which showed a good correspondence to our hierarchical clustering.

The findMarkers function of the scran package was used to find the best markers across annotated cell types (parameters direction=“up” and pval.type=“any”). The top 10 markers for each cell type were extracted and pooled to from a list of 68 markers.

The Bioconductor package SingleR (version 1.2.4) was used for cell-type annotation of the cells [73] using as references (i) a public bulk RNA-seq dataset of sorted immune cell types from human PBMC samples [16], available through the celldex Bioconductor package; (ii) a bulk RNA-seq dataset of sorted immune cell types from the tumor microenvironment of human gliomas [5] (UMI count matrix and annotation downloaded from https://joycelab.shinyapps.io/braintime/); (iii) a 10X genomics scRNA-seq dataset of TAMs from the tumor micro-environment of glioblastoma tumors from 7 newly diagnosed human patients [9]; (iv) a microglia and a macrophage signature scores were defined by averaging the center and scaled expression levels of gene lists obtained in [7] and (v) a Smartseq2 scRNA-seq dataset of IDH-wild-type glioblastoma tumors [17] (downloaded from GEO accession GSE131928). An endothelial score was defined by averaging the center and scaled expression levels of the genes CDH5, VWF, CD34 and PECAM1.

The SingleR high-quality assignments (pruned scores) from the comparisons to the multiple references (some being better at enlightening some immune subsets than others) were used to manually derive a consensus cell type annotation for each cluster. This manual annotation was done jointly with the input from the signature scores and relative expression patterns of known marker genes and of cluster-specific genes.

When we judged that the resolution was not sufficient to discriminate well the subtypes of cells within a subset (in particular for the CD45neg and myeloid cells subsets), we isolated these cells and re-performed selection of hyper-variable genes, dimensionality reduction and clustering on the isolated subset.

Differential abundance analysis of the identified cell types between tumor sites was performed using limma-voom [74]. This method, implemented in the diffcyt package, was shown to be able to handle well the overdispersion in the cell type proportions estimates typically observed with single cell technologies [75]. Differential abundance of cell types was considered to be significant at a false discovery rate (FDR) lower than 5 %. To provide increased resolution into local trends within the tumor innate immune subset we used the miloR bioconductor package [76] (version 1.4.0) to test for differential abundance on neighborhoods of a k-nearest neighbor graph. Testing was done with limma-voom similarly to above. Cell type labels were assigned to each neighborhood by finding the most abundant cell type across cells and requiring that at least 70% of the cells come from one cell type.

Differential expression between tumor sites, or between PBMC cells and tumor periphery cells, stratified by annotated cell type, was performed using a pseudo-bulk approach, which was shown to be the best-performing approach for differential expression testing in recent benchmarks [77]. UMI counts were summed across cells from each cell type in each sample when at least 20 cells could be aggregated. The aggregated samples were then treated as bulk RNA-seq samples [78] and for each pairwise comparison genes were filtered to keep genes detected in at least 5% of the cells aggregated. The package edgeR (version 3.30.3) [79] was used to perform TMM normalization [80] and to test for differential expression with the Generalized Linear Model (GLM) framework, using a model accounting for patient-specific effects. Genes with a FDR lower than 5 % were considered differentially expressed. To validate our differential expression analysis, we rerun a pseudo-bulk expression analysis using DESeq2 (version 1.38.3) [81] for the main immune populations in the periphery versus center and periphery versus PBMC comparison. Although, the results were not identical, they overall agreed with edgeR (Supplementary File 7). Gene set enrichment analysis was performed with the function camera [82] on gene sets from the Molecular Signature Database (MSigDB, version 7.4) [27, 83]. We retained only sets containing more than 5 genes, and gene sets with a FDR lower than 5% were considered as significant.

Cell chat analysis

The R package CellChat (1.1.3) [50] was used to analyze cell-cell interactions in our dataset (with previously annotated 9 cell types). We followed the recommended workflow to infer the cell state-specific communications (using identifyOverExpressedGenes, identifyOver-ExpressedInteractions and projectData with the default parameters). We performed 3 separate analyses, on the center and the periphery subsets and a comparison analysis as described in the official workflow. We visualized the significant interactions for the microglia cluster using netVisual_chord_gene and used plotGeneExpression to display of the expression of all genes involved SPP1 signaling pathway in the cell populations. Finally, netAnalysis_signalin-gRole_scatter was used to calculate and visualize incoming and outgoing signaling strength.

Flow cytometry analysis of paired PBMC and periphery samples

Cryopreserved samples were thawed and washed with excess ice-cold 1xPBS and spun down at 350xg for 5 min. Cells were resuspended in FACS buffer (PBS plus 2% FBS) and blocked with monoclonal antibody to CD16/32 (Human TruStain FcX, #422302, Biolegend) for 10 min at 4°C before staining with surface antibodies: CD45RA (PE, clone HI100, #304108), CD45RO (APC, clone UCHL1, #304210), CD3e (BV650, clone UCHT1, #300468), CD8a (BV421, clone RPA-T8, #301036), CCR7 (FITC, clone G043H7, #353216), CD62L (AF700, clone DREG-56, #304820), CD69 (APC-Cy7, clone FN50, #310914), CD103 (PerCP/Cy5.5, clone Ber-ACT8, #350226) and CX3CR1 (PE/Cy7, clone 2A9-1, #341612). All antibodies were purchased from BioLegend, USA. Cells were stained for 30 min at 4°C, and subsequently washed with FACS buffer. To exclude dead cells Zombie Aqua Fixable Viability Kit (#423102, 1:100, BioLegend) was added. Acquisition was performed on a CytoFLEX (Beckman). Data was analyzed using FlowJo software, version 10.8.1 (TreeStar). Gates were drawn by using Fluorescent Minus One (FMO) controls.

Statistical analysis of flow cytometry data

Data analysis and graph generation was performed using GraphPad Prism 9 (GraphPad Prism Software Inc.). Paired comparisons between two groups were performed using Wilcoxon matched-pairs signed rank test. Differences of more than two paired groups were assessed using repeated measures one-way ANOVA test, followed by post-hoc Šidák’s multiple comparisons correction. A p value < 0.05 was considered statistically significant. *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001.

Graphical illustrations

All graphical illustrations were created with BioRender.com.

Data availability

The UMI count matrix and cell metadata from the scRNA-seq dataset are available on GEO under accession number GSE197543. The remaining data are available within the Article and Supplementary Information.

Acknowledgements

We are grateful to the patients and their families for their consent to donate tissue to our brain tumor biobank. We thank Tamara Hüssen, Alison Riberio and Florian Limani for experimental support and Cecile Buenter for critical revision of the manuscript. Calculations were performed at sciCORE (http://scicore.unibas.ch/) scientific computing center at the University of Basel.

Author contribution

G.H., S. Herter, M.B., S.D. and P.S., conceived and planned the project. J.R. and S. Hogan analyzed the scRNA-sequencing data and performed statistical analysis. P.S. and S.D. performed experiments and interpreted the results. T.S., T.A.M. and M.-F.R. helped coordinating experiments. P.S. wrote the manuscript. G.H. supervised and coordinated the study and critically revised the manuscript. All authors reviewed the paper and approved its final version.

Competing interests

G.H. has equity in, and is a cofounder of Incephalo Inc. S.D., S.H. and M.B. are affiliated with Roche Pharmaceutical Research and Early Development. The authors have no financial interests to declare.