https://orcid.org/0000-0003-1058-4604
Abstract
Pediatric high-grade gliomas (PHGG) are aggressive brain tumors with 5-year survival rates ranging from <2% to 20% depending upon subtype. PHGG presents differently from patient to patient and is intratumorally heterogeneous, posing challenges in designing therapies. We hypothesized that heterogeneity occurs because PHGG comprises multiple distinct tumor and immune cell types in varying proportions, each of which may influence tumor characteristics.
We obtained 19 PHGG samples from our institution’s pediatric brain tumor bank. We constructed a comprehensive transcriptomic dataset at the single-cell level using single-cell RNA-Seq (scRNA-Seq), identified known glial and immune cell types, and performed differential gene expression and gene set enrichment analysis. We conducted multi-channel immunofluorescence (IF) staining to confirm the transcriptomic results.
Our PHGG samples included 3 principal predicted tumor cell types: astrocytes, oligodendrocyte progenitors (OPCs), and mesenchymal-like cells (Mes). These cell types differed in their gene expression profiles, pathway enrichment, and mesenchymal character. We identified a macrophage population enriched in mesenchymal and inflammatory gene expression as a possible source of mesenchymal tumor characteristics. We found evidence of T-cell exhaustion and suppression.
PHGG comprises multiple distinct proliferating tumor cell types. Microglia-derived macrophages may drive mesenchymal gene expression in PHGG. The predicted Mes tumor cell population likely derives from OPCs. The variable tumor cell populations rely on different oncogenic pathways and are thus likely to vary in their responses to therapy.
Analysis of scRNA-Seq data from PHGG patient samples enables identification of tumor and immune cell types.
Tumor cell types include astrocytes, oligodendrocyte progenitor cells (OPC), and OPC-derived mesenchymal-like cells.
Immune cells include a macrophage population that may induce the transformation of OPCs into mesenchymal-like cells.
PHGG is diverse intratumorally and from patient to patient, posing challenges for designing therapies and leading to therapy resistance. Using single-cell transcriptomics, we identified multiple tumor cell types in a diverse group of 19 PHGG patient samples. We also found that CNS-resident macrophages likely induce the emergence of mesenchymal-like tumor cells. Differential gene expression and gene set enrichment analyses revealed distinct oncogenic processes driving proliferation in the different tumor cell types. These distinct processes were not resolved with bulk RNA-Seq analysis. This work paves the way for investigation and development of treatment strategies targeting distinct oncogenic processes in differing tumor cell types of PHGG.
Pediatric high-grade glioma (PHGG) features extensive inter- and intratumoral heterogeneity.1,2 This presents challenges in characterizing and understanding basic tumor biology, as well as in treating these deadly tumors. PHGG are highly invasive and grow diffusely among normal cells, limiting surgery as a therapeutic option.3,4 Radiation therapy is transiently effective, but the tumors nearly always recur.4 Despite hundreds of clinical trials of cytotoxic and targeted chemotherapy, radiation remains the only accepted standard-of-care therapy.
The specific types and origins of PHGG cells have not been fully characterized. Recent studies suggest oligodendrocyte progenitors as one PHGG cell of origin.5,6 A study of primarily adult HGG that included a few pediatric samples proposed differing cell states and whether cells had undergone a mesenchymal transition as agents of tumor heterogeneity.7 This study, however, did not separately analyze the pediatric tumors and the resulting model is attuned to the majority adult HGG sample population rather than the limited number of pediatric samples. Differences in tumor biology limit the application of adult studies to PHGG. PHGG and adult HGG arise from different embryonic layers.8 Adult HGGs frequently include EGFR amplification, TERT promoter mutations, IDH alterations, and chromosome 7 gain/chromosome 10 losses that rarely occur in PHGG.7,8 These differences require separate investigation of PHGG biology, using exclusively pediatric samples.
Here, we identify normal CNS cell type analogs to PHGG cells, the CNS-resident and systemic immune cells present in PHGG, and the potential interactions between the tumor and immune cells. We hypothesized that PHGG comprises several distinct glial lineage cell types in varying proportions from case to case and that the tumor cell phenotypes may be modified through the activity of systemic or CNS-resident (microglial) immune cells. To test our hypotheses, we performed single-cell RNA-Seq (scRNA-Seq), bulk methylation, and bulk RNA-Seq on 19 patient-derived PHGG samples. Our cell capture included immune cells, which to date have received little attention in single-cell studies of PHGG.5,7 We separately analyzed the single-cell transcriptomic data from tumor and immune cells using unsupervised clustering to identify gene expression patterns and supervised clustering to map cells onto analogous normal cells. Using differential gene expression and pathway analysis, we elucidated roles and identified specialized populations related to tumor growth. We conducted multi-channel immunofluorescence (IF) staining to validate the results of transcriptomic analyses at the protein level and to investigate spatial relationships among tumor and immune cells.
Materials and Methods
Methods are summarized below and more fully described in the Supplementary Methods.
Sample Collection
Patients were enrolled on an institutional review board–approved protocol for the use of additional tumors and other tissue for research (Colorado Multi-Institutional Review Board protocol 95-500) after they or their families provided written consent. Sample inclusion in the study was based on clinical diagnosis of a Grade 3 or 4 glial tumor using the classification system in use at the time of diagnosis (Supplementary Table S1). We screened the samples for IDH1R132 variants (of which there were none), and histone 3 variants H3K27M and H3G34R/V.
Single-Cell RNA-Seq Sample Preparation and Sequencing
Frozen disaggregated tumor cells were thawed and flow-sorted to yield viable single cells. The sequencing library was prepared using a Chromium Controller in combination with Chromium Single Cell V3 Chemistry Library Kits, Gel Bead & Multiplex Kit, and Chip Kit (10X Genomics). Samples were sequenced on a HiSeq 2500 or NovaSeq6000 (Illumina) to obtain approximately 50 thousand reads per cell.
Bulk RNA-Seq Sample Preparation and Sequencing
DNA and RNA were extracted from frozen tumor samples using the AllPrep DNA/RNA Mini Kit (Qiagen). RNA library preparation was performed using poly A selected total RNA. Sequencing was conducted using paired end reads with a sequencing length of 50 bp and 40 million reads per sample.
Bulk Methylation Profiling Sample Preparation and Sequencing
DNA extracted from patient tumor tissue samples underwent bisulphite conversion. Samples were scanned using the 850K Infinium Methylation Epic BeadChip Kit (Illumina).
scRNA-Seq Data Analysis
Raw sequencing reads were demultiplexed, mapped to the human reference genome (build GRCh38) and gene-expression matrices generated using CellRanger (version 3.0.1). Count matrices were filtered in Seurat 3.1.0.9,10 We applied Harmony alignment to correct for inter-sample variation due to experimental or sequencing batch effects.11 We analyzed the data using the Seurat, Azimuth, and scCustomize packages.12
Large-Scale CNV Inference of scRNA-Seq Data
Chromosomal CNVs of single cells were inferred on the basis of average relative expression in variable genomic windows using InferCNV (https://github.com/broadinstitute/inferCNV). We deemed cells with large-scale CNVs to be tumor cells and cells lacking them to be immune cells. Results were validated by clustering predicted tumor and immune cells separately to determine the overlap of clusters.
Cell Mapping to Normal Cell Populations
We used the Azimuth package (0.4.3) (https://github.com/satijalab/azimuth) running under R (4.1.2). Tumor cells were mapped to a human motor cortex reference and immune cells to a PBMC reference.13,14 Validation included a comparison of mapping to unsupervised clustering followed by manual annotation and unsupervised clustering of mapped clusters to determine homogeneity and separation from other cell types.
DNA Methylation Data Analysis
DNA methylation results were uploaded to the DKFZ classifier to determine tumor classification (Supplementary Table S1, https://www.molecularneuropathology.org/mnp).15
Bulk RNA-Seq Data Analysis
Raw sequencing reads were demultiplexed and mapped to the human reference genome (build GRCh38). Expression was derived for approximately 66 000 transcripts and differential expression was analyzed using ANOVA.16–18 Following filtering, approximately 30 000 transcripts remained. Pathway analysis was performed using geneset enrichment analysis (GSEA),19,20 and single-sample GSEA (ssGSEA).21
Immunofluorescence Staining
Formalin-fixed paraffin embedded slides were obtained for 8 patient samples on which scRNA-Seq was performed. Primary antibody incubation was performed overnight at 4°C followed for unconjugated primaries by incubation with an Alexa Fluor (Invitrogen) secondary antibody for 1 h at room temperature. Imaging was performed at 200X using a Zeiss 780 LSM confocal microscope. Images were obtained using Z-stacking and quantified in ImageJ using stack depth sufficient to obtain an acceptable signal to noise ratio without blurring the image. Image quantification was performed in ImageJ and analyzed in R (4.2.2).
Multiplexed Immunofluorescence Staining
Seven-channel whole slide and nine-channel region of interest staining and imaging were performed to identify immune cell types (Supplementary Data S9–10). Slides were imaged at 200x on a Vectra Polaris system. We used Phenochart (version 1.1.0, Akoya Biosciences) to perform initial quality control and spectral unmixing, followed by Inform (version 2.6, Akoya Biosciences) to perform phenotyping and scoring. We analyzed, organized, and displayed the data in R (version 4.1.2) using phenoptr (version 0.3.2) and phenoptrReports (0.3.3) (Akoya Biosciences) along with routines from tidyverse, including ggplot2, and pheatmap.
Results
Patient and Disease Characteristics
We classified the samples into 6 subgroups: H3K27M diffuse midline glioma (5 samples), PHGG with epithelioid/giant cell characteristics (e-GBM, 4 samples), radiation-induced high-grade glioma (RIG, 5 samples), IDH wild-type hemispheric glioblastoma (GBM, 3 samples), diffuse hemispheric glioma H3 G34-mutant (G34-GBM, 1 sample), and PHGG otherwise unspecified (1 sample; Supplementary Table S1). Inclusion in the e-GBM category was based on clinicopathological identification of epithelioid/giant cell features. Identification of the RIG cases was based upon clinical application of the Cahan criteria for definition of radiation-induced tumors.22,23 We hypothesized that variation in cell type composition is a general characteristic of PHGG that explains its phenotypic variations and thus conducted our analysis and report most of our results on the basis of the PHGG population as a whole. Where meaningful differences by subtype were identified, they are reported in the results.
Tumor Cell Populations Comprise Distinct Cell Types
Unsupervised clustering of the complete scRNA-Seq dataset distinguished tumors from immune cells and identified multiple subtypes of each (Figure 1A and B, Supplementary Figure S1A–C). We used large-scale copy number variations (CNV) in the scRNA-Seq data to segregate tumor and immune cells for analysis (Supplementary Figure 2A). Due in part to the diversity among tumor samples, we found that unsupervised clustering of tumor cells produced multiple patient-specific subclusters of the same cell type. To obtain cell groupings that were easier to analyze, we decided to map by cell type to a normal brain atlas with supervised clustering and then perform unsupervised clustering to analyze each cell type (Figure 1C).13 We similarly mapped the tumor immune cell population to an atlas of normal PBMCs.14 The supervised clustering verified that tumor and immune cell types clustered separately (Supplementary Figure S2B and C). The total cell population of 16 452 cells segregated into 10 160 tumor and 6292 immune cells. Following the removal of 2812 low-expressing tumor cells that we attributed to poor cell condition or difficulties in library preparation, the populations for analysis comprised 7348 tumors and 6292 immune cells (Figure 1D).
CNV analysis and bioinformatic clustering identify normal CNS analogs of tumor cells. a. and b. UMAP projections of tumor cells show Seurat unsupervised cluster membership (a) is consistent with cell type (b), but Seurat clusters include additional subgroups due to phenotypic differences such as expression of cell cycle, angiogenic, metastasic, or inflammatory genes; c. Analysis strategy for scRNA-Seq data; d. Predicted nonimmune cell types from supervised clustering; e. Predicted tumor cell type percentages by sample (one sample is omitted from this display due to cell count <10 following elimination of low-expressing cells).
The tumor cells were mapped to 5 principal predicted cell types: oligodendrocyte progenitors (OPCs), mesenchymal-like cells (Mes), astrocytes, neurons, and oligodendrocytes (Figure 1D). Individual samples varied in their proportions of the primary cell types (Figure 1E). Analysis by PHGG subtype showed the RIGs are enriched in Mes cells and epitheloid tumors are depleted in OPCs. (Supplementary Figure S3A and B). We used unsupervised clustering to verify the primary cell types identified using supervised clustering (Supplementary Figure S4A). Primary cell types were distributed among multiple clusters in the unsupervised dataset because cellular processes such as cell division, invasion, or angiogenesis that occur in a portion of the cells of a particular type produce sufficiently differential gene expression to cluster separately (Figure 1B and C, Supplementary Figure S4A). Accordingly, supervised clustering, which considers only cell type characteristics in assigning cells to clusters, was preferred for most of the analyses performed.
Astrocytes, OPCs, and Mes cells comprised the primary tumorigenic cell populations. To validate the assignment of tumor cells to these 3 primary cell types, we performed unsupervised clustering by cell type (Supplementary Figure S4B). Astrocyte and OPC subclusters separated, validating their cell type assignment. Mes cell subclusters, which as discussed below are likely derived from the OPCs, clustered primarily with outlying OPC subclusters (Supplementary Figure S4B). We omitted the predicted neurons and oligodendrocytes from the tumor cell analyses because they did not demonstrate tumor cell characteristics in their gene expression and may be nontumor cells that were entangled in the resected samples (Supplementary Figure S5A).
Differential gene expression.—
We analyzed differential gene expression of the unsupervised clusters and predicted tumor cell types to better understand their character and potential origins. We also conducted gene set enrichment analysis (GSEA) and single-sample gene set enrichment analysis (ssGSEA) to identify enriched and depleted expression pathways.10,14,19,24 K means clustering and non-negative matrix factorization of the predicted cell type clusters showed the cell types form distinct entities based on gene expression (Figure 2A and B). The OPCs and Mes clusters were more closely related than either was to the astrocyte population (Figure 2B).
Differential gene expression among predicted tumor cell types. a. Heatmap by gene expression shows segregation of the 3 principal tumor cell types; b. Non-negative matrix factorization (NMF) of principal cell types shows OPC and Mes cells are similar but differ markedly from astrocytes; c. Marker gene expression confirms cell types differentially express known characteristic markers; d. OPCs are enriched in ASCL1, SOX2, and SOX10, which control stemness characteristics, and the stem cell marker NES; e. OPCs are enriched in transcription factors whose expression precedes emergence of a mesenchymal phenotype and Mes cells are enriched in genes characteristic of mesenchymal cells, suggesting Mes cells may develop from OPCs during tumor growth; f. OPC and Mes marker genes in predicted cell types; g. Ki67 staining shows OPCs and astrocytes are more proliferative than the average of all cells; Mes cells have a proliferative subpopulation but on average are less proliferative than the average cell population, representative IF images in Supplementary Figure S3a–b.
The OPCs highly expressed the oligodendrocyte lineage transcription factors OLIG1 and OLIG2 (Figure 2C, Supplementary Figure S5A, Supplementary Data S1). The OPCs were differentially enriched in transcription factors ASCL1, SOX2, and SOX10, which are important in stem cell maintenance and CNS development, and the stemness marker NES (Figure 2D). The Mes cells were also enriched in PDGFRA, as well as IGFBP2 and IGFBP5, which are associated with a mesenchymal phenotype as well as increased migration and invasion in adult HGG (Figure 2C).25–29 Most of the most highly differentially enriched genes in the Mes cell type, including GNAS, IGFBP2, IGFBP5, PTN, S100B, CD63, and VIM play a role in metastatic processes such as loss of adhesion, cell motility, and invasiveness (Figure 2C, Supplementary Figure S5A, Supplementary Data S1).30–33
We noticed that genes such as OLIG1, OLIG2, and PTPRZ1, which were highly expressed in the OPCs, were also expressed in the Mes population, and that highly expressed Mes genes such as IGFBP2, GNAS, and PDGFRA, were also expressed in the OPCs (Supplementary Figure S5B). The OPCs also expressed PTN and S100B, which are involved in inducing the emergence of the mesenchymal phenotype (Figure 2C). We wondered whether the Mes cells could be derived during oncogenesis or progression from OPCs. We plotted expression of transcription factors that incept expression of mesenchymal genes, cadherins that remodel cells undergoing transformation to a mesenchymal phenotype, and genes expressed in mesenchymal cells (Figure 2E). The plots show a progression from mesenchymal-inducing transcription factor and cadherin expression in OPCs but not Mes cells, to mesenchymal gene expression in Mes cells but not OPCs. ZEB1 and ZEB2 gene expression in OPCs and VIM expression in the Mes cells was particularly noteworthy because ZEB transcription factors induce VIM expression and oncogenic progression.34 We plotted the expression of OPC markers and genes involved in mesenchymal cell differentiation, which also distinguished the OPCs from the Mes cells (Figure 2F). These results suggest Mes cells may be derived from OPCs. We investigate a possible mechanism below.
The predicted astrocyte population was highly enriched in expression of the characteristic astrocytic genes GFAP (log2FC = 19.0, P = 0) and APOE (log2FC = 25.5, P = 0), as well as metallothionein family genes (Figure 2C, Supplementary Figure S5A, Supplementary Data S1). Metallothioneins promote cell proliferation, inhibit apoptosis, and promote angiogenesis.35 High metallothionein expression correlates with poor survival in adult HGG.36 In contrast to adult HGG, the astrocytes were only slightly enriched in EGFR (FC = 0.58, adj. p = 8.31E−12). The measured gene expression of the predicted astrocytes confirmed their identification as an astrocytic cell population.
To determine whether astrocyte, OPC and Mes populations contributed to tumor expansion, we assessed whether they are enriched in expression of the proliferation marker MKI67. While all three cell types included a fraction of cells that express MKI67 (Supplementary Data S1), it was unclear whether they all contributed to tumor expansion. Accordingly, we performed immunofluorescence (IF) co-staining for Ki67 with the astrocyte, OPC, and Mes cell markers GFAP, Olig2, and IGFBP2, respectively. These results confirmed that all 3 cell types are proliferative and contribute to tumor expansion (Figure 2G, Supplementary Figure S6A, B).
Differential gene expression analysis established the separate identities of the astrocyte, OPC, and Mes populations, which comprise 3 principal tumor cell types in our PHGG patient samples. Staining for Ki67 confirmed that all 3 cell types contribute to tumor expansion. Gene expression patterns confirmed that OPCs and Mes cells are closely related and suggested the possibility, further investigated below, that Mes cells are derived from OPCs.
Pathway analysis.—
To further probe distinctions among the predicted tumor cell types, we conducted GSEA to identify enriched signaling and functional pathways. We started with gene sets used to describe expression subtypes of adult HGG.37 GSEA using bulk RNA-Seq data showed samples varied in their overall enrichment of genes that characterize the mesenchymal, classical, and proneural adult HGG subtypes, confirming heterogeneity within our sample population (Figure 3A, Supplementary Data S3). We wondered whether differing gene expression among specific cell types could explain the observed heterogeneity.
Tumor cell types differ in pathway enrichment. a. Enrichment in gene sets from adult HGG shows intertumoral heterogeneity based on bulk RNA-seq data (FDR < 1.6E−4 for all gene sets); b. Tumor cell types and microglia vary in enrichment of adult HGG gene sets; c. Pathway analysis in tumor cell types; d. Histone methylation mark gene sets predicting stemness/differentiation (greater enrichment in these gene sets suggests greater cell differentiation); e. Heatmap summarizing pathway enrichment shows differences in oncogenic processes and targetable signaling pathways among principal tumor cell types.
GSEA of the scRNA-Seq data revealed that tumor cell types varied in their expression of adult HGG gene sets: OPCs were primarily proneural, Mes cells displayed a mix of classical and proneural expression, and astrocytes were primarily classical (Figure 3B). None of the 3 tumor cell types was enriched in the expression of the adult HGG mesenchymal geneset, but we found an immune cell population with this enrichment (Figure 3B). These immune cells were enriched in the microglial marker, AIF1, so we tentatively assigned them an identity of CNS-resident immune cells. We further studied these cells and describe our findings below. We note at the outset that the adult HGG mesenchymal geneset refers to a specific tumor characteristic that includes inflammatory gene expression but does not purport to measure whether cells have undergone a mesenchymal transition.37 (Supplementary Data S11). For example, key mesenchymal transcription factors and genes, including VIM, are not included in the adult GBM mesenchymal geneset (Supplementary Data S11).
The astrocytes exhibited multiple differences in pathway expression versus the OPC and Mes cell types (Figure 3B–D). On average, the OPCs and Mes cells were enriched in self-renewal, cell cycle, protein translation, and DNA repair gene sets (Figure 3C, Supplementary Data S2). They were metabolically active and enriched in the expression of targets of the MYC transcription factor and, as noted above, genes associated with proneural adult HGG (Figure 3B, C, Supplementary Data S2). The OPCs and Mes cells were also enriched in core embryonic stem cell gene expression. (Figure 3D, Supplementary Data S2). The astrocytes, in contrast, were enriched in epigenetic histone marks associated with lineage commitment (Figure 3D).38,39 One astrocyte subpopulation, however, was enriched in genes expressed in basal radial glia—a precursor to both neuronal and glial lineage cells during normal development—suggesting that a portion of them were highly immature (Figure 3C, Supplementary Data S2).40 The OPCs were depleted and the astrocytes enriched in expression of inflammatory and immune response pathways (Figure 3E). As compared to the OPCs, the Mes cells were enriched in inflammatory and immune response pathways, but less so than the astrocytes (Figure 3E).
We cataloged by cell type the expression levels of gene sets representing common cancer signaling pathways (Figure 3E, Supplementary Data S3). Notably, the astrocytes and Mes cells were enriched in expression of inflammatory pathways, astrocytes were depleted in DNA repair gene expression, and variability existed in stem cell gene enrichment. The cell types differed substantially in expression of cancer pathway genes, but all were enriched in at least one of the cancer signaling pathways that we assessed (Figure 3E).
In summary, the astrocyte, OPC, and Mes cell populations had gene expression programs that confirmed their identity as tumor cells. The variability in observed gene expression suggests different oncogenic drivers for the 3 tumor cell types (Figure 3E). All 3 populations included undifferentiated cell populations, suggesting the presence of progenitors for each cell type.
Different Tumor Cell Types Rely on Distinct Oncogenic Pathways
To evaluate the heterogeneity and understand the gene expression enrichment and consequent activation of tumorigenic pathways within each predicted cell type, we analyzed differential gene expression and conducted ssGSEA comparing subclusters of each predicted cell type.19,21 Unsupervised clustering of the astrocytes yielded 9 subclusters. Six had enriched mesenchymal gene expression (P = .0033, Figure 4A, Supplementary Data S4). Among the mesenchymal astrocyte clusters, 5 of 6 were enriched in NF-κB pathway expression versus one of 3 nonmesenchymal clusters (p = .002, Supplementary Data S4). NF-κB pathway gene expression was highly correlated with VEGF pathway activation (r = 0.85, P = .0036) (Figure 4B, Supplementary Data S4), suggesting involvement of the mesenchymal astrocyte clusters in angiogenesis. The astrocytes also included subpopulations with expression suggesting stem or progenitor states (Figure 4C, Supplementary Figure S5C).
Pathway and differential gene expression analysis reveal differences by cell type in angiogenesis, metastasis, and stem cell gene expression. a. Expression of mesenchymal versus epithelial pathway genes for astrocyte subclusters; b. Astrocyte subcluster expression of VEGF pathway genes strongly correlates with NF-κB pathway expression (r = 0.85, P = .0036); c. OPC cluster 6 expresses apical radial glial cell genes (“Multi” box), suggesting neural stem cell identity, and OPC cluster 7 expresses committed progenitor genes suggesting a hierarchical relationship within the OPC population (left panel); astrocyte subpopulation 5 expresses apical radial glial cell genes (“Multi” box) and subpopulation 2 expresses committed astrocyte progenitor genes (“As-IPC” boxes) (right panel); d. Expression of mesenchymal versus epithelial genes for OPC subclusters; e. OPC subcluster has low expression of VEGF pathway genes that is uncorrelated with also low NF-κB pathway expression; f. Expression of mesenchymal versus epithelial genes for Mes subclusters; g. Mes subcluster expression of VEGF pathway genes strongly correlates with NF-κB pathway expression (r = 0.84, P = .0042).
Among the 8 predicted OPC subclusters, 3 were mesenchymal (Figure 4D, Supplementary Data S4). NF-κB pathway expression was not correlated with mesenchymal character or VEGF pathway expression (Figure 4E, Supplementary Data S4). Two OPC subclusters, OPC-6 and OPC-7, were enriched in expression of genes indicating stem or progenitor states (Figure 4C, Supplementary Figure S5C). OPC-6 was enriched in genes expressed by multipotent apical radial glial cells, a progenitor capable of differentiating into oligodendrocyte or astrocyte progenitors (Figure 4E, Supplementary Figure S5D).41 OPC-7 expressed committed progenitor genes suggesting a potential lineage hierarchy within the OPC subcluster (Figure 4E, Supplementary Figure S5D).
Four of the 5 Mes cell subclusters were mesenchymal (P = .0286, Figure 4F, Supplementary Data S4). Expression of NF-κB pathway genes was highly correlated with VEGF pathway gene expression (r = 0.84, P = .0042) (Figure 4G). The Mes cells included a subpopulation (Mes-2) that differentially expressed proliferative genes (Supplementary Figure S5E). The Mes cells did not include a stemlike cluster with parallel expression to that of OPC-6.
Analysis of heterogeneity within tumor cell types identified potentially important subpopulations. The OPCs included possible stemlike and transit amplifying subpopulations. The astrocyte and Mes cells included predominant subpopulations with a mesenchymal phenotype that were enriched in expression of angiogenic pathways. In the astrocytes, the NF-κB pathway was enriched in the mesenchymal subpopulations, suggesting a possible driver of the mesenchymal phenotype and angiogenesis.
Tumor-Resident Immune Cells Modify the Gene Expression of Tumor Cells
Tumor-resident immune cells were mapped to a peripheral blood mononuclear cell (PBMC) atlas using supervised clustering of the scRNA-Seq data (Supplementary Figure S7A, B). To validate the mapping of immune cells to the PBMC reference and identify interactions between CNS-derived immune cells and tumor cells, we analyzed immune cell gene expression pathways using the scRNA-Seq data. We also conducted immunofluorescence staining to validate key markers and identify spatial relationships among cell types.
Immune cell population.—
Tumor-resident immune cells included monocytes (including macrophages), microglia, and T cells, with smaller populations of dendritic and hematopoietic stem/progenitor cells (HSPCs) (Figure 5A, B). The number of immune cells and their relative proportions varied widely among tumor samples (Figure 5C). To validate the immune cell mapping, we analyzed differential gene expression and performed GSEA on each predicted cell type using gene sets derived from immune cells of known type (Supplementary Data S5–6). For gene sets matching predicted cell type, GSEA yielded consistently positive normalized enrichment scores (NES) with a statistically significant false discovery rate (Figure 5D). This validates that the supervised clustering accurately mapped the tumor-resident immune cells to the appropriate immune cell type. To better understand the composition of the microglia population, we co-stained for CD163, CD74, and IBA1, from which we estimated that on average 54.8% of the total microglia are type 2 macrophages (CD163+/CD74−/IBA1+, SEM = 0.042, Figure 5E, Supplementary Figure S6C).
Immune cell classification. a. Counts and proportions of cell types found in tumor-associated immune cells; b. Key marker gene expression by immune cell type; c. Immune cell totals and proportions by sample; d. Normalized enrichment score of gene sets indicating identity of immune cells, false discovery rate (FDR) shown by symbol color, error bars show standard error of the mean; e. Proportion by tumor sample of microglial population (CD74−/IBA1+) that expresses type 2 (tumor-associated) macrophage marker CD163 (mean with standard error of the mean shown), representative IF images in Supplementary Figure S3c.
We performed a correlation analysis to determine whether tumor cell composition varied with the immune cells present. Of 15 potential combinations of tumor and immune cells, we identified 3 instances where immune cells were depleted and 2 in which they were enriched versus the null hypothesis of no systematic variation with tumor cell composition (Figure 6A, Supplementary Data S7). The most notable change was that monocytes (including monocyte-derived macrophages) decreased, and microglia (including CNS-resident macrophages) increased with increasing Mes tumor cell population size.
Microglial derived macrophages may induce mesenchymal transformation and evidence for T-cell exhaustion and tumor permissive immune microenvironment. a. Correlations between tumor and immune cell types show microglia-associated macrophages are enriched and monocyte-associated macrophages depleted in Mes subtype; FDR (q < 0.05) for all data points; b. Heatmap of immune cell types shows microglial cells are enriched in mesenchymal adult HGG gene set; FDR < 0.05; c. Nearest neighbor analysis using immunofluorescence shows Mes (IGFBP2+) cells are on average closer to microglial derived type 2 macrophages (CD163+/CD74−) than to OPCs (Olig2+); data points show statistical mode by sample; total number of cells analyzed exceeds 290 for each sample, representative IF images in Supplementary Figure S3d; d. Microglia-derived macrophage (micro) population has greater proportion of M2 macrophages than CD14+ or CD16+ monocyte-associated macrophages; e. M1 (HLA-DPA1, C3) and M2 (CD163) markers validate macrophage classification by gene expression; f. Marker gene expression suggests T-cell exhaustion; exhaustion markers (CD44, SPN, CD69, PRDM1 are enriched), while enrichment markers (IL7R, SELL, GZMB) are depleted; g. GSEA showing negative NES for gene sets upregulated during tumor rejection suggests tumor permissive immunosuppressed microenvironment.
We wondered how some PHGG samples could have bulk enriched expression of the adult HGG mesenchymal gene set given the lack of enrichment of these genes in the principal tumor cell types (Figure 3A, B). To answer this question, we turned to the microglia, which were highly enriched in adult HGG mesenchymal gene set expression and depleted in the adult proneural gene set (Figure 6B, Supplementary Data S6). The microglia were also depleted in embryonic stem cell genes and were nonproliferative (Supplementary Data S5). The microglial population size correlated positively only with the Mes cells (Supplementary Data S7). Among the immune cells, the only strong enrichment in any adult HGG subgroup was the microglia/Mes tumor cell relationship (Figure 6B). A similarly strong depletion existed between the HSPCs and the mesenchymal adult HGG gene set (Figure 6B).
IF staining for the type 2 macrophage marker CD163, along with monocyte, Mes, and OPC markers showed that, on average, type 2 macrophages were situated next to Mes cells more often than they were to OPCs (Figure 6C, Supplementary Figure S6D, Supplementary Data 12–13). The Mes-macrophage distance, measured from nuclear center to nuclear center, is on the order of a cell diameter, suggesting interactions such as receptor-ligand or paracrine signaling could occur. Using the scRNA-Seq data, we compared expression of type 1 and type 2 macrophage polarity in the microglial and monocyte-derived macrophage populations (Figure 6D, Supplementary Data S8). We utilized gene sets representing types 1 and 2 polarity states from a breast cancer model, which we applied to assign polarity to the macrophage populations.42 The microglial-associated macrophage (micro) population had a greater proportion of type 2 macrophages than either of the monocyte-associated macrophage populations (Figure 6D, Supplementary Figure 7C, D). We used the markers HLA-DPA1, C3, and CD163 to verify identification of types 1 and 2 polarity (Figure 6E).
To obtain an overall picture of macrophage gene expression we iteratively reclustered the macrophage population, removing nonmacrophage cells after each iteration (Supplementary Figure S8A). We verified the final population was enriched in monocyte (CD74), microglial (AIF1 – IBA1 protein), and macrophage markers (Supplementary Figure S8B). Turning to specific microglial markers, we found variable results that suggested the possibility the tumor microenvironment may alter microglial and monocyte gene expression in a manner that ordinary tissue does not (Supplementary Figure S8C). The macrophage population included cells that express the CCL2 receptor CXCR2, suggesting immune suppression could be occurring through myeloid-derived suppressor cell infiltration.43 Macrophage subpopulations were also enriched in TREM2 and CD84, which are each associated with immune suppression in the tumor microenvironment.44,45
We investigated gene expression differences that might reflect T-cell exhaustion. CD4 and CD8 positive T cells, NK, and dendritic cells were all enriched in genes associated with an exhaustion phenotype and depleted in genes associated with a nonexhausted phenotype (Figure 6F). All 4 T-lineage cell types were also depleted in genes upregulated during tumor rejection, suggesting a tumor immune microenvironment consistent with exhaustion (Figure 6G).
To better understand the tumor immune microenvironment and validate the scRNA-Seq analyses of immune cells, we performed multi-channel immunofluorescence staining using a panel of immune cell markers. We stained for markers of T cells, monocytes, dendritic cells, HSPCs, macrophages, and microglia. We also stained for markers of T-cell activation, exhaustion, and suppression, and of types 1 and 2 macrophage polarization. Like the scRNA-Seq data, the multi-channel IF staining also showed that total immune cells and the cell type mix varied by sample and PHGG subtype (Supplementary Figure 9A, B, Supplementary Data 9–10). The IF data showed larger dendritic cell population sizes and smaller microglia populations than the scRNA-Seq data (Figure 5C, Supplementary Figure 9A, B). IF staining also showed generally larger type 2 macrophage fractions than the scRNA-Seq data (Figure 6D; Supplementary Figure 9C). The methods of identification were different, however, with the scRNA-Seq analysis relying on expression of genes shown to be upregulated in macrophages with a type 2 phenotype, while the IF staining classified CD68+ macrophages as type 1 and CD163+ as type 2. The IF staining enabled estimation of overall T-cell density in the stained samples as well as proportions of exhausted T cells (Tim3+) and exhausted suppressed T cells (Tim3+ within one cell radius of CD33+ cells; Supplementary Figure S9D). Like the scRNA-Seq data, the IF data showed a sizable proportion of T cells are exhausted/suppressed. Finally, we analyzed whether there were correlations between the tumor and immune cell type proportions in the IF data but found none of significance (Supplementary Data S7,10).
In summary, the immune cell analysis most prominently depicted the likely immunosuppressive role of tumor-associated macrophages and potential that microglia-associated macrophages drive the mesenchymal cell state. We also identified gene expression markers and IF staining patterns consistent with T-cell exhaustion and suppression.
Discussion
Our approach of grouping tumor cells into recognized cell types using supervised clustering and then using unsupervised clustering to characterize and probe differences within cell types facilitated identification of common subpopulations within a diverse set of tumors. We concluded that astrocytes, OPCs, and Mes cells related to OPCs are the primary contributors to oncogenesis and tumor growth, including angiogenesis and metastasis in PHGG. OPCs have previously been identified as a likely PHGG cell of origin.5–7 The OPCs in our study included a proliferating subpopulation with a stemlike signature and a small, slowly dividing stemlike subpopulation enriched in embryonic stem cell gene expression. These subpopulations may represent a neural stem cell-like population of tumor-initiating cells and a separate transit amplifying cell population. To our knowledge, ours is the first direct identification of likely stemlike populations in the OPC component of PHGG.
Astrocyte and mesenchymal cell populations also contributed to PHGG growth and progression. Astrocytes frequently depend on epidermal growth factor receptor (EGFR) signaling for cell growth in adult HGG7 but rarely do so in PHGG, where EFGR is not typically overexpressed. We found only a slight enrichment of EGFR in our astrocyte population. Metallothioneins, NF-κB pathway, and VEGF pathway genes were enriched in the astrocyte populations, suggesting participation in angiogenesis and metastasis. The Mes cell population also expressed NF-κB pathway and VEGF pathway genes, suggesting a potential role in angiogenesis.
The presence of OPC-like cells within the predicted Mes subpopulations and cells expressing mesenchymal genes within some of the OPC subtypes further suggested that the predicted Mes cells may derive from OPCs. In particular, the OPCs express the transcription factor ZEB2, an early transcription factor driving EMT that is required to initiate expression of vimentin, a key mesenchymal marker enriched in Mes cells.34
Interaction between immune and tumor cell populations suggests a role for immune cells in the development of mesenchymal gene expression. The microglial-associated macrophages expressed mesenchymal genes, they occurred more frequently as the Mes tumor cell population size increased, and they also occurred preferentially as nearest neighbors to Mes tumor cells. These associations suggest a causative link, possibly through induction of changes to the tumor microenvironment mediated by receptor-ligand or paracrine signaling interactions.
Intertumoral heterogeneity due to varying proportions of tumor cells with different oncogenic drivers suggests a potential explanation for inconsistent therapeutic results in treating PHGG. A solution may be to use single-cell transcriptomics to identify and design personalized treatments to address key oncogenic pathways that drive tumor growth and metastasis.46,47
Limitations of the current study include the number of samples available for analysis given the heterogeneous nature of PHGG. While we had at least four samples of most PHGG subtypes studied, these numbers proved to be limiting in identifying statistically meaningful differences in gene expression between subtypes. Diversity within each subtype further complicated identification of patterns of variation within the current PHGG subtype structure. A strength of the study design was the ability to identify predicted cell types across all PHGG subtypes, enabling us to propose an overall structure to identify and explain sources of heterogeneity within PHGG.
Key directions for future investigation include whether and how perivascular macrophages induce a mesenchymal tumor microenvironment and the mechanism by which tumor cells react to those changes. An additional key future question is how the cell types we identified respond to radiotherapy and chemotherapies. Therapeutic resistance is well documented in PHGG, but the mechanisms have so far been difficult to decipher.9,48 We hope the advances we have made in understanding tumor biology will accelerate identification of resistance mechanisms and methods to combat them.
Acknowledgments
The Anschutz Medical Campus Genomics Core performed bulk and single-cell RNA sequencing and DNA methylation analysis. The Anschutz Medical Campus Human Immune Monitoring Shared Resource performed multi-channel immunofluorescence staining and imaging.
Conflict of interest statement
The authors declare no competing interests.
Funding
Support was provided by grants to ALG from Alex’s Lemonade Stand Foundation and the Morgan Adams Foundation. The Anschutz Medical Campus Genomics Shared Resource is supported by National Cancer Institute University of Colorado Cancer Center grant [P30CA046934].
Data availability
Sequencing data are deposited in the NCBI Gene Expression Omnibus (GEO) database (GSE231860). The processed scRNA-Seq dataset is also available at Browse Projects - ScPCA Portal (alexslemonade.org). A browsable web page of the scRNA-Seq data is available at: https://www.pneuroonccellatlas.org/.
Authorship statement
JD and ALG designed the project; JD wrote the manuscript; JD and AD prepared samples for analysis; AMG assisted with design and interpretation of immune cell analyses; JD, RF, and KR performed bioinformatics analyses; KR constructed the cell browser; JAS analyzed and suggested analyses of the microglia and macrophages; ALG, JML, JAS, NKF and RV provided overall guidance for the project.
