The landscape of tumor cell states and spatial organization in h3-k27m mutant diffuse midline glioma across age and location

ABSTRACT Histone 3 lysine27-to-methionine (H3-K27M) mutations most frequently occur in diffuse midline gliomas (DMGs) of the childhood pons but are also increasingly recognized in adults.

Their potential heterogeneity at different ages and midline locations is vastly understudied. Here, through dissecting the single-cell transcriptomic, epigenomic and spatial architectures of

a comprehensive cohort of patient H3-K27M DMGs, we delineate how age and anatomical location shape glioma cell-intrinsic and -extrinsic features in light of the shared driver mutation. We

show that stem-like oligodendroglial precursor-like cells, present across all clinico-anatomical groups, display varying levels of maturation dependent on location. We reveal a previously

underappreciated relationship between mesenchymal cancer cell states and age, linked to age-dependent differences in the immune microenvironment. Further, we resolve the spatial organization

of H3-K27M DMG cell populations and identify a mitotic oligodendroglial-lineage niche. Collectively, our study provides a powerful framework for rational modeling and therapeutic

interventions. SIMILAR CONTENT BEING VIEWED BY OTHERS K27M IN CANONICAL AND NONCANONICAL H3 VARIANTS OCCURS IN DISTINCT OLIGODENDROGLIAL CELL LINEAGES IN BRAIN MIDLINE GLIOMAS Article 05

December 2022 EVOLVING CELL STATES AND ONCOGENIC DRIVERS DURING THE PROGRESSION OF IDH-MUTANT GLIOMAS Article 21 November 2024 SPATIAL TRANSCRIPTOMICS REVEALS NICHE-SPECIFIC ENRICHMENT AND

VULNERABILITIES OF RADIAL GLIAL STEM-LIKE CELLS IN MALIGNANT GLIOMAS Article Open access 23 February 2023 MAIN Diffuse midline gliomas (DMG) driven by a lysine27-to-methionine (K27M)

mutation in histone 3 (H3) are among the most lethal brain tumors1,2,3,4,5. Primarily identified in younger children (<10 years), the same oncohistone mutation is also recurrently

observed in midline gliomas in adults6,7,8. In children, the spatiotemporal pattern of H3-K27M DMG incidence, peaking at 6–9 years of age in the brainstem pontine region, has shaped the

hypothesis that the cell-intrinsic and -extrinsic context in which the K27M mutation occurs and elicits oncogenic transformation is developmental stage specific9. Indeed, previous studies

have hinted at precursor cells in the pons10 and an early neurodevelopmental window11 as spatiotemporal correlates in K27M mutation-mediated gliomagenesis. Cell-intrinsically, the K27M

mutation leads to broad epigenetic dysregulation and thus transformation of a developmentally restricted cell to a tumorigenic stem-like state12,13,14,15,16,17,18. The resulting active

chromatin landscape reflects an early oligodendroglial lineage19,20. Single-cell RNA-sequencing (scRNA-seq) of pediatric, predominantly pontine H3-K27M tumors, further demonstrated that most

glioma cells are stalled in a cancer stem cell-like oligodendrocyte precursor cell (OPC)-like state that is capable of self-renewal and tumor initiation21,22. In contrast, more

differentiated noncycling glia-like cells were shown to have lost their tumorigenic capacity21. Together, this indicates OPC-like cells to be at the core of K27M mutation-mediated

tumorigenesis, and hence, may present a strategic therapeutic target in pediatric pontine H3-K27M DMGs. However, it remains incompletely understood whether H3-K27M DMGs of different midline

locations—such as thalamus, pons or spinal cord—as well as different age groups and different morphological features at presentation, have similar cellular compositions. In particular, the

more recently recognized group of adolescent (10–19 years) and adult (≥20 years) H3-K27M DMGs remains understudied. In addition to cell-intrinsic modes of dysregulation, mounting evidence

indicates that microenvironmental factors critically contribute to glioma growth23,24,25,26,27,28, and it has been suggested that the developing brain provides a permissive environment that

can be exploited for pediatric brain tumor growth29,30. However, the interplay between age- and region-specific tissue environments and the varying clinico-anatomical characteristics of

H3-K27M DMGs, and its contribution to tumor pathology remain unexplored. To address these questions, we have utilized single-cell multi-omics and spatial transcriptomic approaches to profile

an extended cohort of H3-K27M DMGs encompassing a broad range of age groups and anatomical locations. We thereby identify how age- and location-dependent contexts underlie cell-intrinsic

and -extrinsic features that together determine variation in glioma spatial and cellular architecture in light of the common K27M mutation. RESULTS COHORT OF H3-K27M DMGS ACROSS AGE GROUPS

AND LOCATIONS We conducted multi-omic profiling of 50 H3-K27M mutant patient tumors, selected only by criteria of the oncohistone mutation, spanning pontine (_n_ = 27), thalamic (_n_ = 20),

lower brainstem (_n_ = 1) and spinal (_n_ = 2) locations (Fig. 1a,b and Supplementary Table 1). The median age was 12 (2.5–68) years, encompassing 36 pediatric (18 early childhood (0–9

years), 18 adolescent (10–19 years)) and 14 adult (20–68 years) tumors. Samples were obtained pre-treatment (_n_ = 30) and post-treatment (_n_ = 20) from 29 female and 21 male patients. We

performed deep full-length Smart-seq2 fresh single-cell (_n_ = 18) or frozen single-nucleus (_n_ = 25) RNA-sequencing (scRNA-seq/snRNA-seq) of 43 tumors (Fig. 1a–c). We additionally analyzed

the open chromatin profiles of eight tumors utilizing the single-cell/single-nucleus assay for transposase-accessible chromatin using sequencing (scATAC-seq/snATAC-seq), as well as the

single-cell spatial transcriptomic architecture of 14 tumors by in situ sequencing (Fig. 1a,b). To identify other mutations, we performed whole or targeted exome sequencing in 43 of 50

tumors (Fig. 1b). Recurrent mutations in _TP53, PDGFRA_ and _PIK3CA_ were broadly observed across all clinico-anatomical groups stratified by age and location, while alterations in

_HIST1H3B_ and _BRAF_ were only rarely detected in childhood tumors, which is in line with previous reports of H3-K27M DMGs1,5,8,31,32. Overall, our cohort covers a representative

clinico-molecular range of H3-K27M DMGs. Interestingly, we did not detect significant differences in co-mutational profiles between different groups, and next set out to investigate

non-genetic features and heterogeneity of H3-K27M DMGs across different spatiotemporal contexts. H3-K27M DMG CELL COMPOSITION ACROSS AGE AND LOCATION We aimed at delineating and comparing

transcriptional heterogeneity within our cohort stratified by age and location (Fig. 1c,d and Extended Data Fig. 1a–d). Complementary approaches assessing inter- and intratumoral

heterogeneity concordantly identified tumor cells differentially expressing actively cycling, OPC-like, ‘astrocyte-likeʼ (AC-like), ‘oligodendrocyte-likeʼ (OC-like) and ‘mesenchymal-likeʼ

(MES-like) signatures (Fig. 2a–d, Extended Data Fig. 2a–g and Supplementary Table 2). OPC-like cells were further resolved into three subpopulations (OPC-like-1, OPC-like-2 and OPC-like-3)

(Fig. 2c,d and Extended Data Fig. 2b,c). Interestingly, the MES-like signature, which has been described in glioblastoma (GBM)33,34, has not been identified in H3-K27M DMGs before, hinting

at unique properties uncovered from previously understudied clinico-anatomical groups within our extended cohort. OPC-like cells were ubiquitously present in all tumors independent of age or

location (Fig. 2a,d and Extended Data Fig. 2e,f). Interestingly, even in this expanded cohort, we did not detect any neuronal lineage tumor cells, placing this in contrast to all other

high-grade glioma types and isocitrate dehydrogenase (IDH)-mutant glioma33,35,36,37. To investigate whether this may be a phenomenon specific to the midline location, we single-cell profiled

two location- and age-matched IDH-mutant midline gliomas (Supplementary Table 1), revealing that neuronal lineage programs are present within rare midline IDH-mutant tumors (Extended Data

Fig. 2h). Hence, this comparison of primary gliomas of the same location and age groups, but different genotypes, supports a direct cell-intrinsic effect of the K27M mutation to skew tumor

cells toward a glial/OPC-like instead of a neuron-like identity. We next reconstructed networks of active transcription factors (TFs) and their downstream gene targets (gene regulatory

networks (GRNs)) (Fig. 2e and Supplementary Table 3) by single-cell regulatory network inference and analysis (SCENIC)38. We indeed found key GRNs known from normal glial specification (for

example, SOX10 in OPC-like cells, TFEB in OC-like cells, SOX9 in AC-like cells) to be likewise active in respective H3-K27M DMG tumor cell counterparts, highlighting parallels between normal

developmental and glioma cell fate determination. Moreover, we identify GRNs (for example, GLI2 and NFATC4) that have not yet been implicated in normal development and may hence present

glioma-specific regulatory aberrations. We next compared cellular compositions across tumor locations and age groups (Fig. 2f,g and Extended Data Fig. 2i,j). Interestingly, the MES-like

metaprogram was substantially enriched in adult tumors (Fig. 2f), which persisted when we controlled for location as a potential confounding factor (Extended Data Fig. 2i). This was

validated by RNA in situ hybridization (Fig. 2h). Except for one NF1_-_mutated pediatric tumor (Fig. 1b), which was associated with a stronger MES-like signature as previously reported33,39,

we did not detect any additional recurring genetic mutations in coding gene regions in tumors enriched for MES-like cells, suggesting either non-coding mutations and/or non-genetic

determinants may underlie the observed age-specificity. As such, this age-related difference points toward the emerging role of the tumor microenvironment in shaping the MES-like signature,

as has been illustrated in recent studies27,28,40. Together, we demonstrate that H3-K27M DMGs are biased toward an OPC-like cell identity independent of age or midline location, which

suggests cell-intrinsic effects of the K27M oncohistone mutation itself rather than environmental determinants to underlie this cellular state. Contrastingly, an association with age is

observed for the MES-like signature (Fig. 2i), potentially linking this cellular state to cell-extrinsic/environmental drivers. LOCATION SPECIFICITY OF OPC-LIKE SUBPOPULATIONS We next

examined the three OPC-like subpopulations uniquely detected in our extended scRNA-seq dataset, termed OPC-like-1, OPC-like-2 and OPC-like-3 (Figs. 2c,d and 3a). While all OPC-like

subpopulations were defined by high expression of canonical OPC markers (for example, _PDGFRA, SOX10_ and _OLIG1/2_), these markers together with other known marker genes of committed OPCs

(for example, _CSPG4, GPR17_ and _EPN2_) were most highly expressed by OPC-like-1 cells (Fig. 3a–c)41,42. In contrast, OPC-like-2 and −3 cells depicted higher expression of marker genes

linked to more immature oligodendrocyte precursors of the developing brain, also termed pre-OPCs—a state of oligodendroglial lineage differentiation between less differentiated neural stem

cell and more differentiated OPC (for example, _ASCL1, HES6, BTG2, DLL1_ and _EGFR_) (Fig. 3c)41,42,43. Additionally, OPC-like-2 cells highly expressed genes encoding ribosomal proteins (for

example, _RPL17_ and _RPS18_), and OPC-like-3 cells exhibited higher expression of immediate early response genes (for example, _JUNB_ and _EGR1_) (Fig. 3a), which have been previously

described as markers of different normal (pre-)OPC subpopulations44,45. When we projected these OPC-like subpopulations onto scRNA-seq atlases of the human telencephalon and mouse

cortex41,43,46, the OPC-like-1 subpopulation indeed mapped to committed/maturing OPCs, whereas OPC-like-2 and OPC-like-3 cells were more similar to pre-OPCs (Fig. 3d–f and Extended Data Fig.

3a–c). Comparison with cell populations from other glioma types and trajectory analyses (Extended Data Fig. 2c; 3d,e) also pointed toward a more immature state of OPC-like-2 and OPC-like-3

cells, and stronger lineage commitment of OPC-like-1 cells. Analysis of OPC-like subpopulation-specific GRNs using our scRNA-seq dataset identified TFs such as SHOX2 and OTX2 to be most

specifically active in OPC-like-1 cells (Fig. 3g and Supplementary Table 3). GRNs specific to OPC-like-2 cells included Notch signaling regulator HES6 as well as multiple patterning TFs of

the HOX family, and GRN characteristics of OPC-like-3 cells were linked to the AP-1 TF family (Fig. 3g). Of note, HOX patterning TFs have been demonstrated to be expressed in mice embryonal

pre-OPCs while being downregulated in postnatal OPCs44. Moreover, immediate early response regulators have been implicated as specific to human pre-OPCs compared to committed OPCs45, further

hinting at a more immature and pre-OPC-like state of DMG OPC-like-2 and OPC-like-3 cells. We next compared proportions of these OPC-like subpopulations across our spatiotemporally

stratified cohort and observed a remarkable enrichment of pre-OPC-like (OPC-like-2 and OPC-like-3) cells in pontine compared to thalamic tumors. Conversely, OPC-like-1 cells were enriched in

thalamic tumors (Fig. 3h). These differences remained when stratifying for age groups as potential confounders (Extended Data Fig. 2j). Therefore, we identify tumor location as a contextual

determinant of OPC-like states, with immature pre-OPC-like progenitors enriched in pontine, and more committed OPC-like cells enriched in thalamic tumors. THE OPEN CHROMATIN LANDSCAPE OF

H3-K27M DMG CELL POPULATIONS To resolve how H3-K27M DMG cellular heterogeneity is governed at the chromatin level, we probed single-nucleus accessible chromatin profiles by snATAC-seq of

eight tumors complementing their single-cell transcriptomes. De novo annotation of malignant cell clusters proved largely concordant with scRNA-seq-derived cell populations and included an

additional group of AC-like (AC-like-alternative) cells with increased gene activity scores for synaptic marker genes (for example, _GABBR2, GRIA1_ and _CAMK2B_) (Fig. 4a,b, Extended Data

Fig. 4a–f and Supplementary Table 4; Supplementary Note). Cross-modality integration with scRNA-seq data further demonstrated overall congruence between chromatin- and transcriptome-defined

cell states (Extended Data Fig. 4g). Notably, this also revealed distinct clusters of OPC-like-1, OPC-like-2 and OPC-like-3 cells in snATAC-seq space (Extended Data Fig. 4h). Concordant with

our scRNA-seq findings, OPC-like-2 and OPC-like-3 cells also exhibited similarities with pre-OPCs at open chromatin level, whereas OPC-like-1 cells depicted higher chromatin accessibility

for genes also described in healthy committed OPCs45 (Extended Data Fig. 4h–j). Thus, our finding of different OPC-like subpopulations is represented at both transcriptome and accessible

chromatin levels. As snATAC-seq resolves gene-distal and intragenic accessible chromatin regions containing potential _cis_-regulatory DNA elements (CREs) that underlie gene expression, we

next inferred putative CREs integrating snATAC-seq and scRNA-seq modalities. By correlating snATAC-seq-derived accessible chromatin regions/peaks to scRNA-seq measured expression levels of

their nearest associated gene (Fig. 4c)47,48, we identified 13,632 potential peak-gene links of CREs and their target genes (Supplementary Table 4 and Extended Data Fig. 4k). Among these,

287 genes exhibited more than eight (top 5%) linked CREs, denoting high regulatory locus complexity that has been described as ‘predictiveʼ chromatin and thereby a determinant of key lineage

marker genes (Fig. 4d,e)47,48. We identified a higher number of genes linked with predictive chromatin (termed ‘GPCsʼ) specific to OPC/OC-like as compared to AC-like/MES-like cells,

indicating highly cooperative regulation of the oligodendroglial lineage pervasively underlying H3-K27M DMGs (Fig. 4e; Methods). Because large groups of CREs are related to the concept of

‘super-enhancersʼ47,48, we overlayed our candidate GPCs with H3-K27ac ChIP-seq derived super-enhancer profiles of H3-K27M primary tumors19. This demonstrated a significant overlap of GPCs

with H3-K27M DMG super-enhancer regulated genes (Fig. 4f), and further points toward a key role of these multimodally derived marker genes in orchestrating H3-K27M tumor cell identities. We

next sought to reconstruct and refine interdependent circuits of gene regulation by integrating expressions and activities of TFs inferred from scRNA-seq and enrichment of TF binding motifs

in CREs derived from snATAC-seq (Methods). We identified 65 putative cell state-specific TFs that our analysis indicated to be (1) expressed at sufficient levels, (2) binding to

characteristic motifs substantially enriched in CREs and (3) altering expressions of downstream target genes in a cell-type-specific manner (Fig. 4g). Moreover, we examined which TFs

potentially regulate GPCs, focusing on TFs predicted to regulate expressions of GPCs and having binding sites detected within GPC-linked CREs. For example, the OPC-like marker gene _SEZ6L_

is differentially expressed and accessible in OPC-like cells, and is linked to 16 CREs containing TF binding sites of SOX8, which is again predicted to be differentially active in OPC-like

cells (Fig. 4h). We describe the same interdependencies between gene expression, chromatin accessibility and enrichment of cell state-specific TFs in CREs for GPCs of all tumor cell states,

such as for AC-like marker gene _ITPKB_ (Fig. 4i), which is linked to 11 CREs that harbor TF binding sites for SOX9, NFATC4 and RFX3, whose regulons are predicted to govern the expression of

_ITPKB_. Together, our data further corroborate the closely interwoven and cell state-specific loops of chromatin regulation and gene expression identified at multiple levels. In summary,

we show that single-cell chromatin accessibility independently recapitulates the main cellular lineages identified in corresponding single-cell transcriptomes of H3-K27M DMG tumors. Our

multimodal analysis reveals putative cell state-specific CREs as building blocks of larger GPC-associated regulatory complexes. These GPCs are enriched in OPC-like/OC-like cells, reinforcing

the central role of the oligodendroglial lineage in H3-K27M DMGs. These results can be leveraged to more deeply investigate select key intrinsic regulators of H3-K27M DMG cell identities.

THE AGE-SPECIFIC MYELOID CELL LANDSCAPE IN H3-K27M DMGS Various cellular and structural components constitute the glioma microenvironment and extrinsically influence glioma cell

identities49,50. It remains to be elucidated whether these components are characteristic of their respective location or age-related brain environments. Here our age- and location-stratified

H3-K27M glioma cohort uniquely lends itself to dissecting such context-specific differences largely independent of tumor subtype and genetic drivers. Because glioma- or tumor-associated

myeloid cells (GAMs/TAMs) presented the largest proportion of nonmalignant cells within our scRNA-seq dataset (Fig. 1c), we focused on characterizing and comparing this microenvironmental

component across our clinico-anatomical patient groups. We classified TAMs into brain-resident microglia or monocyte-derived macrophages using reported sets of canonical marker genes35 (Fig.

5a–c). Overall TAM proportions were not different between adult and pediatric samples (Extended Data Fig. 5a). However, comparison of microglia versus macrophage proportions across age

groups revealed a higher rate of microglia in pediatric DMGs, while adult DMGs contained higher rates of macrophages (Fig. 5d). Tumor location did not seem to influence these proportions

(Extended Data Fig. 5b). Mounting evidence suggests a causal role of TAMs in establishing a mesenchymal cell state in GBM through TAM-secreted ligands binding to receptors on glioma cells,

such as between ligand-receptor pair OSM-OSMR, or via chemokine signaling27,28,40,51. Given the significant enrichment of MES-like cells in adults compared to pediatric H3-K27M DMGs in our

cohort, we hypothesized that this may be driven by differences in such tumor–immune interactions. We indeed detected higher expression of _OSM_ in adult TAMs, and the corresponding receptor

_OSMR_ in adult tumor cells (Fig. 5e,f), indicating immune-mediated engagement of a previously validated pathway27 in inducing the MES-like phenotype in adult tumors. Moreover, we observed

increased expression of MES-like marker genes in adults compared to pediatric TAMs, which were shown to be increased in mesenchymally enriched gliomas27 (Fig. 5g). To assess whether these

transcriptional differences of MES-like state marker genes and inducing ligands may be inherent to normal brain myeloid cells during temporal development and aging, we analyzed gene

expressions across age in a normal mouse brain myeloid cell atlas. Indeed, we observed an increase of ligands such as OSM and of mesenchymal marker genes with age (Fig. 5h)52, supporting

that the increase of the H3-K27M DMG tumor MES-like state with age is linked to changes of the brain myeloid compartment that also occur during normal development and aging processes. Last,

we interrogated receptor–ligand interactions between TAMs and OPC-like subpopulations, revealing shared OPC-wide (for example, SEMA3E-PLXND1) and subpopulation-specific interactions

(Extended Data Fig. 5c–f). This may point toward a harnessing of microenvironmental factors in reinforcing the OPC-like lineage and further determining their varying maturation, which

provides the basis for follow-up investigations to better understand the contributions of cell-extrinsic regulators to the different OPC-like states. In summary, we reveal that adult H3-K27M

DMGs harbor higher proportions of monocyte-derived macrophages, while pediatric tumors are enriched for brain-resident microglia. We also show that H3-K27M DMG-associated TAMs upregulate

ligands and marker genes that can induce tumor cell MES-like states with increasing age, thereby linking the age-specific tumor immune microenvironment to the observed increase of MES-like

tumor cells in adult H3-K27M DMGs. This illustrates how age-related microenvironmental factors can differentially shape tumor cellular states. CHARTING THE SINGLE-CELL SPATIAL ARCHITECTURE

OF H3-K27M DMG To map our scRNA-seq/snATAC-seq derived cell populations to their spatial positions within intact H3-K27M DMG tissues, we performed hybridization-based in-situ sequencing

(HybISS)53 in 16 patient H3-K27M DMG tissue sections (14 different tumors, 2 tumors with multi-region sampling), using a panel of 116 cell-type-specific combinatorial marker genes curated

from our scRNA-seq dataset (Fig. 1b, 6a–c, Extended Data Fig. 6a–e and Supplementary Table 5). We analyzed spatial cell state compositions by probabilistic cell typing by in situ sequencing

(pciSeq). Here we interestingly observed AC-like cells to constitute the major malignant cell compartment (Fig. 6c), which is in contrast to the predominance of OPC-like cancer cells

observed by scRNA-seq. This held true across tumor sections of different sizes, cell densities and qualities. Our spatial analysis also identified larger numbers and diversity of

nonmalignant cell types, that were either not detected or showed only low representation in scRNA-seq (Extended Data Fig. 6e). Because larger numbers of cells are assessed on average, and

processing-associated biases are reduced in intact tissues, spatial transcriptomics is likely more representative of true cell state compositions than conventional scRNA-seq. Stratification

within our spatially profiled cohort again revealed that adult H3-K27M DMG sections harbor substantially higher proportions of MES-like tumor cells relative to pediatric tumors (Fig. 6d),

orthogonally underscoring the association of age with the MES-like state. We next performed neighborhood enrichment analyses to investigate spatial relationships between individual cell

populations. Here we observed marked variability in neighborhood structures, highlighting overall intertumoral spatial heterogeneity (Supplementary Fig. 3). Global analysis of malignant cell

neighborhoods indicated higher colocalization of OPC-like/cycling and OC-like cells (Fig. 6e). We validated these findings on the protein level by multiplexed immunofluorescence (IF)

imaging (codetection by indexing (CODEX) system) in four H3-K27M gliomas (Fig. 6f and Extended Data Fig. 6f–h). Concordantly, this approach indicated a preferred mitotic niche of

proliferating OPC-like and OC-like cells, encircled by more differentiated, nonproliferating AC-like cells (Fig. 6f). Neighborhood analysis between cancer and noncancer cells revealed closer

proximities between vascular cells and MES-like tumor cells (Extended Data Fig. 6i), pointing toward increased vascularization that has been associated with the mesenchymal state54. Within

a subset of samples (7 of 12 with >1,000 cells profiled), we also observe increased colocalization of microglia/macrophages with MES-like, OC-like and AC-like cancer cells (Supplementary

Fig. 3). Further, we assessed the tendencies of each cell population to either form their own homogeneous cluster, by calculating their clustering coefficient (that is, degree to which

members of a cell population favor clustering together), or to cluster heterogeneously with other populations, as represented by their degree centrality (that is, ratio of nonmembers

connected to members of a cell population). Here we observed that AC-like cells, nonmalignant astrocytes and TAMs depicted the highest tendency to cluster with other cell types/states,

hinting at their more diffuse distribution rather than localization within a restricted spatial compartment. In comparison, vascular cells, neurons, and cycling OPC-like cells exhibited a

higher tendency to cluster with members of the same cell population, which is further indicative of a propensity to form specific structures/niches (Fig. 6g). In summary, we resolved the

spatial architecture of scRNA-seq–defined H3-K27M DMG cell populations directly within the native tumor tissue. Our results shed light on global and heterogeneous cellular relationships and

neighborhoods, notably suggesting the presence of mitotic stem-like niches in which H3-K27M tumor cells of oligodendroglial lineage (OPC-like and OC-like cancer cells) colocalize. These

findings lend themselves to further investigation of potential therapeutic avenues directed at regional and temporal perturbation of H3-K27M DMG tumor cell populations and their associated

niches. DISCUSSION We previously demonstrated the preponderance of OPC-like tumor cells in seven pediatric H3-K27M DMGs through scRNA-seq. However, it remained unknown whether the same

cellular composition—proposed to arise as a function of early pontine development—holds true across multiple spatiotemporal environments in which these tumors occur. To address these

questions, we generated a multi-omic single-cell atlas of H3-K27M DMGs, comprising various midline locations and ages ranging from 2 to 68 years. Our data shed light on understudied thalamic

locations and adolescent/adult age groups and provide a blueprint for the spatiotemporal context-specificity of tumor cell-intrinsic properties, spatial tissue architectures, and

microenvironmental interactions that co-orchestrate cellular identity against the backdrop of the shared K27M driver mutation. Our study reveals a ubiquitous presence of OPC-like and more

differentiated glia-like cells across all clinico-anatomical groups. Concomitantly, neuronal-like tumor cells are absent, which is independent of age and location and stands in contrast to

other glioma types. Thus, this likely presents direct consequences of the K27M mutation universally skewing tumor cells toward an OPC-like and away from a neuronal-like state, decoupled from

spatiotemporal influences. We identify two major variable features as a function of regional or temporal context, respectively (Fig. 7): First, we resolve pontine H3-K27M DMGs to harbor

more immature pre-OPC-like tumor cells than their thalamic counterparts. This raises the question of whether this diversity reflects region-specific cell-intrinsic features or it is driven

by local environmental interactions. While normal murine OPCs have been shown to lack heterogeneity across different brain regions44,55, it is possible that region-specific microenvironments

provide distinct cues to differentially foster OPC differentiation. This has been observed in the gray matter where OPC differentiation takes place more slowly compared to white

matter56,57. In glioblastoma, the white matter has likewise been suggested as a prodifferentiative niche for oligodendroglial lineage stem-like cells58. It will be of interest to explore in

future studies what extrinsic factors in the pons relative to the thalamus may contribute to preserving healthy and aberrant OPC(-like cell)s in a less committed pre-OPC(-like) state and how

these specific microenvironmental contexts could be perturbed by targeting such factors. The finding of a more immature precursor-like cell is accordant with previous modeling studies

postulating embryonic neural stem/progenitor cells instead of OPCs as the H3-K27M DMG cell of origin11,59,60,61,62. While the K27M mutation could occur in such an earlier state, it

subsequently induces a cellular arrest in a self-renewing OPC-like state59, and the hypothesized original cell of mutation may become diluted and eliminated from fully transformed tumors9.

Taken together, the literature supports the idea that the cell state of transformation is an oligodendroglial lineage precursor, whose precise state may vary from pre-OPC to more mature OPC

with different histone variants19, anatomical locations and ages. Second, we observe the mesenchymal signature to increase with higher age, which we link to age-related differences in TAMs

that have been illustrated to induce this myeloid-affiliated tumor signature27,28,40. As the mesenchymal state has been associated with a more aggressive phenotype in a broad range of solid

tumors63,64, and mesenchymal- and myeloid-directed therapies are under active investigation, it will be of interest to investigate such an age and outcome association in H3-K27M gliomas and

other tumors. Lastly, we reconstructed the single-cell spatial architecture of patient H3-K27M tumors, identifying a niche of proliferating OPC-like/OC-like tumor cells, surrounded by

AC-like cells, which constitute the major tumor cell population in situ. This finding contrasts the predominance of OPC-like cells observed by conventional and especially fresh scRNA-seq and

may arise due to technical and biological reasons. As AC-like glioma cells have been shown to be interconnected through tumor microtubes19,25,65, we speculate that they may be less viable

and more sensitive to tumor dissociation, thereby biasing toward capturing more aggressive OPC-like cells in scRNA-seq. By contrast, AC-like cells may be better preserved in frozen snRNA-seq

and spatial approaches. Such a potential predominance of AC-like cells instead of OPC-like cells does not stand in contrast to the proposed role of OPC-like cells as the stem-like drivers

of H3-K27M DMGs and would align with a more traditional model in which cancer stem cells present the minority of tumor cells66. With the emergence of spatial technologies, it will be

relevant to assess whether similar differences are observed throughout other tumor types and biological systems, pinpointing the importance of multimodal profiling to further refine models

derived primarily through the lens of a single modality. Altogether, we provide an extensive resource of H3-K27M DMG cellular heterogeneity across space and time that lends itself to

delineating the multi-faceted interplay between spatiotemporal context-specific cellular properties and microenvironmental niches for the design of rational modeling studies and therapeutic

frameworks tailored to the different clinico-anatomical groups of this lethal glioma. METHODS HUMAN SUBJECTS AND ETHICAL CONSIDERATIONS All samples used in this study were deidentified and

obtained with properly informed consent of patients and/or their legal representatives, who did not receive compensation. The study was approved by the Institutional Review Board at Boston

Children’s Hospital/Dana-Farber Cancer Institute (DFCI 10-417) and at affiliated research hospitals or via waiver of consent as appropriate. Clinical information (age, sex and location) and

mutation status are presented in Fig. 1b and Supplementary Table 1. TUMOR TISSUE COLLECTION AND DISSOCIATION Fresh tumor tissue acquired at the time of surgery was immediately mechanically

and enzymatically dissociated for 30 min at 37 °C using the Brain Tumor Dissociation Kit (Miltenyi Biotec). Single-cell suspensions were filtered through a 70 µm strainer, centrifuged at

500_g_ for 5 min, and resuspended in PBS/1% BSA for fluorescence-activated cell sorting (FACS). To extract single nuclei from frozen tissues for snRNA-seq, snap-frozen or OCT-embedded tumor

tissue was disaggregated on ice in 1 ml 0.49% CHAPS detergent-based nuclear extraction buffer67, aided by mild chopping. Single-nuclei suspensions were filtered using a 40 µm strainer and

centrifuged at 500_g_ for 5 min. All steps were performed at 4 °C. To prepare single-nuclei suspensions for snATAC-seq, snap-frozen DMG tissue was lysed on ice in lysis buffer (10 mM

Tris–HCl, 10 mM NaCl, 3 mM MgCl2, 1% BSA, 0.01% Tween-20, 0.01% NP-40, 0.001% digitonin) under mild chopping for 5 min, followed by ten times mixing using a wide-bore pipette tip and 10 min

incubation on ice. Wash buffer (10 mM Tris–HCl, 10 mM NaCl, 3 mM MgCl2, 1% BSA, 0.1% Tween-20) was added and mixed five times before filtering through 70 and 40 µm Flowmi cell strainers.

Single-nuclei suspensions were then centrifuged at 500_g_ for 5 min at 4 °C, resuspended in 1× diluted nuclei buffer and counted. SCRNA-SEQ/SNRNA-SEQ DATA GENERATION Whole transcriptome

amplification, library preparation and sequencing of single cells/nuclei were performed using the Smart-seq2 modified protocol21,33,35,68,69. RNA was purified with Agencourt RNAClean XP

beads (Beckman Coulter). Oligo-dT primed reverse transcription (RT) was performed using Maxima H Minus reverse transcriptase (Life Technologies) and a template-switching oligonucleotide

(TSO; Qiagen). PCR amplification (20 cycles for scRNA-seq and 22 cycles for snRNA-seq) was performed using KAPA HiFi HotStart ReadyMix (KAPA Biosystems), followed by Agencourt AMPure XP bead

(Beckman Coulter) purification. Libraries were generated using the Nextera XT Library Prep kit (Illumina). Libraries from 768 cells with unique barcodes were combined and sequenced using a

NextSeq 500/550 High Output Kit v2.5 (Illumina). SCATAC-SEQ DATA GENERATION scATAC-seq libraries were generated using the 10X Chromium Controller and Chromium Next GEM Single Cell ATAC &

Library Gel Bead Kit v1.1 kit according to the manufacturer’s instructions (Document CG000209). Briefly, 7,000–10,000 nuclei were tagmented at 37 °C for 60 min and loaded on a Chromium Next

GEM Chip H and Chromium Controller for generation of single-cell Gel Bead-In-Emulsions, followed by linear amplification of barcoded tagmented DNA. GEMs were then broken up, DNA fragments

were purified using Dynabeads MyOne SILANE (10X 2000048) and SPRIselect Reagent (Beckman Coulter, B23318), and further PCR-amplified for 10–11 cycles undergoing sample indexing. Libraries

were sequenced using a NextSeq 500/550 High Output Kit v2.5 (Illumina) at targeted 25,000 reads per cell. GENE SELECTION FOR TARGETED HYBISS Gene panel selection was based on the scRNA-seq

data from ten H3-K27M DMG patient tumors spanning multiple clinico-anatomical groups and on published datasets of normal brain-resident cell types70 (Supplementary Table 5). Genes were

prioritized based on differential expression between cell types, followed by manual filtering of genes with likely high background expression levels being strongly expressed in all cell

types. A total of 618 probes were designed for 116 genes encompassing malignant (OPC-like, AC-like, OC-like and MES-like) and nonmalignant cells (oligodendrocytes, astrocytes, neurons,

macrophages, microglia, T cells, endothelia, pericytes and ependymal cells) (Supplementary Note 1). HYBISS After fixation with 3% PFA for 30 min, sections were permeabilized with 0.1 M HCl

and washed with PBS. After rehydration for 1 min in 100% ethanol, 1 min in 75% ethanol and 1 min in PBS, cDNA was synthesized overnight with reverse transcriptase (BLIRT), RNase inhibitor,

and primed with random decamers. Sections were postfixed before padlock probe (PLP) hybridization and ligation at a final concentration of 10 nM/PLP, with Tth Ligase and RNaseH (BLIRT). This

was performed at 37 °C for 30 min and then 45 °C for 1 h. Sections were washed with PBS, followed by rolling circle amplification (RCA) with phi29 polymerase (Monserate) and Exonuclease I

(Thermo Fisher Scientific) overnight at 30 °C. Bridge probes (10 nM) (Supplementary Table 6) were hybridized at RT for 1 h in hybridization buffer (2× saline- sodiumcitrate buffer (SSC), 20%

formamide), followed by hybridization of readout detection probes (100 nM) and DAPI (Biotium) in hybridization buffer for 1 h at RT. The sections were washed with PBS and mounted with

SlowFade Gold Antifade Mountant (Thermo Fisher Scientific). After each imaging round, coverslips were removed and sections were washed five times with 2× SSC. Bridge probe/detection

oligonucleotides were then stripped with 65% formamide and 2× SSC for 30 min at 30 °C, followed by five washes with 2× SSC. The above procedure was repeated for cycles 1 through 5, leading

to hybridization of cycle-specific individual bridge probes (for imaging, see Supplementary Note 1). CODEX FFPE tissue sections were collected onto poly(l-lysine)-coated coverslips and

prepared according to the Akoya Biosciences CODEX protocol71. Sections were then deparaffinized and rehydrated. Antigen retrieval was performed using a pressure cooker and 1× citrate buffer,

pH 6.0. Sections were then quenched for autofluorescence72, and subsequently stained and imaged using the Akoya Biosciences CODEX staining kit (7000008). Tissue was stained using the

following preconjugated antibodies purchased from Akoya: DAPI (7000003), Ki67-BX047 (B56)—Atto 550-RX047 (4250019) 1:200, CD44-BX005 (IM7)—Atto 550-RX005 (4250002) 1:50. The following

antibodies were custom conjugated using the Akoya Biosciences conjugation kit (7000009) and indicated barcodes: anti-PDGFRα antibody (Abcam, ab234965) Barcode BX002—Atto 550-RX002 (5450023)

1:50, anti-BCAS1 antibody (Santa Cruz Biotechnology, sc-136342) Barcode BX027—Cy5-RX027 (5350004) 1:50, anti-GFAP antibody (Invitrogen, 13-0300) Barcode BX030—Cy5-RX030 (5350005) 1:50,

antihistone H3 (mutated K27M) antibody (Abcam, ab240310) Barcode BX004—Alexa FluorTM 488-RX004 (5450014) 1:100, anti-IBA1 antibody (Thermo Fisher Scientific, GT10312) Barcode BX020—Atto

550-RX020 (5250002) 1:50, anti-CD63 antibody (353039, Biolegend) Barcode BX029—Atto 550-RX029 (5250005) 1:50. Imaging was performed using a Keyence BZ-X800E fluorescent microscope equipped

with a BZ Nikon Objective Lens (×20). Images were processed using the CODEX processor software (Akoya) and visualized using the ImageJ plugin CODEX Multiplex Analysis Viewer. STATISTICS AND

REPRODUCIBILITY No statistical method was used to predetermine the sample size. No data were excluded from the analyses. The experiments were not randomized. Data collection and analysis

were not performed blind to the conditions of the experiments. Statistical analysis was performed in R v.4.0.3. A Bayesian statistical framework scCODA (v0.1.4) was used to identify changes

in the proportion of different cell populations between age groups and anatomical departments. Comparisons of numerical variables between different conditions were carried out using Wilcoxon

rank-sum test and Kolmogorov–Smirnov test, as appropriate. Overlap between groups of genes was assessed using a Hypergeometric test. Single-cell sequencing for each tumor was performed in

one experimental replicate. This is typical for human studies because tissues are usually limited and cannot be analyzed more than once. At least three samples per age and anatomical group

were collected to verify reproducibility. The ISS and IF experiments for each tumor sample were performed in one experimental replicate, where the entire section was imaged. For RNAish

experiments, two to three slides were stained per sample and approximately 10–15 fields of view were captured per slide. Further information on research design is available in the Nature

Research Reporting Summary. SCRNA-SEQ DATA PROCESSING We aligned raw sequencing reads to hg19 genome by hisat2 (v2.1.0) and quantified and normalized gene counts using RSEM (v1.3.0) as

transcript-per-million/TPM73. For snRNA-seq data, we modified the gene annotation files to count introns74. We calculated expression levels as _Ei_,_j_ = log2(TPM_i_,_j_/10 + 1) for gene _i_

in sample _j_. To filter out low-quality cells in fresh samples, we removed cells with <2,000 genes or an average housekeeping gene expression of <2.5. For frozen tumors, a filtering

threshold of <1,000 genes and an alignment rate of <0.4 were employed. In sum, 9,911 high-quality cells were retained. We also removed genes with TPM > 16 in <10 cells. For the

remaining cells and genes, we computed the aggregate expression of each gene as Ea(_i_) = log2(average(TPM_i_,1…_n_) + 1) and defined relative expression as centered expression levels,

Er_i_,_j_ = E_i_,_j_ − average(E_i_,1…_n_). On average, we detected 6,866 uniquely expressed genes per cell in fresh, and 4,432 uniquely expressed genes in frozen tumors. DATA HARMONIZATION,

LOUVAIN CLUSTERING AND IDENTIFICATION OF DIFFERENTIALLY EXPRESSED GENES Graph-based clustering with data integration was adapted for independent identification of cellular clusters and gene

signatures. We selected highly variable genes (HVGs) using Seurat (v3.2.2)75 and used the relative expression values of these HVGs for PCA. To disentangle sample-specific biological

variations (that is, tumor-specific genetic and epigenetic alterations) from cell subpopulation-specific variations and to integrate multiple samples, we applied a linear adjustment method

(Harmony v1.0) to the first 100 PCs with default parameters to generate a corrected embedding76. We chose the first 20 Harmony-corrected dimensions for uniform manifold approximation and

projection embedding (UMAP) embeddings, and clustered cells by Seurat’s Louvain algorithm-based FindClusters function. Cells from different samples expressing similar gene programs were well

mixed (Extended Data Fig. 2a). We next identified differentially expressed genes by Seurat’s FindAllMarkers function. We tested genes that were detected in a minimum of 30% of the cells

within each cluster and that showed at least a 0.5-fold mean log difference. We utilized Wilcoxon rank-sum test with Bonferroni correction for multiple testing and only kept genes with

adjusted _P_ value < 0.05. NONNEGATIVE MATRIX FACTORIZATION (NMF) METAPROGRAM ANALYSIS NMF was used to assemble transcriptional programs from relative expressions (with negative values

converted to zero)21,68,69. We derived NMF programs for malignant cells from each sample using the top 10,000 over-dispersed genes, as determined by PAGODA2 (v0.1.4)77. The number of factors

was set to six for each sample. Because redundant NMF programs were merged into a single metaprogram, the final metaprogram was not sensitive to the initially chosen number of factors. We

selected the top 30 genes with the highest NMF weights from each NMF factor and scored all malignant cells with these NMF programs. We then clustered NMF programs by hierarchical clustering

(distance metric: 1 − Pearson correlation; linkage: Wardʼs linkage) on the scores for each NMF program (Extended Data Fig. 2b). This revealed eight highly correlated sets of programs in

fresh tumors and nine in frozen tumors. We merged these correlated programs into metaprograms by selecting the top 30 genes with the highest average NMF weight within each correlated program

set (Supplementary Table 2 and Supplementary Note). COMPARISON AND INTEGRATION OF FRESH AND FROZEN TUMOR METAPROGRAMS We compared transcriptional metaprograms independently derived from

fresh and frozen tumors by pairwise correlation analysis, showing high correlations between the cycling, fresh OPC-like-1/frozen OPC-like-a, OC-like, AC-like and MES-like signatures

(Extended Data Fig. 2d). Even though ribosomal protein-encoding genes marking the fresh OPC-like-2 metaprogram were filtered out in the frozen dataset to exclude potential technical

artifacts from random capture of nuclei-associated ribosomes67, the frozen OPC-like-b program showed high correlation with the fresh OPC-like-2 signature (Extended Data Fig. 2d) and showed

higher expression of pre-OPC markers, such as _DLL1, HES6_ and _EGFR_ (Supplementary Fig. 2). Therefore, we independently identified pre-OPC-like cells in our fresh and frozen

scRNA-seq/snRNA-seq data. We consequently scored frozen nuclei for all fresh metaprograms, only exchanging fresh OPC-like-2 with frozen OPC-like-b to avoid artifacts due to the filtering of

ribosomal protein genes. If the resulting maximum expression score was <0.2, single nuclei were classified as ‘score_too_lowʼ; if ≥0.2, nuclei were assigned according to the

highest-scored metaprogram (Extended Data Fig. 2e,f). ANALYSIS OF CELL TYPE COMPOSITIONS We applied the Bayesian model-based single-cell compositional data analysis (scCODA v0.1.4) framework

to identify associations of cell compositions with different clinical covariates78. scCODA employs hierarchical Dirichlet-multinomial distribution that accounts for the uncertainty and

negative correlative bias in compositional analysis of cell type proportions. The model uses a logit-normal spike-and-slab prior with a log-link function and Hamiltonian Monte Carlo sampling

to estimate the effects of covariates on cell type proportions. The sample level counts of cell annotations and clinical covariates were used as inputs for scCODA. The default parameter was

used with AC-like cells selected as the reference cell type. Locations and ages were included as covariates in the model. The statistical significance of changes in cell compositions was

assessed using credible effects with a 5% false discovery rate. SNATAC-SEQ DATA PROCESSING Cell Ranger ATAC (v1.0.1) was used to process 10X Chromium snATAC-seq data. We used cellranger-atac

counts to generate single-cell accessibility counts and cellranger-aggr to aggregate multiple samples without setting any normalization. The resulting peak-cell matrix and metadata were

then analyzed in Signac (v1.1.0)79. We removed nuclei with <200 detected peaks and peaks detected in <10 nuclei. We further kept nuclei with the following: (1) total number of

fragments in peaks (peak_region_fragments) between 1,500 and 15,000, (2) percent of reads in peaks (pct_reads_in_peaks) >15, (3) ratio of reads in genomic blacklist regions

(blacklist_ratio) <0.02, (4) approximate ratio of mononucleosomal to nucleosome-free fragments (nucleosome_signal) <2 and (5) ratio of fragments centered at the transcription start

site (TSS) to fragments in TSS-flanking regions (TSS_enrichment) >4. After quality control and filtering, a dataset comprising 211,096 peaks and 9,797 nuclei was used for downstream

analysis. We normalized data using term frequency-inverse document frequency (RunTFIDF) and conducted dimensionality reduction using singular value decomposition and top 25% of features. We

calculated k-nearest neighbors using FindNeighbours (reduction = ‘lsiʼ, dims = 2:30) and omitted the first latent semantic indexing (LSI) component as it exhibited a strong correlation with

sequencing depth. We then identified cell clusters by shared nearest neighbor modularity optimization-based clustering algorithm and ran the FindClusters function (algorithm = 3/SLM and

resolution = 0.8), and generated a UMAP embedding using the RunUMAP function with 2–30 LSI components. We calculated gene activities for each gene in each nucleus by summing the peak counts

in the gene body + promoter region (2 kb upstream of TSS). We then normalized gene activities to the median of total gene activities and performed log transformation. Genes with differential

activities (DAGs) were identified by running FindAllMarkers on normalized gene activities. We tested genes that were detected in a minimum of 20% of the cells within each cluster by

Wilcoxon rank-sum test with Bonferroni multiple test correction and only kept genes with log fold change >0.1 and adjusted _P_ < 0.05. Top DAGs were used for initial annotation of each

cell cluster. Putative nonmalignant clusters with highly accessible canonical marker genes were identified, including microglia (for example, _CD14_, _CSF1R_ and _SPP1_), T cells (for

example, _CD2_, _CD3D_ and _RHOH_) and tumor-associated oligodendrocytes (for example, _BCAS1_, _SOX10_ and _SIRT2_). SCRNA-SEQ/SNATAC-SEQ DATA INTEGRATION We applied canonical correlation

analysis as implemented in Seurat to integrate log normalized gene activity scores of ATAC-seq data and gene expression scores of RNA-seq data. We used Seurat’s ‘FindTransferAnchorsʼ

function for integration, specified the union of the 2,764 and 2,000 most variable genes in scRNA-seq and snATAC-seq respectively as input features, ‘ccaʼ as the reduction method, and

default values for the rest of the parameters. For each cell profiled by snATAC-seq, we identified the nearest neighbor cell in those profiled by scRNA-seq with a nearest-neighbor search in

the joint canonical correlation (CCA) L2 space. Nearest neighbors were determined by the ‘FNNʼ R package with the ‘kd_treeʼ algorithm. LINKING GENE REGULATORY ELEMENTS AND GENE EXPRESSION

ACROSS ALL CELL TYPES Because RNA expressions and chromatin accessibilities were measured in different cells, we applied a correlation-based approach to pseudobulk samples aggregating

snATAC-seq and scRNA-seq counts from computationally matched cells to identify peak-to-gene links as putative CREs. We defined pseudobulk samples by randomly sampling 200 cells from the

snATAC-seq dataset and combined each of these 200 seed cells with their respective 99 nearest neighbor cells in the Harmony-corrected ATAC-LSI space. Hence, each of the resulting pseudobulk

sample comprised 100 cells. We computed pseudobulk peak counts by summing peak counts across respective counts of all 100 cells within each pseudobulk sample. Within each pseudobulk, we

matched 100 ATAC cells with 100 RNA cells as their nearest neighbors in CCA L2 space and obtained pseudobulk RNA gene counts by summing gene counts across the respective counts of all 100

cells within each pseudobulk sample. Pseudobulk gene counts were normalized as TPM. We then defined putative peak-gene pairs by associating peaks with a genomic distance within 250 kb of the

TSS of genes profiled by scRNA-seq. Each peak is only linked to its nearest gene. For each candidate peak-gene pair, we determined the Pearson correlation coefficient of peak counts

(normalized as CPM) and gene expression (TPM), and adjusted _P_ values for these coefficients from a _t_-statistic using Benjamini-Hochberg (BH) procedure. We identified a set of 13,632

high-confidence peak-to-gene links by only retaining pairs with |PCC|>0.2 and BH-adjusted _P_ < 0.05. INTEGRATIVE TF ANALYSIS We integrated scRNA-seq and scATAC-seq data to identify

putative regulatory networks of TF-target pairs. For each TF documented in the JASPAR (2020) TF motif database, we computed its mean expression (TPM) and examined the frequency of its

motif(s) within the CREs located in the TSS ± 10 kb region of its predicted target genes by SCENIC80. We then kept TFs with mean TPM > 4 and over-represented binding motifs in CREs. Next,

we kept TFs that were among the top 30 TF regulons with the highest specificity score of any cell type. This resulted in a total of 65 TFs (Supplementary Table 4). Of these TFs, 19 were

specific to OPC-like cells (for example, EGR1, JUN, HES6), 10 were specific to OC-like cells (for example, SOX4, SOX10), 21 were specific to AC-like cells (for example, GLI2, STAT3 and SOX9)

and 15 were specific to MES-like cells (for example, FOSL2, CEBPD and ELK3). For each GPC, we leveraged two complementary approaches to identify core TFs that may regulate expressions of

this gene. First, we selected TFs that were predicted to regulate expressions of the target GPC by SCENIC analysis. Second, we examined if TFs identified above possess binding motifs that

are over-represented in the CREs linked to the target GPC using a hypergeometric test. We kept TFs that are predicted to govern the expression of a target GPC and harbor binding motifs

substantially enriched in CREs linked to the target GPC (Supplementary Table 4). ANALYSIS OF HYBISS DATA Image processing and decoding. Each field of view (FOV) image was maximum intensity

projected to obtain a flattened two-dimensional image. These images were then analyzed using in-house custom software that handles image processing and gene calling based on the python

package Starfish v0.2.1 (ref. 81). Each two-dimensional FOV was exported, and preprocessed including alignment between cycles, and stitched together using the MIST algorithm. Stitching was

followed by retiling to create smaller nonoverlapping 6,000 × 6,000 pixel images that were then used for decoding. The decoding pipeline can be found at

https://github.com/Moldia/iss_starfish/. Using Starfish, images were initially filtered by applying a white top hat filter. The filtered images were subsequently normalized, and spots were

then detected using the FindSpots module from Starfish and decoded using MetricDistance decoding. MALIGNANT VERSUS NONMALIGNANT CELL TYPING To distinguish between malignant (H3-K27M

positive) and nonmalignant (H3-K27M negative) cells, ISS expression maps were aligned to IF images, both taken from the same tissue section, and the mean IF intensity of each cell was

calculated. All IF H3-K27M positive cells were categorized as malignant based on a minimum IF threshold in each sample, while DAPI positive and H3-K27M negative cells were categorized as

nonmalignant based on a maximum IF threshold. Cells with IF intensities between the two thresholds were considered ambiguous and excluded from the analysis. We obtained spatial

transcriptomic profiles of a total of 125,801 high-quality cells (56,664 malignant cells, 69,137 nonmalignant cells). PCISEQ To identify the cellular identity of nonmalignant and cancer

cells, two different methods were applied. Probabilistic cell maps of malignant cells were created using pciSeq v0.0.45. The pciSeq pipeline assigns the spatial coordinates of genes from the

ISS maps to DAPI-stained nuclei based on the proximity and assigns individual cells to cell type definitions defined by our H3-K27M DMG scRNA-seq dataset. The pciSeq pipeline is publicly

available (https://github.com/acycliq/pciSeq)82. In contrast, due to the presence of uniquely expressed markers in the panel, nonmalignant cell types were identified by the expression of key

marker genes in each sample. Here we assigned nonmalignant cell types by lack of H3-K27M signal in IF staining and concomitant expression of key markers, such as MBP for oligodendrocytes,

ESAM for endothelial, MYL9 for pericytes, GFAP for astrocytes, CD74 for TAMs, DLG4 for neurons. T cells were excluded from downstream analyses due to very low numbers identified. SPATIAL

ENRICHMENT AND NEIGHBORS ANALYSIS To explore proximities between the different cell types, neighborhood enrichment analysis was performed using Squidpy v1.1.2 (ref. 83). Briefly, the spatial

coordinates of the mapped cells were used to identify spatial enrichment of cell types at a specific radius, and an enrichment score for each defined cell type was calculated based on the

number of connections for each cell cluster. The number of observed connection events was compared against 100 permutations, and a _Z_ score was computed for each cell type that can be

positive (indicating positive colocalization) or negative (indicating negative colocalization). Centrality scores and clustering coefficients were calculated for all samples and each

individual sample as previously indicated83. Degree centrality represents the fraction of nongroup members, establishing each cell type as a group, connected to the cells assigned to the

cell type analyzed. The clustering coefficient represents the degree to which nodes in the graph tend to cluster together. It is formulated as the number of closed triplets, defining a

triplet as three connected nodes, over the total number of triplets. Calculation of scores was implemented in SquidPy v1.1.2 (ref. 83). REPORTING SUMMARY Further information on research

design is available in the Nature Portfolio Reporting Summary linked to this article. DATA AVAILABILITY ScRNA-seq and scATAC-seq data of primary patient DMGs have been submitted to GEO

(GSE184357). ISS data are available at Zenodo under ID 6805729. Previously published scRNA-seq data reanalyzed in this study are available under accession codes GSE102130 (ref. 21),

GSE122871 (ref. 43), GSE144462 (ref. 41), GSE131258 (ref. 46) and GSE123030 (ref. 52). WES data generated in this study are deposited in EGA (EGAS00001006431). For targeted exome-sequencing

data, the majority of which was generated as part of routine clinical care, variant data have been included as Supplementary Table 7 for all samples except for A21–238 and AAA010043 as these

were generated by external care providers with restricted data access. Previously published WGS data of tumors ICGC-GBM27, ICGC-GBM96 and ICGC-GBM60 are deposited at EGA00001001139, and WGS

data for BT836 and BT869 have been published under dbGaP accession number phs002380.v1.p1 (ref. 84). H3-K27M DMG ChIP-seq data were utilized from GSE126319 (ref. 19). CODE AVAILABILITY

Custom scripts and code used in data processing and figure creation are available at https://doi.org/10.5281/zenodo.7073167 (ref. 85). REFERENCES * Schwartzentruber, J. et al. Driver

mutations in histone H3.3 and chromatin remodelling genes in paediatric glioblastoma. _Nature_ 482, 226–231 (2012). Article  CAS  Google Scholar  * Wu, G. et al. Somatic histone H3

alterations in pediatric diffuse intrinsic pontine gliomas and non-brainstem glioblastomas. _Nat. Genet._ 44, 251–253 (2012). Article  CAS  Google Scholar  * Khuong-Quang, D. A. et al. K27M

mutation in histone H3.3 defines clinically and biologically distinct subgroups of pediatric diffuse intrinsic pontine gliomas. _Acta Neuropathol._ 124, 439–447 (2012). Article  CAS  Google

Scholar  * Sturm, D. et al. Hotspot mutations in H3F3A and IDH1 define distinct epigenetic and biological subgroups of glioblastoma. _Cancer Cell_ 22, 425–437 (2012). Article  CAS  Google

Scholar  * Mackay, A. et al. Integrated molecular meta-analysis of 1,000 pediatric high-grade and diffuse intrinsic pontine glioma. _Cancer Cell_ 32, 520–537 (2017). Article  CAS  Google

Scholar  * Louis, D. N. et al. The 2016 World Health Organization classification of tumors of the central nervous system: a summary. _Acta Neuropathol._ 131, 803–820 (2016). Article  Google

Scholar  * Meyronet, D. et al. Characteristics of H3 K27M-mutant gliomas in adults. _Neuro Oncol._ 19, 1127–1134 (2017). Article  CAS  Google Scholar  * Schulte, J. D. et al. Clinical,

radiologic, and genetic characteristics of histone H3 K27M-mutant diffuse midline gliomas in adults. _Neurooncol. Adv._ 2, vdaa142 (2020). Google Scholar  * Filbin, M. & Monje, M.

Developmental origins and emerging therapeutic opportunities for childhood cancer. _Nat. Med._ 25, 367–376 (2019). Article  CAS  Google Scholar  * Monje, M. et al. Hedgehog-responsive

candidate cell of origin for diffuse intrinsic pontine glioma. _Proc. Natl Acad. Sci. USA_ 108, 4453–4458 (2011). Article  CAS  Google Scholar  * Pathania, M. et al. H3.3(K27M) cooperates

with Trp53 loss and PDGFRA gain in mouse embryonic neural progenitor cells to induce invasive high-grade gliomas. _Cancer Cell_ 32, 684–700 (2017). Article  CAS  Google Scholar  * Lewis, P.

W. et al. Inhibition of PRC2 activity by a gain-of-function H3 mutation found in pediatric glioblastoma. _Science_ 340, 857–861 (2013). Article  CAS  Google Scholar  * Bender, S. et al.

Reduced H3K27me3 and DNA hypomethylation are major drivers of gene expression in K27M mutant pediatric high-grade gliomas. _Cancer Cell_ 24, 660–672 (2013). Article  CAS  Google Scholar  *

Venneti, S. et al. Evaluation of histone 3 lysine 27 trimethylation (H3K27me3) and enhancer of Zest 2 (EZH2) in pediatric glial and glioneuronal tumors shows decreased H3K27me3 in H3F3A K27M

mutant glioblastomas. _Brain Pathol._ 23, 558–564 (2013). Article  CAS  Google Scholar  * Piunti, A. et al. Therapeutic targeting of polycomb and BET bromodomain proteins in diffuse

intrinsic pontine gliomas. _Nat. Med._ 23, 493–500 (2017). Article  CAS  Google Scholar  * Stafford, J. M. et al. Multiple modes of PRC2 inhibition elicit global chromatin alterations in

H3K27M pediatric glioma. _Sci. Adv._ 4, eaau5935 (2018). Article  CAS  Google Scholar  * Harutyunyan, A. S. et al. H3K27M induces defective chromatin spread of PRC2-mediated repressive

H3K27me2/me3 and is essential for glioma tumorigenesis. _Nat. Commun._ 10, 1262 (2019). Article  Google Scholar  * Brien, G. L. et al. Simultaneous disruption of PRC2 and enhancer function

underlies histone H3.3-K27M oncogenic activity in human hindbrain neural stem cells. _Nat. Genet._ 53, 1221–1232 (2021). Article  CAS  Google Scholar  * Nagaraja, S. et al. Histone variant

and cell context determine H3K27M reprogramming of the enhancer landscape and oncogenic state. _Mol. Cell_ 76, 965–980 (2019). Article  CAS  Google Scholar  * Krug, B. et al. Pervasive H3K27

acetylation leads to ERV expression and a therapeutic vulnerability in H3K27M gliomas. _Cancer Cell_ 36, 338–339 (2019). Article  CAS  Google Scholar  * Filbin, M. G. et al. Developmental

and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq. _Science_ 360, 331–335 (2018). Article  CAS  Google Scholar  * Jessa, S. et al. Stalled developmental programs at

the root of pediatric brain tumors. _Nat. Genet._ 51, 1702–1713 (2019). Article  CAS  Google Scholar  * Venkatesh, H. S. et al. Neuronal activity promotes glioma growth through neuroligin-3

secretion. _Cell_ 161, 803–816 (2015). Article  CAS  Google Scholar  * Venkatesh, H. S. et al. Targeting neuronal activity-regulated neuroligin-3 dependency in high-grade glioma. _Nature_

549, 533–537 (2017). Article  Google Scholar  * Venkatesh, H. S. et al. Electrical and synaptic integration of glioma into neural circuits. _Nature_ 573, 539–545 (2019). Article  CAS  Google

Scholar  * Venkataramani, V. et al. Glutamatergic synaptic input to glioma cells drives brain tumour progression. _Nature_ 573, 532–538 (2019). Article  CAS  Google Scholar  * Hara, T. et

al. Interactions between cancer cells and immune cells drive transitions to mesenchymal-like states in glioblastoma. _Cancer Cell_ 39, 779–792 (2021). Article  CAS  Google Scholar  *

Gangoso, E. et al. Glioblastomas acquire myeloid-affiliated transcriptional programs via epigenetic immunoediting to elicit immune evasion. _Cell_ 184, 2454–2470 (2021). Article  CAS  Google

Scholar  * Gibson, E. M. et al. Neuronal activity promotes oligodendrogenesis and adaptive myelination in the mammalian brain. _Science_ 344, 1252304 (2014). Article  Google Scholar  *

Michealraj, K. A. et al. Metabolic regulation of the epigenome drives lethal infantile ependymoma. _Cell_ 181, 1329–1345 (2020). Article  CAS  Google Scholar  * Fontebasso, A. M. et al.

Recurrent somatic mutations in ACVR1 in pediatric midline high-grade astrocytoma. _Nat. Genet._ 46, 462–466 (2014). Article  CAS  Google Scholar  * Ebrahimi, A. et al. High frequency of H3

K27M mutations in adult midline gliomas. _J. Cancer Res. Clin. Oncol._ 145, 839–850 (2019). Article  CAS  Google Scholar  * Neftel, C. et al. An integrative model of cellular states,

plasticity, and genetics for glioblastoma. _Cell_ 178, 835–849 (2019). Article  CAS  Google Scholar  * Wang, L. et al. The phenotypes of proliferating glioblastoma cells reside on a single

axis of variation. _Cancer Discov._ 9, 1708–1719 (2019). Article  CAS  Google Scholar  * Venteicher, A. S. et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by

single-cell RNA-seq. _Science_ 355, eaai8478 (2017). Article  Google Scholar  * Tirosh, I. et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. _Nature_

539, 309–313 (2016). Article  Google Scholar  * Chen, C. C. L. et al. Histone H3.3G34-mutant interneuron progenitors Co-opt PDGFRA for gliomagenesis. _Cell_ 183, 1617–1633 (2020). Article 

CAS  Google Scholar  * Van de Sande, B. et al. A scalable SCENIC workflow for single-cell gene regulatory network analysis. _Nat. Protoc._ 15, 2247–2276 (2020). Article  Google Scholar  *

Verhaak, R. G. et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. _Cancer Cell_ 17,

98–110 (2010). Article  CAS  Google Scholar  * Schmitt, M. J. et al. Phenotypic mapping of pathologic cross-talk between glioblastoma and innate immune cells by synthetic genetic tracing.

_Cancer Discov._ 11, 754–777 (2021). Article  CAS  Google Scholar  * Fu, Y. et al. Heterogeneity of glial progenitor cells during the neurogenesis-to-gliogenesis switch in the developing

human cerebral cortex. _Cell Rep._ 34, 108788 (2021). Article  CAS  Google Scholar  * Huang, W. et al. Origins and proliferative states of human oligodendrocyte precursor cells. _Cell_ 182,

594–608 (2020). Article  CAS  Google Scholar  * Weng, Q. et al. Single-cell transcriptomics uncovers glial progenitor diversity and cell fate determinants during development and

gliomagenesis. _Cell Stem Cell_ 24, 707–723 (2019). Article  CAS  Google Scholar  * Marques, S. et al. Transcriptional convergence of oligodendrocyte lineage progenitors during development.

_Dev. Cell_ 46, 504–517 (2018). Article  CAS  Google Scholar  * van Bruggen, D. et al. Developmental landscape of human forebrain at a single-cell level identifies early waves of

oligodendrogenesis. _Dev. Cell_ 57, 1421–1436 (2022). Article  Google Scholar  * Zhong, S. et al. Decoding the development of the human hippocampus. _Nature_ 577, 531–536 (2020). Article 

CAS  Google Scholar  * Ma, S. et al. Chromatin potential identified by shared single-cell profiling of RNA and chromatin. _Cell_ 183, 1103–1116 (2020). Article  CAS  Google Scholar  *

Trevino, A. E. et al. Chromatin and gene-regulatory dynamics of the developing human cerebral cortex at single-cell resolution. _Cell_ 184, 5053–5069 (2021). Article  CAS  Google Scholar  *

Gilbertson, R. J. & Rich, J. N. Making a tumour’s bed: glioblastoma stem cells and the vascular niche. _Nat. Rev. Cancer_ 7, 733–736 (2007). Article  CAS  Google Scholar  * Gutmann, D.

H. & Kettenmann, H. Microglia/brain macrophages as central drivers of brain tumor pathobiology. _Neuron_ 104, 442–449 (2019). Article  CAS  Google Scholar  * Sa, J. K. et al.

Transcriptional regulatory networks of tumor-associated macrophages that drive malignancy in mesenchymal glioblastoma. _Genome Biol._ 21, 216 (2020). Article  CAS  Google Scholar  * Li, Q.

et al. Developmental heterogeneity of microglia and brain myeloid cells revealed by deep single-cell RNA sequencing. _Neuron_ 101, 207–223 (2019). Article  CAS  Google Scholar  * Gyllborg,

D. et al. Hybridization-based in situ sequencing (HybISS) for spatially resolved transcriptomics in human and mouse brain tissue. _Nucleic Acids Res._ 48, e112 (2020). Article  CAS  Google

Scholar  * Kim, Y. et al. Perspective of mesenchymal transformation in glioblastoma. _Acta Neuropathol. Commun._ 9, 50 (2021). Article  CAS  Google Scholar  * Marques, S. et al.

Oligodendrocyte heterogeneity in the mouse juvenile and adult central nervous system. _Science_ 352, 1326–1329 (2016). Article  CAS  Google Scholar  * Vigano, F., Mobius, W., Gotz, M. &

Dimou, L. Transplantation reveals regional differences in oligodendrocyte differentiation in the adult brain. _Nat. Neurosci._ 16, 1370–1372 (2013). Article  CAS  Google Scholar  * Boshans,

L. L., Sherafat, A. & Nishiyama, A. The effects of developmental and current niches on oligodendrocyte precursor dynamics and fate. _Neurosci. Lett._ 715, 134593 (2020). Article  CAS 

Google Scholar  * Brooks, L. J. et al. The white matter is a pro-differentiative niche for glioblastoma. _Nat. Commun._ 12, 2184 (2021). Article  CAS  Google Scholar  * Haag, D. et al.

H3.3-K27M drives neural stem cell-specific gliomagenesis in a human iPSC-derived model. _Cancer Cell_ 39, 407–422 (2021). Article  CAS  Google Scholar  * Larson, J. D. et al. Histone H3.3

K27M accelerates spontaneous brainstem glioma and drives restricted changes in bivalent gene expression. _Cancer Cell_ 35, 140–155 (2019). Article  CAS  Google Scholar  * Cordero, F. J. et

al. Histone H3.3K27M represses p16 to accelerate gliomagenesis in a murine model of DIPG. _Mol. Cancer Res._ 15, 1243–1254 (2017). Article  CAS  Google Scholar  * Funato, K., Major, T.,

Lewis, P. W., Allis, C. D. & Tabar, V. Use of human embryonic stem cells to model pediatric gliomas with H3.3K27M histone mutation. _Science_ 346, 1529–1533 (2014). Article  CAS  Google

Scholar  * Wang, Q. et al. Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. _Cancer Cell_ 32, 42–56 (2017). Article

  Google Scholar  * Tam, W. L. & Weinberg, R. A. The epigenetics of epithelial-mesenchymal plasticity in cancer. _Nat. Med._ 19, 1438–1449 (2013). Article  CAS  Google Scholar  *

Osswald, M. et al. Brain tumour cells interconnect to a functional and resistant network. _Nature_ 528, 93–98 (2015). Article  CAS  Google Scholar  * Batlle, E. & Clevers, H. Cancer stem

cells revisited. _Nat. Med._ 23, 1124–1134 (2017). Article  CAS  Google Scholar  * Drokhlyansky, E. et al. The Human and mouse enteric nervous system at single-cell resolution. _Cell_ 182,

1606–1622 (2020). Article  CAS  Google Scholar  * Gojo, J. et al. Single-cell RNA-Seq reveals cellular hierarchies and impaired developmental trajectories in pediatric ependymoma. _Cancer

Cell_ 38, 44–59 (2020). Article  CAS  Google Scholar  * Hovestadt, V. et al. Resolving medulloblastoma cellular architecture by single-cell genomics. _Nature_ 572, 74–79 (2019). Article  CAS

  Google Scholar  * La Manno, G. et al. Molecular diversity of midbrain development in mouse, human, and stem cells. _Cell_ 167, 566–580 (2016). Article  Google Scholar  * Goltsev, Y. et al.

Deep profiling of mouse splenic architecture with CODEX multiplexed imaging. _Cell_ 174, 968–981 (2018). Article  CAS  Google Scholar  * Du, Z. et al. Qualifying antibodies for image-based

immune profiling and multiplexed tissue imaging. _Nat. Protoc._ 14, 2900–2930 (2019). Article  CAS  Google Scholar  * Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from

RNA-Seq data with or without a reference genome. _BMC Bioinf._ 12, 323 (2011). Article  CAS  Google Scholar  * Ding, J. et al. Systematic comparison of single-cell and single-nucleus

RNA-sequencing methods. _Nat. Biotechnol._ 38, 737–746 (2020). Article  CAS  Google Scholar  * Stuart, T. et al. Comprehensive integration of single-cell data. _Cell_ 177, 1888–1902 (2019).

Article  CAS  Google Scholar  * Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. _Nat. Methods_ 16, 1289–1296 (2019). Article  CAS  Google

Scholar  * Fan, J. et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. _Nat. Methods_ 13, 241–244 (2016). Article  CAS  Google Scholar 

* Buttner, M., Ostner, J., Muller, C. L., Theis, F. J. & Schubert, B. scCODA is a Bayesian model for compositional single-cell data analysis. _Nat. Commun._ 12, 6876 (2021). Article  CAS

  Google Scholar  * Stuart, T., Srivastava, A., Madad, S., Lareau, C. A. & Satija, R. Single-cell chromatin state analysis with Signac. _Nat. Methods_ 18, 1333–1341 (2021). Article  CAS

  Google Scholar  * Aibar, S. et al. SCENIC: single-cell regulatory network inference and clustering. _Nat. Methods_ 14, 1083–1086 (2017). Article  CAS  Google Scholar  * Axelrod, S. et al.

Starfish: scalable pipelines for image-based transcriptomics. _J. Open Source Softw._ 6, 2440 (2021). Article  Google Scholar  * Qian, X. et al. Probabilistic cell typing enables fine

mapping of closely related cell types in situ. _Nat. Methods_ 17, 101–106 (2020). Article  CAS  Google Scholar  * Palla, G. et al. Squidpy: a scalable framework for spatial omics analysis.

_Nat. Methods_ 19, 171–178 (2022). Article  CAS  Google Scholar  * Dubois, F. P. B. et al. Structural variants shape driver combinations and outcomes in pediatric high-grade glioma. _Nat.

Cancer_ 3, 994–1011 (2022). Article  CAS  Google Scholar  * Liu, I. et al. The landscape of tumor cell states and spatial organization in H3-K27M mutant diffuse midline glioma across age and

location. Custom scripts v1.0.0. https://doi.org/10.5281/zenodo.7073167 (2022). Download references ACKNOWLEDGEMENTS This work was supported by generous funding from the Hope/Care project

NIH CCSG cancer center (grant P30CA124435 to M.G.F., M.M., A.D., A.R., W.K.A.Y. and M.L. Suvà), the Sajni Fund (M.G.F.), the Claudia Adams Barr Program in Innovative Cancer Research (DFCI)

(M.G.F.), the Cuming Family Fund for Pediatric Brain Tumor Research (M.G.F.), Andruzzi Foundation (M.G.F.), the Anita, Sophia and Athena Fund to Advance DIPG Research and Care (M.G.F.),

Prabal Chakrabarti & Vanessa Ruget (M.G.F.), Hyundai Hope on Wheels (M.G.F.), Liv Like A Unicorn (M.G.F.), Alex’s Lemonade Stand Foundation Crazy 8 Initiative (M.G.F., M.M.) and Solving

Kids’ Cancer/The Bibi Fund (M.G.F.). M.G.F. holds an NIH director’s New Innovator (award DP2NS127705), a Career Award for Medical Scientist from the Burroughs Wellcome Fund, the

Distinguished Scientist Award from the Sontag Foundation and the A-Award from the Alex’s Lemonade Stand Foundation. M.G.F. was also supported by National Cancer Institute SPORE (grant

2P50CA165962). M.N. received funding from the Knut and Alice Wallenberg Foundation (KAW 2018.0172), the Erling Persson Foundation, the Chan Zuckerberg Initiative (SVCF 2017-173964),

Cancerfonden (CAN 2018/604), EU H2020 Marie Skłodowska-Curie Actions project AiPBAND (grant agreement 764281) and the Swedish Research Council (2019-01238). M.M. was supported by the Swifty

Foundation, McKenna Claire Foundation, NIH Director’s Pioneer Award (DP1NS111132 to M.M.), National Cancer Institute (P50CA165962, R01CA258384 and U19CA264504), Robert J. Kleberg, Jr. and

Helen C. Kleberg Foundation (to M.M.) and Cancer Research UK (to M.M.). I.L. was supported by the German Research Foundation (DFG, LI-3486/1-1). B.E. was supported by the Erwin Schrödinger

Fellowship of the Austrian Science Fund (J-4311, B.E.). P.P. was supported by the Ministry of Health of the Czech Republic (grant NU20-03-00240). O.S. received funding from the project

National Institute for Cancer Research (Programme EXCELES, Project ID LX22NPO5102)—Funded by the European Union—Next Generation EU. K.L.L. received funding support from NCI (R01 CA219943 and

P50CA165962). M.D.D. received funding from the Australian National Health and Medical Research Council (NHMRC), RUN DIPG, Tour de Cure, and Kiriwina Investments. The work was further

supported by the ‘Verein unser Kind’ (J.G.). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. We thank the CHLA

Pediatric Research Biorepository supported by the USC Norris Comprehensive Cancer Center (P30 CA014089) for providing tissue resources. We thank Angela Halfmann for assistance with FACS, and

the Molecular Pathology Core Laboratory at Dana-Farber Cancer Institute for help with tissue sectioning. AUTHOR INFORMATION Author notes * These authors contributed equally: Ilon Liu, Li

Jiang, Erik R. Samuelsson. * These authors jointly supervised this work: Michelle Monje, Mats Nilsson, Mariella G. Filbin. AUTHORS AND AFFILIATIONS * Department of Pediatric Oncology,

Dana-Farber Boston Children’s Cancer and Blood Disorders Center, Boston, MA, USA Ilon Liu, Li Jiang, Olivia A. Hack, Daeun Jeong, McKenzie L. Shaw, Bernhard Englinger, Jenna LaBelle, Hafsa

M. Mire, Maria Trissal, Eshini Panditharatna, Johannes Gojo & Mariella G. Filbin * Broad Institute of MIT and Harvard, Cambridge, MA, USA Ilon Liu, Li Jiang, Olivia A. Hack, Daeun Jeong,

 McKenzie L. Shaw, Bernhard Englinger, Jenna LaBelle, Hafsa M. Mire, Maria Trissal, Eshini Panditharatna, Mario L. Suvà, Keith L. Ligon & Mariella G. Filbin * Science for Life

Laboratory, Department of Biochemistry and Biophysics, Stockholm University, Stockholm, Sweden Erik R. Samuelsson, Sergio Marco Salas, Jessica Svedlund & Mats Nilsson * Center for

Neuropathology, Ludwig-Maximilians-University, Munich, Germany Alexander Beck * Department of Urology, Comprehensive Cancer Center, Medical University of Vienna, Vienna, Austria Bernhard

Englinger * Department of Pediatrics and Adolescent Medicine, Comprehensive Center for Pediatrics and Comprehensive Cancer Center, Medical University of Vienna, Vienna, Austria Sibylle

Madlener, Lisa Mayr, Rene Geyeregger, Irene Slavc & Johannes Gojo * Department of Neurology and Neurological Sciences, Stanford University School of Medicine, Stanford, CA, USA Michael

A. Quezada & Michelle Monje * Hopp Children’s Cancer Center Heidelberg (KiTZ), Division of Pediatric Glioma Research, German Cancer Research Center (DKFZ), Heidelberg, Germany Kati J.

Ernst & David T. W. Jones * Department of Oncologic Pathology, Dana-Farber Cancer Institute, Boston, MA, USA Jayne Vogelzang & Keith L. Ligon * Department of Neurological Surgery,

University of Pittsburgh School of Medicine, Pittsburgh, PA, USA Taylor A. Gatesman, Matthew E. Halbert & Sameer Agnihotri * John G. Rangos Sr. Research Center, Children’s Hospital of

Pittsburgh, Pittsburgh, PA, USA Taylor A. Gatesman, Matthew E. Halbert & Sameer Agnihotri * Central European Institute of Technology, Masaryk University, Brno, Czech Republic Hana

Palova, Petra Pokorna & Ondrej Slaby * Pediatric Oncology Department, University Hospital Brno, Faculty of Medicine, Masaryk University, ICRC, Brno, Czech Republic Jaroslav Sterba *

Department of Biology, Faculty of Medicine, Masaryk University, Brno, Czech Republic Ondrej Slaby * Department of Clinical Cell Biology and FACS Core Unit, St. Anna Children’s Cancer

Research Institute (CCRI), Vienna, Austria Rene Geyeregger * Department of Neurological Surgery, University of California San Francisco, San Francisco, CA, USA Aaron Diaz * Cancer Signalling

Research Group, School of Biomedical Sciences and Pharmacy, College of Health, Medicine and Wellbeing, University of Newcastle, Callaghan, New South Wales, Australia Izac J. Findlay & 

Matthew D. Dun * Precision Medicine Program, Hunter Medical Research Institute, New Lambton Heights, New South Wales, Australia Izac J. Findlay & Matthew D. Dun * Center for Data Driven

Discovery in Biomedicine, Children’s Hospital of Philadelphia, Philadelphia, PA, USA Adam Resnick * Department of Pathology, Center for Cancer Research, Massachusetts General Hospital,

Boston, MA, USA Mario L. Suvà * Division of Pediatric Hematology/Oncology, Department of Pediatrics, Michigan Medicine, Ann Arbor, MI, USA Carl Koschmann * Division of Neuropathology and

Neurochemistry, Department of Neurology, Medical University of Vienna, Vienna, Austria Christine Haberler * Department of Neurosurgery, Medical University of Vienna, Vienna, Austria Thomas

Czech * Department of Pathology and Laboratory Medicine, Children’s Hospital Los Angeles, Keck School of Medicine of University of Southern California, Los Angeles, CA, USA Jennifer A.

Cotter * Department of Pathology, Brigham and Women’s Hospital, Boston, MA, USA Keith L. Ligon * Department of Pathology, Boston Children’s Hospital, Boston, MA, USA Keith L. Ligon & 

Sanda Alexandrescu * Department of Neuro-Oncology, Brain Tumor Center, The University of Texas MD Anderson Cancer Center, Houston, TX, USA W. K. Alfred Yung * Massachusetts General Hospital,

Cancer Center, Boston, MA, USA Isabel Arrillaga-Romany * Howard Hughes Medical Institute, Stanford, CA, USA Michelle Monje Authors * Ilon Liu View author publications You can also search

for this author inPubMed Google Scholar * Li Jiang View author publications You can also search for this author inPubMed Google Scholar * Erik R. Samuelsson View author publications You can

also search for this author inPubMed Google Scholar * Sergio Marco Salas View author publications You can also search for this author inPubMed Google Scholar * Alexander Beck View author

publications You can also search for this author inPubMed Google Scholar * Olivia A. Hack View author publications You can also search for this author inPubMed Google Scholar * Daeun Jeong

View author publications You can also search for this author inPubMed Google Scholar * McKenzie L. Shaw View author publications You can also search for this author inPubMed Google Scholar *

Bernhard Englinger View author publications You can also search for this author inPubMed Google Scholar * Jenna LaBelle View author publications You can also search for this author inPubMed

 Google Scholar * Hafsa M. Mire View author publications You can also search for this author inPubMed Google Scholar * Sibylle Madlener View author publications You can also search for this

author inPubMed Google Scholar * Lisa Mayr View author publications You can also search for this author inPubMed Google Scholar * Michael A. Quezada View author publications You can also

search for this author inPubMed Google Scholar * Maria Trissal View author publications You can also search for this author inPubMed Google Scholar * Eshini Panditharatna View author

publications You can also search for this author inPubMed Google Scholar * Kati J. Ernst View author publications You can also search for this author inPubMed Google Scholar * Jayne

Vogelzang View author publications You can also search for this author inPubMed Google Scholar * Taylor A. Gatesman View author publications You can also search for this author inPubMed 

Google Scholar * Matthew E. Halbert View author publications You can also search for this author inPubMed Google Scholar * Hana Palova View author publications You can also search for this

author inPubMed Google Scholar * Petra Pokorna View author publications You can also search for this author inPubMed Google Scholar * Jaroslav Sterba View author publications You can also

search for this author inPubMed Google Scholar * Ondrej Slaby View author publications You can also search for this author inPubMed Google Scholar * Rene Geyeregger View author publications

You can also search for this author inPubMed Google Scholar * Aaron Diaz View author publications You can also search for this author inPubMed Google Scholar * Izac J. Findlay View author

publications You can also search for this author inPubMed Google Scholar * Matthew D. Dun View author publications You can also search for this author inPubMed Google Scholar * Adam Resnick

View author publications You can also search for this author inPubMed Google Scholar * Mario L. Suvà View author publications You can also search for this author inPubMed Google Scholar *

David T. W. Jones View author publications You can also search for this author inPubMed Google Scholar * Sameer Agnihotri View author publications You can also search for this author

inPubMed Google Scholar * Jessica Svedlund View author publications You can also search for this author inPubMed Google Scholar * Carl Koschmann View author publications You can also search

for this author inPubMed Google Scholar * Christine Haberler View author publications You can also search for this author inPubMed Google Scholar * Thomas Czech View author publications You

can also search for this author inPubMed Google Scholar * Irene Slavc View author publications You can also search for this author inPubMed Google Scholar * Jennifer A. Cotter View author

publications You can also search for this author inPubMed Google Scholar * Keith L. Ligon View author publications You can also search for this author inPubMed Google Scholar * Sanda

Alexandrescu View author publications You can also search for this author inPubMed Google Scholar * W. K. Alfred Yung View author publications You can also search for this author inPubMed 

Google Scholar * Isabel Arrillaga-Romany View author publications You can also search for this author inPubMed Google Scholar * Johannes Gojo View author publications You can also search for

this author inPubMed Google Scholar * Michelle Monje View author publications You can also search for this author inPubMed Google Scholar * Mats Nilsson View author publications You can

also search for this author inPubMed Google Scholar * Mariella G. Filbin View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS I.L., L.J.,

E.R.S., M.M., M.N. and M.G.F. conceived the study, designed the experiments, interpreted results and wrote the manuscript with the input of all co-authors. I.L., O.A.H., M.L. Shaw, B.E. and

S.M. performed glioma tissue processing and FACS, with contributions from M.T., E.P., K.J.E., T.A.G. and M.E.H. I.L., O.A.H., M.L. Shaw and H.M.M. generated scRNA-seq data. L.J. and I.L.

conducted glioma scRNA-seq analysis. D.J. performed RNAish experiments, with contributions from B.E. I.L. generated snATAC-seq data, which was analyzed by L.J. I.L., M.A.Q., H.P., O.S.,

P.P., I.J.F., M.D.D. and J.S. generated WES data, which was analyzed by J.L., I.J.F., M.D.D. and P.P. E.R.S. generated ISS data, which was analyzed by S.M.S., supervised by J.S. and M.N.

A.B. performed and analyzed CODEX experiments. Primary tissue resources and pathology consultation were provided by J.V., A.D., A.R., M.L. Suvà, D.T.W.J., S.A., C.K., C.H., T.C., I.S.,

J.A.C., K.L.L., S.A., W.K.A.Y., I.AR., J.G. and M.M. M.N. supervised ISS data generation and analysis. M.G.F. supervised all aspects of the study. CORRESPONDING AUTHORS Correspondence to

Ilon Liu or Mariella G. Filbin. ETHICS DECLARATIONS COMPETING INTERESTS M.G.F. is a consultant for Twentyeight-Seven Therapeutics and Blueprint Medicines. M.N. is Scientific Advisor for 10X

Genomics. M.M. is a SAB member for Cygnal Therapeutics. M.L. Suvà is an equity holder, scientific cofounder and advisory board member of Immunitas Therapeutics. K.L.L. is the founder and

equity holder of Travera and receives consulting fees from BMS, Integragen, Rarecyte and research support from Lilly, BMS and Amgen. J.S. is now (but not when contributing to this

manuscript) an employee of 10X Genomics. The remaining authors declare no competing interests. PEER REVIEW PEER REVIEW INFORMATION _Nature Genetics_ thanks Xiao-nan Li and the other,

anonymous, reviewer(s) for their contribution to the peer review of this work. ADDITIONAL INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in

published maps and institutional affiliations. EXTENDED DATA EXTENDED DATA FIG. 1 NON-MALIGNANT CELL POPULATIONS. UMAP projections highlighting non-malignant cell clusters by expression of

canonical markers of (a) Tumor-associated myeloid cells. (b) T cells. (c) Oligodendrocytes. (d) Endothelial cells. EXTENDED DATA FIG. 2 INTRATUMORAL TRANSCRIPTIONAL HETEROGENEITY OF H3-K27M

DMGS. (a) UMAP of fresh tumor cells, highlighting original samples (color legend) after batch effect correction. (b) Pairwise Pearson correlations (color scale) between NMF factors derived

from each fresh tumor sample (x-axis). Highly correlated NMF factors were combined as metaprograms. (c) Pairwise Pearson correlations (color scale) between metaprograms derived from fresh

H3-K27M DMGs, GBM33, IDH-mutant glioma35. (d) Pairwise Pearson correlations (color scale) between metaprograms independently derived from fresh and frozen tumors. (e) Proportions (y-axis) of

projected fresh tumor derived metaprograms (color legend), that were highly correlated to respective frozen metaprograms, and of fresh OPC-like-3, across frozen tumor nuclei (x-axis).

Instead of fresh OPC-like-2, correlated frozen OPC-like-b was scored to minimize technical ariefacts (see methods). Nuclei with scores <0.2 are denoted as ‘score too low’. (f) UMAP of

frozen tumor nuclei after batch effect correction, with color legend depicting annotation based on single-cell scores of all fresh metaprograms and frozen OPC-like-b (see methods). (g)

Proportion of all cells/nuclei assigned as cycling vs. non-cycling (color legend) across metaprograms. (h) UMAP of location matched IDH-mutant midline tumors, highlighting independently

derived metaprograms. (i) Boxplots depicting metaprogram proportions in all tumors compared by adult vs. pediatric age groups, controlled for pontine (left) or thalamic (right) locations

(Thalamic: adult (N = 6), pediatric (N = 8); Pontine: adult (N = 4), pediatric (N = 15)). (j) Boxplots depicting metaprogram proportions in all tumors compared by pontine and thalamic

locations, controlled for pediatric (left) or adult (right) age groups (Adult: thalamic (N = 6), pontine (N = 4); Pediatric: thalamic (N = 8), pontine (N = 15)). In (i) and (j) The median is

marked by the thick line within the boxplot, the first and third quartiles by the upper and lower limits, and the 1.5x interquartile range by the whiskers. *** denotes credible statistical

changes as assessed by a Bayesian scCODA model, with FDR < 0.05, without multiple test correction. EXTENDED DATA FIG. 3 REGION-SPECIFIC STATES OF OPC-LIKE CELLS. (a) Projection of fresh

tumor-derived metaprograms (x-axis) onto scRNA-seq derived normal cell types (y-axis) of the human hippocampus43. Color scale presents expression scores of normal cell signatures in tumor

cells, while symbol sizes depict expression scores of tumor cell signatures in normal cells. Symbol shape denotes Pearson correlation of expressions, with circle denoting r > =0.5, and

square denoting r < 0.5. (b) Projection of fresh tumor-derived metaprograms (x-axis) onto scRNA-seq derived normal cell types (y-axis) of the developing human cortex40. Color scale

presents expression scores of normal cell signatures in tumor cells, while symbol sizes depict expression scores of tumor cell signatures in normal cells. Symbol shape denotes Pearson

correlation of expressions, with circle denoting r > =0.5, and square denoting r < 0. (c) Projection of fresh tumor-derived metaprograms (x-axis) onto scRNA-seq derived normal cell

types (y-axis) of the neonatal mouse cortex42. Color scale presents expression scores of normal cell signatures in tumor cells, while symbol sizes depict expression scores of tumor cell

signatures in normal cells. Symbol shape denotes Pearson correlation of expressions, with circle denoting r > =0.5, and square denoting r < 0. (d) Diffusion map embedding of single

OPC-like subpopulation transcriptomes (left) and pseudotime analysis by Slingshot where the color scale represents the relative pseudotime (right). (e) Heatmap representing Z-scored

expression levels (color scale) of pre-OPC and OPC marker genes (rows) in tumor OPC-like subpopulations ordered along pseudotime (columns). EXTENDED DATA FIG. 4 CHARACTERISTIC CHROMATIN

PROFILES OF H3-K27M DMG CELL POPULATIONS. (a) UMAP of all nuclei profiled by snATAC-seq from 8 samples (color legend). (b) UMAP of all nuclei profiled by snATAC-seq, highlighting assignments

as malignant or different nonmalignant cell types. (c) UMAP of all snATAC-seq derived tumor nuclei, highlighting sample of origin after batch effect correction. (d) Sample level clustering

analyses and de novo cell type annotations (color legends). (e) Dotplot representation of gene activities (color scale) and proportion of nuclei accessible (dot size) in snATAC-seq profiles

of AC-like-alt., AC-like and OPC-like cells (y-axis) for canonical marker genes of AC-like, OPC-like, NPC-like (as identified in Neftel et al., 201933), and glutamatergic (as described to be

enriched in OPC-like cells by Venkatesh et al., 201925) tumor cells. (f) ScRNA-seq derived log transformed expression levels of synapse-associated genes differentially accessible in

AC-like-alt. cells. (g) Cell state annotations of all snATAC-seq tumor nuclei based on scRNA-seq data following canonical correlation (CCA) and label transfer analyses. (h) UMAP of chromatin

accessibility profiles of all OPC-like subpopulations (color legend). (i) & (j) Dotplot representation of gene activities (color scale) and proportion of nuclei accessible (dot size) in

snATAC-seq profiles of different tumor OPC-like subpopulations (x-axis) for top differentially accessible marker genes (i) and TFs (j) derived from studies of normal pre-OPCs and OPCs45.

(k) Venn diagrams depicting the intersection of differentially accessible chromatin sites with CREs that are linked to GPCs for each cell type. p-values calculated from a two-sided

hypergeometric test are shown. EXTENDED DATA FIG. 5 THE MYELOID CELL LANDSCAPE OF H3-K27M DMGS. (a) Boxplot depicting TAM proportions in all tumor and normal cells profiled by scRNA-seq and

grouped by adult and pediatric sample groups across N = 16 biologically independent samples. The median is marked by the thick line within the boxplot, the first and third quartiles by the

upper and lower limits, and the 1.5x interquartile range by the whiskers. (b) Distributions (mean values + /− 2xSEM) of macrophage and microglia proportions within TAMs across N = 16 pontine

and thalamic tumors. (c) & (e) Venn diagram depicting shared and specific OPC-like-to-myeloid (c) and myeloid-to-OPC-like (e) ligand-receptor interactions between different OPC-like

subpopulations. (d) & (f) Ligand-receptor interactions assessed for each OPC-like subpopulation for OPC-like-to-myeloid (d) and myeloid-to-OPC-like (f) interactions. Color scale depicts

probabilities of interaction, while dot size denotes Benjamini-Hochberg (BH)-corrected p-values from a two-sided permutation test. EXTENDED DATA FIG. 6 THE SINGLE-CELL SPATIAL TRANSCRIPTOMIC

ARCHITECTURE OF H3-K27M DMGS. (a) Representative HybISS gene maps for 16 H3-K27M tumors (1 experiment/tumor over the entire image section with 100-20,000 cells profiled/tumor). Scale bar

corresponds to 100 µm in all panels. (b) Confusion matrix of pciSeq derived tumor cell state scores for all samples. The color scale represents the mean probability assigned to a cell when a

specific cell state is predicted. Higher values indicate a more probable prediction. (c) Scatter plot representing numbers of malignant cells assigned to a cell state (color scale) for each

sample (dot), as inferred from pciSeq based on 116 marker genes (y-axis) or on the 4 best markers (x-axis). The Pearson correlation coefficient between both marker sets is shown in red. (d)

Sample-level proportions (x-axis) of malignant and non-malignant cells (color legend) across 16 tumors (y-axis) profiled by HybISS as assessed by anti-H3.3K27M IF. (e) Sample-level

proportions (x-axis) of non-malignant cell types (color legend) assigned by HybISS for the 16 H3-K27M DMGs (y-axis). (f) Scatter plot representing numbers of malignant cells assigned to a

specific cell state (color scale) for each sample profiled (dot), as inferred from pciSeq based on 116 marker genes (y-axis) or on selected IF markers (PDGFRA, BCAS1, GFAP, CD44/CD63)

(x-axis). The Pearson correlation coefficient between both marker sets is shown in red. (g) & (h) Representative multiplexed IF (CODEX) images, showing spatially distinct subpopulations

of malignant (marker: H3-K27M) OPC-like (marker: PDGFRA), OC-like (marker: BCAS1), AC-like (marker: GFAP), and proliferating cells (marker: Ki67) in (g), and of MES-like (marker: CD44/CD63)

and myeloid cells (marker: IBA1) in (h). For each tumor, one experiment was performed with ~70,000-1.2 million individual cells profiled per sample over the entire tissue section. (i)

Neighborhood enrichment analysis between all malignant and non-malignant cell populations, identified at 50 μm. The color scale denotes the probability of finding a cell when a second cell

type is present divided by the probability of finding the second cell type. SUPPLEMENTARY INFORMATION SUPPLEMENTARY INFORMATION Supplementary Figs. 1–3 and Supplementary Note. REPORTING

SUMMARY SUPPLEMENTARY TABLE Supplementary Tables 1–7. RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which

permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to

the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless

indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or

exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints

and permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Liu, I., Jiang, L., Samuelsson, E.R. _et al._ The landscape of tumor cell states and spatial organization in H3-K27M mutant diffuse

midline glioma across age and location. _Nat Genet_ 54, 1881–1894 (2022). https://doi.org/10.1038/s41588-022-01236-3 Download citation * Received: 28 October 2021 * Accepted: 20 October 2022

* Published: 05 December 2022 * Issue Date: December 2022 * DOI: https://doi.org/10.1038/s41588-022-01236-3 SHARE THIS ARTICLE Anyone you share the following link with will be able to read

this content: Get shareable link Sorry, a shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative