Functional profiling of murine glioma models highlights targetable immune evasion phenotypes

doi:10.21203/rs.3.rs-4946878/v1

Download PDF

Article

Functional profiling of murine glioma models highlights targetable immune evasion phenotypes

https://doi.org/10.21203/rs.3.rs-4946878/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Cancer intrinsic immune evasion mechanisms and pleiotropy represent a barrier to effective translation of cancer immunotherapy. This is acutely apparent for certain highly fatal cancers such as high-grade gliomas and glioblastomas. In this study, we use functional genetic screens, single-cell transcriptomics and machine-learning approaches to deeply characterize murine syngeneic glioma models in vitro and in vivo, and compare-and-contrast their value as preclinical models for human glioblastoma (GBM). Systematic genome-wide co-culture killing screens with cytotoxic T cells, natural killer cells or macrophages established NFkB signaling, autophagy/endosome machinery, and chromatin remodeling as pan-immune cancer intrinsic evasion mechanisms. Additional fitness screens identified dependencies in murine gliomas that partially recapitulated those seen in human GBM (e.g., UFMylation). Different models associated with contrasting immune infiltrates including macrophages and microglia, and both models recapitulate hallmark immune gene programs seen in human GBM, including hypoxia, interferon and TNF signaling. Moreover, in vivo orthotopic tumor engraftment is associated with phenotypic shifts and changes in proliferative capacity, with models recapitulating the intratumoral heterogeneity observed in human GBM, exhibiting propensities for developmental- and mesenchymal-like phenotypes. Notably, we observed common transcription factors and cofactors shared with human GBM, including developmental (Nfia, Tcf4), mesenchymal (Prrx1 and Wwtr1), as well as cycling-associated genes (Bub3, Cenpa, Bard1, Brca1, and Mis18bp1). Perturbation of these genes led to reciprocal phenotypic shifts suggesting intrinsic feedback mechanisms that balance in vivo cellular states. Finally, we used a machine-learning approach to identify evasion genes that revealed two gene programs, one of which represents a clinically relevant phenotype and delineates a subpopulation of stem-like glioma cells that predict response to immune checkpoint inhibition in human patients. This study offers relevant insights and serves to bridge the knowledge gap between murine glioma models and human GBM.

glioma

glioblastoma

murine

human

CT2A

GL261

scRNA-seq

genome-wide CRIPSR screen

Pre-clinical models of cancer that faithfully recapitulate the complexity of cancer cells and their interactions with the tumor microenvironment are essential for investigating therapeutic strategies. Commonly used in vitro systems for cancer research, such as conventional 2D cell culture or 3D organoids, are unsatisfactory for modeling the complexity of tumor microenvironments, particularly for complex tissues such as the brain¹. Gliomas are the most common primary brain tumor in adults and vary in severity and aggressiveness with subtypes that include astrocytomas, oligodendrogliomas and glioblastomas^2,3. Glioblastomas (GBM) are malignant grade 4 gliomas with no evidence of a lower-grade precursor and are predominantly made up of abnormal astrocytic cells that are diffusely infiltrative and invasive^2,3. Drug development efforts for GBM are largely impeded by substantial financial risks associated with historical translational failures⁴, and the lack of preclinical models that recapitulate the complexity and heterogeneity of the disease in patients⁵. Intentional functional discovery and drug development efforts are vital for GBM given its unique therapy challenges including heterogeneity, an immunosuppressive tumor micro-environment, and the presence of the blood-brain or blood-tumor barrier⁵.

Murine syngeneic GL261 and CT2A glioma models are commonly used to study glioma biology in immunocompetent C57BL/6 mice⁶. Both models recapitulate various hallmarks of GBM, and their unique mutational^7-13 and transcriptomic⁹ profiles result in distinct responses to immune checkpoint inhibitors (ICI)^9,11,14,15, radiation^16,17, and chemotherapy^16-18 (Table 1). Histologically, both tumors have characteristic pseudopalisading necrotic cores, along with extensive angiogenesis. Unlike GL261 tumors, CT2A tumors are known to resemble a mesenchymal phenotype that correlates with its immune-suppressive tumor-immune microenvironment (TIME) and resistance to immunotherapies. Given that preclinical findings in GL261 and CT2A models often fail to correlate with clinical findings, it behoves us to better understand how these syngeneic glioma cell lines deviate from each other and from human GBM.

In the current study, we used experimental and computational approaches and evaluated murine glioma models using genome-wide CRISPR-Cas9 screens to identify cancer intrinsic immune evasion genes to cytotoxic T cells, natural killer cells and macrophages, cancer gene dependencies, and single cell transcriptomics to evaluate the TIME. We also compared murine glioma models using public single cell RNA-sequencing (scRNA-seq) and CRISPR data to determine how these models recapitulate what is known about human GBM. We report on numerous relevant findings, including in vitro-to-in vivo changes that are acquired upon tumor engraftment, intratumoral heterogeneity, phenotypic regulators, the TIME, and intrinsic mechanisms of immune evasion with relevance to immunotherapy.

Table 1. Literature review of murine glioma models evaluated in this study¹⁹
Feature	CT2A	GL261
Origin	C57BL/6 mouse	C57BL/6 mouse
Tumorigenesis	Methylcholanthrene-induced²⁰	Methylcholanthrene-induced²¹
Histology	High-grade astrocytoma¹⁰ Microvascular proliferation, angiogenic^10,22 Pseudopalisading necrosis¹⁰ Spindled cells, fascicular tissue⁹	Ependymoblastoma²³ Poorly differentiated⁹ Pseudopalisading necrosis²⁴ Angiogenic²⁴
Transcriptional Profile	Mesenchymal/angiogenic/WNT⁹	Interferon signaling⁹
Genetic Alterations	¯ PTEN/TSC2^7,8 N-ras* mutation⁹ Cdkn2a/b heterozygous deletion⁹ mutation load⁹ Idh1/2, Trp53** wild-type^9,10	K-Ras/Trp53 mutations¹¹ C-MYC expression¹¹ mutation load^9,12 Idh1/2 wild-type^9,13
Growth	CT2A tumors are more aggressive than GL261 tumors⁸
Immune Response	Low immunogenicity^9,14 Interferon/antigen presentation deficits⁹ myeloid infiltration^8,9	Moderate immunogenicity^9,11,14,15 MHC I expression¹² microglial activation⁸
Radiation Response	+++Radiosensitive^16,17	++Radiosensitive^16,17
TMZ Response	TMZ resistant^16,17	TMZ sensitive^16-18
Based on protein expression; not mutated⁹. *Based on nuclear immunohistochemical staining of p53 protein CT2A tumor mass Abbreviations: TMZ; temozolomide.

Autophagy mediates 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 ²⁵. 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 S1).

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 S1). Additionally, we identified 69 sensitizers and 13 resisters involved in antibody-dependent cellular phagocytosis (ADCP; J774.1 + CT2A opsonized with anti-CD29). NFkB 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).

Like the myeloid screens, lymphoid screens revealed that perturbation of NFkB signaling and autophagy sensitized CT2A cells to CTL and NK killing (Fig. 1D, Fig S3, Table S1). 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²⁶, 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²⁵. In CT2A cells, we found that these core CTL sensitizers, including genes involved in NFkB 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 TNFa/NFkB signaling, but not IFNg 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 coculture screens establish NFkB 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 in CT2A and GL261 cells and essential fitness genes identified using BAGEL (Fig. 2A-B, Table S2; BF > 5 threshold). 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 (G7 program; 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²⁷. 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²⁸. We profiled 159,270 single cells with a median of 1,786 UMI/cell and 1,055 gene/cells (Fig. 3B, Fig S6). Glioma and resident brain cells were identified using a combination of differential-expression analyses (Fig S7, Table S3) and label-transfers from reference atlases (Fig S8)^29–32. Anatomical information was also assigned by means of label-transfer of the spatially-resolved brain atlas (10x Genomics, Adult Mouse Brain FFPE dataset; Fig S9).

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). Inferred anatomical 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 S8D).

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 (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³³; 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)^34,35, was the top upregulated transcript in in vivo GL261 and CT2A cells, whereas Vim, a mesenchymal marker, was the top downregulated transcript in vivo (Fig. 4C, Fig S10B, Table S4). These transcriptomic changes were associated with a relative loss of 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 downregulation of cell cycle, hypoxia and MYC-associated signaling (Fig. 4D, Fig S10D-G).

We also 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 (see Methods). 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³⁶ (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 IV primary GBM tumors than Grade II (low grade glioma; LGG) and recurrent Grade IV recurrent GBMs (Fig S11)^37,38. Given the resemblance to human GBM, we sought to determine whether murine gliomas recapitulate canonical GBM expression programs^39,40. We performed unsupervised gene program discovery using non-negative matrix factorization (NMF) in in vivo CT2A and GL261 cells (Fig. 5B-C, Fig. 12, Table S5) 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. 13).

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 determined 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 S13). 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α/NFkB 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 S13). 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³⁹. Differential expression and pathway analysis corroborated these findings (Fig S12).

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 towards 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, see Methods for description and Table S6) 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 S7).

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 S12B-C).

Lastly, to validate the cycling-related GTRs, we analyzed pooled loss-of-function genetic screens in CT2A and GL261 cells (discussed below; Fig. 7, Table S4), as well as human GBM (Project Score database)^41,42. 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 S14D, Table S6).

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 beneficiary to understand the TIME in preclinical CT2A and GL261 models that partly served as the basis for these clinical trials^43,44. 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. 8A, B).

To functionally annotate the immune populations, we performed NMF-based gene program discovery and resolved 9 murine immune programs (denoted IM1-9; Fig S15A-C; Table S8). The macrophage population was the most abundant and heterogeneous, consisting of distinct subpopulations involved in TNFα signaling (Mp-4 cells; IM3 program), IFN signaling (Mp-9 cells; IM5 program) and antigen presentation (Mp-8, Mp-9 and Mp-14 cells; IM9 program) (Fig. 7A, top heatmap). There were also macrophages under hypoxic stress (Mp-4, Mp-5 and Mp-14 cells; IM8 program). Among dendritic cells (DC), DC-12 was associated with higher MHC class II-mediated antigen presentation (Cd74, H2-Ab1; IM9 program) whereas DC-15 was more pro-inflammatory, consisting of higher TNFα (IM3) and IFN (IM5) signaling. Unlike the other myeloid lineages, microglial cells (P2ry12) were relatively homogenous. T cells varied in processes related to cell cycle (IM6) and differentiation (IM1) and included a T-regulatory subpopulation (TC-11; Ctla4).

To determine the relevance of the TIME in GL261 and CT2A gliomas to human GBM, we evaluated the TIME across three human GBM cohorts^37–39. Notably, biased sampling of immune cells from human samples (enrichment sorting) precluded direct compositional comparison of the murine and human TIMEs. Instead, we performed unsupervised gene program discovery and annotation using digitally sorted immune cells from human GBM samples and resolved 10 human immune programs (denoted IH1-10; Fig S15D-E, Table S8). We found that all murine immune programs were recapitulated in the human TIME (Fig. 7C, Fig S15F), suggesting that the TIME in glioma-engrafted mice recapitulates hallmark features of the human immune response in GBM.

To elucidate the cytokines driving these diverse immune gene programs, we leveraged a scRNA-seq-derived cytokine response dictionary to calculate immune response enrichment scores (IRES) for each gene program (Fig. 7D-E)⁴⁵. Select cytokines were non-specifically linked to multiple immune programs, like IL1a/b in the DC-specific response (IM4), T-cell differentiation (IM1), hypoxia (IM8), and antigen-presentation (IM9), suggesting pleiotropic and cell-type specific cytokine responses. Conversely, other programs, like the microglial-specific response, was linked to multiple cytokines, including SCF, Noggin, and IL30. Notably, inflammatory TNFa (IM3) and Interferon- (IM5) programs mapped to their cognate cytokines. Of the inferred cytokines, many were significantly enriched in CSF sampled from glioma patients⁴⁶, including IFNg (p = 0.015), TNFa (p = 0.052), EGF (p = 0.082), IL1b (p = 0.022), but not IL1a (p = 0.84), thereby linking the extracellular presence of select cytokines to the downstream responses observed in the TIME.

Earlier we reported distinct inflammatory responses in CT2A and GL261 tumor cells; CT2A tumors were associated with TNFα/NFkB signaling (G4 program) whereas GL261 tumors had higher levels of IFN signaling (G7 program) (Fig. 5D-F). This led us to hypothesize that CT2A and GL261 have distinct intrinsic sensitivities to different cytokines. To test this, we performed cell viability assays in CT2A and GL261 cells treated with IFNγ and TNFα (Fig. 7G). GL261 cells (IC50_IFNγ < 5 ng/mL) were > 20-fold more sensitive to IFNγ than CT2A cells (IC50_IFNγ = 120 ng/mL), whereas the inverse was observed in TNFα treated GL261 (IC50_TNFα = 85 ng/mL) and CT2A (IC50_TNFα < 5 ng/mL) cells.

Finally, we evaluated the composition of the TIME in CT2A and GL261-engrafted mice. The relative abundance of immune cells in CT2A (5.3%) and GL261 (2.1%) was significantly higher than sham (1.1%) (Fig. 7H). Regardless of glioma model, tumor engraftment was associated with significant infiltration of pro-inflammatory (IFN-signaling) macrophages (Mp-9), cycling T-cells [TC-6 and TC-11 (Ctla-positive T regulatory cells)], and antigen-presenting (DC-12) and pro-inflammatory (DC-15) dendritic cells (Fig. 7I). Comparison of CT2A and GL261 immune infiltrates showed that Mp-4 and Mp-8 macrophage populations were unique to CT2A tumors, and otherwise absent in sham control and GL261-engrafted mice (Fig. 7I, J). In contrast, Mg-3 and to a lesser extent Mg-2 microglia were over-represented in the GL261 TIME compared to CT2A. T-cells and dendritic cells were equally represented in both tumor models.

To summarize, here we characterized the TIME in GL261 and CT2A tumors, including the active immune programs and their associated cytokines, and TIME composition.

Intrinsic immune evasion of stem-like glioma cells

Given the prevalent heterogeneity in GBM, we hypothesized there to be intra-tumoral variation in immune evasion activity. We performed NMF to identify evasion genes with common patterns of expression that were reproducible across six independent CT2A sci-RNA-seq3 experiments (Fig. 8A). Of the 265 evasion genes identified (Fig. 7), 114 genes clustered into two gene programs, denoted as evasion phenotype E1 (57 genes) and E2 (57 genes) (Fig. 8B, Table S9). E1 and E2 activities were negatively correlated (r = -0.33) and delineated two mutually exclusive glioma subpopulations (Fig. 8C). Interestingly, each evasion phenotype was comprised of a mixture of resister and sensitizer hits (Fig. 8D), with little overlap with Lawson’s core evasion genes (Fig. 8E)²⁵, possibly reflecting glioma specificity. E1 was associated with sensitizer genes implicated in Toll-like receptor signaling, and resister genes implicated in TGFb/SMAD signaling and Runx1/Chromatin regulation resisters (Fig. 8F). Alternatively, E2 was associated with sensitizer genes involved in chromatin remodeling and resister genes involved in apoptosis and GPI anchor activity (Fig. 8F). We also explored which glioma state each evasion phenotype was associated with. Using a curated list of tumor/glioma-associated gene programs, intratumoral correlations were calculated and pooled across 69 human tumors from three independent human cohorts^37,47. This revealed that the E1 and E2 phenotype activities were correlated with neurodevelopmental-like and stem-like GBM states, respectively (Fig. 8G, H).

We next evaluated clinical correlates for each evasion phenotype. Meta-analysis of survival data from several human glioma cohorts (973 LGG patients across 2 cohorts and 638 GBM patients across 5 cohorts) revealed that high E2 activity was associated with worse survival outcomes in LGG [HR (95% CI) = 4954 (577, 42498); p < 0.0001, I² = 0%] and GBM [HR (95% CI) = 10.9 (1.7, 71.5); p = 0.013, I² = 64%] (Fig. 8I). In other words, this corresponded to a 50.4-month and 5.4-month survival difference between LGG and GBM patients, when stratified by E2 activity levels (Fig. 8J). E2 activity was also positively associated with WHO grade gliomas in the TCGA (Fig. 8K, p = 1.65×10^− 29) and CCGA (Fig. 8L, p = 1.23×10^− 14) cohorts, and with GBM recurrence across three additional cohorts profiled by scRNA-seq (Fig. 8M, p = 0.046). Unlike E2, the E1 phenotype had no significant clinical correlates (data not shown). Finally, we hypothesized that tumor intrinsic evasion phenotypes can be leveraged to predict response to ICIs. Using RNA-seq profiles obtained from patients at baseline (i.e., prior to immunotherapy)⁴⁸, we found that E2 was significantly upregulated in ICI (PD-1 inhibitor) responders (Fig. 8N, p = 0.024). We benchmarked the evasion phenotype signatures against other immunotherapy response indices that have been proposed, including CD274 expression, interferon signature (IFNG)⁴⁹, Estimation of STromal and Immune cells in Malignant Tumors (ESTIMATE)⁵⁰, Tumor Immune Dysfunction and Exclusion (TIDE)⁴⁹ and Immune Escape-Related Gene Prognosis Index (IEGPI)⁵¹, and found that the E2 signature outperformed all, with an AUROC of 0.74 (Fig. 8O, P).

Overall, these findings show that intrinsic immune evasion genes have distinct patterns of expression, E1 and E2. E2 represents the more clinically relevant phenotype and delineates a subpopulation of stem-like GBM cells that are associated with worse prognosis, higher WHO grade and tumor recurrence. Importantly, high E2 activity is predictive of ICI response, and outperforms all other predictive indices.

In the current study we evaluated data from single-cell RNA sequencing and genome-wide pooled CRISPR screening approaches to compare-and-contrast functional dependencies present in two syngeneic murine models of glioma. Not surprisingly, comparison of in vitro and in vivo snRNA-seq profiles demonstrated a profound influence of the in vivo microenvironment on cell state, resulting in increased tumor heterogeneity, downregulation of the mesenchymal and stress response, and lower proliferative capacity secondary to Tcf4 upregulation. Unsupervised gene program discovery in in vivo tumor cells further revealed that CT2A and GL261 cells are mesenchymal- and developmental-like tumors, respectively. The gene programs in human GBM and murine glioma were regulated by common GTFs, and we experimentally validated several developmental (Nfia, Tcf4), mesenchymal (Prrx1 and Wwtr1) and cycling-associated (Bub3, Cenpa, Bard1, Brca1, and Mis18bp1) GTFs in CT2A glioma cells. Genome-wide CRISPR-Cas9 screens revealed distinct genetic dependencies in CT2A (epigenetic and post-translational regulation) and GL261 (metabolic) cells and demonstrated that murine gliomas recapitulate various GBM-specific genetic dependencies (e.g., UFMylation in GL261). Moreover, the murine TIME was found to be macrophage-dominant in CT2A tumors, and microglial-dominant in GL261 tumors. Immune-glioma co-culture screens established NFkB signaling, autophagy/endosome machinery, and chromatin remodeling as the predominant mechanisms of CT2A-intrinsic immune evasion. Lastly, we discovered that cancer intrinsic immune evasion genes have heterogeneous patterns of expression. Specifically, the E2 evasion phenotype was associated with a stem-like GBM subpopulation that was correlated with prognosis, WHO grade and recurrence, and predicted response to ICI.

Glioma cells engrafted orthotopically are subject to selective in vivo pressures that include nutrient limitations, hypoxia (21% oxygen in vitro vs. 0.3–7.4% in vivo⁵²), and immune cell interactions, while in vitro cells are exposed to artificial culture substrates and media. The transcriptomic changes observed upon in vivo CT2A and GL261 tumor engraftment revealed significant downregulation of proliferation and mesenchymal programs, thereby confirming reports by others in GL261 and 4T1 murine mammary carcinoma cells^53,54. Notably, 13% of genes in the mesenchymal signature (Neftel MES2³⁹) overlapped with the in vitro stress signature⁵⁵ (including DDIT3, HSPA5 and HSPA9), and these were downregulated upon in vivo engraftment. This raises the distinct possibility that the mesenchymal state in vitro is an artifact of in vitro stress, as suggested by others⁵⁶. However, the TCF4 gene product has also been proposed as a metabolic sensor that is upregulated in response to metabolic demands⁵⁷. High glucose concentrations in DMEM have been associated with epithelial-to-mesenchymal (EMT; e.g., Vim, Cd44) upregulation in MDA-MB-231 breast cancer cells⁵⁸. Similarly, hyperglycemia has induced EMT and HIF1a/hypoxia signaling in other cell lines^59–61. While this supports a model in which limited in vivo glucose availability induces mesenchymal downregulation in a Tcf4-dependent manner, further experimental investigations were out of the scope of the current study. Instead, we found that perturbation of Tcf4, which was upregulated in vivo, was able to revert many of the in vivo-induced changes back to an in vitro-like state, specifically leading to increased cell cycle and mesenchymal-like program activities. TCF4:TCF12 dimer activity in a primary GBM line was previously shown to inhibit cellular proliferation and suppress POSTN expression, while TCF4 knockdown led to increased POSTN expression, similar to that seen in our experiments³⁶. Periostin, encoded by Postn, promotes EMT, invasion and integrin-mediated adhesion in glioma cells⁶², and its negative regulation by Tcf4 may in part explain the mesenchymal shift observed in vitro and Tcf4 perturbed glioma cells.

Intratumoral heterogeneity is a hallmark feature of human GBMs, and the subtypes of GBM have been exhaustively characterized and include the Verhaak subtypes (classical, mesenchymal, proneural and neural)⁴⁰, Neftel subtypes (MES; mesenchymal, OPC; oligodendrocyte progenitor-like cells, NPC; neural progenitor-like cells, and AC; astrocyte-like cells)³⁹, and Richards subtypes (developmental and injury-response)⁶³. Our unsupervised gene program discovery revealed that CT2A and GL261 gliomas are predominantly mesenchymal- and developmental-like tumors, respectfully, as suggested by others⁹. Importantly, perturbation of developmental (Nfia, Tcf4) and mesenchymal (Wwtr1, Prrx1) GTRs led to reciprocal phenotypic shifts. Many of the GTRs predicted here have been implicated in GBM by others in shaping GBM identity, including EZH2^64,65, FOXM1^66,67, NFIA/B^68–72, OLIG1/2^73–76, SOX2^77–79, SOX4^80–82, SOX6⁸³, SOX8^84,85, TCF4^36,86, TCF12^87,88, and ZEB1^68,89.

CT2A and GL261 engraftment resulted in significant immune recruitment. Among the myeloid compartment, macrophages were overrepresented in CT2A tumors whereas microglia were overrepresented in GL261 tumors, thereby confirming earlier reports^8,14. Although tumor-associated dendritic cell infiltrates were similar in both models, they were still higher than that observed in human GBM⁸. CT2A cells secrete significantly higher levels of myeloid chemokines (including CCL-2, CCL-5 and CCL-22) than GL261s⁹, possibly contributing to the differential pattern of myeloid recruitment. Owing to the limited recovery of rare cell types by sci-RNA-seq3⁹⁰, we were unable to detect NK, NKT, and B-cells, which are known to be present in GL261 and CT2A TIMEs^8,32,91,92. Using mass cytometry, Khalsa et al. reported no significant differences in the abundance of T cell, NK, NKT or B cell fractions between GL261 and CT2A tumors⁸. We also observed no differences in the T cell populations infiltrating the two models, however Khan et al. showed that CT2A tumors were enriched for exhausted CD8 + and regulator CD4 + T cells, whereas GL261 tumors were enriched for progenitor exhausted CD8 + T cells⁹³. Although the TIME composition in GL261 more closely resembles human GBMs than CT2A⁸, the murine TIME shares key gene programs that are characteristic of human GBMs, including TNF, IFN and hypoxia.

Gene perturbations conferring resistance to immune cell killing revealed common and distinct anti-tumor mechanisms in different immune populations. CTLs mediated CT2A killing via TNF (Tnfrsf1a and Tnfrsf1b) and IFN (Ifngr1 and Ifngr2), whereas non-phagocytic myeloid cells predominantly relied on TNF (Tnfrsf1a). Notably, NK-mediated cytotoxicity was dependent on GPI anchor-mediated signaling (e.g., Piga, Pigh, Pigm, etc.) and UL16-binding protein (Raet1e). GPI anchors are essential for the surface expression of UL16-binding proteins (i.e., ligands for the activating NK cell receptor NKG2D⁹⁴) and subsequent secretion of cytolytic granzyme and perforin from NK cells⁹⁵. Although other NK screens have implicated TNF, IFN or antigen-presentation in NK-mediated cytotoxicity, we were unable to reproduce these findings. However, there was a report that impaired TNF-mediated signaling (e.g., TNFRSF10A, TNFRSF10B) conferred resistance to NK killing in GPI-deficient HAP1 cells thereby suggesting that GPI anchor resistance may mask the effects of other antitumor mediators^96–98. Phagocytic myeloid evasion was dependent on known phagocytosis inhibitors Cd47 and Apmap⁹⁹.

Our sci-RNA-seq3 analyses supported the abundance of TNF and IFN in the TIME, albeit with noteworthy differences between CT2A and GL261. CT2A tumors (G4-MES1 program) and macrophages (Mp-4 population) had higher levels of TNF signaling compared to GL261. Conversely, IFN signaling was similar among macrophages (Mp-9 population) infiltrating CT2A and GL261 tumors, but only observed in GL261 tumor cells (G5-Inflammatory program). The attenuated IFN response in CT2A tumor cells, confirmed by our IFN dose-response experiment, was due to intrinsic resistance of CT2A to IFN secondary to known single-copy deletions of chromosomal regions encompassing type I IFN genes, and Stat2, Stat6, and Ifng⁹. However, despite this reduced sensitivity, impaired IFN signaling stilled conferred resistance to CTL killing in CT2A cells. Defects in IFN response have been associated with resistance to checkpoint immunotherapy^100,101, thus positioning the CT2A line as a relevant preclinical model for optimizing immunotherapies.

The autophagy-NFkB axis has previously been implicated in the immune evasion of various cancers to CTLs (REFs). Our work expands this to encompass additional immune cells, including myeloid and NK cells, thereby establishing the autophagy-NFkB axis as pan-immune evasion pathway. NFkB signaling deficits are known to sensitize tumor cells to immune cell killing by promoting Caspase 8-mediated apoptosis downstream of TNF signaling (e.g., Rnf31)^25,102 and upregulating expression of MHC1 (e.g., Tnip)¹⁰³. Similarly, autophagy mediates tumor immune evasion through multiple mechanisms, including granzyme clearance (e.g., Vps4bi, Atg5)¹⁰² and TNF resistance (e.g., Atg12)²⁵. Consistent with this, we showed that autophagy deficient CT2A tumors (ΔAtg12) were associated with elevated TNFa and apoptotic signaling, and favorable survival outcomes.

Contrary to other tumors²⁵, we found that CTL-mediated cytotoxicity in CT2A cells was independent of MHC1 antigen presentation (e.g., B2m, H2-K1, Psmb9, Tap2). Iorgulescu et al. reported that CT2A cells harbor multiple mutations in antigen presentation machinery (e.g., Tap1 p.Y488C, Psmb8 p.A275P), however that the defect in MHC1 expression could be overcome with IFNg treatment⁹. Consistent with this, we found that perturbation of Ptpn2, a protein tyrosine phosphatase that negatively regulates IFNg-mediated effects on antigen presentation¹⁰⁴, sensitized CT2A cells to CTL killing. Human GBMs are among the most common tumors associated with loss of MHC1 expression (e.g., HLA, b2M), and MHCI loss is independently associated with unfavorable clinical outcomes¹⁰⁰. In some 1p/19q intact IDH-mutant gliomas, recurrence is associated with loss of heterozygosity in HLA genes¹⁰⁵. These findings position CT2A as a relevant preclinical model of MHC1-deficiency in which to study immunotherapies.

We showed that immune evasion genes were heterogeneously expressed, with a subset (E2) enriched among stem-like GBM cells. This goes to suggest that mechanisms of immune evasion can vary even within individual tumors. Moreover, it has been hypothesized that immunosuppressive factors are enriched among stem-like tumor cells to establish an immune-privileged microenvironment in which clonal diversification and expansion can occur¹⁰⁶. Perturbation of Sox2, a known stemness marker, sensitized CT2A cells to J774- and CTL-mediated cytotoxicity (Table S1). Others have also shown that glioma stem cells escape immune surveillance by downregulating antigen presentation¹⁰⁷. Notable evasion mechanisms captured in the E2 signature included the phagocytic inhibition (Apmap), TNF (Tnfrsf1a) and interferon (Ifngr2) signalling, GPI anchorage (Pigt, Gpaa1, Dpm3), and NFkB signaling (Nfkb2), thereby reflecting a multimodal evasion phenotype. The E2 signature outperformed other immune indices in predicting response to immune checkpoint inhibitors in GBM patients with a competitive AUROC of 0.74. This highlights the translational capacity of CT2A gliomas in identifying mechanisms of immune evasion, with direct implications for identifying patients that will response to ICI therapy.

This study was not without limitations: i) We were unable to recover rare cell populations using sci-RNA-seq3. Furthermore, sci-RNA-seq3 experiments involving Atg12 and GTR-perturbed CT2A cells did not recover any immune cells, thereby precluding the evaluation of the TIME in these tumors. ii) Immune coculture screens were only performed in CT2A cells. iii) While our GTR inference was able to reliably identify GBM phenotype-associated GTRs, the direction of regulation at times could not always be reliably inferred by our approach. For example, Tcf12 and Zeb1 were predicted to promote a developmental phenotype, where in fact the associations are negative^88,89.

Our study offers relevant insights into the biology of GL261 and CT2A and serves to bridge the gap between murine models and human GBM.

Resources and Data sources

Software. Endnote X9 (Thomson Reuters) was used to manage references, R (version 4.2.2) was used for data analysis, Excel for Office 365 (Microsoft) was used for data storage, and CorelDRAW X8 (Corel) was used for Figure preparation.

Computational Resources. Analyses were run on a desktop computer with an Intel Core i9-10900L CPU (3.70 GHz, 10 cores, 20 threads) with 120 GB RAM running Windows 10 Pro (v21H2).

Public data sources. scRNA-seq data (Table S10) from Yu et al. were obtained from Gene Expression Omnibus (GSE117891)³⁸; Neftel et al. from GSE131928³⁹; Abdelfattah et al. from GSE182109³⁷; Qazi et al. from GSE196583¹⁰⁸. Richards et al. data was obtained from the Broad Institute Single-Cell Portal (SCP503)⁶³. Spatial transcriptomics data of coronal section from 8-week male C57Bl/6 mice was from 10x Genomics (Adult Mouse Brain FFPE dataset). Bulk RNA-seq data from the Cancer Genome Atlas (TCGA) RNA-seq V2 data from the TCGA PanCancer Atlas were retrieved from the National Cancer Institute (NCI) Genomic Data Commons using TCGAbiolinks R package v.2.16.0; Chinese Glioma Genome Atlas (CCGA) from http://www.cgga.org.cn/download.jsp¹⁰⁹; Glioma Longitudinal Analysis (GLASS) consortium from Synapse (http://synapse.org/glass)¹¹⁰; Zhao et al. PDL1 responder vs. non-responder data from SRA (PRJNA482620)⁴⁸; and Cloughesy et al. from GEO (GSE121810)¹¹¹.

Animal studies

Animal Studies. Animal use and studies were performed in accordance with guidelines outlined by Animal Use Protocols within the Division of Comparative Medicine at the University of Toronto and the McMaster University Central Animal Facility. Intracranial injections were conducted in 6-8-week-old male C57BL/6 mice as previously described using 10⁶ cells per mouse¹¹². Briefly, a small burr hole was drilled 2 mm behind the coronal suture and 3 mm to the right of the sagittal suture. Cells suspended in 10 µL PBS were injected intracranially using a Hamilton syringe (Hamilton, #7635-01) into the right frontal lobes. Animals were sacrificed at humane endpoint, defined loss of 15-20% of bodyweight from date of tumor engraftment, and a combination of hunched posture, decreased response when picked up, lack of grooming and loss of skin elasticity. Brains were removed, and right frontal pieces were snap frozen for sci-RNA-seq3 processing.

Cell lines

CT2A and GL261 glioma cells. CT2A cells were purchased from Millipore/Sigma-Aldrich (#SCC194) and GL261 cells were purchased from DSMZ (#ACC802). These lines were maintained in Dulbecco’s modified eagle medium (DMEM; Wisent, #319-005-en) supplemented with 10% fetal bovine serum (FBS; Gibco, #12483-020) and 0.1% penicillin-streptomycin (Gibco, #15140122). Cells were cultured at 37 °C, 5% CO₂. Mycoplasma testing was routinely performed. To generate Cas9 knock-in CT2A and GL261 cell lines, cells were transduced with Lenti-Cas9-2A-Blast (Addgene, #73310) and selected with Blasticidin S HCl (Gibco, #A1113903) as previously described²⁵. CT2A-Ova cell lines were generated via transduction with lenti-Ova, and sorted for high expression.

NK and CTL primary cells. Primary NK and CTL cells were isolated from OT-1 (C57BL/6-Tg(TcraTcrb)1100Mjb/J mice (Jackson Laboratory). Spleens were harvested, minced, and strained to obtain splenocytes. NK cells were then isolated through negative selection using the NK Cell Isolation Kit (Miltenyi Biotec, #130-115-818). NK cells were cultured in 10% RPMI with 1% penicillin-streptomycin, 100 ng/mL IL-2 and 55 mM 2-mercaptoethanol. Primary CTL cells were isolated, cultured, and activated as described previously²⁵. In short, OT-1 CD8+ T cells were isolated using an antibody-based magnetic separation kit (Miltenyi, #130-096-543), and activated and expanded with CD3/CD28 beads (Myltenyi, #130-093-627).

Myeloid cell lines. RAW264.7 cells were purchased from MilliporeSigma, BV-2 cells from AcceGen Biotech, and J774A.1 cells from ATCC. RAW264.1 and J774A.1 were cultured in DMEM (Wisent, #319-005-CL) supplemented with 10% FBS (Gibco, #12483-020). BV-2 cells were cultured in RPMI-1640 (Wisent, #350-000-CL) supplemented with 10% FBS. Cell lines were maintained in humidified incubators at 37 °C and 5% CO₂ and were routinely tested for mycoplasma contamination.

CRISPR-Cas9-edited cell lines. CRISPR-mediated gene knockouts in CT2A cell lines were generated by electroporation using the Neon Transfection System (Invitrogen, MPK10096) following the manufacturer’s instructions. In brief, Cas9 ribonucleoproteins (RNP) were prepared by combining 20 pmol single-guide (sg) RNA (Synthego) consisting of 3 sgRNA’s targeting the same gene with 20 pmol s.p. Cas9 nuclease (IDT #1081059) in 5 μL buffer R (Invitrogen) for 15 min at room temperature. sgRNA sequences are listed in Table S11. CT2A cells were lifted with trypsin (Gibco) and washed twice with Dulbecco’s phosphate buffered saline (DPBS, Wisent) and resuspended with buffer R. 200,000 cells were combined with 20 pmol Cas9 RNP in 12 μL buffer R and electroporated at 1200 V, for 2 pulses of 30 ms with 10 μL tips. 72 h after electroporation, cells were subjected to limiting dilution and single-cell clonal expansion. Genomic DNA from selected clones was extracted using the DNA Fast Extract kit (Wisent, #801-200-DR) and sgRNA target regions were PCR amplified with DNA 2X HS-Red Taq PCR mastermix (Wisent, #801-200-DM). Primer sequences are provided in Table S12. The PCR product was sequenced by Sanger sequencing. Confirmation of gene knockout was performed using ICE (https://ice.synthego.com/#/) to identify out-of-frame insertion-and-deletion mutations. CRISPR-mediated gene knockouts were also verified by western blot. Antibodies used are listed in Table S13.

Immunoblot analysis. To confirm gene knockouts in CT2A cells, cells were cultured to 90% confluency in 10 cm dishes. The cells were washed with DPBS and lysed with RIPA buffer (Thermo Scientific, #89901) supplemented with 1X protease inhibitor (Thermo Scientific, #78420) at 4°C. Protein quantification was done by Pierce BCA Assay (Life Technologies, #23225). Lysates were loaded on precast SDS-PAGE gels (Invitrogen, NP0321PK2) and subsequently transferred onto nitrocellulose membrane for detection. All primary antibodies were probed overnight at 4°C, and membranes were washed with Tris-buffered saline with Tween-20 (TBST; Cell Signaling Technology, #9997) and incubated with appropriate HRP-conjugated secondary antibodies for 1 h. Subsequently membranes were washed with TBST and the signal was detected with chemiluminescent substrate (Thermo Scientific, #34579) on an iBright Imaging system (Thermo Fisher Scientific).

Cell Proliferation Assay. Single cells were plated in a 96-well plate at a density of 1,000 or 200 cells/200 μL per well and incubated at 37°C and 5% CO2 for 5 days. 20 μL of Presto Blue (ThermoFisher, Cat.A13262), a fluorescent cell metabolism indicator, was added to each well four hours prior to reading out the assy. Fluorescence was measured using a FLUOstar Omega Fluorescence 556 Microplate reader (BMG LABTECH) with excitation and emission wavelengths of 544 nm and 590 nm, respectively. Readings were then analyzed using Omega analysis software. Proliferation was calculated for each well by subtracting the average RFI of blank wells from the RFI of individual wells. Normalized RFI was plotted for each well being tested as a side-by-side comparison of proliferation.

Immune cell killing assays. The killing efficiency of immune cells cocultured with CT2A cells were assessed across a range of effector-to-tumor ratios (E:T) and quantified as a percentages killing relative to untreated conditions. 24 h prior to killing assay, CT2A cells were incubated in 1 mM carboxyfluorescein succinimidyl ester (CFSE) dye (37 ºC × 10-20 min) and then cultured in complete medium overnight. On day of killing assay, CT2A cells were re-plated with immune cells at various E:T ratios. At 24 h endpoint, immune cells and dead CT2A cells were removed by gentle PBS wash. The remaining viable and adherent CT2A cells were assessed visually and by counting using a Coulter counter. For ADCP-mediated killing assays, killing efficiency was assessed by flow cytometry as before²⁵. CFSE and anti-CD11a were used as CT2A and J774.1 cell markers, and double-positive cells were interpreted as phagocytosed CT2A cells. E:T ratios that achieved between 20-50% killing efficiency were selected for screening conditions.

Single cell transcriptomic analysis using sci-RNA-seq3

Sample processing, sci-RNA-seq3 library generation, and sequencing. Cells were harvested with 0.25% typsin-EDTA and neuron dissociation solution¹¹³, respectively. Cell pellets were immediately snap-frozen in liquid nitrogen and then stored at -80°C for sci-RNA-seq3 based single-nucleus RNA-Seq processing. Samples from all genotypes were processed together to minimize batch effects. Nuclei extraction and fixation were performed as previously described¹¹³, except for the use of a modified CST lysis buffer¹¹⁴ plus 1% of SUPERase In RNase Inhibitor (AM2696). Nuclei quality was checked with DAPI and Wheat Germ Agglutinin staining. Sci-RNA-seq3 libraries were generated as previously described using three-level combinatorial indexing¹¹³. The final libraries were sequenced on Illumina NovaSeq as follows: read 1: 34bp, read 2: 69bp, index 1: 10bp, index 2: 10bp.

Demultiplexing and read alignments. Raw sequencing reads were first demultiplexed based on i5/i7 PCR barcodes. FASTQ files were then processed using the sci-RNA-Seq3 pipeline¹¹³. After barcodes and unique molecular identifiers (UMIs) were extracted from the read1 of FASTQ files, read alignment was performed using STAR short-read aligner (v2.5.2b) with the mouse genome (mm10) and Gencode vM12 gene annotations. After removing duplicate reads based on UMI, barcode, chromosome and alignment position, reads were summarized into a count matrix of M genes × N nuclei.

Filtering. Raw count matrices were loaded into a Seurat object (version 4.0.1) and filtered to retain cells with (i) 200 – 9000 recovered genes per cell, (ii) less than 60% mitochondrial content, and (iii) unmatched rate between 0.11 to 0.27 (median ± 3 median absolute deviations).

Normalization. To normalize expression values, we adopted the modeling framework previously described and implemented in sctransform (R Package, version 0.3.2)¹¹⁵. In brief, count data were modelled by regularized negative binomial regression, using sequencing depth as a model covariate to regress out the influence of technical effects, and Pearson residuals were used as the normalized and variance stabilized biological signal for downstream analysis.

Integration. Data from each timepoint and replicate were integrated with the reciprocal PCA method (Seurat) using the top 2000 variable features.

Dimensional reduction. PCA was performed on the integrated dataset, and the top N components that accounted for 90% of the observed variance were used for UMAP embedding, RunUMAP(max_components = 2, n_neighbours = 50, min_dist = 01, metric = cosine).

Clustering and annotation. To identify cellular sub-populations, we performed clustering using the Louvain algorithm implemented in Seurat (resolution = 1.5). Cluster were assigned numerical identifiers that were ranked in descending order according to cluster size, such that the smallest identifier (i.e., 0) corresponding to the largest cluster size, and vice versa. Clustered populations were then annotated using spatial and scRNA-seq reference atlases. Using Seurat’s label transfer pipeline, we identiﬁed transfer anchors for each query-reference pair using FindTransferAnchors(normalization.method = ‘‘SCT’’, reference.reudction = ‘‘pca’’, dims = 1:50) and then mapped the query samples onto the reference atlas using MapQuery(reference.reduction = ‘‘pca’’, reduction.model = ‘‘umap’’). The resulting prediction scores were cross-validated with cluster-speciﬁc markers obtained from differential-expression analyses to inform cluster annotation.

Differential-expression analysis. Differential expression analyses were performed using the Wilcoxon [wilcoxauc(…) function, (presto R package, v 1.0.0)] and co-dependency index (CDI) methods [FindCDIMarkers(…) function, (scMiko R package, v. 1.0.0)]. Differentially expressed genes (thresholded at 5% FDR) were ranked by area under receiver operating characteristic curve (AUROC) or normalized CDI (nCDI) scores and the top 2-4 markers for each cell type cluster were used for visualization.

NMF gene program discovery. For each sample, NMF was performed across multiple rank parameters and gene programs that were consistently resolved within and between tumors across multiple ranks were retained as robust NMF programs. NMF workflow parameters are summarized in Table S11. The workflow described here was modified from work by Gavish and colleagues:⁵⁵

Expression matrix preparation: For each sample, genes expressed in 0.5-5% of cells were used for NMF analysis. Expression matrices were normalized and scaled using the sctransform workflow (see above), and negative values were truncated at zero to obtain a non-negative scaled expression matrix (NSM).
Run NMF analysis: NMF analyses was performed for each NSM using nnmf(…, k = c(2,...,15), loss = “mse”, rel.tol = 1e-4, max.iter = 50) (NNLM R package, v 0.4.4). NMF was performed for rank parameters between k = 2-15 (Table S14). For each NMF run, we defined the resulting gene programs as the top 70-150 genes, ordered by factor loading.
Identify component programs: We reasoned that certain NMF runs will yield non-informative decompositions, and that only a subset of NMF programs will be reproducible within and between samples. To ensure within-sample robustness, the Jaccard similarity between NMF programs was computed across all intra-sample ranks, and programs with similarities exceeding 0.2-0.7 across more than 2 ranks were retained. Next, to ensure between-sample robustness, the Jaccard similarities between the remaining NMF programs were computed between samples, and programs with similarities exceeding 0.15-0.5 in at least 2 samples were retained. The resulting NMF programs that passed the intra- and inter-sample filtering steps were termed component programs.
Identify consensus programs: To obtain a non-redundant set of gene programs from the component programs, a Jaccard similarity matrix was computed for the remaining component programs, and hierarchically-clustered using Pearson correlation as the distance metric. Clustering parameters that separated the component programs were determined by visual inspection of the hierarchically-clustered Jaccard similarity heatmap. For each cluster, the consensus program was defined as the set of genes that were represented in at least 25-50% of component programs (Table S14).

Gene set enrichment analysis. To functionally-annotated gene sets, we performed hypergeometric overrepresentation analysis using the fora(…) (fgsea R package, v 1.14.0). Annotated gene sets used for enrichment analyses included GO ontology (biological processes, cellular components, molecular function) and those curated by the Bader Lab (http://baderlab.org/GeneSets).

Enrichment Network. Jaccard similarity between enriched pathways was used to account for gene set redundancies. Functional enrichment co-similarities were visualized using emapplot(…, edge.params = list(min = 0.2)) (enrichplot R package, v1.18.1).

Gene program activity. For each gene set, cell-level gene program activities were calculated using AddModuleScore(…) in Seurat. Where indicated, the programs were classified as “on” or “off” states by performing gaussian decomposition of gene program activity values. This was implemented using MClust(…, G = 1:5) (mclust R package, v 6.0.0) where parameterised finite gaussian mixture models were fit to a numeric vector of activity values and models were estimated by expectation-maximization algorithm initialized by hierarchical model-based agglomerative clustering. For each program, the optimal model was selected according to Bayesian information criterion (BIC).

Glioma-associated transcriptomic regulators. For each GBM scRNAseq dataset (overview of datasets in Table S10), cell-level gene program activities were calculated using i) a curated list of GBM- and tumor-associated gene panels and ii) TF-specific target gene sets consolidated from TRRUST v2¹¹⁶ and AnimalTFDB 3.0¹¹⁷. Machine-learning-based random forest regression models were then trained to identify which transcription regulators were associated with each GBM subtype. This was repeated across each GBM dataset, and the top transcription factors that were consistently associated with a GBM subtype were nominated as GBM glioma-associated transcription regulators (GTRs, 5% FDR). For each GBM transcriptomic phenotype, the top regulatory transcription factors were visualized in a bipartite network, where one layer of nodes represented GBM-specific gene programs and the other layer of nodes represented associated GTRs.

Relative abundance analysis. To compare the relative abundance of each cell type across conditions, counts for each cell type were tallied within each condition, and divided by the total number of cells profiled.

Differential abundance analysis. To evaluate regional differential abundances of cells in UMAP space across genotypes, we adopted the Milo method¹¹⁸. In brief, for each comparison between the CT2A and GL261-engrafted mouse samples, cells were first resampled to normalize cell counts, and then a K-nearest neighbor (KNN) graph representing higher-dimensional relationships between single cells was constructed. The KNN graph was then used to define neighborhoods of cells using the refined sampling scheme¹¹⁸. Finally, the number of cells belonging to each condition within each neighborhood was counted and the differential abundance was computed using a negative binomial generalized linear model. The differential abundance estimates for each neighborhood were visualized in UMAP space, with each node representing a given neighborhood (comprised of 20-80 cells each), and the color representing the differential abundance expressed as log fold-change between GL261 and CT2A samples. Non-significant differentials (FDR > 0.1) were truncated at zero.

Transcriptomic similarity. To evaluate the relative transcriptomic similarity of various human gliomas to murine gliomas, we adopted the transfer learning approach implemented in Seurat¹¹⁹. Specifically, murine glioma cells were mapped onto human GBM transcriptomic space by projecting the query (murine glioma cells) PCA structure onto the reference (human GBM) PCA structure. This enabled identification of corresponding cells (i.e., anchors), thereby allowing the mapping of murine gliomas cells onto the human GBM space and inference of the corresponding human glioma label in murine cells. The transfer scores computed using the FindTransferAnchor(…, normalization.method = “SCT”, k.anchor = 5, n.trees = 100, npcs = 40) and TransferData(…) functions in Seurat were used as surrogate measures of transcriptomic similarity.

Survival analyses

Murine glioma programs survival associations. The prognostic value of NMF gene programs identified in murine glioma cells was assessed by performing survival analyses using TCGA RNA-seq profiles of GBM patient samples (N=162). For each sample, NMF program activity was calculated by ssGSEA, and survival associations were assessed constructing Kaplan-Meier curves (survminer R package, v 0.4.9) in which samples were stratified into low (score <0) and High (score ≥0) groups, split based on median score. Corresponding hazard ratios were estimated by fitting proportional hazards regression models using the coxph(…) function (survival R package, v 3.4-0).

Meta-analysis of E1 and E2 survival associations. Using RNA-seq data from each GBM cohort, E1 and E2 program activities were scored using ssGSEA. Then, E1 and E2 hazard ratios (HR) were calculated using proportional hazards regression. Meta-analytic HR estimates were estimated by fitting a random effects model using the rma(…, method = “REML”) function (metafor R package, v 4.0) and forest plots were generated using forest.rma(…). I², a measure of heterogeneity, was quantified to estimate the percentage of total variance that is attributed to interstudy variance.

Genome-wide CRISPR screens

Pooled genome-wide CRISPR screens in CT2A and GL261 cells. CRISPR screens in sTable CT2A-Cas9 and GL261-Cas9 cells were performed as previously described. In brief, 150×10⁶ cells were infected with the mTKO lentiviral library (Addgene #159393) at an MOI of around 0.3. 24 h after infection, culture medium was changed to puromycin-supplemented medium (5 µg/mL). 72 h after infection, 1×10⁸ cells were cryo-banked whereas 9×10⁷ cells were split into three replicates of 3×10⁷ cells and further passaged every 3–4 d while maintaining 200-fold library coverage. 3×10⁷ cells were collected for genomic DNA extraction at day 0 (T₀) post-selection and at every subsequent passage until days 19 (T₁₉, CT2A) or 15 (T₁₅, GL261) post-selection. Genomic DNA extraction, library preparation, and sequencing were performed as described below.

CD29-sort screen in CT2A cells. Infection, selection and passaging of CT2A-Cas9 cells for the CD29 sort screen was performed as described above. At T₁₉, cells were stained with anti-CD29-FITC antibody (20 mg/1000 mL) and FACS sorting was performed as before¹²⁰. The unsorted T₁₉ sample was used as a reference. Genomic DNA extraction, library preparation, sequencing, and data processing were performed as described below.

Genome-wide CRISPR loss-of-function immune killing screens. Cas9-expressing CT2A-Ova cells were infected with the mTKO lentiviral library at a multiplicity of infection of 0.3 and maintained at 200-fold-coverage for each gRNA in the library. Infected cells were selected with 5 µg/mL of puromycin for 48 h. After selection cells were split in technical triplicates and maintained in culture for 72 h. For myeloid killing screens, CT2A cells (3×10⁶)were harvested and co-cultured with LPS-polarized (100 ng/mL LPS × 24 h) RAW264.7 (1.5×10⁵; 0.05 E:T), J774A.1 (3×10⁵; 0.1 E:T), and BV-2 (6×10⁵; 0.2 E:T) myeloid cell lines for 24 h to achieve killing efficiencies ranging between 22.9% to 44.1% (Fig S12A). For the J774A.1 coculture screen, an additional treatment arm was included in which CT2A cells were opsonized with anti-CD29 antibody (37ºC × 20 min) followed by anti-Armenian Hamster IgG2b antibody (37ºC × 20 min, then room temperature × 40 min) prior to coculturing with J774A.1 cells to facilitate ADCP. Anti-CD29 antibody was confirmed to not affect CT2A viability (Fig S12B). ADCP-dependent killing efficiency was 33.7% across replicates (Fig S12C). After 24 h of selective pressure, myeloid cells were eliminated using puromycin (8 µg/mL × 48 h) and confirmed by CD11b staining. For CTL killing screens, CT2A cells were co-cultured with preactivated CD8+ T cells at a 0.2 E:T ratio to achieve 46.2% killing efficiency (Fig S12D, see Lawson et al. for further details²⁵). For NK killing screens, a 0.5 E:T ratio was used to achieve 43.8% killing efficiency (Fig S12E). Untreated CT-2A cells were kept in parallel as a control. Across all conditions, killing efficiencies was calculated as the number of CT-2A cells in co-culture relative to the number of unchallenged cells. Cell pellets at 200-fold library coverage were collected for genomic DNA extraction at the end of screen. Genomic DNA extraction, library preparation, sequencing, and data processing were performed as described below.

Genomic DNA extraction, library preparation, sequencing, and data processing. For all screens, genomic DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega, #A1120)¹²¹. Sequencing libraries were prepared, sequenced using Illumina HiSeq2500, and processed as before^25,122. Bayes factor (BF) scores were calculated using BAGEL for CT2A and GL261 dropout screens, such that genes with a BF > 5 were classified as essential genes¹²³. Scaled BF was calculated as scaled BF = BF - 5. For immune coculture screens, fold changes between immune cell-treated and untreated CT2A cells were calculated using limma¹²⁴.

Human fitness scores. Processed CRIPSR-Cas9 genome screen data was retrieved from Project Score database for 41 human GBM cell-lines and 1031 human non-CNS cell lines^41,42. For each cell-line, gene-level fitness was represented as a scaled Bayes factor (BF), which corresponded to BAGEL2-derived BF subtracted by the BF at the 5% FDR threshold³¹. Scaled BFs were then pooled as averages for GBM and non-CNS cell lines and genes with BF > 0 were classified as essential within each respective group. If genes were essential in both groups, they were additionally termed common essentials, and if genes were essential in one group but not the other, they were termed either GBM-specific or Non-CNS-specific. All other genes were non-essentials.

Data visualization and statistics

Data visualization. Unless otherwise specified, the ggplot2 R package (v 3.3.5) was used for data visualization. scRNA-seq gene expression was visualized using FeaturePlot function (Seurat) or DotPlot function (Seurat). Venn diagrams were generated using either ssvFeatureEuler (seqsetvis R package, v 1.8.0) or ggVennDiagram (ggVennDiagram R package, v 1.1.4).

Statistics and reproducibility. All pairwise comparisons were performed using the signed Wilcoxon rank sum test, and p values were adjusted for multiple comparisons using the Benjamini-Hochberg procedure, as indicated. In cases where methods were compared across a common set of data, paired Wilcoxon tests were performed.

Data and Code Availability

R scripts used to perform the analyses are provided on GIT repository (https://github.com/NMikolajewicz/Mikolajewicz-2024). Sci-RNA-seq3 data is available on FigShare (10.6084/m9.Figshare.25685523).

CONFLICT OF INTEREST STATEMENT

The research was conducted in the absence of any commercial/financial relationships that could be construed as a conflict of interest.

AUTHOR CONTRIBUTIONS

Conceptualization: N.M., H.H., S.S., J.M.; Methodology: N.M., J.W., H.H., V.D., N.T., J.M.; Data Acquisition: J.W., N.S., A.G.F., V.D., K.D., M.A.U., Y.X., S.Y.L., N.T., C.V., H.H., V.D.; Analysis and Interpretation: N.M., V.D., D.C., Z.Z., K.B., H.H.; Visualization: N.M., D.C.; Drafting of Manuscript: N.M.; Funding and Supervision: S.S., J.M. All authors contributed to the critical revision and approval of the final manuscript.

FUNDING

This research work was supported by the 2020 William Donald Nash Brain Tumor Research Fellowship awarded to N.M. and the Canadian Institutes for Health Research (PJT438232 to J.M.). S.S. is a Tier 1 Canada Research Chair in Human Cancer Stem Cell Biology. J.M. is the GlaxoSmithKline Chair in Genetics and Genome Biology at the Hospital for Sick Children.

Jubelin C, Munoz-Garcia J, Griscom L, et al. Three-dimensional in vitro culture models in oncology research. Cell Biosci. 2022; 12(1):155.
Louis DN, Perry A, Wesseling P, et al. The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro Oncol. 2021; 23(8):1231-1251.
Horbinski C, Berger T, Packer RJ, Wen PY. Clinical implications of the 2021 edition of the WHO classification of central nervous system tumours. Nat Rev Neurol. 2022; 18(9):515-529.
Siah KW, Xu Q, Tanner K, Futer O, Frishkopf JJ, Lo AW. Accelerating glioblastoma therapeutics via venture philanthropy. Drug Discov Today. 2021; 26(7):1744-1749.
Singh K, Hotchkiss KM, Parney IF, et al. Correcting the drug development paradigm for glioblastoma requires serial tissue sampling. Nat Med. 2023; 29(10):2402-2405.
Ren AL, Wu JY, Lee SY, Lim M. Translational Models in Glioma Immunotherapy Research. Curr Oncol. 2023; 30(6):5704-5718.
Marsh J, Mukherjee P, Seyfried TN. Akt-dependent proapoptotic effects of dietary restriction on late-stage management of a phosphatase and tensin homologue/tuberous sclerosis complex 2-deficient mouse astrocytoma. Clin Cancer Res. 2008; 14(23):7751-7762.
Khalsa JK, Cheng N, Keegan J, et al. Immune phenotyping of diverse syngeneic murine brain tumors identifies immunologically distinct types. Nat Commun. 2020; 11(1):3912.
Iorgulescu JB, Ruthen N, Ahn R, et al. Antigen presentation deficiency, mesenchymal differentiation, and resistance to immunotherapy in the murine syngeneic CT2A tumor model. Front Immunol. 2023; 14:1297932.
Martinez-Murillo R, Martinez A. Standardization of an orthotopic mouse brain tumor model following transplantation of CT-2A astrocytoma cells. Histol Histopathol. 2007; 22(12):1309-1326.
Szatmari T, Lumniczky K, Desaknai S, et al. Detailed characterization of the mouse glioma 261 tumor model for experimental glioblastoma therapy. Cancer Sci. 2006; 97(6):546-553.
Johanns TM, Ward JP, Miller CA, et al. Endogenous Neoantigen-Specific CD8 T Cells Identified in Two Glioblastoma Models Using a Cancer Immunogenomics Approach. Cancer Immunol Res. 2016; 4(12):1007-1015.
Pellegatta S, Valletta L, Corbetta C, et al. Effective immuno-targeting of the IDH1 mutation R132H in a murine model of intracranial glioma. Acta Neuropathol Commun. 2015; 3:4.
Liu CJ, Schaettler M, Blaha DT, et al. Treatment of an aggressive orthotopic murine glioblastoma model with combination checkpoint blockade and a multivalent neoantigen vaccine. Neuro Oncol. 2020; 22(9):1276-1288.
Wainwright DA, Chang AL, Dey M, et al. Durable therapeutic efficacy utilizing combinatorial blockade against IDO, CTLA-4, and PD-L1 in mice with brain tumors. Clin Cancer Res. 2014; 20(20):5290-5301.
McKelvey KJ, Hudson AL, Donaghy H, et al. Differential effects of radiation fractionation regimens on glioblastoma. Radiat Oncol. 2022; 17(1):17.
McKelvey KJ, Wilson EB, Short S, et al. Glycolysis and Fatty Acid Oxidation Inhibition Improves Survival in Glioblastoma. Front Oncol. 2021; 11:633210.
Wu S, Calero-Perez P, Arus C, Candiota AP. Anti-PD-1 Immunotherapy in Preclinical GL261 Glioblastoma: Influence of Therapeutic Parameters and Non-Invasive Response Biomarker Assessment with MRSI-Based Approaches. Int J Mol Sci. 2020; 21(22).
Kijima N, Kanemura Y. Mouse Models of Glioblastoma. In: De Vleeschouwer S, ed. Glioblastoma. Brisbane (AU)2017.
Seyfried TN, el-Abbadi M, Roy ML. Ganglioside distribution in murine neural tumors. Mol Chem Neuropathol. 1992; 17(2):147-167.
Seligman AM, Shear MJ, Alexander L. Studies in Carcinogenesis: VIII. Experimental Production of Brain Tumors in Mice with Methylcholanthrene1. The American Journal of Cancer. 1939; 37(3):364-395.
Mukherjee P, Abate LE, Seyfried TN. Antiangiogenic and proapoptotic effects of dietary restriction on experimental mouse and human brain tumors. Clin Cancer Res. 2004; 10(16):5622-5629.
Ausman JI, Shapiro WR, Rall DP. Studies on the chemotherapy of experimental brain tumors: development of an experimental model. Cancer Res. 1970; 30(9):2394-2400.
Zagzag D, Amirnovin R, Greco MA, et al. Vascular apoptosis and involution in gliomas precede neovascularization: a novel concept for glioma growth and angiogenesis. Lab Invest. 2000; 80(6):837-849.
Lawson KA, Sousa CM, Zhang X, et al. Functional genomic landscape of cancer-intrinsic evasion of killing by T cells. Nature. 2020; 586(7827):120-126.
Cao W, Xi X, Hao Z, et al. RAET1E2, a soluble isoform of the UL16-binding protein RAET1E produced by tumor cells, inhibits NKG2D-mediated NK cytotoxicity. J Biol Chem. 2007; 282(26):18922-18928.
MacLeod G, Bozek DA, Rajakulendran N, et al. Genome-Wide CRISPR-Cas9 Screens Expose Genetic Vulnerabilities and Mechanisms of Temozolomide Sensitivity in Glioblastoma Stem Cells. Cell Rep. 2019; 27(3):971-986 e979.
Martin BK, Qiu C, Nichols E, et al. Optimized single-nucleus transcriptional profiling by combinatorial indexing. Nat Protoc. 2023; 18(1):188-207.
Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell. 2021; 184(13):3573-3587 e3529.
Zeisel A, Hochgerner H, Lonnerberg P, et al. Molecular Architecture of the Mouse Nervous System. Cell. 2018; 174(4):999-1014 e1022.
La Manno G, Siletti K, Furlan A, et al. Molecular architecture of the developing mouse brain. Nature. 2021; 596(7870):92-96.
Ochocka N, Segit P, Walentynowicz KA, et al. Single-cell RNA sequencing reveals functional heterogeneity of glioma-associated brain macrophages. Nat Commun. 2021; 12(1):1151.
Liu B, Li C, Li Z, Wang D, Ren X, Zhang Z. An entropy-based metric for assessing the purity of single cell populations. Nat Commun. 2020; 11(1):3155.
Yang J, Horton JR, Li J, et al. Structural basis for preferential binding of human TCF4 to DNA containing 5-carboxylcytosine. Nucleic Acids Res. 2019; 47(16):8375-8387.
Wittmann MT, Katada S, Sock E, et al. scRNA sequencing uncovers a TCF4-dependent transcription factor network regulating commissure development in mouse. Development. 2021; 148(14).
Mikheeva SA, Funk CC, Horner PJ, Rostomily RC, Mikheev AM. Novel TCF4:TCF12 heterodimer inhibits glioblastoma growth. Mol Oncol. 2024; 18(3):517-527.
Abdelfattah N, Kumar P, Wang C, et al. Single-cell analysis of human glioma and immune cells identifies S100A4 as an immunotherapy target. Nat Commun. 2022; 13(1):767.
Yu K, Hu Y, Wu F, et al. Surveying brain tumor heterogeneity by single-cell RNA-sequencing of multi-sector biopsies. Natl Sci Rev. 2020; 7(8):1306-1318.
Neftel C, Laffy J, Filbin MG, et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell. 2019; 178(4):835-849 e821.
Verhaak RG, Hoadley KA, Purdom E, et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010; 17(1):98-110.
Dwane L, Behan FM, Goncalves E, et al. Project Score database: a resource for investigating cancer cell dependencies and prioritizing therapeutic targets. Nucleic Acids Res. 2021; 49(D1):D1365-D1372.
Pacini C, Dempster JM, Boyle I, et al. Integrated cross-study datasets of genetic dependencies in cancer. Nat Commun. 2021; 12(1):1661.
McGranahan T, Therkelsen KE, Ahmad S, Nagpal S. Current State of Immunotherapy for Treatment of Glioblastoma. Curr Treat Options Oncol. 2019; 20(3):24.
Noffsinger B, Witter A, Sheybani N, et al. Technical choices significantly alter the adaptive immune response against immunocompetent murine gliomas in a model-dependent manner. J Neurooncol. 2021; 154(2):145-157.
Cui A, Huang T, Li S, et al. Dictionary of immune responses to cytokines at single-cell resolution. Nature. 2024; 625(7994):377-384.
Fortuna D, Hooper DC, Roberts AL, Harshyne LA, Nagurney M, Curtis MT. Potential role of CSF cytokine profiles in discriminating infectious from non-infectious CNS disorders. PLoS One. 2018; 13(10):e0205501.
Wang L, Jung J, Babikir H, et al. A single-cell atlas of glioblastoma evolution under therapy reveals cell-intrinsic and cell-extrinsic therapeutic targets. Nat Cancer. 2022; 3(12):1534-1552.
Zhao J, Chen AX, Gartrell RD, et al. Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med. 2019; 25(3):462-469.
Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24(10):1550-1558.
Yoshihara K, Shahmoradgoli M, Martinez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612.
Lu H, Zheng LY, Wu LY, Chen J, Xu N, Mi SC. The immune escape signature predicts the prognosis and immunotherapy sensitivity for pancreatic ductal adenocarcinoma. Front Oncol. 2022; 12:978921.
McKeown SR. Defining normoxia, physoxia and hypoxia in tumours-implications for treatment response. Br J Radiol. 2014; 87(1035):20130676.
Hum NR, Sebastian A, Gilmore SF, et al. Comparative Molecular Analysis of Cancer Behavior Cultured In Vitro, In Vivo, and Ex Vivo. Cancers (Basel). 2020; 12(3).
García-Vicente L, Borja M, Tran V, et al. Single-nucleus RNA sequencing provides insights into the GL261-GSC syngeneic mouse model of glioblastoma. bioRxiv. 2023:2023.2010.2026.564166.
Gavish A, Tyler M, Simkin D, et al. The transcriptional hallmarks of intra-tumor heterogeneity across a thousand tumors. bioRxiv. 2021:2021.2012.2019.473368.
Kinker GS, Greenwald AC, Tal R, et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity. Nat Genet. 2020; 52(11):1208-1218.
Boj SF, van Es JH, Huch M, et al. Diabetes risk gene and Wnt effector Tcf7l2/TCF4 controls hepatic response to perinatal and adult metabolic demand. Cell. 2012; 151(7):1595-1607.
Kim SW, Kim SJ, Langley RR, Fidler IJ. Modulation of the cancer cell transcriptome by culture media formulations and cell density. Int J Oncol. 2015; 46(5):2067-2075.
Tang L, Li H, Gou R, et al. Endothelin-1 mediated high glucose-induced epithelial-mesenchymal transition in renal tubular cells. Diabetes Res Clin Pract. 2014; 104(1):176-182.
Iacobini C, Vitale M, Pugliese G, Menini S. Normalizing HIF-1alpha Signaling Improves Cellular Glucose Metabolism and Blocks the Pathological Pathways of Hyperglycemic Damage. Biomedicines. 2021; 9(9).
Rogalska A, Forma E, Brys M, Sliwinska A, Marczak A. Hyperglycemia-Associated Dysregulation of O-GlcNAcylation and HIF1A Reduces Anticancer Action of Metformin in Ovarian Cancer Cells (SKOV-3). Int J Mol Sci. 2018; 19(9).
Mikheev AM, Mikheeva SA, Trister AD, et al. Periostin is a novel therapeutic target that predicts and regulates glioma malignancy. Neuro Oncol. 2015; 17(3):372-382.
Richards LM, Whitley OKN, MacLeod G, et al. Gradient of Developmental and Injury Response transcriptional states defines functional vulnerabilities underpinning glioblastoma heterogeneity. Nat Cancer. 2021; 2(2):157-173.
Zhang J, Chen L, Han L, et al. EZH2 is a negative prognostic factor and exhibits pro-oncogenic activity in glioblastoma. Cancer Lett. 2015; 356(2 Pt B):929-936.
Ma L, Lin K, Chang G, et al. Aberrant Activation of beta-Catenin Signaling Drives Glioma Tumorigenesis via USP1-Mediated Stabilization of EZH2. Cancer Res. 2019; 79(1):72-85.
Gouaze-Andersson V, Gherardi MJ, Lemarie A, et al. FGFR1/FOXM1 pathway: a key regulator of glioblastoma stem cells radioresistance and a prognosis biomarker. Oncotarget. 2018; 9(60):31637-31649.
Lee Y, Kim KH, Kim DG, et al. FoxM1 Promotes Stemness and Radio-Resistance of Glioblastoma by Regulating the Master Stem Cell Regulator Sox2. PLoS One. 2015; 10(10):e0137703.
Chandra A, Jahangiri A, Chen W, et al. Clonal ZEB1-Driven Mesenchymal Transition Promotes Targetable Oncologic Antiangiogenic Therapy Resistance. Cancer Res. 2020; 80(7):1498-1511.
Glasgow SM, Zhu W, Stolt CC, et al. Mutual antagonism between Sox10 and NFIA regulates diversification of glial lineages and glioma subtypes. Nat Neurosci. 2014; 17(10):1322-1329.
Lee J, Hoxha E, Song HR. A novel NFIA-NFkappaB feed-forward loop contributes to glioblastoma cell survival. Neuro Oncol. 2017; 19(4):524-534.
Yu X, Wang M, Zuo J, et al. Nuclear factor I A promotes temozolomide resistance in glioblastoma via activation of nuclear factor kappaB pathway. Life Sci. 2019; 236:116917.
Mathur R, Wang Q, Schupp PG, et al. Glioblastoma evolution and heterogeneity from a 3D whole-tumor perspective. Cell. 2024; 187(2):446-463 e416.
Azzarelli B, Miravalle L, Vidal R. Immunolocalization of the oligodendrocyte transcription factor 1 (Olig1) in brain tumors. J Neuropathol Exp Neurol. 2004; 63(2):170-179.
Ohnishi A, Sawa H, Tsuda M, et al. Expression of the oligodendroglial lineage-associated markers Olig1 and Olig2 in different types of human gliomas. J Neuropathol Exp Neurol. 2003; 62(10):1052-1059.
Riemenschneider MJ, Koy TH, Reifenberger G. Expression of oligodendrocyte lineage genes in oligodendroglial and astrocytic gliomas. Acta Neuropathol. 2004; 107(3):277-282.
Mokhtari K, Paris S, Aguirre-Cruz L, et al. Olig2 expression, GFAP, p53 and 1p loss analysis contribute to glioma subclassification. Neuropathol Appl Neurobiol. 2005; 31(1):62-69.
Bulstrode H, Johnstone E, Marques-Torrejon MA, et al. Elevated FOXG1 and SOX2 in glioblastoma enforces neural stem cell identity through transcriptional control of cell cycle and epigenetic regulators. Genes Dev. 2017; 31(8):757-773.
Song WS, Yang YP, Huang CS, et al. Sox2, a stemness gene, regulates tumor-initiating and drug-resistant properties in CD133-positive glioblastoma stem cells. J Chin Med Assoc. 2016; 79(10):538-545.
Berezovsky AD, Poisson LM, Cherba D, et al. Sox2 promotes malignancy in glioblastoma by regulating plasticity and astrocytic differentiation. Neoplasia. 2014; 16(3):193-206, 206 e119-125.
Lin B, Madan A, Yoon JG, et al. Massively parallel signature sequencing and bioinformatics analysis identifies up-regulation of TGFBI and SOX4 in human glioblastoma. PLoS One. 2010; 5(4):e10210.
Luo C, Quan Z, Zhong B, et al. lncRNA XIST promotes glioma proliferation and metastasis through miR-133a/SOX4. Exp Ther Med. 2020; 19(3):1641-1648.
Wu J, Li R, Li L, et al. MYC-activated lncRNA HNF1A-AS1 overexpression facilitates glioma progression via cooperating with miR-32-5p/SOX4 axis. Cancer Med. 2020; 9(17):6387-6398.
Ueda R, Yoshida K, Kawakami Y, Kawase T, Toda M. Immunohistochemical analysis of SOX6 expression in human brain tumors. Brain Tumor Pathol. 2004; 21(3):117-120.
Schlierf B, Friedrich RP, Roerig P, Felsberg J, Reifenberger G, Wegner M. Expression of SoxE and SoxD genes in human gliomas. Neuropathol Appl Neurobiol. 2007; 33(6):621-630.
Jiang YW, Wang R, Zhuang YD, Chen CM. Identification and validation of potential novel prognostic biomarkers for patients with glioma based on a gene co-expression network. Transl Cancer Res. 2020; 9(10):6444-6454.
Fu H, Cai J, Clevers H, et al. A genome-wide screen for spatially restricted expression patterns identifies transcription factors that regulate glial development. J Neurosci. 2009; 29(36):11399-11408.
Zhu G, Yang S, Wang R, et al. P53/miR-154 Pathway Regulates the Epithelial-Mesenchymal Transition in Glioblastoma Multiforme Cells by Targeting TCF12. Neuropsychiatr Dis Treat. 2021; 17:681-693.
Pang Y, Zhou S, Zumbo P, Betel D, Cisse B. TCF12 Deficiency Impairs the Proliferation of Glioblastoma Tumor Cells and Improves Survival. Cancers (Basel). 2023; 15(7).
Joseph JV, Conroy S, Tomar T, et al. TGF-beta is an inducer of ZEB1-dependent mesenchymal transdifferentiation in glioblastoma that is associated with tumor invasion. Cell Death Dis. 2014; 5(10):e1443.
Ding J, Adiconis X, Simmons SK, et al. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods. Nat Biotechnol. 2020; 38(6):737-746.
McKelvey KJ, Hudson AL, Prasanna Kumar R, et al. Temporal and spatial modulation of the tumor and systemic immune response in the murine Gl261 glioma model. PLoS One. 2020; 15(4):e0226444.
Kirschenbaum D, Xie K, Ingelfinger F, et al. Time-resolved single-cell transcriptomics defines immune trajectories in glioblastoma. Cell. 2024; 187(1):149-165 e123.
Khan SM, Desai R, Coxon A, et al. Impact of CD4 T cells on intratumoral CD8 T-cell exhaustion and responsiveness to PD-1 blockade therapy in mouse brain tumors. J Immunother Cancer. 2022; 10(12).
Duan S, Guo W, Xu Z, et al. Natural killer group 2D receptor and its ligands in cancer immune escape. Mol Cancer. 2019; 18(1):29.
Menasche BL, Davis EM, Wang S, et al. PBRM1 and the glycosylphosphatidylinositol biosynthetic pathway promote tumor killing mediated by MHC-unrestricted cytotoxic lymphocytes. Sci Adv. 2020; 6(48).
Sheffer M, Lowry E, Beelen N, et al. Genome-scale screens identify factors regulating tumor cell responses to natural killer cells. Nat Genet. 2021; 53(8):1196-1206.
Bernareggi D, Xie Q, Prager BC, et al. CHMP2A regulates tumor sensitivity to natural killer cell-mediated cytotoxicity. Nat Commun. 2022; 13(1):1899.
Kearney CJ, Vervoort SJ, Hogg SJ, et al. Tumor immune evasion arises through loss of TNF sensitivity. Sci Immunol. 2018; 3(23).
Kamber RA, Nishiga Y, Morton B, et al. Inter-cellular CRISPR screens reveal regulators of cancer cell phagocytosis. Nature. 2021; 597(7877):549-554.
Dhatchinamoorthy K, Colbert JD, Rock KL. Cancer Immune Evasion Through Loss of MHC Class I Antigen Presentation. Front Immunol. 2021; 12:636568.
Dubrot J, Du PP, Lane-Reticker SK, et al. In vivo CRISPR screens reveal the landscape of immune evasion pathways across cancer. Nat Immunol. 2022; 23(10):1495-1506.
Frey N, Tortola L, Egli D, et al. Loss of Rnf31 and Vps4b sensitizes pancreatic cancer to T cell-mediated killing. Nat Commun. 2022; 13(1):1804.
Spel L, Nieuwenhuis J, Haarsma R, et al. Nedd4-Binding Protein 1 and TNFAIP3-Interacting Protein 1 Control MHC-1 Display in Neuroblastoma. Cancer Res. 2018; 78(23):6621-6631.
Manguso RT, Pope HW, Zimmer MD, et al. In vivo CRISPR screening identifies Ptpn2 as a cancer immunotherapy target. Nature. 2017; 547(7664):413-418.
Varn FS, Johnson KC, Martinek J, et al. Glioma progression is shaped by genetic evolution and microenvironment interactions. Cell. 2022; 185(12):2184-2199 e2116.
Miranda A, Hamilton PT, Zhang AW, et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci U S A. 2019; 116(18):9020-9029.
Sese B, Iniguez-Munoz S, Ensenyat-Mendez M, et al. Glioblastoma Embryonic-like Stem Cells Exhibit Immune-Evasive Phenotype. Cancers (Basel). 2022; 14(9).
Qazi MA, Salim SK, Brown KR, et al. Characterization of the minimal residual disease state reveals distinct evolutionary trajectories of human glioblastoma. Cell Rep. 2022; 40(13):111420.
Zhao Z, Zhang KN, Wang Q, et al. Chinese Glioma Genome Atlas (CGGA): A Comprehensive Resource with Functional Genomic Data from Chinese Glioma Patients. Genomics Proteomics Bioinformatics. 2021; 19(1):1-12.
Barthel FP, Johnson KC, Varn FS, et al. Longitudinal molecular trajectories of diffuse glioma in adults. Nature. 2019; 576(7785):112-120.
Cloughesy TF, Mochizuki AY, Orpilla JR, et al. Neoadjuvant anti-PD-1 immunotherapy promotes a survival benefit with intratumoral and systemic immune responses in recurrent glioblastoma. Nat Med. 2019; 25(3):477-486.
Chokshi CR, Savage N, Venugopal C, Singh SK. A Patient-Derived Xenograft Model of Glioblastoma. STAR Protoc. 2020; 1(3):100179.
Cao J, Spielmann M, Qiu X, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature. 2019; 566(7745):496-502.
Slyper M, Porter CBM, Ashenberg O, et al. A single-cell and single-nucleus RNA-Seq toolbox for fresh and frozen human tumors. Nat Med. 2020; 26(5):792-802.
Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 2019; 20(1):296.
Han H, Cho JW, Lee S, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res. 2018; 46(D1):D380-D386.
Hu H, Miao YR, Jia LH, Yu QY, Zhang Q, Guo AY. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 2019; 47(D1):D33-D38.
Dann E, Henderson NC, Teichmann SA, Morgan MD, Marioni JC. Differential abundance testing on single-cell data using k-nearest neighbor graphs. Nat Biotechnol. 2022; 40(2):245-253.
Stuart T, Butler A, Hoffman P, et al. Comprehensive Integration of Single-Cell Data. Cell. 2019; 177(7):1888-1902 e1821.
Mair B, Aldridge PM, Atwal RS, et al. High-throughput genome-wide phenotypic screening via immunomagnetic cell sorting. Nat Biomed Eng. 2019; 3(10):796-805.
Ngo W, Wu JLY, Lin ZP, et al. Identifying cell receptors for the nanoparticle protein corona using genome screens. Nat Chem Biol. 2022; 18(9):1023-1031.
Chan K, Tong AHY, Brown KR, Mero P, Moffat J. Pooled CRISPR-Based Genetic Screens in Mammalian Cells. J Vis Exp. 2019(151).
Hart T, Moffat J. BAGEL: a computational framework for identifying essential genes from pooled library screens. BMC Bioinformatics. 2016; 17:164.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43(7):e47.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Functional profiling of murine glioma models highlights targetable immune evasion phenotypes

Status:

Version 1

Abstract

Figures

INTRODUCTION

RESULTS

Autophagy mediates intrinsic pan-immune evasion

Murine glioma cells partially recapitulate human genetic dependencies

Unbiased transcriptomic profiling of murine brain tumors

Glioma cells in vitro and in vivo have distinct transcriptomic signatures

GL261 and CT2A tumors recapitulate canonical GBM transcriptomic phenotypes

Human and murine gliomas have common transcriptional regulators

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

Intrinsic immune evasion of stem-like glioma cells

DISCUSSION

MATERIALS AND METHODS

Resources and Data sources

Animal studies

Cell lines

Single cell transcriptomic analysis using sci-RNA-seq3

Survival analyses

Genome-wide CRISPR screens

Data visualization and statistics

Declarations

Data and Code Availability

CONFLICT OF INTEREST STATEMENT

AUTHOR CONTRIBUTIONS

FUNDING

References

Additional Declarations

Supplementary Files

Status:

Version 1