Spatial immune profiling of glioblastoma identifies an inflammatory, perivascular phenotype associated with longer survival

The historic term “glioblastoma multiforme” reflects the heterogeneity of histopathologic compartments in malignant gliomas [1]. Prominent and disease-defining histopathologic features of glioblastoma, the most common type of glioma, include microvascular proliferation and necrosis, the latter often with perinecrotic palisading of tumor cells [3]. These features are thought to result from a vicious cycle between dysfunctional, clotting blood vessels and perinecrotic areas where hypoxia and low pH drive an aberrant pro-angiogenic response, resulting in more dysfunctional blood vessels [6]. The interconversion between both compartments results in spatially resolved niches, where metabolically driven myeloid cell and tumor cell subpopulations interact [5]. Here, we have analyzed spatial immune profiles from clinically and molecularly well-annotated glioblastoma patient samples to explore clinical implications of spatial immune phenotyping, including the abundance of immunotherapy targets, and associations of spatial immune profiles with outcome. Patients and methods are detailed in supplementary Note S1. First, we analyzed 360 anatomically defined glioblastoma regions of interest by spatial immune profiling utilizing a panel of 28 DNA bar-coded antibodies to quantify immune cell types and immunotherapy targets in perivascular, perinecrotic and cell-dense tumor regions (Fig. 1a, Table S1). Regions of interest were defined in 30 tissue sections from 20 patients diagnosed with glioblastoma according to the 2021 World Health Organization classification of central nervous system tumors [3], including paired samples from untreated primary tumors at diagnosis and recurrent tumors after standard chemoradiotherapy with temozolomide (paired sample cohort, PSC) of 10 patients from the central nervous systems tumor tissue bank Dusseldorf, and 10 samples from untreated primary tumors of patients included in the EORTC 1419 ETERNITY study, who survived for more than 5 years after diagnosis (long-term survivors, LTS, Figure S1, Table S2, [2]). Enrichment of antibody-linked DNA barcodes in perivascular and perinecrotic compartments was analyzed utilizing cellular tumor regions as reference (Fig. 1b). In the perivascular zone of LTS, there was enrichment of CD8 (p = 0.002) and of the innate immunity activator, stimulator of interferon genes (STING, p = 0.019). By contrast, the immunosuppressive macrophage marker CD163 was enriched in the perivascular zone among newly diagnosed tumors from non-LTS patients (p = 0.016). In


Patient tissue samples
Formalin-fixed and paraffin-embedded tissue sections with confirmed diagnosis of an IDH wild-type glioblastoma according to the 2021 World Health Organization classification of central nervous system tumors [8] were utilized for GeoMx digital spatial immune profiling [9].Sections were selected by a board-certified neuropathologist (JF) based on the presence of perivascular and perinecrotic compartments.Patients provided written informed consent for the use of their tissue samples and clinical data for research purposes.Matched primary and recurrent tissue samples were derived from the central nervous system tumor tissue bank Düsseldorf from a previous research project that was approved by the Ethics Committee of the Medical Faculty at Heinrich Heine University Düsseldorf (study number 4940, [7]).Tissue samples of long-term surviving patients were derived from the EORTC 1419 study, which was approved by the ethics board of the Canton of Zurich (study number KEK 2014-0555) and locally at each participating site [4].
Slides were then scanned utilizing a GeoMx Digital Spatial Profiler (Nanostring Inc) and regions of interest selected remotely by a board-certified neuropathologist (JF).

RNA sequencing datasets
The Ivy glioblastoma atlas project (GAP) spatial RNA sequencing core dataset, comprising 120 spatial transcriptomes from 10 patients with IDH wild-type glioblastoma [13] was utilized for gene ontology (GO) network analyses of differentially expressed immunity genesets (Figure 1c) and for the development of a spatial immunity classifier (Figure 1d).The spatial immunity classifier was then applied to bulk RNA sequencing datasets of IDH wild-type gliobastoma samples, including N=139 newly diagnosed glioblastomas from The Cancer Genome Atlas pan-cancer study [5] (Figure 1e), a merged dataset of N=72 recurrent glioblastomas [14] (Figure 1f, Figure S6a) as well as two smaller studies comprising N=36 [7] (Figure S6b) and N=29 [2] (Figure S6c) recurrent glioblastomas, respectively.

Geneset network analyses
Geneset network analyses were conducted in the Ivy GAP dataset to identify differentially expressed genesets in the perivascular zone versus cellular tumor and in the perinecrotic zone versus cellular tumor, respectively (Figure 1c [11,12]).First, differential expression of GO genesets [1,3] with a size < 200 genes was analyzed between two groups of interest.Genesets with a false discovery rate < 0.05 were then utilized for network discovery by depicting differentially expressed genesets as nodes that were connected by edges if they shared at least one gene.Genesets that were up-regulated in the perivascular or perinecrotic zone compared to cellular tumor were colored red, genesets that were down-regulated were colored blue, and genesets in which some genes were up-and others down-regulated were colored gray.Doublets and community-connecting edges were removed to split large clusters into smaller ones based on edge betweenness.Gene set names were split into individual words and counted to automatically construct cluster labels.

Classification of spatial immunity gene expression patterns
A random forest algorithm to discriminate spatial gene expression patterns based on N=72 pan-caner immunity-associated gene sets [13] was trained on N=80 randomly selected spatial RNAseq samples from the Ivy GAP dataset and was tested on N=40 samples from the same dataset [12].Weighted correlation network analysis delivered 3 coexpression clusters of immune genesets characterizing the different glioma regions (eigengenesets).Clustering of all N=120 spatial RNAseq samples based on the expression of these eigengenesets revealed associations with the perivascular, perinecrotic and infiltration zone samples, whereas cellular tumor was characterized by lack of expression of these clusters (Figure S5).Eigengeneset clustering indicated overlap of the cellular tumor expression pattern with the other three patterns.The GO terms constituting eigengeneset clusters were imputed to the CIBERSORT online tool to classify bulk RNAseq datasets (Figure S4).On note, the lack of genesets specific for classifying cellular tumor led to an intended focus on the relative distribution of perivascular, perinecrotic and infiltration zone immunity patterns.No further filtering methods were applied for this purpose.

Survival analyses
Survival analyses were performed in R using the coxph function from the survival package (Version 3.4.0).Scores from the random forest-based classifier were used to automatically devise a cutoff using the surv_cutpoint function from the survminer package (Version 0.4.9) and then used as input for the coxph function.For the newly diagnosed glioblastoma cohort [3], MGMT promoter methylation status was also used as input, whereas limited availability of MGMT promoter methylation status precluded respective adjustments in the recurrent glioblastoma cohorts [2, 7,14].

Biostatistics tools
All statistical analysis was performed using the R environment for statistical computing (R version 3.5.1,2018-07-02) with Bioconductor [6] and dedicated packages.All analysis and this report is produced programatically using R markdown (https://CRAN.R-project.org/package=rmarkdown) in Rstudio (RStudio Inc, Boston, MA, USA), compliant with the principles of Reproducible Research [10].

Note S2. Funding sources.
This study was funded by grants from the OPO foundation, the Desirée und Niels Yde foundation, Charlotte and Nelly Dornacher Foundation, the Swiss National Science   CIBERSORT was utilized for classifier application and a cut-off of 0.5 was employed to define infiltration zone (IFZ) high vs low scores; the logrank test was utilized for curve comparison; datasets: a, ref. [14], N=72; b, ref. [7], N=36; c, ref.

Figure S1 .Figure S2 .
Figure S1.Overall survival of spatial immune profiling cohorts of patients with

Figure S3 .
Figure S3.Spatial digital cytometry.a, immune cell type abundance in the Ivy GAP