Multi-sector biopsy-based single cell RNA sequencing data (12) of glioma samples were downloaded from the Chinese Glioma Genome Atlas (CGGA) database.¹ These data were from 6,148 cells, involving 73 regions from 14 patients. We also downloaded an independent collection of single cell RNA sequencing data (13) from the Gene Expression Omnibus (GEO) database (GSE141383) for validation, which included data from 29,415 cells from nine glioma patients.

We collected publicly available glioma gene expression data with full clinical information from two databases: the CGGA (China) (14) and The Cancer Genome Atlas (TCGA) (the United States²). To remove the impact of different tumor origins, patients diagnosed with oligodendroglioma, anaplastic oligodendroglioma, oligoastrocytomas and anaplastic oligoastrocytomas were excluded. We obtained 378 samples with clinical and molecular information from the CGGA database. A further 810 astrocytic glioma (grade II–III astrocytoma and glioblastoma) samples with clinical and molecular information (copy number, mutation, gene expression data) were downloaded from TCGA database. The details of patients are described in Table 1.

Expression profiling of samples from 28 recurrent glioblastoma patients who received neoadjuvant anti-PD-1 immunotherapy were also downloaded from the GEO database (GSE121810) to explore the predictive value of TARA infiltration in immune checkpoint inhibition.

To quantify the tissue-infiltrating immune and stromal cell populations in glioma samples, we applied the Microenvironment Cell Populations-counter (MCP-counter) algorithm, which allows for sensitive and specific discrimination of eight immune and two stromal cell populations, including CD3+ T cells, CD8+ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, cells originating from monocytes (monocytic lineage), myeloid dendritic cells, neutrophils, as well as endothelial cells and fibroblasts. MCP-counter is an objective method based on a methodological framework of transcriptome markers that outputs an abundance estimate per cell population that enables an inter-sample comparison. This procedure was performed using the MCPcounter R package (15).

Chromosome 1p/19q status was inferred from gene expression on chromosome 1p and 19q from start to end using a Gaussian window smoothing function (16).

The infiltration scores of TARAs were calculated for 1,079 astrocytic glioma patients for whom mRNA expression data (RPKM values) were available using the reactive-astrocyte representative genes determined by Heiland et al. (11). Single sample gene set enrichment analysis (ssGSEA) was used to calculate the enrichment scores (11).

Gene annotation enrichment analysis of TARA-related genes was executed using the BP function in clusterProfiler R package (17). A p value <0.05 and a q value <0.05 were used as cutoff values. Differentially expressed genes between TARA low and TARA high groups were calculated using limma in R package (adjusted p value < 0.05). The gene set enrichment analysis (GSEA) results of RAS-low and RAS-high groups were visualized by enrichplot in R package. The corrplot package was used to explore the correlation between RAS and other relevant biological processes curated by Mariathasan et al. (18).

The single-cell sequencing data were analyzed using a standard Seurat V4 pipeline (19). Associations between continuous variables were calculated using the Pearson correlation method. Differences between groups were estimated using unpaired Student’s t tests, one-way ANOVA, or chi-square tests. The Kaplan–Meier survival curve and log-rank test were applied using survminer R package to estimate the distribution of overall patient survival. The forestplot package was applied to depict the efficacy of RAS in subgroups. The timeROC package was used to generate time-dependent ROC curves to evaluate the 1 year predictive accuracy of RAS and PD-1 expression. Maftools package was used to depict the mutation landscape of glioma samples in TCGA dataset. All statistical analyses were executed using R software (version 3.6.2). A two-sided p value < 0.05 was considered statistically significant.

To assess the possible infiltration of TARAs in gliomas, we first analyzed multi-sector biopsy-based single-cell RNA sequencing data of 6,148 cells. The type of the cells were first auto-annotated by SingleR package (Supplementary Figure S1) and the copy number of the cells were analyzed by inferCNV package (Supplementary Figure S2). The astrocytes showed glioma copy number features: 1p and 19q codeletion, 7 gain and 10 loss. To avoid potential contamination of malignant cells in the following analysis, the astrocytes were excluded. The remaining cells (N = 2,286) were well clustered into 13 clusters (Figure 1A, 0 ~ 12). To annotate TARAs among these cells, we used a previously published TARA signature, and found high expression of this signature in clusters 1, 3, 6, 8, 9, and 11 (Figures 1B,C). These cells accounted for approximately 43.4% (993/2,286) of non-malignant cells and 15.7% (993/6,148) of the total number of cells.

To validate these results, we analyzed GEO data of 29,415 cells from nine glioma patients (GSE141383). Similarly, the non-malignant cells (N = 7,243) were clustered into 16 clusters (0 ~ 15, Figure 1D). Clusters 0, 2, 7, 11, and 13 (2,685/7243, 37.1%; 9.1% (2,685/29,415) of total cells) showed high expression of the TARA signature (Figures 1E,F). These results indicated that TARAs abundantly infiltrated the glioma micro-environment.

To characterize and understand the biological and clinical features of reactive astrocyte infiltration, we quantified the relative abundance of TARAs in 1079 astrocytic glioma samples using the mRNA-based gene signature determined by Heiland et al. (11). The glioma samples were ranked according to increasing TARA infiltration (Figure 2A). Analysis of TCGA dataset revealed that WHO grade IV samples (Figure 2B), mesenchymal samples (Figure 2C) and IDH wild-type GBM samples (Figure 2D) were more likely to have more TARA infiltration. Meanwhile, correlations between the age at diagnosis and gender and reactive astrocytic score were not statistically significant (Figure 2A). In terms of well-known molecular markers, IDH mutation, 1p/19q codeletion, MGMT promotor methylation, wild-type TERT promoter, and ATRX mutation glioma samples correlated with lower infiltration of reactive astrocytes. Similar results were found in the CGGA dataset (Supplementary Figures S3A–D). In addition, recurrent glioma samples had a higher reactive astrocyte score than primary glioma samples. Moreover, when WHO grade (Supplementary Figure S3E), IDH status (Supplementary Figure S3F) and 1p/19q status (Supplementary Figure S3G) were added as sub-classifiers, the reactive astrocyte scores of recurrent glioma samples were significantly higher than those of primary glioma samples in the same subset. These results indicated that TARA infiltration was tightly associated with clinical and pathological features of glioma patients and might be a new diagnostic biomarker.

To uncover genomic alterations associated with TARA infiltration, we analyzed somatic alterations and copy-number variations in TCGA dataset. Based on the increasing reactive astrocyte scores, we stratified TCGA glioma samples by defining tumors in the top 25% as the high TARA infiltration group (HRA) and the bottom 25% as the low TARA infiltration group (LRA), and we found mutated genes to occur significantly in each group. LRA gliomas had more mutations in TP53, IDH1, and ATRX (Figure 3A). Mutations in PTEN, TTN, EGFR, and NF1 occurred more frequently in HRA samples (Figure 3A). We also observed a significantly different frequency of mutations in MUC16, FLG, and PICK3A.

Glioma is characterized by increased genomic instability with extensive copy number alterations (20, 21). We therefore assessed the copy number variations between the two groups. As shown in Figure 3B, the incidence of chromosome 7, 19, and 20 amplification, and chromosome 10 and 22 deletion were increased with increasing TARA infiltration. The HRA samples possessed more copy number variations, including deletions at 9p21.3 (CDKN2A/CDKN2B), 10q23.3 (PTEN/ATAD1) and 13q14.2 (KPNA3), and amplification of 7p11.2 (EGFR) (Figure 3C). Meanwhile, the recurrent copy number losses at 9p21.3 encompassing the CDKN2A/CDKN2B locus and at 10q26.3 encompassing the MIR202 locus, and amplification at 7p11.2 encompassing the EGFR locus were the most common copy-number variations in LRA samples (Figure 3C). Globally, LRA samples were more stable and had significantly fewer somatic copy number alterations than HRA samples. These results indicated that TARAs tended to infiltrate gliomas with distinct genomic alterations.

To reveal the functional characteristics of TARA infiltration, we identified the genes that were significantly positively (Pearson coefficient > 0.5) or negatively (Pearson coefficient < −0.5) correlated with reactive astrocyte score and annotated the gene sets using R package clusterProfiler. The top 20 significant biological processes are summarized in Figures 4A,B. Genes positively correlated with TARA scores were mainly involved in neutrophil activation, leukocyte migration, response to interferon gamma, T cell activation and other immune-related activations. The negatively correlated genes were enriched in synapse organization, modulation of chemical synaptic transmission, regulation of trans-synaptic signaling, cell morphogenesis involved in neuron differentiation, and other synaptic-signal-transduction-related biological processes.

Furthermore, GSEA of differentially expressed genes ranked by the log2-transformed fold change between HRA and LRA samples revealed significantly enriched hallmarks of malignant tumors in HRA samples to include interferon gamma-mediated signaling pathway, inflammatory response, positive regulation of the JAK–STAT cascade, positive regulation of NIK/NF-kappa B signaling, and the tumor necrosis factor biosynthetic process. Meanwhile, functions, including carboxylic acid transport, detection of calcium ion, acid receptor signaling pathway, L-glutamate import and neuron cell–cell adhesion were related to the low RAS group (Figures 4C,D). Therefore, TARA infiltration might activate oncogenic pathways in glioma to promote malignant transformation.

To determine the characteristics of the immune micro-environment associated with TARA infiltration, we deconvolved the immune cells with mRNA sequencing data using the MCPcounter method. The immune score (Cor = 0.68, p < 0.0001), stromal score (Cor = 0.65, p < 0.0001), fibroblasts (Cor = 0.71, p < 0.0001), monocytic lineage (Cor = 0.37, p < 0.0001), neutrophils (Cor = 0.37, p < 0.0001), CD8⁺ T cells (Cor = 0.68, p < 0.0001), T cells (Cor = 0.68, p < 0.0001) and myeloid dendritic cells (Cor = 0.68, p < 0.0001) all showed significantly positive correlation with reactive astrocyte score. Glioma purity showed significant negative correlation with the reactive astrocyte score (Figure 5A). Meanwhile, the HRA samples were characterized by increases in the infiltration of T cells, CD8⁺ T cells, monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells and fibroblasts (Figure 5B).

We then analyzed the differential expression of checkpoint relevant genes and cytokine and chemokine transcripts (22) between HRA and LRA groups (Figures 5C,D). HRA cases were associated with high expression levels of IDO1, PD1, CTLA4, PD-L1, PD-L2, HAVCR2, and LAG3. Moreover, the expression levels of cytokine and chemokine-related transcripts (CD8A, CXCL10, CXCL9, GZMA, GZMB, PRF1, and TBX2) were all significantly upregulated in the high RAS group. We also tested known signatures within TCGA dataset to better describe the functionality of the reactive astrocyte score. The results indicated that the reactive astrocyte score was significantly associated with mismatch repair, cell cycle regulation, epithelial-mesenchymal transition, CD8+ T effector, immune checkpoint, antigen processing machinery, and TGFb response signal (Figures 5A–E). These analyses indicated that patients with more reactive astrocyte infiltration are accompanied by immune infiltration and high activation of oncogenic signaling pathways.

We stratified glioma patients into distinct subgroups based on gender, WHO grades, IDH status, 1p/19q status, MGMT promotor status, tumor stage, chemotherapy and radiotherapy, and observed the prognostic value of the reactive astrocyte score as a continuous variable in TCGA and CGGA datasets (Figures 6A,B). The forest plot indicated that the reactive astrocyte score was a prognostic predictor in nearly all subgroups, except for WHO grade II and 1p/19q codeletion groups. Meanwhile, we used dichotomization to separate cases to depict the survival curves according to the median value and found that patients with higher reactive astrocyte scores had a significantly shorter overall survival than their counterparts in nearly all distinct subgroups in both datasets (Supplementary Figures S4, S5). The PD-1 immune checkpoint inhibitor can induce lasting tumor remission in patients with diverse advanced cancers (23). We therefore evaluated the predictive value of recurrent glioblastoma patients with neoadjuvant anti-PD-1 immunotherapy (GSE121810) (24). The clinical characteristics of recurrent glioblastoma patients are shown in Figure 6C. The time-dependent AUC indicated that the reactive astrocyte score had a high accuracy in predicting 1 year survival (Figure 6D). Furthermore, the reactive astrocyte score (Figure 6E) demonstrated a stronger power than PD-1 expression (Figure 6F) in predicting recurrent glioblastoma patient overall survival with anti-PD-1 treatment. Patients with lower reactive astrocyte score were more likely to benefit from immune checkpoint treatment. Therefore, TARA infiltration might be a prognostic and a predictive marker for immune checkpoint inhibition.