Autophagy mediates cancer intrinsic pan-immune evasion

To identify the underlying genes regulating glioma-intrinsic immune evasion across a spectrum of immune cell pressures, we performed genome-scale pooled CRISPR loss-of-function screens in a murine glioma model using the mTKO library [52]. CRISPR-mutagenized CT2A cells were propagated in the presence or absence of various immune cell lines [microglia; BV-2, non-phagocytic macrophages; Raw 264.7 and J774.1, phagocytic macrophages; J774.1 with anti-CD29-opsonized CT2A cells, cytotoxic T-lymphocytes (CTLs), or natural killer (NK) cells] (see Methods for details). Following a period of co-culture (i.e., selective pressure; Fig. S1), CT2A cells were subjected to deep sequencing of gRNA barcodes to identify genes that were enriched or depleted, i.e., genetic perturbations that conferred resistance or sensitivity to immune cell killing, respectively (Fig. 1a, b; Table S6).

Fig. 1

Mechanisms of CT2A-intrinsic immune evasion. a Workflow for mTKO genome-scale pooled CRISPR screens to identify immune-evasion genes. CRISPR-mutagenized CT2A cells were propagated in the present or absence of various immune cell lines (microglia; BV-2, macrophages; Raw 264.7 and J774.1, phagocytes; J774.1 treated with anti-CD29, cytotoxic T-lymphocytes, or natural killer cells) to apply selective pressure and CT2A cells were subjected to deep sequencing to identify sgRNA that were enriched (i.e., resister genes) or depleted (i.e., sensitizer genes) relative to untreated cells. b Rank-ordered z-score of sgRNA enriched/depleted in mutagenized CT2A cells after exposure to immune cells. Hits at FDR < 5% are highlighted in yellow (resistor genes) and blue (sensitizer genes). Point size is inversely scaled by FDR. c–e STRING network analysis of myeloid c and lymphoid d sensitizer genes, and resister genes (e). Clusters determined by Markov clustering. Nodes represents genes, and solid and broken edges represented intra- and inter-cluster connectivity, respectively. f Precision-recall (top) and ROC analysis (bottom) illustrating recovery of core CTL sensitizers and resisters identified by Lawson et al.[52] g Enrichment maps comparing CTL resisters (yellow) and resisters (blue) between CT2A and core sets. Nodes represent gene sets, and edges represent Jaccard similarities between gene sets. h GSEA for select pathways in in vivo ΔAtg12 CT2A tumors, compared to parental tumors, using snRNA-seq data. i Survival of C57BL/6 mice orthotopically engrafted with parental and ΔAtg12 CT2A cells. AUPRC, area under precision-recall curve; AUROC, area under receiver operating characteristic curve; CTL, cytotoxic T-lymphocytes; GSEA, gene set enrichment analysis; NK, natural killer cells; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins

CT2A-intrinsic non-phagocytic myeloid-evasion genes were defined as sensitizer or resister genes significant (5% FDR) across at least two myeloid cell lines (BV-2, Raw 264.7, or J774.1). This yielded 54 sensitizer and 8 resister hits (Fig. 1c, Table S6). In addition, we identified 69 sensitizers and 13 resisters involved in antibody-dependent cellular phagocytosis (ADCP; J774.1 + CT2A opsonized with anti-CD29). NFκB signaling (e.g., Traf2/3/6, Tnfaip3, Bcl2l1, Ikbkg and Ikbkb) and autophagy (e.g., Wipi2, Atg12) were shared sensitizing hits in phagocytic and non-phagocytic myeloid cells (Fig. 1c), likely due to residual non-phagocytic effects in the ADCP conditions (Fig. S1c). To clarify which genes were directly involved in ADCP evasion, we performed a CD29 sort screen in CT2A cells to identify regulators of the antibody target CD29 (encoded by Itgb1, Fig. S2a-b). Among the ADCP evasion genes not involved in regulating CD29 expression were Apmap and Cd47, both known inhibitors of phagocytosis (Fig. S2c), and several cytoskeletal (Mbtd1, Rab11a and Itgb5) and mediator complex (Cdk8) genes (Fig. S2d). Gene perturbations conferring resistance to myeloid-mediated killing were sparse, but included Tnfrsf1a, indicating the role of TNF in non-phagocytic myeloid-mediated cytotoxicity (Fig. 2e).

Fig. 2

Genetic dependencies in murine and human glioblastoma. a Workflow for mTKO genome-scale pooled CRISPR screens to identify fitness genes in CT2A and GL261 cells. b Distribution of gene-level differential logFC of sgRNAs in CT2A and GL261, stratified by essentiality. Gene fitness was scored using BAGEL. c Comparison of CT2A and GL261 gene-level fitness. Scatter plot shows CT2A and GL261 scaled BFs. Scaled BF was calculated as BF – 5 such that scaled BF > 0 represents essential genes. d Ranked differential fitness between GL261 and CT2A. Y-axis for differential fitness is signed log10(FDR) derived from difference between scaled BF scores. e Enrichment map illustrating CT2A and GL261-specific dependencies. Nodes represent gene sets, and edges represent Jaccard similarities between gene sets. f Scatter plot of scaled BF scores for human GBM cells and non-CNS cells. Scores were retrieved from Project Score Database (see methods). g Ranked differential fitness between human GBM and non-CNS cell lines. Genes were ranked by signed log10(FDR) derived from difference between scaled BF scores. h Venn diagram of human (GBM and non-CNS) and murine (CT2A and GL261) essential genes (scaled BF > 0). i Boxplot of scaled BFs from CT2A and GL261 screens grouped by human essentiality gene sets (as defined in f). j Dot plot of GBM-specific fitness genes that are common to human GBM and murine gliomas. BAGEL, Bayesian analysis of gene essentiality; BF, Bayes factor; CNS, central nervous system; ETC, electron transport chain; logFC, log fold-change

Like the myeloid screens, lymphoid screens revealed that perturbation of NFκB signaling and autophagy sensitized CT2A cells to CTL and NK killing (Fig. 1d, Fig. S3, Table S6). Furthermore, components of the chromatin remodeling pathway, including those seen in non-phagocytic myeloid evasion (i.e., Hdac2, Eed, and Dot1l) sensitized against CTL but not NK cells. CTL and NK cells had distinct resisters which reflected the unique mechanisms of anti-tumor immunity by each effector cell. Perturbation of GPI-anchor components (e.g., Pigu, Pigk, Dpm1, Dpm3, etc.) and Raet1e, which encodes the UL16-binding protein that serves as a NKG2D ligand [10], conferred resistance to NK-mediated cytotoxicity. Conversely, interferon (Ifngr1 and Ifngr2) and TNF (Tnfrsf1a and Tnfrsf1b) were required for CTL-mediated cytotoxicity (Fig. 1e).

We previously performed CTL coculture screens across six diverse syngeneic murine cancer cell models [colon; CT26 and MC38, kidney; Renca, breast; 4T1 and EMT6, melanoma; B16] to identify 182 core CTL evasion genes [52]. In CT2A cells, we found that these core CTL sensitizers, including genes involved in NFκB signaling and autophagy, were recovered with an AUROC of 0.81 and AUPRC of 0.30 (Fig. 1f). In contrast to other cancer cell lines, UFMylation, gene silencing, and GPI-anchored pathways were not involved in CT2A CTL evasion (Fig. 1g). Next, we found that core CTL resister genes were recovered in CT2A cells with an AUROC of 0.70 and AUPRC of 0.13, suggesting some contextual divergence in CT2A cells (Fig. 1f).

Given that autophagy was involved in pan-immune evasion including in CT2A cells, we sought to characterize the survival effect on mice engrafted intracranially with CT2A cells that have been genetically engineering with the autophagy pathway perturbed. Thus, clonal ΔAtg12 CT2A cell lines were engineered and engrafted orthotopically (Fig. S4). Single-cell transcriptome profiling of the engrafted ΔAtg12 CT2A tumors revealed significant downregulation of the autophagy pathway compared to parental controls (Fig. 1h). Decreased autophagy was associated with significant increases in apoptotic and TNFα/NFκB signaling, but not IFNγ signaling (Fig. 1h). Indeed, this also corresponded with a significant survival advantage (p < 0.01, Fig. 1i). Given that Atg12 was not an essential gene, we attributed this survival benefit to immune sensitization, rather than intrinsic impairment of tumor growth.

Taken together, the immune-glioma co-culture screens establish NFκB signaling, autophagy/endosome machinery, and chromatin remodeling as the predominant mechanisms of CT2A-intrinsic immune evasion.

Murine glioma cells partially recapitulate human genetic dependencies

In addition to identifying intrinsic immune evasion genes, we also defined fitness genes in CT2A cells and another murine glioma model GL261, then compared the results to similar screens performed in human GBM models. Pooled loss-of-function genetic screens were performed in CT2A and GL261 cells and essential fitness genes were identified using BAGEL (Fig. 2a, b, Table S7; BF > 5 threshold) and 1392 genes were deemed essential by this criterion (i.e., BF > 5) in both murine models, while 408 genes were GL261-specific and 250 were CT2A-specific (Fig. 2c). Notably, among the GL261-specific hits, Kras and Sox6 were top differential fitness genes, consistent with Kras being a known GL261 oncogene and Sox6 being a transcriptional regulator of the OPC-like GBM phenotype (Fig. 2d). Functional annotation of CT2A- and GL261-specific fitness genes further revealed that CT2A-specific fitness genes were enriched for processes involved in cell division and epigenetic and post-translational regulation of gene expression (e.g., RNA processing, spliceosome, cell division, histone modification) whereas GL261-specific genes were associated with metabolic processes (e.g., TCA cycle/ETC, nucleotide/flavin/cholesterol biosynthesis) (Fig. 2e).

We next evaluated the fitness landscape in human GBM cells. Comparison of gene essentiality profiles from 41 human GBM cell lines and 1031 human non-central nervous system (CNS) cell lines (Project Score database) identified 1625 common essential genes and 124 GBM-specific genes; notable GBM-specific human fitness genes included JUN, FERTM2, FGFR1, WWTR1 and ADAR (Fig. 2f, g). Of the 124 GBM-specific fitness genes identified in human cell lines, 44 (35%) and 54 (44%) genes were essential in CT2A and GL261 cells, respectively (Fig. 2h). However, by comparison, 51/123 (41%) and 68/123 (55%) of non-CNS-specific fitness genes were also essential in CT2A and GL261 cells. This suggests that CT2A and GL261 have unique genetic dependency profiles that resemble GBM in some ways, but not others. These findings were consistent across different essentiality thresholds and supported by precision-recall analysis (Fig. 2i, Fig. S5a). Among the human GBM-specific fitness genes that were recovered by GL261 were UFMylation-related genes, including Ufc1, Ube2g2 and Ufl1; these are known essential regulators of cell stress in human glioma stem cells [63]. Conversely, CT2A shared dependencies with human GBM cells related to epigenetic regulation (Dnmt1, Ttf1) and DNA damage response (Brat1, Rnf8) (Fig. 2j, Fig. S5b). Together our analyses provide insight into the genetic fitness landscape in CT2A and GL261 glioma models and highlight dependencies that are uniquely shared with human GBM.

Unbiased transcriptomic profiling of murine brain tumors

To further characterize CT2A and GL261 murine glioma models, each cell line was orthotopically engrafted into the right frontal hemisphere of immunocompetent C57BL/6 mice (Fig. 3a). PBS-injected mice were included as sham controls. Brain samples were collected at humane endpoint for sci-RNA-seq3 profiling [67]. We profiled 159,270 single cells with a median of 1786 UMI/cell and 1055 gene/cells (Fig. 3b, Fig. S6). Glioma and resident brain cells were identified using a combination of differential-expression analyses (Fig. S7, Table S8) and label-transfers from reference atlases (Fig. S8) [31, 51, 85, 123]. Anatomic information was also assigned by means of label-transfer of the spatially-resolved brain atlas (10 × Genomics, Adult Mouse Brain FFPE dataset; Fig. S9).

Fig. 3

Unbiased snRNA-seq profiling of glioma-engrafted mouse brains. a Workflow of snRNA-seq profiling of murine glioma models. CT2A and GL261 cells were expanded in vitro and orthotopically engrafted into the frontal right hemisphere of C57Bl/6 mice. At humane end point, brain tissue was sampled and nuclei profiled by sci-RNA-seq3. b UMAP of in vivo samples obtained from sham, GL261- and CT2A-engrafted mice. Neuronal populations are annotated using inferred anatomic (cerebellar, cerebral nuclei, cortical, hippocampal, hypothalamic and thalamic) and neurotransmitter (glutaminergic, GABAergic, glycinergic, dopaminergic, cholinergic) labels (see Methods). Numerical suffix corresponds to unique cluster identifier for each subpopulation

Owing to the lack of enrichment sorting prior to sci-RNA-seq3 analysis, our sci-RNA-seq3 profiles represent an unbiased snapshot of the intracranial milieu. This was reflected by the diverse representation of cells, including excitatory and inhibitory neurons, oligodendrocytes, astrocytes, lymphoid and myeloid cells, ependymal and meningeal cells, and CT2A or GL261 glioma cells (Fig. 3b, Fig. S7a, b, Fig. S8d). Inferred anatomic labels further reaffirm this diversity, with cell types arising from the cortex (CTX), cerebral nuclei (CNU), cerebellum (CB), hippocampus (HIP), hypothalamus (HY), thalamus (TH) and ventricles (VEN) (Fig. 3b, Fig. S9d).

Several markers distinguished GL261 and CT2A from non-malignant populations, including Hmga2, Piezo2, and B2m (Fig. S7a). In designing our experiments, we had intentionally used male mice, hypothesizing that sex-specific markers would discriminate the female-derived glioma lines from male host cells. Contrary to expectations, we found that Xist, a female-specific transcript, was only upregulated in CT2A and not GL261 cells (Fig. S7a). Bnc2 and Moxd1 were among the most sensitive and selective CT2A and GL261 markers, respectively, thereby representing gene markers than can be used to differentiate tumor cells from the surrounding microenvironment (Fig. S7a, b).

Glioma cells in vitro and in vivo have distinct transcriptomic signatures

We next assessed how the in vivo environment affects murine glioma biology (Fig. 4a). In vivo glioma cell engraftment led to increased transcriptomic dissimilarity (Fig. S10a) and decreased population purity (ROGUE score [56]; Fig. 4b) compared to in vitro conditions. Differential gene expression analysis revealed significant differences between in vitro and in vivo conditions in both glioma lines. Tcf4, a basic helix–loop–helix transcription factor that binds to specific DNA regulatory sequences (CANNTG) known as Ephrussi boxes (E-boxes) [114, 118], was the top up-regulated transcript in in vivo GL261 and CT2A cells, whereas Vim, a mesenchymal marker, was the top down-regulated transcript in vivo (Fig. 4c, Fig. S10b, Table S9). These transcriptomic changes were associated with a relative decrease (but not complete loss) of the mesenchymal-like phenotype and acquisition of oligodendrocyte progenitor-like (OPC) and neural progenitor-like (NPC) phenotype in both glioma models (Fig. 4d, e, Fig. S10c). Moreover, in vivo engraftment was associated with down-regulation of cell cycle, hypoxia and MYC-associated signaling (Fig. 4d, Fig. S10d–g).

Fig. 4

In vitro vs. in vivo comparison of syngeneic glioma models. a UMAPs of in vitro and in vivo GL261 and CT2A glioma cells. b In vitro vs. in vivo population purity (i.e., homogeneity), quantified by ROGUE [56] and compared by Wilcoxon test. c Differential gene expression between in vitro and in vivo GL261 and CT2A glioma cells. Log fold changes (logFCs) are compared between cell lines in sectored scatter plot. d, e Differential pathway activities between in vitro and in vivo GL261 and CT2A glioma cell. Differential activities are compared between cell lines in scatter plot (d) and representative GSEA plots are shown (e). f Volcano plot of differential expression between in vivo ΔTcf4 and parental CT2A cells. g Functional annotation of genes upregulated in in vivo ΔTcf4 cells, by hypergeometric gene set enrichment analysis. h, i Comparison of ΔTcf4 signature activity (h) and GSEA enrichment (i) in parental in vivo vs. in vitro GL261 and CT2A cells. j Proliferation assay in parental and ΔTcf4 CT2A clones. AC, astrocyte-like; GSEA, gene set enrichment analysis; MES, mesenchymal-like; NES, normalized enrichment score; NPC, neural progenitor-like; OPC, oligodendrocyte progenitor-like

We evaluated whether the differences acquired in vivo could be explained by Tcf4. We generated a clonal Tcf4 knockout CT2A cell line (ΔTcf4) using CRISPR-Cas9 and orthotopically engrafted these ΔTcf4 cells into the right frontal hemisphere of immunocompetent C57BL/6 mice. Tcf4 knockout led to upregulation of mesenchymal markers, including Col1a1, Col3a1, and Vim (Fig. 4f), and pathway analyses demonstrated significant enrichment for mesenchymal and MYC-related signaling (Fig. 4g). The ΔTcf4-associated signature effectively mimicked the gene expression profile in Tcf4-low in vitro glioma cells and was significantly depleted among the genes that were upregulated in the in vivo setting (Fig. 4h, i). Consistent with the high cell cycle signature observed in Tcf4-low in vitro glioma cells and previous reports [77] (Fig. 4d, Fig. S10g), ΔTcf4 CT2A cells proliferated significantly faster than parental cells (Fig. 4j).

Together these data demonstrated that GL261 and CT2A biology are influenced by environmental factors. Specifically, (i) in vivo glioma cells are more phenotypically heterogeneous than in vitro cells, (ii) in vivo engraftment of glioma cells impacts up-regulation of the NPC/OPC-like phenotype, and down-regulation of the mesenchymal-like phenotype, and (iii) in vitro cultures are more proliferative than in vivo glioma cells. Mechanistically, we found that the in vivo environment induces Tcf4 upregulation (or alleviates in vitro suppression), which in turn mediates these transcriptomic changes.

GL261 and CT2A tumors recapitulate canonical GBM transcriptomic phenotypes

The extent to which CT2A and GL261 tumors recapitulate human gliomas was next examined at the transcriptomic level. Using a transfer-learning-based approach (see Methods), we found that in vivo CT2A and GL261 tumor cells (Fig. 5a) had a higher degree of transcriptomic similarity to human Grade 4 primary GBM tumors than Grade 1 (low grade glioma; LGG) and recurrent Grade 4 recurrent GBMs (Fig. S11) [1, 120]. Given the resemblance to human GBM, we sought to determine whether murine gliomas recapitulate canonical GBM expression programs [82, 111]. We performed unsupervised gene program discovery using non-negative matrix factorization (NMF) in in vivo CT2A and GL261 cells (Fig. 5b, c, Table S10) and compared each program to established GBM and tumor-associated gene signatures (Fig. 5d). For each program we compared activity levels between CT2A and GL261 tumors (Fig. 5e, f), and evaluated the prognostic value using human survival data from The Cancer Genome Atlas (TCGA) program (Fig. S12).

Fig. 5

In vivo characterization of intrinsic GL261 and CT2A tumor biology. a UMAPs of in vivo GL261 and CT2A glioma cells. b Flowchart for NMF-based gene program discovery and annotation. c Heatmap of Jaccard similarity between component NMF programs used to derive consensus NMF programs in murine glioma models. d–f GL261- and CT2A-intrinsic gene programs were discovered using unsupervised NMF algorithm and characterized using hypergeometric gene set enrichment (d), gene program activity visualization on UMAPs (e), and differential gene program activity between CT2A and GL261 glioma cells (f). A, activity; H0, null hypothesis; NMF, non-negative matrix factorization

Altogether, we identified 8 gene programs, G1-G8, representing CT2A and GL261 intrinsic processes; three were GL261-biased (G5, G7, G8), four were CT2A-biased (G2, G3, G4, G6) and one was non-specific (G1; cell-cycle, 75 genes, e.g., Top2a) (Fig. 5c–e). G7 and G8 were associated with favorable survival in human glioma patients. G7 (87 genes, e.g., Sox6 and Ptprz1) represented a developmental-like program whereas G8 (99 genes, e.g., Met) had sparse functional annotations and was interpreted to be a GL261-specific signature. G5 [100 genes, e.g., Cd274 (encodes PD-L1), Irf1-2, Jak2, Tap1-2 and Stat1-3] was a GL261-biased inflammatory program associated with unfavorable survival outcomes in human gliomas (Fig. S12). Among the CT2A-biased programs, 3 of 4 were associated with mesenchymal processes. G4 (MES1; 88 genes, e.g., Fos/Fosb, Cd44, Nfkbiz and Vim) was associated with TNFα/NFκB signaling and epithelial-to-mesenchymal transition (EMT). G6 (MES2; 88 genes, e.g., Hk2 and Mxi1) was associated with glycolytic and hypoxic signaling. Finally, G2 (MES3; 79 genes, e.g., Prrx1, Pdgfra/Pdgfrb, Tfgb2 and Col1a1) was associated with angiogenesis, EMT and invasion. Among these only G4 was associated with unfavorable survival in glioma patients (Fig. S12). Finally, G3 (93 genes, e.g., Mast4) was a CT2A-enriched program with no known functional associations, and was interpreted as a CT2A-specific signature, akin to its GL261 counterpart G8. CT2A-specific G3 and GL261-specific G8 likely represent cell-line intrinsic programs with uncertain relevance to human GBM biology.

CT2A-enriched G4 and G6 programs directly mapped to the mesenchymal MES1 and MES2 GBM programs described by Neftel et al., whereas GL261-enriched G7 program mapped to Neftel’s OPC- and AC-like GBM programs [82]. Differential expression and pathway analysis corroborated these findings (Fig. S13).

In summary, CT2A and GL261 murine models recapitulate the canonical transcriptomic phenotypes of human GBM, and position GL261 and CT2A as developmental- and mesenchymal-like glioma models, respectively.

Human and murine gliomas have common transcriptional regulators

GBM is a notoriously heterogeneous and plastic tumor so understanding the transcriptomic regulators that govern different states may expose opportunities to bias tumors toward more therapeutically vulnerable states. Having established that GL261 and CT2A recapitulate canonical GBM phenotypes, we sought to define the transcription regulators responsible for these states. We bioinformatically identified GBM-associated transcriptional regulators (GTRs, Table S11) using a random forest machine-learning based strategy implemented across seven independent human GBM cohorts (N = 146 tumors, Fig. S14a). Three phenotypic axes were identified encompassing developmental (23 GTRs), mesenchymal (12 GTRs), and cycling-related processes (22 GTRs, Fig. 6a, Fig. S14b–d). In addition to Tcf4, which was identified as a developmental transcription factor, we selected 3 additional GTFs for experimental validation, including mesenchymal Wwtr1 and Prrx1, and developmental Nfia. For each candidate GTF, clonal GTR-perturbed CT2A lines were generated using CRISPR-Cas9 and engrafted into murine brains (Fig. S14e–h). At humane end point, mice were sacrificed, and brain tissue was sampled and subject to sci-RNA-seq3 profiling to evaluate the effect of each GTF perturbation on glioma biology (Fig. S14i–l, Table S12).

Fig. 6

Glioma transcriptional regulators. a Bipartite network illustrating relationship between GBM phenotypes (red nodes) and GTR activities (blue nodes). Edges represent random forest regression-derived feature importance scores, pooled across all human GBM datasets (Fig S14). b GSEA plots showing effect of Wwtr1 and Prrx1 perturbation in CT2A cells on developmental (G7-Dev) and mesenchymal (G4-MES1) phenotypes. (c) GTR essentiality scores (scaled BF) in CT2A, GL261, and human GBMs. Essential genes were defined as scaled BF > 0, where scaled BF = BF − 5. Bolded GTRs represent cycling associated GTRs that are essential across all glioma models. Differences (p values) between phenotypes were determined by ANOVA. BF, Bayes factor; Dev, developmental; GSEA, gene set enrichment analysis; GTR, glioma transcriptional regulators; MES, mesenchymal

As predicted bioinformatically, perturbation of mesenchymal GTRs Wwtr1 and Prrx1 resulted in developmental phenotypic shifts (Fig. 6b), whereas perturbation of developmental GTRs Nfia and Tcf4 resulted in mesenchymal shifts (Fig. 4g, Fig. S14i–l). GTR perturbations also resulted in the differential expression of other GTRs in patterns expected based on their inferred phenotypes (Fig. S14i–l). Finally, Sox6, although not interrogated here, was abundantly expressed in GL261—but not CT2A—thereby supporting its role as a developmental GTR (Fig. S13b, c).

Lastly, to validate the cycling-related GTRs, we analyzed pooled loss-of-function genetic screens in CT2A and GL261 cells (Fig. 2, Table S7), as well as human GBM (Project Score database) [21, 87]. We reasoned that GTRs implicated in the cycling-related phenotypic axis could be associated with glioma fitness in vitro. Of the 22 predicted cycling GTRs, 9, 8, and 11 were essential genes in CT2A, GL261 and human GBM cell lines, respectively, and five were essential across all models (i.e., Bub3, Cenpa, Bard1, Brca1, and Mis18bp1; Fig. 6c). By contrast, developmental and mesenchymal GTRs were overwhelmingly non-essential for cellular fitness, except for mesenchymal Eno1 in GL261 and CT2A, developmental Sox2/4/6 in GL261, and mesenchymal WWTR1 in human GBM lines (Fig. 4c). Among these off-target hits, Eno1 and Sox2/4/6 were inferred to have some cycling activity, thereby explaining their essentiality (Fig. S14b, Table S11).

These data represent a catalog of high-yield candidate GTRs and showcase the utility of CT2A in modeling GTR-associated phenotypic shifts. Furthermore, we provide experimental evidence supporting Wwtr1 and Prrx1 as mesenchymal GTRs, Nfia and Tcf4 as developmental GTRs, and Bub3, Cenpa, Bard1, Brca1, and Mis18bp1 as cycling GTRs.

Myeloid recruitment and cytokine signaling patterns distinguish the CT2A and GL261 tumor immune microenvironments

Human GBM is regarded as an immunosuppressive tumor, and as immunotherapies emerge to address this challenge, it is beneficial to understand the TIME in preclinical CT2A and GL261 models that partly served as the basis for these clinical trials [70, 84]. Thus, we characterized the immune microenvironment in CT2A and GL261 glioma models. We digitally sorted lymphoid and myeloid immune populations from sham, CT2A- and GL261-engrafted mice brains, and resolved 4 main types of immune cells with several distinct subpopulations observed at higher clustering resolutions: Macrophages (Mp; 7 subtypes), microglia (Mg; 2 subtypes), dendritic cells (DC; 2 subtypes), and T cells (TC; 3 subtypes) (Fig. 7a, b).

Fig. 7

Immune microenvironment in CT2A and GL261 tumors. a Gene program activity (top heatmap) and marker gene expression (bottom dot plot) in immune cells types. b UMAP of immune cells recovered from sham, GL261, and CT2A-engrafted brains. c Comparison of murine and human immune gene programs. Size of dots reflect degree of enrichment of murine gene sets in human gene sets, and color reflects correlation between murine and human gene program activities scored in murine immune population. d, e Inferred cytokine activities for each immune program. Immune response enrichment scores (IRES) were computed (IM_5 program shown as example) (d) and scores aggregated across each cell type were used to infer upstream cytokines activities (e). f Cytokine abundance in CSF from glioma patients. Significance determined by t test. Data from Fortuna et al. [