Identification of a unique tumor cell subset employing myeloid transcriptional circuits to create an immunomodulatory microenvironment in glioblastoma.

Glioblastoma (GBM) is an aggressive primary brain tumor with unique immunity predominated by myeloid cells. GBM cells have been implicated to evade immune attack through hijacking myeloid-affiliated transcriptional programs to establish an immunosuppressive microenvironment. However, molecular features of immune-evading GBM cells in heterogeneous GBMs and their interactions with immune cells remain unclear. Herein, we employed single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq data to develop an in silico method for delineating GBM immune signature and identifying new molecular subsets for immunotherapy. We identified a new GBM cell subset, termed TC-6, that harbored immune-invading signature and actively interacted with tumor-associated macrophages (TAMs) to orchestrate an immune-suppressive niche. Proinflammatory transcriptional factors STAT1, STAT2, IRF1, IRF2, IRF3, and IRF7 were identified as the core regulons defining TC-6 subsets. Further immune transcriptome analyses revealed three immune subtypes (C1, C2, and C3). C3 subtype GBMs were enriched with TC-6 cells and immunosuppressive TAMs, and exhibited an immunomodulatory signature that associated with reduced efficacy of anti-PD-1 treatment. Interferon-related DNA damage resistance signaling was upregulated in C3 GBMs, predicting shortened survival of GBM patients who received chemo-radiation treatment. Treatment of OSI-930 as a molecular agent targeting c-kit and VEGFR2 tyrosine kinases may compromise the immunomodulatory signature of C3 GBMs and synergize with chemo-radiation therapy. We further developed a simplified 11-gene set for defining C3 GBMs. Our work identified TC-6 subset as an immune-evading hub that creates an immunomodulatory signature of C3 GBMs, gaining insights into the heterogeneity of GBM immune microenvironment and holding promise for optimized anti-GBM immunotherapy.

Introduction

Glioblastoma (GBM) is the most lethal and prevalent malignant brain tumor with marked cellular and molecular heterogenicity. GBM exhibits a dismal response to standard-of-care therapeutics including maximal safe surgical resection, radiation and chemo-treatment.[1] Genetic profiling of GBM has identified distinct transcriptional subtypes (termed proneural/classical/mesenchymal subtypes), reflecting the distinct states of GBM tumor cells.[2,3] Whether and how GBM heterogenicity is driven by the intrinsic evolution of tumor cells or the extrinsic microenvironmental cues remains largely unknown. A recent work by Gangoso etal.[4] reveals that GBM cells employ myeloid transcriptional circuits to establish a tumor-supportive, immunosuppressive microenvironment. However, molecular features of immune-evading GBM cells in heterogeneous GBMs and their interactions with immune cells remain unclear.

Immune checkpoint blockade (ICB) treatment reshapes the scope of anti-tumortherapeutics.[5] In GBMs, although preclinical results of ICB treatment are promising,[6,7] clinical trials are under evaluation with discrepant responsiveness.[8,9] Genetic changes of PTEN and MAPK signaling molecules, as well as the activation of WNT-β-catenin pathway in GBM cells, have been proposed as intrinsic immune evasion mechanisms.[8,10] Adaptive changes of tumor microenvironment, including induction of PD-1 expression and sustained immune activation, also contribute to neoadjuvant anti-PD-1 treatment inefficacy.[11,12] Likewise, a growing number of data imply that sustained interferon(IFN) signaling arose from PD-1 blockade is a major driving force of anti-tumor adaptive immune resistance.[13,14] In clinical settings, comprehensive analyses of the interactions between tumor-infiltrating immune cells and tumor cells could predict immunotherapy efficiency.[4,15] However, whether there are GBM cell populations that could drive immune evasion and create a specific immune signature associated with ICB ineffectiveness remains largely unknown.

Herein, we employed single-cell RNA-seqencing (scRNA-seq) and bulk RNA-seq data to develop an in silico method for delineating GBM immune signature, and identified a unique tumor subset (TC-6) as a hub for creating an immunomodulatory signature of human GBMs. TC-6 actively interacts with tumor microenvironmental (TME) cells and acquires TAM-affiliated transcriptional programs to elicit immunosuppression. To delineate specific transcriptional programs, we performed single-cell regulatory network inference and clustering (SCENIC) based on the scRNA-seq data to determine the core gene regulatory networks for immune subtyping. Three immunological subtypes were identified with distinct genomic alterations, immunotherapeutic response, and patient survival. Our work provides insights into the transcriptional modules in GBMs, with potential value for developing and optimizing therapeutic strategies against this lethal tumor.

Materials and methods

Human GBM cohorts

The Cancer Genome Atlas (TCGA) GBM gene expression data (Affymetrix Human Genome U133A microarray platform, n = 529) and the corresponding patient clinical information data were from the UCSC data portal (https://xenabrowser.net). Quantile normalization was performed on gene expression data using the limma R package. RNA sequencing (RNA-seq) raw data of primary GBMs (n = 154) and recurrent GBMs (n = 5) were obtained from the TCGA portal (http://tcga-data.nci.nih.gov/tcga/) and were normalized using Deseq2 R package. GBM gene expression profiles of Chinese patients were retrieved from the Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn/). The batch effect between mRNAseq-325 series and mRNAseq-693 series was corrected using the “ComBat” function of sva R package. RNA-seq data of recurrent GBM patients treated with adjuvant, post-surgical anti-PD-1 blockade were retrieved from GEO portal (GSE121810, https://www.ncbi.nlm.nih.gov/geo/).

Enrichment analysis

Gene-set enrichment analysis (GSEA) was performed using GSEA software (v4.0.1) (https://www.gsea-msigdb.org/gsea). Enrichment results with a criterion of p-value < 0.05 and false discovery rate (FDR) < 0.25 were considered as significant. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler R package. The files of oncogenic tumor biology were interpreted and visualized using the enrichplot R package. Single-sample GSEA (ssGSEA) algorithm was employed for analyzing the enrichment scores usingGSVA package. Gene sets were derived from the Broad Institute’s Molecular Signatures Database (MsigDB) or literature resources.

Immune signature analysis

CIBERSORT platform (https://cibersort.stanford.edu/) was used to estimate immune cell proportions in GBM dataset based on the LM22 reference profile. xCell (https://xcell.ucsf.edu/) gene signature-based method was used for calculating the proportions of Th1 and Th2 cells. Computing gene set scores for representative immune signatures (immune checkpoint, immunosuppression modulators for TC-6, immune response to tumor cells, et al.) were performed through ssGSEA algorithm. Immune-related indices and an IFN-related DNA damage signature have been defined in previous studies.[16,17] Immune-related functional modules used for GSEA analyses were summarized in Supplementary Table 1.

SynergySeq

SynergySeq platform (https://schurerlab.shinyapps.io/synergyseq/) was used for identifying compounds associated with inflammatory signature and with potential to combine with chemo-radiotherapy.[18] The LINCS compounds displaying high reversal scores of disease signature and high dissimilarity of cytotoxic therapeutics (chemoagent temozolomide (TMZ) and radiation) were considered as potential molecular candidates for combined treatment. The candidates with signature discordance score ≥ 0.5 and TMZ orthogonality score < 0.02 were selected for in silico compound library screening.

Single-cell RNA sequencing data analyses

GBM scRNA-seq dataset was downloaded from Sequence Read Archive (SRA, PRJNA579593). The CellRanger software (version 2.1.0) was performed to process single-cell data to demultiplex, align reads against the GRCh38 human reference genome, and generate feature-barcode matrices. Raw matrix was preprocessed using computational methods deposited in the Seurat R package. Quality control was performed to further ensure that only high-quality single-cell data was processed further, and cells with fewer than 4000 genes/cell and fewer than 8000 UMIs/cell were removed, as well were cells with greater than 20% of their transcriptome represented in mitochondrial transcripts. We applied the Canonical Correlation Analysis (CCA) algorithm to obtain the batch-corrected expression matrix to remove the batch effects from different donors. To cluster single cells by their expression profiles, we chose 20 principal components and used an unsupervised graph-based clustering algorithm Louvain (resolutions: 0.15). Cell type annotations were based on the expression of cellular entity gene marker visualized on the t-SNE plot. Gene signatures for different cell subsets were identified as log2 fold-change value > 0.5 and expression percent > 50% in all qualified cells (Supplementary Table 2). SCENIC (https://aertslab.org/) protocol was used for single-cell regulatory network inference.[19] The inference of regulons was performed as follows: (1) Identification of co-expression modules between transcription factors and the potential targets based on the gene-expression matrix through GENIE3; (2) Screening for the significantly enriched motif of the correct upstream regulators to identify truly functional regulons through cis-regulatory motif analysis from RcisTarget; (3) Calculation of the activity score of each regulon using AUCell algorithm. The regulon-activity matrix for each cell was used for dimensionality reduction for visualization via the t-SNE technique. For cell-type specific regulons identification, calcRSS function was performed based on regulon specificity score. For trajectory inference, we used the spliced and unspliced transcript reads to calculate the RNA velocity using velocity.Rscript with the Cell Ranger output. All the tumor cell subsets were analyzed and the evolution-based trajectory was projected onto the UMAP plot.

Cell-cell communication inference

The intracellular interactions through receptor-ligand complexes were analyzed through CellphoneDB (www.cellphonedb.org).[20] The average expression level of ligands and receptors in each cell type was generated and presented as a null distribution of each ligand-receptor pair. The likelihood of cell-specificity of a given receptor-ligand complex was determined.

GBM immunological subtyping

Genes in the metagenes (n = 112) were used for GBM immunological subtyping using R package(parameter: pItem = 0.8, pFeature = 1, reps = 1000) with Euclidean clustering method. The best cluster number (k) was determined using Consensus Cumulative Distribution Function (CDF) Plot. The t-SNE-based approach was employed for subtype assignment validation using the Rtsne R package. The molecular similarity between the predicted subtypes from the TCGA dataset and those from the CGGA dataset was compared using Subclass Mapping on the platform of GenePattern.

Generation of the C3 gene classifier

Prediction Analysis for Microarrays (PAM)[21] was employed for identifying subset-associated genes using the nearest shrunken centroids through the PAMR package. GBM samples from the TCGA database were randomly stratified into the training and validation cohorts with a ratio of 2:1. The PAM fits the model for 100 different threshold values, 10-fold cross-validations were used to choose the optimal threshold value for centroid shrinkage. The accuracy of the 11-gene classifier was calculated and validated through external GBM datasets.

Genomic mutation analysis

The mutation annotation format (MAF) file with aggregated mutation information of GBM cases was obtained from the TCGA portal. Candidate somatic single nucleotide variations were identified using Mutect algorithms. The mutational summary and landscape plots were performed by the maftools R package. The “mafCompare” function was applied to compare two different cohorts to identify differentially mutated genes, and results were visualized using forest plots. The “clinicalenrichment” function was used for various groupwise and pairwise comparisons to identify enriched mutations and copy number alterations (CNAs) in each subtype. The “plotEnrichmentResults” function was used to plot these results. MutSigCV1.41 on the GenePattern platform was used to identify the mutated genes highly relevant to cancer (q-value <0.25).

Immunohistochemical (IHC) and immunofluorescent (IF) staining

Human GBM specimens were from the Biobank of Southwest Hospital, Army Medical University, with informed consent from patients or their guardians. The usage of human specimens was approved by the Ethics Committees of Southwest Hospital and was in accordance with the principles of the Helsinki Declaration. IHC and IF staining were performed on formalin-fixed, paraffin-embedded GBM sections as previously described.[22] Primary antibodies used in this study were listed in Supplementary Table 3. Dako REAL EnVision Detection System was used for visualization following the manufacturer’s protocol. IHC staining was evaluated independently by two neuropathologists blinded to C3 status of each tumor. The expression level of each molecule was semi-quantitatively scored based on the staining area and staining intensity of positive tumor cells as previously described.[23] The proportion of positive tumor cells (staining area) was scored as follows: 0 (no positive tumor cells), 1 (<5% positive tumor cells), 2 (5%–19% positive tumor cells), 3 (20%–50% positive tumor cells), 4 (>50% positive tumor cells). The staining intensity was evaluated according to the following criteria: 0 (no staining), 1 (weak staining = light yellow), 2 (moderate staining = yellow brown), and 3 (strong staining = brown). The staining index (SI) was calculated as staining intensity score × staining area score for each field of microscope.

RNA-sequencing (RNA-seq)

The RNA of GBM specimens from Southwest Hospital was extracted using TRIzol (Invitrogen). The cDNA libraries were constructed according to the TruSeqRNA Sample Preparation v2 protocol, and sequenced using an Illumina HiSeq X sequencer in paired-end sequencing mode (2 × 150 bp reads). The HISAT alignment tool was used to align all sequencing raw reads to the human genome (hg38) reference sequence. Only uniquely and properly mapped read pairs were used for further analysis. Cufflinks program was used to assemble individual transcripts from reads aligned to the human genome.

Statistical analyses

Statistical analyses were performed using R software (3.6.3) or SPSS (19.0 Inc, Chicago, IL). Kruskal–Wallis H test was used for comparisons of variables with abnormal distribution among three independent groups. The Mann–Whitney U test was used for the comparison between two independent groups when the dependent variable was continuous but was not normally distributed. The Chi-Square test was used for comparison of categorical variables. Kaplan–Meier survival curve, Log-rank test and multivariate Cox regression analyses were used for survival analyses. The forest plot was used for the visualization of covariate effects using the forest model R package. P< 0.05 or 0.01 was considered statistical significance.

Results

Single-cell RNA sequencing identifies a unique tumor cell subset orchestrating immunosuppression

To delineate GBM immune TME feature at the single-cell level, we used scRNA-seq data from 5 cases of IDH wild-type human GBMs (SRA, PRJNA579593) and performed unsupervised graph-based clustering on 36789 cells. Sixteen subsets were identified and were annotated based on CNAs and canonical cell markers (Figure 1a). Glioma cell lineages were characterized by the gain of chromosome 7 and the loss of chromosome 10.[24] Normalization of CNA profiles revealed chromosomal aberrations in each cell and 6 tumor cell subsets (termed TC-1 to TC-6) were identified in all the samples analyzed (Supplementary Figure S1a-b). TC-1 to TC-5 subsets recapitulated the features of neural progenitor-like (NPC-like), oligodendrocyte-progenitor-like (OPC-like), astrocyte-like (AC-like), and mesenchymal-like (MES-like) GBM cells identified in previous studies,[2] whereas TC-6 did not match the feature of abovementioned subsets. Of note, TC-6 showed high expression of EGFR and IFN-stimulated genes (CXCL10, STAT1, ISG15, IFI44L) (Figure 1b). Functional enrichment analyses of hallmark collection indicated that TC-6 was associated with signatures involving IFN response, epithelial-mesenchymal transition, hypoxia, and inflammatory response (Figure 1c). To address the interplay between these tumor cell subsets and immune cells, we performed an unbiased ligand-receptor interaction analysis among these TC subsets and niche cells (Supplementary Figure S2a, Figure 1d). In comparison with other TC subsets, TC-6 subset actively interacted with TAMs (TAM-1, -2) through predicted pathways involving TAM recruitment and phenotype reprogramming (CX3CR1-CX3CL1 axis; CSF1-CSF1R axis), growth factor-initiated tumor growth (MIF-EGFR axis, Supplementary Figure S2b) and “Do not Eat Me” immune surveillance (LILRB1-HLA-F axis).[25,26] TC-6 subset was associated with abundant inflammatory stimuli (IL1β, IL1RN, MIF and Oncostatin-M) mediating the TAM reprogramming (Figure 1e). Further analyses of TCGA GBM cohort indicated that TAM signature was positively correlated with TC-6 signature (Figure 1f). In silico cellular interaction analyses also suggest TC-6 subset may interact with tumor-infiltrating T cells in GBMs to regulate T cell recruitment (CXCR3-CXCL10 axis; SPN-ICAM1 axis), expansion and cytotoxic activity (CD94:NKG2A-HLA-E axis and TIGIT-NECTIN2 axis).[27-29] Intriguingly, in comparison with other tumor cell subsets, TC-6 exhibited an upregulation of immune checkpoints PD-L1, IDO1, NT5E (CD73), and LGALS9 for immune evasion[30-32] (Figure 1g). RNA velocity estimation indicated that TC-6 was the precursor cells of other GBM cell populations (Supplementary Figure S3a). Altogether, these data suggest TC-6 subset as a unique tumor cell subset that could actively interact with TAMs and T cells to create an immunosuppressive microenvironment.

Figure 1.

Single-cell RNA sequencing identifies a novel tumor cell subset for orchestrating immunosuppression. (a) t-SNE plot showing 16 major cell types using the single-cell RNA sequencing data from 5 cases of GBMs. TAM, Tumor-associated macrophages; RBC, Red blood cell; OPC, Oligodendrocyte progenitor cell. (b) Dot plots showing the distinct cell state markers in the indicated tumor cell subsets (TC-1, -2, -3, -4, -5, and -6). Dot size was proportional to the fraction of cells expressing specific genes. Color intensity corresponds to the relative expression of specific genes. TC, Tumor cell; NPC, Neural progenitor cell; AC, Astrocyte; MES, Mesenchymal. (c) Enrichment analyses of the differential expressed genes in TC-6 versus other subsets (log2 fold-change > 0.5, min.pct > 0.10). The bubble size indicates the number of genes in each term, and different colors correspond to different adjusted p-values. (d) Heatmap showing the total number of interactions between cell types inferred from whole clusters of GBMs using CellphoneDB. Oligo., Oligodendrocyte. (e) Overview of representative ligand–receptor interactions between tumor cells and immune cells. P-values are indicated by circle size, the average expression of interacting molecules in cluster 1 and 2 are indicated by colors. (f) Correlation plot between TC-6 subset and TAM subset in TCGA GBM cohort, correlation coefficient was listed on the top. (g) Dot plots showing the expression of known immune checkpoints (PD-L1, IDO1, NT5E, and LGALS9) across the TC subsets.

Identification of core regulons coordinating the TC-6 phenotype

To identify transcriptional regulatory networks determining the state of each cell subset, we applied single-cell regulatory network inference and clustering (SCENIC) to evaluate the activity of the gene regulatory networks (GRNs, termed regulons). Eighty regulons were recognized across different cell subsets (Figure 2a). To evaluate the performance of SCENIC, annotations inferred from the expression profiling were mapped to the t-SNE clustering of regulon activity matrix (Supplementary Figure S3b). Cell entity clustering was similar to those defined by transcriptome profiling. Among the 88 regulons, IRF1, IRF2, IRF3, IRF7, STAT1, and STAT2 were identified as specific regulators for TC-6 based on the Regulon Specificity Score (RSS) (Figure 2b).[33] TC-6 also showed high gene set activities in these regulons (Figure 2c) that may mediate unique immune transcriptional programs.

Figure 2.

Identification of core regulons delineating the TC-6 phenotype. (a) Heatmap showing the normalized activity score of representative regulons across different subsets of GBM tissues. (b) Identification of TC-6 subset specific regulators based on the Regulon Specificity Score (RSS, y-axis). (c) Activities of IRF1, IRF2, IRF3, IRF7, STAT1 and STAT2 regulons on t-SNE clustering, cells were colored according to the corresponding activity score of regulon (upper panels). Dark blue dots representative of TC-6 in the t-SNE single-cell map were shown in the bottom panel.

Novel GBM subtyping based on TC-6 core regulons

We next investigated whether the core regulons of TC-6 could be used to define GBM immune subtypes. The workflow was shown in Figure 3a. The top 30 genes with high Genie3weight were screened and 112 candidate genes were identified as metagenes (Supplementary Table 4). TCGA GBM dataset was clustered based on the metagenes using consensus clustering. The results showed that three-subtype-clustering was the most optimized clustering method (Figure 3b, Supplementary Figure S4a, b), and was consistent with the clustering of spectral patterns through t-SNE dimensionality reduction (Figure 3c). The rationality and reproducibility of the three-subtype-clustering method were validated through the CGGA GBM cohort (Supplementary Figure S4c,d). Most metagenes associated with TC-6 subtype were upregulated in C3 subtype (Supplementary Figure S4e), and TC6 signature score was significantly higher in C3 subtype of GBMs as compared with C1 or C2 subtype of GBMs (Supplementary Figure S4f), indicating C3 subtype of GBMs are enriched with TC-6 subset. To determine the prognostic significance of immunological subtypes, we employed TCGA and CGGA GBM cohorts and found that C3 subtype GBMs displayed shortened survival relative to non-C3 subtype GBMs (Figure 3d,e).

Figure 3.

Novel GBM subtyping based on TC-6 core regulons. (a) The workflow of TC-6 regulons-based subtyping of human GBMs. (b) Consensus clustering based on metagenes from core regulons identified three subgroups. (c) t-SNE analysis of the three subtypes. (d-e) Kaplan–Meier survival analyses of overall survival in GBMs with C3 and non-C3 subtypes using the TCGA mRNA cohort (d) and CGGA RNA-seq cohort (e). **, p < 0.01; *, p < 0.05.

IFN-related DNA damage resistance signature informs therapeutic resistance to chemo/radio-therapy in C3-subtype GBMs

We further delineate the molecular features and clinical characteristics of the three subtypes. C3 subtype of GBMs was associated with wild-type IDH1, non-glioma-CpG Island methylator phenotype (non-G-CIMP), and classical subtype features of glioma samples (Figure 4a) . ssGSEA analysis revealed distinct functions of C1-C3 subtypes. Intriguingly, C3 subtype was associated with IFN-related DNA damage resistance (IRDS), cytosolic DNA sensing pathway, and inflammatory response (Figure 4a). Meanwhile, C1 subtype was enriched with G2/M checkpoints and Notch signaling gene signature. C2 subtype was associated with pathways regulating hypoxia, angiogenesis, and glycolysis. To address how C3 molecular feature underlies worse clinical outcomes of GBM patients who received chemo/radiation, we analyzed IRDS gene set associated with chemo/radio-resistance.[16,34] GSEA analysis revealed that IRDS was significantly enriched in C3 subtype relative to non-C3 subtypes (Figure 4b). Twenty-four of 44 IRDS genes were among the top 10% genes differentially expressed in C3 versus non-C3 subtype (Figure 4b). Consistently, the majority of the C3-associated leading-edge genes in IRDS were upregulated in recurrent GBMs relative to the paired primary GBMs (Figure 4b). Besides, C3 GBMs overexpressed another geneset associated with chemo/radio-resistance (Figure 4c). To investigate whether targeting C3 signature could reverse the resistance to chemoagent TMZ and radiation thus improve therapeutic efficacy, we used the SynergySeq platform and screened for the potential agents targeting C3-subtype GBMs. Orthogonal plot analyses identified potential candidate agents that targeted C3 signature and could combine with TMZ/radiation (Figure 4d, Supplementary Table 5). OSI-930, a brain-penetrable molecular inhibitor of c-kit, VEGFR2 and CSF1 R,[35] was identified as a potential C3-targeting agent. We also identified EX-527, a SIRT1 inhibitor exhibiting anti-tumor effects on glioma cells,[36] as a promising chemo-agent for C2-subtype GBMs (Supplementary Table 6).

Figure 4.

IFN-related DNA damage resistance signature informs therapeutic resistance to chemo/radio-therapy treatment in C3-subtype GBMs. (a) The landscape of clinical and molecular characteristics among different subtypes (C1, C2, C3). The normalized ssGSEA score for functional gene sets was plotted through a heatmap. G-CIMP, Glioma-CpG island methylator phenotype; MGMT, O6-methylguanine-DNA methyltransferase; ECM, Extracellular matrix. (b) GSEA plot of the enrichment of IFN-related DNA damage resistance signature in the C3 subtype relative to the non-C3 subtypes (C1 and C2). Leading-edge genes were labeled as the red dotted line (left panel). Heatmap in the right panel showed the expression of these leading-edge genes in five pairs of primary-recurrent GBM samples from the TCGA database. (c) GSEA plot of C3 over-expressed geneset (log2 fold-change >1) in temozolomide (TMZ)/radiation-resistant group relative to the TMZ/radiation-sensitive group. (d) Rankings of the library of integrated network-based cellular signatures (LINCS) compounds based on their concordance with TMZ/radiation therapy and discordance with C3 signature. The x-axis value suggests the concordant degree with TMZ/radiation and the y-axis value suggests the discordance with C3 signature.

In the revised Figure 3, “10-gene classifier” has been replaced by “11-gene classifier” in panel a.

In the revised Supplementary tables, “LINCS IID” has been replaced by “LINCS ID” in Table S5 and Table S6.

Genomic landscape of the three immunological subtypes

To determine whether the three immunological subtypes were associated with genomic changes of GBM cells, we evaluated the distribution of somatic gene alterations among the three subtypes using MutSig algorithms. As shown in Supplementary Figure S5a–c, gene mutations associated with GBM progression (EGFR, PTEN, TP53, PDGFRA, NF1, RB1, PIK3CA, and IDH1) varied across C1-C3 subtypes. Of note, C3 subtype had the highest genomic alteration (mutation and regional DNA amplification) rate of EGFR and PIK3CA (70%, 16%) relative to C1 (45%, 6%) and C2 (52%, 12%) subtypes. HYDIN mutation, a recently identified cancer-associated antigen involving adaptive immune responses,[37] was associated with poor prognosis of GBM patients and showed higher mutation frequency in C3 subtype GBMs relative to C1 and C2 subtypes (Supplementary Figure S5b–d). EGFR amplification was also significantly higher in C3 (Supplementary Figure S5e). We determined the protein level and kinase phosphorylation of EGFR across all three subtypes using TCGA datasets. As expected, C3 subtype GBMs exhibited a higher expression of EGFR and phosphorylated EGFR relative to C1 and C2 subtypes (Supplementary Figure S5f), which may sustain cell survival and compromise chemo-radiation-induced cytotoxicity.

In the revised Figure 4, “CIMP” has been replaced by “G-CIMP” in the top left corner of panel a.

The label of y-axis has been added in the revised Figure S5c.

Distinct immunological characteristics of the immunological subtypes

We further profiled the immune milieu and found that C1-C3 subsets manifested distinct immune signatures (Figure 5a) and different immune cell fractions (Supplementary Figure S6a–c). C1 subtype was defined as a lymphocyte-dominant subtype, harboring increased Th1 lymphocytes (Figure 5a) and CD8+ T lymphocytes (Supplementary Figure S6a–c). High indel-derived neoantigens and immune response to tumor cells in C1 subtype informed increased probability of immune responsiveness to immune checkpoint inhibitors (Figure 5a). C2 subtype was defined as a myeloid-predominant subtype, showing an increased fraction of monocytes, macrophages and neutrophils (Figure 5a, Supplementary Figure S6a–c). C3 subtype was featured by the enrichment of polarized macrophages (i.e., tumor-associated macrophages) and NK cells (Supplementary Figure S6a–c), reduced indel neoantigens and cancertestis antigens, increased expressions of immune checkpoint molecules and immunosuppressive modulators overexpressed in TC-6 (Figure 5a). Correlation analyses of TCGA GBM cohort revealed an association between immunosuppressive TAM signature and TC-2, -3 and -6 signature (Supplementary Figure S7a,b). The immunological characteristics of C1-C3 subtype GBMs were summarized in Figure 5b.

Figure 5.

Immunological and clinical characteristics of GBMs with distinct immunological subtypes. (a) Scoring or fraction of immune cell components and indices in C1, C2 and C3 subtypes. (b) The key immune signatures among C1, C2 and C3 subtypes. IFN, Interferon; TCR, T cell receptor; Tfh, T follicular helper.

C3 subtype GBMs manifest reduced efficacy to anti-PD-1 therapy

Recent clinical trials of neoadjuvant anti-PD-1 treatment in GBMs showed a modest effect in improving patient survival.[7,8] To determine whether the efficacy of anti-PD-1 treatment is associated with C1-C3 immunological subtypes, we investigated a GBM cohort (GSE121810) treated with anti-PD-1[12] and assigned the GBM cases to immunological subtypes. As expected, C3-subtype GBMs showed a shortened overall survival (median: 87.5 days) relative to C2-subtype GBMs (median: 236.9 days) or C1-subtype GBMs (median: 329.0 days). Kaplan–Meier survival analyses confirmed that C3-subtype GBMs exhibited worst patient outcome relative to C1- or C2-subtype GBMs (Figure 6a; C3 vs. C1: p = 0.0067; C3 vs. C2: p = 0.002; C2 vs. C1: p = 0.1875). C3 subtype GBMs were enriched for pathways regulating IFN-γ response, T cell exhaustion, TGF-β signaling and TC-6-associated immunosuppressive modulators (Figure 6b,c), thus informing a reduced efficacy to anti-PD-1 therapy.

Figure 6.

C3 subtype predicts poor immunotherapeutic response. (a) Kaplan–Meier analyses of patient overall survival of anti-PD-1-treated GBMs with different subtypes (GSE121810). (b) GSEA analysis of enriched immune modules with statistical significance in C3 subtype relative to the non-C3 subtypes. The column plot (left panel) showed the normalized enrichment score of each module. Concordance index (C-index) was presented for the indicated characteristic immune module (right panel). Predictive power (high to low) was indicated as bar color (red to blue). (c) GSEA plots of the representative enriched functional modules in C3 subtype. NES, Normalized enrichment score.

Development of an 11-gene C3 classifier for GBM molecular diagnosis

To simplify a biomarker set defining C3 GBMs for molecular diagnosis, we used the PAM method to develop a multigene classifier. TCGA GBMs were randomly stratified into training/testing cohorts at the ratio of 2:1, and 10-fold cross-validation was used to compute the cross-validation error (Supplementary Figure S8a). An 11-gene classifier was identified as the best optimal set, showing the best sensitivity (0.978), accuracy (0.989), specificity (0.992) and the highest discriminant performance using the testing cohort (Figure 7a, Supplementary Figure S8b,c). The efficiency of the 11-gene C3 GBM classifier was confirmed by testing cohort, external TCGA RNA-seq testing dataset and CGGA RNA-seq testing dataset (Figure 7b,c). The immune functional modules associated with IFN-γ response, antigen processing and presentation were enriched in C3-subtype GBMs using the 11-gene-classifier (Supplementary Figure S8d). C3-subtype GBMs identified through the 11-gene classifier showed worse outcomes relative to the non-C3 GBMs using the TCGA and CGGA GBM cohorts (Figure 7d,e). Multivariate Cox regression analysis confirmed that the 11-gene C3 classifier was a promising and independent biomarker set for predicting GBM patient survival (Supplementary Figure S9a). Intriguingly, most of the classifier genes were upregulated in TC-6 subsets with good consistency, suggesting that TC-6 subset is a determinant for defining C3 immune signature (Supplementary Figure S9b). We analyzed bulk RNA-seq data from six cases of human GBMs and found a consistent upregulation of the 11 classifier genes in C3 subtype GBMs (Figure 7f). Moreover, the TC-6 and IFN-related DNA damage resistance signatures were significantly enriched in C3-subtype GBMs relative to non-C3-subtype GBMs (Figure 7g). IHC staining confirmed upregulated expression of 11 classifier genes (Figure 8a,b), TAM-1 marker (CD163), TAM-2 markers (CTSD and GPNMB) and EGFR (Supplementary Figure S10a,b) in C3-subtype GBMs relative to non-C3-subtype GBMs. These data indicate that the 11-gene classifier is sufficient for defining C3-subtype GBMs with potential value for molecular diagnosis and prognostic determination.

Figure 7.

Development of an 11-gene C3 classifier for GBM molecular diagnosis. (a) Plots of sensitivity (red) and specificity (black) of the shrinkage parameter (thresholds) computed by cross-validation. The value (7.9, indicaded in blue line) yielded a subset of 11 genes with the most optimized efficiency. (b) Percentage of classified samples using the 11-gene classifier for C3 subtype in the indicated GBM datasets. The TCGA mRNA (training set) was employed as the predictor model. TCGA mRNA dataset, TCGA RNA-seq dataset, CGGA RNA-seq dataset were used as the testing sets. (c) Unsupervised hierarchical clustering of the TCGA mRNA testing set based on the 11-gene classifier. Metagenes-based subtypes are indicated with different colors at the top column. (d-e) Kaplan–Meier analyses of GBM patients from the TCGA dataset (d) and CGGA dataset (e) based on the 11-gene classifier of C3 subtype. (f) Unsupervised hierarchical clustering of six human GBM samples of RNA-seq data based on the 11-gene classifier of C3 subtype. (g) Enrichment score of TC-6 signature and IFN-related DNA damage resistance in human GBMs with non-C3 and C3 subtypes. The data was presented as the mean ± SEM. **, p < 0.01; *, p < 0.05.

Figure 8.

IHC staining of 11 classifier genes in C3-subtype GBMs relative to non-C3-subtypes. (a) Representative IHC images showing the expressions of 11 classifier genes in human GBMs with C3 and non-C3 subtypes. Scale bar: 50 μm. (b) Semi-quantification of IHC staining of 11 classifier genes in human GBMs with C3 and non-C3 subtypes. ***, p < 0.001; **, p< 0.01; *, p < 0.05.

In the revised Figure S8b and c, “C1&C2” has been replaced by“Non-C3”.

In the revised Figure 7, “C1&C2” has been replaced by“Non-C3”in panel d and panel e.

Discussion

The limited understandings of GBM tumor cell heterogeneity and the complex interaction with immune components impede the development of effective immune therapeutics. Herein, we employed scRNA-seq and bulk RNA-seq data to develop an in silico method for delineating GBM immune signature. We identified TC-6 subset as a hub for creating an immunomodulatory signature of C3 subtype GBMs, holding promise for advancing molecular diagnosis of GBM immune feature and subsequent immunotherapy development. We speculate that the generation of TC-6 subset is not elicited via intrinsic genetic evolution, but through immunoediting reprogramming by immune cells and niche factors such as IFN. A recent work by Gangoso et al.[4] indicated that IFN-γ treatment induced IRF8 expression in GBM cells to contribute to immune evasion. In agreement with this, our results showed that IRF1, IRF2, IRF3, IRF7, STAT1, and STAT2 downstream of IFN-γ signaling were highly activated in TC-6 subset to launch a myeloid-affiliated transcriptional program. While the sources of IFN-γ in GBMs warrant further exploration, emerging data have demonstrated that tumor-infiltrating macrophages especially the proinflammatory macrophages and brain-residing microglia secreted abundant IFN-γ.[38] We speculate that IFN-γ signaling may represent a sustaining circuit that enables the retrospective interaction between TC-6 cells and tumor-infiltrating myeloid cells, thus rendering tumor cells undergo myeloid-affiliated transcriptional reconfiguration and creating an immunosuppressive environment in C3 subtype GBMs. Apart from immune cells, our in silico analyses also suggest a possible interplay between TC-6 and vascular cells through EFGR and VEGF signaling. Vascular components are important participants of GBM immunity. Besides, glioma-vascular interactions facilitate vasculogenic mimicry, vessel co-option mediated invasion, and glioma stem cell stemness.[39-41] Future studies are warranted to investigate whether and how TC-6 could interact with vascular cells to reshape immune cell recruitment and functions.

GBMs are typically T cell-deprived immune cold tumors that may not respond to ICB therapy. The molecular signatures related to unfavorable ICB response in GBMs remain largely unknown. The enriched IFN signature in C3 subtype GBMs prompts us to investigate its association with anti-PD-1 efficiency. IFN signaling has been implicated in orchestrating antitumor responses, whereas conflicting data exist regarding the role of IFN signaling in restraining ICB treatment.[42,43] Recent evidence supports that induction of IFN-γ at the low tumor burden stage could promote apoptosis of tumor-cytotoxic CD8+ T cells, resulting in a compromised long-term response of memory T cells and impaired ICB therapeutic efficacy.[44] Besides, tumors may occupy a natural defensive mechanism from immune response.[45,46] IFN-γ could also exert an immune-suppressive function by inducing the expression of immune inhibitory factors like PD-L1/2, IDO1 and TIM3 to antagonize innate and adaptive immune responses in solid tumors.[10,13,45,47] Our work identified upregulated PD-L1 and IDO1 in TC-6 subsets, which might be attributed to sustained inflammatory stimuli like IFNs in C3 subtype GBM microenvironment. Delineating immune features in GBMs received neoadjuvant anti-PD-1 are warranted in future works to confirm our in silico findings of the role of C3 signature in predicting anti-PD-1 response.

Previous Transcriptomic analyses revealed proneural, classical, and mesenchymal subtypes as the major GBM molecular clusters.[48] One of the major limitations of this molecular classification system is the limited involvement of environmental cues, especially the immune-competent signature. Our classification scheme integrates immune cell entity and tumor-immune cell interaction analyses to make a supplement to Verhaak’s GBM subtyping to delineate GBM molecular features. We also integrate tumor genomic alternations such as IDH1, NF1 and EGFRvIII in our analyses, as genomic mutation patterns are fundamental for inflammatory factor production, peripheral immune cell recruitment and activation.[8,49,50] EGFR amplification and mutation were frequently observed in C3 subtype GBMs and were predominately presented in TC-6 cells. Previous studies have demonstrated that EGFR activation impaired anti-tumoral immunity in lung carcinoma and melanoma.[51-53] GBM shares many of the molecular drivers used by other solid tumors, but possesses a unique immune environment that renders immunotherapy failure.[46] It is worthwhile to investigate the unique mutational patterns in C3 subtype GBMs and their association with immune reprogramming.

In conclusion, we developed a novel, reliable, and feasible subtype classification system for delineating GBM immune signatures. A simplified 11-gene classifier defining C3-subtype GBMs shows the potential value for prognostic determination and anti-PD-1 efficacy prediction. Our findings of a unique GBM cell subset employing myeloid transcriptional circuits to create an immunomodulatory signature delineate the interplays of GBM cell subsets and infiltrating myeloid cells, with potential value for GBM molecular diagnosis and immunotherapy development.

Click here for additional data file.