OTME-23. Single-cell transcriptomic and epigenomic immune landscape of isocitrate dehydrogenase stratified human gliomas

Abstract The brain tumor immune microenvironment (TIME) continuously evolves during glioma progression, but a comprehensive characterization of the glioma-centric immune cell repertoire beyond a priori cell states is uncharted. In this study, we performed single-cell RNA-sequencing (scRNA-seq) and single cell- Assay for Transposase-Accessible Chromatin using sequencing (sc-ATAC-seq) on ~100,000 tumor-associated immune cells from seventeen isocitrate dehydrogenase (IDH) mutation classified primary and recurrent human gliomas and non-glioma brains (NGBs). Our analyses revealed sixty-two transcriptionally distinct myeloid and lymphoid cell states within and across glioma subtypes and we noted microglial attrition with increasing disease severity concomitant with invading monocyte-derived cells and lymphocytes. Specifically, certain microglial and monocyte-derived subpopulations were associated with antigen presentation gene modules, akin to cross-presenting dendritic cells (DCs). We identified cytotoxic T cells with poly-functional cytolytic states mostly in recurrent IDH-wt gliomas. Furthermore, ligand-receptor interactome analyses showed a preponderance of antigen presentation and phagocytosis over the checkpoint axis in IDH-wt compared to IDH-mut gliomas. Additionally, our sc-ATAC-seq analyses revealed differences in regulatory networks in NGBs, IDH-mut and IDH-wt glioma associated immune cells. In particular, we noted abundant usage of inflammatory transcription factors (TFs) as exemplified by Nuclear factor kappa B and Activator Protein-1 TF family in IDH-wt microglia when compared with microglia from IDH-mut and NGBs. Unique features such as amplification of 11- Zinc Finger Protein accessibility were restricted to monocyte derived cells and were not observed in microglia. Finally, sc-ATAC-seq profiles of CD8+ exhausted T cells from IDH-wt showed strong enhancer accessibility on Cytotoxic T-lymphocyte-associated protein 4, Layilin and Hepatitis A Virus Cellular Receptor 2 but no enrichment on PDCD1 (gene encoding Programmed cell death protein 1) was seen. In summary, our study provides unprecedented granular detail of transcriptionally defined glioma- specific immune contexture that can be exploited for immunotherapy applications. This study in K.B. laboratory was supported by the generous philanthropic contributions to The University of Texas (UT) MD Anderson Cancer Center (MDACC) Moon Shots Program™, Marnie Rose Foundation, NIH grants: R21 CA222992 and R01CA225963. This study was partly supported by the UT MDACC start-up research fund to L.W. and CPRIT Single Core grant RP180684 to N. E. N.

progression, but only a limited view of a highly complex glioma associated immune contexture 23 across isocitrate dehydrogenase mutation (IDH) classified gliomas is known. Herein, we present 24 an unprecedentedly comprehensive view of myeloid and lymphoid cell type diversity based on 25 our single cell RNA sequencing and spectral cytometry-based interrogation of tumor-associated 26 leukocytes from fifty-five IDH stratified primary and recurrent human gliomas and three non-27 glioma brains. Our analyses revealed twenty-two myeloid and lymphoid cell types within and 28 across glioma subtypes. Glioma severity correlated with microglial attrition concomitant with a 29 continuum of invading monocyte-derived microglia-like and macrophages amongst other 30 infiltrating conventional T and NK lymphocytes and unconventional mucosa associated invariant 31 T (MAIT) cells. Specifically, certain microglial and monocyte-derived subpopulations were 32 associated with antigen presentation gene modules, akin to cross-presenting dendritic cells 33 (DCs). Furthermore, we identified phagocytosis and antigen presentation gene modules 34 enriched in Triggering receptor expressed on myeloid (TREM)-2 + cells as a putative anti-glioma 35 axis. Accelerated glioma growth was observed in Trem2 deficient mice implanted with CT2A 36 glioma cells affirming the anti-glioma role of TREM2 + myeloid cells. In addition to providing a 37 comprehensive landscape of glioma-specific immune contexture, our investigations discover 38 TREM2 as a novel immunotherapy target for brain malignancies.

INTRODUCTION 40
The discovery of meningeal and parenchymal access of immune cells 1 and the presence of 41 meningeal 2-4 and dural lymphatics 5 in the central nervous system (CNS) has led to redefinition 42 of the brain being immunologically distinct rather than immune privileged; a notion held for 43 several decades 6 . In brain pathologies such as traumatic injury or neurodegenerative disorders, 44 phagocytic cells comprising brain resident microglia and CNS associated 45 monocytes/macrophages are the first responders to resolve associated inflammation. For 46 instance, in Alzheimer's disease the protective role of microglia in clearing amyloid plaques has 47 been established to prevent disease progression 7 . In contrast, tumor associated macrophages 48 have been linked to poor prognosis in brain neoplasms such as gliomas that arise from 49 captured the leukocyte diversity of the brain TIME in primary gliomas and brain metastasis 30,31 68 to a certain extent albeit based on a priori markers. 69 To address these knowledge gaps and delineate the glioma associated leukocyte 70 diversity in the TIME, we performed single cell (sc)-and bulk RNA Sequencing (RNA-seq) and 71 spectral cytometry analyses on tumor-associated leukocytes from fifty-five IDH-stratified primary (treatment naïve) and standard of care treated recurrent glioma subtypes to define their immune 73 landscape. In addition to corroborating established myeloid dominant glioma characteristics, we 74 redefine glioma TIME by superimposing our advanced findings across glioma subtypes with 75 largely IDH-wt glioma restrictive studies 17,18,22-24 . Our major findings include the following: i) We 76 observed significant attrition of microglia (MG) accompanied by increased infiltration of classical 77 monocytes (c-Mo), monocyte-derived-microglia-like (Mo-MG or MG-like), -macrophages (MDM) 78 and conventional dendritic cells (cDC)-2 in recurrent IDH-wild type gliomas relative to other 79 glioma subtypes; ii) We demonstrate eleven transcriptionally distinct glioma associated MG 80 states inclusive of tumoricidal, inflammatory and metabolic phenotypes; ii) Infiltration of Tregs, 81 NK cells and mucosa associated invariant T (MAIT) were significantly abundant in recurrent IDH-82 wild type gliomas; and iv) We identified glioma associated myeloid cells with triggering receptor 83 expressed on myeloid (TREM)-2 cells enriched for phagocytosis and antigen-presentation gene 84 modules as putative anti-glioma axis and demonstrate their anti-glioma functions using a 85 xenograft mouse model. In summary, our reverse translational glioma immunophenotyping 86 investigations reveal an unprecedently advanced landscape of glioma TIME that can be 87 exploited for future immunotherapy applications. We further uncover TREM2 as a novel glioma 88 specific immunomodulatory target with likely implications in other brain malignancies. 89

Spectrum of invading non-MG myeloid cells in human gliomas 197
In order to understand non-MG myeloid cell diversity, which have been a subject of intense 198 investigation in IDH-wt gliomas 12,13,18,20,21,33 , we performed a comprehensive analyses of 199 invading non-MG myeloid subpopulations across IDH-classified gliomas and identified sixteen 200 non-myeloid cell states (Fig. 3D). We uncovered six MAC and two MDM clusters. Based on DEG 201 and GO analysis we annotated glioma associated macrophages as IL10 + MAC_Anti-Inflam, 202 multiple metabolic phenotypes such as clusters enriched with hypoxia genes (e.g., SDS, 203 HMOX1) MAC_Metab/Hypoxia 1, MAC_Metab/Hypoxia 2 and LIPA + MAC_Lipid Metab 204 subpopulations. A cluster of LYVE1 expressing MAC associated with vasculature defined as 205 MAC_Perivascular was observed across glioma subtypes (Fig. 3E, Supplementary Fig. S2F-206 G). Amongst the MDMs, C1QA + MDM_Phagocytic 1 was associated with receptor mediated 207 endocytic process and JUN + SELENOP + MDM_Inflam subset correlated with positive regulation 208 of inflammation GO term (Fig. 3E, Supplementary Fig. S1M, S2F-G). Although MDMs showed 209 marked increase in IWR glioma subtype, abundance of MAC was proportionately similar across 210 IMP, IMR and IWP gliomas (Fig. 2D). Monocytic infiltrates in glioma TIME included significantly 211 abundant FCN1 + CD14 + c-Mo in IWR compared to IMP gliomas whereas FCGR3A + TCF7L2 + nc-212 Mo did not show any noticeable differences across other glioma subtypes (Fig. 3D, 2C, 213 Supplementary Fig. S1M, S2F). Glioma associated inflammation led to DC infiltration in 214 contrast to negligible DCs in NGB CD1C + cDC2 was significantly abundant in IWR glioma 215 subtypes, while a proportionally similar levels of infiltration of CLEC9A + cDC1 and IL3RA + pDC 216 were observed across other glioma subtypes (Fig. 3D, 2E, Supplementary Fig. S1M, S2F-G). 217 Taken together, our data provides advanced insight of infiltrating myeloid cell subsets in the 218 glioma TIME and highlights their similarities and differences with MG. 219 220

Identification of Trem2 as an anti-glioma modulator 221
To gain a deeper understanding of pathways altered in MG, we inspected genes that regulate 222 phagocytosis and antigen presentation in myeloid cells, especially MG_APC-Like 1 and 223 MG_APC-like 2 clusters (Fig. 4A). These clusters were enriched for HLA-DR (Fig. 4A), a major 224 histocompatibility complex (MHC) class II gene. In addition, we found higher expression of 225 TREM2 and its regulator MS4A6A (Fig. 4B) reduced survival compared to syngeneic WT mice with intracranial gliomas (Fig. 4C). These anti-glioma properties. Although TREM2 was shown to be immunosuppressive in other cancers, 283 TREM2 genetic ablation induced glioma growth in mice indicating differential functions of 284 TREM2 in brain tumors. Further investigation of TREM2 expression in brain resident and 285 infiltrating immune cells is warranted. 286 In summary, our unbiased high dimensional studies have paved the way for an advanced 287 understanding of immune landscape across gliomas that can be exploited for novel 288 immunotherapy strategies for these cancer types.

ACKNOWLEDGEMENTS 290
This study in K.B. laboratory was supported by the generous philanthropic contributions to The 291

University of Texas (UT) MD Anderson Cancer Center (MDACC) Moon Shots Program™, 292
Marnie Rose Foundation, NIH grants: R21 CA222992 and R01CA225963. This study was partly 293 supported by the UT MDACC start-up research fund to L.W. and CPRIT Single Core 294 grant RP180684 to N. E. N. We are thankful to Cynthia Kassab and Martina Ott for their 295 contribution to tissue logistics for scRNA-seq pipeline. In part, the patient specimens analyzed 296 with spectral cytometry in this study were contributed by CNS tumor Analysis Stream 297 (CATALYST) program at UT MDACC and we extend our gratitude to entire CATALYST team. 298 We thank Sanaalarab Al-Enazy for schematic illustrations. We would like to thank Haidee D. 299 Chancoco for her contribution to RNA isolation for bulk mRNA-seq at Biospecimen Extraction 300 Resource of MD Anderson Cancer Center; Advanced Technology Genomics Core supported by 301 NCI grant CA016672 for library preparation and sequencing; and David W. Dwyer and Karen C. All data and code available upon request 320

Conflict Statement 321
All authors declare no competing or financial interests METHODS 323 324

Human brain tumor and tissue collection 325
The brain tumor/tissue samples were collected from seventeen-patients post appropriate 326 informed consent during neurosurgery with detailed information pertaining to gender, age, 327 glioma grade, subtype and, brain site of tumor extraction etc. mentioned in Supplementary 328 Table S1. The brain tumor/tissue samples were collected as per MD Anderson internal review to the 10X protocol. Cleanup cDNA was checked via a tape station (Agilent 42000) HSD5000 379 (cat# part # 5067-5593) for cDNA traces. 25% of the cDNA was used to generate the library, 380 and the Chromium i7 multiplex kit (part #120262) was used to identify each sample. Library 381 cleanup was performed using AMPure beads and QC was done again with tape station D1000 382 tapes (Part # 5067-5583). Ten libraries of equal amount were pooled to give a final concentration 383 of 10nM and submitted for sequencing with the NovaSeq6000 S2 sequencer, 28 cycles for 384 read1, 8 cycles for i7 index, and 91 cycles for read 2 through the ATGC core at MD Anderson. 385 Sequence data was then put through the 10X genomic cell ranger 3.0 pipeline. QC and Fastq 386 files were obtained and checked for data quality, and Fastq files were used to do further analysis.
RNeasy Mini Kit (Cat. #74104, Qiagen) was used to extract sorted cells and achieve efficient 389 purification of total RNA from small amounts of starting material. The technology simplifies total 390 RNA isolation by combining the stringency of guanidine-isothiocyanate lysis with the speed and 391 purity of silica-membrane purification to ensure highest-quality RNA with minimum copurification 392 of DNA. With the RNeasy Mini Kit, total RNA can be purified from 10 to 1 x 10 7 animal or human 393 cells, 0.5-30 mg human tissues. Briefly, samples are first lysed and then homogenized. Ethanol 394 is added to the lysate to provide ideal binding conditions. The lysate is then loaded onto the 395 separately with default parameters. Each cell type was down sampled to 500 cells. We chose 445 the k input value from 1% to 100% of the sample size. In each run, the number of tested 446 neighborhoods was 10% of the sample size. The mean and maximal rejection rates were then 447 calculated based on a total of 100 repeated k-BET runs. Following estimation of sample 448 processing-or sequencing-related batch effects using k-BET, we employed Harmony for actual 449 batch effect correction 53 . Harmony was run with default parameters to remove batch effects in 450 the PCA space when clustering of major cell lineages before any clustering analysis or cell type 451 identification/annotation was performed. We carefully evaluated the performance of Harmony in 452 terms of its ability to integrate batches while maintaining cell type separation. Harmony was run 453 on all cells to firstly identify major cell types. It was also run on each of the major cell types for 454 subclustering analysis to further identify different cell states. To quantify the performance of 455 Harmony, we further used k-BET and compared the rejection rate (reflecting batch effect) before 456 and after Harmony. The data after Harmony showed a low rejection rate, indicating an excellent 457 performance of batch effect correction in this study. 458 Moreover, we also applied the local inverse Simpson's Index (LISI) to assess the performance 459 of Harmony. As described previously 53 , the 'integration LISI' (iLISI) measures the degree of 460 mixing among datasets (batches), ranging from 1 in an unmixed space to the number of datasets 461 (batches) in a well-mixed space. And the 'cell-type LISI' (cLISI) measures integration accuracy 462 using the same formulation but computed on cell-type labels instead. An accurate embedding 463 has a cLISI close to 1 for every neighborhood, reflecting separation of different cell types. Before 464 batch correction with Harmony, cells were mainly grouped by dataset (iLISI is around 1) and 465 cells from different cell types were mixed (cLISI is far from 1). After batch correction with 466 Harmony, iLISI and cLISI were re-computed in the Harmony embedding. iLISI is around 3.5, 467 indicating a high degree of mixing among different datasets, and cLISI is very close to 1, 468 reflecting excellent separation of different cell types while remain the well-mixed space.

Cytometry Analyses 550
The acquired data was analyzed by Cytek SpectroFlo and Becton Dickinson FlowJo 10.8.1.

Mouse experiments and survival analysis 552
Six-week-old WT and TREM2 -/mice were obtained from Jackson laboratories. One week after 553 guide screw implantation, 10,000 CT-2A cells were injected intracranially. Animals were 554 monitored for tumor growth using bioluminescence imaging. For survival analysis, we used the 555 log-rank test to calculate P values between groups, and the Kaplan-Meier method to plot survival 556 curves using GraphPad Prism9 version 9.2.0 software. 557