Tumor type and cell type-specific gene expression alterations in diverse pediatric central nervous system tumors identified using single nuclei RNA-seq

Central nervous system (CNS) tumors are the leading cause of pediatric cancer death, and these patients have an increased risk for developing secondary neoplasms. Due to the low prevalence of pediatric CNS tumors, major advances in targeted therapies have been lagging compared to other adult tumors. We collected single nuclei RNA-seq data from 35 pediatric CNS tumors and three non-tumoral pediatric brain tissues (84,700 nuclei) and characterized tumor heterogeneity and transcriptomic alterations. We distinguished cell subpopulations associated with specific tumor types including radial glial cells in ependymomas and oligodendrocyte precursor cells in astrocytomas. In tumors, we observed pathways important in neural stem cell-like populations, a cell type previously associated with therapy resistance. Lastly, we identified transcriptomic alterations among pediatric CNS tumor types compared to non-tumor tissues, while accounting for cell type effects on gene expression. Our results suggest potential tumor type and cell type-specific targets for pediatric CNS tumor treatment. In this study, we address current gaps in understanding single nuclei gene expression profiles of previously uninvestigated tumor types and enhance current knowledge of gene expression profiles of single cells of various pediatric CNS tumors.


Introduction
Central nervous system (CNS) tumors account for ~25% of pediatric cancer cases and are the leading cause of cancer death in children and adolescents in the United States 1 . Incident pediatric CNS tumors are comprised of many histologically distinct tumor types including pilocytic astrocytomas (15.2%), embryonal tumors (9.4%), and neuronal/mixed neuronal-glial tumors (7.9%) 2 . Survival rates vary widely among tumor types, with a good 10-year survival of 95.4% for pilocytic astrocytomas and a poor 10-year survival of 15.9% for pediatric high-grade gliomas 2 . Pediatric CNS tumor patients are at risk of developing secondary neoplasms, with a 30-year cumulative incidence of malignant secondary neoplasms ranging from 4.7 -7.8% 3,4 . The standard of care treatments for primary CNS tumors include surgery, radiotherapy, and chemotherapy with relatively limited options for targeted therapy compared to tumors in other anatomic regions.
Recent advances in identifying molecular subtypes in various pediatric CNS tumor types have been made utilizing genomic, transcriptomic and epigenomic data as reflected in the 2021 World Health Organization classification of CNS tumors 5 . For example, medulloblastoma can be classified into four separate molecularly defined subtypes: WNT-activated, SHH-activated and TP53-wildtype, SHHactivated and TP53-mutant, and non-WNT/non-SHH [6][7][8][9][10] . In addition, supratentorial ependymoma can be categorized into ZFTA fusion-positive or YAP1 fusion-positive 11,12 . A better understanding of the molecular variations that exist even among each tumor type has led to novel treatment options. For example, Larotrectinib and entrectinib, targeted therapies for NTRK fusion, which has been found in brain tumors, have been approved by the Food and Drug Administration to treat some brain tumors that are metastatic or unresectable with surgery 13,14 .
In addition to the molecular characterization of bulk pediatric CNS tumor tissue, emerging work has begun to investigate the transcriptome and cellular states that exist in these tumors at the single cell level. One of the first single cell transcriptomics contributions focused on H3K27M -altered pediatric gliomas (n=6, and 3,300 cells) showed that tumors are mainly composed of progenitor cell-like oligodendrocyte populations, rather than differentiated malignant cells 15 . Later, Gojo et al. identified that cellular hierarchies in primary ependymomas (n=28) reflect impaired neurodevelopment and that undifferentiated programs can infer prognosis 16  While these single cell and single nucleus transcriptomics studies in 85 total primary CNS tumors to date have improved our understanding of cell states in pediatric CNS tumors, there is still much to be investigated to advance optimal therapeutic options for both primary cancer treatment and reduction of secondary neoplasms. Due to limited sample availability for these rare pediatric CNS tumors, progress in single cell level characterization of these tumors has been relatively slow. Here, we characterized single nuclei gene expression profiles of 35 pediatric CNS tumors and 3 non-tumor pediatric brain tissues. Our study augments previous studies by incorporating single nuclei gene expression profiles of additional pediatric CNS tumor types (dysembryoplastic neuroepithelial tumors, gangliogliomas, etc.) and non-tumor pediatric brain tissue which have not yet been published to our knowledge.

Results
Samples from pediatric central nervous system tumors and non-tumor pediatric brain tissue were obtained from patients being treated at Dartmouth-Hitchcock Medical Center and Dartmouth Cancer Center from 1993 to 2017. Non-tumor pediatric brain tissues from the supratentorial regions were collected from patients undergoing surgical resection for epilepsy. Patient characteristics are described in Table 1. Pathological re-review for histopathologic tumor type and grade were done according to the 2021 World Health Organization CNS tumor classification system and categorized into the broader tumor types to balance sample size per tumor type 5 . Specific diagnoses for each sample can be found in Supplementary Table 1.
Genetic variants were identified using bulk tissue RNA-seq data for all tumors except for two tumors due to low bulk RNA-seq data quality. Copy number variations (CNV) were determined using bisulfite treated DNA methylation array data. Genetic and cytogenic variations varied among tumors and tumor types (Figure 1). Interestingly, across all but one tumor sample, tumors had genetic variants in MALAT1. Many of the genetic variants detected within the pediatric CNS tumors were associated with epigenetic processes. For example, almost half of the tumors, across tumor types, had genetic variants in HIST1H1E (14/33). CNV patterns in some tumor types were as expected from previous literature. For instance, 5 out of the 9 ependymoma had chromosome 1q gain, which has been considered to be an early tumorigenic event in ependymoma 18,19 .
Integrated de-multiplexing method to increase single nuclei RNA-seq data yield Using lipid-tagged hashtag oligonucleotides (HTO), 34 samples (out of 38 total samples) were multiplexed in 17 pools to collect 10X genomics snRNA-seq data 20 . The distribution of samples across sequencing runs and pools is provided in Supplementary Table 1. As many nuclei were not tagged with sufficient HTO to be efficiently demultiplexed in downstream analyses, we aimed to augment demultiplexing by analyzing sequencing-derived genotype data from each nucleus together with HTO information and assign additional nuclei to specific samples (Figure 2A). To summarize our demultiplexing process, we first used HTOs to assign the nuclei to their respective samples. For samples that were assigned confidently with the HTO, we filtered to keep only the nuclei that were assigned to the same sample concordantly using genotype information. For nuclei that were either unassigned to a sample or assigned as a doublet with HTO, we assigned nuclei to samples using genotype information (detailed in the methods section). The final set of nuclei per sample were comprised of the filtered nuclei from HTO and genotype identified nuclei. An example of the single nucleotide variants identified per pool along with their assigned sample can be found in Figure 2B. An example of how many nuclei were obtained for one pool, during each step, is shown in Figure 2A   nuclei) over the genotype-based method alone ( Figure 2C). Gene expression profiles for a total of 84,700 nuclei were used for downstream analyses.

Cell type heterogeneity in pediatric central nervous system tumors and non-tumor pediatric brains
Out of 84,700 nuclei, 67,249 nuclei (79%) were from pediatric CNS tumors and 17,451 nuclei (21%) were from non-tumor tissue ( Figure 3A). Across all samples, snRNA-seq data revealed 58 clusters that were grouped into 16 Figure   1A, Figure 3B). The clusters were classified into cell types using classical markers for cell types found in the brain ( Table 2). The following gene markers for cell types were used: GFAP and AQP4 for astrocytes; FN1 and COL4A1 for endothelial cells; CSF1R and PTPRC for macrophage/microglia;  (Supplementary Figure 2). Interestingly, cell types expected to be more differentiated, like astrocytes, had relatively high levels of CD44 and VIM, and these genes were expressed in many of the cell types. In addition, although the NSC-like cluster had a high stemness score, the expression of cancer stem cell markers was minimal. Unexpectedly, the UBC-like clusters also had elevated stemness scores. While gene expression levels may not always correlate with protein expression, our results indicate cell types identified using classical stem cell markers may not capture all tumor cells with stemness features.
Next, we tested for potential associations of clinical variables with tumor stemness scores. We first assessed the distribution of stemness scores among nuclei in each sample and determined the median stemness score (Supplementary Figure 3). We found that the stemness scores were higher in embryonal compared to other tumor types or non-tumor tissue (Supplementary Figure 4A).
Compared with low grade tumors, high grade tumors had higher stemness scores (P-value = 0.008,

Cell type-specific pathway enrichment in pediatric CNS tumors
First, to determine cell type-specific pathway enrichment in each nuclei of the tumor samples, we conducted a pathways analysis at the single cell level using the Variance-adjusted Mahalanobis (VAM) method, which computes cell-level pathway scores that account for the technical noise and inflated zero counts of single cell RNA-seq data 21 . We used 196 pathways from the MSigDB Pathway Interaction Database (PID) collection for our enrichment testing 23,24,31 . The cell-level enrichment pvalues generated by VAM were corrected for false discovery rate using the Benjamini-Hochberg method and classified to be significantly enriched in each nucleus if the FDR adjusted p-value was less than 0.1 as binary classifications (enriched or not enriched).
Next, we determined any pathways that were more specific for each cell type to determine pathways important in each cell type. The PID pathways were considered to be important/specific to the cell type under adjusted p-value < 0.05 threshold in the differential enrichment test. For cell types with a limited presence in tumor tissues, like many of the excitatory neurons and A1, we observed no pathways that were specific to the clusters (Supplementary Figure 6A, Supplementary Table 2). The immune-related cells (MG1, MG2, and TC), which were present in tumor tissue at slightly higher levels than in non-tumor tissue, had more than 44% of the PID pathways specific to these cell types. The high percentage of PID pathways that were important in the immune-related cell types is likely due to the relatively greater number of cytokine and other immune-associated pathways are included in the PID database.
All NSC clusters, except for NSC6 (3.6% pathways), had more than 10% of PID pathways that

Transcriptomic alterations in tumors compared to non-tumor at the single cell level
We next aimed to determine transcriptomic alterations in pediatric CNS tumors compared to non-tumor pediatric brain tissue. In bulk differential gene expression analyses, it is typically not possible to account for the impact of cell composition differences on gene expression levels [32][33][34][35] . Here, using single nuclei level data, we compared expression of the 4,000 most variable genes in nuclei from each tumor type to the gene expression of nuclei in non-tumor tissue, controlling for cell-type composition differences ( Figure 6A). Genes were considered differentially expressed if they met the FDR < 0.05 threshold.
As expected, adjusting for cell type proportions reduced the number of significantly differentially expressed genes compared with cell-type-unadjusted analyses. However, importantly, cell-type- Of the 4,000 most variable genes that were used in differential gene expression analysis, there were 558 genes that were differentially expressed in all six of the tumor types, 717 in five of the tumor types, and 596 in four of the tumor types compared with non-tumor tissue ( Figure 7A, Table 3). There were differentially expressed genes specific to a single tumor type: 43 genes for astrocytomas, 61 for embryonal tumors, 52 for ependymomas, 68 for glioneuronal/neuronal tumors, 98 for glioblastoma, and 57 for Schwannoma. While 60.9% (340/558) of the differentially expressed genes shared among all the tumor types were either increased or decreased the same direction, the remainder of genes varied in the direction of change based on tumor type compared to non-tumor tissue. The proportion of genes that either increased or decreased in the same direction for the shared significantly differentially expressed among all tumor types were significantly higher than expected (P-value < 2.2x10 -16 ). Proteincoding genes with increased expression across all tumor types included E2F7, ETS1, EZH2, ID3/4, MKI67, PIK3R3, and TOP2A. We conducted a pathways analysis of the genes with increased expression across all tumor types, and genes with decreased expression across all tumor types with Reactome pathways 37 . Interestingly, translation or nonsense mediated decay related processes having increased expression across all tumor types compared to non-tumor tissues (Figure 7B). Shared decreased protein-coding genes across all tumor types included FOXP2, GABRA1/2/4/5, NRGN, SST, and SYNPR. Even when differential gene expression analyses were adjusted for cell type, across all tumor types, there was decreased expression in genes associated with neuronal system such as transmission across chemical synapses and activation of NMDA or GABA receptors ( Figure 7C).
Hierarchical clustering of the differentially expressed genes revealed that transcriptomic alterations were similar in ependymomas and glioneuronal/neuronal tumors and likewise in astrocytomas and embryonal tumors (Figure 7A).

Discussion
In this study, we characterized gene expression profiles of 84,700 nuclei from snRNA-seq of 35 pediatric CNS tumors and 3 pediatric non-tumor brain tissues. We utilized an integrated hashtag oligonucleotide and genotype-based methods to maximize the number of sample-assigned nuclei from our multiplexed snRNA-seq experiment. Although the original MULTI-seq 20 work showed that multiplexing nuclei was feasible, some difficultly encountered with the approach in our study may have been attributable to use of fresh frozen samples that had been stored in the freezer for decades. In our study, we detail a novel approach to increase the number of cells assigned to a specific sample from pooled sequencing runs by integrating a genotype-based approach to demultiplex snRNA-seq data.
Future studies are expected to benefit from our integrated demultiplexing method to maximize data usage while decreasing the cost of snRNA-seq experiments.
Our study incorporates pediatric CNS tumor types that have not yet been characterized with single cell or single nuclei RNA-seq such as gangliogliomas. Moreover, we incorporated non-tumor pediatric tissues in our experiment, which to our knowledge have not been included in previous pediatric CNS tumor single cell RNA-seq studies. We describe changes in cell type proportions specific to each tumor type and use this information to identify the gene expression profiles and pathways enriched across tumor and normal samples through a cell type-adjusted analysis. We identified the pathways enriched in varying cell types, with a focus on neural stem like cells.
Since NSCs have been shown to be associated with therapy resistance, metastasis, and tumor malignancy, it is important to specifically consider NSCs when treating pediatric CNS tumors and reducing risk for secondary neoplasms [41][42][43][44][45][46][47] . We determined potential targetable NSC-specific pathways. While some commonly enriched pathways like MYC and FOXM1 in NSCs may be considered very difficult to target as MYC and transcription factors are considered to be less druggable, there were more easily targetable pathways enriched in NSCs like Aurora-B kinase and retinoic acid pathway.
With our cell type-adjusted approach, we addressed a critical confounder in differential gene expression analyses to identify transcriptomic alterations that exist in tumors compared to non-tumor tissue. Although the number of significantly differentially expressed genes decreased in the cell typeadjusted model compared to the cell type-unadjusted model, the adjusted model identified novel genes associated with tumors that would not have been uncovered in the unadjusted model. Moreover, the significantly differentially expressed genes exclusive to the unadjusted model likely stem from variations in cell type proportions, rather than from the underlying tumor biology that would be necessary for discovering effective therapeutic targets.
The pathways associated with the differentially expressed genes across the multiple tumor types in the cell type-adjusted model (translation associated processes like peptide chain elongation and translation initiation/termination along with nonsense mediated decay (NMD) processes) suggest the importance of these pathways commonly being dysregulated in pediatric central nervous system tumors. Previous studies have suggested the importance of downregulation of NMD responses in the differentiation of neural stem cells [48][49][50] . Moreover, high levels of NMD factors were sufficient to keep the stemness of neural stem cells 48 . Interestingly, our results indicate upregulation of NMD associated genes across all pediatric CNS tumor types in comparison to non-tumor pediatric brain which suggest the potential mechanism of upregulation of NMD maintaining more stem-like cells in these tumors. As more stem-like cells contributes to therapy resistance and recurrence, further studies investigating the NMD pathways and how they can be exploited to be potential therapeutic targets in pediatric CNS tumors are necessary.
Our study characterizes the heterogeneity that exists across pediatric CNS tumor types in comparison to non-tumoral pediatric brain tissue at the single cell level. We also identify potential tumor type and cell type-specific molecular characteristics that may be used therapeutic targets for the various pediatric CNS tumors from primary tissue samples. Although there were very limited samples for Schwannomas and glioblastoma, our study included thousands of nuclei from these tumor types to gain a better understanding of cells that exist in these tumor types that previous studies have not investigated yet. From our results, complementary preclinical in vitro and in vivo experiments are needed to validate these targets to advance these potential targets as therapeutic options in the clinic.

Study population
This study of pediatric central nervous system tumors was approved by the Institutional Review Board Study #00030211. Tumor and non-tumor tissues were collected from patients treated at Dartmouth Hitchcock Medical Center from 1993 to 2017. Patients consented to use of tissues for research purposes. Histopathologic tumor type and grade for each sample were re-reviewed according to the 2021 WHO classification of CNS tumors and categorized into the major tumor types 5 . Tumor types included in this study are astrocytoma, embryonal tumors, ependymoma, glioneuronal/neuronal tumors, glioblastoma, and Schwannoma. The average age at diagnosis of subjects from whom the tumor tissues were derived from in this study was 9.3 (range: 0. 75 -18). Male subjects accounted for 62.9% of the tumor samples and female subjects accounted for 37.1% of the tumor samples. Nontumor brain tissues were obtained from pediatric patients with epilepsy who underwent surgical resection. The average age at diagnosis of subjects from whom the non-tumor samples were derived from was 6.2 (0.58 -11). Male subjects accounted for 33.3% of the non-tumor samples and female subjects accounted for 66.7% of the non-tumor samples. Specific demographic characteristics of patients for the study are provided in Table 1 and sample information for each subject are provided in Supplementary Table 1.

Identification of genetic variation with bulk tissue RNA-seq
RNA was collected using Qiagen RNeasy plus kit (Catalog ID: 74034, Qiagen, Hilden, Germany). RNA-seq libraries were prepared following the Takara Pico v3 low input protocol and sequenced on Illumina NextSeq500.
Raw RNA-seq data were trimmed for polyA sequences and low-quality bases using cutadapt (v2.4) 51 . Reads were aligned to human genome hg38 using STAR (v 2.7.2b) 52 . Duplicate read identification and other quality control checks for read alignment were performed using CollectRNASeqMetrics and MarkDuplicates in Picard Tools. 53 Reads containing N were split using SplitNCigarReads function in the Genome Analysis Toolkit (GATK) 54,55 . Bases quality scores were recalibrated using known variants from the GATK resource bundle and with the BaseRecalibrator and ApplyBQSR functions in GATK 54,55 . Somatic SNV and indels were called with Mutect2 in tumor-only mode 54,55 . Only variants with at least read depth of 10, 5% allele frequency, read depth of 5 for the alternate allele were kept for analysis. The variants were then filtered for variants in sex or mitochondrial chromosomes, RNA editing sites, repeat masker regions, and variants in Panel of Normal (from GATK) references. Variants were then annotated using the Funcotator function in GATK 54,55 .

Identification of copy number variation with DNA methylation arrays
DNA were treated with sodium bisulfite following the TrueMethyl® oxBS Module (Tecan Genomics Inc, Redwood City, CA). Converted DNA were hybridized to Infinium HumanMethylationEPIC BeadChips. Raw idat files from the EPIC arrays were processed using preprocessNoob function in minfi in R 56  Lawrence, MA). We aimed for 2,500 -5,000 nuclei per sample to be sequenced.
Libraries for single nuclei RNA-seq were prepared following the 10x Genomics Single Cell Gene Expression workflows (10x Genomics, Pleasanton, CA) and were sequenced on Illumina NextSeq500 to average 45,000 reads per cell. 10X Cell Ranger software was used to align sequences to GRch38 pre-mRNA reference genome and generate feature-barcode matrices for downstream analyses.

Pre-processing snRNA-seq data
To filter low quality nuclei, only those with greater than 200 and less than 10,000 features and less than 5% of reads that map to the mitochondrial genes were used in downstream analyses. Pooled nuclei were demultiplexed by hashtag oligonucleotides using HTODemux function in Seurat v4 [58][59][60][61] .
Pooled samples were also demultiplexed using Vireo, a genotype based demultiplexing method 62 . We performed genetic demultiplexing analysis using genotype data following the methods described in Weber et al. 63 , implemented in a Nextflow workflow 64 . Briefly, bulk RNA-seq reads from each sample were mapped to the reference genome (GRCh38.p13) using STAR 52 . Pooled single-nuclei RNA-seq reads were mapped to the reference genome using STARsolo 65 . Variants among the samples within each pool were identified and genotyped with bcftools mpileup 66 using the mapped bulk reads.
Individual cells were then genotyped only at the sites identified using the bulk RNA using cellsnplite (mode 1a) 67 . Cell genotypes were used to identify the sample of origin for each cell using Vireo 62 . Code for the genetic demultiplexing workflow can be found at https://github.com/AlexsLemonade/alsf-scpca/tree/main/workflows/genetic-demux.
To integrate the methods, we first used sample identity assigned from the hashtag oligonucleotides. If the nuclei were confidently assigned a sample, it was compared to the genotypebased sample assignment. Those that did not match the same sample were filtered out. If the nuclei were assigned as a doublet or to none of the samples, the nuclei were assigned to a sample based on the genotype-based approach. 84,700 nuclei with confident sample assignment were used in analysis.
As our dataset included a very large number of nuclei to be integrated and was expected to have certain cell types only present in certain samples, we used the reciprocal PCA integration approach on the 2,000 most variable features to combine the nuclei from each sample. We first found the integration anchors with the FindIntegrationAnchors function then used the IntegrateData function in Seurat v4 to integrate all our filtered nuclei [59][60][61] .

Dimension reduction and clustering of snRNA-seq data
The integrated dataset was scaled using the ScaleData function in Seurat. First, PCA dimensionality reduction analyses were done to identify 100 principal components (PCs). To further reduce the dimensionality and cluster our nuclei by their gene expression profile, we conducted UMAP analyses on the 50 PCs with highest standard deviation with RunUMAP function in Seurat 58,68 . Then, we clustered our cells using FindNeighbors (n_neighbors = 30) and FindClusters (resolution = 1.0) function in Seurat 58 .

Gene set enrichment testing
Gene set enrichment tests at the single cell level were conducted using the Variance-Adjusted Mahalanobis (VAM) method 21 . The vamForSeurat function from the VAM R package was used to calculate enrichment scores for each nucleus. Brain cell type specific gene sets from the Molecular Signatures Database (MSigDB) v7.5.1 were used to validate our single cell identities [23][24][25][26][27][28] . For identifying cell types, p-values were calculated from the cumulative distribution function values generated by VAM. Nuclei were considered to be associated with a specific brain cell type-pathway if the VAM-generated p-value was ≤ 0.05. Nucleus-level pathway scoring was also conducted using VAM for pathways in the MsigDB Pathways Interaction Database (PID) collection 31 . PID Pathways were considered to be enriched in each nucleus at the FDR adjusted p-value threshold of 0.1 for the VAMgenerated p-values.
Stemness scores for each nucleus were calculated using the stemness-associated gene list from Tirosh et al 29 and the AddModuleScores function in Seurat.

Differential gene expression and pathways
Differential expression analysis between tumor nuclei and non-tumor nuclei were conducted using monocle3 [69][70][71][72] . Differential expression analyses were conducted only on the top 4,000 most variable features identified from the FindVariableFeatures Seurat function. The unadjusted differential expression testing was done using the fit_models function in monocle3 (v1.0.0) R package with the quasi-poisson distribution with the non-tumor nuclei being the referent gene expression profile 69,71,72 .
The adjusted differential expression testing was done with the same quasi-poisson distribution with non-tumor nuclei being the referent but including the major cell type identity in the model. Gene types for each gene used in the differential expression testing were annotated using the org.Hs.eg.db 73 , Human genome annotation package, and mapIds function in the AnnotationDbi R package 74 . Pathways associated with the differentially expressed genes were identified using the Reactome pathways and ReactomePA R package 37 .
Pathways important for each cell cluster were identified using FindAllMarkers function in Seurat v4 with the Wilcoxon rank sum test in Seurat on the binary classification of PID pathways enrichment for each nuclei 75 . Log fold change and minimum percentage of cells enriched in each pathway were both set to 0. To identify the pathways with greater number of nuclei with enriched pathway per cluster, we selected pathways that were only positive in direction in the FindAllMarkers options.

Statistical testing
Observed proportion of genes that were either increased or decreased in the same direction in the shared differentially expressed genes among all the tumor types (60.9%) were compared to expected proportion of genes that would be increased or decreased in the same direction across all the tumor types (3.13%) using a one-sample proportion test. The expected proportions were determined based on the permutations of direction of change compared to non-tumor for the six tumor types.

Availability of data and materials
The single nuclei RNA-seq datasets used in the current study is available in GSE211362. A token for access by reviewers is available by contacting the corresponding author. Detailed annotation of which pool each sample was multiplexed in can be found in Supplementary Table 1.

Funding
This work was supported by a Prouty Pilot award from the Dartmouth Cancer Center and a Single-cell Pediatric Cancer Atlas (ScPCA) grant from the Alex's Lemonade Stand Foundation. MKL was supported by the Burroughs-Wellcome Fund: Big Data in the Life Sciences at Dartmouth. NA was supported by the S.M. Tenney Fellowship at Dartmouth. This work was also supported by R01CA216265, R01CA253976, and P20GM104416 -6369 to BCC and P20GM130454 -7702 and R21CA253408 to HRF. Single nuclei RNA-seq experiments were conducted in the Genomics and Molecular Biology Shared Resource (GMBSR) at Dartmouth, which is supported by NCI Cancer Center Support Grant 5P30CA023108 and NIH S10 (1S10OD030242) awards. Single nuclei RNA-experiments were also supported through the Dartmouth Center for Quantitative in collaboration with the GMBSR with support from NIGMS (P20GM130454) and NIH S10 (S10OD025235) awards.