A spatially resolved human glioblastoma atlas reveals distinct cellular and molecular patterns of anatomical niches

Domain: chemistry
Published: 2026-02-20
Editorial outcome: minor_revision

Abstract
Glioblastoma is an aggressive brain cancer with limited treatment options and poor patient survival, driven in part by cellular diversity within tumors. While individual cell types have been catalogued, how malignant, vascular, and immune cells are spatially organized inside human tumors remains incompletely understood. Here we show a spatially resolved, multi-modal atlas of human glioblastoma that integrates gene expression profiling across tissue sections with matched single-cell and protein measurements at subcellular resolution. Using a targeted 348 gene panel enriched for vascular and stromal markers, we identify less well-characterized endothelial, perivascular, and fibroblast-like cell states and define their spatial associations with malignant and immune compartments. We further identify a distinct oligodendrocyte population restricted to tumor core and perivascular regions that exhibits gene expression patterns associated with tumor recurrence and poor clinical outcome. This publicly accessible atlas provides a high-resolution framework for studying the spatial organization of glioblastoma and highlights region-specific cellular interactions that may represent therapeutically actionable vulnerabilities.

Introduction
Glioblastoma (GBM) is the most lethal form of primary brain tumor in adults, with a median survival of only 15 months despite intensive treatment regimens, including surgical resection, radiation, and chemotherapy1. The inevitable recurrence of GBM and its resistance to therapies are largely driven by profound cellular and molecular heterogeneity2. Understanding how the tumor heterogeneity arises and contributes to GBM progression is crucial for improving patient outcomes and developing more effective treatments.

Early efforts in comprehensive genomic profiling suggested three molecular subtypes of GBM: classical, mesenchymal, and proneural3. Subsequently, single-cell RNA sequencing (scRNA-seq) has further elucidated the molecular diversity of GBM, revealing the coexistence of multiple malignant cell states such as mesenchymal-like (MES-like), neural progenitor cell-like (NPC-like), astrocyte-like (AC-like), and oligodendrocyte precursor cell-like (OPC-like) within individual tumors4–7. Transcriptional heterogeneity in tumor tissue is complex and shaped by cell-intrinsic as well as cell-extrinsic factors. Although scRNA-seq has advanced our understanding of GBM cell-intrinsic heterogeneity, it remains unclear whether tumor cell states are randomly dispersed throughout the tissue or instead shaped by regional microenvironmental pressures.

It has long been recognized that GBM exhibits significant histological heterogeneity, with distinct regions of the tumor displaying different patterns of cellular organization8. Hans-Joachim Scherer’s pioneering work on the spatial architecture of gliomas highlighted the importance of these spatially organized patterns in the progression and invasiveness of gliomas9. In our previous work, we reported the Ivy Glioblastoma Atlas (Ivy GAP) study, a comprehensive molecular-histology map of GBM designed to address the tumor’s spatial and molecular heterogeneity10. Building on the Ivy Glioblastoma Atlas (Ivy GAP), which linked molecular programs to GBM histological features, recent advances such as 10x Visium, Xenium, and GeoMx now provide finer resolution, enabling spatial mapping of tumor–microenvironment interactions at the cellular scale11–15. These studies showed that the tumor microenvironment, in particular hypoxia and immune cell infiltrates, reshapes the transcriptional and genetic heterogeneity of glioblastoma cells.

A major finding of recent scRNA-seq analyses was that unmutated oligodendrocytes are a part of the glioma microenvironment, infiltrating gliomas at substantial and varying frequencies across glioma subtypes and patient ages. Although the most frequent type of non-tumor cells in some adult high-grade glioma, surprisingly little is known about their spatial distribution. Ravi et al. have documented spatial segregation of neurodevelopmental neighborhoods from hypoxia-associated mesenchymal and inflammation-associated neighborhoods13. Greenwald et al. documented 14 metaprograms, including 6 non-malignant programs, which feature Oligo (oligodendrocytes) alongside Mac (macrophages/microglia), Reactive-AST (reactive astrocytes)15. The Oligo metaprogram was localized in Layer 5, corresponding to the infiltrating area of the brain parenchyma. Whether OL/OPCs exhibit temporally distinct and functional states in glioma, similar to normal OPCs which display spatiotemporal heterogeneity, depending on age and their anatomic location, remains unknown. While past studies have significantly enhanced our understanding of the microenvironmental factors driving glioma heterogeneity, their limited resolution reduces the capacity to fully capture the cellular composition and complex patterns of cell-cell interactions within GBM anatomical subdomains.

In this study, we present a comprehensive, high-resolution mapping of the cellular architecture within the anatomical niches of GBM. We generated whole spatial transcriptome data using Visium CytAssist, integrated with patient-matched single cell multi-omic data, including from matched tumor core and periphery. In addition, for better analysis of tissue architecture, we acquired high-resolution H&E images with spatial transcriptome at 40x in contrast to the more commonly used 5x or 10x magnifications. Additionally, we generated high-depth Spatial Transcriptome (ST) data, leading to a median of 4694 genes detected per section, compared to 2159 in a recent study and identifying an additional extra 482 genes, providing a more comprehensive view of the molecular landscape13. Furthermore, previous studies have primarily used fresh frozen (FF) sections for ST data generation13,15 however, we generated ST data from not only FF (n = 5), but also FFPE sections (n = 27), which are more commonly used in clinical samples, offering broader applicability while preserving tissue architecture for comprehensive analysis. As a result, we obtained 115,914 spatial transcriptomes across 32 genotyped glioma sections, mapping 56 distinct cellular components, including oligodendrocyte subtypes not emphasized in prior studies, some localized to distinctive spatial domains in the tumor core. We identified a neutrophil population from the freshly processed samples, a cell type missing from most human GBM scRNA analyses16. In addition, we employed Xenium in situ technology with subcellular resolution to refine cell type annotations and validate the existence of the rare oligodendrocyte subtype. This study underscores the extensive spatial and cellular heterogeneity of GBM and emphasizes the intricate interplay among diverse cellular components. To facilitate broader access and utilization, our analytical results are presented as a public resource in a user-friendly format, available for all at https://github.com/nameetas/TSKGA.

Results
Spatially resolved GBM atlas
While considerable progress has been made in exploring molecular GBM tumor heterogeneity with spatial transcriptomics (ST) technologies, high-resolution cellular organization and cell–cell interaction analysis is still lacking. To generate a high-resolution spatial map of the cellular structure across the anatomical regions of GBM, we employed single cell, whole transcriptome ST (Visium CytAssist) as well as subcellular resolution in situ ST (Xenium) technologies on both fresh frozen (FF) and formalin-fixed paraffin embedded (FFPE) preserved blocks (Fig. 1A–E). To comprehensively capture GBM heterogeneity, we analyzed patients across different age groups with diverse genomic alterations (Supplementary Data 1 and Fig. 1F) and profiled ST on at least two tissue sections prepared from different locations with distinct anatomical features within the same tumor.

We investigated the cellular architecture of GBM in 28 patients (detailed patient and specimen information is provided in Supplementary Data 1 and Fig. 1F), including 18 with GBM, 5 with astrocytoma (IDH-mutant), 1 with oligodendroglioma (IDH-mutant, 1p19q co-deleted), and 4 with miscellaneous CNS pathologies (Fig. 1A). Complementary single-cell RNA sequencing (scRNA-seq) and cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) data were acquired from 25 patients, resulting in a collective analysis encompassing 223,113 cells (Fig. 1B). In addition, we cataloged the GBM single-cell reference atlas by compiling the newly generated single-cell data with publicly available sc/snRNA-seq datasets, containing 10 major cell types, which were further classified into 58 transcriptional states (Fig. 1B). Furthermore, we collected ST data for 17 cases, comprising 16 brain tumor patients (13 GBM, 3 Astrocytoma) and one healthy brain tissue control using the 10x Genomics Visium ST platform (Fig. 1C). As a result, we successfully generated a comprehensive dataset comprising 115,914 individual ST spots from 32 tissue slices. Among these, 14 cases have both matched ST and single-cell data. By integrating the reference single cell atlas with ST data, we identified recurrent molecular features associated with anatomical features (AF)-based niches including cell-type composition, and co-occurrence (Fig. 1D). Additionally, we generated subcellular resolution ST data using recently developed Xenium in situ technology (Fig. 1C) for 4 GBM patients, including a custom panel of 348 genes. From this, we annotated 729,307 cells within an area of 155mm2. Collectively, we generated an integrated spatially resolved GBM atlas which is publicly available through an interactive website at https://github.com/nameetas/TSKGA to enhance user accessibility (Fig. 1E).

GBM cellular complexity revealed through matched core and periphery single-cell analysis
To investigate cell states in both the core and periphery (peri) of GBM, we generated single-cell transcriptomes from six matched core–peri fresh samples, processed on the same day (Supplementary Data 1). After quality filtering, 223,113 cells were pooled for single-cell transcriptomics, batch-corrected, and clustered using Recursive Consensus Clustering (RCC)16. The level 1 clusters of RCC were used to assign major cell types (microglia, oligodendrocyte, neoplastic, myeloid, lymphocyte) by canonical marker gene expression (Supplementary Data 2). The major cell types were then recursively clustered up to three levels and markers were identified for each cell subtype or cellular state (Fig. S1A–S1F). The cluster name indicates the major cell type followed by the cluster numbers separated by dashes for each level or an abbreviation based on previously described cellular states (see Methods, Fig. 2A, Supplementary Data 3). Notably, we were able to detect a granulocyte population from the freshly processed samples (4.6% vs 1% in freshly frozen), a cell type missing from most human glioblastoma scRNA-seq datasets unless captured using alternate protocols21,26. Since, we integrated these cell types in our reference dataset we were able to see that Chromatin.reg (malig.) program correlates to AstroSERPINI2, Prolif.Metab (malig. and shown to be strongly associated with vascular program in Greenwald et al.) to gbmEndoPeri_2, and MES.Ast (malig. and shown to be associated with immune activity in Greenwald et al.) to gbmFib_5 (Supplementary Data 6, Supplementary Data 13).

To investigate the spatial organization of GBM anatomical features (AFs), we generated ST data comprising 32 tumors from 17 patients, including 13 patients diagnosed with IDH-wildtype GBM, 3 patients with IDH-mutant high-grade astrocytoma, and a normal brain region serving as a comparative control (Supplementary Data 1, Fig. 1F). The ST specimens were carefully chosen to encompass diverse AFs of GBM identifiable by H&E staining; leading edge (LE), cellular tumor (CT), microvascular proliferation (MVP) and pseudopalisading cells around necrosis (PAN)10. We first defined AFs of each ST spot. We utilized previously established marker gene sets for all AFs from Ivy GAP10,13 (Fig. 3B; Supplementary Data 7). Besides the histologically visible AFs, we also used glioma stem-like cell (GSC) niches-vascular (HBV_ccGSC - hyperplastic blood vessels with community curated glioma stem-like cell markers (ccGSC)) and hypoxic (PNZ_ccGSC - perinecrotic zone with ccGSC expression)-defined by integrating histological traits and transcriptional profiles of laser-captured microdissected regions expressing ccGSC markers10,13. Tumor regions lacking expression of any ccGSC markers were categorized as CT_control. We further divided the LE area into LE_WM (white matter) and LE_GM (gray matter) as these areas have distinct cellular profiles.

Each spot was scored for the 56 reference cell states using averaged marker expression, followed by binary presence/absence assignment (Fig. 3C, Supplementary Data 8 see Methods, Supplementary File 14). Visualization of Ivy GAP AFs over high-resolution H&E images indicated that further refinement was needed for AF assignment as HBV, MVP did not always overlap with hyperplastic/glomerulus blood vessels and PAN did not show any visible PAN on H&E. This can largely be attributed to overlapping immune and vascular signature across HBV, MVP, PNZ and overlapping hypoxic signature between PNZ and PAN. We were able to use our cell type assignment and H&E features to further stratify AFs into 10 AFs with additional introduction of BV (Blood Vessel), and IGN (Immune-Glial Niche) (Fig. 3D–E, Fig. S3). Notably, normal sample spots (SNU38) were classified into LE_GM, LE_WM, and BV AFs, while spots from patients with astrocytoma (Ast), IDH mutant (IDHm) (SNU19, SNU22) were classified into LE_GM, LE_WM, BV, IT, and CT (Fig. 3F, Fig. S3). In stark contrast, GBM patient spots exhibited classification across all AFs (Fig. S3). We performed copy number analysis as an alternate method to estimate the presence of tumor cells for each spot. We also calculated the number of major cell types in each spot using our binary assignment. Furthermore, cell counts per spot were determined by identifying nuclei on H&E images using the DINO object detection algorithm. As expected, CT spots showed maximum CNV scores, highest cell density but low cell type diversity; IGN, HBV, and MVP followed by PNZ show lower CNV scores due to presence of nonmalignant cell types but highest cell type diversity (Fig. 3G, Fig. 4A).

Mapping GBM tumor architecture by anatomical features specific cell states
In order to assess spatial AF preference for a given cell subtype, we compared the distribution of normalized expression values for that cell subtype in real data vs. AF randomized data. If the area under the curve for actual data is higher than randomized data then the cell subtype is preferentially present in that AF (Fig. 4B, Fig. S4). The distribution of the 56 cellular subtypes in 115,914 ST spots from 32 specimens, as well as detailed specimen analysis, including matched single cell UMAP, AF annotation, and spatial CNV analysis, are provided on the interactive website, https://github.com/nameetas/TSKGA (Supplementary Data 15).

Spatial annotation of cellular subtypes revealed varying degrees of spatial and cellular heterogeneity within the GBM tumor microenvironment (TME) (Fig. 4C and Fig. S4). BV are enriched with Endo_cap cells and consist mainly of normal vasculature, predominantly in LE_GM, followed by LE_WM, IT, and CT regions. LE_GM is enriched in neurons, AstroPLCG1 astrocytes, and Art2 cells, while LE_WM contains more oligodendrocytes, particularly Oligo_2_1. As the region transitions from LE_WM to IT, there is increased presence of other oligodendrocytes (Oligo_2_3_1), reactive astrocytes (AstroAQP1), microglia with an emergence of non-mesenchymal tumor cells. CT shows a clear enrichment of non-mesenchymal tumor cells. The Immune-Glial Niche (IGN) contains elevated levels of GBM-associated vascular cells, reactive astrocytes, fibroblasts, lymphoid and myeloid cells, as well as mesenchymal tumor cells (TC_mesnh). The identified Oligo_2_3_2 subtype is also found in IGN. HBV and MVP regions show the highest enrichment of vascular cells, with MVP notably containing more GBM-associated pericytes (gbmPeri, gbmEndoPeri). The PNZ is a hypoxic zone rich in mesenchymal tumor cells (TC_mesnh, TC_mesh), immune cells, and vascular cells. In contrast, the PAN-also hypoxic-features TC_mesh, TC_NPC, and hypoxic myeloid cells, and appears as pseudopalisades on H&E staining. The immune cell profile changes from IGN to HBV to MVP: microglia decrease while macrophages and lymphoid cells increase. In the hypoxic niches (PNZ and PAN), immune diversity drops further, with only macrophages remaining (Fig. 4C).

GBM landscape visualized at sub-cellular resolution
The spatial resolution of Visium spatial transcriptomics (55µm) enables the detection of gene expression patterns across tissue sections but has inherent limitations in resolving individual cells. To complement this, we employed the Xenium platform, which provides subcellular-resolution transcriptomic data, enabling spatial profiling at the level of individual cells. We designed a custom panel of 348 genes representing all major cell types (Methods, Supplementary Data 10, Fig. S5A) and used it to generate Xenium data from four GBM patients (SNU21, SNU33, SNU18, and SNU25) within our Visium cohort (Fig. 5A). This panel enabled spatial validation of 10 major cell types and their associated states, including previously uncharacterized subtypes of endothelial cells, pericytes, and oligodendrocytes.

We first annotated all cells using a semi-automated process, resulting in 729,307 annotated cells across a total area of 155mm² (Fig. 5B–D). Gene expression profiles from the spatial single-cell data showed excellent correlation with our reference single-cell dataset using the 348 panel genes (Fig. 5E).

The Xenium data enabled analysis of anatomic features (AFs) across continuous regions spanning up to 1.2cm. We selected 11 representative trajectories to capture transitions observed in the data (Fig. S5B). Grids were defined to match the area of a Visium spot (2376μm²). A subset of 40 genes (Supplementary Data 11) was chosen to assign AF labels to these grids based on known AF-specific cellular compositions (Fig. S5C). Spatial gene expression profiles were used to assign AFs to trajectory grids, with validation by H&E image visualization (Fig. 6 and Fig. S6). Representative areas for each AF with spatial transcripts overlaid on H&E image are shown in Fig. 6B.

Among these, IGN showed two distinct gene expression profiles. A subset of IGN grids displayed high expression of fibroblast and immune cell markers but lacked glial cell signatures. These were labeled as IFN (Immune Fibroblast Niche), corresponding to fibrotic scar. Fibrotic scarring was found to contribute to an immune-suppressive and pro-tumorigenic microenvironment that facilitates glioblastoma regrowth after therapy27. Remaining grids were assigned AF labels based on their shortest distance from the annotated trajectory grids (Fig. 7A).

We then assessed the cell type composition of AF-assigned grids (Fig. 7B–C). As expected, the cellular composition of AFs was consistent with findings from the Visium data (Fig. 7B–C). AF neighborhood analysis also revealed similar spatial relationships to those observed in the Visium dataset (Fig. 7D).

Distinct spatial and molecular characteristics of oligodendrocytes in GBM
Emerging evidence indicates the pivotal involvement of oligodenglia in various neurodegenerative conditions21,22. However, the precise extent of oligodendrocyte diversity and their role in GBM pathology remains elusive. In this study, we identified four distinct oligodendrocyte cell states within GBM (Fig. 2F) and delineated markers characterizing these states (Fig. S1A). Notably, Oligo_2_3_2 was uniquely detected in GBM, but not in IDH-mutant high-grade astrocytoma, particularly within the tumor core rather than the periphery (Fig. 8A). We observed the same in another dataset (Fig. S7A)13. To investigate the potential roles of Oligo_2_3_2 in GBM, we first examined the differentially expressed genes (DEGs) between Oligo_2_3_2 and the other oligodendrocyte subtypes using both single-cell and pseudo-bulk transcriptome datasets (Fig. 8B, Fig. S7B; Supplementary Data 12). Key upregulated genes in Oligo_2_3_2 were various immune-related genes including GSN, TUBB2B, HLA-A, ALDOA, CLU, TIMP1, S100A1, SERPINA3, and NGFR (collectively referred to as the Oligo_2_3_2 gene set). On the other hand, genes significantly downregulated in Oligo_2_3_2 were myelin-related genes such as OPALIN, KIF19, ALDOC, PCDH9, and DOCK9 (referred to as the Oligodendrocytes gene set). Gene ontology (GO) analysis highlighted enrichment in antigen processing and presentation and terms related to metabolic shift towards glycolysis in Oligo_2_3_2 genes, while myelination and neural projection development in Oligodendrocyte genes (Fig. 8C). We further analyzed enrichment of these two gene signatures using ssGSEA28 in bulk RNA-seq dataset from TCGA29 and CGGA30. In both datasets, we observed the trend of increased enrichment of Oligo_2_3_2 gene set from IDHmut, 1p19q co-deleted gliomas to IDHmut, 1p19q non-co-deleted gliomas to IDHwt gliomas. Similarly, we observed the opposite trend for Oligodendrocytes gene set (Fig. S7C).

To investigate the spatial organization of Oligo_2_3_2 associated marker genes at single cell resolution, we analyzed the 10x Xenium dataset. Consistent with the single-cell data, we found oligodendrocyte populations with high expression of SERPINA3, HLA-A, TUBB2B, TIMP1, CLU, HLA-DRA, and S100A1 and low OPALIN transcripts exclusively in the tumor core in Xenium data consistent with single-cell data (Fig. 8D–E). To identify cell types co-localized with Oligo_2_3_2 gene sets we analyzed both Visium and Xenium dataset. Correlation analyses of Visium dataset show that besides immune and stromal cell signatures, TC_mesnh signature is highly co-expressed with Oligo_2_3_2 genes in comparison to other subtypes of tumor cells (Fig. 8F). Using higher resolution Xenium data, we see that TC_mesnh and Mac_5_1_2 cells are found within a 20µm radius of the 90% of Oligo_2_3_2 cells (Fig. S7D). In contrast Mg_1_3 and normal Endo_cap are in close proximity of Oligo_2_1 cells (Fig. S7D). In summary, Oligo_2_1 exhibited closer spatial association with normal brain cell types, including Mg_1_3 (inactive microglia), pericytes, and Endo_cap (normal endothelial cells). In contrast, Oligo_2_3_2 was more proximally located to TC_mesnh (neoplastic mesenchymal non-hypoxic cells), Mg_1_1 (active microglia), Mac_5_1_2, and glioblastoma-associated stromal and endothelial cell types (Fig. S7D). Interestingly, Oligo_2_3_2 does not show spatial proximity to TC_mesh (neoplastic mesenchymal hypoxic cells) cells. Our Visium ST data showed a mutually exclusive expression of TC_mesh and Oligo_2_3_2 signatures (Fig. S7E). Greenwald et al.15 used a 40-marker antibody panel (CODEX) to characterize spatial organization of tumor cell states especially in the context of hypoxia at single-cell protein level. They annotated 428,395 cells from 12 samples. This analysis identified all major differentiated non-malignant cell types, including astrocytes, oligodendrocytes, neurons, vascular cells, T cells, B cells, and macrophages/microglia and key malignant cell states such as MES-like (TC_mesnh), MES-Hyp (TC_mesh), Chromatin-Reg, OPC-like (TC_OPC), NPC-like (TC_NPC), and AC-like (TC_AC). Their marker panel did not have markers to identify different oligodendrocyte subtypes but we were able to see that the GBM core samples with necrosis and high number of TC_mesh cells had almost none to few oligodendrocytes (Fig. S7F) validating the depletion of oligodendrocytes in hypoxic regions.

In order to validate Oligo_2_3_2 gene markers at protein level, we performed IHC for myelin marker Mag, actin cytoskeleton remodeling marker Gsn, and neuronal marker Ngfr on matching tissue from patients SNU33, and SNU21 used for Visium and Xenium samples. Based on Visium and Xenium spatial transcriptomic data Oligo_2_3_2 co-express these three markers. The gene and protein expression levels show high concordance for all three markers, as can be seen in the side-by-side Visium spatial gene profile and IHC protein expression where Visium equivalent tissue is highlighted in blue rectangle (Fig. S8). SNU33 showed highly heterogeneous and AF specific patterns. We analyzed three distinct regions (LE/IT, IGN, and CT) in SNU33 sample using the Oligo_2_3_2 markers, MAG, GSN, and NGFR. This was done for both Xenium and equivalent tissue areas in IHC dataset. Major cell type distributions within each region revealed a strong enrichment of neurons, and astrocytes in the LE/IT region, oligodendrocytes, macrophages, and endothelial cells in the IGN region, and neoplastic cells in the CT region (Fig. S9A). Notably, spatial overlay analysis confirmed the differential presence of MAG, GSN, and NGFR transcripts within these regions, aligning with expected cell type distributions, highest expression of all three markers in oligodendrocytes within the IGN region (Fig. S9B–C). The same pattern is observed in IHC data (Fig. S9D–E). We then performed IHC on an additional set of 11 patient samples (three grade 3 oligodendrogliomas, one grade 3 astrocytoma, two grade 4 astrocytomas and 5 glioblastomas). Four out of five glioblastoma cases showed same IGN pattern whereas none of the six IDH mutant cases exhibited such pattern (Fig. S10).

We further validated the expression of protein key markers for Oligo_2_3_2 such as S100a1, Serpina3n, B2M (MHC class I subunit), Opalin, and Sox10 as a marker for oligodendrocyte lineage cells and CRE as a marker for tumor cells, using immunofluorescence of murine high-grade glioma in situ31. Glioma-associated oligodendrocytes - identified as Sox10+ CRE- cells - were found within the tumor core and periphery, including infiltrating tumors. Quantifications of marker-positive cells revealed higher proportions of Sox10+ cells positive for S100a1 and Serpina3n in the tumor core, compared to normal mouse brain and the tumor periphery, including infiltrative regions. Sox10+ cells with Opalin expressions show a decreasing trend in the tumor core without reaching significant differences (Fig. S11).

Clinical implications of GBM associated oligodendrocytes
To explore the clinical significance of the GBM associated oligodendrocyte (Oligo_2_3_2) markers, we examined multiple publicly available datasets, including single-cell and bulk transcriptomic data with associated clinical outcomes. To assess the prognostic value of the Oligodendrocyte and Oligo_2_3_2 gene set scores in the context of established clinical and molecular parameters, a multivariable Cox proportional hazards regression analysis was performed using both the TCGA (n=603) and CGGA (n=1018) glioma cohorts. The models incorporated age, IDH mutation and 1p/19q co-deletion status, and MGMT promoter methylation status alongside the two transcriptional signatures (Fig. 9A). In both TCGA and CGGA cohorts, the Oligodendrocyte gene set score was significantly associated with favorable prognosis, whereas the Oligo_2_3_2 gene set score showed a strong association with poor overall survival. These associations remained consistent and independent of key clinical variables, including age, IDH mutation/1p19q co-deletion status, and MGMT promoter methylation, highlighting the distinct prognostic relevance of these transcriptional programs in glioma. While multivariable Cox regression strengthens the prognostic analysis, it remains challenging to completely disentangle IDH-associated transcriptional effects from subtype-specific gene programs in bulk RNA datasets, and residual confounding cannot be completely excluded.

We then evaluated the prognostic significance of the Oligodendrocytes and Oligo_2_3_2 gene signatures using TCGA and CGGA, by classifying patients into four groups based on the expression of the respective gene signatures. The analysis revealed that patients with high scores for the Oligo_2_3_2 gene set and low scores for Oligodendrocytes gene set exhibited significantly poorer survival outcomes, whereas those with high scores for the Oligodendrocytes gene set and low scores for the Oligo_2_3_2 gene set demonstrated the most favorable prognosis (Fig. 9B). Consistent trends were observed when patients were further stratified by IDH mutation status and 1p/19q co-deletion status, underscoring the prognostic relevance of these gene sets across molecularly defined gliomas subgroups (Fig. 9B). The bulk transcriptomic data (GLASS dataset)32 from matched primary (TP) and recurrent (R1) gliomas, was analyzed to determine the dynamic changes in Oligo_2_3_2 gene signatures at recurrence. ssGSEA scores demonstrated that genes upregulated in Oligo_2_3_2 exhibited increased expression in recurrent gliomas in both IDH-mutant, 1p/19q-codeleted and IDH-wildtype samples (Fig. 9C). Additionally, within IDH-wildtype tumors, recurrent samples with unmethylated MGMT promoters showed elevated expression of the Oligo_2_3_2 gene set, suggesting a potential association with treatment resistance and tumor recurrence. At single-cell level this was observed in Abdelfatteh et al. dataset where we see higher prevalence of Oligo_2_3_2 in an unmethylated recurrent tumor (Fig. S7A).

Methods
Patients samples
This study was conducted in accordance with all relevant ethical regulations. Tissue samples for this study were obtained from patients undergoing surgical resection with informed consent for their usage for research purposes in accordance with the guidelines of the Institutional Review Board of Seoul National University Hospital, which approved this study (IRB Nos. H-0507-509-153 and H-2010-123-1166). No individually identifiable information is reported in this manuscript. All clinical and demographic data are presented in anonymized and aggregated form, such that individual participants cannot be identified. Sex and gender information were collected via patient self-reports. Disaggregated data for all individual experiments are included in the supplementary data submitted with the study. Overall distribution by sex is reported in the Supplementary Data 1. Consent for reporting and sharing individual-level disaggregated data was obtained from all participants prior to study commencement. No sex-based stratified analyses were performed, as the study was not designed or powered to assess sex-specific effects and analyses focused on tumor-intrinsic transcriptional programs. No financial or material incentives were provided to patients who donated their samples. These samples were entirely consumed during the process of genetic sequencing and data acquisition; consequently, no biological material remains. All resulting genomic data, however, are available in public database or upon request.

Animals
All animal procedures were performed in accordance with institutional and federal guidelines (Stanford University APLAC protocol #33633). Mice were housed in groups of three to five and maintained on a 12-h light/dark schedule. Food and water were available ad libitum. Mice were fed a standard chow diet (Inotiv, #2018, Tekland Global 18% protein rodent diet) starting at weaning (p21). Compound transgenic mice were generated by crossing Cre-dependent BRAF CA with Ink4a/Arf flox/flox mice, as previously described (Xing et al 2025). BRAF CA/WT; Ink4a/Arf flox/flox mice were bred on a pure C57Bl/6 genetic background in our animal facility to produce sufficient mice for experiments.

No formal statistical methods were used to predetermine sample sizes. Sample sizes were determined by the availability of high-quality human tumor specimens and experimental material generated in this study, consistent with standard practice in exploratory spatial transcriptomic and histopathological analyses. Sample sizes for spatial transcriptomics and histological experiments reflect the number of biologically independent patient samples meeting pre-defined quality-control criteria. For experimental validation and animal studies, sample sizes were chosen based on prior studies in the field and feasibility considerations and were sufficient to observe reproducible effects across biologically independent samples. No data generated in this study were excluded from analysis on the basis of sample size.

Intracranial injections of adenovirus-Cre
BRAF CA/WT; Ink4a/Arf flox/flox mice were injected with 2 µl of adenovirus expressing cre recombinase under the control of the ubiquitous CMV promoter (Ad-Cre; Titer = 1×10^9 PFU/ml; Vector Biolabs, Malvern, PA) to induce BRAF V600E expression and deletion of tumor suppressor genes p16 Ink and p19 Arf. The stereotaxic coordinates selected to target the cortical layer, relative to bregma, are anterior-posterior –1.00 mm, medial-lateral –1.80 mm, and dorsal-ventral –1.5 mm. Experimental mice were monitored and humanely sacrificed in accordance with institutional animal ethics committee-approved humane endpoints, which defined the maximal allowable tumour burden as the loss of more than 15% of original body weight and/or the development of clinical signs consistent with intracranial tumor formation, such as hydrocephalus, lethargy, and hunched posture, and these limits were not exceeded in any animal. Brains of symptomatic mice were collected, and sectioned and tumor-bearing brains were subjected to immunofluorescence staining.

Spatial transcriptomics (ST)
ST experiments were performed using the 10X Visium Spatial Gene Expression kit, following manufacturer’s instructions. Here, we briefly describe the methods.

Sample preparation (FFPE patient samples)
A total of 2 FFPE sections (5um each) per sample were used to determine RNA integrity. Total RNA was extracted using RNA Isolation Kit (Qiagen, 73504) according to the manufacturer’s protocol. DV200 (the percentage of total RNA fragments >200) were determined using TapeStation (HS RNA ScreenTape, Agilent, 5067-5579) according to the manufacturer’s protocol. Only samples with DV200 ≥ 50% were qualified and used.

Deparaffinization, staining, imaging, and FFPE library construction
5µm thick sections were immediately placed on Visium array slide (Visium Spatial Gene Expression slides, 10× Genomics). Slides were deparaffinized and stained with hematoxylin and eosin. Brightfield imaging was carried out at 100x magnification with Nikon Eclipse Ti2, and post-processing was performed using NIS software. After imaging, spatially barcoded, ligated probe products were amplified using KAPA SYBR FAST qPCR Master Mix (Roche, KK4600). Fragments in the size of interest were selected using SPRIselect reagent (Beckman Coulter, B23318). Library quality check was performed using TapeStation (HS D1000 ScreenTape, Agilent, 5067-5584).

Sample preparation (Fresh frozen patient samples)
Fresh GBM tissue was collected immediately after resection and embedded in cryomolds using Tissue Tek OCT (Sakura, 4583) using a bath of isopentane and liquid nitrogen. Following freezing in OCT, blocks were stored at –80°C until processing. A total of 2 sections (10um each) per sample were used to determine RNA integrity. Total RNA was extracted using RNA Isolation Kit (Qiagen, 74104) according to the manufacturer’s protocol. RIN values were determined using TapeStation (HS RNA ScreenTape, Agilent, 5067-5579) according to the manufacturer’s protocol. Only samples with an RNA integrity value > 7 were qualified and used.

Sample preparation (Fixed frozen mouse samples)
Mice were intraperitoneally injected with ketamine (100 mg/kg) and xylazine (10 mg/kg) to induce deep anesthesia and subsequently transcardiac perfusion was performed with phosphate-buffered saline (PBS) followed by 4% paraformaldehyde (PFA)/PBS. Brains perfused with 4% PFA were then resected, post-fixed in 4% PFA/PBS for 2 h on ice, cryoprotected in 15% sucrose/PBS for 24 h and then 30% sucrose/PBS for at least 48 h, followed by embedding in Tissue-Tek® O.C.T. Compound (Sakura Finetek USA, Inc., Torrance, CA). The embedded tissues were stored at –80°C until sectioned. 10µm-thick coronal sections of the brain were cut on a cryostat (Leica, Wetzlar, Germany), collected onto Superfrost Plus slides (Menzel Glaser, VWR) and air dried for 1 h before storing at –80°C until stained.

Fixation, staining, imaging, and construction of fresh frozen cDNA libraries
10µm thick sections were immediately placed on Visium array slide (Visium Spatial Gene Expression slides, 10× Genomics). Slides were fixed in 100% methanol, stained with hematoxylin and eosin. Brightfield imaging was carried out at 100x magnification with Nikon Eclipse Ti2, and post-processing was performed using NIS software. After imaging, sections were permeabilized to release and capture mRNA to oligonucleotides on the capture areas. Following permeabilization, the on-slide reverse transcription (RT) reaction and second strand synthesis was performed. Spatially barcoded cDNA was amplified using KAPA SYBR FAST qPCR Master Mix (Roche, KK4600). Fragments in the size of interest were selected using SPRIselect reagent (Beckman Coulter, B23318). Library quality check was performed using TapeStation (HS D1000 ScreenTape, Agilent, 5067-5584). To optimize cDNA fragments for Illumina NextSeq, fragmentation and double-sided size selection was performed using SPRIselect reagent.

Immunofluorescence staining on frozen mouse brain sections
Frozen 10µm mouse brain sections on Fisherbrand™ Superfrost™ Plus Microscope Slides (Thermo Fisher Scientific) slides were washed briefly in PBS and blocked with PBS containing 0.3% Triton X-100 and 10% normal donkey serum for 1 hr at room temperature. The sections were then incubated with primary antibodies at room temperature overnight in a humidified chamber, washed briefly in PBS three times next day, followed by 1 h incubation with secondary antibodies, the slides were then washed three times with PBS, and coverslipped with a DAPI-containing mounting solution (VectaShield H20002; Vector Laboratories, Newark, CA) and subjected to confocal microscopic analysis. Some primary antibodies were incubated simultaneously. The following primary antibodies were used: mouse anti-Cre recombinase (1:500; MAB3120, Millipore Sigma, Burlington, MA), rabbit anti-Serpina3n (1:200; MA5-41236, Thermo Fisher Scientific), rabbit anti-Opalin (1:300; ab121425, Abcam), rabbit anti-S100a1 (1:300; PA1-932, Thermo Fisher Scientific), rabbit anti-B2M (1:500; PI701250, Thermo Fisher Scientific) and goat anti-Sox10 (1:100; AF2864, R&D Systems). Secondary antibodies raised in donkey and conjugated to Alexa Fluor 488, Alexa Fluor 594 or Alexa Fluor 647 were purchased from Jackson ImmunoResearch Labs (West Grove, PA) or Invitrogen (now Thermo Fisher Scientific) and used at 1:200 dilution.

Confocal imaging, image processing and statistical analysis
Stained 10µm thick coronal sections were imaged by laser scanning confocal microscopy (Zeiss LSM700) which was used to detect up to four fluorophores by laser excitation at 405, 488, 561, and 639 wavelengths. Confocal images were acquired using a 20X objective. Images were imported into Fiji/ImageJ software for quantification of cellular density in the regions of interest. All analyses were performed in a blinded fashion. Statistical analysis of quantitative results was performed using GraphPad Prism software. Statistical significance was determined using two-way ANOVA with Sidak’s multiple-comparison test. Quantitative data are reported as mean ± SEM, and any values of P < 0.05 were considered statistically significant.

Immunohistochemical staining on tumor tissues from human patients
Formalin-fixed paraffin-embedded tumor tissues were sectioned at 4μm thickness for IHC. IHC staining was performed using the BenchMark ULTRA system (Ventana-Roche, Manheim, Germany) with antibodies for GSN (PA5-29910, Invitrogen), MAG (D4G3, Cell Signaling), and NGFR (MA5-31968, Thermo Fisher Scientific).

ST sequencing
The average length of the final libraries was quantified using TapeStation (HS D1000 ScreenTape, Agilent, 5067-5584). The paired-end sequencing was performed on the Illumina NextSeq 550 using the following read protocol; 28 cycles for read 1; 10 cycles for i7 index read: 10 cycles for i5 index read; 90 cycles for read 2.

Single-cell RNA-sequencing
Tissue dissociation. Tissue specimens were collected and immersed in RPMI media at room temperature immediately after surgical resection. The MACS brain tumor dissociation kit (Miltenyi Biotec, Auburn, CA, USA) and gentle MACS™ Dissociator (Miltenyi Biotec) were used to dissociate the tissue samples within 2 h after collection, and debris was removed using MACS Debris Removal Solution (Miltenyi Biotec) according to the manufacturer’s protocol. Only the samples with cell viability > 80% after dissociation were used in this study. Isolated single cells were cryopreserved until further use.

Construction of cDNA libraries and sequencing. The library preparation was performed using Chromium Next Gem Single-cell 5’ kit v2 (10x Genomics, USA, 1000263) with a cell recovery target of 10,000, following the manufacturer’s instruction. The libraries were processed according to manufacturer’s recommendation (10x Genomics, USA). The libraries were sequenced on a HiSeq 2500 system using 100-bp paired-end sequencing, and the raw sequencing data was processed using Cell Ranger v5.0.1 with default setting. The reads were aligned to the human reference genome (hg38).

CITE-seq. Approximately 1 million cells were incubated with Fc receptor blocking solution (Human Trustain FcX™, BioLegend, 422301) at 4’C for 10 min and were stained with TotalSeq™-C Universal Cocktail, V1.0 (BioLegend, 399905) at 4’C for 30 min. After a centrifuge (600g, 4’C, 5 min), cells were resuspended with 3 ml of Cell Staining Buffer (BioLegend, 420201). Cells were washed twice more and resuspended in 500 ul of Cell Staining Buffer. Antibody-stained cells were processed by Chromium controller (10x Genomics, USA) with a cell recovery target of 5000 for the generation of gene expression libraries and Feature Barcode libraries. Gene expression libraries and Feature Barcode libraries were prepared using Chromium Next Gem Single-cell 5’ kit v2 (10x Genomics, USA, 1000263) and 5’ Feature Barcode kit (10x Genomics, USA, 1000256) following the manufacturer’s instruction. The libraries were sequenced on HiSeq 2500 system using 100-bp paired-end sequencing yielding 25K reads per cell for gene expression libraries and 5K reads per cell Feature Barcode libraries, respectively. The raw sequencing data was processed using Cell Ranger v5.0.1 with default setting. The reads were aligned to the human reference genome (hg38).

Quantification and statistical analysis
Preprocessing of 10x files. The 10x Visium data was aligned using SpaceRanger, and the Cell-seq/CITE-seq data was aligned using the CellRanger software. The following are the commands used,
CellRanger command for CITE-seq:
cellranger count –id=[sampleid] citeseq --libraries=library.csv --transcriptome=refdata-gex-GRCh38-2020-A --feature-ref=TotalSeqC_universal_GBM_feature_ref.csv
Cellranger command for Cell-Seq:
cellranger count --id=[sampleid] --fastqs=[fastq file location] --sample=[sample name] --transcriptome = refdata-gex-GRCh38-2020-A
Spaceranger command for Visium:
spaceranger count --id=[sampleid] --transcriptome=refdata-gex-GRCh38-2020-A/ --probe-set = Visium_Human_Transcriptome_Probe_Set_v1.0_GRCh38-2020-A.csv --fastqs=[fastq file location] --sample=[sample name] --image=image.tif --slide=[slide name] --area=A1 --localcores=32 --localmem=64

Data normalization of single-cell and spatial transcriptome datasets. Raw count matrices (raw_feature_bc_matrix) were used for downstream analysis. The data was normalized using the counts per million (CPM) method, followed by log-transformation. To prevent undefined values during log transformation, a pseudocount of 1 was added. For single-cell RNA-seq data, cells with a total RNA count (nCount_RNA) below 1000 were excluded. Additionally, to enrich for granulocytes, cells were retained if the mean expression of IFITM2, S100A12, and S100A8 exceeded a threshold of 2. For 10x Genomics Visium spatial transcriptomics data, spots with nCount_RNA less than 2000 were filtered out.

Data normalization of CITE-seq. mRNA expression data were normalized using the same approach applied to the single-cell RNA-seq data, involving CPM normalization followed by log-transformation with a pseudocount. Protein expression data were further normalized to enable comparison across samples, as protein counts exhibited greater inter-sample variability compared to mRNA expression. To address this, log-transformed CPM protein expression values were z-scored within each sample using the mean and standard deviation derived specifically from oligodendrocyte cells. Oligodendrocytes were selected as internal controls for normalization due to their consistently low levels of cell surface marker protein expression.

Clustering of single-cell data. A modified version of Recursive Consensus Clustering (RCC) was employed for cell clustering and marker gene identification16. Briefly, RCC identifies primary clusters at the first level and then recursively applies clustering within each cluster to reveal hierarchical substructures across multiple levels. In our implementation, the major modification involved replacing the original ConsensusClusterPlus algorithm34 with Seurat v4.2.035 for improved scalability and integration with single-cell analysis workflows. Before clustering, batch correction was performed using Harmony v0.1.036 to integrate multiple samples. Additionally, UMAP coordinates were recalculated after recursive clustering to enhance the visualization of hierarchical cluster relationships.

Cell-type assignment in single-cell. Level 1 RCC clusters were annotated as major cell types based on the expression of canonical marker genes. These major cell types were then subjected to recursive clustering up to three levels to resolve finer cellular subtypes. Marker genes were identified for each sub-cluster. Cluster names reflect the hierarchical structure, with the major cell type followed by cluster numbers separated by dashes for each level, or by abbreviations corresponding to previously characterized cellular states. Cell type assignments were validated using CITE-seq protein expression data. Neoplastic cells were further confirmed based on their inferCNV-derived copy number variation (CNV) scores.

CNV calculation. Copy number variation (CNV) analysis of spatial transcriptomics (ST) data was performed using the SPATA2 package in R13. CNV profiles were inferred by smoothing gene expression signals across spatially adjacent spots and comparing relative expression along chromosomal coordinates to reference non-malignant cell populations. Spatial spots exhibiting large-scale, concordant chromosomal gains and losses were classified as malignant, whereas spots lacking such patterns were classified as non-malignant. For single-cell RNA-seq data, CNV inference was conducted using inferCNV5. The single-cell data was split by major cell types for CNV analysis. Cells labeled as ‘Oligo_2_1’, ‘Mg_1_2’, ‘Mg_1_3’, and ‘Mg_1_1’ from normal brain samples were used as reference controls for CNV inference.

Reference single-cell dataset. The reference single-cell dataset was constructed using the datasets listed in Supplementary Data 3:
1. Allen Institute: Single-nucleus transcriptomes from 76,533 nuclei obtained from two post-mortem human brain specimens were used to survey cell type diversity in the primary motor cortex (M1 Cor M1)20. Cell type annotations from the original publication were retained.
2. Xie et al.: Single cell RNA-seq data from tumor cores and matched peripheral tissue of four glioblastoma patients22 were used, particularly for endothelial cell identification. RCC was applied to identify clusters, and cell types were annotated using scMAGIC37.
3. Winkler et al.: Single cell data from dissociated vascular cells of the adult human brain and arteriovenous malformations21 were incorporated. The original cell type annotations were used.
4. Abdelfattah et al.: Single-cell RNA-seq data from 44 glioma samples (newly diagnosed and recurrent GBM and LGG)23 were included. Clustering was conducted using the RCC algorithm, and cell types were subsequently annotated using scMAGIC based on canonical markers.
5. Pombo Antunes et al.: CITE-seq data from myeloid cells in newly diagnosed and recurrent GBM samples were incorporated25, using the cell type annotations provided in the original study.
6. Matthewson et al.: Single-cell RNA-seq data of glioma-infiltrating T cells isolated from fresh tumor samples in 31 adult patients with either IDH-wildtype GBM or IDH-mutant glioma24 were included. Cell type annotations from the original publication were retained.
To capture further biologically significant subtypes identified in our cohort, we also incorporated selected cells from our single-cell dataset into the referenced dataset.

Cell state assignment for reference single-cell data. Cell states were assigned using four complementary approaches:
1. Original publication annotations: Cell state annotations from the original publications were adopted without modification. These include a range of cell types and subtypes such as Artery subtypes (Art1, Art2, Art3), Astrocyte subtypes based on single-gene enrichment (AstroAQP1, AstroPLCG1, AstroSERPINI2), DC, migrating DC, FBMC, Fibroblast, non-classical Monocyte, Excitatory Neurons, Inhibitory Neurons, OPC, Pericyte, SMC, proliferating T cells, Venous, and Venule.
2. Recursive Consensus Clustering (RCC) combined with canonical markers: Marker genes were identified across multiple clustering levels using RCC. These markers were then reviewed and matched with well-established cell type signatures to assign cell identities. This approach allowed classification of immune and stromal subsets including B cell, B cell (plasma), CD4 T cell, CD8 T cell, CD8 tissue-resident memory T cell (TRM), capillary endothelial cell (Endo_cap), glioblastoma-perivascular cells (gbmPeri), granulocyte, proliferating macrophage (Mac_prolif), proliferating microglia (Mg_prolif), monocyte, natural killer T cell (NKT), red blood cell (RBC), proliferating T cell (TC_prolif), and regulatory T cell (Treg).
3. Literature-based marker analysis for tumor cell states: To assign states to tumor cells from the SNUH single-cell dataset, we integrated clustering with gene signature analysis based on studies by Neftel et al. (2019) and Garofano et al. (2021). Identified tumor cell subtypes included TC_AC, TC_mesh, TC_mesnh, TC_MTC, TC_NPC, TC_oligo, and TC_OPC. Other signatures, such as glycolytic/plurimetabolic (GPM), neuronal (NEU), and proliferative/progenitor (PRR), were excluded due to substantial overlap with the Neftel classification. For each neoplastic cell, we computed:
• Mmax = max(TC_log2CPM average)
• If Mmax < 0, the cell was labeled as misc; otherwise, it was assigned to the corresponding TC subtype with the highest Mmax. The gene signatures used are listed in Supplementary Data 4.
4. RCC-based de novo stromal subtype identification: RCC enabled the identification of robust subclusters within glioblastoma stromal cells, microglia, macrophages, and oligodendrocytes. These subclusters were labeled using a hierarchical cluster naming scheme (e.g., gbmEndo_9_1, Mac_5_1_3, Oligo_2_3_1), reflecting their recursive clustering structure. As the functional roles of these subpopulations remain to be fully elucidated, we retained cluster-based labels to preserve analytical granularity for future studies.

Marker calculation. Marker gene identification was performed using the scran package38. Genes were considered highly specific markers if they met a false discovery rate (FDR) threshold of < 0.01 and a minimum log2 fold change > 0. Differential expression was calculated as the log2 fold change of each gene in the target cluster (annotated feature/cell type) relative to all other clusters. The minimum log2 fold change represents the smallest log2 fold change across all pairwise comparisons, ensuring that selected markers consistently exhibit higher average expression within the target cluster.

Differential gene expression calculation for single-cell data. Differentially expressed genes were identified using the clusterDE39 R package, which corrects for the double-dipping bias that occurs when clustering and differential expression testing are performed on the same dataset. clusterDE was applied to the normalized Seurat object to identify genes whose expression patterns robustly defined each cluster while controlling the false discovery rate (FDR). Genes with FDR < 0.05 were considered significant and retained as differentially expressed genes. Log2 fold change values were calculated by subtracting the average expression of control cells from that of test cells. Alternatively, differential gene expression analysis was also performed using pseudo-bulk data derived from the single-cell transcriptomic dataset. DESeq240 was used to normalize and identify differentially expressed genes. For the pseudo-bulk differential expression analysis, DESeq2 was run with a design formula that included both patient identity and cell type as covariates to account for inter-patient variability and differences in cellular composition. Genes with adjusted p < 0.01 and |log₂FC| > 0.8 were considered significant, with log₂FC > 0.8 classified as up-regulated and log₂FC < –0.8 as down-regulated. Volcano plots were generated using the EnhancedVolcano package in R41, which enables the creation of publication-ready visualizations with enhanced coloring and annotation features.

Cell type proportions calculation for core vs periphery tissue. To evaluate whether specific cell states were enriched in the tumor core compared to the periphery, we calculated cell type frequencies for each sample in both regions. Statistical differences between core and peripheral tissues were assessed using a paired t-test on the cell proportions across all cellular states.

Integration of single-cell and spatial transcriptome data. Cell state scores for each spatial transcriptomics spot were calculated using the marker gene sets listed in Supplementary Data 6. For a given cell state C, the score was defined as the average of the log2 CPM values of all marker genes associated with C. To enable binary classification of cell states, we visualized histograms of cell state scores across all spots. These distributions were typically multimodal, allowing manual selection of score thresholds for each cell state. Supplementary Data 14 includes the histograms with corresponding cut-off values, spatial plots of cell state scores, and binary cell state assignments.

Anatomic feature (AF) assignment for spatial transcriptome data. Based on binary cell type assignments and expression profiles, spot identities were refined using a hierarchical rule-based approach:
1. LE_GM (Leading Edge - Gray Matter): Spots positive for either Neurons_Ex or Neurons_In were assigned as LE_GM.
2. MVP (Microvascular Proliferation): Spots co-expressing both gbmEndo_9_1 and gbmPeri markers were assigned as MVP.
3. HBV (Hyperplastic Blood Vessel): Spots expressing either gbmEndo_9_1 or gbmPeri, and not already assigned MVP, were assigned as HBV.
4. PAN (Pseudopalisading cells Around Necrosis): Spots expressing TC_mesh but lacking TC_mesnh expression, were assigned as PAN.
5. PNZ (PeriNecrotic Zone):
Spots previously labeled as PAN and additionally positive for Mac_5_1_2 were reassigned as PNZ.
Spots co-expressing TC_mesh and TC_mesnh, while still unassigned, were also assigned to PNZ.
6. IGN (Immune-Glial Niche):
Spots positive for TC_mesnh and negative for TC_mesh, while still unassigned, were annotated as IGN.
Additionally, spots positive for either Mac_5_1_2 or gbmFib_6, and originally labeled LE_GM or were yet unassigned, were also classified as IGN.
7. CT (Cellular Tumor): Spots with the presence of any of the non-mesenchymal neoplastic cell populations (TC_oligo, TC_prolif, TC_OPC, TC_NPC, or TC_AC), and yet unassigned, were designated as CT.
8. IT (Infiltrating Tumor): Spots expressing Oligo_2_1, Oligo_2_2, or Oligo_2_3_1, and previously classified as CT, were relabeled as IT.
9. LE_WM (Leading Edge - White Matter): Spots positive for Oligo_2_1, Oligo_2_2, Oligo_2_3_1, or OPC, and still unassigned, were assigned as LE_WM.
10. BV (Blood Vessels): Spots expressing Endo_cap or SMC and previously labeled LE_GM or yet unassigned were designated as BV.
11. PNZ (Refined): Among spots labeled as PAN, those with average gene set expression of TC_mesh < 6 were reclassified as PNZ.
This rule-based spatial reannotation enabled consistent, biologically informed classification of microenvironmental regions across the tumor tissue.

Xenium panel development. Gene selection for the custom glioma Xenium panel was performed using a reference single-cell dataset. Marker genes significantly enriched in both major and minor cell types were identified using the findMarkers function from the scran package. A total of 2689 cell type-enriched genes were selected for initial analysis. To eliminate functionally redundant genes, recursive feature elimination (RFE) was applied using the caret package in R (Kuhn, 2008)42. This process yielded 383 genes as significant discriminative markers. To validate their cell type specificity and utility, cell type assignment was independently reassessed using only these 383 genes. Based on assignment accuracy and inter-gene correlation analyses, a final set of 348 genes was selected for inclusion in the glioma Xenium custom panel.

Xenium data analysis. Xenium data was generated for four patients, namely SNU18, SNU21, SNU25, and SNU33.

Data normalization. The cell_feature_matrix.h5 files generated by 10x Genomics Xenium were used for downstream analysis. Gene expression data were normalized using the square root transformation method implemented in the stlearn module43, which has been shown to provide robust normalization for spatial transcriptomics data.

Cell type assignment
Major cell type. Based on the single-cell reference dataset and previously identified marker genes, a subset of genes with exclusive enrichment in major cell types was selected. From the Xenium transcript file provided by 10x Genomics, only cells with nucleus_overlap were retained for downstream analysis to ensure high-confidence cell assignments. For each cell, the sum of expression values for marker genes corresponding to each major cell type was computed. An initial cell type was assigned based on the gene set with the highest cumulative expression. To support this, a Xenium reference matrix was generated, consisting of 348 selected genes (rows) and their average expression across individual major cell types (columns). Cells that remained unassigned in the initial round were subsequently annotated by calculating their Pearson correlation with the reference matrix and assigned to the cell type exhibiting the highest correlation coefficient.

Minor cell type. Major cell types were further resolved into minor subtypes by subsetting and reclustering each major population using Seurat. Dimensionality reduction was performed using PCA with 100 principal components. Clustering was carried out using the Louvain algorithm, with the resolution parameter adjusted between 0.2 and 1, depending on the heterogeneity observed within each major cell type. Minor cell type identities were assigned based on marker gene expression, using reference markers derived from the reference single-cell dataset (see Supplementary Data 5).

Anatomical feature assignment. To investigate anatomic features (AFs) within the high-resolution Xenium spatial transcriptomics data, we defined spatial grids (2376 μm²; area-matched to Visium spots) across tissue sections from four GBM patients (SNU21, SNU33, SNU18, and SNU25).

Trajectory selection. We selected 11 representative spatial trajectories across the four samples to capture key biological transitions and microenvironmental heterogeneity observed in our dataset. Trajectories were chosen based on anatomical landmarks, visible transitions in gene expression, and areas of interest identified from H&E morphology and Visium correspondence. Each trajectory was designed to traverse multiple relevant AFs, such as perivascular regions, tumor parenchyma, infiltrative margins, and fibrotic or immune-infiltrated zones. Each spatial trajectory comprised a variable number of spatial grids, reflecting differences in tissue architecture and trajectory length (range: 85–235 grids per trajectory; exact counts provided in the Source Data files).

40-gene selection for AF assignment. A curated panel of 40 genes was used to assign AF identities to each grid. Selection was based on established markers of canonical GBM anatomical structures—incorporating markers that distinguish normal versus GBM-associated vasculature26, as well as our in-house oligodendrocyte markers and established microglia/macrophage/monocyte25.

Manual AF assignment. Initially, AF labels were assigned to grids along selected trajectories using a combination of spatial gene expression profiles (from the 40 selected genes), focused cell type composition, and direct comparison to H&E images. Special attention was given to subregions such as the immune fibroblast niche (IFN), which was defined by the presence of fibroblast and immune markers and the absence of glial signatures within the Immune-glial Niche (IGN).

Distance-based AF label propagation. For the remaining grids outside the manually annotated regions, AF labels were assigned based on the shortest Euclidean distance to an AF-annotated grid within the same trajectory. This approach allowed continuous annotation of AFs across tissue while minimizing manual bias and ensuring spatial coherence.

Gene ontology analysis
Functional enrichment was assessed using ToppFun44, based on the top 100 markers from the upregulated and downregulated gene sets.

Survival analysis
Survival analysis was performed using the survival package in R45. Patient stratification was based on Oligodendrocytes and Oligo_2_3_2 gene set scores computed by ssGSEA. Stratification thresholds were selected such that samples with the highest Oligodendrocytes scores were enriched for IDH-mutant tumors, whereas those with the highest Oligo_2_3_2 scores predominantly corresponded to IDH-wildtype cases. This approach was motivated by consistent observations across multiple datasets indicating that the Oligo_2_3_2 signature is largely, though not exclusively, associated with IDH-wildtype gliomas. Kaplan–Meier survival curves were generated using the survminer46 package to visualize survival distributions across groups. In addition, multivariable Cox proportional hazards regression models were constructed to evaluate the independent prognostic value of molecular and clinical attributes. Hazard ratios (HR) with 95% confidence intervals were estimated, and forest plots were generated using the ggforest function to visualize the relative effect of each variable on survival. The Oligodendrocytes and Oligo_2_3_2 gene signatures identified in this study were incorporated as continuous covariates in the Cox models alongside patient age, IDH mutation status, 1p/19q co-deletion status and MGMT status to assess their prognostic relevance.

Single sample gene set enrichment analysis (SSGSEA)
To assess specific cell type enrichment patterns in both single-cell and spatial transcriptomics datasets, the Single Sample Gene Set Enrichment Analysis (ssGSEA) method was applied. Gene sets were derived from the marker identification approaches described above. The resulting ssGSEA scores were subsequently z-scored and used for downstream analyses.

Cell-type enrichment in anatomical features
We first calculated the average log2 CPM values of each cell type gene set for every spot. Using the density function in R, we computed the distribution of each gene set score across all anatomical features (AFs). To assess the significance of these distributions, AF assignments were randomized 1000 times, and density curves were recalculated for each iteration. For both the original and shuffled datasets, the area under the curve (AUC) was computed for gene set scores greater than 2. The enrichment score was defined as the ratio of the AUC from the original data to the mean AUC obtained across the 1000 shuffled iterations, thereby minimizing bias from random sampling. For each AF-cell type pair, empirical p-values were computed as the AUC of shuffled statistics greater than or equal to the observed AUC, following the standard permutation formula:
p = (count of permuted values ≥ observed) + 1 / Nperm + 1
False discovery rate (FDR) correction across all comparisons was then applied using the Benjamini–Hochberg method.

Neighborhood analysis
Neighborhood analysis was performed to characterize spatial associations between anatomical features (AFs) in both the Visium and Xenium transcriptomic datasets. For each dataset, spatial proximity was defined based on immediate neighbors - spots in the Visium dataset and grids (each grid area equivalent to a Visium spot) in the Xenium dataset. The immediate neighbors for each spot/grid were first identified using their spatial coordinates, and neighboring units were then aggregated according to their assigned anatomical features. This process yielded an AF × AF adjacency matrix, where each entry represented the number of neighboring interactions between two anatomical features. To correct for differences in feature abundance, the raw interaction counts were normalized by the total number of neighboring pairs, resulting in a normalized neighborhood frequency matrix.

To assess the statistical significance of spatial associations, a permutation-based null model was generated. The anatomical feature assignments were randomly shuffled 1000 times, and the neighborhood analysis was recomputed for each iteration to establish a background distribution of random neighbor counts. For each pairwise interaction, the ratio of the observed normalized neighborhood count to the mean randomized count was computed, providing an enrichment score that quantified the strength of spatial association. Ratios greater than 1 indicated enrichment (features co-localizing more frequently than expected by chance), whereas ratios below 1 indicated spatial exclusion. For each spatial pair, empirical p-values were computed as the area under the curve (AUC) of shuffled statistics greater than or equal to the observed AUC, following the standard permutation formula:
p = (count of permuted values ≥ observed) + 1 / Nperm + 1
where Nperm = 1000 corresponds to the number of randomizations performed. False discovery rate (FDR) correction across all pairwise comparisons was applied using the Benjamini–Hochberg method to control for multiple testing. These resulting ratios were used as edge weights, while the anatomical features served as nodes, to construct an interaction network, which was visualized using Cytoscape47.

DETR model for nuclei detection
We employed DINO-DETR48 a variant of the DETR49 model, for object detection tasks. The training dataset consisted of 27 regions of interest (ROIs) derived from 12 whole slide images (WSIs) from the 10x Genomics Visium dataset. Within these ROIs, nuclei were manually annotated, including areas with substantial background noise to ensure robustness. For model training, each ROI was subdivided into smaller image patches of size 256 × 256 pixels. In total, 352 patches were used for training and 89 for validation.

Mutual exclusivity test
The Rediscover package50 in R used to identify mutually exclusive patterns in ST data.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Park et al. analyzed different anatomic regions of GBM through ST and paired scRNA-seq, adding more data to the fast-growing GBM ST datasets. However, the results are mostly descriptive offering limited biological insights, lack comparison with existing datasets, and many of the claims are not substantiated by data. Also, the writing and figure presentation needs to be substantially revised to improve clarity.\n1. The introduction is poorly written. The first paragraph has no reference. The logic flow from the second to third paragraph is mixed up. Many recent publications using spatial transcriptomics to investigate GBM have not been properly discussed and compared with this study to highlight the knowledge gap.\n2. The sampling and sequencing strategy appears confusing. The authors should cite and describe their Table S1 in the first paragraph of Results, and summarize the mutational profiles of GBM samples analyzed, especially the IDH mutation status.\n3. Fig.2, it is unclear how the authors distinguish neoplastic vs TME cells.\n4.“We were able to detect a granulocyte population from the freshly processed samples, a cell type missing, from most human GBM scRNA-seq public datasets” Which granulocyte population? Neutrophils? The authors need to back up their claim with specific markers and compare with other scRNA-seq datasets.\n5. What is the rationale for dividing Oligo into 4 subtypes? Have these subtypes been previously shown? Is there any disease relevance? Also, I don’t understand why the authors put the comparison in Fig. S7A. They should immediately show it in the main figure.\n6. Fig.3D, I am not sure what the authors are showing here.\n7. The web link provided does not appear to allow users to explore the spatial gene expression.\n8. Fig.6A, OPALIN and MOBP are not classical markers for myelinating OCs.\n9. What is the cellular frequency of Oligo_2_3_2 in the tumor core and periphery? How does it impact on the tumor cells? I am not fully convinced this is an important cell subpopulation.", "reviewer_2": "Park et al. present a comprehensive single-cell and spatial transcriptomic dataset of glioblastoma samples, complementing previous studies by Ravi (2022) and Greenwald (2024). While this dataset, which employs CITE-seq, sc/snRNA-seq, Visium, and Xenium technologies, is a valuable resource, it offers limited biological insights.\nThe authors identified a potential novel cell type, Oligo_2_3_2, interacting with another new cell type, gbmEndoPeri, but further characterization is needed to understand their phenotypes and clinical implications. The proposed link between neuroinflammation and the transformation of oligodendrocytes into pathogenic states associated with impaired remyelination and neurodegenerative processes remains unsupported by current analysis.\nMajor Points:\nThe authors provide a rich dataset of spatial transcriptomics. However, there have been previously published spatial transcriptomic datasets:\n- Ravi, V. M. et al. Spatially resolved multi-omics deciphers bidirectional tumor-host interdependence in glioblastoma.\nCancer Cell 40, 639-655.e13 (2022).\n- Greenwald, A. C. et al. Integrative spatial analysis reveals a multi-layered organization of glioblastoma. Cell (2024).\nThe analysis of spatial organization leads to similar conclusions made in the studies above. The conclusion and the knowledge contribution to the field did not come across strongly. Can the authors clarify specifically how this dataset can enable novel discovery for the community compared to the other datasets?\nTo increase value to the community and further validate findings, increasing the Xenium dataset would be beneficial.\nMore attention is required in the figures and text. Some legends and labels are unclear (see section below).\nThe AUC analysis for spatial enrichment of cell types in AFs is missing in the manuscript. More evidence and analysis are required to support the claims about spatial co-localization.\nThe pseudotime analysis in Figure 6C is not supported. These are non-malignant cells that are terminally differentiated and thus, the pseudotime framework to model trajectories should not apply.\nThe main finding of this manuscript is the novel oligodendrocyte subtype. It would be good to validate the existence of this subtype with protein marker staining. Although the authors have identified Oligo_2_3_2, it is not clear from the analysis how the novel subtype contributes to the pathology of GBM. Is there an association between the proportion of Oligo_2_3_2 and overall survival?\nMinor Points:\nFigure 1B: It is unclear where the cell count per cell type is coming from. Why is the microglia count much higher than tumour cells?\nFigure 2A/B: Why does the UMAP look different between different levels of labels?\nFigure 2C: Need to briefly describe what the states stand for and how they align with previously reported cell states.\nFigure 2C-D: Add a statistical test to support the claim that core samples have a higher proportion of tumor cells compared to peri samples.\nFigure 3D: This plot is unclear. Please make the message clearer, either through text or an alternative visualization.\nFigure 5B: Can be moved to supplementary.\nFigure 6A: “We scrutinized the differentially expressed genes (DEGs) between oligodendrocytes from the tumor core and periphery”: Oligo_2_1 is not enriched towards the periphery. What is the rationale for specifically comparing Oligo_2_3_2 to Oligo_2_1?\nFigure 6 Caption: Ensure captions are consistent with the figure.\nFigure S1: Legend missing for dot size and colour.\nFigure S3: Legend missing for Visium spot annotations.\nFigure S6B: Needs a legend. If the enrichment of Oligo_2_3_2 towards the core is significant.\nFigure S7A: The Y-axis is missing or too small.\nFigure S7B: Change the colour of the colour bar to be more easily interpretable. Report FDR instead of p-value.\nIn the text, it should be Neftel and not Verhaak.\n\"TC_mesnh, resembling Mes1-like cells in Verhaak classification, demonstrate an association with immune infiltration, whereas TC_mesh (Mes2-like cells in Verhaak classification) are linked with chronic hypoxia.\"\n“Spots from patients with astrocytoma (Ast), IDH mutant (IDHm) (SNU19, SNU22) were classified into LE_GM, LE_WM, IT, and CT”: Does it make sense to label IDHmut cases with GBM anatomical features?\n“This suggests that despite sharing mesenchymal gene signatures, these two tumor cell states may have arisen under different environmental pressures.”: How does this connect with Sten Linnarsson’s GW-like and FB-like states?", "reviewer_3": "In this manuscript (NCOMMS-24-33567-T, “A spatially resolved human glioblastoma atlas reveals distinct cellular and molecular patterns of anatomical niches”) Park*, Sonpatki* et al., perform spatial transcriptomics with matched single cell RNA sequencing and cellular indexing of transcriptomes and epitopes (CITE-Seq) on FFPE embedded glioma samples. They present a spatial map of 56 cellular compartment based on 115,914 spatial transcriptomes from 32 tissues sections that could be of use to the research community. The authors suggest that their study provides novel insight into the cellular architecture of GBM compared to similar studies that have been published in Nature Communications and other journals (PMCID: PMC9974099, PMCID: PMC9950149, PMID: 35700707, PMCID: PMC10336135, PMCID: PMC9767870).\nAlthough the authors show some unique data, including some previously unreported diversity in oligodendrocytes, there are some issues ranging from major to minor that reduce my enthusiasm for the manuscript in its present form.\n1. In general, in its current form I fail to understand the overall questions or hypothesis that is trying to be answered through this analysis. I understand the authors are claiming they have a unique dataset, which I do not want to diminish as it could be useful to the research community. However, multiple recent papers, some listed above, have extensively conducted spatial transcriptomics with matched RNA sequencing.\n2. The authors highlight their unique oligodendrocyte cell type that they have identified throughout their analysis, but that is not discussed in the abstract or the title.\n3. It is current form, it is not clear how the Oligo2_3_2 population marked by expression of HLA-A is novel and/or clinically relevant. What is the function of this population? Many manuscripts have extensively categorized and identified the OPC lineage. If the authors want to focus on this, I suggest making the paper about this and providing more evidence as to how this population is unique and why it is important (linked to any patient outcome, genetic mutation, or any other characteristic (treatment response, biological sex, etc.)).\n4. The authors suggest that this shared resource will help develop new and effective combinatorial therapies but do not give an example as to how. What makes this analysis and dataset different than what is already out there? How can this dataset provide information that no other dataset can when looking for combination therapies. Not to belabor the point, this gets back to the authors not answering a clearly delineated question or hypothesis.\n5. The patient classification should be based on sex (not gender) as this is the biological variable.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nPark et al. analyzed different anatomic regions of GBM through ST and paired scRNA-seq,\nadding more data to the fast-growing GBM ST datasets. However, the results are mostly\ndescriptive, offering limited biological insights, lack comparison with existing datasets, and\nmany of the claims are not substantiated by data. Also, the writing and figure presentation needs\nto be substantially revised to improve clarity.\n- We thank the reviewer for the constructive feedback. As suggested, we have revised the\nmanuscript accordingly and addressed each of the points raised. Specifically, we have expanded\nthe biological insights through additional analyses, incorporated systematic comparisons with\nexisting GBM ST datasets, carefully revised or removed claims to ensure they are fully\nsupported by data, and substantially improved the writing and figure presentation to enhance\nclarity and readability. A detailed point-by-point response is provided below.\nComment #1.1. The introduction is poorly written. The first paragraph has no reference. The\nlogic flow from the second to third paragraph is mixed up. Many recent publications using\nspatial transcriptomics to investigate GBM have not been properly discussed and compared with\nthis study to highlight the knowledge gap.\nResponse #1.1. We thank the reviewer for this valuable comment. In the revised manuscript, we\nadded appropriate references to the opening paragraph and have thoroughly rewritten the\nIntroduction to improve clarity and logical flow. In the updated manuscript, we included detailed\ncomparisons with recent spatial transcriptomics studies, mapping our AFs to modules and layers\ndefined by Ravi et al. and Greenwald et al. and clarified how our study builds upon and differs\nfrom them (Figure 4C, Supplementary File 1, Table S6). Specifically, our work uniquely\nintegrates both fresh frozen and FFPE samples with matched single-cell multi-omics and\nsubcellular-resolution Xenium profiling, enabling us to identify novel vascular and glial niches,\nincluding the clinically relevant Oligo_2_3_2 subtype. These revisions highlight the knowledge\ngaps left by prior studies and frame the distinct contributions of our atlas.\nComment #1.2. The sampling and sequencing strategy appears confusing. The authors should\ncite and describe their Table S1 in the first paragraph of Results, and summarize the mutational\nprofiles of GBM samples analyzed, especially the IDH mutation status.\nResponse#1.2. We agree that the description in the submitted version may not have adequately\nconveyed the deliberate, stepwise nature of our experimental design.\nOur sampling and sequencing strategy followed a discovery-to-validation progression:\n1. Single-cell RNA-seq and CITE-seq were performed first, and integrated with public\ndomain datasets to establish a high-resolution cellular atlas of the glioma\nmicroenvironment and to define reference cell states.\n2. Visium spatial transcriptomics was then applied to FFPE and frozen tissue sections\nfrom the same cohort and some additional samples to map those cellular programs back\ninto their tissue architecture and spatial niches.\nXenium in situ transcriptomics was conducted lastly, as a cross validation platform, enabling\nhigher-resolution single-molecule detection of key marker sets within the tissue context.\nThis tiered design of single-cell discovery, spatial mapping, and targeted high-resolution\nvalidation was intentional. While not every assay was possible for every patient due to tissue\navailability, the overall sequencing order was consistent and conceptually aligned. Each modality\nwas chosen for its complementary strengths, with Xenium providing crucial independent\nvalidation of our single-cell–derived findings in situ.\nTo further clarify the approach, we have expanded the description of the cohort in both the\nIntroduction and Results sections. In addition, we now provide the comprehensive mutational\nprofiles of all patients in Figure 1F and Supplementary Table S1.\nComment #1.3. Fig.2, it is unclear how the authors distinguish neoplastic vs TME cells.\nResponse #1.3. We appreciate the reviewer raising this important point.\nTo distinguish neoplastic from tumor microenvironment (TME) populations, we combined\nrecursive clustering, marker-based annotation, and CNV analysis.\n1. Level-1 clustering at low resolution:\nWe applied a modified Recursive Consensus Clustering (RCC) with a resolution of 0.1\nat Level-1 to generate only broad, stable partitions and avoid over-fragmentation. This\nyielded seven macro-lineages.\n2. Marker-based annotation of Level-1 clusters:\nCluster identities were assigned using differential gene expression and canonical lineage\nmarkers:\n Cluster 1 – Microglia: P2RY12, CSF1R, TREM2, CX3CR1\n Cluster 2 – Oligodendrocytes: MOG, MAG, PLP1, MBP\n Cluster 3 – Endothelial/stromal + neoplastic admixture: vascular/stromal\n(CLU, TNR, MARCKS) together with neoplastic (SOX2, PTPRZ1)\n Cluster 4 – Microglia + granulocytes: CD74, CCL3, FCGR2B, TYROBP\n Cluster 5 – Myeloid lineage: S100A8, S100A9, S100A11, SPI1, CD68\n Cluster 6 – Lymphocytes: HLA-B, CD3E, IL7R, CXCR4, CD52\n Cluster 7 – Immune + neoplastic admixture: proliferative/neoplastic (STMN1,\nTOP2A, MKI67, UBE2C) with immune (HLA-A, HLA-C)\nThe full gene lists per cluster are provided in Table S2.\nRecursive refinement through three levels:\nEach of the seven macro-lineages was recursively subdivided up to three levels, progressively\nresolving finer cell states. The final level of clusters was used for cell type assignments.\n3. CNV validation after full cell-type assignment:\nInferCNV was then applied after the final, finer assignments. This allowed CNV\nanalysis at the level of minor cell types, as shown in Figure 2B. Neoplastic subclusters\ndisplayed consistently elevated CNV scores, while immune and stromal populations\nremained low, thus providing validation of the malignant lineages.\nThus, neoplastic vs. TME assignment was based on transcriptomic clustering and canonical\nmarkers, with inferCNV providing a validation across the entire dataset. The graphical\nrepresentation of this validation is included in Figure 2B, and the Methods and figure legend\nhave been revised accordingly.\nComment #1.4. “We were able to detect a granulocyte population from the freshly processed\nsamples, a cell type missing, from most human GBM scRNA-seq public datasets” Which\ngranulocyte population? Neutrophils? The authors need to back up their claim with specific\nmarkers and compare with other scRNA-seq datasets.\nResponse #1.4. We appreciate the reviewer’s suggestion to clarify the identity of the granulocyte\npopulation. As detailed in the CITE-seq section, the granulocytes we detected show high CD16\nand low CD14 expression, consistent with neutrophils, in line with Maas et al., Cell 2023\n(Highlighted in green, Figure 2C). We also compared our findings with the GBmap dataset\n(Ruiz-Moreno et al., Neuro Oncol 2025), where granulocytes are rare (~0.45%), underscoring\nthat our freshly processed samples enabled detection of this population often missed in other\nGBM scRNA-seq studies (Figure S2A). The corresponding description can be found in the\nrevised manuscript (page 7, lines 161–175).\nComment #1.5. What is the rationale for dividing Oligo into 4 subtypes? Have these subtypes\nbeen previously shown? Is there any disease relevance? Also, I don’t understand why the authors\nput the comparison in Fig. S7A. They should immediately show it in the main figure.\nResponse # 1.5. We thank the reviewer for this comment. Recursive clustering of\noligodendrocytes in our dataset resulted in four distinct clusters. To the best of our knowledge,\nthese subtypes have not been described previously. Their disease relevance is demonstrated in\nFigure 9, where the Oligo_2_3_2 gene signatures stratify patient survival in TCGA/CGGA\ncohorts and increased in recurrent gliomas, particularly in IDH-wildtype and MGMT-\nunmethylated tumors. In addition, the comparison previously shown in Figure S7A has been\nmoved to the main text as Figure 8A for clarity.\nComment #1.6. Fig.3D, I am not sure what the authors are showing here.\nResponse #1.6. We thank the reviewer for pointing this out. Figure 3D has been removed and is\nno longer part of the revised manuscript.\nComment #1.7. The web link provided does not appear to allow users to explore the spatial gene\nexpression.\nResponse #1.7. Thank you for bringing this to our attention. We have resolved the issue with the\nprovided web link, which now allows users to fully explore the spatial gene expression data. The\nusers can now access the data at https://www.amaranth-insights.com/ with the following\ncredentials:\nEmail: amaranthinsights@gmail.com\nPassword: password@123\nThey can view all the spatial datasets under the Spatial Glioblastoma Atlas tenant and Glioma\nmanuscript project.\nComment #1.8. Fig.6A, OPALIN and MOBP are not classical markers for myelinating OCs.\nResponse #1.8. Thank you for your comment regarding the use of OPALIN and MOBP as\nmarkers in our study. While we acknowledge that OPALIN and MOBP are not traditionally\nconsidered classical markers for myelinating oligodendrocytes (OCs), both are well-documented\nin the literature as being highly expressed during oligodendrocyte maturation and myelination.\nWe appreciate the reviewer’s suggestion and have now expanded our discussion of these markers\nin the revised text.\nFunctional studies indicate that Opalin (also known as TMEM10) plays a potentially redundant,\nnon-essential role in oligodendrocyte differentiation and maintenance of oligodendrocyte states.\nTMEM10/Opalin was first identified as a novel type I transmembrane protein and the most\nenriched transcript upon OPC differentiation into oligodendrocytes, as shown by microarray\nprofiling and genetic ablation approaches (PMCID: PMC2396630; PMC2830273).\nTMEM10/Opalin expression increases during early OPC differentiation and continues to rise as\noligodendrocytes myelinate axons in vivo. Although conventional TMEM10/Opalin knockout\nmice show normal oligodendrocyte maturation and myelination, they display learning deficits\n(PMCID: PMC5113975). Overexpression of TMEM10/Opalin in Oli-neu cells promotes myelin\ngene expression (MAG, CNP, CGT), while siRNA knock-down impairs differentiation and\nbranching (PMCID: PMC6400977). Furthermore, TMEM10/Opalin functions as a receptor for\nneuronal LGI1, orchestrating oligodendrocyte differentiation and myelination (PMCID:\nPMC11317624). TMEM10/Opalin-positive oligodendrocytes have also been detected in\nremyelinating MS plaques and EAE models (PMCID: PMC6400977), implicating them in\ndisease contexts. Human OPALIN is listed in the Human Protein Atlas as an oligodendrocytic\nmyelin paranodal protein, selectively expressed in CNS oligodendrocytes and highly expressed\nin glioblastoma, consistent with our findings.\nSimilarly, MOBP is a known myelin gene enriched in newly-formed oligodendrocytes,\nexpressed at higher levels than CNP during OPC differentiation (PMCID: PMC5221728), and\nhas been identified alongside MOG and other myelin-associated genes in single-nucleus\nsequencing studies of the human cortex and substantia nigra (PMCID: PMC7442652).\nTaken together, these data support the use of OPALIN/TMEM10 and MOBP as relevant\noligodendrocyte markers. We have revised the manuscript to include this discussion and the\ncorresponding references. In addition, for comparison and validation, we also analyzed classical\nmarkers such as PLP1, MBP, and MAG in our dataset.\nComment #1.9. What is the cellular frequency of Oligo_2_3_2 in the tumor core and periphery?\nHow does it impact the tumor cells? I am not fully convinced this is an important cell\nsubpopulation.\nResponse #1.9. We thank the reviewer for this insightful comment and for prompting us to\nelaborate on the significance of the Oligo_2_3_2 population.\nSpatial distribution\nIn our analysis, Oligo_2_3_2 was found almost exclusively in the tumor core, accounting for\n~1.9% of total cells and 24.6% within oligodendrocytes, whereas in the periphery its frequency\nwas negligible (~0.2% of total cells and 1.7% within oligodendrocytes) (Figure 2F, Figure 8A).\nWithin the core, this population was predominantly localized to perivascular regions, where it\nshowed a close spatial association with immune cells, particularly macrophages.\nClinical implications\nWe further demonstrate that enrichment of Oligo_2_3_2 is significantly associated with adverse\nclinical outcomes in two independent cohorts (TCGA and CGGA). Patients with high\nenrichment of the Oligo_2_3_2 up signature and low enrichment of the Oligo_2_3_2 down\nsignature had the poorest prognosis, whereas the reverse pattern predicted the most favorable\nsurvival (Figure 9A). In addition, analysis of the GLASS [PMCID: PMC6280138] longitudinal\ndataset revealed that Oligo_2_3_2 up signatures were further increased at recurrence, particularly\nin IDH-wildtype and MGMT-unmethylated tumors (Figure 9B), suggesting a role in therapy\nresistance and tumor relapse. Consistently, this population was primarily observed in\nglioblastoma, IDH-wildtype, but not in IDH-mutant astrocytoma or oligodendroglioma, further\nunderscoring its disease specificity and potential clinical relevance.\nWe have revised the manuscript to better contextualize these findings and emphasize their\nimplications. In the revised manuscript, we added the following passages: “The analysis revealed\nthat patients with high scores for the Oligo_2_3_2 up gene set and low scores for the\nOligo_2_3_2 down gene set exhibited significantly poorer survival outcomes… recurrent\nsamples with unmethylated MGMT promoters showed elevated expression of the Oligo_2_3_2\ngene set, suggesting a potential association with treatment resistance and tumor recurrence.” on\npage 16 corresponding to lines 420-434.\nReviewer #2 (Remarks to the Author):\nPark et al. present a comprehensive single-cell and spatial transcriptomic dataset of glioblastoma\nsamples, complementing previous studies by Ravi (2022) and Greenwald (2024). While this\ndataset, which employs CITE-seq, sc/snRNA-seq, Visium, and Xenium technologies, is a\nvaluable resource, it offers limited biological insights. The authors identified a potential novel\ncell type, Oligo_2_3_2, interacting with another new cell type, gbmEndoPeri, but further\ncharacterization is needed to understand their phenotypes and clinical implications. The proposed\nlink between neuroinflammation and the transformation of oligodendrocytes into pathogenic\nstates associated with impaired remyelination and neurodegenerative processes remains\nunsupported by current analysis.\n-We thank the reviewer for this thoughtful comment. We agree that our current data do not fully\nestablish a mechanistic link between neuroinflammation and the transformation of\noligodendrocytes into pathogenic states. We have revised the Discussion to explicitly\nacknowledge this as a limitation and to note that further functional validation will be required to\nconfirm these findings.\nDespite these limitations, our study provides several important advances. We identify and\nspatially localize a previously unrecognized oligodendrocyte subtype (Oligo_2_3_2) enriched in\nperivascular and immune-glial niches of the GBM core; we show that this population displays an\nimmune-activated profile with reduced myelin gene expression; and we demonstrate its potential\nclinical relevance through associations with patient survival and recurrence (Figure 9). By\nintegrating multiple modalities (CITE-seq, sc/snRNA-seq, Visium, and Xenium) across both FF\nand FFPE tissues, our study provides a resource that can guide future functional investigations\ninto this population and its role in GBM.\nMajor Points:\nComment #2.1. The authors provide a rich dataset of spatial transcriptomics. However, there\nhave been previously published spatial transcriptomic datasets:\n–Ravi, V. M. et al. Spatially resolved multi-omics deciphers bidirectional tumor-host\ninterdependence in glioblastoma. Cancer Cell 40, 639-655.e13 (2022).\n–Greenwald, A. C. et al. Integrative spatial analysis reveals a multi-layered organization of\nglioblastoma. Cell (2024).\nThe analysis of spatial organization leads to similar conclusions made in the studies above. The\nconclusion and the knowledge contribution to the field did not come across strongly. Can the\nauthors clarify specifically how this dataset can enable novel discovery for the community\ncompared to the other datasets?\nResponse #2.1. We appreciate the reviewer’s acknowledgment of the value of our dataset and\nshare the concern regarding explaining how our work builds upon existing studies. In the revised\nmanuscript, we now provide a clearer description of what distinguishes our dataset, for example,\nsequencing depth, the inclusion of both FF and FFPE samples, and the use of complementary\nplatforms such as Visium and Xenium, and how these features together create a resource that can\nfacilitate future discovery for the GBM research community.\nDataset Strengths and Novel Contributions\nWhile we acknowledge the value of prior studies in the field, including those mentioned in your\ncomments, our dataset provides several critical advancements:\n1) Scale and Resolution:\no Our dataset comprises 115,914 spatial transcriptomes from 32 glioma sections,\nrepresenting one of the largest GBM spatial resources to date.\no Deeper sequencing enabled detection of more genes and finer cell type resolution,\nallowing us to map 56 distinct cellular components, including novel\noligodendrocyte and vascular subtypes. By contrast, Greenwald et al. analyzed\n70,618 spots from 19 cases and identified 14 metaprograms broadly overlapping\nwith single-cell states. This difference is largely attributable to lower transcript\nand feature counts in earlier studies, which may be explained by technical\ndifferences such as the exclusive use of frozen sections (as in Greenwald,\nRavi, and Linnarson), thereby limiting the ability to reliably assign cell types.\nA detailed comparison of these datasets is provided in Supplementary File 1.\n2) Application of Xenium technology:\no We use custom gene panel with more even marker distribution to include more\nmarkers for endothelial cells, fibroblasts, and pericytes compared to other\nstudies (Data provided in Supplementary File 1, sheet 8). The number of markers\nidentified for each major cell type was compared across three Xenium panels; our\ncohort (SNUH), the Moffet dataset [PMC10699851], and the GBM-space dataset\n[https://doi.org/10.1101/2025.05.13.653495]. While all datasets captured\ncanonical neoplastic, oligodendrocyte, and immune cell markers, our cohort\n(SNUH) showed a more balanced distribution across lineages, including robust\nrepresentation of stromal (fibroblast, pericyte, endothelial) and immune (myeloid,\nlymphoid) markers. In contrast, GBM-space and Moffet displayed strong\nenrichment for neuronal and astrocytic markers but under-represented vascular\nand myeloid populations. Thus our design choice ensured broader coverage of the\ntumor microenvironment in our analyses.\n90\n80\n70\n60\n50\n40\nn\nA str o c yt\ne\nO P d\nC\ne n dr o c y N t\ne\ne o pl a sti\nc\nMi cr o gli\na\nM y el oi L\\d\ny m p h oi F\\d\ni br o bl a\\st\nP eri c y E t\\e\nn d ot h eli\\al\no\ng\nOli\nMajor cell types\nsrekram\nfo\nrebmuN\no We leverage Xenium in situ technology for subcellular spatial resolution. By\nutilizing Xenium, we demonstrate a level of granularity and cell-type resolution\nconfirming the accuracy of cell-type assignments, particularly for complex tumor\nniches. We also show that our findings, such as hypoxic regions devoid of\nvasculature and the exclusion of immune cells and oligodendrocytes in these\nregions, align with findings from Greenwald et al.'s CODEX data, further\nvalidating our conclusions.\n3) Advantages of FFPE Data:\no Preservation of Tissue Architecture: FFPE tissue fixation preserves the tissue's\nhistopathological structure, enabling us to directly correlate molecular\ndata with tissue morphology in a way that frozen sections cannot always\nachieve.\no Broad Applicability: We utilized FFPE tissue samples, the gold standard in\nclinical pathology. FFPE samples offer a much wider applicability for\nretrospective clinical studies and offer the potential to re-analyze rare, historically\narchived samples for further discovery. This makes our dataset especially\nvaluable for expanding spatial transcriptomics into clinical research.\no We did high-resolution imaging so that our data is compatible with AI-driven\nimage analysis based on H&E staining. This integration bridges molecular data\nwith histopathological context, offering significant utility for future computational\nanalyses.\n4) Integration with Genotypic Data:\no Although our cohort size is not yet sufficient for a full genotype-stratified\nanalysis, all tumor samples were comprehensively genotyped, providing a\nframework for linking spatial and molecular backgrounds. This annotation enhances the long-term value of the dataset, as it can be\nintegrated with future studies to investigate genotype-phenotype-spatial\nrelationships in GBM. We also include outcome data for future analysis.\n5) Niche-Specific Insights:\no Unlike previous spatial transcriptomic studies that largely focus on the tumor core\nand hypoxic regions, we carefully selected tumor regions to include all known\nanatomical and histological niches. This broader sampling approach provides a\nmore complete understanding of the tumor microenvironment.\no By systematically sampling across diverse histological regions, we were able to\ncharacterize all canonical GBM niches and introduce two additional\ncompartments, the immune-glial niche (IGN) and the immune-fibroblast niche\n(IFN), which were not captured in earlier atlases. These refinements provide a\nmore complete view of how stromal, immune, and glial elements organize around\nthe tumor.\nThese advancements collectively distinguish our study and dataset, providing a unique resource\nthat integrates high-resolution spatial and molecular data. This level of granularity and\ncompatibility with computational tools makes the dataset well-suited for advanced analyses,\nincluding AI-based models for predicting tumor behavior and patient outcomes. We have\nupdated the manuscript to more explicitly highlight these technical strengths and their\nimplications for advancing GBM research.\nComment #2.2. To increase value to the community and further validate findings, increasing the\nXenium dataset would be beneficial.\nResponse #2.2. Thank you for your suggestion. We expanded the Xenium dataset by profiling\nfour GBM patients (SNU21, SNU33, SNU18, and SNU25) using our custom 348-gene panel,\nresulting in the annotation of 729,307 cells across a total area of 155 mm².\nThe tissue cores\nsampled from these tumors encompassed cellular tumor regions, infiltrative tumor margins, and\nperipheral tissue, thereby enabling us to capture the full spectrum of tumor-associated\nmicroenvironments. This expansion provided subcellular-resolution validation of major and\nminor cell types and allowed us to refine niche-specific organization within the tumor.\nWe have\nupdated the manuscript accordingly, with revisions in the Methods (Xenium data analysis),\nResults (GBM landscape visualized at subcellular resolution), and Discussion sections. The\ncellular and spatial landscape of the samples can be viewed in Figures 5 and 7.\nComment #2.3. More attention is required in the figures and text. Some legends and labels are\nunclear (see section below).\nResponse #2.3. Thank you for highlighting the need for clearer legends and labels in the figures\nand text. We have carefully reviewed all figures and their corresponding legends to ensure clarity\nand accuracy.\nComment #2.4. The AUC analysis for spatial enrichment of cell types in AFs is missing in the\nmanuscript. More evidence and analysis are required to support the claims about spatial co-\nlocalization.\nResponse #2.4. We appreciate the reviewer’s comment.\nIn the revised manuscript, we have now:\n Included a description of the AUC analysis for cell-type enrichment in anatomical\nfeatures in the Methods section (under Cell-type enrichment in anatomical features).\n Performed an additional neighborhood analysis in both the Visium and Xenium datasets\nto more rigorously evaluate spatial co-localization. The methodology for this has been\ndetailed in the Methods section (under Neighborhood analysis).\n Presented the results of these analyses in the revised Figures 4 and 7, which provide\nfurther evidence supporting the spatial enrichment and co-localization of cell types within\nspecific anatomical features.\nComment #2.5. The pseudotime analysis in Figure 6C is not supported. These are non-malignant\ncells that are terminally differentiated and thus, the pseudotime framework to model trajectories\nshould not apply.\nResponse #2.5. We thank the reviewer for this important comment. In response, we have\nremoved the pseudotime analysis previously shown in Figure 6C from the revised manuscript.\nComment #2.6. The main finding of this manuscript is the novel oligodendrocyte subtype. It\nwould be good to validate the existence of this subtype with protein marker staining.\nResponse #2.6. Thank you for highlighting the importance of validating our main finding. We\nappreciate the suggestion to confirm the novel oligodendrocyte subtype Oligo_2_3_2 with\nprotein marker staining. In the revised manuscript, we now include immunohistochemical\nstaining of human GBM samples, which confirmed the presence of Oligo_2_3_2 in the tumor\ncore (Figures S8 and S9). Additionally, we validated its existence using a murine glioblastoma\nmodel, where immunofluorescence demonstrated Sox10⁺/Cre⁻ oligodendrocytes co-expressing\nOligo_2_3_2 markers (S100a1, Serpina3n, B2M) within the tumor core but not in normal brain.\nThese validation experiments are presented in Figure S9 (human GBM IHC) and Supplemental\nFigure S10 (mouse model immunofluorescence), with both images and quantitative analyses\nsupporting the robustness of this subtype identification.\nComment #2.7. Although the authors have identified Oligo_2_3_2, it is not clear from the\nanalysis how the novel subtype contributes to the pathology of GBM. Is there an association\nbetween the proportion of Oligo_2_3_2 and overall survival?\nResponse #2.7. We thank the reviewer for this question. Our analysis demonstrate that\nenrichment of Oligo_2_3_2 in the tumor core is significantly correlated with worse overall\nsurvival in two independent datasets (TCGA and CGGA), highlighting its potential contribution\nto GBM pathology (revised manuscript Figure 9A). Importantly, Oligo_2_3_2 is predominantly\nobserved in glioblastoma, IDH-wildtype, which is associated with a more aggressive clinical\ncourse, further underscoring its relevance to disease severity. Moreover, this population is more\nabundant in recurrent glioblastomas compared to primary tumors (revised manuscript Figure\n9B), suggesting a role in tumor progression and recurrence. We acknowledge that a limitation of\nour manuscript is that we do not functionally validate these potential roles of Oligo_2_3_2.\nNevertheless, its unique spatial distribution and strong association with adverse outcomes\nprovide a compelling rationale for additional functional studies in the future. The manuscript has\nbeen revised to better articulate these findings and their clinical implications.\nMinor Points:\nComment #2.8. Figure 1B: It is unclear where the cell count per cell type is coming from. Why\nis the microglia count much higher than tumour cells?\nResponse #2.8. The cell type counts shown in Figure 1B represent the major cell types included\nin our Glioma reference atlas. The earlier values were incorrect and have now been corrected in\nthe updated figure.\nComment #2.9. Figure 2A/B: Why does the UMAP look different between different levels of\nlabels?\nResponse #2.9. We applied recursive consensus clustering up to three levels. From the second\nlevel onward, each level-1 cluster was re-clustered independently, generating a new UMAP for\nthe corresponding subpopulation. At every level, UMAP coordinates were recalculated to\nrepresent the full cellular population, which explains why UMAPs at different levels appear\ndifferent. In the revised figure 2A/B, we now present only the level-3 UMAPs to avoid potential\nconfusion.\nComment #2.10. Figure 2C: Need to briefly describe what the states stand for and how they\nalign with previously reported cell states.\nResponse #2.10. Some of these cellular states are novel and have not been described in any\nprevious publications (e.g., Oligo_2_3_2); the rest are described in the results section now. We\nhave also included a detailed description of each cellular state and correspondence to cell states\ndescribed in other studies in Table S3.\nComment #2.11. Figure 2C-D: Add a statistical test to support the claim that core samples have\na higher proportion of tumor cells compared to peri samples.\nResponse #2.11. We now mark cell types with significant core–peri differences using asterisks.\nSignificance was assessed with a two-sided paired t-test across patients on per-sample\nproportions (tumor-cell fraction per type), with Benjamini–Hochberg correction for multiple\ntesting (q < 0.05). This is stated in the figure and detailed in Methods.\nComment #2.12. Figure 3D: This plot is unclear. Please make the message clearer, either\nthrough text or an alternative visualization.\nResponse #2.12. We have removed the figure from the revised manuscript.\nComment #2.13. Figure 5B: Can be moved to supplementary.\nResponse #2.13. This panel has been retained in the main figures (now Fig. 5E). It provides the\nprimary cross-modality validation, demonstrating concordance between cell-type assignments\nfrom single-cell RNA-seq and Xenium spatial profiling. This evidence is central to our claim that\nthe single-cell atlas transfers to spatial data and underpins downstream analyses\nComment #2.14. Figure 6A: “We scrutinized the differentially expressed genes (DEGs) between\noligodendrocytes from the tumor core and periphery”: Oligo_2_1 is not enriched towards the\nperiphery. What is the rationale for specifically comparing Oligo_2_3_2 to Oligo_2_1?\nResponse #2.14. In the revised manuscript, we have addressed the concern by performing a\ndifferential gene expression analysis between Oligo_2_3_2 and all other oligodendrocyte\nsubclusters (Oligo_2_1, Oligo_2_2, and Oligo_2_3_1), rather than limiting the comparison to\nOligo_2_1. This broader analysis provides a more comprehensive assessment of differences\nbetween the tumor core–enriched subcluster and the remaining oligodendrocyte populations.\nComment #2.15. Figure 6 Caption: Ensure captions are consistent with the figure.\nResponse #2.15. We have corrected the legends to align with the figures to avoid\ninconsistencies.\nComment #2.16. Figure S1: Legend missing for dot size and colour.\nResponse #2.16. The legends have been added to the respective figure.\nComment #2.17. Figure S3: Legend missing for Visium spot annotations.\nResponse #2.17. The legend has been added to the respective figure.\nComment #2.18. Figure S6B: Needs a legend. If the enrichment of Oligo_2_3_2 towards the\ncore is significant.\nResponse #2.18. We have removed the figure from the revised manuscript.\nComment # 2.19. Figure S7A: The Y-axis is missing or too small.\nResponse # 2.19. The figure is now moved to Figure 8 panel A and the axis labels are adjusted.\nComment #2.20. Figure S7B: Change the colour of the colour bar to be more easily\ninterpretable. Report FDR instead of p-value.\nResponse #2.20. We have redone the figure and are now reporting -log10 FDR values instead of\np-values. The color scheme is adjusted to show better distinction in the significance values.\nComment #2.21. In the text, it should be Neftel and not Verhaak.\n\"TC_mesnh, resembling Mes1-like cells in Verhaak classification, demonstrate an association\nwith immune infiltration, whereas TC_mesh (Mes2-like cells in Verhaak classification) are\nlinked with chronic hypoxia.\"\nResponse #2.21. We thank the reviewer for catching this error. The text in the Methods section,\n“Cell state assignment for reference single-cell data” (page 30, line 755 in the revised\nmanuscript) has been corrected to read “Neftel” instead of “Verhaak.”\nComment #2.22. “Spots from patients with astrocytoma (Ast), IDH mutant (IDHm) (SNU19,\nSNU22) were classified into LE_GM, LE_WM, IT, and CT”: Does it make sense to label\nIDHmut cases with GBM anatomical features?\nResponse #2.22. Thank you for raising this important point. We recognize that anatomical\nfeatures (AFs) were introduced in GBM. We apply them to IDH-mutant gliomas as a neutral\nspatial vocabulary to enable cross-dataset comparison, while restricting labels to BV, LE_GM,\nLE_WM, IT, and CT and explicitly excluding GBM-defining features such as necrosis and\nmicrovascular proliferation. Two observations support this choice: (i) Greenwald et al. report\nrecurrent cell-state programs in IDH-mutant gliomas (vascular, neuronal, oligodendroglial,\nastrocytic/OPC-like, myeloid) that align with AF-linked niches (vascular ↔ perivascular regions,\nneuronal ↔ LE_GM, oligodendroglial ↔ LE_WM/IT, astrocytic/OPC-like ↔ CT, myeloid ↔\nimmune-enriched zones); and (ii) the foundational anatomic transcriptional atlas of human\nglioblastoma (doi:10.1126/science.aaf2666) included a small subset of IDH-mutant cases,\nindicating that AF-oriented analyses have been considered alongside IDH status. We have\nclarified this rationale and the exclusion of GBM-specific labels in the revised text.\nComment #2.23. “This suggests that despite sharing mesenchymal gene signatures, these two\ntumor cell states may have arisen under different environmental pressures.”: How does this\nconnect with Sten Linnarsson’s GW-like and FB-like states?\nResponse # 2.23. Thank you for this insightful comment. Upon assigning our cell types to the\nLinnarsson data we found that FB-like 1 and FB-like 2 align with the mesenchymal, non-hypoxic\nmalignant state, whereas FB-like 3 corresponds to GBM-associated fibroblasts. Most GW-like\nprograms likewise map to our mesenchymal, non-hypoxic malignant state; the exceptions are\nGW-like 4 and GW-like 6, which correspond to our mesenchymal, hypoxia-associated malignant\nstate. This cross-mapping separates stromal CAFs from tumor programs and clarifies how\nsimilar mesenchymal signatures can arise under distinct microenvironmental pressures (ECM\nremodeling versus hypoxia/vascular stress). We have not included this analysis in the manuscript\nas the paper is still not peer-reviewed.\nReviewer #3 (Remarks to the Author):\nIn this manuscript (NCOMMS-24-33567-T, “A spatially resolved human glioblastoma atlas\nreveals distinct cellular and molecular patterns of anatomical niches”) Park*, Sonpatki* et al.,\nperform spatial transcriptomics with matched single cell RNA sequencing and cellular indexing\nof transcriptomes and epitopes (CITE-Seq) on FFPE embedded glioma samples. They present a\nspatial map of 56 cellular compartments based on 115,914 spatial transcriptomes from 32 tissue\nsections that could be of use to the research community. The authors suggest that their study\nprovides novel insight into the cellular architecture of GBM compared to similar studies that\nhave been published in Nature Communications and other journals (PMCID: PMC9974099,\nPMCID: PMC9950149, PMID: 35700707, PMCID: PMC10336135, PMCID: PMC9767870).\nAlthough the authors show some unique data, including some previously unreported diversity in\noligodendrocytes, there are some issues ranging from major to minor that reduce my enthusiasm\nfor the manuscript in its present form.\nThank you for your insightful feedback. Below, we provide a point-by-point response to your\ncomments.\nComment # 3.1. In general, in its current form I fail to understand the overall questions or\nhypothesis that is trying to be answered through this analysis. I understand the authors are\nclaiming they have a unique dataset, which I do not want to diminish as it could be useful to the\nresearch community. However, multiple recent papers, some listed above, have extensively\nconducted spatial transcriptomics with matched RNA sequencing.\nResponse # 3.1. We thank the reviewer for raising this important point and for allowing us to\nclarify the overall motivation of our study.\nWhile we acknowledge that several recent publications have successfully applied spatial\ntranscriptomics in glioblastoma, our single-cell sequencing and Visium profiling was initiated\nbefore the other studies were published. The central aim of our study was to understand how the\nremarkable cellular heterogeneity of GBM is organized across distinct anatomical niches and\nhow this organization contributes to cellular states. Given the other datasets and analysis focused\non tumor and immune cell populations and hypoxic niches, we set out to determine whether non-\nmalignant cell populations, in particular oligodendrocytes, endothelial cells, fibroblasts, and\npericytes, harbor novel disease-associated subtypes within the GBM microenvironment. Through\nrecursive clustering of single-cell data, we identified 21 cellular subtypes/states for these major\ncell types. We identified Oligo_2_3_2, a previously unrecognized oligodendrocyte state enriched\nin the tumor core. This population displayed an immune-activated transcriptomic profile and loss\nof canonical myelin genes, features that distinguish it from oligodendrocytes in normal brain or\nIDH-mutant gliomas. Importantly, enrichment of the Oligo_2_3_2 signature correlated with poor\nsurvival in TCGA and CGGA cohorts and was further increased at recurrence in the GLASS\ndataset, suggesting a potential role in GBM aggressiveness and therapy resistance. We further\nvalidated our findings with subcellular-resolution Xenium profiling using a fully custom gene\npanel and were able to elucidate two additional niches; Immune-glial and Immune-fibroblast.\nA further goal of our work was to extend the applicability of spatial transcriptomics to clinically\nrelevant settings. By generating data from both fresh frozen and FFPE tissues and integrating\nthese with matched single-cell multi-omic data, we demonstrate that high-resolution spatial\nmapping can be performed on routinely archived clinical material. This provides a framework for\nfuture studies to interrogate rare or historically collected specimens that were not accessible to\nearlier frozen-tissue–based atlases.\nTaken together, our study was not only intended to introduce a large and technically diverse\ndataset, but also to address fundamental questions regarding the spatial organization of malignant\nand non-malignant states in GBM and to highlight the clinical significance of previously\nunrecognized cell populations. We have revised the Introduction and Discussion to state these\noverarching questions more explicitly.\nComment # 3.2. The authors highlight their unique oligodendrocyte cell type that they have\nidentified throughout their analysis, but that is not discussed in the abstract or the title.\nResponse # 3.2. We thank the reviewer for this helpful suggestion. In response, we have revised\nthe abstract to explicitly include our main finding. The revised abstract now states: “Notably, we\nidentified a distinct oligodendrocyte subset confined to core and perivascular regions,\ncharacterized by transcriptional signatures predictive of recurrence and adverse clinical\noutcome.” (see Abstract, line 46-48)\nComment # 3.3. It is current form, it is not clear how the Olig2_3_2 population marked by\nexpression of HLA-A is novel and/or clinically relevant. What is the function of this population?\nMany manuscripts have extensively categorized and identified the OPC lineage. If the authors\nwant to focus on this, I suggest making the paper about this and providing more evidence as to\nhow this population is unique and why it is important (linked to any patient outcome, genetic\nmutation, or any other characteristic (treatment response, biological sex, etc.)).\nResponse # 3.3. Thank you for raising this important point regarding the Oligo_2_3_2\npopulation and its clinical relevance. This point is similar to comment #2.7, so we would like to\nrefer the reviewer to our response #2.7. Briefly, we have included additional analyses to explore\npotential links between the Oligo_2_3_2 population and clinical characteristics such as disease\nstatus, genetic mutations, and patient outcomes. As mentioned above, we demonstrate that the\nOligo_2_3_2 population is significantly associated with worse clinical outcomes in two\nindependent datasets (TCGA and CGGA). Furthermore, this subpopulation is more abundant in\nrecurrent glioblastomas compared to primary tumors and is predominantly observed in IDH\nwild-type glioblastomas, which are associated with a more aggressive clinical course. These\nfindings underscore its potential role in tumor progression and recurrence. These data are shown\nin the revised manuscript in Figures 9A and 9B.\nThe reviewer also asks about the potential function of the Oligo_2_3_2 subpopulation. It is\ncharacterized by high expression of immune-modulating genes including Serpina3n, S100a1, and\nGSN. While the precise function of Oligo_2_3_2 remains to be fully defined, its selective\nenrichment in the tumor core, particularly in perivascular regions, and its close association with\nmacrophages suggest a potential role in modulating the immune landscape of the tumor\nmicroenvironment. We acknowledge that further experimental investigation will be required to\ndelineate the specific mechanisms by which Oligo_2_3_2 contributes to GBM pathology, but its\nunique spatial distribution and strong association with adverse outcomes provide a compelling\nrationale for additional studies.\nAgain, we appreciate your feedback and hope these additional analyses clarify the significance of\nthe Oligo_2_3_2 population in the tumor microenvironment. We have revised the manuscript to\nemphasize these points.\nComment # 3.4. The authors suggest that this shared resource will help develop new and\neffective combinatorial therapies but do not give an example as to how. What makes this analysis\nand dataset different than what is already out there? How can this dataset provide information\nthat no other dataset can when looking for combination therapies. Not to belabor the point, this\ngets back to the authors not answering a clearly delineated question or hypothesis.\nResponse # 3.4. We thank the reviewer for this valuable comment. We agree that our initial\nmanuscript did not provide sufficient context regarding how our resource may inform the\ndevelopment of combinatorial therapies. The primary aim of our study was to identify the\nspecific combinations of cellular populations that define distinct anatomical niches in GBM. By\nresolving the cellular composition and interactions within each niche, we demonstrate that\ncritical tumor-stroma-immune interactions vary in a context-dependent manner. While we do not\nprovide direct functional validation, our hypothesis is that targeting these key interacting cell\npopulations simultaneously, rather than individually, may offer greater therapeutic benefit. In\nthis way, our atlas highlights cellular combinations that could serve as a foundation for\nformulating rational combinatorial therapeutic strategies. We have replaced the text,\n“This valuable and openly shared resource will help develop new and effective combinatorial\ntherapies to target multiple niches simultaneously to improve patient outcomes.”\nWith\n“Our publicly accessible atlas provides a high-resolution framework for investigating GBM\nspatial biology and offers a resource for identifying spatially restricted therapeutic\nvulnerabilities.” (see Abstract, line 48-49)\nWe added the following paragraphs to Discussion section including the limitations of our study\nand future work that can be done based on the findings from our atlas,\n\"Translational implications....them within their spatial and cellular contexts.\" (page 19, line 466-483\nin the revised manuscript)\nComment # 3.5. The patient classification should be based on sex (not gender) as this is the\nbiological variable.\nResponse # 3.5. We thank the reviewer for pointing this out. We agree that patient classification\nshould be based on sex as the appropriate biological variable. The manuscript has been revised\naccordingly to replace “gender” with “sex” throughout the text and figures where applicable.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nRESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have properly addressed my concerns. I have no further comments.\nReviewer #2 (Remarks to the Author):\nThis revised manuscript presents a multi-modal glioblastoma atlas integrating sc/snRNA-seq,\nCITE-seq, Visium profiling of both frozen and FFPE samples, and subcellular-resolution Xenium\ndata. The dataset is undeniably large, encompassing over 115,000 Visium spots, more than\n220,000 single cells, and nearly three-quarters of a million Xenium profiles generated with a\ncustom 348-gene panel. The authors have also genotyped all tumors, incorporated\nhigh-resolution imaging, and provided a public browser. In terms of scale, technical breadth, and\npotential clinical applicability through FFPE inclusion, the resource is both large and diverse.\nThe revision addresses most of the concerns raised in the previous round. The authors now\nmore clearly articulate how their dataset differs from and extends prior atlases such as Ravi\n(2022) and Greenwald (2024). They highlight the deeper sequencing, broader histological\nsampling across canonical niches, the advantages of FFPE for translational research, and the\nmore balanced Xenium panel design that better captures stromal and vascular components.\nThe Xenium dataset has been extended to four GBM patients, providing some subcellular\nvalidation of cell-type assignments and niche architecture. The presentation is cleaner, and the\ndiscussion has been moderated; speculative links have been either removed or clearly marked\nas limitations.\nThe central biological claim remains the identification of a distinct oligodendrocyte subtype,\nOligo_2_3_2, enriched in the tumor core and perivascular regions, characterized by\nimmune-associated transcripts and loss of myelin genes. In this revision, the authors show that\nthe Oligo_2_3_2 signature stratifies survival in TCGA and CGGA, is enriched at recurrence in\nGLASS, and associates with IDH-wildtype and MGMT-unmethylated tumors. They supplement\nthese findings with Xenium validation and orthogonal evidence from human IHC and mouse IF.\nThese additions reduce the likelihood that Oligo_2_3_2 is a computational artifact.\nThat said, the advance remains primarily descriptive. The dataset is broad and technically\nimpressive, but the biological insights are incremental rather than transformative, and several of\nthe supporting analyses remain statistically underpowered or methodologically weak. The work’s\nvalue lies chiefly in the scale and accessibility of the resource, not in mechanistic discovery. The\nfollowing issues must be addressed for the manuscript to meet Nature Communications\nstandards and for the atlas to serve as a credible and durable community resource.\n1. The differential expression analysis strategy is unclear and methodologically weak. The\nmanuscript describes per-cell t-tests with random subsampling to equalize group sizes and, in\nsome cases, “paired = TRUE wherever applicable”. It is not evident what “pairing” means at the\nsingle-cell level, and treating individual cells as independent replicates inflates degrees of\nfreedom and increases the risk of false positives. This approach also assumes normally\ndistributed data, which is inappropriate for sparse, zero-inflated single-cell counts. Subsampling\nfurther introduces unnecessary stochasticity. For a flagship resource, key claims (particularly the\ndefinition of the Oligo_2_3_2 population) should be supported by a patient-aware framework.\nBest practice would involve pseudobulk aggregation per patient per state, analyzed with edgeR\nor DESeq2, or mixed-effects models with patient identity included as a random effect.\nResponse #1: We thank the reviewer for this valuable comment and agree that our initial\ndifferential expression strategy had methodological limitations. In the revised analysis, we have\nreplaced the earlier per-cell t-test–based approach with clusterDE, a method specifically\ndesigned to address the double-dipping issue inherent to differential expression testing following\nclustering. clusterDE statistically controls for false discoveries and identifies robust, cell type\nspecific markers even when cluster boundaries are uncertain. Using this updated framework, we\nre-evaluated the two oligodendrocyte-associated cell states and redefined their marker genes.\nAs a result, the Oligodendrocyte gene set was slightly modified, with two genes (MOBP and\nPLP1) removed from the list of differentially expressed genes. The Oligo_2_3_2 gene set\nremained unchanged. These revisions have been incorporated into the manuscript, and Figure\n8B has been updated accordingly. Alternatively, differential gene expression analysis was also\nperformed using pseudo-bulk data derived from the single-cell transcriptomic dataset. DESeq2\nwas used to normalize and identify differentially expressed genes. Genes with p-value < 0.01\nand |log₂FC| > 0.8 were considered significant, with log₂FC > 0.8 classified as up-regulated and\nlog₂FC < –0.8 as down-regulated. The results for the same are included in Figure S7B.\n2. Although the prognostic analyses suggest that high expression of the Oligo_2_3_2_up gene\nset and low expression of the Oligo_2_3_2_down gene set correlate with poorer outcomes\nacross glioma subgroups, it is difficult to reconcile these results with the single-cell finding that\nOligo_2_3_2 is absent from IDH-mutant high-grade cases. If this population does not exist in\nIDH-mutant tumors, it is unclear how its signatures could have prognostic impact in these cases.\nThis raises the concern that the observed associations may be driven by confounding factors,\nsuch as general stress-related transcriptional programs, rather than by Oligo_2_3_2 as a unique\ncell type. Moreover, the current survival analyses are descriptive rather than inferential:\nKaplan–Meier curves are shown, but without hazard ratios, confidence intervals, or multivariable\nCox regression, it is impossible to determine whether the Oligo_2_3_2 signal is independent of\nestablished clinical covariates. With TCGA and CGGA data available, the authors should adjust,\nat minimum, for age, IDH status, 1p/19q codeletion, and MGMT methylation, and should state\nexplicitly how cut-points were defined (median, quartiles, or optimized). Reporting hazard ratios\nwith 95% confidence intervals, global p-values, and cohort or batch terms will be essential to\nclarify whether the observed associations reflect a truly GBM-specific oligodendrocyte\npopulation or instead broader gene expression patterns that correlate with tumor\naggressiveness.\nResponse #2: We thank the reviewer for this thoughtful and important observation. Ideally, one\nwould not expect to see Oligo_2_3_2_up signature to be high in IDH-mutant cases but given the\nlimitation of bulk RNAseq data, we cannot distinguish between co-localized expression of the\ngeneset vs. non-overlapping expression. Also, given the small IDHmutant single-cell data cohort\nand AF associated heterogeneity of this cell population it is possible that a small fraction of\nIDHmutant cases may have Oligo_2_3_2 population. Limitations of this analysis is mentioned in\nDiscussion section. We stratified the high and low groups using 2 methods with similar results.\nWe agree that the previous stratification method, which used the median ssGSEA score as the\nthreshold, led to more IDH-mutant samples being classified as Oligo_2_3_2-high. To address\nthis, we have updated our approach to define patient groups based on Oligodendrocyte and\nOligo_2_3_2 gene set ssGSEA scores using optimized thresholds. These thresholds were\nselected such that the Oligodendrocyte-high group was enriched for IDH-mutant tumors, while\nthe Oligo_2_3_2-high group predominantly represented IDH-wildtype cases (81% in TCGA and\n73% in CGGA).\nFurthermore, we have strengthened our prognostic analysis by performing multivariable Cox\nproportional hazards regression using TCGA and CGGA datasets. The models included\nclinical covariates such as age, IDH mutation status, 1p/19q codeletion, and MGMT promoter\nmethylation, along with the Oligodendrocyte and Oligo_2_3_2 ssGSEA scores as continuous\nvariables. The corresponding forest plots with hazard ratios, 95% confidence intervals, and\nglobal p-values are presented in Figure 9A, while the updated Kaplan–Meier curves based on\nthe revised stratification are shown in Figure 9B.\n3. The new AUC and neighbourhood analyses are useful additions but currently lack\ntransparency. The manuscript does not state the number of permutations performed, the criteria\nused to define significance, or whether multiple-testing correction was applied across the many\nAF×cell-type comparisons. In the neighbourhood analysis, a fixed 50-µm radius is used without\njustification or sensitivity analysis. The authors should specify permutation counts (at least\n1,000), clearly define how p-values were derived, apply appropriate FDR correction, and\ndemonstrate that conclusions are robust across a range of neighbourhood radii.\nResponse #3: We thank the reviewer for these valuable suggestions, which have been\nimplemented in the revised manuscript. Both the AUC-based and neighbourhood analyses were\nrepeated using 1,000 permutation iterations to ensure statistical robustness. The\nneighborhood/enrichment scores were calculated by taking the ratios of the original values and\nthe average across 1000 iterations of the shuffled values. For each test, empirical p-values\nwere computed as the proportion of shuffled statistics greater than or equal to the observed\nstatistic, following the standard permutation formula:\np = (count of permuted values ≥ observed) + 1/N + 1\nperm\nFalse discovery rate (FDR) correction across all comparisons was then applied using the\nBenjamini–Hochberg method.\nIn the neighborhood analysis of the Xenium dataset, grids with an area of 2,376 μm² were\ngenerated to match the spatial coverage of a single spot in the Visium dataset, ensuring\ncomparability between platforms [Page 35, Line 877-879]. The analysis pipeline is updated to\nrefine the estimation of spatial proximity, now considering only immediate neighboring grids\nwhen quantifying local spatial interactions. All methodological details, permutation parameters,\nand significance criteria have been updated in the revised Methods section and corresponding\nfigures 4E and 7D.\n4. Data availability is inadequate and falls below Nature Communications standards. The\nmanuscript states only that raw data will be made “available upon request” and provides links to\nmultiple ad hoc portals for processed data. This is unacceptable for a resource paper.\nAccession numbers must be provided for all RAW datasets together with complete metadata.\nRaw data should encompass FASTQs, Visium brightfield and IF images, Xenium raw outputs\nand segmentation masks, probe design files, and per-sample annotations. In addition, all code\nshould be placed in a versioned public repository with a DOI (e.g., via Zenodo) and must include\nexecutable scripts or notebooks that can reproduce every main figure. Without these changes,\nthe resource will not be reproducible or durable.\nResponse #4: We thank the reviewer for this constructive suggestion and have now made all\ndata publicly available in compliance with Nature Communications data-sharing standards. The\nraw single-cell RNA-seq data have been deposited in the NCBI Sequence Read Archive\n(SRA) under submission SUB15673052, and the Visium spatial transcriptomics data under\nsubmission SUB15686159. The processed spatial transcriptomic data, corresponding\nimages, and associated metadata and codes are available on Zenodo\n(10.5281/zenodo.17572905). Both datasets would be accessible upon publication of the\nmanuscript.\nFor convenient exploration and visualization, the datasets can also be accessed through the\nfollowing interactive portals:\n● SNUH-GliomaAtlas (GBM Visium): [https://gbmvisium.snu.ac.kr/]\nA spatially resolved glioblastoma atlas providing transcriptomic maps across tissue\nsections.\n● UCSC Multi-Omic GBM Browser: [https://multiomic-gbm.cells.ucsc.edu/]\nAn integrated portal for multi-omic glioblastoma datasets.\n● Amaranth Insights: [https://www.amaranth-insights.com/]\n○ (Login: amaranthinsights@gmail.com | Password: password@123)\n○ Tenant: Spatial Glioblastoma Atlas Project: Glioma manuscript\n5. Several clarifications are needed to strengthen the central conclusion around the\nOligo_2_3_2 claim. First, the per-patient distribution of Oligo_2_3_2 frequencies in core and\nperiphery should be shown, not just pooled fractions, and the defining DEGs confirmed with\npatient-aware analysis. Second, given that only four GBM patients were profiled with Xenium,\nan independent FFPE validation (either with the same panel or via orthogonal ISH/IF) would\nsubstantially improve confidence in the generalizability of this signature. Third, because\nOligo_2_3_2 is proposed as a malignant oligodendrocyte-like state, CNV score distributions\nshould be provided comparing this cluster to bona fide neoplastic and non-neoplastic\npopulations to rule out misclassification. More broadly, several of the reported cell type\nannotations are not fully convincing. For example, in Figure 4B, the BV population does not\nshow strong enrichment for the vascular signature described by Greenwald et al., whereas\nneuronal and NPC populations unexpectedly display higher enrichment. Similarly, the PNZ and\nPAN populations appear depleted of tumor cells, which is difficult to reconcile with expectations.\nThese observations raise the possibility of discrepancies in cell-type labelling that could\ncontribute to downstream inconsistencies in interpretation. To strengthen confidence in the\nannotations, the authors should provide additional validation, for instance, a heatmap of CNV\nprofiles with AF and cell type labels overlaid to assess biological coherence, as well as a dot\nplot or heatmap showing expression of key marker genes for each annotated cell type (similar to\nFigure S1 for single-cell data).\nWith the revisions outlined above, this study will provide a credible and durable resource for the\nGBM research community.\nResponse #5: We thank the reviewer for this comment.\na. The per-patient distribution of Oligodendrocyte sub-populations in core and peri is\nprovided in Figure 2F.\nb. We have now added more IHC data of NGFR, MAG, and GSN expression in\nrepresentative Oligodendroglioma, IDH-mutant and 1p/19q-codeleted (n = 3),\nAstrocytoma, IDH-mutant (n = 3), and Glioblastoma, IDH-wildtype (n = 5) samples. The\nsame is represented in Figure S10, and the images are part of the data repository.\nc. The CNV score for all the cell types and their sub-populations is provided in Figure 2C.\nAlso, the gene expression profile of all the cellular states is provided in Supplementary\nFigure 5C.\nd. The vascular signature reported by Greenwald et al. primarily captures aberrant,\nglioblastoma-associated vasculature rather than normal blood vessels (BV). This\ndistinction is supported by the enrichment pattern of the Endo_cap signature, which\nmore accurately reflects canonical endothelial expression (Figure 4C). The apparent\nenrichment of the Greenwald signature in neuronal and NPC clusters arises from the fact\nthat these vascular structures are predominantly located within gray matter regions. To\nfurther clarify this point, we have now included representative snapshots of these\nvascular features in Figure 4B of the revised manuscript. The PNZ and PAN regions\nare depleted specifically of non-mesenchymal tumor cells, while mesenchymal tumor\ncells remain present within these zones. This observation is consistent with prior reports\ndescribing mesenchymal enrichment in perinecrotic and perivascular niches. Supporting\nthis, the MES-HYP metaprogram defined by Greenwald et al. also shows strong\nenrichment in these regions, further validating the mesenchymal identity of the tumor\ncomponent within PNZ and PAN (Figure 4C). In the updated figure, the overlapping lines\ndepicting enrichment from both the Greenwald and Ravi datasets have been made\nthinner to improve visual clarity and highlight these changes more effectively. To further\nillustrate the morphological distinctions among the different vascular subtypes,\nrepresentative H&E snapshots of BV, HBV, MVP, and IGN have been included in\nFigure 4A.\ne. As suggested by the reviewer, we have included the heatmap of aggregated gene set\nexpression profiles of cellular states across the anatomical features in Figure S4H.\nReviewer #3 (Remarks to the Author):\nIn general, the authors have been responsive to the previous review and have made\nsatisfactory modifications to their manuscript to increase clarity. I have no additional issues.\nI would like to point out that part of the response to point 3.1 is not logical, the claims that work\nwas started before other papers were published has no bearing on the clarity of the message,\nwhich was the original issue.\nResponse: We appreciate the reviewer’s overall positive feedback and their acknowledgment of\nour efforts to improve the manuscript. We accept the point made. Our intention was to provide\ncontext, but we agree that this did not address the core issue raised.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nRESPONSE TO REVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nOverall, the authors have been highly responsive, and the manuscript is substantially\nimproved. I view the manuscript as acceptable pending a few minor additions/clarifications.\n- The pseudobulk DESeq2 analysis is appreciated. For transparency, it would be helpful if\nthe authors could explicitly state the design formula used (e.g., whether patient identity was\nincluded as a covariate).\nResponse #1:\nWe thank the reviewer for their comment. We have now explicitly stated the design formula\nused for the pseudobulk differential expression analysis in the Methods section.\nSpecifically, DESeq2 was run using a design formula that included both patient identity and\ncell type as covariates to account for inter-patient variability and differences in cellular\ncomposition. This clarification has been added to improve transparency and reproducibility\n(Lines 841-843).\n- The multivariable Cox models strengthen the prognostic analysis. I suggest tempering the\ninterpretation to acknowledge that some IDH-related confounding is difficult to eliminate\nfully in bulk RNA datasets.\nResponse #2:\nWe thank the reviewer for this thoughtful suggestion. We agree that, while the multivariable\nCox models adjust for key clinical and molecular covariates, residual confounding related to\nIDH-associated transcriptional programs cannot be fully excluded in bulk RNA-seq\ndatasets. We have therefore tempered the interpretation accordingly in the Results section\n(Lines 446-448).\n- While raw data deposition now meets Nature Communications standards, reproducibility\nwould be further enhanced if Visium annotations (spot-level metadata or cell-type\nassignments, when available) were documented more explicitly in the repository.\nResponse #3:\nWe thank the reviewer for this suggestion. To enhance reproducibility, we have now\nexplicitly provided Visium spot-level annotations, including anatomical feature labels and\ncell-type assignments where available, in the Source Data files deposited alongside the\nmanuscript (10.5281/zenodo.18093125). These files document spot-level metadata used for\nall downstream analyses and figure generation.\n- While the single-cell CNV profiles are clearly presented in Fig. 2B, the approach used to\nassign CNV status in the spatial datasets is less explicitly documented. If feasible, providing\na brief visualization or description analogous to Fig. 2B for the spatial data would\nstrengthen confidence in the malignant vs. non-malignant spatial assignments.\n- A few small presentation items (gene name capitalization, missing scale bars in select\nsupplemental figures) can be addressed during final revision.\nResponse #4:\nWe thank the reviewer for this helpful suggestion. We agree that explicit documentation of\nCNV assignment in spatial datasets is important for transparency. To address this, we have\nprovided spatial CNV profiles in Supplementary File 3, showcasing inferred CNV scores for\nall spatial transcriptomics spots. In addition, spatial overlays of CNV scores mapped onto\ntissue sections are available through the UCSC Cell Browser\n(https://cells.ucsc.edu/?ds=multiomic-gbm+spatial&meta=inferCNV), enabling direct\nvisualization of malignant versus non-malignant spatial regions, analogous to the single-cell\nCNV presentation in Fig. 2B.\nWe have also addressed the minor presentation issues noted by the reviewer, including\nstandardizing gene name capitalization and adding missing scale bars to the relevant\nsupplementary figures."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have properly addressed my concerns. I have no further comments.", "reviewer_2": "This revised manuscript presents a multi-modal glioblastoma atlas integrating sc/snRNA-seq, CITE-seq, Visium profiling of both frozen and FFPE samples, and subcellular-resolution Xenium data. The dataset is undeniably large, encompassing over 115,000 Visium spots, more than 220,000 single cells, and nearly three-quarters of a million Xenium profiles generated with a custom 348-gene panel. The authors have also genotyped all tumors, incorporated high-resolution imaging, and provided a public browser. In terms of scale, technical breadth, and potential clinical applicability through FFPE inclusion, the resource is both large and diverse.\nThe revision addresses most of the concerns raised in the previous round. The authors now more clearly articulate how their dataset differs from and extends prior atlases such as Ravi (2022) and Greenwald (2024). They highlight the deeper sequencing, broader histological sampling across canonical niches, the advantages of FFPE for translational research, and the more balanced Xenium panel design that better captures stromal and vascular components. The Xenium dataset has been extended to four GBM patients, providing some subcellular validation of cell-type assignments and niche architecture. The presentation is cleaner, and the discussion has been moderated; speculative links have been either removed or clearly marked as limitations.\nThe central biological claim remains the identification of a distinct oligodendrocyte subtype, Oligo_2_3_2, enriched in the tumor core and perivascular regions, characterized by immune-associated transcripts and loss of myelin genes. In this revision, the authors show that the Oligo_2_3_2 signature stratifies survival in TCGA and CGGA, is enriched at recurrence in GLASS, and associates with IDH-wildtype and MGMT-unmethylated tumors. They supplement these findings with Xenium validation and orthogonal evidence from human IHC and mouse IF. These additions reduce the likelihood that Oligo_2_3_2 is a computational artifact.\nThat said, the advance remains primarily descriptive. The dataset is broad and technically impressive, but the biological insights are incremental rather than transformative, and several of the supporting analyses remain statistically underpowered or methodologically weak. The work’s value lies chiefly in the scale and accessibility of the resource, not in mechanistic discovery. The following issues must be addressed for the manuscript to meet Nature Communications standards and for the atlas to serve as a credible and durable community resource.\n1. The differential expression analysis strategy is unclear and methodologically weak. The manuscript describes per-cell t-tests with random subsampling to equalize group sizes and, in some cases, “paired = TRUE wherever applicable”. It is not evident what “pairing” means at the single-cell level, and treating individual cells as independent replicates inflates degrees of freedom and increases the risk of false positives. This approach also assumes normally distributed data, which is inappropriate for sparse, zero-inflated single-cell counts. Subsampling further introduces unnecessary stochasticity. For a flagship resource, key claims (particularly the definition of the Oligo_2_3_2 population) should be supported by a patient-aware framework. Best practice would involve pseudobulk aggregation per patient per state analyzed with edgeR or DESeq2, or mixed-effects models with patient identity included as a random effect.\n2. Although the prognostic analyses suggest that high expression of the Oligo_2_3_2_up gene set and low expression of the Oligo_2_3_2_down gene set correlate with poorer outcomes across glioma subgroups, it is difficult to reconcile these results with the single-cell finding that Oligo_2_3_2 is absent from IDH-mutant high-grade cases. If this population does not exist in IDH-mutant tumors, it is unclear how its signatures could have prognostic impact in these cases. This raises the concern that the observed associations may be driven by confounding factors, such as general stress-related transcriptional programs, rather than by Oligo_2_3_2 as a unique cell type. Moreover, the current survival analyses are descriptive rather than inferential: Kaplan–Meier curves are shown, but without hazard ratios, confidence intervals, or multivariable Cox regression, it is impossible to determine whether the Oligo_2_3_2 signal is independent of established clinical covariates. With TCGA and CGGA data available, the authors should adjust, at minimum, for age, IDH status, 1p/19q codeletion, and MGMT methylation, and should state explicitly how cut-points were defined (median, quartiles, or optimized). Reporting hazard ratios with 95% confidence intervals, global p-values, and cohort or batch terms will be essential to clarify whether the observed associations reflect a truly GBM-specific oligodendrocyte population or instead broader gene expression patterns that correlate with tumor aggressiveness.\n3. The new AUC and neighbourhood analyses are useful additions but currently lack transparency. The manuscript does not state the number of permutations performed, the criteria used to define significance, or whether multiple-testing correction was applied across the many AF×cell-type comparisons. In the neighbourhood analysis, a fixed 50-µm radius is used without justification or sensitivity analysis. The authors should specify permutation counts (at least 1,000), clearly define how p-values were derived, apply appropriate FDR correction, and demonstrate that conclusions are robust across a range of neighbourhood radii.\n4. Data availability is inadequate and falls below Nature Communications standards. The manuscript states only that raw data will be made “available upon request” and provides links to multiple ad hoc portals for processed data. This is unacceptable for a resource paper. Accession numbers must be provided for all RAW datasets together with complete metadata. Raw data should encompass FASTQs, Visium brightfield and IF images, Xenium raw outputs and segmentation masks, probe design files, and per-sample annotations. In addition, all code should be placed in a versioned public repository with a DOI (e.g., via Zenodo) and must include executable scripts or notebooks that can reproduce every main figure. Without these changes, the resource will not be reproducible or durable.\n5. Several clarifications are needed to strengthen the central conclusion around the Oligo_2_3_2 claim. First, the per-patient distribution of Oligo_2_3_2 frequencies in core and periphery should be shown, not just pooled fractions, and the defining DEGs confirmed with patient-aware analysis. Second, given that only four GBM patients were profiled with Xenium, an independent FFPE validation (either with the same panel or via orthogonal ISH/IF) would substantially improve confidence in the generalizability of this signature. Third, because Oligo_2_3_2 is proposed as a malignant oligodendrocyte-like state, CNV score distributions should be provided comparing this cluster to bona fide neoplastic and non-neoplastic populations to rule out misclassification. More broadly, several of the reported cell type annotations are not fully convincing. For example, in Figure 4B, the BV population does not show strong enrichment for the vascular signature described by Greenwald et al., whereas neuronal and NPC populations unexpectedly display higher enrichment. Similarly, the PNZ and PAN populations appear depleted of tumor cells, which is difficult to reconcile with expectations. These observations raise the possibility of discrepancies in cell-type labelling that could contribute to downstream inconsistencies in interpretation. To strengthen confidence in the annotations, the authors should provide additional validation, for instance, a heatmap of CNV profiles with AF and cell type labels overlaid to assess biological coherence, as well as a dot plot or heatmap showing expression of key marker genes for each annotated cell type (similar to Figure S1 for single-cell data).\nWith the revisions outlined above, this study will provide a credible and durable resource for the GBM research community.", "reviewer_3": "In general, the authors have been responsive to the previous review and have made satisfactory modifications to their manuscript to increase clarity. I have no additional issues.\nI would like to point out that part of the response to point 3.1 is not logical, the claims that work was started before other papers were published has no bearing on the clarity of the message, which was the original issue.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nOverall, the authors have been highly responsive, and the manuscript is substantially improved. I view the manuscript as acceptable pending a few minor additions/clarifications.\n- The pseudobulk DESeq2 analysis is appreciated. For transparency, it would be helpful if the authors could explicitly state the design formula used (e.g., whether patient identity was included as a covariate).\nResponse #1:\nWe thank the reviewer for their comment. We have now explicitly stated the design formula used for the pseudobulk differential expression analysis in the Methods section.\nSpecifically, DESeq2 was run using a design formula that included both patient identity and cell type as covariates to account for inter-patient variability and differences in cellular composition. This clarification has been added to improve transparency and reproducibility (Lines 841-843).\n- The multivariable Cox models strengthen the prognostic analysis. I suggest tempering the interpretation to acknowledge that some IDH-related confounding is difficult to eliminate fully in bulk RNA datasets.\nResponse #2:\nWe thank the reviewer for this thoughtful suggestion. We agree that, while the multivariable Cox models adjust for key clinical and molecular covariates, residual confounding related to IDH-associated transcriptional programs cannot be fully excluded in bulk RNA-seq datasets. We have therefore tempered the interpretation accordingly in the Results section (Lines 446-448).\n- While raw data deposition now meets Nature Communications standards, reproducibility would be further enhanced if Visium annotations (spot-level metadata or cell-type assignments, when available) were documented more explicitly in the repository.\nResponse #3:\nWe thank the reviewer for this suggestion. To enhance reproducibility, we have now explicitly provided Visium spot-level annotations, including anatomical feature labels and cell-type assignments where available, in the Source Data files deposited alongside the manuscript (10.5281/zenodo.18093125). These files document spot-level metadata used for all downstream analyses and figure generation.\n- While the single-cell CNV profiles are clearly presented in Fig. 2B, the approach used to assign CNV status in the spatial datasets is less explicitly documented. If feasible, providing a brief visualization or description analogous to Fig. 2B for the spatial data would strengthen confidence in the malignant vs. non-malignant spatial assignments.\n- A few small presentation items (gene name capitalization, missing scale bars in select supplemental figures) can be addressed during final revision.\nResponse #4:\nWe thank the reviewer for this helpful suggestion. We agree that explicit documentation of CNV assignment in spatial datasets is important for transparency. To address this, we have provided spatial CNV profiles in Supplementary File 3, showcasing inferred CNV scores for all spatial transcriptomics spots. In addition, spatial overlays of CNV scores mapped onto tissue sections are available through the UCSC Cell Browser (https://cells.ucsc.edu/?ds=multiomic-gbm+spatial&meta=inferCNV), enabling direct visualization of malignant versus non-malignant spatial regions, analogous to the single-cell CNV presentation in Fig. 2B.\nWe have also addressed the minor presentation issues noted by the reviewer, including standardizing gene name capitalization and adding missing scale bars to the relevant supplementary figures."}]

Decision letter