A blueprint for local and distal invasion programs in glioblastoma

Domain: chemistry
Published: 2026-03-17
Editorial outcome: major_revision

Abstract
Glioblastoma invasion into brain parenchyma presents significant challenges for treatment but remains
poorly understood. In this study, we combine single-cell RNA sequencing, spatial transcriptomics, and
multiplexed imaging of orthotopic xenograft models to investigate glioblastoma invasion. We first screen
20 patient-derived gliomasphere models for their distal (i.e., extending to the contralateral hemisphere)
and local invasive potential in mice. We show that models with distal invasion potential are enriched with
oligodendrocyte progenitor-like cells, while models with only local invasion potential are enriched with
mesenchymal-like cells. These patterns reflect predominantly peri-axonal vs peri-vascular invasion
routes, respectively. Next, we analyze the transcriptomes of invading cells within models (compared to
tumor core) and identify programs associated with distal and local invasion. Thus, we decouple
transcriptional features associated with invasion potential from those associated with the process of
invasion. We validate our findings by spatial transcriptomics and multiplexed imaging, further describing
the spatial niche of invasive cells. Taken together, our results provide a blueprint for the invasive potential
of glioblastoma cell states and of the programs associated with invasion across different scales.

Introduction
Glioblastoma, isocitrate dehydrogenase IDH-wildtype (hereafter GBM), is an incurable primary brain
tumor that inevitably recurs despite maximal surgery and adjuvant chemoradiotherapy1–4. A major
obstacle to treatment is the topographically diffuse nature of GBM cells5,6. Gross total resection of tumors
(defined as resection without detectable residual tumor by imaging) extends recurrence-free survival7,8;
however, clinically undetectable GBM cells that spread beyond the tumor cavity invariably lead to tumor
recurrence. Importantly, invasion of the brain parenchyma by GBM cells occurs both locally and distally,
with the inter-hemispheric white matter tracts (e.g., corpus callosum) allowing GBM to extend their reach
to the contralateral hemisphere9,10. Landmark histopathology and molecular studies highlighted that
invasion of the brain parenchyma by malignant cells tends to occur along preexisting brain anatomical
structures, primarily white matter tracts (peri-axonal) and blood vessels (peri-vascular), which constitute
the so-called Scherer secondary structures9,11–14.

Dissecting the molecular features of invading GBM cells is of critical importance and may contribute to
identifying therapeutic opportunities to limit GBM recurrence. Single-cell RNA-sequencing (scRNA-seq)
analyses of GBM have identified malignant cell states recapitulating distorted versions of normal
neurodevelopment cell types, including oligodendrocyte progenitor (OPC)-like, neural progenitor (NPC)-
like, and astrocyte (AC)-like states, as well as the mesenchymal (MES)-like state, which is regulated by
the tumor microenvironment (TME)15–19. MES-like programs in GBM have been suggested to underlie
GBM invasion, given their similarity to the epithelial-to-mesenchymal transition (EMT) programs
associated with invasion and metastasis in epithelial cancers20,21. Consistently, in vitro studies have
indicated that MES-like GBM cells may exhibit migratory and invasive phenotypes22–25. In contrast, other
studies suggest the enrichment of progenitor-like transcriptional states at the periphery of GBM26,27 and
of MES-like programs at the tumor cores18,28. A recent study leveraging PDX models coupled with a gap
junction-permeable dye and in vivo microscopy proposed that GBM invasion is primarily driven by
“unconnected” GBM cells that are enriched for neurodevelopmental progenitor states (OPC/NPC-like)29.
Thus, it is plausible that more than one GBM state may be responsible for invasion, and that each invasion-
associated state may have differing interactions with non-malignant cell types and unique TME
associations.

Prior studies aiming to profile invasive cells by scRNA-seq were limited to cells from the adjacent
infiltrating edge or “rim” of bulk tumors30,31. We hypothesized that distance from primary site and route
of invasion may influence the invasion phenotype. Accordingly, we sought to analyze both local
(ipsilateral) and distal (contralateral) invasive GBM cells and to separately analyze peri-vascular invasion
and peri-axonal/neuronal invasion. To establish a controllable system for these variables, we relied on
GBM xenografts that recapitulate various patterns of GBM invasion29,32,33.

In this work, we directly profiled invasive GBM cells and their matched tumor core cells in mouse
intracranial xenografts by scRNA-seq. We first identified expression differences between distally invasive
and locally invasive GBM models, uncovering associations of cellular states with the propensity for
distinct infiltration patterns. We then analyzed the states of invading cells within invasive models and
uncovered transcriptional programs associated with route-specific invasion. We validated our findings by
spatial transcriptomics and antibody-based multiplexed imaging, and additionally explored the TME
associations of invasive cells. Taken together, we delineate the invasive potential and the programs
associated with GBM invasion through diverse anatomical routes.

Results
A pooled in vivo PDX system reveals differential capacity for distal invasion in GBM
models
We first aimed to investigate the capacity of GBM cells for distal invasion in patient-derived xenograft
(PDX) models. Previous work showed that only about 20% of GBM PDX models invade the contralateral
hemisphere32. We therefore sought to test the distal invasive capacity of different GBM patient-derived
gliomasphere models and compare the expression profiles of distally invasive and non-distally invasive
models.

We developed a pooled in vivo system to simultaneously analyze the distal invasion capacity of multiple
patient-derived GBM models. In this assay, we grew 20 GBM patient-derived gliomasphere models
separately and pooled them for stereotaxic injection into the brains of mice (10 models per pool – [Pool
#1, Pool #2], each with 2 mice replicates) (Fig. 1A). At symptom onset, the mice were euthanized, and
brains were dissected under a stereomicroscope and separated into the primary injection site and the
contralateral hemispheric invasion fractions (Fig. 1B-C). Malignant cells from each site were isolated by
flow cytometry based on GFP expression and profiled by scRNA-seq using the 10X Genomics 3’ gene
expression platform (Fig. S1A-B). Since each GBM model has a unique genetic and expression profile,
the identity of each cell could be robustly assigned to one model. For assigning cell line identity within
each pool, we utilized a combined approach, including calling SNPs34,35 and comparing overall
transcriptomes to those of gliomasphere models profiled in vitro, as previously described36 (Fig. S1C-E,
Methods). In each pool, a single model made up the majority (~70%) of detected malignant cells within
the resulting tumor core, while about half of the models were not detected, suggesting a competition
between models in their expansion capabilities (Fig. 1D). Distally invasive cells, defined as those found
in the contralateral hemisphere, were primarily derived from three highly invasive GBM models (L0627,
L0104, GS143), with 53-83% of cells from each of those models detected in the contralateral hemisphere
(Fig. 1E). A fourth GBM model (MGG75) was moderately invasive, with 11% of its cells in the
contralateral hemisphere, while the remaining models had minimal to no distally invasive cells (Fig. 1E
and S1F).

To validate our analyses and avoid the confounding effects of our multiplexing strategy, we generated
stereotactic xenografts of each of the individual four distally invasive models described above, as well as
two noninvasive models (MGG72, MGG65), and MGG123, a model with local ipsilateral invasion37.
Consistent with the pooled results, we observed distal invasion in L0627, L0104, GS143, and to a lesser
extent in MGG75, while no malignant cells were detected in the contralateral hemisphere in the other
models (Fig. 1F and S1F-G). Histological analysis confirmed tumor cell infiltration into the brain
parenchyma in each of the three top distal-invasion models (Fig. 1G). L0104 appeared to have the most
robust invasion phenotype (>80% of invasive tumor cells, Fig. 1E-F), and developed peri-neuronal
satellitosis, as observed in human GBM38 (Fig. 1G).

Between-models analysis: Distal invasion potential is correlated with an OPC-like transcriptional
program
To identify expression features of distal invasion models, we integrated our dataset with a previously
published compendium of 29 GBM PDXs assessed for their invasiveness across the hemispheres (by a
score of I-IV) and profiled by bulk RNAseq32. Distally invasive models consistently upregulated 204
genes (Supplementary Data 1), which were significantly enriched (FDR<0.05) for the OPC-like (e.g.
OLIG2, TNR) and NPC1-like metaprograms (e.g. SOX11, DLL1) (Fig. 2A and S2A). Distally invasive
models also consistently downregulated 171 genes (Supplementary Data 1), which were enriched for
MES-related genes (e.g. TIMP1, ANXA2, and EMP1) (Fig. 2A and S2A).

When examining the distribution of previously defined GBM cellular states in our models, we found that
the tumor cores of distally invasive models were depleted of MES-like cells and enriched with OPC-like,
NPC-like and AC-like cells compared to the non-invasive models (Fig. S2B). However, only the OPC-
like enrichment was statistically significant (Fig. S2B, p= 0.033, Wilcoxon rank test). Notably, out of the
ten models with sufficient data, OPC-like cells were robustly present (>5% of cells) in the tumor cores of
five models and were largely absent (<1% of cells) in the remaining five models, and this dichotomy was
consistent with the invasive capacity of the models, with only one exception (Fig. 1H). We thus conclude
that the capacity for distal invasion is linked primarily to abundance of the OPC-like state, perhaps relating
to the extensive migratory capability of normal OPCs across the developing brain39 and that distal invasion
capacity is, in fact, negatively associated with MES-like and AC-like states.

Our observed pattern is consistent with the knowledge that distal invasion primarily occurs by crossing
inter-hemispheric white matter tracts such as the corpus callosum, thus primarily reflecting peri-axonal
invasion. In line with this, antibody-based multiplexed imaging (CODEX) of a distally invasive model
(L0627) demonstrated that SOX2+OLIG2+ malignant cells invade along white matter tracts (MOG+),
consistent with peri-axonal invasion (Fig. 2B).

While distal invasion proceeds along axonal tracts and is associated with OPC-like states, cancer cells can
also infiltrate the brain via alternative routes. We hypothesized that some of the tumor-forming models
that did not have detectable distal invasion may still show local invasion along other routes, such as along
blood vessels9. Accordingly, we found direct evidence of peri-vascular invasion in two models (MGG65
and MGG123) within the ipsilateral hemisphere37 (Fig. 2C-D and S2C-E).

Within-models analysis: A consensus transcriptional program defines distally invading cells
While inter-model comparisons revealed features associated with the capacity for distal invasion, we next
sought to identify transcriptional programs specifically enriched in cells undergoing invasion. For this,
we leveraged the paired single-cell RNA-seq profiles of malignant cells from both the tumor core and the
contralateral hemisphere within each model. This allowed us to identify a set of genes upregulated or
downregulated in distally invasive cells, (i.e. found in the contralateral hemisphere). By comparing such
genes across models, we then identified a set of consensus genes that were consistently upregulated or
downregulated in distally-invading glioma (DiG) across models (Fig. 3A, Supplementary Data 2).

The DiG downregulated signature was significantly enriched (FDR<0.05) with AC-like and MES-like
genes (Fig. 3A and S3A). In contrast, the DiG upregulated signature was significantly enriched with
NPC-like (e.g. DLL3, DCX), OPC-like (e.g. OLIG1, CADM2) and cell cycle (e.g. MKI67) genes, and
included invasion-related regulators not found in the core GBM cell states, such as ZEB1—a transcription
factor recently implicated in GBM infiltration40 (Fig. 3A and S3A). Mapping individual cells to the DiG
invasion signature revealed that this signature highlights a subset of cells with limited similarity to the
OPC-like and NPC-like states and to other signatures involved in GBM invasion29 (Fig. 3B and S3B-C).
Thus, the DiG program represents a distinct transcriptional state induced during distal invasion, rather
than a simple subset of pre-existing lineages.

The enrichment of cell cycle genes in the DiG upregulated signature prompted us to examine the
distribution of cycling cells across samples. In all three distal invasion models, cycling genes were
significantly enriched in cells that upregulated the DiG signature (Fig. S3D; p < 0.005, hypergeometric
test), consistent with recent studies that highlighted the proliferation of invasive and progenitor-like GBM
cells41. This result was observed within NPC/OPC-like cells and within AC/MES-like cells, again
highlighting that the DiG signature is distinct from the previously defined Neftel states16 (Fig. S3E).

To investigate the dynamics and spatial distribution of the DiG program during the invasion process, we
turned to spatial transcriptomics using the 10x Visium platform. We profiled both hemispheres of GS143
(by both Visium and Visium HD) and integrated our data with an external Visium dataset42 including
three additional distally invasive models (n = 5 mice total; Supplementary Data 3). In PDX models,
malignant cells are of human origin, while non-malignant TME cells are of mouse origin, enabling us to
confidently distinguish malignant reads from TME reads such that spots could be annotated as TME-
dominated, cancer-dominated, or mixed (Fig. 3C).

To examine the spatial distribution of the DiG program across locations during the distal invasion process,
we segregated PDX samples of whole brain into four regions: tumor core, tumor periphery, medial, and
contralateral (Fig. 3D, Methods). Across all models, the DiG program and cell cycle were upregulated by
malignant cells as they are distancing from the tumor core (Fig. 3D and S3F-H, p<0.008 for each model,
Wilcoxon rank test), suggesting an active invasion state rather than a program expressed at the culmination
of invasion. In parallel, DiG downregulated genes gradually decrease in expression along this spatial axis.
To further examine its temporal dynamics, we analyzed an external Visium PDX gliomasphere dataset in
which the same model was profiled at different stages of the invasion process (Fig. 3D). We found that
the DiG program was highly expressed in the medial region and contralateral hemisphere during the early
stages of invasion. In contrast, at a later stage, the expression of DiG was subsequently lower within the
medial region and mostly confined to the contralateral hemisphere, suggesting that cells that have already
undertaken invasion begin to downregulate the DiG program.

Trajectory inference of single-cell transcriptomes revealed a pseudo-temporal progression from tumor
core cells through a DiG-high state and into a cycling DiG state, consistent with a model in which cells
first upregulate DiG and later also upregulate cell cycle genes (Fig. S3I-K). Accordingly, ST analysis
identified a spatial DiG and periphery DiG program, the latter reflecting malignant cells that invaded
furthest, reaching the cortical layers of the brain (Fig. S3L, Supplementary Data 4). The periphery DiG
cells were highly proliferative (17% cycling, compared to 4.8% cycling for DiG), consistent with the
neuronal microenvironment promoting the proliferation of glioma cells (Fig. S3L)43.

We next asked how the Neftel cell states16 are spatially distributed within the tumor core and across the
contralateral hemisphere. We found that the NPC-like and OPC-like states were consistently enriched in
the contralateral hemisphere in three out of four models, while the AC-like state was consistently depleted
in the contralateral hemisphere of all models (p=0.002) (Fig. 3E and S3M-N). In the GS143 model in
particular, the interhemispheric white matter marked an abrupt transition point for AC-like depletion (Fig.
S3N), suggesting that the physical constraints and/or cell type composition of white matter tracts may
result in a rapid state transition, consistent with the notion that OPC/NPC-like cells are the states endowed
with peri-axonal invasion potential44.

A signature of local peri-vascular invasion in MES-like tumors
As noted above, models lacking distal spread may instead invade locally along blood vessels, possibly
upregulating a distinct transcriptional state. To characterize the program expressed by locally-invading
glioma (LiG), we employed fluorescent-guided microdissection and isolated around 200 malignant cells
associated with vasculature within the ipsilateral hemisphere of two non-distally invasive models
(MGG65 and MGG123) (Fig. 2C-D and S2C-E). These cells were profiled by scRNA-seq, along with
around 200 cells of the tumor core of each model. As described for the DiG program, we defined
signatures consisting of genes recurrently upregulated and downregulated in LiG compared to cells from
their matched tumor cores (Fig. 4A and S4A-B, Supplementary Data 2).

The LiG up-regulated program was enriched (FDR<0.05) with AC-like genes (e.g., CLU) and OPC-like
genes (e.g., ALCAM, NLGN3) (Fig. S4B-C). However, as observed with DiG, the LiG up-regulated
program was distinct from the previously defined states and was specifically enriched in the context of
peri-vascular invasion (Fig. 4B). Additionally, cells from the tumor core in LiG models up-regulated
mesenchymal genes, including ECM genes (e.g., FN1, COL11A1, LAMB1) and ECM-remodeling
enzymes (MMP2, LOXL2). Multiplexed imaging confirmed the presence of cancer cells exhibiting local
invasion along blood vessels (CD31) upregulating AC-OPC markers (i.e., EGFR, OLIG2) as well as
ECM-remodeling (glyco-)proteins (e.g., tenascin-C, TNC) (Fig 4C-D), suggesting that LiG may rely on
ECM interactions for peri-vascular invasion. Furthermore, ST analysis revealed a local blood vessel
invasion signature (spatial-LiG) that was localized to vascular spots of the ipsilateral hemisphere at the
periphery of the tumor core (Fig. 4E, Supplementary Data 4) and correlated with the single-cell LiG
program (r=0.43, p<2.2x10-16 for MGG123, r=0.32, p=7.8x10-12 for MGG65) (Fig. S4D). Similar to
the LiG program, the spatial LiG was significantly enriched for AC-like genes (Supplementary Data 4).

Distinct cell states reflect the propensity for invasion vs. the process of invasion
While GBM invasion is often described as a single phenotype, our analyses show multiple associations
of cellular states (i.e., expression signatures) with invasion-related phenotypes, highlighting two
important distinctions. First, different routes of invasion may rely on distinct mechanisms and hence are
associated with different cellular states. Indeed, while peri-axonal invasion is positively associated with
OPC-like cells and negatively associated with MES-like cells, while inverse associations to those states
are observed in peri-vascular invasion.

Second, even within a single route of invasion, the basal states associated with the propensity for a specific
invasion route may not be the same as the states that are observed during the process of invasion, when
the cells detach from the bulk of the tumor and move into completely different niches, with an altered
milieu of interacting cells and other tissue components. Accordingly, DiG and LiG reflect unique
signatures, with only partial associations to predefined cell states. Therefore, the previously defined Neftel
cell states can consistently distinguish models favoring distal vs. local invasion (Fig. 4F left and 4G top),
but the unique invasion states more effectively distinguish the invasive vs. non-invasive cells within a
model (Fig. 4F right and 4G bottom). While the spatial DiG state captures cancer cells actively invading
along white matter tracts, once cancer cells have traversed the corpus callosum or the anterior commissure,
they may also invade along BVs where they upregulate ECM markers (TNC, HPLN1) in addition to
neurodevelopmental lineage markers (SOX2, OLIG2) (Fig. S4E), consistent with the role of ECM
remodeling in vascular invasion45.

Mapping the TME of GBM invasion
Having defined distinct transcriptional programs associated with GBM invasion, we next sought to
characterize the TME and spatial context of invasion along these routes. Using both spatial
transcriptomics and CODEX, we explored how malignant invasion states are embedded within – and
potentially shaped by – their local environment. We first examined the anatomical distribution of the
malignant invasion states. Several location-dependent associations between malignant invasion states and
anatomically-defined neuronal subsets emerged: for example, the core spatial DiG state was most
associated with subsets of inhibitory neurons in close proximity to white matter tracts, whereas distant
DiG was most associated with cortical excitatory neurons (Fig. S5A-D, Methods).

Next, we asked whether non-neuronal components of the TME differ between the ipsilateral (tumor core)
and contralateral (distal invasion) hemispheres. We observed two consistent differences between the
ipsilateral the contralateral TMEs: a shift in glial cell composition from more astrocyte-abundant
ipsilateral TME to more oligodendrocyte-abundant contralateral TME (p=0.03) and a shift in myeloid
composition from a canonical GBM TAME program (Mac) (e.g., Apoe, Trem2) to an inflammatory
macrophage/microglia program (Inf-Mac) (e.g., Ccl2, Ccl4) (p=0.02, Fig. 5A). In distally invasive
tumors, we observed recruitment of F4/80+, P2y12+ myeloid cells to cancer cells invading along MOG+
white matter tracts (MOG) (Fig. 5B). In contrast, myeloid cells of the tumor core were F4/80+ while
myeloid cells in the grey matter were P2y12+ only. This aligns with recent work showing that microglia
are mobilized ahead of invasion, guiding GBM migration along invasion routes46. Accordingly, DiG
cancer cells were significantly colocalized with Inf-Mac, while LiG cancer cells were significantly
associated with endothelial cells and pericytes (Vasc) (Fig. 4E and 5C).

The enrichment of Inf-Mac in the contralateral hemisphere and its spatial association with DiG cancer
cells led us to examine further the differences between cancer-associated macrophages in the two
hemispheres. Apart from Inf-Mac genes, contralateral macrophages also upregulated OXPHOS-related
genes, while downregulating canonical TAM genes, with no changes in the expression of microglial genes
(e.g., Cx3cr1, Selplg), indicating they are distinct from the dominant myeloid populations of the tumor
core and the normal brain (Fig. 5D and S5E-F). Interestingly, cancer cells in the contralateral hemisphere
that were coupled to Inf-Mac also upregulated OXPHOS-related genes, as well as GPX4, a member of
the glutathione peroxidase family, which protects cells from oxidative stress47 (Fig. S5F-G).

We additionally inferred interactions between distally invasive malignant cells with neurons, and peri-
vascular invading malignant cells with endothelial cells (Methods), identifying ligand-receptor pairs
upregulated by invading cancer cells and their interacting TME counterparts (Fig. S5H-J). For example,
interacting endothelial cells and invading cancer cells upregulated ECM-related genes such as Ptn-SDC3,
Fn1-ITGAV, and Bcan-EGFR.

Collectively, we resolved GBM invasion at three levels of granularity (Fig. 5E): First, malignant features
associated with overall invasion capacity between models; second, malignant features associated with
specific modes of invasion within models; and third, malignant features within those modes of invasion
that reflect specific associations with the TME. Importantly, the expression patterns associated with each
level of invasion are distinct.

Methods
Ethical statement
In vivo experimental procedures were approved by the Institutional Care and Use Committee (IACUC) at
Massachusetts General Hospital (2019NL000229) and the University of Michigan (PRO00011658) and
were performed in accordance with institutional regulations.

GBM spheroid culture.
Patient-derived human glioblastoma cells19,55,56 derived under IRB #2005P001609 were maintained on a
low-attachment surface (Corning) with Neurobasal Medium (Thermo Fisher Scientific) supplemented
with 1/2 x N2 and 1 x B27 (Thermo Fisher Scientific), 1% Penicillin/Streptomycin (Thermo Fisher
Scientific), 1.5 x Glutamax (Thermo Fisher Scientific), 20 ng/mL of EGF and 20 ng/mL of FGF2
(Shenandoah Biotechnology). Prior to multiplexed in vivo xenograft experiments, GBM spheroid models
were cultivated in T75 flasks separately and were profiled by scRNA-seq as described below.
Subsequently, 10 GBM spheroid models were mixed (1.5 x 105 cells each) in 6 well plates, grown for 48
hours, snap labeled with CMV-N16-GFP lentiviral vector for 24 hours before harvested for stereotaxic
injection. A portion of cells was further maintained in culture for 48 hours, and analyzed by flow
cytometry to determine the efficiency of the viral labeling and cell viability.

Experimental mouse model
Healthy, naive female NOD.Cg-Prkdc < scid > IL2rg < tm1Wjl > /SzJ (NSG) mice aged 6-8 weeks were
purchased from the Jackson Laboratory (Bar Harbor, ME). All mice were housed in specific pathogen-
free conditions. Housing rooms employ centrally controlled and monitored light cycles that utilize a 12-
hour light / 12-hour dark photoperiod. Individual room has remote temperature setpoint adjustment and
maintained within plus or minus 2 degrees throughout a range of 18 - 26 C. Relative humidity within the
animal suite is maintained at 30 - 70%. Stereotaxic injections were performed as described
previously19,37,57. Briefly, cells (1.5 x 105 cells for multiplexed experiment; 105 for single-model injection)
resuspended in 1 uL PBS were loaded on a Hamilton syringe with 26-gauge needle, and injected at a
speed of 0.1uL/30s-1min using the following coordinate: 0 mm posterior, 2.5 mm lateral, and 2.5 mm
dorsal to the bregma. Upon completing injection, the needle was left in place for 3 min, then withdrawn
slowly in 2 min to help reduce virus reflux. In this study, tumor size or burden was not used as a criterion
for endpoints because the tumors were grown intracranially and were not externally visible. Endpoints
were instead determined based on IACUC-approved (2019N000229) clinical symptoms caused by the
growth of intracerebral tumors, such as significant body weight loss (20% or more weight loss from
baseline weight), seizures lasting more than 10 seconds (seizures involving blank stares or shaking will
both be considered), or a persistent hunched posture with diminished activity.

Statistics & Reproducibility
No statistical method was used to determine sample sizes. The number of cells injected and profiled was
determined empirically. GBM spheroids were pooled for the above-described multiplexed in vivo assays
without regard to spheroid identity, but experiments were not randomized. No data were excluded from
the analysis, except for individual cell and spatial spot data points that did not pass quality control. Mouse
brain tissue collection, processing, and single-cell data generation were performed without the knowledge
of model identity. During data analysis, the identities of GBM spheroids used to generate transcriptome
data were not blinded.

Isolating invading cells and tumor cores in GBM PDXs
Mice were euthanized upon reaching pre-determined clinical endpoints such as significant body weight
loss, seizures, and a persistent hunched posture with diminished activity. The collected mice brains were
imaged with a fluorescent stereomicroscope with 0.63 objective magnification (Olympus MVX10, Leica
M205), which covers both the entire brain and cell morphology at single-cell resolution. Under the
microscope, the brains were observed from the ipsilateral hemisphere for tumor-peripheral invasion to the
contralateral hemisphere for distal invasion. For distal invasion, the tumor cores were dissected from the
ipsilateral hemispheres, centered around the injection site, identified by the hole created during stereotaxic
injection, and based on GFP intensity, ensuring the removal of GFP-negative outer surface of the brains.
Contralateral samples were collected from the hemisphere opposite the injection site, avoiding the midline
border regions (approximately 500 μm). During dissections of the contralateral hemispheres, we observed
lower and more restricted GFP signals compared to the ipsilateral hemisphere. For local blood vessel-
associated invasion, tube-structured GFP+ cell clusters were dissected under a stereomicroscope without
introducing any bulk tumor mass. Then, samples from each location were mechanically and enzymatically
dissociated using a papain-based neural tissue dissociation kit (Miltenyi Biotec) supplemented with 0.1%
Type I collagenase (Thermo Fisher Scientific)/PBS. The dissociated cells were first stained with calcein
Blue AM (BD Biosciences) and Zombie NIR (BioLegend) for 25 min at 4C, and with anti-mouse
CD16/32 (BD Biosciences) for 5 min. After washing cells with ice-cold 1% BSA/PBS, the cells were then
stained with PerCP anti-mouse CD45 (clone 30-F11, BD Biosciences) for 30 min at 4C. Sorting was
performed with Becton Dickinson Influx cytometer (Becton Dickinson). Side scatter (SSC) width versus
forward scatter (FSC) area and Trigger Pulse Width versus FSC criteria were used to discriminate doublets
and gate-only singleton cells. Viable single cells were identified as calcein blue AM positive and Zombie
NIR negative to low cells. We sorted viable CD45 negative/GFP positive cells from blood-vessel
associated invasion samples into 96-well plates containing 10uL of TCL buffer (QIAGEN) with 1% beta-
mercaptoethanol (Sigma Aldrich). Plates were frozen immediately after sorting and stored at -80 degree
prior to whole transcriptome amplification, library preparation, and sequencing. For multiplexed samples,
CD45 negative/GFP positive cells were sorted into 15mL tumors for subsequent 10X scRNA-seq
processing.

Single cell RNA-sequencing
Smart-seq2 whole transcriptome amplification, library construction, and sequencing were performed as
previously published16. The spheroid models in vitro were processed as a pool using TotalSeq-B/Hashtag
anti-bodies (BioLegend), described previously19. Each of the dissociated cells was first stained with
HumanTruStain FcX (BioLegend) in Cell Staining Buffer (BioLegend) for 10 min at 4C, and then stained
with barcoded TotalSeq-B anti-bodies targeting human CD298 and b2 microglobulin in a separated tube
for 30 min at 4C. We generated a single cell suspension of the pool in 0.04% PBS-BSA and loaded
approximately 9,900 single cells to the Chromium Controller with a targeted recovery of 6,000 cells.
Single cells, reagents and single gel beads containing barcoded oligonucleotides were encapsulated into
nanoliter-sized droplets and subjected to reverse transcription. Droplets were broken and the barcoded
cDNAs were purified with DynaBeads and amplified by 12 cycles of PCR (98C for 3 min; [98C for 15 s,
63C for 20 s, 72C for 1 min] x 12; 72C for 1 min). Gene expression cDNA and TotalSeq antibody derived
barcode cDNA were size-selected and separated with SPRIselect Reagent (Beckman Coulter), and then
fragmented, end-repaired, ligated with index adaptors. Quality control of the resulting libraries was
performed with the Bio-analyzer High Sensitivity DNA Analysis (Agilent). The constructed gene
expression library and barcode library were combined at 4:1 ratio and sequenced with paired-end, 28 and
55-base reads, using a NextSeq 2000 sequencer (Illumina).

Single cell RNA-sequencing data processing
Data was used in transcript-per-million (TPM) or UMI count matrices. As quality control, we excluded
cells with fewer than 2,000 detected genes. For TPM matrices, expression levels were quantified as
E(i,j) = log2(TPM(i,j)/10+1), where TPM(i,j) refers to transcript-per-million for gene i in cell j, as
quantified by RSEM58. The average number of transcripts detected per cell was less than 100,000, thus
TPM values were divided by 10, to avoid inflating the differences between detected (E(i,j) > 0) and non-
detected (E(i,j) = 0) genes, as previously described59. For UMI count matrices, expression levels were
quantified as E(i,j) = log2(1 + CPM(i,j)/10), where CPM(i,j) refers to 106*UMI(i,j)/sum[UMI 1..n,j], for
gene i in sample j, with n being the total number of analyzed genes. CPM values were divided by 10, as
described above for TPM values. We defined relative expression over the remaining cells for each study
separately, by centering the expression levels per gene, Erel(i,j) = Ei,j– mean[Ei,1...n]. Next, we
calculated the Eavg(i) = log2(mean(CPMi,1...n)+1) for each gene, and excluded genes with Eavg < 4. The
gene-cell matrix underwent dimension reduction using UMAP and Louvain clustering using the R
package Seurat (https://satijalab.org/seurat/).

Definition of gene signature scores
Cells or bulk tumors were scored for a gene signature as previously described, using the R package scalop
(https://github.com/jlaffy/scalop). Given a set of genes (Gj) reflecting an expression signature of a specific
cell type or biological function, we calculate for each cell i, a score, SCj (i), quantifying the relative
expression of Gj in cell i, as the average relative expression (Er) of the genes in Gj, compared to the
average relative expression of a control gene-set (Gj cont): SCj (i) = average[Er(Gj,i)] – average[Er(Gj
cont,i)]. The control gene-set is defined by first binning all analyzed genes into 30 bins of aggregate
expression levels (Ea) and then, for each gene in the gene-set Gj, randomly selecting 100 genes from the
same expression bin. In this way, the control gene-set has a comparable distribution of expression levels
to that of Gj, and the control gene set is 100-fold larger, such that its average expression is analogous to
averaging over 100 randomly selected gene-sets of the same size as the considered gene-set. Cells were
scored for each study separately.

Determining the identity of GBM models in multiplex experiment
First, we applied an expression-based method to assign cells from each multiplexed PDX to the individual
spheroids we profiled by scRNA-seq separately, similar to the method previously described36. In this
approach, we generated cell line signatures for each spheroid using in vitro scRNA-seq data from all of
the injected cell lines. The signatures were defined as the top 50 most differentially expressed genes when
comparing each cell line to the remaining cell lines. Next, we clustered the cells from each multiplex pool
group using Louvain clustering. We scored each cluster for the cell line signatures, where clusters with
high overall similarity (>0.5) for the cell line scores were combined. Clusters were assigned by their
cluster’s average score for the signatures.

Next, we applied the SNV-based classification using souporcell34 to identify the unique genotypes for
each cell line grown in vitro. We then applied this same method to our multiplex PDX pools. The
algorithm was used with the parameter to identify 10 clusters because we had originally injected a pool
of 10 cell lines into each mouse. After grouping individual genotype clusters with high overall correlation
(>0.5), we correlated clusters with the in vitro genotypes using Demuxafy35 to classify the cell-line identity
of each cluster. Finally, a consensus classification approach, combining the expression-based and SNV-
based methods, was used to identify PDX models.

Comparing invasive and non-invasive models
We determined each cell line’s invasion capability based on the presence of cells from that cell line in the
invasive site (contralateral hemisphere) of the multiplex experiment. Next, we generated pseudo-bulks of
the tumor cores for each cell line from the scRNA-seq multiplex experiment. In addition, we combined
these pseduobulks with a previously published dataset of GBM PDXs profiled by RNA-seq and annotated
for their invasion capability (Invasion Score between I and IV32) and compared the expression of tumors
that were invasive (Invasion Score of III-IV or models with invasion detected in our experiment) to those
that were noninvasive (Invasion Score of I-II or models with no invasion detected in our experiment).

Characterization and comparison of invasion programs to human metaprogram states
The invasion programs were defined as the top 75 differentially expressed genes between invasive cells
and tumor core cells that were conserved across models in either distal invasion (DiG) or local invasion
(LiG). For both invasion programs identified, we calculated the Jaccard similarity with each of the
previously defined human metaprograms16 and other invasion related signatures. In addition, we
calculated for each gene in each program, the correlation of its expression with the average expression of
the canonical human GBM metaprograms defined in Neftel et al 2019.

Ligand-Receptor Interaction Analysis
To investigate potential cell-cell communication pathways, we performed a ligand-receptor interaction
analysis as previously described.19 We focused on identifying interactions between invading malignant
cells (DiG or LiG populations) and the specific non-malignant cells that constitute their respective tumor
microenvironment niches: neurons for DiG cells and endothelial cells for LiG cells. Ligand and receptor
expression was assessed in our invading cell populations. We then queried for corresponding receptor and
ligand expression within a publicly available single-cell RNA-sequencing dataset of non-malignant cells
from human GBM tumors61-62. A ligand-receptor pair was considered a potential interaction if the ligand
was expressed in one cell population and its cognate receptor was expressed in the other, as previously
explained in Hara et al 2021. This approach allowed us to map the potential communication networks
between invading cancer cells and their niche-specific stromal partners.

Pseudotime Trajectory Analysis
To infer developmental trajectories and order cells along a continuous path, we performed pseudotime
analysis using the slingshot R package63. A low-dimensional embedding of the cells from each invasive
models separately was first generated via principal component analysis (PCA) on the normalized gene
expression data for each model. The slingshot algorithm was then applied to this PCA embedding to
identify the principal cellular trajectory and assign a pseudotime value to each cell, representing its
progress along the inferred developmental continuum. To identify gene programs associated with this
trajectory, we used the R package tradeSeq64. From this analysis, we identified three distinct programs in
each model that were significantly correlated with pseudotime.

Tissue Immunostaining
Collected PDX samples were fixed with 2% PFA in PBS for 16 hours at 4°C and then stored in PBS until
further processing. Thick 100um brain floating sections were prepared with a vibrating blade microtome
(Leica VT1000 S) and stored in tissue freezing media (25% Glycerol, 30% Ethylene glycol, 1.38 g/L
NaH2PO4, and 5.48 g/L Na2HPO4) at -20°C. Sections were blocked with 10% Goat serum and 0.3%
TritonX-100 in PBS for 2h and immunostained for brain endothelium using anti-CD31 (MEC13.3; BD
Biosciences) antibody overnight at 4°C, followed by secondary antibody AF647 anti-rat IgG (Invitrogen)
for 2h at room temperature and DAPI for nucleus staining. The images were acquired with laser-scanning
confocal microscopy (Zeiss).

CODEX (multiplexed imaging)
The antibody panel was designed to capture the human and murine cell types and states of interest
(Supplementary Data S5). With the exception of the pre-conjugated CD31 antibody from Akoya, all
antibodies were validated by conventional immunofluorescence on control tissues, including human
glioblastoma and mouse control brain samples. For validation, fresh-frozen sections were fixed in 4%
paraformaldehyde (20 min), blocked with 4% BSA and 0.25% Triton X-100 (30 min), incubated with
primary antibodies overnight at 4 °C, and subsequently with secondary antibodies for 2 h at room
temperature. All antibodies were further validated post-conjugation in control samples using multiple
CODEX runs, allowing adjustment of antibody concentrations and exposure times to optimize signal-to-
noise ratios. CODEX experiments were performed on the PhenoCycler-Fusion and imaged on the
PhenoImager (Akoya Biosciences) following the manufacturer’s instructions.

Visium Sample Preparation
Samples were flash frozen by liquid nitrogen and embedded in cold OCT on dry ice (Scigen OCT
Compound, #4586). The RNA quality of each sample was evaluated by Tapestation (Tapestation RNA
Screen Tape, Agilent) following RNA isolation (Zymo Quick RNA MicroPrep Kit, #ZR-R1051). Fresh
frozen samples were sectioned at 10um thickness with a cryostat onto 10X Visium Spatial
Transcriptomics slides (Visium Spatial Gene Expression Slide and Reagent Kit, PN-1000184) according
to the manufacturer’s instructions.

Visium H&E staining and imaging
Tissue sections on Visium slides were first fixed in methanol (Millipore Sigma #34860) followed by an
aqueous eosin-based H&E protocol according to manufacturer’s instructions (10X Visium Methanol
Fixation, H&E Staining, and Imaging Protocol CG000160). Brightfield imaging was performed using a
wide-field Leica DMIi8 inverted microscope (Leica-microsystems CMS GmbH Germany) equipped with
a DFC310FX color camera. Images were acquired with a 10x/0.25 dry objective and stitched by Leica
Application Suite X software. Image post-processing was performed with Fiji version 2.3.1.

Visium cDNA synthesis and library generation
Following imaging of H&E staining, permeabilization (9 minutes, as determined by Visium tissue
optimization experiment) was carried out on the Visium slide to capture mRNA released from the tissue.
cDNA synthesis and library generation were performed with the Visium Spatial Gene Expression Slide
and Reagent Kit (10X Genomics). Quality control and quantification of the resulting dual-indexed
barcoded libraries were performed with Agilent TapeStation and by qPCR (NEBNext Library Quant Kit
for Illumina, New England Biolabs).

Tissue preparation and sectioning for Visium HD
Brains from GS143 PDX models were extracted and fixed with 10% formalin for 16-24 hours at room
temperature, then switched to 70% ethanol. Paraffin blocks were generated and adjacent or near-adjacent
5 µm sections were taken from the FFPE tissue blocks following the Visium HD FFPE Tissue Preparation
Handbook (CG000684, Revision B).

RNA quality assessment for Visium HD
RNA quality was assessed by scraping off two adjacent 5 µm sections, performing tissue slide
deparaffinization (Qiagen, cat #19093), extracting RNA following the AllPrep DNA/RNA FFPE
Handbook (Qiagen cat #80234) and calculating the DV200 score. Briefly, mean RNA fragment size is a
reliable indicator of RNA quality. As such, the percentage of total RNA fragments >200 nucleotides
(DV200) were measured for RNA quality assessment upstream of library preparation. Following the
recommendation of the Visium HD FFPE Tissue Preparation Handbook (CG000684, Revision B), tissues
were used with the Visium HD Spatial Gene Expression if they had a DV200 of >30%.

Visium HD CytAssist Spatial Gene Expression for FFPE
5 µm sections were cut from previously mentioned formalin-fixed paraffin-embedded blocks and mounted
onto Visium HD Cassettes. H&E staining and imaging were performed following the Visium HD FFPE
Tissue Preparation Handbook (CG000684, Revision B). Samples were then processed and sequenced
following the Visium HD Spatial Gene Expression Reagent Kits User Guide (CG000685, Revision B).
Briefly, the region of interest in the slides was subjected to heat- and solvent-mediated deparaffinization,
decrosslinking, and tissue permeabilization. Probe hybridization, ligation, and probe extension were
carried out on-slide to capture transcriptomic information while preserving spatial architecture. Following
second-strand synthesis and cDNA amplification, libraries were prepared according to the Visium HD
protocol, indexed, and sequenced on an Illumina platform following the Visium HD Spatial Gene
Expression Reagent Kits User Guide (CG000685, Revision B) to generate spatially resolved
transcriptomic profiles.

Visium Sequencing
Final library QC was performed by Tapestation. Pair-end dual-indexed final libraries were diluted to
1.8nM, pooled, and denatured prior to sequencing on Novaseq (Illumina) using the Novaseq SP 100 cycles
sequencing kit (Illumina) with 1% PhiX and the following sequencing parameters: Read 1 – 28 cycles,
Read 2 – 90 cycles, Index 1 – 10 cycles, Index 2 – 10 cycles.

Visium Alignment and data pre-processing
Alignment of FASTQ files to a concatenated human (GRCh38) and mouse (mm10) transcriptome
(GRCh38-mm10-2020A), UMI counting, and spot barcode filtering were performed using SpaceRanger
(version 1.1, 10X Genomics). Alignment between positionally barcoded Visium spots and tissue images
to obtain spatial coordinates necessary to generate spatial maps was performed using Loupe Browser
(version 5.0.1, 10x Genomics). Expression levels were quantified (separately for mouse and human genes)
as Ei,j = log2 (1 + CPMi,j/10), in which counts per million (CPM)i,j refers to
106 × UMIi,j/sum(UMI1...n,j), for gene i in sample j, with n being the total number of analyzed genes.
The average number of UMIs detected per spot was less than 100,000; thus, CPM values were divided by
10 to avoid inflating the differences between detected (Ei,j > 0) and undetected (Ei,j = 0) genes as
previously described59.

For each spot, the number of counts was used as a proxy for sample quality. Spot QC filtering was
performed separately for mouse and human reads. For each species, spots with fewer than 500 counts
and/or expressing more than 10% mitochondrial genes, another proxy for low quality, were filtered out.
For mouse and human separately, the top ~7,000 most highly expressed genes were retained, and centering
was performed per sample in order to define relative expression values by subtracting the average
expression of each gene i across all k spots: Eri,j = Ei,j − average(Ei,1...k), where Er represents relative
expression values.

Five GBM PDX Visium mouse brains (spanning four invasion models) that were profiled by Manoharan
et al.42 were included in our analysis (Supplementary Data 4). The same QC filtering and pre-processing
were performed. No batch integration was required since our analysis is generally performed per sample,
and only the results of various analyses are combined across samples.

Segmentation by hemisphere
In some cases, an entire brain was sectioned onto a single Visium capture area. Therefore, in order to
perform analysis comparing the ipsilateral (tumor core) and contralateral (invading) hemisphere, it was
necessary to assign all tissue-covered spots to a hemisphere. We performed segmentation and annotation
of the hemispheres using the “createSegmentation” function in the SPATA2 R package, which provides
a GUI for manual segmentation and then returns a vector of spot assignments based on the
segmentation60.

Scoring and annotating Visium data by malignancy fraction and by single-cell invasion programs
First, the malignancy fraction of each spot was inferred by the ratio of human/mouse reads. Second, spots
could be further binned as “Cancer,” “TME,” and “Mixed” based on the distribution of human/mouse
reads. Spots were then scored for the single-cell invasion programs (using only human reads) using the
scalop R package (https://github.com/jlaffy/scalop) as previously described16. Given a set of genes (Gj)
reflecting an MP corresponding to a cell state or cell type, we calculate for each spot i, a score, SCj(i),
quantifying the relative expression of Gj in spot i, as the average relative expression (Er) of the genes in
Gj, compared to the average relative expression of a control gene set (Gj cont): SCj (i) = average[Er(Gj,i)]
– average[Er(Gj cont,i)]. The control gene-set is defined by first binning all analyzed genes into 30 bins
of aggregate expression levels (Ea) and then, for each gene in the gene set Gj, randomly selecting 100
genes from the same expression bin. In this way, the control gene set has a comparable distribution of
expression levels to that of Gj, and the control gene set is 100-fold larger, such that its average expression
is analogous to averaging over 100 randomly selected gene sets of the same size as the considered gene
set.

Spots were annotated for the spatial DiG program based on the difference between the spatial DiG
program and the DiG program. Annotation thresholds were then set based on the distribution of the
differences, where a spot was classified as “spatial DiG” if the difference between the scores was >=0.5.
For spot annotation to Neftel states and other malignant programs, spots were assigned to the program for
which it scored most highly with a minimum difference of at least 0.2 between the highest-scoring
program and the second highest-scoring program.

Defining malignant signatures from spatial data
In addition to annotating the Visium PDX data by the single-cell derived invasion programs, we also
performed unsupervised analysis to derive malignant signatures directly from the spatial data in several
different ways using only the human (i.e., malignant) genes and expression data. First, we considered
differentially expressed genes (DEGs) between the ipsilateral and contralateral hemispheres per model.
We then derived consensus signatures for ipsilateral and contralateral malignant programs based on the
recurring DEGs across models. Second, we defined gene signatures based on DEGs between malignant
spots colocalized with a specific TME cell type (for example, vascular or macrophage) and malignant
spots not colocalized with a specific TME cell type. Third, we derived consensus gene signatures as
previously described28. Briefly, for each sample (using only human reads), Leiden clustering was
performed on the SNN graph (implemented with Seurat version 4.3.0). Per sample gene programs were
defined per Leiden cluster by differential expression analysis based on the top 50 most differentially
expressed genes by the Wilcoxon Rank Sum Test with a p value of <0.005. The per sample gene programs
were then clustered based on their Jaccard similarity and consensus gene signatures were defined from
the top 50 recurring genes per cluster as defined by Jaccard similarity.

Defining nonmalignant cell type signatures from Visium data
In order to annotate spots by their dominant TME cell type, we first redefined single-cell and spatial
nonmalignant TME cell type signatures from the Visium PDX data as follows: spots were scored for the
original signatures (mouse reads only). Genes were retained in the cell type signature if they had a
correlation >0.2 with the overall program score for that signature.

Next, we derived consensus signatures for neuronal subsets/anatomical locations. First, we generated per
sample Leiden-derived gene programs as described above. Second, we performed NMF on each sample
separately to capture continuous patterns of gene expression variation, as previously described28. Negative
values in each centered expression matrix were transformed to zero. To minimize the influence of
selection of an individual k parameter, we ran NMF with multiple k values ranging from 2-11. Each NMF
program was summarized by the top 50 genes based on NMF score. Consensus signatures corresponding
primarily to neuronal subsets and anatomical locations were derived by jointly clustering per sample NMF
and Leiden gene programs collected from all samples by their Jaccard overlap and selecting the top 50
most recurrent genes across samples per cluster. Finally, in several cases where we observed clear
subclusters within a cluster, we then performed hierarchical clustering, which resulted in further splitting
neuronal clusters to additional subclusters.

Spots were given two TME annotations: one based on the dominant non-neuronal TME score (i.e., Mac,
Oligo, Astrocyte, Inf-Mac, or Vasc) and an additional annotation based on neuronal subtype/anatomical
location). For each annotation category, in order for a spot to be assigned to a given program, the spot
needed to score at least 0.25 higher for that program than any other program

Colocalization and spot enrichment analysis
For colocalization between malignant states and nonmalignant cell types, only mixed composition spots
(defined as at least 20% malignant and 20% nonmalignant) were considered. First, p-values (calculated
by hypergeometric test) and fold-enrichment values (defined as observed colocalization/expected
colocalization) were calculated per sample. Finally, a mean fold-enrichment value for each pair of cell
states/cell types across all samples and effect size (defined as the proportion of samples for which a given
colocalization was significant) were calculated.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Utilizing single-cell and spatial transcriptomics of PDX glioma sphere systems in mice, the authors determined the cell states and invasion signatures of glioma cells in the ipsilateral and contralateral hemispheres. Distant glioma cells with a peri-axonal migration route (CLAT-CAI) predominantly resembled OPC and NPC-like cells that were also proliferative whereas locally invading cells with a mesenchymal gene signature preferred to migrate along blood vessels (BV-CAI). Conceivably, glioma cells with the specific CAI invasion signatures differ from those of the primary tumor, likely due to under-representation within the mass and/or selective induction during the invading process.\n\nThe use of 20 different PDX gliomas are a substantial strength. It is, however, not clear how a pool of 10 different gliomas can be clearly segregated by RNA expression; specifically, since glioma cells may change their expression profile during growth and invasion. Although not necessarily critical for this analysis here, but barcoding of the glioma populations would have been more convincing, specifically since the authors conducted scRNAseq and ST.\n\nThe observation that glioma cells can move along basement membranes, including vascular basement membranes and neurons, has been known for a long time, but that this correlates with a specific cell state is less well understood. There exist publications demonstrating that mesenchymal GBM cells prefer a perivascular invasive mode. The analysis overall is quite informative determining CLAT-CAI and BV-CAI invasion signatures. It is still not evident though how and why gliomas cells select for a specific invasion program if all the different tumor subtypes are present in the primary tumor. It is also less clear why perivascular invasion is only local and not also contralateral given that the brain is highly vascularized. Is it then possible to change the invasion mode by changing the cell state? Would OPC/NPC-like gliomas in general be more infiltrative than mesenchymal GBM; would these spread more upon recurrence?\n\nOverall, the authors demonstrate that the gene signatures are differently enriched in the distinct invasion programs but the specific mechanisms that lead to and cause these invasion modes still need to be determined. This will be instrumental in blocking these invasive programs.", "reviewer_2": "Chanoch-Myers et al\nThis manuscript by Chanoch-Myers et al describes experiments aimed at understand the gene expression patterns associated with two distinct types of invasion of gliomas in vivo. They use a series of 20 human PDX/organoid models with different propensities to invaded into the mouse brain. Through a tour de force of technologies they have isolated tumor cells that invade from the central bulk of the tumor across to the other hemisphere either along the perineuronal pattern, or adjacent to the blood vessels in a perivascular pattern. Using a combination of bulk and spatial single cell sequencing, they have identified expression for the cells that invade in both of these patterns. The cells invading along neurons show a oligo progenitor expression pattern, while the ones that invade in a perivascular manner show a mesenchymal expression pattern.\n\nIt is certainly true that brain invasion is one of the reasons that gliomas are so difficult to treat, certainly surgically. The results are also very believable, and similar to those of other glioma modeling systems. Not that it would be necessary, but one would presume that these two types of invasion would be seen in mouse, rat and swine glioma models and likely have largely overlapping gene expression signatures.\n\nI guess what might help the manuscript would be more insight of how this connects with the behavior of the human disease.\nFor example:\n1. For those 20 lines that were used and had different ratios of MES- vs OPC-like invasion, what did the MRI scans of the patients look like with respect to the actual invasion pattern of the glioma?\n2. Or since the ratio of OPC vs MES is well studied with single cell sequencing, are there correlations between MRI scan patterns of invasion and the ratio of OPC vs MES cells in the tumor\n3. Also, is there any thought around whether either type of invasion is more or less responsive to therapy? XRT for example. If that were the case, it might suggest that post radiation, the invasion type might switch to a more MES types if the OPC invading cells were more likely to die.", "reviewer_3": "Chanoch-Myers, Hara, Greenwald et al present a detailed transcriptomic characterization of glioblastoma invasion. The authors demonstrate differential contralateral and ipsilateral invasive potential of various patient-derived gliomaspheres, and characterize the molecular identities of models primed for metastasis using single cell RNA-Seq. They find that distal invasion potential is associated with OPC-like cells and gene expression, while local invasion is enriched for MES-like cells, but interestingly downregulates MES-like genes. The authors further confirm these results using pixel-based spatially-resolved RNA-Seq and delineate changes in cellular composition of the neighboring tumor TME.\n\nThe underlying question in this study is an important and intriguing one, and the methodologies employed were impressive. In particular, I found the use of a gliomasphere screen to identify models primed for distant invasion to be innovative. It provides a useful resource to the community and identified novel expression programs which govern peri-axonal vs. peri-vascular invasion routes. However, many of the conclusions were based on shallowly tested hypotheses, and several claims were largely speculative and overdrawn. In general, the manuscript would also benefit from a more coherent organization and writing. I include specific comments below.\n\nNote: The authors provided neither line numbers nor page numbers. I therefore quote the text and page number as it appears to me, where applicable. Please include this critical information in future revisions.\n\n<Major Concerns:>\n1. One of the main novelties the authors stress is the decoupling of transcriptional features associated with invasion potential (i.e. basal features of the primary tumor) vs. cells actively undergoing invasion. However, this conclusion is simply untrue and not tested in the current set of experiments, at least in the case of contralateral invasion. I see two major points that need to be addressed.\n\na. There are several interesting spatial-temporal dimensions of GBM invasion from a biological perspective. The primary tumor presumably first establishes itself at the injection site, which constitutes the basal cellular composition, as the authors suggest. In a second phase, subpopulations of GBM migrates to distant regions on the contralateral or ipsilateral side, while the remaining tumor corresponds to an altered composition. Lastly, invaded cells proliferate (as the authors identify) and may take on additional states. All the authors did was to compare gene expression signatures between the primary tumor site with metastases at a timepoint after invasion has already occurred. This in no way represents cells actively undergoing invasion to the contralateral side. The ipsilateral BV-CAI cells were actively undergoing invasion which was interesting to see. But this made comparisons more confusing as CLAT-CAI and BV-CAI were in fact two completely different populations in two stages of invasion.\n\nRelated to this, it is unclear what was actually sequenced in this set of experiments. The results section note that dissections were separated into the primary injection site and invasion fractions. But from the methods section it seems that the two hemispheres were simply separated with the midline removed. This needs to be clarified. In either case, it seems unlikely that the cells actively undergoing invasion would even have been sequenced. In the case of the latter, how would primary tumor cells be differentiated from ipsilaterally migrated cells?\n\nb. The authors should perform pseudotime trajectory analysis to identify stages of invasion and its associated molecular features, and if there are in fact multiple lineages of cells during this process.\n\n2. The authors need to perform single molecule FISH (i.e. RNA-Scope or HCR) to confirm the key findings from spatial transcriptomics with single cell resolution. The authors should confirm a selection of disrupted genes and their associated cell-types. In particular, it would be interesting to see if MES-like genes are downregulated in MES-like cells or other fractions.\n\n3. The authors interestingly propose peri-axonal vs. peri-vascular interactions (suggested by previous literature), but do not go on to characterize or confirm this. The last section delves into the tumor TME, but in alterations in immune cell compositions post-invasion. It would be interesting to perform cell-cell interaction inference to probe key players involved in the process of invasion, both from the side of GBM cells and host cells. This would be more relevant to the current study in my opinion. Additionally, previous studies identified neuronal activity to be a critical source for driving GBM invasion. Do the authors observe similar synaptic alterations in the tumour TME, and are there changes in GBM composition at different stages of invasion? This would be another interesting question that could be answered with a pseudotime analysis but which is completely missing from the current study.\n\n<Minor Concerns:>\n1. On page 5: \"Invasive cells...primarily derived from three highly invasive GBM models (L0627, L0104, GS143), with 54-78% of cells from each model being invasive\". It is unclear where these numbers are coming from. The highest range noted here is 78% but L0104 in Figure 1E,H is clearly over 80% in the [% invading cells] axis. Additionally, these numbers do not seem to correlate to the number of primary tumor vs. invasion cells presented in Figure S1E, where L0104 seem to have almost entirely invasion cells, while L0627 also appear to have a lot more than 50% invading cells. Please clarify where the discrepancies are coming from, and if it is from overlapping datapoints in the UMAP, this needs to be corrected in all figures.\n\n2. Please comb through the figures to include all necessary information for interpretation. For example, Figure S1B is missing cluster labels which is critical for interpretation and Figure 2C is missing scale bars.\n\n3. Page 8: \"Thus, the distinction between primary tumor cells and invasive cells lies in the frequency of the CLAT-CAI signature rather than a coherent difference across all cells\" - it is unclear what this sentence means.\n\n4. Page 9: \"...our analyses show multiple associations of cellular states (i.e., expression signatures) with invasion-related phenotypes...(Fig 3C). How does Figure 3C show multiple cellular states?\n\n5. If possible, the authors should provide higher resolution confocal images of vessels and the corpus callosum during invasion, instead of the blurry and over-saturated images in Figure 2C where it is difficult to tell what one is looking at.\n\n6. There are some inaccuracies in terminology. Please correct spatial genomics to spatial transcriptomics as the manuscript does not utilize the former. Please also provide the spot sizes clearly in the results and figure legends, and include the corresponding H&E counterstain images in the supplementary figures."}, "author_response": "Overview\nWe thank the Editor and Reviewers for their thoughtful and thorough evaluation of our manuscript.\nTheir insightful feedback raised important questions and offered valuable suggestions to enhance\nthe impact of our study. In response, we have incorporated extensive additional data and analyses\nthat substantively reinforce our conclusions, as detailed in the point-by-point responses below.\nKey highlights in our revision are:\n New spatial data, including spatial proteomics (CODEX) and spatial transcriptomics\n(Visium HD), supporting our prior conclusions.\n Barcoding experiment and analyses combining spatial transcriptomics and trajectory\nanalysis of scRNA-seq data to validate that the Distally invading Glioma (DiG) signature\nis expressed by invading cells.\n Niche-labeling experiment and analysis of microenvironment interactions of peri-\naxonal/peri-vascular invasive GBM cells.\n Analysis and discussion on the relevance of the invasion states for recurrence.\n Analysis and discussion on glioma selection for a specific invasion program.\nOverall, the key contributions of this work are:\n Models with distal invasion potential are enriched with oligodendrocyte progenitor-like\n(OPC-like) cells, while models with local perivascular invasion potential are enriched with\nmesenchymal-like (MES-like) cells.\n Identified novel programs associated with distal and local invasion.\n Validated our findings by spatial transcriptomics and proteomics.\n Analyzed interactions between invasive GBM cells and their microenvironment.\nReviewer #1 (Remarks to the Author): Expert in glioblastoma, tumor microenvironment and\nmetastasis\nUtilizing single-cell and spatial transcriptomics of PDX glioma sphere systems in mice, the\nauthors determined the cell states and invasion signatures of glioma cells in the ipsilateral and\ncontralateral hemispheres. Distant glioma cells with a peri-axonal migration route (CLAT-CAI)\npredominantly resembled OPC and NPC-like cells that were also proliferative whereas locally\ninvading cells with a mesenchymal gene signature preferred to migrate along blood vessels (BV-\nCAI). Conceivably, glioma cells with the specific CAI invasion signatures differ from those of the\nprimary tumor, likely due to under-representation within the mass and/or selective induction\nduring the invading process.\n1-1. The use of 20 different PDX gliomas are a substantial strength. It is, however, not clear how\na pool of 10 different gliomas can be clearly segregated by RNA expression; specifically, since\nglioma cells may change their expression profile during growth and invasion. Although not\nnecessarily critical for this analysis here, but barcoding of the glioma populations would\nhave been more convincing, specifically since the authors conducted scRNAseq and ST.\nWe appreciate the reviewer’s comment and have now strengthened our methodology to segregate\nglioma populations, ensuring more rigorous identification of glioma PDX models and programs.\nIn our previous manuscript, we identified pooled glioma models by correlating their overall\nexpression profiles with those of the same spheroids profiled in vitro, a computational method\nused previously and demonstrated to be robust1. However, we understand the reviewer’s concern\nthat cells may alter their expression profile when grown or invading in vivo. Therefore, we have\nadded another approach to improve the identification and separation of PDXs by combining our\ninitial transcriptome analysis with a genetic-based method that analyzes single-nucleotide variant\n(SNV) profiles from the scRNA-seq data (Fig. R1A).\nFirst, we applied the previous expression-based method to assign cells from each multiplexed PDX\nto the individual spheroids we profiled by scRNA-seq separately. In this approach, we generated\nsignatures for each spheroid using in vitro scRNA-seq data from all of the injected models. The\nsignatures were defined as the top 50 most differentially expressed genes when comparing each\nspheroid model to the remaining spheroid models. Next, we clustered the cells from each multiplex\npool group using Louvain clustering. We scored each cluster for the model-specific signatures,\nwhere clusters with high overall similarity (>0.5) for the spheroid model scores were combined.\nClusters were assigned by their cluster’s average score for the signatures (Fig. R1B).\nNext, we applied the SNV-based classification using souporcell2 to identify the unique genotypes\nof each spheroid model grown in vitro. We then applied this same method to our multiplex PDX\npools. The algorithm was initially used with parameters set to identify 10 genetic clusters, because\nwe had originally injected a pool of 10 spheroid models into each mouse. However, genotype\nclusters with high overall correlation (>0.5) were merged, since some models were not recovered\n2\nin the final scRNA-seq data. We then correlated clusters with the in vitro genotypes using\nDemuxafy3 (Fig. R1C) to assign each genetic cluster into a cell-line identity.\nFinally, a consensus classification approach, combining the expression-based and SNV-based\nmethods, was used to identify PDX models (Fig. R1D). Both methods yielded the same\nclassification, with one exception: in Pool 1 Mouse 2, the SNV-based method only detected three\ngenetic clusters, but one of them was clearly separated further by gene expression into two clusters,\neach with high similarity to a different in vitro spheroid model (MGG 75 and MGG72). We thus\nclassified them based on expression data as distinct models.\nImportantly, regarding the reviewer’s concern that invasion and growth may dynamically affect\nexpression such that expression-based assignment might be unreliable, in both PDX pools,\ngenotypes identified from clusters in the invasion sites and primary tumor sites were highly\nconsistent in the case of the three invasive models (Fig. R1D, GS143, L0627, L0104). This\ndemonstrates the robustness of our previous RNA-based demultiplexing strategy and its ability to\nreliably segregate glioma populations regardless of invasion. Thus, while invasion does affect gene\nexpression, there is a core expression profile of each model that persists during invasion, growth\nand other perturbations. We have now updated the methods and figures of the text to the consensus\nclassification approach (Figs. S1B-D in the updated manuscript). While we acknowledge that\nbarcoding would have been an alternative approach, our genetic-based demultiplexing provides\nstrong evidence for the clear segregation of the glioma populations in our study.\nglobal gene expression approach and SNV-based approach. B) Heatmaps show the scores of each Louvain\ncluster (y axis) from each mouse (2 pools, 2 mice) for each in vitro cell line signature (x axis). C) Heatmaps\nshow the correlation of the genotypes for each souporcell cluster (y axis) from each mouse with the\ngenotypes of the in vitro cell lines (x axis). D) UMAPs in each row of panels show scRNA-seq dataset of\neach mouse, where cells are colored by (left to right): location isolated, expression-based classification,\nSNV-based classification, or consensus classification.\n4\n1-2: The observation that glioma cells can move along basement membranes, including vascular\nbasement membranes and neurons, has been known for a long time, but that this correlates with a\nspecific cell state is less well understood. There exist publications demonstrating that\nmesenchymal GBM cells prefer a perivascular invasive mode. The analysis overall is quite\ninformative determining CLAT-CAI and BV-CAI invasion signatures. It is still not evident though\nhow and why gliomas cells select for a specific invasion program if all the different tumor subtypes\nare present in the primary tumor. It is also less clear why perivascular invasion is only local and\nnot also contralateral given that the brain is highly vascularized.\nWe agree that these are interesting questions to explore further. Regarding how and why gliomas\nselect for a specific invasion program, since a single primary tumor could have cells migrating\nalong multiple routes, each tumor does not necessarily correspond to a single mode of invasion.\nNevertheless, as we and others have described previously, primary tumors typically have a bias\ntowards one cell state4,5; this bias tends to be even stronger in glioma models than in patient\ntumors6. This bias allows samples to be classified into subtypes according to their dominant state.\nConsequently, if certain cell states are predisposed to follow specific invasion routes, as our\nfindings indicate, it is likely that one invasion route would appear dominant regionally or even at\nthe level of the entire tumor.\nThe mechanisms underlying the preference of certain cell states for specific migratory routes are\nbeyond the scope of the current study. However, our new results add some insights into these\nprocesses. We performed multiplexed immunofluorescence imaging of an invasive PDX model\n(L0627) in biological replicates using the CODEX platform. Our results confirmed the presence\nof cancer cells exhibiting local invasion along blood vessels (CD31) upregulating AC-OPC\nmarkers (i.e., EGFR, OLIG2) as well as ECM-remodeling glycoproteins (e.g., TNC) (Fig. R2A-\nB, and Fig. 4C-D in the revised manuscript), suggesting that locally invading glioma (LiG) cells\nmay rely on ECM interactions for peri-vascular invasion. Supporting this, we previously found\nsignificant upregulation of ECM genes (e.g., FN1, COL11A1, COL5A2, COL6A1, COL1A1,\nCOL1A2, LAMB1) and ECM-remodeling enzymes (MMP2, LOXL2) in models showing local\nperivascular invasion but lacking distal invasion (Table S2 of the manuscript).\nConsistent with a role for ECM remodeling in vascular invasion, our spatial transcriptomics\nanalysis of cancer cells colocalized with blood vessels revealed upregulation of multiple ECM-\nrelated secreted factors (e.g. ADAM19, PDGFA, MDK, and S100A10) (Fig. R2C). Furthermore,\nligand-receptor analysis of vascular-associated cancer cells with endothelial cells showed\nincreased expression of ECM-related cancer ligands (i.e., TNC, SERPINE2, TIMP1) and their\nendothelial cognate receptors (SDC3, LRP1, CD63) (Fig. R2D, Fig. S5H in the revised\nmanuscript) as well as ECM-related vascular ligands (BCAN, FN1) and their corresponding cancer\ncell receptors (EGFR, ITGAV, ITGB8) (Fig. R2E, Fig. S5H in the revised manuscript).\nRegarding the local nature of perivascular invasion, although the brain is highly vascularized and\ninvasion along blood vessels may also occur in the contralateral hemisphere (Fig. R2F and Fig.\n2B in the revised manuscript), anatomical constraints likely limit the use of perivascular routes\nfor contralateral migration, as blood vessels rarely cross the midline, especially within the\nhemispheres. In contrast, white matter tracts - especially in the corpus callosum, which provides a\n5\nmajor inter-hemispheric route - may be exploited for contralateral dissemination (Fig. R2G). This\nlikely explains why isolating cells that actively invade across the midline selects for those capable\nof migrating along axonal or non-vascular invasion routes. While the spatial DiG state likely\ncaptures cancer cells actively invading along white matter tracts, once cancer cells have traversed\nthe corpus callosum or the smaller anterior commissure, they may also invade along BVs where\nthey upregulate ECM markers (TNC, HPLN1) in addition to neurodevelopmental lineage markers\n(SOX2, OLIG2) (Fig R2F-H and Fig. S4E in the revised manuscript), consistent with the role\nof ECM remodeling in vascular invasion.\nWe hope our expanded description, along with the new analyses and multiplexed spatial imaging,\nprovides further insights into the mechanisms and logic underlying glioma invasion route\nselection.\nFigure R2: Associations with perivascular invading cells. A) Representative CODEX image of the\nipsilateral hemisphere showing glioma cells (SOX2+) with local invasion along blood vessels (CD31+) and\nupregulation of AC-OPC markers (EGFR+, OLIG2+). B) Representative CODEX image of the ipsilateral\nhemisphere showing glioma cells (SOX2+) with local invasion along blood vessels (CD31+) and deposition\nof the glioma specific ECM glycoprotein tenascin-C (TNC). C) Volcano plot shows the differentially\nexpressed genes between spots with cancer cells and vascular and spots without vascular in the visium data.\nD) Ligand-receptor analysis: dot plots show the expression of genes in spots with vascular (x-axis) and in\nspots with malignant cells associated with vascular (y-axis). Ligands in red are highly expressed in vascular\nwith corresponding receptors (blue) expressed highly in malignant cells. E) Same as D, but in the reverse\ndirection. F-H) CODEX image showing peri-vascular (CD31+HPLN1+TNC+) and peri-neuronal\n(NeuN+) infiltration in the contralateral hemisphere by OPC glioma cells (SOX2+OLIG2+) along\nthe white matter tracts (MOG+). Scale bar: 50µm (F and H) and 800µm (G).\n6\n1-3: Is it then possible to change the invasion mode by changing the cell state?\nOur findings suggest that inducing an OPC/NPC-like state in PDXs may promote distal invasion,\nwhereas inducing MES-like states may favor local perivascular invasion. However, we did not\ndirectly test this possibility, partly because it is challenging to experimentally manipulate the\nhighly plastic GBM cell states in PDX models. Even if we were to apply perturbations that do\naffect cell states, it would be difficult to distinguish the effect of the state change from the other\neffects of these perturbations.\n1-4: Would OPC/NPC-like gliomas in general be more infiltrative than mesenchymal GBM; would\nthese spread more upon recurrence?\nGBM infiltration occurs along several pre-existing brain structures with different anatomical\nreaches. We observed that NPC/OPC cell states are associated with distal invasion via the corpus\ncallosum, which we term Distally invading Glioma (DiG) in the revised manuscript (previously\ndenoted as CLAT-CAI). In contrast, MES-like cells favor local invasion, which we term Locally\ninvading Glioma (LiG) in the revised manuscript (previously denoted as BV-CAI). Thus, these\nresults, along with findings from other studies7–9, suggest an adaptive model in which both an\noriginal GBM cell state and its favored brain structure collectively generate an infiltrative\nphenotype, and the anatomical route may determine the distance for spread. This contrasts with a\nmodel in which a single GBM cell state is associated with intrinsic capacity for invasion along\nmultiple routes.\nWe think that infiltration and recurrence are related but different biological processes, in which\nthe features of recurrent tumors are influenced not only by the infiltration of tumor cells, but also\nby treatment history and the location of tumor cells that recur. Additionally, recurrent GBM tumor\nsamples are often collected from the tumor core, within the region of contrast enhancement.\nIn our recent publications, we investigated differences between primary and recurrent tumor pairs\nby leveraging a single-cell dataset of 121 primary and recurrent GBM samples from 59 patients\n(i.e. one pre-treatment primary sample and one post-treatment recurrent sample for the majority of\npatients10,11). We found that differences in the malignant cell state distributions across time-points\nwere highly patient-specific, i.e. there was no consistent enrichment across patients of either the\nOPC/NPC state or MES state in the matched primary or recurrent tumors.\nFurthermore, this dataset includes 4 patients who had a distant recurrence (i.e., in a brain area\ndifferent than the location of the primary tumor). We performed an additional analysis using this\ndataset and found that even when comparing these 4 patients to the patients with local recurrence\nwe did not detect any major shifts or changes in the distributions of the malignant cell states\n(Figure R3).\nand recurrence location (local or distant).\n1-5 Overall, the authors demonstrate that the gene signatures are differently enriched in the\ndistinct invasion programs but the specific mechanisms that lead to and cause these invasion\nmodes still need to be determined. This will be instrumental in blocking these invasive programs.\nWe fully agree that uncovering the specific mechanisms driving distinct modes of invasion will be\nessential for developing strategies to block glioma invasion. Our study provides a foundation for\npredicting regulators of invasion by establishing robust and reproducible associations between\ntranscriptional cell states and distinct invasion routes. We have included an additional analysis to\nfurther strengthen these aspects.\nWe examined the correlation between expression of transcription factors (TFs) and our defined\ninvasion programs - DiG in the distally invading models, (GS143, L0627, L0104) (Fig. R4A), and\nLiG in locally invading models, (MGG123, MGG65) (Fig. R4B). Notably, ZEB1 emerged as the\nTF most strongly correlated with the DiG program, consistent with its recent identification as a\nkey regulator of glioma invasion8. Importantly, ZEB1 is also a part of our DiG invasion program,\nfurther validating the relevance of our transcriptional programs and highlighting ZEB1 as a\npotential therapeutic target in glioma invasion. Additionally, we identified several candidate TFs\nimplicated in regulation of LiG, including ZEB2, a known master regulator of EMT and invasion.\naverage expression, defined as log2(mean(CPM)+1) (x-axis), against the correlation of its expression with\neither A) the DiG program or B) the LiG program (y-axis).\nThis manuscript by Chanoch-Myers et al describes experiments aimed at understand the gene\nexpression patterns associated with two distinct types of invasion of gliomas in vivo. They use a\nseries of 20 human PDX/organoid models with different propensities to invaded into the mouse\nbrain. Through a tour de force of technologies they have isolated tumor cells that invade from the\ncentral bulk of the tumor across to the other hemisphere either along the perineuronal pattern, or\nadjacent to the blood vessels in a perivascular pattern. Using a combination of bulk and spatial\nsingle cell sequencing, they have identified expression for the cells that invade in both of these\npatterns. The cells invading along neurons show a oligo progenitor expression pattern, while the\nones that invade in a perivascular manner show a mesenchymal expression pattern.\nIt is certainly true that brain invasion is one of the reasons that gliomas are so difficult to treat,\ncertainly surgically. The results are also very believable, and similar to those of other glioma\nmodeling systems. Not that it would be necessary, but one would presume that these two types of\ninvasion would be seen in mouse, rat and swine glioma models and likely have largely overlapping\ngene expression signatures.\nI guess what might help the manuscript would be more insight of how this connects with the\nbehavior of the human disease.\nFor example:\n2-1. For those 20 lines that were used and had different ratios of MES- vs OPC-like invasion, what\ndid the MRI scans of the patients look like with respect to the actual invasion pattern of the glioma?\nWe agree with the reviewer that it would be interesting to evaluate the MRI scans and try to make\nvisual associations between the tumors with more MES-like or OPC-like states and their associated\ninvasion patterns. However, we do not have access to MRI scans of these patients. Instead, we\nhave added an analysis examining the enrichment of GBM cell states, as well as our invasion\nprograms, in clinical cases where GBMs recurred at a site distant from the initial resection cavity.\nSpecifically, in our recent study of single-nucleus RNA-seq (snRNA-seq) profiling 121 primary\nand recurrent GBM samples from 59 patients10,11, we found and analyzed 4 samples with a distal\nrecurrence. However, we performed an additional analysis using this dataset and found that we did\nnot detect any major shifts or changes in the distributions of the malignant cell states between\ndistant and local recurrence (Figure R3).\n2-2. Or since the ratio of OPC vs MES is well studied with single cell sequencing, are there\ncorrelations between MRI scan patterns of invasion and the ratio of OPC vs MES cells in the\ntumor.\nWe thank the reviewer’s comment, and please see our responses in Comment #2-1.\n2-3. Also, is there any thought around whether either type of invasion is more or less responsive\nto therapy? XRT for example. If that were the case, it might suggest that post radiation, the\ninvasion type might switch to a more MES types if the OPC invading cells were more likely to die.\nconnection to invasion phenotypes. While our current study does not directly assess the effect of\nXRT on specific invasive programs, in a recent publication from our lab we addressed this question\nby leveraging a single-cell dataset of 121 primary and recurrent GBM samples from 59 patients\n(i.e. one treatment unexposed primary sample and one treatment exposed recurrent sample for the\nmajority of patients10,11).\nIn this study we showed that differences in the malignant cell state distribution across time-points\nare highly patient-specific, i.e. there is no post-treatment enrichment or depletion of specific states\nthat is consistent across patients (at least as measured at the time of the second resection).\nFurthermore, we performed an additional analysis using this dataset and computed the expression\nlevel of the Distally invading Glioma (DiG; previously CLAT-CAI) and Locally invading Glioma\n(LiG; previously BV-CAI) programs in the malignant compartment in each sample and assigned\nthe samples based on the expression levels to DiG-high, LiG-high or intermediate. We did not find\nsignificant differences between the groups in the time interval between resections (a proxy for\nresponse to treatment) or in overall survival in a multivariate Cox regression analysis adjusted for\nage, gender and MGMT methylation status (Figure R5A-B).\nFigure R5 - Kaplan Meier curves for time between resections (A) and overall survival (B).\nneuroscience\nChanoch-Myers, Hara, Greenwald et al present a detailed transcriptomic characterization of\nglioblastoma invasion. The authors demonstrate differential contralateral and ipsilateral invasive\npotential of various patient-derived gliomaspheres, and characterize the molecular identities of\nmodels primed for metastasis using single cell RNA-Seq. They find that distal invasion potential is\nassociated with OPC-like cells and gene expression, while local invasion is enriched for MES-like\ncells, but interestingly downregulates MES-like genes. The authors further confirm these results\nusing pixel-based spatially-resolved RNA-Seq and delineate changes in cellular composition of\nthe neighboring tumor TME.\nThe underlying question in this study is an important and intriguing one, and the methodologies\nemployed were impressive. In particular, I found the use of a gliomasphere screen to identify\nmodels primed for distant invasion to be innovative. It provides a useful resource to the community\nand identified novel expression programs which govern peri-axonal vs. peri-vascular invasion\nroutes. However, many of the conclusions were based on shallowly tested hypotheses, and several\nclaims were largely speculative and overdrawn. In general, the manuscript would also benefit\nfrom a more coherent organization and writing. I include specific comments below.\nNote: The authors provided neither line numbers nor page numbers. I therefore quote the text and\npage number as it appears to me, where applicable. Please include this critical information in\nfuture revisions.\nWe apologize for the inconvenience and have updated the text to include page numbers.\n<Major Concerns:>\n3-1. One of the main novelties the authors stress is the decoupling of transcriptional features\nassociated with invasion potential (i.e. basal features of the primary tumor) vs. cells actively\nundergoing invasion. However, this conclusion is simply untrue and not tested in the current set\nof experiments, at least in the case of contralateral invasion. I see two major points that need\nto be addressed.\n3-1-1. There are several interesting spatial-temporal dimensions of GBM invasion from a\nbiological perspective. The primary tumor presumably first establishes itself at the injection site,\nwhich constitutes the basal cellular composition, as the authors suggest. In a second phase,\nsubpopulations of GBM migrates to distant regions on the contralateral or ipsilateral side, while\nthe remaining tumor corresponds to an altered composition. Lastly, invaded cells proliferate (as\nthe authors identify) and may take on additional states. All the authors did was to compare gene\n12\nexpression signatures between the primary tumor site with metastases at a timepoint after invasion\nhas already occurred. This in no way represents cells actively undergoing invasion to the\ncontralateral side. The ipsilateral BV-CAI cells were actively undergoing invasion which was\ninteresting to see. But this made comparisons more confusing as CLAT-CAI and BV-CAI were in\nfact two completely different populations in two stages of invasion.\nWe understand the reviewer’s concern that the CLAT-CAI program, as originally defined through\ndifferential expression between malignant cells at the contralateral and primary tumor site, may\nnot directly capture an active invasion program, but rather the end result of invasion. We have\ntherefore toned down our claim of active invasion, and renamed this program as Distally invading\nGlioma (DiG). Nevertheless, we believe that our work provides multiple lines of evidence to\nsupport the claim that DiG does represent active invasion, and we have added additional analysis\nto address this concern.\nFirst, our spatial transcriptomics analysis demonstrated that this program was consistently\nupregulated by cells across multiple time points and locations during the invasion process,\nextending from the midline and through the contralateral hemisphere. We have added additional\nanalysis to demonstrate this more clearly. We segregated the spatial map of each PDX sample into\nfour regions: bulk tumor, bulk tumor periphery, medial, and contralateral (Fig. R6A, Fig. 3D in\nthe revised manuscript, Methods). We found that the DiG and cell cycle programs are increasingly\nupregulated by malignant cells as they are distancing from the bulk tumor (Fig. R6B, Fig. 3D in\nthe revised manuscript), supporting an active invasion state rather than a program expressed at the\nculmination of invasion. In parallel, DiG downregulated genes gradually decrease in expression\nalong this spatial axis.\nIn addition, we examined the expression and spatial distribution of the DiG program at different\nstages of GBM invasion within the same model using an external Visium PDX dataset in which\nthe same model was profiled at multiple time points during the invasion process (Fig R6B, Fig.\n3D in the revised manuscript). We found that the DiG program was highly expressed in the medial\nregion and contralateral hemisphere during the early stages of invasion. In contrast, at a later stage,\nthe expression of DiG was subsequently lower within the medial region and mostly confined to\nthe contralateral hemisphere, suggesting that cells that have already undergone invasion begin to\ndownregulate it.\nFurthermore, we performed multiplexed imaging by CODEX of an invasive PDX model (L0627)\nto visualize the invasion process across spatial locations at single-cell resolution. Cells of the bulk\ntumor upregulate AC and MES state in the spatial transcriptomics analysis (Fig. R6C), and\nconsistent with this, they also express AC/MES-like markers (PDPN) by CODEX (Fig. R6D).\nHowever, while shifting patterns across spatial locations may imply a temporal axis, it does not\nprove one. To this end, we performed trajectory analysis as later suggested in Reviewer 3\nvessel invasion, while others favor distal invasion, these invasion patterns are, of course, not\nmutually exclusive. Therefore, we agree with the reviewer that it is possible that the DiG and LiG\nprograms represent not only different spatial modes of invasion but also different time points of\nthe invasion process. Ongoing work in our labs to comprehensively characterize patterns of\ninvasion in patient samples spatially and temporally will enable us to further explore these\nquestions.\nFigure R6: DiG is an active invasion program. A) Spatial map of sample BT134 at early invasion time\npoint, segregated into 4 regions: bulk, periphery bulk, medial, and contralateral. B) Top: Spatial map\ncolored by segregated zones in A, with the malignant spots colored by their scores for DiG upregulated or\ndownregulated, and cell cycle program. Bottom: Distribution of scores for the programs, left to right: DiG\ndownregulated genes, DiG upregulated genes, and cell cycle genes at each of the 4 regions. C) Spatial map\nwith distribution of Neftel states corresponding to each model in B. D) Representative CODEX image of\nipsilateral hemisphere with AC/MES-like state (PDPN+) enriched at bulk tumor.\n14\n3-1-2. Related to this, it is unclear what was actually sequenced in this set of experiments. The\nresults section note that dissections were separated into the primary injection site and invasion\nfractions. But from the methods section it seems that the two hemispheres were simply separated\nwith the midline removed. This needs to be clarified. In either case, it seems unlikely that the cells\nactively undergoing invasion would even have been sequenced. In the case of the latter, how would\nprimary tumor cells be differentiated from ipsilaterally migrated cells?\nWe agree that the previous version of the text lacked sufficient detail regarding what was profiled\nand we have modified the Methods and Results sections to improve clarity. Specifically, to isolate\ntumor cells from the multiplexed experiments, we used a stereomicroscope that allows observation\nof GFP+ tumor cells at single-cell resolution throughout the whole brain. The primary tumors were\ndissected from the ipsilateral hemispheres, centered around the injection site, identified by the hole\ncreated during stereotaxic injection, and based on GFP intensity, ensuring the removal of GFP-\nnegative outer surface of the brains. Contralateral samples were collected from the hemisphere\nopposite the injection site, avoiding the midline border regions (approximately 500 μm). During\ndissections of the contralateral hemispheres, we observed lower and restricted GFP signals when\ncompared with the ipsilateral hemisphere; GFP+ tumor cells were diffusely present in the brain\ntissue, primarily at the single-cell level, and did not extend to the entire brain regions. This\nindicates that malignant cells in the contralateral hemispheres are still in the process of active\ninvasion. In contrast, the ipsilateral samples displayed a dense GFP signal from the tissue pieces.\nBoth ipsilateral and contralateral samples were then dissociated and underwent sorting for GFP+\nmalignant cells, followed by scRNA-seq. Importantly, we found that a subset of the tumor cells in\nthe ipsilateral hemisphere upregulated DiG (Fig 2G in the original manuscript), and these\nipsilateral DiG expressing cells were found in the periphery of the tumors (Fig 4C in the original\nmanuscript). To isolate malignant cells along blood vessels, GFP+ tubular structures observed at\nthe periphery of primary tumors were manually dissected using stereomicroscopy. These vascular\ninvasion structures were not observed clearly in the multiplexed experiments, likely because highly\ninvasive tumor models take over the tumor-normal border.\nAlthough this physical dissection approach was not initially intended to capture the temporal\ndynamics of invasion, the defined invasion programs and generated single-cell data appear to\nresolve such information (see trajectory analysis in Reviewer 3 Comment -1-3), which was\nfurther supported by our spatial profiling, as described in Reviewer 3 Comment -1-1.\n3-1-3. The authors should perform pseudotime trajectory analysis to identify stages of invasion\nand its associated molecular features, and if there are in fact multiple lineages of cells during this\nprocess.\nWe performed a trajectory analysis with the three invasive models (GS143, L0627, L0104) from\nboth PDX multiplex pools. We inferred a trajectory of the cells for each model using the slingshot\nmethod12. As expected, the inferred trajectory began with cells isolated from the primary tumor\n15\ncells and ended with those isolated from the contralateral hemisphere (Fig. R7A, Fig. S3I in the\nrevised manuscript).\nIn all three models, we identified three subsequent transcriptional programs associated with\npseudotime progression that correlated with our defined invasion programs: Time-1 (early\npseudotime), which was highly correlated with DiG downregulated genes, and Time-2\n(intermediate) and Time-3 (late pseudotime) which both highly correlated with DiG upregulated\ngenes. Time-3 also showed the strongest correlation with cell cycle genes (Fig. R7A-C, Fig. S3J-\nK in the revised manuscript). Overall, the consistency of pseudotime signatures with our previous\ninvasion signatures demonstrates that our signatures were robust and could be detected by an\northogonal approach.\nNext, we examined how the pseudotime programs correspond to spatial location in the GS143\nmodel, for which we had high-resolution spatial transcriptomics data (Visium HD). We defined\nthe malignant cell populations by unsupervised spatial clustering (Singhal et al. 2024). We found\nthree malignant programs that were correlated with the pseudotime programs (Early timepoint,\nR=0.4; Intermediate timepoint, R=0.72; Late timepoint, R=0.56). Spatial mapping of these\nprograms generated a spatial continuum from Time 1-3 (Fig. R7B, Fig. S3K in the revised\nmanuscript). Taken together, the pseudotime and spatial analysis support a temporal and spatial\ntrajectory progressing from the tumor core to proneural states (DiG) to cell cycle.\nWhile trajectory analysis can suggest potential progression of cellular states, it is an inference\nmethod, and, therefore, doesn’t prove cell transitions between states. To examine if transitions\noccur within a lineage, we employed lentiviral genetic barcoding to track clonal relationships\nbetween primary and invasive cells. This approach ensures that observed transitions in pseudotime\nreflect actual within-lineage invasion events, rather than being technical artifacts or unrelated\ntranscriptional variation. We injected GS143 cells with a barcoding lentivirus into mice (n=2) and\nisolated single cells from both the primary tumor and contralateral invasion sites. Cells were\nprofiled by single-cell RNA-seq using SS2, with 484 cells that passed QC.\nWe detected 14 unique barcodes that were found in at least 3 cells, of which 5 barcodes had\nrepresentation in both the primary tumor and contralateral invasion site (Fig. R7D). We used the\nsame method as described above to infer a trajectory in this sample. First, we found that cells from\nthe primary tumor had lower pseudotimes than those in the contralateral invasion sites (Fig. R7E,\np<2.2e-16). This was also true for individual clones, but the difference was not statistically\nsignificant due to the low number of cells. We found similar time-associated patterns in the\nbarcoding experiment. Importantly, cells from individual clones mapped along the inferred\ntrajectory (Fig. R7F), demonstrating that cells from the same lineage transition through distinct\ninvasion-associated transcriptional states.\nCollectively, these results support the notion that DiG captures a program dynamically linked to\ninvasion and proliferation, rather than simply reflecting a static difference between primary and\ninvasive sites. By integrating barcoding with trajectory analysis, we provide a direct lineage-based\nvalidation of invasion-associated transcriptional states, strengthening our conclusions.\n(GS143, L0104, L0627), separated by mouse (two rows). Plots are colored by Location (Primary, Invasion),\nPseudotime, Time-associated program (early, intermediate, late) assignment. B) Spatial map of GS143.\nMalignant cells are colored by their assignment to time-associated programs. C) Correlation plot shows the\ncorrelation of expression of each time associated program generated per mouse, as well as the invasion\nprograms and cell cycle programs. D) Table shows each barcode and the number of cells that it was detected\nin. E) Boxplot shows the pseudotimes of cells in primary or invasion sites, separated by cells with the same\nbarcode. F) Scatterplots show the pseudotime (x-axis) against the score for time associated programs (y-\naxis). Each row colors cells from an individual barcode\n17\n3-2 The authors need to perform single molecule FISH (i.e. RNA-Scope or HCR) to confirm the\nkey findings from spatial transcriptomics with single cell resolution. The authors should confirm\na selection of disrupted genes and their associated cell-types. In particular, it would be interesting\nto see if MES-like genes are downregulated in MES-like cells or other fractions.\nWe thank the reviewer for this important suggestion and agree that single-cell resolution spatial\ndata would strengthen our claims. We used two different approaches to address this. First, we\nperformed spatial transcriptomics with Visium HD on an invasive PDX sample (GS143). This\nhigh-resolution platform enables near-single-cell spatial profiling. However, a significant\nlimitation of this approach is the use of a human probe set, enabling us to only robustly capture\nthe malignant cells. Using this approach, we could clearly identify invasive tumor cells expressing\nthe DiG signature, thus validating the presence of this program at single-cell resolution (Fig R8A).\nSecond, as an orthogonal approach, we performed CODEX multiplexed imaging on an additional\ninvasive model (L0627), employing a 12-marker antibody panel that spans the major cell states\nand cell types as well as ECM markers. CODEX enabled us to validate several major claims from\nthe single-cell and spatial transcriptomics analysis at single-cell resolution: (1) Local invasion of\ncancer cells (EGFR, OLIG2) with upregulation of ECM markers (HPLN1, TNC) along blood\nvessels (CD31) (Fig. R2B, and Fig. 4D in the updated manuscript) (2) Distant invasion of cancer\ncells with pro-neural features (SOX2, OLIG2) along white matter tracts/axons (MOG, NeuN) (Fig.\nR2F-H and Fig. 2B in the updated manuscript) (3) Recruitment of myeloid cells (CD45,F4/80,\nP2Y12) to invading cancer cells along white matter tracts (MOG) (Fig. R8B and Fig. 5B in the\nupdated manuscript). These findings strengthen the specificity of our invasion programs and\nhighlight the dynamic reprogramming of cell states during invasion.\nFigure R8: Invasion programs at single-cell resolution with Visium HD/CODEX. A) Spatial\nmap of Visium HD data for GS143. Malignant spots colored by DiG upregulated or downregulated\ngenes. Nonmalignant spots colored in green. Mixed (malignant and nonmalignant) spots colored in pink. B)\nCODEX image showing myeloid recruitment (F4/80+ P2y12+ CD45+) to the white matter tracts\nalongside invading SOX2+ cancer cells in the contralateral hemisphere in PDX model L627.\nNeurons (NEuN) are in green.\n18\n3-3. The authors interestingly propose peri-axonal vs. peri-vascular interactions (suggested by\nprevious literature), but do not go on to characterize or confirm this. The last section delves into\nthe tumor TME, but in alterations in immune cell compositions post-invasion. It would be\ninteresting to perform cell-cell interaction inference to probe key players involved in the process\nof invasion, both from the side of GBM cells and host cells. This would be more relevant to the\ncurrent study in my opinion.\nWe thank the reviewer for this suggestion and have performed additional experiments and analysis\nto examine interactions between invading cells and their TME, with a focus on\nmacrophage/microglia (henceforth referred to as macrophages), endothelial cells, and neurons.\nFirst, we incorporated a niche-labeling system into our invasive GBM model to dissect the\ncomposition of non-malignant cells within the invasive niche at single-cell resolution. Specifically,\nGFP+ GBM cells were engineered to express a cell-penetrating mCherry reporter13, enabling its\nuptake by adjacent cells. This allowed for isolation of GFP⁻mCherry⁺ nonmalignant TME cells\nand GFP⁺ malignant cells by FACS (Fig. R9A-B), which were then profiled by Smart-seq2\nscRNA-seq. Among the nonmalignant, niche-associated populations, we identified macrophages,\nastrocytes, OPCs, and oligodendrocytes (Fig. R9B).\nWe focused on macrophage–tumor interactions by inferring ligand-receptor interactions. This\nanalysis revealed several interactions, with invasive macrophages upregulating ligands (e.g., MIF,\nTGFBI, TGM2, SERPINF1), paired with their corresponding receptors on GBM cells. Conversely,\ninvasive GBM cells expressed ligands (e.g., MDK, ANGPTL1, LGALS1, MIF, PDGFA),\ncorresponding to cognate receptors on macrophages (Fig. R9C-D). For validation, we inferred\nligand-receptor interactions between colocalized invasive macrophages and DiG cancer cells in\nthe spatial transcriptomics data, with MDK-LRP1 (Fig. R9E) and MIF-CD44/LRP1(Fig. R9F)\nemerging as the most robust interactions. Upregulation of MDK by invading cancer cells and its\nreceptor LRP1 on macrophages is consistent with recent work showing that GBM cell-secreted\nMDK mobilizes macrophages/microglia at the invasion front by upregulating plexin-B214. In both\nthe niche labeling experiment and the spatial transcriptomics analysis, MIF was expressed by both\ninvasive macrophages and GBM cells, suggesting potential autocrine and paracrine signaling. The\nenriched expression of MIF in the invasive niche supports its possible role as a mediator of TME–\ntumor crosstalk in facilitating invasion.\nNext, we used the niche-profiling experiment to validate our previous finding from the spatial\ntranscriptomics and found that upregulated genes by contralateral macrophages were enriched for\nOXPHOS-related genes (Fig. R9G). Interestingly, a greater enrichment in the contralateral\nmacrophages was observed when comparing them to those from normal brain, suggesting that\nOXPHOS upregulation is linked to their interaction with invading malignant cells.\nSince the niche-labeling experiment did not robustly capture endothelial cells and neurons, cell\ntypes relevant to the invasion phenotypes, we leveraged the spatial transcriptomics data of the\ninvasion models as well as scRNA-seq (our PDX models) and snRNA-seq datasets10,11, to\ncharacterize these populations.\nsecreted factors (e.g., ADAM19, PDGFA, MDK, and S100A10), suggesting ECM remodeling by\nLiG cancer cells (Fig. R2C).\nSecond, spatially-dictated ligand-receptor analysis of vascular-associated cancer cells and\nendothelial cells showed upregulation of ECM-related cancer ligands (i.e., TNC, SERPINE2,\nTIMP1) and their endothelial cognate receptors as well as ECM-related vascular ligands (BCAN,\nFN1) and cognate cancer cell receptors (EGFR, ITGAV, ITGB8) (Fig. R2D-E, Fig. S5H in the\nrevised manuscript). Furthermore, ligand-receptor interactions inferred between LiG cancer cells\n(PDX models MGG123, MGG65) and endothelial cells were enriched for ECM-related markers\n(e.g., ligands ANXA1, FBLN1, and MDK on LiG cancer cells and VEGFC, FN1, COL4A1, and\nADAM10 on endothelial cells (Fig. R10A-B, Fig. S5I in the revised manuscript). To corroborate\nthese findings at the protein level, we used CODEX imaging, confirming that locally invasive and\ncontralateral vessel-associated tumor cells upregulate ECM proteins such as TNC and HAPLN1\n(Fig. R2A-B). Collectively, these findings suggest ECM remodeling as a key mechanism of\nvascular invasion, driven by. By both tumor and adjacent vasculature.\nLastly, we inferred ligand-receptor interactions between neurons and DiG cancer cells in invasive\nmodels (GS143, L104, L627). We identified ligands expressed by neurons and their cognate\nreceptors on DiG cancer cells that may provide guidance cues for invading cancer cells (e.g.,\nSLIT3-ROBO1) (Fig. R10C, Fig. S5J in the revised manuscript). We also identified potential\nsynaptic interactions between DiG cancer cells and neurons, suggesting neuronal activation by\ninvading cancer cells (e.g., NLGN4Y-NRXN1). (Fig. R10D, Fig. S5J in the revised manuscript).\nprofiling experiment: GFP pos GBM cells secrete mCherry dye which can be uptaken by neighboring GFP-\nneg nonmalignant cells. B) UMAPS of niche experiment: cells colored by organism origin, GFP/mCherry,\nLocation or identified cell type. C) Ligand receptor analysis. Ligands expressed by Macrophages (Top) and\ntheir corresponding receptors in Malignant cells (Bottom). Colored by expression in normal mouse brain,\nInvasion, or Primary tumor. D) Same as C) except Ligands expressed by Malignant cells and Receptors\nexpressed by macrophages. E) Ligands expressed by Inflammatory Macrophages (y-axis) and their\ncorresponding receptors in Inflammatory-macrophage associated malignant cells (x-axis). F) Same as E)\nexcept Ligands expressed by Inflammatory-macrophage associated malignant cells and Receptors\nexpressed by Inflammatory-macrophages. G) Enrichment of OXPHOS genes by contralateral macrophages\nvs ipsilateral macrophages (top) or contralateral macrophages and normal mouse brain macrophages\n(bottom).\n21\nreceptors in blood-vessel invading cells (bottom). B) Ligands expressed in blood-vessel invading cells (top)\nwith receptors expressed in endothelial cells from (REF) (bottom). C) Ligand Receptor analysis: Ligands\nexpressed in neurons from (REF) (top) with corresponding receptors in contralateral invading cells\n(bottom). D) Ligands expressed in contralateral invading cells (top) with receptors expressed in neurons\nfrom (REF) (bottom).\n3-4. Additionally, previous studies identified neuronal activity to be a critical source for driving\nGBM invasion. Do the authors observe similar synaptic alterations in the tumour TME, and are\nthere changes in GBM composition at different stages of invasion? This would be another\ninteresting question that could be answered with a pseudotime analysis but which is completely\nmissing from the current study.\nWe thank the reviewer for this insightful comment. Indeed, prior studies, including a recently\npublished manuscript8, have shown that neuronal activity and neuron–glioma interactions play an\nimportant role in GBM invasion. Notably, the DiG invasion program includes many OPC/NPC-\nrelated genes, making it challenging to disentangle whether this shift reflects engagement with\nneuronal circuits (e.g., through synaptic interactions) or simply a transition to a progenitor-like\nstate.\nWhile we do not directly assess synaptic gene expression or functional neuronal signaling, the\nemergence of NPC/OPC-like genes in the DiG program, as well as the pseudotime programs\n(Comment 3-1-3), supports the idea that GBM invasion may involve a dynamic reprogramming\nprocess, potentially influenced by neuronal cues.\nIn addition, we have performed additional analysis of neuronal interactions, as described in\nresponse to the previous comment (Comment 3-3), revealing upregulation of several markers\nsuggestive of synaptic activation (e.g, NLGN4Y, CADM2, ALCAM). We have now clarified these\nfindings and their interpretation in the revised manuscript (Fig. S5K), acknowledging both the\nstrengths and the limitations of our current data in addressing neuron–glioma interactions.\n<Minor Concerns:>\n3-5 <invasive cell %>\nOn page 5: \"Invasive cells...primarily derived from three highly invasive GBM models (L0627,\nL0104, GS143), with 54-78% of cells from each model being invasive\". It is unclear where these\nnumbers are coming from. The highest range noted here is 78% but L0104 in Figure 1E,H is\nclearly over 80% in the [% invading cells] axis. Additionally, these numbers do not seem to\ncorrelate to the number of primary tumor vs. invasion cells presented in Figure S1E, where\nL0104 seem to have almost entirely invasion cells, while L0627 also appear to have a lot more\nthan 50% invading cells. Please clarify where the discrepancies are coming from, and if it is\nfrom overlapping datapoints in the UMAP, this needs to be corrected in all figures.\nWe appreciate the reviewer pointing this out. Regarding the values of invasion in the text, they\nwere previously incorrect. We have updated the text with the correct values for average invasion\nper cell line (Fig R11A). Regarding L0627, we agree that the current coloring may give the\nimpression that more than 50% of the cells were invading, but this is because of overlapping points\nas the reviewer pointed out. We have now updated this UMAP to fix this issue (Fig R11B, and\nFig S1D in the updated manuscript).\nFigure R11: Percent invading cells per model. A) Table shows the percent cells found in the\ncontralateral hemisphere for each model. B) UMAPS fixed to show random sampling of cells\n(right) compared to previous version which plotted the invasive cells on top (left).\n3-6. Please comb through the figures to include all necessary information for interpretation. For\nexample, Figure S1B is missing cluster labels which is critical for interpretation and Figure 2C is\nmissing scale bars.\nWe appreciate the reviewer pointing this out. We have added missing labels and scale bars\naccordingly.\n3-7. Page 8: \"Thus, the distinction between primary tumor cells and invasive cells lies in the\nfrequency of the CLAT-CAI signature rather than a coherent difference across all cells\" - it is\nunclear what this sentence means.\nIn this sentence, we aimed to explain that the difference between primary tumor cells and invasive\ncells was in the proportion of cells that express the CLAT-CAI program and not that this program\nis upregulated in all of the invasive cells and downregulated in all of the primary cells. We have\nremoved this section of the text because we address the same point again later with an additional\nspatial transcriptomics analysis (see Reviewer 3 Comment 1-1).\n3-8. Page 9: \"...our analyses show multiple associations of cellular states (i.e., expression\nsignatures) with invasion-related phenotypes...(Fig 3C). How does Figure 3C show multiple\ncellular states?\nthat LiG is associated with AC/OPC states. We removed the reference to the figure since it is a\nclaim we make from the earlier sections.\n3-9. If possible, the authors should provide higher resolution confocal images of vessels and the\ncorpus callosum during invasion, instead of the blurry and over-saturated images in Figure 2C\nwhere it is difficult to tell what one is looking at.\nWe strongly agree with this comment. We have now included confocal scanning images of GBM\nmodels showing local perivascular invasion (Fig. R12, Fig. 2D in the revised manuscript). From\nsuch areas, we isolated GFP+ cells using microdissection and dissociated them for isolation and\nscRNA-seq profiling. As addressed in Reviewer 3 Comment 2 and Fig. R8, we have now added\nCODEX images of local (LiG) invasion and distal (DiG) invasion as well as macrophage\nrecruitment to DiG cancer cells along white matter tracts.\nFigure R12: Confocal imaging of perivascular invasion. Representative confocal imaging of GFP+\ntumor cells invading along blood vessels (CD31+).\n3-10. There are some inaccuracies in terminology. Please correct spatial genomics to spatial\ntranscriptomics as the manuscript does not utilize the former. Please also provide the spot sizes\nclearly in the results and figure legends, and include the corresponding H&E counterstain images\nin the supplementary figures.\nWe have updated the text accordingly.\n1. Kinker, G. S. et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular\nheterogeneity. Nat. Genet. 52, 1208–1218 (2020).\n2. Heaton, H. et al. Souporcell: robust clustering of single-cell RNA-seq data by genotype\nwithout reference genotypes. Nat. Methods 17, 615–620 (2020).\n3. Neavin, D. et al. Demuxafy: improvement in droplet assignment by integrating multiple\nsingle-cell demultiplexing and doublet detection methods. Genome Biol. 25, 94 (2024).\n4. Neftel, C. et al. An Integrative Model of Cellular States, Plasticity, and Genetics for\nGlioblastoma. Cell 178, 835-849.e21 (2019).\n5. Couturier, C. P. et al. Single-cell RNA-seq reveals that glioblastoma recapitulates a normal\nneurodevelopmental hierarchy. Nat. Commun. 11, 3406 (2020).\n6. Hara, T. et al. Interactions between cancer cells and immune cells drive transitions to\nmesenchymal-like states in glioblastoma. Cancer Cell 39, 779-792.e11 (2021).\n7. Venkataramani, V. et al. Glioblastoma hijacks neuronal mechanisms for brain invasion. Cell\n185, 2899-2917.e31 (2022).\n8. Wu, Y. et al. Neurodevelopmental hijacking of oligodendrocyte lineage programs drives\nglioblastoma infiltration. Dev. Cell (2025) doi:10.1016/j.devcel.2025.04.022.\n9. Manoharan, V. T. Spatiotemporal modeling reveals high-resolution invasion states in\nglioblastoma.\n10. Spitzer, A. et al. Deciphering the longitudinal trajectories of glioblastoma ecosystems by\nintegrative single-cell genomics. Nat. Genet. 57, 1168–1178 (2025).\n11. Nomura, M. et al. The multilayered transcriptional architecture of glioblastoma ecosystems.\nNat. Genet. 57, 1155–1167 (2025).\n12. Street, K. et al. Slingshot: cell lineage and pseudotime inference for single-cell\ntranscriptomics. BMC Genom. 19, 477 (2018).\n13. Ombrato, L. et al. Metastatic-niche labelling reveals parenchymal cells with stem features.\nNature 572, 603–608 (2019).\n14. Kang, S. et al. Glioblastoma shift from bulk to infiltrative growth is guided by plexin-B2-\nmediated microglia alignment in invasive niches. Nat. Cancer 1–19 (2025) doi:10.1038/s43018-\n025-00985-4.\n26"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors responded adequately to the reviewers’ comments and incorporated substantial new data to strengthen their findings. Overall, the inclusion of the latest data is very impressive.\nIt would be, however, important to discuss that some of the observed local and distal invasion programs may still have similarities, as they can also represent different time points of glioma invasion. An example would be that the authors list ZEB1 as a TF for distal and ZEB2 as a TF for local invasion, although both TFs drive EMT/invasion.", "reviewer_2": "there is a lot of analysis done here, and even more now after the revisions. it is unfortunate that the authors do not have access to the MRI scans for these patients since that is actually how invasion would be meaningfully measured in patients. in general, my questions were not really answered by the rebuttal. however, the other reviewer's questions seemed to be. and further work is not likely to solve the problems i raised. the paper deserves to be published in my opinion.", "reviewer_3": "The authors have significantly improved their manuscript through additional experiments and data analyses, and addressed all my concerns. In particular, I applaud their efforts in extending to additional lentiviral barcoding-based lineage tracing and mCherry niche labeling. The quality of the revision experiments are quite high."}, "author_response": "Response to Reviewers:\nReviewer #1\nComment: The authors responded adequately to the reviewers’ comments and incorporated\nsubstantial new data to strengthen their findings. Overall, the inclusion of the latest data is very\nimpressive. It would be, however, important to discuss that some of the observed local and distal\ninvasion programs may still have similarities, as they can also represent different time points of\nglioma invasion. An example would be that the authors list ZEB1 as a TF for distal and ZEB2 as\na TF for local invasion, although both TFs drive EMT/invasion.\nResponse: We thank the reviewer for this insightful comment and for recognizing the strength of\nour revised data. We agree that while we have delineated distinct transcriptional programs for\nLocal (LiG) and Distal (DiG) invasion, these cellular states may exist on a dynamic continuum\nand share fundamental mechanisms, such as EMT/invasion pathways driven by transcription\nfactors like ZEB1 and ZEB2.\nTo address this, we have added the following sentence to the Discussion section to provide this\nnecessary context:\n\"While our findings delineate distinct invasion route-specific programs (LiG vs. DiG), we\nacknowledge these states may also represent different points along a dynamic temporal\nprogression, potentially sharing core mechanisms that are adapted to the immediate anatomical\nniche.\"\nThis addition clarifies that our findings represent snapshot states optimized for specific routes\nbut do not preclude shared or transitional features between them.\nReviewer #2\nComment: there is a lot of analysis done here, and even more now after the revisions. it is\nunfortunate that the authors do not have access to the MRI scans for these patients since that is\nactually how invasion would be meaningfully measured in patients. in general, my questions were\nnot really answered by the rebuttal. however, the other reviewer's questions seemed to be. and\nfurther work is not likely to solve the problems i raised. the paper deserves to be published in\nmy opinion.\nResponse: We appreciate the reviewer's comprehensive assessment of the manuscript and their\nultimate support for publication. We acknowledge the reviewer's valuable point regarding the\nclinical relevance and the limitation of not having patient-matched, pre-operative MRI data to\ncorrelate with our in vitro and in vivo model findings. Clinically, gross invasion is indeed\nmeasured by MRI, and correlating our molecular data with this macroscopic radiological\nphenotype would be highly informative.\nReviewer #3\nComment: The authors have significantly improved their manuscript through additional\nexperiments and data analyses, and addressed all my concerns. In particular, I applaud their\nefforts in extending to additional lentiviral barcoding-based lineage tracing and mCherry niche\nlabeling. The quality of the revision experiments are quite high.\nResponse: We are very grateful for the reviewer's positive assessment and their recognition of\nthe significant effort dedicated to the revision. We are pleased that the additional experiments,\nparticularly the lentiviral barcoding-based lineage tracing and the mCherry niche labeling, have\nsuccessfully addressed all of their initial concerns and enhanced the manuscript's overall quality\nand impact."}]

Decision letter