Cryptic developmental events determine medulloblastoma radiosensitivity and cellular heterogeneity without altering transcriptomic profile

It is unclear why medulloblastoma patients receiving similar treatments experience different outcomes. Transcriptomic profiling identified subgroups with different prognoses, but in each subgroup, individuals remain at risk of incurable recurrence. To investigate why similar-appearing tumors produce variable outcomes, we analyzed medulloblastomas triggered in transgenic mice by a common driver mutation expressed at different points in brain development. We genetically engineered mice to express oncogenic SmoM2, starting in multipotent glio-neuronal stem cells, or committed neural progenitors. Both groups developed medulloblastomas with similar transcriptomic profiles. We compared medulloblastoma progression, radiosensitivity, and cellular heterogeneity, determined by single-cell transcriptomic analysis (scRNA-seq). Stem cell-triggered medulloblastomas progressed faster, contained more OLIG2-expressing stem-like cells, and consistently showed radioresistance. In contrast, progenitor-triggered MBs progressed slower, down-regulated stem-like cells and were curable with radiation. Progenitor-triggered medulloblastomas also contained more diverse stromal populations, with more Ccr2+ macrophages and fewer Igf1+ microglia, indicating that developmental events affected the subsequent tumor microenvironment. Reduced mTORC1 activity in M-Smo tumors suggests that differential Igf1 contributed to differences in phenotype. Developmental events in tumorigenesis that were obscure in transcriptomic profiles thus remained cryptic determinants of tumor composition and outcome. Precise understanding of medulloblastoma pathogenesis and prognosis requires supplementing transcriptomic/methylomic studies with analyses that resolve cellular heterogeneity. Malawsky, Weir et al. compare medulloblastomas that form in mice engineered to express SmoM2, a constitutively active component of the Sonic hedgehog signaling pathway, starting in either neural progenitor or stem cells. Medulloblastomas in the two models show similar bulk transcriptomic profiles, but contain different fractions of tumor stem cells and different myeloid populations, and demonstrate different therapeutic responses.

M edulloblastoma, the most common malignant pediatric brain tumor, is typically treated with surgery, radiation, and chemotherapy, which are effective in 80-90% of patients. However, individual patients face a significant risk of treatment failure despite uniform treatment, and the causes of treatment failure are incompletely understood. Transcriptomic studies have identified four major subgroups of medulloblastoma: WNT, SHH, group 3, and group 4 1 . Each subgroup has a different prognosis, but within each subgroup, outcomes for individuals are heterogeneous 2 . The factors that determine the variable outcomes for patients with similar-appearing medulloblastomas are unclear. We tested the possibility that different outcomes of medulloblastoma therapy can be determined by developmental events that are cryptic at the time of presentation.
SHH-subgroup medulloblastomas are grouped together because they show similar patterns of gene expression in bulk transcriptomic studies, indicating SHH pathway activation. Despite sharing a common oncogenic pathway, patients with SHH-subgroup medulloblastomas show different responses to treatment, with~20% developing incurable recurrence. It is not clear whether differences in outcome are stochastic or driven by determinants that remain to be identified. Age of onset, however, is a factor that clearly influences prognosis within each subgroup 2 , suggesting a potential effect of developmental timing on tumor behavior.
Prior studies in genetically engineered mice show that cerebellar granule cell progenitors (CGNPs) are proximal cells of origin for SHH-driven medulloblastoma 3,4 . CGNPs are a population of SHH-responsive, committed neural progenitors that derive from the rhombic-lip, migrate to the external granule layer (EGL) of the cerebellum, and then proliferate rapidly in response to SHH ligand secreted by the Purkinje neurons 5,6 . CGNPs proliferate in the first two weeks of life in mice, or the first year of life in humans, to generate the cerebellar granule neurons (CGNs) the largest neuronal population in the brain 7 . Mutations that hyperactivate SHH signaling in CGNPs cause familial medulloblastoma in humans and recapitulate medulloblastoma formation in mice, providing genetically faithful primary tumor models 3,4,8 .
Importantly, CGNPs are not a homogeneous population. While the EGL is predominantly populated by ATOH1 (AKA MATH1)-expressing progenitors, a small subset of NESTIN+/ ATOH1-cells reside in the EGL 9 . These EGL cells are typically quiescent in vivo but proliferate in response to SHH pathway activation and can give rise to SHH-driven medulloblastoma 9 . Moreover, the ATOH1+ cells of the EGL comprise different subsets 10,11 , including a transient subpopulation that expresses the stem cell marker SOX2 and may be particularly vulnerable to SHH-driven tumorigenesis 11 . CGNPs are thus a heterogeneous population with varying oncogenic potential.
The developmental origins of medulloblastoma can be analyzed by pairing different Cre drivers with conditional mutations of the SHH receptor components Ptc and Smo directing SHH hyperactivation to broad lineages that include CGNPs or more narrow lineages within the CGNP population. For example, hGFAP-Cre targets a lineage of neuroglial stem cells throughout the brain that includes both ATOH1+ CGNPs and NESTIN +/ATOH1− cells of the EGL 12,13 . Math1-Cre, in contrast, targets the ATOH1-expressing CGNPs [14][15][16] , including the SOX2+ subset 11 . A prior study compared the effects of deleting Ptc either in neural stem cells in hGFAP-Cre/Ptc loxP/loxP mice or in CGNPs in Math1-Cre/Ptc loxP/loxP mice 17 . Both genotypes developed medulloblastoma with 100% penetrance and in hGFAP-Cre/ Ptc loxP/loxP where SHH was hyperactivated throughout the brain, tumors developed only in the cerebellum. Similarly, inducing a Cre-dependent, constitutively active allele of Smo (SmoM2) in either stem cells or CGNPs, using respectively hGFAP-Cre or Math1-Cre, resulted in medulloblastoma with 100% penetrance and no other brain tumors 18 . These studies show that the EGL population is uniquely competent to undergo SHH-mediated tumorigenesis and that cells with hyperactivation of SHH signaling prior to the formation of the EGL must advance along the CGNP developmental trajectory by migrating to the cerebellar surface before giving rise to tumors 17,18 .
Medulloblastomas initiated by hGFAP-Cre in stem cells progress faster than medulloblastomas initiated by Math1-Cre, producing a shorter EFS despite occurring in the same location and showing similar pathology and gene expression profile 17,18 . Similarly, medulloblastomas initiated by conditional activation of SmoM2 either prenatally, using hGFAP-Cre, or postnatally, using tamoxifeninducible Math1-CreER, are histologically and molecularly indistinguishable, but show different propensities for anchorageindependent growth in vitro 19 . In these studies, SHH-driven medulloblastomas triggered early, in stem cells or later, in progenitor cells, show overall similarities but also specific differences.
The different survival times when tumors are triggered with hGFAP-Cre or Math1-Cre suggest that the timing of the oncogenic event can act as a cryptic factor that produces clinically relevant effects that persist throughout the generations of tumor cells. We tested the clinical relevance of this possibility by comparing the responses to therapy of medulloblastomas initiated in either stem cells or neural progenitors, and by subjecting both types of tumors to scRNA-seq analysis. We show that timing of tumor initiation within the lineage trajectory of GFAP+ stem cells to ATOH1+ progenitors influences the cellular heterogeneity within the resulting tumors, without detectably altering average gene expression profiles, producing tumors that appear similar but contain divergent sub-populations with different tumor-stromal interactions and treatment responses.

Results
Similar pathology and gene expression in medulloblastomas from progenitors or stem cells.  18 , the survival times of untreated G-Smo mice were typically shorter than those of M-Smo mice (Fig. 1b). Moreover, while M-Smo mice showed significantly improved survival after radiation therapy (xRT), consistent with our prior studies 20 , xRT did not extend the survival of G-Smo mice (Fig. 1c, d). In all these studies, 8-12 mice of each genotype were studied. The prognosis of G-Smo and M-Smo mice was therefore markedly different, despite their common oncogenic driver, tumor pathology, and similarity in bulk transcriptomic studies.
Apoptotic response after xRT is less uniform in stem cellderived medulloblastoma. To compare the cellular responses to xRT in G-Smo and M-Smo tumors, we analyzed apoptosis after xRT treatment. We have previously shown that a single fraction of xRT induces widespread, synchronous apoptosis after a 3-h latent period 21 , and that this apoptotic response to xRT is required for efficacy 20 . We treated P15 G-Smo (n = 4) and M-Smo (n = 3) mice with 10 Gy xRT delivered in a single dose and then harvested tumors after 3 hours (n = 3). We identified apoptotic cells using immunohistochemistry (IHC) for cleaved caspase-3 (cC3; Fig. 1e) and compared the fraction of apoptotic cells in each genotype. We found that xRT induced apoptosis throughout the tumors in both genotypes. However, radiated G-Smo tumors showed significantly smaller fractions of apoptotic cells compared to radiated M-Smo tumors. While both tumor types showed large radiation-sensitive populations, cells that survived radiation were consistently more numerous in G-Smo tumors.

Similar average gene expression in G-Smo and M-Smo tumors after xRT.
To determine if the differences in cellular responses to xRT produced detectable differences in the transcriptomic profiles in radiated G-Smo and M-Smo tumors, we analyzed 6 tumors of each genotype harvested 2 h after a single 10 Gy dose of xRT using expression microarrays. We used a single fraction of xRT in order to produce synchronized changes in gene expression, and harvested tumors after 2 h, in order to study the latent period prior to the onset of apoptosis. 74 probes sets detected statistically significant signals in radiated tumors versus untreated tumors, representing 73 identified transcripts altered by xRT (Supplementary Data 3). Consistent with prior studies 20 , xRT induced p53-regulated genes, including Cdkn1a, Trp53inp1, and Bbc (aka PUMA). We used 2-way ANOVA to analyze the interaction of genotype and radiation treatment in the combined transcriptomic data from untreated and radiated G-Smo and M-Smo tumors; we found no genes with differential expression with corrected p value of <0.05. Bulk transcriptomic analysis of treated and untreated G-Smo and M-Smo tumors thus identified genes induced by xRT, but did not identify genotype-specific differences in transcriptomic response.
scRNA-seq identifies the difference between medulloblastomas from progenitors or stem cells. To compare differences between G-Smo and M-Smo tumors with cellular resolution, we analyzed both tumor types using scRNA-seq. Transcriptomic analysis at single-cell resolution allowed us to examine tumor subpopulations that might be overlooked in bulk transcriptomic studies. We harvested medulloblastomas from 5 M-Smo and 6 G-Smo at P15, dissociated the tumors, and subjected the cells to bead-based Drop-seq analysis, as previously described 10 . Putative cells identified by bead-specific barcodes were subjected to exclusion criteria described in the "Methods" section, to address the common problems of gene drop out, unintentional cell-cell multiplexing, and premature cell lysis 22,23 10 . We identified 17 principal components that described >78% of the variance and used UMAP to place cells in a 2-dimensional space according to their distances in PCA space, with Louvain clusters color-coded (Fig. 2a). We noted that cells in the same clusters localized close together in the UMAP, supporting the validity of the cluster assignments. To determine the biological relevance of the clusters, we generated cluster-specific differential gene expression profiles (Supplementary Data 4); we compared for each gene the expression by cells within each cluster to the expression by all cells outside the cluster using Wilcoxon rank-sum test. We then used cluster-specific gene expression patterns to infer the type of cells represented by each cluster.
Using these methods, we identified 8 clusters as different types of stromal cells typical of brain tissue, including astrocytes, oligodendrocytes, macrophages/microglia, vascular cells, fibroblasts, and ciliated cells resembling ependymal or choroid plexus cells (Table 1 and Fig. 2b). These 8 clusters localized as discrete single-cluster units on the UMAP projection. The other 14 clusters localized to a multi-cluster complex in which each cluster shared a border with other clusters. We identified these 14 clusters as tumor cells in a range of states that paralleled CGNP development, from proliferative cells at different phases of the cell cycle to non-proliferative cells at different stages of neural differentiation (Table 1). We identified proliferative clusters by expression of proliferation marker Mki67, Cyclin expression, and SHH transcription factor Gli1, and further characterized proliferative cells as quiescent, cycling or M-phase enriched based on cluster-specific gene expression ( Table 1). The nonproliferative clusters showed successive expression of early to later differentiation markers Barhl1, Cntn2, Rbfox3, and Grin2b (Table 1 and Fig. 2c), as in our prior study of M-Smo tumors 10 . We included CGNs as the most differentiated cell type within this group.
To compare the populations within M-Smo and G-Smo tumors, we deconvoluted the UMAP by genotype (Fig. 2d). For quantitative comparison, we determined the number of cells from each replicate animal in each cluster, normalized to the total number of cells from that animal, and then compared the cluster populations from replicate M-Smo and G-Smo mice (Fig. 2d, e). . We excluded stromal cells in order to prevent stromal gene expression patterns from obscuring differences in tumor cell gene expression. We noted that Clusters 1, 2, and 7 showed increased expression of genes associated with stem cell phenotype, including Nes, Vim, Olig1, and Olig2. Feature plots of each of these genes confirmed increased expression in G-Smo tumors, particularly in the region of Clusters 1, 2, and 7 ( Fig. 2f). We selected Olig2 for further study because recent functional genetic studies have shown that Olig2+ tumor cells in medulloblastoma are cancer stem cells that play an important role in medulloblastoma initiation and recurrence 25 .

Different temporal patterns of stem cell behavior in M-Smo
and G-Smo tumors. To confirm the differential expression of Olig2 at the protein level and to compare the temporal course of Olig2 expression efficiently between genotypes, we labeled tumor sections using IHC. Our scRNA-seq data showed that both oligodendrocytes and tumor stem cells express Olig2, and that oligodendrocytes can be distinguished from stem cells by the expression of Sox10 (Fig. 2b, f). Labeling tumor sections with OLIG2 and SOX10 antibodies (Fig. 3a) demonstrated both OLIG2+/SOX10+ cells that we considered to be oligodendrocytes and OLIG2+/SOX10− cells that we considered to be OLIG2-expressing tumor stem cells, equivalent to the Olig2+ cells of clusters 1, 2, and 7. The OLIG2/SOX10 assay allowed us to differentiate tumor stem cells from oligodendrocytes and to compare tumor stem cell populations at different ages and treatment conditions (Fig. 3b) Quantitative analysis confirmed increased OLIG2+ stem cells in P15 G-Smo tumors compared to M-Smo tumors (Fig. 3c). Differences in OLIG2+ stem cell populations could be attributed to differences in stem cell proliferation, or to differences in the tendency to maintain a stem cell phenotype. We compared proliferation in the stem cell populations of P15 G-Smo and M-Smo tumors to determine whether the production of stem cells differed between genotypes. We measured stem cell proliferation using IHC for OLIG2 and the proliferation marker phosphorylated RB (pRB; Supplementary Fig. 1a). Quantitative analysis showed no significant difference in proliferation rate in the Prior studies found that OLIG2+ tumor stem cells were more numerous in the early stages of medulloblastoma tumorigenesis and decreased as tumors enlarged 18 . We, therefore, compared stem cell dynamics in G-Smo and M-Smo tumors at P5, early in tumor formation, and at P15 when tumors were larger, using 3-5 replicates of each genotype at each age. We found that at P5, unlike P15 OLIG2-expressing stem cells have been shown to be treatmentresistant in mouse medulloblastoma models, and to drive recurrence 25 (Fig. 4a). Cells within the fibroblast and macrophage/ microglia clusters did not express Yfp, and the rare Yfp+ cells in the endothelial and ependymal clusters did not indicate a significant trend as they were not observed in more than one replicate of either genotype. Glial cells were therefore the only stromal cell types that derived from the tumor lineage in either stem cell-derived or progenitor-derived tumors.
We compared the expression of specific neural markers to assess differences in tumor cell fates. Both genotypes consisted predominantly of Barhl1+ cells. We detected a significant genotype-specific difference within the Atoh1-lineage, with a higher proportion of Barhl1+ cells in G-Smo tumors (p < 1.0 × 10 −15 , two-proportions z-test) and a  (Fig. 4b). However, we did not detect statistically significant differences in the Ascl1+, Pax+, and Pax2+ GABAergic interneuron-lineage cells derived from the ventricular zone (Fig. 4c). Therefore, although hGFAP-Cre activated   Fig. 2 and Supplementary Data 6). We identified each stromal cell type based on gene expression ( Fig. 5a and Supplementary Data 6), and isolated the endothelial, macrophage/microglial, and fibroblastic populations. We then subjected the cells of each isolated cell type to a new PCA to sub-cluster each cell type.
Endothelial cells show cancer-specific changes without reflecting developmental differences between tumors. Endothelial cells showed significant differences between WT and tumor but did not show statistically significant differences between G-Smo and M-Smo tumors. The unsupervised analysis defined 2 clusters (Fig. 5b). Cluster E 0 included cells from WT cerebella and both tumor genotype, while E 1 was populated predominantly by cells from the tumors, with similar proportions from G-Smo and M-Smo (Fig. 5c, d). Both E 0 and E 1 expressed the endothelial markers Pecam1 and Cldn5, confirming endothelial identity (Fig. 5e).
Each cluster showed cluster-specific gene expression (Fig. 5f, g and Supplementary Data 7). The tumor-specific E 1 cells, showed increased expression of genes likely to contribute to malignancy, including the VEGF receptor Flt1, the p-Glycoprotein Abcb1a (aka Mdr1), and the CXCR4 ligand Cxcl12 (aka Sdf1), which has been shown to promote medulloblastoma growth and glioblastoma-endothelial interactions [31][32][33][34] . Biologically relevant differences between endothelial populations in WT cerebella and tumors were shared between tumor genotypes, consistent with their overall similarity.  (Fig. 6a). Projection of C1qb expression confirmed that clusters M 0-M 3 were populated by myeloid cells (Fig. 6b). In contrast, cluster M 4, which was the least populated, was C1qb-and expressed Cnn3 and Meis1 (Fig. 6b); this marker pattern identified cluster M 4 as choroid plexus epithelial cells that have been noted to cluster with myeloid cells in other scRNA-seq analyses 35 . We identified the types of myeloid cells in each cluster by defining the sets of genes upregulated by cells within the cluster, compared to cells in the other 4 clusters (Supplementary Data 8). We projected key cluster-specific marker genes on a UMAP of the 4 clusters, along with other known markers of myeloid phenotype (Fig. 6c-f). Clusters M 0 and M 1, which localized to clusters on one side of the UMAP, expressed Cx3cr1 which distinguished them as microglia, while clusters M 2 and M 3, opposite in the UMAP, showed minimal Cx3cr1, indicating that they were macrophages (Fig. 6c). M 0 cells showed cluster-specific expression of Sparc (Fig. 6c), identifying them as mature, ramified microglia 36 . M 1 microglia specifically expressed Mrc1 and Wfdc17 (Fig. 6d) which have been linked to an M2-like, anti-inflammatory phenotype [37][38][39] . Both M 0 and M 1 cells expressed Igf1 (Fig. 6d), which has been shown to be a paracrine signal promoting growth in SHH medulloblastoma 28 ; Igf1was not detected in M 2 and M 3 cells. M 2 macrophages specifically expressed MHCII components, including Histocompatibility 2, Class II Antigen E alpha (H2-Ea), and the Invariant Polypeptide of Major Histocompatibility Complex, Class II Antigen-associated (Cd74) as well as Il1b and Ccr2 (Fig. 6e), consistent with a pro-inflammatory M1-like phenotype [40][41][42] . Importantly, Ccr2+ macrophages have previously been shown to exert an anti-tumor effect in a medulloblastoma 42 . M 3 macrophages specifically expressed Cd163 and Mrc1 (Fig. 6f), also consistent with an antiinflammatory M2 phenotype 43,44 . Myeloid cells thus resolved into microglial and macrophage populations, each with M1-like and M2-like subsets ( Table 2). M 0 and M 1 cells expressed Igf1, which has been shown to promote medulloblastoma progression, while M 2 expressed Ccr2, which is associated with tumorinhibiting macrophages.
Each cluster distributed differently across G-Smo and M-Smo tumors and WT cerebella (Fig. 6g, h). The M1-like microglial

Differential cytokine expression in G-Smo and M-Smo tumors.
To consider potential mechanisms for the differences in macrophage populations in the two tumor types, we performed a nonbiased comparison of cytokine expression. To generate a list of known cytokines and chemokines, we used the set of 232 genes tagged with the Gene Ontology term "Cytokine Activity". For each gene, we conducted differential expression testing between  (Fig. 6i). MIF is a ligand for CD74 45 , and increased Mif in G-Smo tumors may contribute to the markedly lower Cd74+ population. Prior studies in human glioblastoma and melanoma associate MIF with cancer stem cells show that intercellular communication through MIF-CD74 is immunosuppressive, and that blocking MIF-CD74 signaling increases tumor-associated M1 macrophages [46][47][48] . Based on these prior studies and the inverse correlation of MIF and CD74 in our tumor model, we propose that MIF functions in medulloblastoma to bias myeloid cells toward an Igf1+ phenotype, and acts more effectively in G-Smo tumors, which have higher Mif expression.  To determine if this difference in Igf1 may be biologically significant, we analyzed mTORC1 activity, which is increased by IGF1 signaling and known to be important in SHH medulloblastoma progression 49 Differential myeloid marker expression in mouse and human SHH medulloblastomas. We confirmed differential myeloid populations using IHC for the MHCII glycoprotein coded by H2-Ea and the pan macrophage/microglial marker IBA1, comparing 4 G-Smo and 3 M-Smo tumors (Fig. 7a). H2-EA+ macrophages were more prevalent in M-Smo tumors (Fig. 7b), consistent with the scRNA-seq data. These results confirm that the protein expression of these markers matched the transcript data from scRNA-seq and demonstrate that immunohistochemical staining for MHCII components was sufficient to distinguish between G-Smo and M-Smo tumors.
To determine whether similar differences in myeloid markers can be used to probe clinical medulloblastoma samples, we used IHC to detect HLA-DR proteins, the human orthologs of mouse H2-EA, as described previously 51 . We analyzed 30 medulloblastoma samples for which SHH, WNT, group 3, and group 4 subgroups and subtypes within each subgroup had been using bulk transcriptomic and methylomic analysis according to published criteria 2,52 . All slides were subjected to blinded review and the frequency of HLA-DR+ cells was scored on a scale of 0-3 (Fig. 7c) 51 . 10 of the medulloblastomas were SHH-subgroup tumors, and these tumors received significantly higher HLA-DR scores compared to tumors of the other types (Fig. 7d). Subtype was determinable for 9/10 SHH medulloblastomas, including 2 SHH-alpha tumors, 4 SHH beta tumors, and 3 SHH delta tumors. HLA-DR staining varied between SHH subtypes, with SHH-alpha tumors showing markedly less HLA-DR compared to beta and delta tumors (Fig. 7e). SHH beta and SHH delta are the infantpredominant subtypes 2 , and the trend toward higher HLA-DR expression in these subtypes suggests that a developmental process affects myeloid subgroups in human tumors, as in our mouse models.
The trends noted in our IHC analysis correlate well with bulk transcriptomic data from prior studies. In the published data first used to establish the subtypes within the 4 medulloblastoma subgroups 2 , SHH-subgroup tumors showed significantly higher HLA-DRA mRNA (Fig. 7f), and SHH-alpha subtype tumors showed significantly lower HLA-DRA mRNA, compared to the other SHH subtypes (Fig. 7g). We also examined SHH-subgroup tumors stratified by age, irrespective of subtype; infant medulloblastomas showed significantly higher HLA-DRA (Fig. 7h), consistent with developmental differences producing differences in the TME, as in our model.
Developmental differences influence tumor fibroblast populations. Genotype-specific differences in the TME were not limited to myeloid cells; we also noted differences in the fibroblast populations. Fibroblasts in G-Smo and M-Smo tumors and P7 WT cerebella assorted into 3 clusters, with significant differences in the types of fibroblasts in each genotype (Fig. 8a-c). Each cluster showed specific gene expression patterns (Fig. 8d-f and Supplementary Data 9), with genotype-specific effects retinoid signaling genes; Fabp5 was WT-specific, and Rbp4 and Crabp2 were specifically down-regulated in G-Smo tumors (Fig. 8g). These differential patterns show that the timing of oncogenesis affects the fibroblastic stroma of the TME as well as the tumorassociated myeloid cells.  enforcing effects; thus, differences in Igf1-expressing myeloid populations correlate with differences in mTORC1 activation in tumor cells Differences in developmental history, therefore did not produce large differences in global gene expression profile, but significantly altered the heterogeneity within tumors, producing differences that can be self-amplifying. Our analysis identified multiple processes that are likely to contribute mechanistically to the differences in the progression and recurrence of G-Smo and M-Smo tumors. Olig2-expressing tumor stem cells, Ccr2+ macrophages, and Igf1+ microglial have all been shown to affect tumor progression and prognosis 25,28,42 . The importance of Olig2+ medulloblastoma stem cells was shown in prior studies where targeting the Olig2+ population, either by conditional ablation of Olig2-expressing cells using HSV TK or conditional genetic deletion of the Olig2 locus, reduced the growth of SHH-driven medulloblastomas in mice 25 . We propose that by preventing p53-dependent apoptosis after xRT, pOLIG2 in tumor stem cells may allow survival long enough for DNA repair, as seen in Bax-deleted medulloblastomas after xRT 20 .

Discussion
Differences in Igf1+ and Ccr2+ myeloid cells may also contribute to poor outcomes in G-Smo mice. Ccr2-expressing macrophages, which are reduced in G-Smo tumors, have been shown to suppress tumor growth 42 . In contrast, Igf1-secreting microglia, which are increased in G-Smo tumors, support medulloblastoma progression 28 . The differences in the polarization of myeloid cells in G-Smo and M-Smo tumors, along with differences in Olig2+ stem cells, are multiple mechanisms that are each sufficient to worsen prognosis, and these mechanisms may act in combination.
Several factors may obscure important features in bulk transcriptomic studies. Transcriptomic signals from small subsets of cells may not be detectable when averaged with larger subsets of cells in bulk tumor lysates. For example, Ccr2+ cells did not produce a detectable signal in our microarray study but were detectable by scRNA-seq. In addition, cell types that are very common, such as proliferating cells, may generate large signals that do not produce statistically detectable variation, as we found that proliferation markers were not statistically different in bulk In contrast, studies with cellular resolution, including scRNA-seq and IHC, allow the detection of both small populations, the counting of proliferative cells identified by multiple markers, and the differentiation of markers by their cellular context. While bulk transcriptomic analysis of HLA-DRA mRNA showed trends that matched our HLA-DR IHC study, the variability between samples in each subtype of the bulk analysis was large; statistical trends could be discerned between subtypes, but each subtype included individuals within a wide range of values. Larger studies are needed to determine if cell-resolved data for markers, including IHC for HLA-DR or OLIG2, may allow stratification into more homogeneous groups than the SHH subtypes identified by bulk methods.
The Our data show that developmental events that would be cryptic in a clinical setting can influence the clinical outcomes by affecting both tumor cells and stroma. While these effects passed undetected in bulk transcriptomic studies, the fractions of OLIG2+ tumor cells and MHCII-expressing macrophages were effective biomarkers that distinguished radioresistant stem cell-derived tumors from radiosensitive tumors originating in CGNPs. These biomarkers succeeded because they were sensitive to differences in cellular heterogeneity. Similarly, Olig2+ populations and myeloid subtypes are readily measurable in clinical samples. While scRNAseq identified differences between tumors, each of these differences can also be probed in clinical samples using IHC or flow cytometry, which preserves cellular information that is lost in transcriptomic analysis. Future studies are needed to determine whether analysis of these parameters provides prognostic information that reduces the heterogeneity of transcriptomic subgroups, improving prognostication and precision therapy. Human subjects. Medulloblastoma samples were obtained from patients consented under the University of Colorado IRB COMIRB 95-500. All patient materials were de-identified prior to their use in this study. Histology, IHC, and western blot. For histology and IHC, mouse brains were processed, immunostained, and quantitatively analyzed as previously described 10,20,53 . For western blot, samples were snap-frozen and lysed, blotted, and stained as previously described 53 . Primary antibodies used were: cC3 1:100 (Biocare, #229), pRB diluted x (x), OLIG2 diluted 1:100 (Cell Marque, # 387R-14), pOLIG2 diluted 1:3000 (x), SOX10 diluted 1:200 (Cell Signaling Technology, #7833S), MHC Class II glycoprotein H2-Ea diluted 1:200 (Novus, # NBP1-43312), IBA1 1:2000 (Wako Chemicals, #019-19741) and HLA-DR diluted 1:50 (Wako NBP2-47670), and p4EBP1 diluted 1:500. HLA-DR studies were performed as previously described 51 . All other stained slides were counterstained with DAPI, digitally imaged using an Aperio Scan Scope XT (Aperio), and subjected to automated cell counting using Tissue Studio (Definiens).
Radiation therapy and survival studies. Medulloblastoma-bearing mice were treated with 10 Gy X-ray irradiation, delivered as 5 fractions of 2 Gy each, as previously described 20 . Briefly, starting P10, G-Smo, and M-Smo mice were irradiated daily for 5 days. Irradiation was performed under general anesthesia with isoflurane, delivered by vaporizer through nose cones, after which mice were allowed to recover and then returned to their dams. Following radiation therapy, radiation-treated mice and untreated littermate controls were observed for symptoms of tumor progression, including movement disorder, ataxia, or sustained weight loss. Mice showing symptoms of progression were euthanized, and the time to progression was considered to be the EFS. Brain pathology was analyzed to confirm tumor progression. Tissue preparation for Drop-seq. Mice were anesthetized using isoflurane and then euthanized via decapitation. The brain was divided along the sagittal midline and one half was processed for histology while a large sample of tumor was dissected from the other half and processed for Drop-seq analysis. This sample was dissociated using the Papain Dissociation System (Worthington Biochemical) following the protocol used in previous studies 10,54 . Briefly, tumor samples were incubated in papain at 37°C for 15 min, then triturated, and the suspended cells were spun through a density gradient of ovomucoid inhibitor.
Pelleted cells were then resuspended in 1 mL HBSS with 6 g/L glucose and diluted in PBS-BSA solution to a concentration of 95-110 cells/µL. Barcoded Seq B Drop-seq beads (ChemGenes) were diluted in Drop-seq lysis buffer to a concentration between 95-110 beads/µL. Tumor cells were co-encapsulated with barcoded beads using FlowJEM brand PDMS devices as previously described 10 . All cells were processed within one hour of tissue dissociation. Droplet breakage and library preparation steps followed Drop-seq protocol V3.1 55 . After PCR, amplified cDNA was subjected Ampure XP cleanup at 0.6× and 1× ratios to eliminate residual PCR primers and debris. found by the bioanalyzer electropherogram. 1.If PCR failed to generate adequate cDNA, the PCR was repeated with the 3rd round increased from 11 to 13 cycles.
For QC purposes, library pools consisting of the tagmented cDNA from 2000 beads/run were prepared and sequenced to low depth (~2.5 M reads/2 K beads). We used the resulting data to assess library efficiency, including total read losses to PolyA regions, nonsense barcodes, and adapter sequences as well as the quality and number of the transcriptomes captured. Passable runs contained 40-60% of reads associated with the top 80-100 barcodes found in 2000 beads.
Drop-seq runs passing QC were then prepared for high-depth sequencing on an Illumina Hi-Seq 4000. Each sample underwent a new generation of bulk cDNA from the stored beads and was prepared with the same ratios as described above. Pools were formulated according to the number of cells/sample to avoid oversampling of each sample and to balance the reads per lane across the Hi-Seq.
Data analysis was performed using the Seurat R package version 3.1.1 56 . Data were subjected to several filtering steps. First, only genes that were detected in at least 30 cells were considered, to prevent misaligned reads appearing as rare transcripts in the data. Cells were then filtered using specific QC criteria to limit the analysis to cells with transcriptomes that were well-characterized and not apoptotic.
We noted that G-Smo cells were sequenced at a greater depth than M-Smo cells which can introduce unwanted batch effects into the analysis. Consistent with best practices 24 , we downsampled the G-Smo cells to 60% of their original depth so as to achieve similar sequencing depth between G-Smo and M-Smo cells prior to further filtering.
Putative cells with fewer than 500 detected RNA molecules (nCount) or 200 different genes (nFeature) were considered to have too little information to be useful, and potentially to contain mostly ambient mRNA reads. Putative cells with >4 standard deviations above the median nCount or nFeature were suspected to be doublets, improperly merged barcodes, or sequencing artifacts and were excluded. As in our previously published work, putative cells with more than 10% mitochondrial transcripts were suspected to be dying cells and also excluded 10 .
In total, 53% of putative cells from G- scRNA-seq Data normalization, clustering, differential gene expression, and visualization. The data were normalized using the SCTransform method as implemented in Seurat. The function then selected the top 3000 most highly variable genes. PCA was performed on the subset of highly variable genes using the RunPCA function. The number of PCs to be used in the downstream analysis was chosen to be 17 based on examining the elbow in the elbow plot as implemented by Seurat.
We used the FindNeighbors and FindClusters functions to identify cell clusters in the data. Briefly, these functions define a graph connecting cells to each other by weighted edges and then identify clusters in the graph that place each cell into a single cluster using the Louvain algorithm. For the FindClusters function, we found that a resolution of 1.2 produced biologically meaningful clusters.
To identify differential genes between clusters of cells, Wilcoxon rank-sum test was used to compare gene expression of cells within the cluster of interest to all cells outside that cluster as implemented by the FindMarkers function. Specific parameters for the genes to be analyzed based on their log fold change between the two compared groups and percent of cells expressing the gene in at least one of the groups are available in the data analysis code. Uniform Manifold Approximation and Projection were used to reduce the PCs to two dimensions for data visualization using the RunUMAP function. For re-iterated analysis of the isolated stromal clusters, the same procedures were used with parameters changed as described in the data analysis code.
Cell-type identification. Following PCA and UMAP, we inspected clusters for expression of indicated markers using the differential gene expression results. Marker genes were plotted using an expression cut off to facilitate the visualization of both high-and low-expression genes on a single plot. Cutoffs are applied so that only cells with expression >cutoff received the color corresponding to that gene. These cutoffs are available in the data analysis code. In feature plots of multiple genes, for individual cells expressing multiple markers, each gene was over-plotted in the order described in the code. Feature plots of genes will be made available upon reasonable request.
Harmony analysis. To merge the previously published WT data set with the tumor data set, we used the Harmony algorithm 30 . First, the WT and tumor data set were analyzed in single SCTransform normalization and PCA steps. The Harmony algorithm then used the cells' PCA coordinates and data set identity to calculate new coordinates for each cell so as to minimize data set dependence when applying clustering to the cells. This algorithm produced a dimensional reduction that was used in place of PCA with the same steps applied to the data as described in the "Data normalization, clustering, differential gene expression, and visualization" and "Cell-type identification" subsections of the "Methods" section.