Distinct immune escape and microenvironment between RG-like and pri-OPC-like glioma revealed by single-cell RNA-seq analysis

Abstract

The association of neurogenesis and gliogenesis with glioma remains unclear. By conducting single-cell RNA-seq analyses on 26 gliomas, we reported their classification into primitive oligodendrocyte precursor cell (pri-OPC)-like and radial glia (RG)-like tumors and validated it in a public cohort and TCGA glioma. The RG-like tumors exhibited wild-type isocitrate dehydrogenase and tended to carry EGFR mutations, and the pri-OPC-like ones were prone to carrying TP53 mutations. Tumor subclones only in pri-OPC-like tumors showed substantially down-regulated MHC-I genes, suggesting their distinct immune evasion programs. Furthermore, the two subgroups appeared to extensively modulate glioma-infiltrating lymphocytes in distinct manners. Some specific genes not expressed in normal immune cells were found in glioma-infiltrating lymphocytes. For example, glial/glioma stem cell markers OLIG1/PTPRZ1 and B cell-specific receptors IGLC2/IGKC were expressed in pri-OPC-like and RG-like glioma-infiltrating lymphocytes, respectively. Their expression was positively correlated with those of immune checkpoint genes (e.g., LGALS3) and poor survivals as validated by the increased expression of LGALS3 upon IGKC overexpression in Jurkat cells. This finding indicated a potential inhibitory role in tumor-infiltrating lymphocytes and could provide a new way of cancer immune evasion.

Keywords

single-cell RNA-seq / glioma / radial glia / primitive oligodendrocyte precursor cell / immune escape

1 Introduction

In the past 20 years, cancer molecular subtypes have been used as basis for prognosis and treatment. More than 80% of brain tumors are gliomas and adult diffuse gliomas that can be classified into the following three molecular subtypes according to the 2021 WHO classification of tumors of the central nervous system [1]: glioblastoma (GBM) with isocitrate dehydrogenase (IDH) wild-type (named as G1 here), astrocytoma with IDH mutation and no chromosome arms 1p and 19q (1p/19q) codeletion (named as G2 here), and oligodendroglioma with IDH mutation and 1p/19q codeletion (named as G3 here). These subtypes have different clinical outcomes [2]. Despite decades of efforts from the glioma community, the understanding of their biology and the resulting development of targeted therapies remain limited.

In the development of the central nervous system, early radial glia (RG) cells arise from neuroepithelial cells lining the ventricles during neurogenesis [3]. RGs cannot only produce intermediate progenitor cells that ultimately generate neurons and provide scaffold to support neuronal migration during human cortical neurogenesis peak [4,5] but also generate oligodendrocyte precursor cells (OPCs) via intermediate primitive OPCs (pri-OPCs) at the onset of human cortical gliogenesis [6]. RG cells express classic markers, such as vimentin (VIM), pri-OPC cells have oligodendrocyte transcription factor 2 (OLIG2), and OPC cells have protocadherin-related 15 (PCDH15), platelet-derived growth factor receptor alpha (PDGFRA), and NK2 homeobox 2 (NKX2.2) as markers [6].

A tumor is a developing organ–tissue system [7]. Abnormal gliogenesis contributes to glioma tumorigenesis [8]. Single-cell analyses of neonatal glial development in a murine model of GBM revealed that pri-OPCs contribute to gliomagenesis and are prominent in all three human glioma subtypes [9]. In addition to gliogenesis, neurogenesis contributes to glioma progression [10,11]. A previous study compared single-cell transcriptomes between adult GBM cells and normal human fetal brain cells covering neurogenesis and gliogenesis and found that most GBM cancer cells were similar to glial progenitor cells (GPCs) [12].

Considering that G1–G3 differ in genomics and clinical outcomes, we intuitively hypothesized that the three subtypes may have different subgroups, though they are transcriptionally similar to one glial cell type, either pri-OPC or GPC. In this work, we analyzed 26 tumors representing all glioma molecular subtypes by conducting single-cell RNA-seq. We found that glioma could be classified into pri-OPC-like and RG-like subgroups with high and low fraction of cancer cells expressing oligo-lineage markers, respectively. This classification was validated in a public cohort of single-cell RNA-seq data set (called Broad cohort hereafter) [13] including 25 primary GBMs and gliomas in TCGA. Further characterization of cancer-driving mutated gene, MHC-I, in tumor subclones and immune microenvironment revealed significant difference between pri-OPC-like and RG-like tumors. We found that they could extensively modulate glioma-infiltrating immune cells in a different manner, for example, glial/glioma stem cell markers OLIG1/PTPRZ1 (i.e., not expressed in immune cells under normal conditions) and B cell-specific receptors IGLC2/IGKC (i.e., only express in B cells in normal conditions) are expressed in pri-OPC-like and RG-like glioma-infiltrating lymphocytes, respectively. Their expression is positively correlated with those of immune checkpoint genes (e.g., LGALS3) and poor survivals, suggesting their potential inhibitory role in tumor-infiltrating lymphocytes that shall provide a new way of cancer immune evasion.

2 Methods

2.1 Patient information

All patients involved in this study were from Fudan University Affiliated Huashan Hospital in China and have given consent preoperatively in accordance with IRB #2018-220 approved by the Ethics Committee of Huashan Hospital. Twenty-six patients were involved and divided into three groups, that is, grade IV glioma with WT IDH (G1, n = 11), astrocytoma with mutant IDH and no 1p/19q co-deletion (G2, n = 10), and oligodendroglioma with mutant IDH and 1p/19q co-deletion (G3, n = 5).

2.2 Glioma single-cell preparation

Tumors were freshly collected and transferred in cold DMEM/neurobasal medium to the laboratory within 30 min after resection. After a brief rinse with sterile DPBS and the removal of blood vessels and meninges, a part of the tumor was fixed for subsequent histology, and the remaining was finely minced with a sterile blade. Owing to the waiting period for pathology results for our grouping, some of the finely minced tumor tissues were cryopreserved in a freezing medium (90% FBS and 10% DMSO). Others were directly processed for single-cell dissociation.

The minced tumor tissues were dissociated into single cells using the method of Zhang et al. [14] with some modifications. The minced tissues were digested with papain (Worthington, LK003176) at 34 °C for 100 min with aeration of 95% O₂ plus 5% CO₂ and occasional agitation. The digested tissues were then washed with 5 mL of Low Ovo [14] (ovomucoid with low concentration, also known as trypsin inhibitor) for four times and gently triturated twice with 4 mL of Low Ovo without introducing any bubble. The supernatant was collected into a new tube, and 10 mL of High Ovo [14] was carefully layered under the supernatant. The cells in the supernatant were spun down through High Ovo at 100× g for 5 min, and the cell pellet was then resuspended with 9 mL of DPBS and 0.02% BSA. After being filtered through a 40 μm cell strainer and spun down, the cells were resuspended with 3 mL of 30% Percoll (Sigma, P1644) and centrifuged at 700× g for 10 min to remove myelin debris and dead cells. The final single-cell suspension was resuspended with 200 μL–1 mL of DPBS, and cell viability was examined using trypan blue exclusion. On average, 6000–8000 single cells per tumor sample were loaded for single-cell library preparation using Chromium Single Cell Gene Expression Solution V2 (10 ⅹ Genomics).

The above method was named as “fresh” in the main text. Meanwhile, the “thawed tissue” method referred to single-cell library preparation from cells dissociated from cryo-preserved tissues. In our study, a combination of these methods was used, and the details are summarized in Table S1.

2.3 Immunohistochemistry and histology

All twenty-six samples were confirmed to be glioma by a clinical pathologist. Hematoxylin and eosin (HE) stain showed a variety of cellularity and cytonuclear pleomorphism. Ki67 immunostaining was performed for evaluating proliferation. IDH1 R132H staining was tested for the presence of the R132H mutation of IDH1. Direct sequencing of IDH (IDH1 and IDH2) was performed as described previously. The sequencing results of IDH status confirmed the results of IDH1 R132H staining. 1p/19q codeletion status was determined in Huashan Hospital either by FISH or next-generation sequencing [15].

Immunofluorescent staining of PDGFRA, TMSB4X, LGALS1, and HLA-A on RG-like and pri-OPC-like tumors was performed as previously described [16]. Primary antibodies included rabbit anti-PDGFRA (SAB Biotech, 38364, 1:200), rabbit anti-TMSBX4 (Proteintech, 19850-1-AP, 1:400), rabbit anti-LGALS1 (Proteintech, 11858-1-AP, 1:250), and mouse anti-HLA-A (Proteintech, 66013-1-Ig, 1:500). Secondary antibodies were CoraLite 488-conjugated goat anti-rabbit IgG (H + L) (Proteintech, SA00013-2, 1:1000) and CoraLite 647-conjugated AffiniPure goat anti-mouse IgG (H + L) (Proteintech, SA00014-10, 1:500). Images were taken under a confocal microscope (FV-1000, Olympus) and processed by ImageJ.

2.4 Cell culture and virus infection

Human Jurkat T cell line (clone E6-1) was purchased from China’s national collection of authenticated cell cultures and cultured at 37 °C in 5% CO₂ in RPMI 1640 medium (Invitrogen, 11875-093) containing 10% fetal bovine serum (GIBCO, 10099141C), 1% L-glutamine (Invitrogen, 35050061), 1% sodium pyruvate (Invitrogen, 11360070), and 1% penicillin-streptomycin (Invitrogen, 15140-122). Human IGKC (genebank AAA58989.1) was constructed into a lentivirus GV492 plasmid (Ubi-MCS-3FLAG-CBh-gcGFP-IRES-puromycin). For IGKC-expressing virus infection, the Jurkat cells were seeded at 1×10⁵ cells/mL in six-well plates, infected with IGKC-expressing virus (20 M.O.I, multiplicity of infection) for 16 h, washed with culture medium, and cultured for another 32 h before infection efficiency observation by EVOS M7000 cell imaging system (Invitrogen) and subsequent RT qPCR.

2.5 RT qPCR

Total RNA was extracted from the infected Jurkat cells using RNAprep Pure Tissue Kit (TIANGEN, DP431). In brief, 1 μg of RNA was reverse-transcribed into cDNA using Evo M-MLV Kit (Agbio, AG11705). Quantitative RT-PCR was performed with SYBR Green Premix Pro Taq HS qPCR Kit (Agbio, AG11701) using Applied Biosystems StepOnePlus Real-Time PCR System. The PCR cycling conditions included an initial denaturation step at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. The sequences of the primers for RT-qPCR were summarized in SourceData.

2.6 Single-cell data processing

The initial data set contained 224 667 cells and 45 068 genes. The 10 ⅹ Genomics single-cell RNA-seq data-processing pipeline was used to count and aggregate cells from all the tumor samples. Mitochondrial genes and the genes expressed in less than 20 cells were removed. The cells expressing less than 1000 unique genes or more that 5000 unique genes were excluded. The cells with UMI mapping to mitochondrial genes of more than 10% were removed. After filtering out these genes and cells, a total of 25 473 genes across 160 288 cells were obtained. The data were log-normalized for further analysis. Seurat package was used for quality control.

2.7 Ambient RNA removal and doublet identification

SoupX [17] was used for the removal of cell free mRNA contamination, and DoubletFinder [18] was applied for identifying doublets. The doublet rate was calculated based on the number of recovered cells as per Chromium 10 ⅹ user guide. Other quality control steps were performed after removing doublets and ambient RNA.

2.8 Cell type identification

Tumor and stromal cells were identified separately for each sample. The cells were clustered by unsupervised clustering based on Louvain algorithm using Seurat’s FindClusters and visualized with UMAP. Cluster-specific markers were generated using Seurat’s FindAllMarkers. The cells in the clusters were annotated according to the marker gene expression. Copy number variations from scRNA-seq data were inferred using infercnv and used to validate the tumor cells. Finally, Seurat’s SCTRansform was used to integrate all samples for downstream analyses.

2.9 Matching cancer cells to developmental glial cell types

CIBERSORTx (version 10) [19] was used for the deconvolution of scRNA-seq glioma data against the scRNA-seq data of mouse developmental cell types identified by Weng et al. [9], who generated two scRNA data sets representing the developmental stages of astroglial (hGFAP⁺) and oligodendroglial (PDGFRα⁺) lineages. Glial cell types in these lineages were subjected to scRNA profiling. Gene signature matrices for cell types in the astroglial (hGFAP⁺) and oligodendroglial (PDGFRα⁺) lineages were generated separately using CIBERSORTx to avoid batch effects. Before generating gene signature matrices, all mouse genes were converted to human genes based on Ensembl annotation (Ensembl Genes 97) and the housekeeping genes were removed following the study of Eisenberg et al. [20]. Given that the genes differ in the hGFAP⁺ and PDGFRα⁺ data sets, only common genes between these two data sets were used for generating gene signature matrices. Average single-cell gene expression levels of the malignant cells in a tumor were used as an input for CIBERSORTx to deconvolute the fraction of the cell types. The housekeeping genes were also removed from these samples. Cell fractions of glial cell types were imputed by CIBERSORTx with the ABSOLUTE mode and 1000 permutations. Given that some glial cell types appeared in both gene signature matrices, the absolute cell fraction values of each tumor were derived from CIBERSORTx. For example, for a given tumor sample, the cell fractions of the cycling neuroblast determined by both gene signature matrices were similar. Thus, its cell fractions derived from the two gene signature matrices were averaged to report the cell fraction for this tumor. Only OPC was exceptional: the cell fractions of a tumor for hGFAP⁺ OPC and PDGFRα⁺ OPC were always different. Thus, these cell fractions were separately inputted in the final report. Finally, the relative fractions of each cell type were generated based on their absolute values for each sample. The gene signature matrices were also appended into a single-gene signature matrix for deconvolution using CIBERSORTx, and similar results were obtained. The findings were validated by deconvoluting the scRNA-seq of 25 primary glioblastoma from Neftel et al. [13].

CIBERSORTx was used to extract the common gene expression features of cancer cells of each tumor, and scmap [21] was applied to map the cancer cells to the developmental cell types cell by cell to further compare the subpopulations of cancer cells within a tumor. Our scRNA seq data of glioma were integrated with the abovementioned two mouse reference data sets separately using Seurat’s IntegrateData. In addition, scmap [21] was used to map the cancer cells onto the reference data sets. An index for the reference mouse data sets was created with indexCluster, and scmapCluster was then used to obtain the assigned cell types. Prior to scmap analysis, the following modifications were conducted: (1) similar cell types from two reference data sets were merged to reduce the number of cell types and increase the number of cells in those cell types, (2) the gene list was refined to have a good separation between cell types, and (3) the threshold for cell type assignments was removed because the reference developmental cell types have a higher degree of similarity than other mature cell types. In summary, scmap [21] was used to map single cancer cells from our and Broad cohorts [13] to the mouse [9] and human [12] developmental glial cell types.

2.10 Gene signature correlation analysis in cancer cells

Most of the cancer cells matched two gene signatures of cell types after CIBERSORTx deconvolution. Correlation analysis was conducted to determine whether the two gene signatures are present in a same cancer cell. For each cell type, its genes were ranked in the hGFAP⁺ gene signature matrix based on gene expression values, and the top 100 genes were taken as a gene signature for the cell type. Pearson correlation analysis between gene expression profiles of a cell type’s gene signature and the cancer cell was performed to quantify the presence of a gene signature in a cancer cell. Pearson correlation was calculated for each malignant cell in a tumor for the two dominant gene signatures identified during deconvolution. The cells with significant correlations (P ≤ 0.05) for at least one gene signature and Pearson correlation coefficient > 0.0 for both gene signatures were collected. Finally, the Pearson correlation coefficients of each from both gene signatures were plotted against each other.

2.11 RNA velocity analysis

Velocyto [22] was applied to our single tumor cells for RNA velocity analysis, and velocyto.py was used to calculate spliced and unspliced counts. Unsupervised clustering was then performed on the tumor cells using Louvain algorithm, and each cluster was annotated based on the scmap results for the mouse developmental cell types [9]. Seurat was used to generate the UMAP embeddings, and velocyto.R for RNA velocity estimation was applied to project the results on embedded graphs.

2.12 Identification of early and late subclones

Early and late clones were identified based on copy number variation (CNV): cells having similar copy number change belong to the same subclone. Somatic CNVs are acquired over time. Hence, late subclones most likely have more CNV events than early subclones. Moreover, the presence of any deletion event in one of the subclones indicates that it is most likely a late subclone. On this basis, early and late subclones were manually defined (see SourceData and supplementary figures).

2.13 HLA-based subclone identification

For each sample, the tumor cells were separated into four HLA-based subclones named extremely low (EL), low (L), slightly low (SL), and high (H) based on the expression of HLA genes. For each HLA-A/-B/-C gene, normalized gene expression (GE) was converted into percentage gene expression:

{\mathrm{G}\mathrm{E}}_{\mathrm{C}\mathrm{e}\mathrm{l}\mathrm{l}|HLA\text{-}X}=({\mathrm{G}\mathrm{E}}_{\mathrm{C}\mathrm{e}\mathrm{l}\mathrm{l}|HLA\text{-}X}/{\mathrm{G}\mathrm{E}}_{\mathrm{M}\mathrm{A}\mathrm{X}|HLA\text{-}X})\times 100,

where GE_MAX|HLA-X = 99.5th percentile of GE_HLA-X and HLA-X refers to HLA-A/-B/-C.

For each gene, every cancer cell was assigned to EL, L, SL, or H group based on the percentage of HLA expression of 0%–10%, 10%–50%, 50%–80%, or 80%–100%. Depending on the subclone group, each cell was assigned to one of the 10 clusters, that is, cluster 1/2/3 as all three genes/at least two genes/at least one gene in EL group, cluster 4/5/6 as all three genes/at least two genes/at least one gene in L group, cluster 7/8/9 as all three genes/at least two genes/at least one gene in SL group, and cluster 10 as all three genes in H group.

2.14 Assignment of pri-OPC-like and RG-like tumors in TCGA

TCGA glioma (GBMLGG) RNA-seq data set was downloaded from Broad Institute Firehose, and associated clinical data were downloaded. Tumor purity data for glioma samples were downloaded from PanCanAtlas. The bulk RNA-seq data from TCGA contained many gene expression signals from normal cells besides cancer cells. We found that considering HLA gene expression and tumor purity gave reasonable results in classifying RG-like and pri-OPC-like tumors. Tumors with purity less than 0.6 were removed from the analysis. The tumor samples were further divided into subgroups based on their purity (i.e., 0.8–0.85, 0.85–0.9, 0.9–0.95, and 0.95–1.0). For each subgroup, averaged HLA-A/-B/-C expression was calculated, and samples were annotated with high/low expression for each gene. If a sample had a high expression of at least two genes, then it was considered as an RG-like tumor; otherwise, it was considered as a pri-OPC-like tumor.

For the analysis of the expression level of oligo-lineage markers PCDH15/PDGFRA/OLIG2/NKX2.2 and RG marker VIM, only tumors with a minimum purity of 80% were considered. Purity groups 0.8–0.85 and 0.85–0.9 were analyzed, but purity groups 0.9–0.95 and 0.95–1.0 were excluded due to their small sample size.

2.15 Mutation frequency

The mutation data of TCGA glioma were extracted from the mutation data annotated from TCGA. Only non-silent mutations were considered for further analyses. Mutation frequencies were calculated for pri-OPC-like and RG-like tumors.

2.16 Network analysis of TCGA glioma samples

TCGA glioma (GBMLGG) RNA-seq data set was downloaded from the Broad Institute Firehose. In this cohort, RG-like (n = 79) and pri-OPC-like (n = 56) tumors were identified according to the gene expression of signature genes (method: assignment of pri-OPC-like and RG-like tumors in TCGA). Network analysis was performed with Multiscale Embedded Gene Co-expression Network Analysis (MEGENA), which identified 566 and 488 significant modules for RG-like and pri-OPC-like tumors, respectively. For each module, gene enrichment analysis was performed with Enrichr. For the RG-like samples, 18 modules were enriched with pathways associated with EGFR, and 2 modules were enriched with pathways associated with LRP1B. For the pri-OPC-like samples, 35 modules were enriched with pathways associated with TP53.

2.17 Survival analysis

Survival analysis was performed using Kaplan–Meier analysis. R packages “survival” (version 2.43.3) and “survminer” (version 0.4.3) were used for estimating and generating survival curves. For survival analysis related to IGKC expression, RG-like tumors were first classified from the TCGA data set, immune cell fractions were obtained using CIBERSORTx, and the samples were sorted by CD8⁺ T cell fraction. After the bottom 20% of the samples were removed, the rest of the samples were further sorted based on IGKC expression. Survival analysis was performed between the top 20% (IGKC high) and bottom 20% (IGKC low) samples. Top/bottom 30% and 50% samples were also examined and showed similar results.

2.18 Immune cell type identification

For the high-confidence annotation of immune cells, the cells with at least one of the key immune gene markers CD45, IL7R, CCR7, CD79A, and KLRB1 were first selected. All immune cells from pri-OPC-like and RG-like samples were integrated using Seurat’s Fast Integration with reciprocal PCA (RPCA). Unsupervised clustering was performed based on Louvain algorithm and plotted with UMAP. Further immune cell type annotations were conducted based on known immune cell markers.

3 Results

To revisit the gene expression similarity of cancer cells in all three glioma subtypes to the developmental cell types, we performed single-cell RNA sequencing on 26 gliomas (i.e., 11, 10, and 5 patients of G1, G2, and G3, respectively). Papain-based digestion [14] was used to dissociate cells, and single-cell 3′ transcriptomes were generated using the droplet-based method (10 ⅹ Genomics, 8986 cells/sample, a total of 224 667 cells before filtering, Fig. S1). Clinical information and representative histology are shown in Table S1 and Fig. S2. A combination of “fresh”, “thawed cells”, and “thawed tissue” methods was used mainly due to the waiting period of pathological reports used for grouping. Pearson’s correlation analysis and linear regression model indicated that the “fresh” and “thaw” methods were comparable (Fig. S3, R² = 0.982 or 0.974). The number of genes per cell (1000–5000), number of transcripts per cell, and percentage of reads mapped to mitochondrial genes (< 10%) were used for quality filtering, and 160 288 single cells were retained (Fig. S4 and Table S2).

3.1 Gliomas could be classified into RG-like and pri-OPC-like subgroups

A recent single-cell study of glial development revealed cell types and their relations in the glial lineage hierarchy [9]. RGs are a heterogeneous population of glial- and neuronal-restricted progenitor cells. Embryonic pluripotent RG and adult neural stem cell may be clonally linked and display stem cell features [23]. RG cells can differentiate into neuroblasts, astrocytes, or pri-OPCs (Fig.1). By conducting a single-cell analysis, Weng et al. [9] showed that pri-OPC-like cells are the only prominent proliferative cell type in all three glioma subtypes. Considering that G1–G3 are distinct genomically and clinically, we hypothesized that they do not belong to one glial cell type. We tested this hypothesis using CIBERSORTx [19], which deconvolutes the global and common gene expression features of all cancer cells within a tumor into cell-type gene expression signatures and has been applied for a similar purpose [24]. We found that the 26 tumors can be mainly classified into RG-like and pri-OPC-like subgroups (Fig.1). The RG-like tumors showed a high fraction of cells transcriptomically similar to RGs and pri-OPCs (Fig.1, blue rectangle), and most of the cells of the pri-OPC-like tumors were similar to pri-OPCs and hGFAP⁺ OPCs (Fig.1, black rectangle). CIBERSORTx extracted the common gene expression features of the cancer cells of each tumor for comparison; however, it could not compare cancer subpopulations within a tumor. Thus, we further investigated cancer subpopulations using scmap [21] (Fig.1), one of the best tools for mapping single cells into known cell types according to a recent study that has benchmarked 17 methods [25]. Three thresholds were set to achieve a consensus assignment of cell types. We found that in the RG-like tumors, most of the cancer cells were transcriptomically similar to RG cells. Meanwhile, most of the cancer cells in the pri-OPC-like tumors were similar to pri-OPCs and OPCs (Fig.1).

To further validate the classification of gliomas into RG-like and pri-OPC-like tumors, we examined PCDH15 and PDGFRA expression in all cancer cells (Fig.1) and cancer subpopulations (Fig.1) in each tumor. PCDH15 and PDGFRA are OPC marker genes [6]. The pri-OPC-like tumors showed high fractions of cancer cells expressing PCDH15 or PDGFRA (Fig.1), and the dominant subpopulations of these cancer cells were pri-OPC-like and hGFAP⁺ OPC-like subpopulations (Fig.1). Thus, PCDH15 and PDGFRA are not only OPC markers to distinguish RG and OPC but also markers to distinguish RG-like and pri-OPC-like tumors. Similar findings were obtained in the Broad cohort (Fig. S5).

We then conducted the same analyses using other oligo-lineage markers (OLIG2/NKX2.2) and RG markers (VIM/SLC1A3) in all cancer cells in the pri-OPC-like and RG-like tumors from our and Broad cohorts (Fig. S6). As expected, the pri-OPC-like tumors showed significantly high fractions of cancer cells expressing oligo-lineage markers, OLIG2/NKX2.2, and the RG-like tumors had significantly high fractions of cancer cells expressing RG markers VIM (P = 0.0026, 0.0148 in our and Broad cohorts, respectively) and SLC1A3 (P = 0.0013, 0.0991 in our and Broad cohort, respectively). We further examined the mean expression level of RG marker genes (PAX6/SLC1A3/GFAP/VIM/HES1/NES/SOX2) between the RG-like and pri-OPC-like tumors from our and Broad cohort, respectively, and found that VIM was the most suitable RG marker to distinguish the two tumor groups (Fig. S6).

In addition to gliogenesis, neurogenesis might play a role in gliomagenesis [11]. To further validate these findings, we extended the same analyses by mapping cancer cells to human fetal brain cell types during neurogenesis and gliogenesis [12] using CIBERSORTx (Fig.2) and scmap (Fig.2) and obtained similar results. We found that most of the RG-like tumors were transcriptionally similar to RG and astrocytes, and most of the pri-OPC-like tumors were matched to OPCs (Fig.2 and 2B). A high fraction of cancer subpopulations expressing PCDH15 and PDGFRA was found in the pri-OPC-like tumors, especially in the OPC-like subpopulation, but not in the RG-like tumors (Fig.2 and 2D). Similar findings were obtained for the Broad cohort (Fig. S7).

Finally, the classification of pri-OPC-like and RG-like tumors was validated by examining the expression level of oligo-lineage markers (PCDH15/PDGFRA/OLIG2/NKX2.2) and RG marker (VIM) between the two subgroups from TCGA gliomas (see Methods for the assignment of the two subgroups for TCGA gliomas). Consistent with the previous findings, the pri-OPC-like tumors showed a high expression of oligo-lineage markers PCDH15/PDGFRA/OLIG2/NKX2.2, and the RG-like tumors had a high expression of RG marker VIM (Fig.2).

IHC validation was performed on PDGFRA, one of the markers discriminating RG-like and pri-OPC-like tumors. Consistent with the previous findings, the pri-OPC-like tumors showed a high fraction of PDGFRA⁺ nuclei per total nuclei (Fig.2, n = 3 pairs, three random fields were counted for each tumor).

To compare with other classification methods that are similarly based on the cellular architecture of tumors, such as the four cellular states of GBM [13], we mapped the pri-OPC-like and RG-like tumors from our cohort to the four cellular states of GBM, that is, astrocyte (AC)-like, mesenchymal cell (MES)-like, neural progenitor cell (NPC)-like, and OPC-like. Each tumor showed a hybrid combination of three or four cellular states. The majority of the pri-OPC-like tumors contained a high fraction of OPC-like and NPC-1-like states, and most of the RG-like tumors showed a high ratio of MES-1-like state. Both tumor groups contained a varied ratio of AC-like state (Fig.2). The NPC-1-like state may reflect the potential of differentiation from NPCs toward OPCs [13], and the MES-1-like state is mainly hypoxia independent [13].

In summary, these results supported that gliomas could be classified into two subgroups: RG-like and pri-OPC-like tumors.

We also noticed that most of the cancer cells in the RG-like tumors were similar to RG and pri-OPCs, and those in the pri-OPC-like ones were similar to pri-OPCs and hGFAP⁺ OPCs (Fig.1). To check whether both signatures existed on the same cell, we performed correlation analysis (Fig. S8). For the pri-OPC-like tumors, the cancer cell gene expression profile was significantly positively correlated with the gene signatures of hGFAP⁺ OPC and pri-OPC cells, suggesting a cellular state of transition from pri-OPC to hGFAP⁺ OPC (Fig. S8). To further predict potential cellular differentiation trajectories, we used RNA velocity [22] to compute for each cancer cell of each pri-OPC-like or RG-like tumor. The cancer cells mapped into the same cell type (color) were aggregated together. The directional flows often pointed from a less differentiated glial cell type to a more differentiated lineage, such as the flow from pri-OPC to OPC for a pri-OPC-like tumor (Fig.1) and the flow from RG to pri-OPC for an RG-like tumor (Fig.1). These results suggested that pri-OPC-like tumors are probably associated with a potential differentiation trajectory of pri-OPC toward OPC, and RG-like tumors may be linked with a putative developmental trajectory of RG to pri-OPC.

Finally, we studied the relationships between glioma molecular subtypes and RG-like and pri-OPC-like subgroups. We found that 78.5% of the IDH mutant samples (G2 and G3) were pri-OPC-like tumors (Fig.1 and 1B). However, the IDH WT samples (G1) had mainly two types: pri-OPC-like or RG-like, with about a 1/2:1/2 ratio (Fig.1 and 1B). Similar results were obtained for the Broad cohort (Fig.1 and Fig. S5A).

3.2 RG-like and pri-OPC-like gliomas have distinct HLA-defined subclones

To understand tumor heterogeneity, we further examined their CNVs, which can be inferred from the single-cell RNA-seq data [26]. One key application of single-cell RNA-seq is to identify CNVs. In cancer, CNVs are genetic fingerprints that can be used to infer the clonal hierarchy of a tumor [27]. In this study, we inferred the early and late subclones in each tumor (see Methods) (Fig. S9). Comparative analysis of gene expression profiles between early and late subclones revealed that the expression of MHC class I (HLA-A/B/C/E/F) and class II (HLA-DM/DO/DP/DQ/DR) was dominantly different between the early and late subclones in the RG-like tumors but not as quite in the pri-OPC-like tumors (Fig.3). In particular, the RG-like gliomas often had a higher expression of MHC class I in early subclones than in late subclones (3 out of 4 compared to 4 out of 11, red). Moreover, only RG-like gliomas had a higher expression of MHC class II in early subclones than in late subclones (2 in 4 compared to 0 in 11, red). These results suggested that during tumor evolution, HLA expression dynamics differ between RG-like and pri-OPC-like gliomas.

To further study subclonal evolution in details, we classified cancer cells into HLA-based subclones based on gene expression levels of HLA-A, -B, and -C. Each gene expression level was normalized and scaled between 0 and 100% (see Methods). We then defined four groups: extremely low (EL), low (L), slightly low (SL), and high (H) groups representing HLA-A, B, or C expression level of 0%–10%, 10%–50%, 50%–80%, and 80%–100%, respectively. A cancer cell was assigned into EL, L, SL, or H group, if any of the three HLA genes has an expression level of EL, L, SL, or H, respectively. As shown in Fig.3, most of the pri-OPC-like cancer cells belonged to the EL group, suggesting that most of them have at least one HLA-A/B/C with extremely low expression. By contrast, most of the RG-like cancer cells belonged to the L and SL groups. These results suggested that pri-OPC-like and RG-like tumors might have distinct underlying mechanisms of subclonal evolution and immune escaping programs.

This finding was statistically robust and significant because the cell number analyzed was between 1000 and 10 000 per patient (Fig.3). The differential reduction in HLA expression between RG-like and pri-OPC-like gliomas was further validated in the Broad cohort (Fig.3, 3E and S10). Pri-OPC-like gliomas had more cells in their HLA EL group, so their cancer cells had significantly reduced HLA-A/B/C (Fig.3). This finding was also statistically robust because the cell number was between 100 and 400 (Fig.3).

3.3 RG-like and pri-OPC-like gliomas have distinct immune escape programs

A reduced HLA expression is one strategy for tumor cells to escape immune surveillance. To reveal the underlying mechanisms of subclonal evolution and immune escaping programs for RG-like gliomas, we conducted pairwise comparative gene expression analyses between RG-like and pri-OPC-like tumors in HLA-based subclonal cancer cells (i.e., L and SL groups) and normal developmental cells [9]. Enrichment analyses of the modulated genes using Enrichr [28,29] showed that the cancer-driving and immune escaping programs were largely different between the RG-like and pri-OPC-like tumors. “Ras signal transduction”, “non-canonical Wnt signaling”, “FGFR signaling”, and “DNA damage response by p53” were more enriched in pri-OPC-like gliomas than in RG-like gliomas when HLA expression was reduced (Fig.4). By contrast, RG-like gliomas had more enriched terms in “regulation of NF-κB”, “MAPK cascade”, “IFN/IFN I signaling”, and “antigen processing and presentation via MHC I/MHC II” compared with pri-OPC-like gliomas (Fig.4).

During HLA reduction from SL to L and then to EL status, the modulated signaling pathways differed between the early transition (SL to L) and late transition (L to EL) status. To investigate such dynamics, we applied the same strategy as mentioned above. Differential expression analyses of L versus SL and EL versus L within pri-OPC-like or RG-like gliomas were performed, followed by GO-BP analysis. In pri-OPC-like glioma, “Wnt signaling” and “DNA damage response by p53” were upregulated first, followed by “positive regulation of Wnt”, “positive regulation of NF-κB”, “FGF signaling”, and “regulation of IFN I” (Fig.4). By contrast, “MAPK cascade” was upregulated during the late stage of HLA expression loss in RG-like glioma (Fig.4).

Given that more than 50% and 20%–30% of the tumor cells were assigned into the EL group in pri-OPC-like and RG-like glioma, respectively (Fig.3), we performed differential expression analysis between the two types for the EL group and subsequent GO-BP analysis. When the HLA expression was extremely low, pri-OPC-like gliomas were enriched in “Wnt signaling”, “regulation of canonical Wnt”, “Ras signal transduction”, “FGFR signaling”, and “Fc-ϵ receptor signaling”. Meanwhile, RG-like gliomas showed an enrichment in “regulation of NF-κB”, “cellular response to IFNγ”, “antigen processing and presentation via MHC II”, and “cellular response to IL12”(Fig.4). For the pri-OPC-like and RG-like cancer cells assigned to EL group, we further classified them into three subclones (C1 to C3) based on the expression loss of one, two and three HLA genes, respectively. Differential expression and GO-BP analyses were performed to understand the dynamic signaling changes during the gene expression loss from one to three HLA genes in pri-OPC-like or RG-like gliomas, respectively. In pri-OPC-like gliomas, “Wnt signaling”, “regulation of IFN I production”, “positive regulation of NF-κB”, and “FGFR signaling” were enriched when the loss occurred in one or two genes, followed by the “DNA damage response by p53” when the loss occurred for three genes (Fig.4). By contrast, in RG-like gliomas, “MAPK cascade”, “positive regulation of NF-κB”, and “cellular response to FGF stimulus” were enriched during early HLA loss (Fig.4).

We selected one of the pathways, NF-κB signaling, from the comparison of RG-like versus pri-OPC-like cancer cells in EL group (Fig.4) and compared the expression of two top DEG genes, LGALS1 (log₂FC = 7.28, P = 1.36E-282) and TMSB4X (log₂FC = 15.89, P = 0), between three RG-like and five pri-OPC-like tumors (Fig.4 and 4I). The ratio of LGALS1⁺HLA-A^-/DAPI⁺ nuclei and TMSB4X⁺HLA-A^-/DAPI⁺ nuclei were significantly higher in the RG-like tumors than in the pri-OPC-like tumors, supporting the enrichment of NF-κB signaling in RG-like tumors.

Immune escaping mechanisms include adaptive and innate immune escape and other oncogenic signaling that inhibit tumor immune infiltration. Our findings shows that only pri-OPC-like tumor cells mainly utilized the loss of HLA expression as one of the adaptive immune escaping mechanisms (Fig.3). Meanwhile, IFN signaling, one of innate immune escape mechanisms, was active in both subgroups (Fig.4). Interestingly, the two subgroups were prone to utilize distinct oncogenic pathways that could contribute to tumor immune escape [30–42]. RG-like gliomas harnessed NF-κB/FGFR/MAPK, whereas pri-OPC-like glioma utilized NF-κB/FGFR/WNT/TP53 immune inhibitory oncogenic pathways to protect tumor cells from immune surveillance. The differences in the immune escape mechanisms between RG- and pri-OPC-like gliomas could provide guidance for individual immunotherapy.

3.4 Distinct driver mutations and survivals between the two subgroups

The distinct HLA-based subclonal distributions and immune escaping strategies between RG-like and pri-OPC-like gliomas encouraged us to propose that RG-like and pri-OPC-like tumors may have evolved from two types of cancer founder cells, which undergo distinct driving mutations for oncogenic transformation from normal RG-like or pri-OPC-like cells (Fig.5). Thus, we used IDH WT gliomas from TCGA data set and separated them into RG-like and pri-OPC-like gliomas mainly based on HLA-A/B/C expression and tumor purity (see Methods). pri-OPC-like gliomas had significantly more reduction on HLA-A/B/C expression as comparing to RG-like ones (Fig.5). Mutation frequency was calculated in both subgroups. By analyzing the differential mutations (i.e., Fisher’s test) for Top200 mutated genes, we found that RG-like gliomas had significantly higher mutation frequency in genes such as EGFR and LRP1B than pri-OPC-like tumors (Fig.5). These results suggested that during malignant transformation, EGFR and LRP1B are the preferential mutating-drivers for transforming RG-like cells into IDH^WT tumors (Fig.5). By contrast, pri-OPC-like gliomas consistently had a higher mutation frequency of TP53 (Fig.5), indicating that TP53 is the dominant-mutating-driver for oncogenic-transforming pri-OPC-like cells into IDH^WT tumors. IDH, TP53, and ATRX are the key mutating-drivers for transforming pri-OPC-like cells into IDH^mut tumors (Fig.5) [2]. We also performed network analysis of RG-like (n = 79) and pri-OPC-like (n = 56) tumors from TCGA with MEGENA [43] and identified 566 and 488 significant modules for RG-like and pri-OPC-like tumors, respectively. Gene enrichment analysis was then performed on each module. For RG-like tumors, we identified 18 modules enriched with EGFR-associated pathways and 2 modules enriched with LRP1B-associated pathways. For pri-OPC-like tumors, 35 modules were enriched with TP53-associated pathways (Table S3).

Considering that RG-like and pri-OPC-like glioma had distinct driving mutations and immune escape mechanisms, we suspected that their survival outcomes might differ from each other. Although patients with 1p/19q codeletion and IDH^mut gliomas (G3) have a longer overall survival than patients with non-codeletion and IDH^mut gliomas (G2) [2], no significant difference was found between pri-OPC-like G2 and G3 tumors (Fig.5, first survival plot, P = 0.21, n = 159 for G2, n = 244 for G3) and between pri-OPC-like and RG-like IDH^WT tumors (Fig.5, second survival plot, P = 0.3, n = 78 for pri-OPC-like tumors, n = 56 for RG-like tumors). Temozolimide (TMZ) is currently the only drug used clinically to treat glioma. MGMT promoter methylation can increase TMZ sensitivity and improved survival [44], which occurs in about 40% of GBM and 80% of IDH^mut LGG [45]. Thus, we investigated whether RG-like or pri-OPC-like tumors affect the survival stratification of MGMT promoter methylation status (Fig.5, last two survival plots). The MGMT promoter methylation status stratified subgroups into poor and good outcomes for pri-OPC-like, but not for RG-like GBM. This finding may provide critical guidance for further improving TMZ sensitivity. Further study with increased sample size must consolidate the value of classification of pri-OPC-like and RG-like tumors in evaluating the survival outcome of MGMT promoter methylation status.

3.5 The two subgroups greatly modulate tumor-infiltrating immune cells

Understanding tumor stroma, especially immune microenvironment, could provide clues for the mechanisms of cancer immune evasion and help improve glioma immunotherapy treatment. Glioma contains a variety of tumor-infiltrating immune cells (TICs), such as microglia, tumor-associated microglia (TAM), dendritic cells (DCs), CD8⁺ T, CD4⁺ T and B cells with a dominant immunosuppressive effect [46,47]. To study whether RG-like and pri-OPC-like glioma had different immune contextures, we integrated the TICs of all RG-like and pri-OPC-like tumors by conducting clustering analysis based on gene expression profiles, and annotated cell clusters using known gene markers of immune cells [48] (Fig.6 and S11). RG-like gliomas contained a significantly higher ratio of proliferative monocytes and macrophage, and pri-OPC-like glioma had a significantly higher ratio of T helper and mast cells. Both subgroups had similar fractions of CD4⁺ T, CD8⁺ T, microglia, and DCs (Fig.6 and 6C). B and NK cells were not included because we did not detect them in most of the samples.

OLIG1 expression was higher in 42.0%, 43.2%, 40.7%, 90.2%, 42.4%, 38.9%, 63.3%, 64.8%, 36.3%, and 63.2% of infiltrating memory CD4⁺ T, recently activated CD4⁺ T, CD8⁺ T, T helper, microglia, proinflammatory microglia, monocytes, proliferative monocytes and macrophages, DCs, and mast cells in pri-OPC-like gliomas, respectively, than in RG-like TICs (Fig.6, upper panel). OLIG1 is an oligodendrocyte lineage-specific marker [49,50] (its paralog OLIG2 is able to reprogram differentiated glioblastoma cells into stem-like cells [51]). Similarly, glioma stem cell marker PTPRZ1 [52] was highly expressed in 66.5%, 64.0%, 54.1%, and 33.3% of infiltrating CD8⁺ T, memory CD4⁺ T, microglia, and DCs in pri-OPC-like tumors but not in RG-like tumors (Fig.6–6H, right panel). OLIG1 and PTPRZ1 are highly expressed in oligodendroglial cancer cells but have not been reported to be expressed in immune cells. These results suggested that pri-OPC-like cancer cells may have molecules mediating cell–cell communications to maintain a high expression of OLIG1 and PTPRZ1, which in turn induce the expression of OLIG1 and PTPRZ1 in most of their TICs. Both markers were not highly expressed in RG-like TICs. However, B cell-specific receptor IGLC2 was highly expressed in 42.8%, 37.8%, 60.6%, 35.9%, 45.3%, 25.3%, 49.2%, 52.6%, and 50.0% of infiltrating CD4⁺ memory T, recently activated CD4⁺ T, CD8⁺ T, T helper, microglia, proinflammatory microglia, monocytes, proliferative monocytes and macrophages, and DCs in RG-like gliomas, respectively. Similarly, B cell-specific receptor IGKC was highly expressed in 50.7%, 42.9%, 73.3%, 44.9%, 57.4%, 34.5%, 53.5%, 61.3%, and 62.8% of the abovementioned cell types, respectively (Fig.6, lower panel). These results suggested that RG-like and pri-OPC-like gliomas could modulate their TICs in a different manner.

OLIG1 and PTPRZ1 are normally expressed in immune cells, and B cell-specific gene markers, IGLC2 and IGKC, are expressed in other immune cells. Thus, we hypothesized that their expression in TICs could be associated with an inhibitory effect on the immune cells and provide a means of cancer immune evasion. To test this hypothesis, we examined expression correlations among OLIG1, PTPRZ1, immune checkpoints (PD1, CTLA4, TIM3, LAG3, and TIGIT), their ligands (PD-L1, PD-L2, CD80, CD86, HMGB1, LGALS9, CEACAM1, FGL1, LGALS3, PVR, NECTIN2, NECTIN3, and NECTIN4), and cytotoxicity genes (GZMB, GZMA, PRF1, GZMK, GZMM, GZMH, and GNLY) in each pri-OPC-like TIC. The expression of OLIG1 and PTPRZ1 was positively correlated with immune checkpoint PD1 in CD8⁺ T cells and TIM3 in microglia/proinflammatory microglia, respectively (Fig.7). Consistently, OLIG1 and PTPRZ1 were positively correlated with CTLA4 ligand CD86 [53] in microglia and proinflammatory microglia; with LAG3 ligand LGALS3 [54] in microglia; with TIM3 ligand HMGB1 [55] in proliferative monocytes and macrophages; and with TIM3 ligand LGALS9 [55] in microglia (Fig.7). All these results suggested that OLIG1 and PTPRZ1 expression in pri-OPC-like glioma TICs is associated mostly with TIM3 pathway (TIM3 and its ligands) as a means of cancer immune evasion. This finding was further confirmed by the positive correlation between NK cell receptor KLRB1 and TIM3 in tumor-infiltrating CD4⁺ T cells.

After extending the same analysis of B cell specific receptors (IGLC2, IGLC3, IGKC, IGHG1, IGHG3) to RG-like TICs, we found a positive correlation of B cell receptors with TIM3 in multiple TICs, with TIGIT in CD4⁺ T cells, and with LAG3 in CD8⁺ T cells (Fig.7). Consistently, positive correlations of B cell receptor expression were observed with TIM3 ligands LGALS9 and HMGB1, CTLA4 ligand CD86, TIGIT ligand NECTIN2 [56], and LAG3 ligand LGALS3 in multiple TICs (Fig.7), suggesting that B cell receptor expression in RG-like glioma TICs was also mainly associated with TIM3 pathway but supplemented with other immune inhibitory pathways such as TIGIT and LAG3 to achieve cancer immunosuppression. Moreover, B cell receptors IGLC2 and IGKC were negatively correlated with the expression of cytotoxicity genes GZMB and GZMA, which should attenuate the cytotoxic CD8⁺ T activity in RG-like tumor. At last, we also found that a higher IGKC expression was marginally associated with poor survival in RG-like tumors using TCGA samples (Fig.7). Expression of IGLC2 and other B cell receptors were very low in TCGA tumors so that we could not correlate their expressions with patients’ survivals. In addition, because OLIG1/PTPRZ1 were expressed in both cancer cells and pri-OPC-like TICs, it was hard to dissect their contributions to either cancer cells or immune cells to patients’ survivals.

We then validated the correlation of IGKC with immune checkpoints (Fig.7) by overexpression of IGKC in a human lymphoma T cell line, Jurkat cell, and examining the upregulation of LAG3 and LGALS3 (Fig.7 and 7E). IGKC was efficiently overexpressed in Jurkat cells, and LGALS3 was significantly upregulated upon IGKC overexpression. This supports our finding that IGKC was positively correlated with LAG3 ligand LGALS3 in tumor infiltrating CD8 T of RG-like glioma (Fig.7). The immunosuppressive molecule LGALS3 was reported to be an unfavorable prognostic factor in diffusely infiltrating glioma [57], which partially explained the poor survival outcome of IGKC high expression (Fig.7).

Next, we investigated differential gene expression in TICs between RG-like and pri-OPC-like tumors. Gene expression profiles of RG-like glioma-infiltrating CD8⁺ T and CD4⁺ T cells, microglia and DC were significantly different from those of pri-OPC-like ones. Further analysis of these modulated genes suggested that they were associated with immune cell activities (Fig.6–6H). For example, expression of the leukocyte trafficking-related genes (e.g., S100A family genes) were significantly lower in pri-OPC-like TICs, indicating lower activities of their immune cells in cell-cell communications.

Furthermore, specific modulation in CD4⁺ and CD8⁺ T cells, microglia and DC was also found in both subgroups (Fig.6–6H). Significant more downregulating of TCR pathway-related genes of tumor-infiltrating CD8⁺ T cells such as IL2, GZMB, GZMH and GZMA was found in pri-OPC-like tumors, suggesting much lower activity of CD8⁺ T cells (Fig.6). Consistently, T cell clone diversity-related genes such as CD3-gamma (CD3G), CD3-delta (CD3D) and CD3-epsilon (CD3E) of tumor-infiltrating CD8⁺ T cells were much more downregulated in pri-OPC-like tumors (Fig.6). In addition, CD8⁺ T cells had significantly lower expression in T cell adhesion-related genes such as CD2 and CD52 in pri-OPC-like tumors, suggesting lower activities in cell-cell communications (Fig.6). Finally, CD8⁺ T cells showed significantly lower expression of T cell costimulatory molecules such as CD27 [58], CD96 [59], CD63 [60] and CD99 [61] in pri-OPC-like tumors (Fig.6). Similar results were obtained in CD4⁺ T cells (Fig.6). Taken together, these results suggested that pri-OPC-like tumor-infiltrating CD8⁺ and CD4⁺ T cells had much lower activity than RG-like ones. As shown in Fig.3, MHC-I genes were largely downregulated in pri-OPC-like tumors so that cancer neoantigens were weakly presented on their cancer cell surface, leading to the low activity of CD8⁺ T cells.

DCs and TAM are antigen-presenting cells that could present neoantigens to other immune cells to kill tumor cells. Interestingly, tumor-infiltrating microglia in RG-like tumors had a higher expression of MHC-II related genes such as HLA-DQ/DM/DP/DR and CD74 compared with those in pri-OPC-like gliomas, suggesting that they were more active in RG-like tumors (Fig.6). TAMs can be subclassified into M1- and M2-like states. M2-like TAMs are tumor-promoting and defined by high expression levels of CD163 and CD14. They also mediate tumor angiogenesis by producing pro-angiogenic factors, such as VEGFA [62]. RG-like tumors expressed higher levels of CD163 (log₂FC = 1.85, P = 0), CD14 (log₂FC = 0.97, P = 2.60E–191), and VEGFA (log₂FC = 0.69, P = 2.75E–104) than pri-OPC-like tumors. Moreover, 65% of tumor-infiltrating microglia expressed CD163 in RG-like tumors, but only 14.7% of tumor-infiltrating microglia expressed CD163 in pri-OPC-like tumors. All these findings suggested that compared with pri-OPC-like gliomas, RG-like tumors had more M2-like TAMs, which are the key cells that create an immunosuppressive tumor microenvironment by producing cytokines, chemokines, growth factors and triggering the release of inhibitory immune checkpoint proteins in T cells.

The tumor-infiltrating DCs in pri-OPC-like glioma had a high expression of MHC-II genes and related chaperone (HLA-DQ/DR/DM/DP and CD74) (Fig.6), suggesting that they were more active in pri-OPC-like tumors. The activated DCs with MHC-II expression are able to stimulate CD8⁺ T cell activity; however, tumor-infiltrating CD8⁺ and CD4⁺ T cells were inhibited by other mechanisms in pri-OPC-like tumors. Therefore, the anti-tumor activity of the TICs was weak in pri-OPC-like tumors. By contrast, the tumor-infiltrating DCs in RG-like tumors expressed T cell-specific marker CD99, which promotes T cell communications, blood–brain barrier, and immune cell transmigration, suggesting that these tumor-infiltrating DCs could activate immune cells in lymph nodes to migrate into RG-like tumors.

In summary, we showed that RG-like and pri-OPC-like tumors extensively modulated their TICs globally (e.g., for a wide range of immune cell types) and specifically (e.g., low expression of TCR pathway genes in specific immune cell types). TIC modulation is most likely to inhibit immune cell activity to provide another layer of immune escaping mechanism for cancer cells. Therefore, RG-like and pri-OPC-like tumors could send distinct signaling molecules into their environments to modulate tumor-infiltrating immune cells (see more details in Discussion section).

4 Discussion

Abnormal gliogenesis contributes to glioma tumorigenesis [8], and neurogenesis might also play a role in gliomagenesis [11]. A recent study using single-cell RNA sequencing suggested that all three glioma subtypes are dominantly similar to pri-OPCs [9]. These results seemed conflict, considering the wide range of prognosis and genomics of the three glioma subtypes with distinct median survivals [2]. By analyzing the single-cell RNA-seq data of twenty-six samples representing the three clinical subtypes, we found that glioma can be generally classified into pri-OPC-like and RG-like tumors. This classification was robust and was confirmed globally by CIBERSORTx and individually by scmap. In addition, the expression of oligo-lineage markers (PCDH15/PDGFRA/OLIG2/NKX2.2) and RG marker (VIM) was validated in our and Broad cohorts. In addition, the expression level of the abovementioned markers was confirmed in the two subgroups from TCGA gliomas. PCDH15 is selectively expressed in OPCs and functionally important for OPC dispersion and proliferation [6]. PDGFRA is a classic OPC marker gene. These genes distinguish between OPCs and RGs [6] and between RG-like and pri-OPC-like tumors.

Our classification of glioma into RG-like and pri-OPC-like tumors is comparable with the current four cellular states of GBM classified by Neftel et al. [13]. Each RG-like or pri-OPC-like tumor showed a combination of three or four cellular states. RG-like tumors showed a higher ratio of hypoxia-independent MES-1-like state. By contrast, pri-OPC-like tumors contained more cells of OPC-like and NPC-1-like states, the latter showing a potential to differentiate toward OPCs. This finding further supports the pri-OPC feature of the pri-OPC-like tumors we classified (Fig.2).

HLA genes encode cell-surface molecules responsible for antigen presentation and immune response modulation. In this study, we found that HLA-A/B/C downregulation is associated with the status of RG- or pri-OPC-like tumors. pri-OPC-like tumor cells showed extremely low expression of HLA class I genes, and RG-like gliomas had a slightly low expression. These findings were further validated using TCGA’s pri-OPC-like and RG-like tumors. Basing on these findings, we hypothesized that pri-OPC-like tumors could have a poorer response to such as anti-PD-1 or anti-CTLA4 treatment than RG-like tumors. In IDH mutant gliomas (e.g., most of them are pri-OPC-like) have a lower immunotherapy response rate than IDH WT gliomas due to the more prominent TIL infiltration and higher PD-L1 expression in the latter [63]. A reduced HLA expression could result in reduced tumor antigen presentation and thus aid in tumor immune evasion, which has been prevalently found across a variety of cancer types and linked with poor survival [64,65].

In addition to HLA loss, innate immune escape and cancer driving pathways that inhibit lymphocyte infiltration within tumors could also lead to tumor immune evasion. Innate immune escape is mainly related to modulating IFN-associated pathways. IFN signaling is constitutively active in glioma cells and contributes to the immune evasion of glioma [66]. Cancer-driving pathways that are immune inhibitory include but not limited to NF-κB, FGFR, MAPK, WNT, and TP53 signaling. The activation of NF-κB pathway is one mechanism that drives immune cell abortive activation [30] and is associated with decreased NK cell recruitment and increased neutrophil attraction [32], all leading to tumor immune escape. b-FGF suppresses lymphocytes, especially CD4⁺ T migration, and leads to cancer immune escape [33]. MAPK is frequently activated in cancer cells, leading to immunosuppressive cascades [34]. Wnt/beta-catenin pathway is one of the important cancer driving pathways contributing to tumor immune escape [35], and its activation is correlated with T cell exclusion [36], associated with innate immune modulation such as DCs [37], and is involved in crosstalk between cancer cells and TAMs [38]. Finally, the canonical p53 mutations are not only vital for tumor escape from apoptosis and senescence but also fuel inflammation and aid tumor immune evasion [39] via their correlation with defective MHC class I expression [40,42], suppressing cGAS-STING signaling [41], and decreasing T cell infiltration [31]. In our study, we found that RG-like and pri-OPC-like tumor cells utilized distinct immune escape mechanisms during tumor subclonal evolution. RG-like gliomas harnessed IFN pathway (innate immune escape), NF-κB/FGFR/MAPK immune inhibitory oncogenic pathways, and mild reduction of HLA class I gene expression (adaptive immune escape) to protect tumor cells from immune surveillance. Meanwhile, pri-OPC-like gliomas utilized IFN pathway, NF-κB/FGFR/WNT/TP53 immune inhibitory oncogenic pathways, and severe loss of HLA class I gene expression to achieve tumor immune escape. The difference in immune escape mechanisms between RG- and pri-OPC-like glioma could provide guidance for individual immunotherapy.

We then suspected single normal RG or pri-OPC cell undergoes distinct malignant transformation programs, that is, distinct mutating-drivers, to become cancer founder cell for subsequent subclonal evolution. In this study, we found that EGFR and LRP1B are preferential drivers for transforming RG-like cells into IDH^WT tumors. This finding on EGFR is consistent with a systematic WES study of 757 patients revealing that IDH^WT is frequently amplified with EGFR [67]. Meanwhile, the finding on LRP1B is different/novel from the reported mutational landscape of 543 IDH^WT tumors [68]. We found that pri-OPC-like IDH^WT gliomas consistently had a high mutation frequency in TP53, indicating that TP53 is the dominant-mutating-driver for oncogenic-transforming pri-OPC-like cells into IDH^WT subtypes. This finding is novel for pri-OPC-like IDH^WT gliomas, although TP53 mutations have been frequently reported for GBM-IDH^mut [69] and LGG-IDH^mut gliomas [67,70].

Given that RG- and pri-OPC-like cells undergo distinct driving mutations for malignant transformation to become cancer founder cell and then undergo distinct immune escape mechanisms during tumor subclonal evolution, their distinct survival outcomes must be determined. We found that glioma with MGMT promoter methylation survived longer than non-methylation status in pri-OPC-like, but not in RG-like GBM. This finding is clinically important because MGMT promoter methylation occurs in about 40% of GBM and 80% of IDH mutant LGG [45], which can increase response to temozolomide [44], currently the only clinical drug to treat glioma. Our finding provided critical guidance for further improving temozolomide sensitivity to specific subtypes of glioma with pri-OPC-like feature.

In addition to tumor cells, studying TICs is critical for understanding immune escape and could eventually improve immunotherapy/clinical outcome. We found that RG-like and pri-OPC-like cancer cells extensively modulated TICs globally and specifically. Their global impact included the following: (1) pri-OPC-like TICs had a consensus oligodendroglial lineage marker OLIG1 and a stem marker PTPRZ1, and RG-like TICs shared B cell markers IGLC2 and IGKC, which normal T cells or microglia are not supposed to highly express. A recent report stated that T cells in type 1 diabetes also express B cell markers, which has triggered a debate [71,72]. Glial/stem markers in pri-OPC-like TICs and B cell receptors in RG-like TICs were positively correlated with immune inhibitory T cell receptors and related ligands, suggesting that RG-like and pri-OPC-like cancer cells could send distinct signals to globally modulate a wide range of TICs to weaken immune cell function and provide a means of immune evasion for cancer cells. (2) The low expression of leukocyte trafficking-related genes (e.g., S100A family genes) in pri-OPC-like TICs indicated low activities in cell–cell communications and high expression of bone marrow T marker SEC61G and commensal T marker MIF in RG-like TICs, suggesting dedifferentiation to some degree and less cell-killing activity. Their specific impact included the following. (1) NK cell-specific receptor KLRB1 and B cell-specific receptors IGLC2/IGLC3/IGKC and IGHG1/IGHG3 were specifically expressed in the tumor-infiltrating CD8⁺ and CD4⁺ T cells, respectively, of pri-OPC-like and RG-like tumors to inhibit their activities. (2) A low activity of tumor-infiltrating CD8⁺ and CD4⁺ T cells was observed in pri-OPC-like tumors. (3) A high activity of tumor-infiltrating microglia was found in RG-like gliomas, and most of them were M2-like TAMs.

RG-like and pri-OPC-like tumors probably send distinct signaling molecules into their environments to modulate respective TICs. pri-OPC-like cancer cells showed high activation for Wnt and TP53 signaling, which could reduce T cell infiltration and immune cell recruitment, respectively [31]. This finding is consistent with the less CD8⁺ T activity found in pri-OPC-like tumors. RG-like tumor cells displayed high activity for MAPK signaling, which could affect the proper development and function of T cells [73]. Additional detailed work should be conducted to illustrate how the two tumor types respectively modulate their TICs.

5 Conclusions

In summary, our study provides new insights into a comprehensive understanding of pairing developmental cell types, tumor subtypes, immune microenvironment, tumor immune escape, cancer drivers, and clinical survival. The findings would trigger new thinking, lead to new directions in the field, and give important hints in the stratification of glioma patients that will be extremely valuable when evaluating personalized therapies.

References

[1] Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, Hawkins C, Ng HK, Pfister SM, Reifenberger G, Soffietti R, von Deimling A, Ellison DW. The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro-oncol 2021; 23(8): 1231–1251

[2] Cancer Genome Atlas Research Network. . Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N Engl J Med 2015; 372(26): 2481–2498

[3] LaMonica BE, Lui JH, Hansen DV, Kriegstein AR. Mitotic spindle orientation predicts outer radial glial cell generation in human neocortex. Nat Commun 2013; 4(1): 1665

[4] Lui JH, Hansen DV, Kriegstein AR. Development and evolution of the human neocortex. Cell 2011; 146(1): 18–36

[5] Gertz CC, Kriegstein AR. Neuronal migration dynamics in the developing ferret cortex. J Neurosci 2015; 35(42): 14307–14315

[6] Huang W, Bhaduri A, Velmeshev D, Wang S, Wang L, Rottkamp CA, Alvarez-Buylla A, Rowitch DH, Kriegstein AR. Origins and proliferative states of human oligodendrocyte precursor cells. Cell 2020; 182(3): 594–608.e11

[7] Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011; 144(5): 646–674

[8] Zong H, Parada LF, Baker SJ. Cell of origin for malignant gliomas and its implication in therapeutic development. Cold Spring Harb Perspect Biol 2015; 7(5): a020610

[9] Weng Q, Wang J, Wang J, He D, Cheng Z, Zhang F, Verma R, Xu L, Dong X, Liao Y, He X, Potter A, Zhang L, Zhao C, Xin M, Zhou Q, Aronow BJ, Blackshear PJ, Rich JN, He Q, Zhou W, Suvà ML, Waclaw RR, Potter SS, Yu G, Lu QR. Single-cell transcriptomics uncovers glial progenitor diversity and cell fate determinants during development and gliomagenesis. Cell Stem Cell 2019; 24(5): 707–723.e8

[10] Zong H, Parada LF, Baker SJ. Cell of origin for malignant gliomas and its implication in therapeutic development. Cold Spring Harb Perspect Biol 2015; 7(5): a020610

[11] Matarredona ER, Zarco N, Castro C, Guerrero-Cazares H. Editorial: neural stem cells of the subventricular zone: from neurogenesis to glioblastoma origin. Front Oncol 2021; 11: 750116

[12] Couturier CP, Ayyadhury S, Le PU, Nadaf J, Monlong J, Riva G, Allache R, Baig S, Yan X, Bourgey M, Lee C, Wang YCD, Wee Yong V, Guiot MC, Najafabadi H, Misic B, Antel J, Bourque G, Ragoussis J, Petrecca K. Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy. Nat Commun 2020; 11(1): 3406

[13] Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, Richman AR, Silverbush D, Shaw MKL, Hebert CM, Dewitt J, Gritsch S, Perez EM, Gonzalez Castro LN, Lan X, Druck N, Rodman C, Dionne D, Kaplan A, Bertalan MS, Small J, Pelton K, Becker S, Bonal D, Nguyen QD, Servis RL, Fung JM, Mylvaganam R, Mayr L, Gojo J, Haberler C, Geyeregger R, Czech T, Slavc I, Nahed BV, Curry WT, Carter BS, Wakimoto H, Brastianos PK, Batchelor TT, Stemmer-Rachamimov A, Martinez-Lage M, Frosch MP, Stamenkovic I, Riggi N, Rheinbay E, Monje M, Rozenblatt-Rosen O, Cahill DP, Patel AP, Hunter T, Verma IM, Ligon KL, Louis DN, Regev A, Bernstein BE, Tirosh I, Suvà ML. An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell 2019; 178(4): 835–849.e21

[14] Zhang Y, Sloan SA, Clarke LE, Caneda C, Plaza CA, Blumenthal PD, Vogel H, Steinberg GK, Edwards MSB, Li G, Duncan JA III, Cheshier SH, Shuer LM, Chang EF, Grant GA, Gephart MGH, Barres BA. Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse. Neuron 2016; 89(1): 37–53

[15] Eckel-Passow JE, Lachance DH, Molinaro AM, Walsh KM, Decker PA, Sicotte H, Pekmezci M, Rice T, Kosel ML, Smirnov IV, Sarkar G, Caron AA, Kollmeyer TM, Praska CE, Chada AR, Halder C, Hansen HM, McCoy LS, Bracci PM, Marshall R, Zheng S, Reis GF, Pico AR, O’Neill BP, Buckner JC, Giannini C, Huse JT, Perry A, Tihan T, Berger MS, Chang SM, Prados MD, Wiemels J, Wiencke JK, Wrensch MR, Jenkins RB. Glioma groups based on 1p/19q, IDH, and TERT promoter mutations in tumors. N Engl J Med 2015; 372(26): 2499–2508

[16] Kim GW, Li L, Gorbani M, You L, Yang XJ. Mice lacking α-tubulin acetyltransferase 1 are viable but display α-tubulin acetylation deficiency and dentate gyrus distortion. J Biol Chem 2013; 288(28): 20334–20350

[17] Young MD, Behjati S. SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data. Gigascience 2020; 9(12): giaa151

[18] McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst 2019; 8(4): 329–337.e4

[19] Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, Khodadoust MS, Esfahani MS, Luca BA, Steiner D, Diehn M, Alizadeh AA. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol 2019; 37(7): 773–782

[20] Eisenberg E, Levanon EY. Human housekeeping genes, revisited. Trends Genet 2013; 29(10): 569–574

[21] Kiselev VY, Yiu A, Hemberg M. scmap: projection of single-cell RNA-seq data across data sets. Nat Methods 2018; 15(5): 359–362

[22] La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, Lidschreiber K, Kastriti ME, Lönnerberg P, Furlan A, Fan J, Borm LE, Liu Z, van Bruggen D, Guo J, He X, Barker R, Sundström E, Castelo-Branco G, Cramer P, Adameyko I, Linnarsson S, Kharchenko PV. RNA velocity of single cells. Nature 2018; 560(7719): 494–498

[23] Malatesta P, Appolloni I, Calzolari F. Radial glia and neural stem cells. Cell Tissue Res 2008; 331(1): 165–178

[24] Vladoiu MC, El-Hamamy I, Donovan LK, Farooq H, Holgado BL, Sundaravadanam Y, Ramaswamy V, Hendrikse LD, Kumar S, Mack SC, Lee JJY, Fong V, Juraschka K, Przelicki D, Michealraj A, Skowron P, Luu B, Suzuki H, Morrissy AS, Cavalli FMG, Garzia L, Daniels C, Wu X, Qazi MA, Singh SK, Chan JA, Marra MA, Malkin D, Dirks P, Heisler L, Pugh T, Ng K, Notta F, Thompson EM, Kleinman CL, Joyner AL, Jabado N, Stein L, Taylor MD. Childhood cerebellar tumours mirror conserved fetal transcriptional programs. Nature 2019; 572(7767): 67–73

[25] Abdelaal T, Michielsen L, Cats D, Hoogduin D, Mei H, Reinders MJT, Mahfouz A. A comparison of automatic cell identification methods for single-cell RNA sequencing data. Genome Biol 2019; 20(1): 194

[26] Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, Hovestadt V, Escalante LE, Shaw MKL, Rodman C, Gillespie SM, Dionne D, Luo CC, Ravichandran H, Mylvaganam R, Mount C, Onozato ML, Nahed BV, Wakimoto H, Curry WT, Iafrate AJ, Rivera MN, Frosch MP, Golub TR, Brastianos PK, Getz G, Patel AP, Monje M, Cahill DP, Rozenblatt-Rosen O, Louis DN, Bernstein BE, Regev A, Suvà ML. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science 2017; 355(6332): eaai8478

[27] Henrichsen CN, Chaignat E, Reymond A. Copy number variants, diseases and gene expression. Hum Mol Genet 2009; 18(R1): R1–R8

[28] Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, McDermott MG, Monteiro CD, Gundersen GW, Ma’ayan A. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res 2016; 44(W1): W90–W97

[29] Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma’ayan A. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics 2013; 14(1): 128

[30] Whiteside TL. The tumor microenvironment and its role in promoting tumor growth. Oncogene 2008; 27(45): 5904–5912

[31] Spranger S, Gajewski TF. Impact of oncogenic pathways on evasion of antitumour immune responses. Nat Rev Cancer 2018; 18(3): 139–147

[32] Fusella F, Seclì L, Busso E, Krepelova A, Moiso E, Rocca S, Conti L, Annaratone L, Rubinetto C, Mello-Grand M, Singh V, Chiorino G, Silengo L, Altruda F, Turco E, Morotti A, Oliviero S, Castellano I, Cavallo F, Provero P, Tarone G, Brancaccio M. The IKK/NF-κB signaling pathway requires Morgana to drive breast cancer metastasis. Nat Commun 2017; 8(1): 1636

[33] Kitayama J, Nagawa H, Yasuhara H, Tsuno N, Kimura W, Shibata Y, Muto T. Suppressive effect of basic fibroblast growth factor on transendothelial emigration of CD4⁺ T-lymphocyte. Cancer Res 1994; 54(17): 4729–4733

[34] Yaguchi T, Sumimoto H, Kudo-Saito C, Tsukamoto N, Ueda R, Iwata-Kajihara T, Nishio H, Kawamura N, Kawakami Y. The mechanisms of cancer immunoescape and development of overcoming strategies. Int J Hematol 2011; 93(3): 294–300

[35] Spranger S, Gajewski TF. A new paradigm for tumor immune escape: β-catenin-driven immune exclusion. J Immunother Cancer 2015; 3(1): 43

[36] Luke JJ, Bao R, Sweis RF, Spranger S, Gajewski TF. WNT/β-catenin pathway activation correlates with immune exclusion across human cancers. Clin Cancer Res 2019; 25(10): 3074–3083

[37] Swafford D, Manicassamy S. Wnt signaling in dendritic cells: its role in regulation of immunity and tolerance. Discov Med 2015; 19(105): 303–310

[38] Kaler P, Augenlicht L, Klampfer L. Activating mutations in β-catenin in colon cancer cells alter their interaction with macrophages; the role of snail. PLoS One 2012; 7(9): e45462

[39] Guo G, Yu M, Xiao W, Celis E, Cui Y. Local activation of p53 in the tumor microenvironment overcomes immune suppression and enhances antitumor immunity. Cancer Res 2017; 77(9): 2292–2305

[40] Sabapathy K, Nam SY. Defective MHC class I antigen surface expression promotes cellular survival through elevated ER stress and modulation of p53 function. Cell Death Differ 2008; 15(9): 1364–1374

[41] Ghosh M, Saha S, Bettke J, Nagar R, Parrales A, Iwakuma T. Mutant p53 aids cancer cells in evading lethal innate immune responses. Cancer Discov 2021; 11(5): OF14

[42] Garancher A, Suzuki H, Haricharan S, Chau LQ, Masihi MB, Rusert JM, Norris PS, Carrette F, Romero MM, Morrissy SA, Skowron P, Cavalli FMG, Farooq H, Ramaswamy V, Jones SJM, Moore RA, Mungall AJ, Ma Y, Thiessen N, Li Y, Morcavallo A, Qi L, Kogiso M, Du Y, Baxter P, Henderson JJ, Crawford JR, Levy ML, Olson JM, Cho YJ, Deshpande AJ, Li XN, Chesler L, Marra MA, Wajant H, Becher OJ, Bradley LM, Ware CF, Taylor MD, Wechsler-Reya RJ. Tumor necrosis factor overcomes immune evasion in p53-mutant medulloblastoma. Nat Neurosci 2020; 23(7): 842–853

[43] Song WM, Zhang B. Multiscale embedded gene co-expression network analysis. PLOS Comput Biol 2015; 11(11): e1004574

[44] Hegi ME, Diserens AC, Gorlia T, Hamou MF, de Tribolet N, Weller M, Kros JM, Hainfellner JA, Mason W, Mariani L, Bromberg JEC, Hau P, Mirimanoff RO, Cairncross JG, Janzer RC, Stupp R. MGMT gene silencing and benefit from temozolomide in glioblastoma. N Engl J Med 2005; 352(10): 997–1003

[45] Cohen AL, Colman H. Glioma biology and molecular markers. Cancer Treat Res 2015; 163: 15–30

[46] Pombo Antunes AR, Scheyltjens I, Duerinck J, Neyns B, Movahedi K, Van Ginderachter JA. Understanding the glioblastoma immune microenvironment as basis for the development of new immunotherapeutic strategies. Elife 2020; 9: e52176

[47] Martinez-Lage M, Lynch TM, Bi Y, Cocito C, Way GP, Pal S, Haller J, Yan RE, Ziober A, Nguyen A, Kandpal M, O’Rourke DM, Greenfield JP, Greene CS, Davuluri RV, Dahmane N. Immune landscapes associated with different glioblastoma molecular subtypes. Acta Neuropathol Commun 2019; 7: 203

[48] Nieto P, Elosua-Bayes M, Trincado JL, Marchese D, Massoni-Badosa R, Salvany M, Henriques A, Nieto J, Aguilar-Fernández S, Mereu E, Moutinho C, Ruiz S, Lorden P, Chin VT, Kaczorowski D, Chan CL, Gallagher R, Chou A, Planas-Rigol E, Rubio-Perez C, Gut I, Piulats JM, Seoane J, Powell JE, Batlle E, Heyn H. A single-cell tumor immune atlas for precision oncology. Genome Res 2021; 31(10): 1913–1926

[49] Dai J, Bercury KK, Ahrendsen JT, Macklin WB. Olig1 function is required for oligodendrocyte differentiation in the mouse brain. J Neurosci 2015; 35(10): 4386–4402

[50] Burton A. Olig1 needed for remyelination. Lancet Neurol 2005; 4(2): 80

[51] Suvà ML, Rheinbay E, Gillespie SM, Patel AP, Wakimoto H, Rabkin SD, Riggi N, Chi AS, Cahill DP, Nahed BV, Curry WT, Martuza RL, Rivera MN, Rossetti N, Kasif S, Beik S, Kadri S, Tirosh I, Wortman I, Shalek AK, Rozenblatt-Rosen O, Regev A, Louis DN, Bernstein BE. Reconstructing and reprogramming the tumor-propagating potential of glioblastoma stem-like cells. Cell 2014; 157(3): 580–594

[52] Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT, Martuza RL, Louis DN, Rozenblatt-Rosen O, Suvà ML, Regev A, Bernstein BE. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science 2014; 344(6190): 1396–1401

[53] Sansom DM. CD28, CTLA-4 and their ligands: who does what and to whom?. Immunology 2000; 101(2): 169–177

[54] Kouo T, Huang L, Pucsek AB, Cao M, Solt S, Armstrong T, Jaffee E. Galectin-3 shapes antitumor immune responses by suppressing CD8 T cells via LAG-3 and inhibiting expansion of plasmacytoid dendritic cells. Cancer Immunol Res 2015; 3(4): 412–423

[55] Wolf Y, Anderson AC, Kuchroo VK. TIM3 comes of age as an inhibitory receptor. Nat Rev Immunol 2020; 20(3): 173–185

[56] Bottino C, Castriconi R, Pende D, Rivera P, Nanni M, Carnemolla B, Cantoni C, Grassi J, Marcenaro S, Reymond N, Vitale M, Moretta L, Lopez M, Moretta A. Identification of PVR (CD155) and Nectin-2 (CD112) as cell surface ligands for the human DNAM-1 (CD226) activating molecule. J Exp Med 2003; 198(4): 557–567

[57] Hu WM, Yang YZ, Zhang TZ, Qin CF, Li XN. LGALS3 is a poor prognostic factor in diffusely infiltrating gliomas and is closely correlated with CD163⁺ tumor-associated macrophages. Front Med (Lausanne) 2020; 7: 182

[58] Song DG, Ye Q, Poussin M, Harms GM, Figini M, Powell DJ Jr. CD27 costimulation augments the survival and antitumor activity of redirected human T cells in vivo. Blood 2012; 119(3): 696–706

[59] Chiang EY, Almeida PE, Almeida Nagata DE, Bowles KH, Du X, Chitre AS, Banta KL, Kwon Y, McKenzie B, Mittman S, Cubas R, Anderson KR, Warming S, Grogan JL. CD96 functions as a co-stimulatory receptor to enhance CD8⁺ T cell activation and effector responses. Eur J Immunol 2020; 50(6): 891–902

[60] Pfistershammer K, Majdic O, Stöckl J, Zlabinger G, Kirchberger S, Steinberger P, Knapp W. CD63 as an activation-linked T cell costimulatory element. J Immunol 2004; 173(10): 6000–6008

[61] Yoon SS, Kim HJ, Chung DH, Kim TJ. CD99 costimulation up-regulates T cell receptor-mediated activation of JNK and AP-1. Mol Cells 2004; 18(2): 186–191

[62] Zhu C, Mustafa D, Zheng P, van der Weiden M, Sacchetti A, Brandt M, Chrifi I, Tempel D, Leenen PJM, Duncker DJ, Cheng C, Kros JM. Activation of CECR1 in M2-like TAMs promotes paracrine stimulation-mediated glial tumor progression. Neuro-oncol 2017; 19(5): 648–659

[63] Berghoff AS, Kiesel B, Widhalm G, Wilhelm D, Rajky O, Kurscheid S, Kresl P, Wöhrer A, Marosi C, Hegi ME, Preusser M. Correlation of immune phenotype with IDH mutation in diffuse glioma. Neuro-oncol 2017; 19(11): 1460–1468

[64] Campoli M, Ferrone S. HLA antigen changes in malignant cells: epigenetic mechanisms and biologic significance. Oncogene 2008; 27(45): 5869–5885

[65] McGranahan N, Rosenthal R, Hiley CT, Rowan AJ, Watkins TBK, Wilson GA, Birkbak NJ, Veeriah S, Van Loo P, Herrero J, Swanton C; TRACERx Consortium. Allele-specific HLA loss and immune escape in lung cancer evolution. Cell 2017; 171(6): 1259–1271.e11

[66] Silginer M, Nagy S, Happold C, Schneider H, Weller M, Roth P. Autocrine activation of the IFN signaling pathway may promote immune escape in glioblastoma. Neuro-oncol 2017; 19(10): 1338–1349

[67] Suzuki H, Aoki K, Chiba K, Sato Y, Shiozawa Y, Shiraishi Y, Shimamura T, Niida A, Motomura K, Ohka F, Yamamoto T, Tanahashi K, Ranjit M, Wakabayashi T, Yoshizato T, Kataoka K, Yoshida K, Nagata Y, Sato-Otsubo A, Tanaka H, Sanada M, Kondo Y, Nakamura H, Mizoguchi M, Abe T, Muragaki Y, Watanabe R, Ito I, Miyano S, Natsume A, Ogawa S. Mutational landscape and clonal architecture in grade II and III gliomas. Nat Genet 2015; 47(5): 458–468

[68] Brennan CW, Verhaak RG, McKenna A, Campos B, Noushmehr H, Salama SR, Zheng S, Chakravarty D, Sanborn JZ, Berman SH, Beroukhim R, Bernard B, Wu CJ, Genovese G, Shmulevich I, Barnholtz-Sloan J, Zou L, Vegesna R, Shukla SA, Ciriello G, Yung WK, Zhang W, Sougnez C, Mikkelsen T, Aldape K, Bigner DD, Van Meir EG, Prados M, Sloan A, Black KL, Eschbacher J, Finocchiaro G, Friedman W, Andrews DW, Guha A, Iacocca M, O'Neill BP, Foltz G, Myers J, Weisenberger DJ, Penny R, Kucherlapati R, Perou CM, Hayes DN, Gibbs R, Marra M, Mills GB, Lander E, Spellman P, Wilson R, Sander C, Weinstein J, Meyerson M, Gabriel S, Laird PW, Haussler D, Getz G, Chin L; TCGA Research Network. The somatic genomic landscape of glioblastoma. Cell 2013; 155(2): 462–477

[69] Ludwig K, Kornblum HI. Molecular markers in glioma. J Neurooncol 2017; 134(3): 505–512

[70] Zhang C, Moore LM, Li X, Yung WKA, Zhang W. IDH1/2 mutations target a key hallmark of cancer by deregulating cellular metabolism in glioma. Neuro-oncol 2013; 15(9): 1114–1126

[71] Ahmed R, Omidian Z, Giwa A, Cornwell B, Majety N, Bell DR, Lee S, Zhang H, Michels A, Desiderio S, Sadegh-Nasseri S, Rabb H, Gritsch S, Suva ML, Cahan P, Zhou R, Jie C, Donner T, Hamad ARA. A public BCR present in a unique dual-receptor-expressing lymphocyte from type 1 Diabetes Patients Encodes a Potent T Cell Autoantigen. Cell 2019; 177(6): 1583–1599.e16

[72] Japp AS, Meng W, Rosenfeld AM, Perry DJ, Thirawatananond P, Bacher RL, Liu C, Gardner JS, Atkinson MA, Kaestner KH, Brusko TM, Naji A, Luning Prak ET, Betts MR. TCR⁺/BCR⁺ dual-expressing cells and their associated public BCR clonotype are not enriched in type 1 diabetes. Cell 2021; 184(3): 827–839.e14

[73] Liu Y, Shepherd EG, Nelin LD. MAPK phosphatases—regulating the immune response. Nat Rev Immunol 2007; 7(3): 202–212

RIGHTS & PERMISSIONS

The Author(s) 2023. This article is published with open access at link.springer.com and journal.hep.com.cn