Gene expression data [read counts and transcripts per million (TPM) values] and corresponding clinical information of glioma patients were downloaded from UCSC Xena browser: Toil RNA-seq recomputed project and TCGA Pan-Cancer cohort¹ (Vivian et al., 2017). Among these, 231 IDH-wt glioma samples with survival data were retrieved for further analyses. Two independent datasets (CGGA and GSE16011) were used for validation (Gravendeel et al., 2009; Bao et al., 2014). RNA sequencing data (FPKM value) of gene expression was downloaded from the Chinese Glioma Genome Atlas (CGGA) database². Then FPKM value was transformed into TPM values. Affymetrix microarray-based gene expression matrix with RMA normalization was downloaded from the Gene Expression Omnibus (GEO) database (accession: GSE16011³) (Gravendeel et al., 2009). The clinical information of these three cohorts is summarized in Supplementary Table 1.

Genes were first filtered to include only protein-coding genes, and gene expressions were omitted when TPM values were less than 0.01. Genes with null expression in more than 75% of samples were removed, followed by removal of genes with median TPM expression <0.1. The filtered TPM values were then log₂ transformed and normalized to the mean and standard deviation among the samples (z score transformation). z score transformation was also performed in microarray-based expression data across samples. Subsequent gene-expression analyses were implemented with the normalized TPM values. Somatic mutation and copy-number variation were estimated and visualized by using the TCGA OncoGrid algorithm. Tumor mutation burden (TMB) was represented by the sum of the number of non-synonymous mutations. Mutation rate, aneuploidy score, and neoantigen count of glioma samples were obtained from Thorsson’s study (Thorsson et al., 2019).

Principal component analyses was implemented by using software R package ade4 with the normalized TPM values. Unsupervised clustering analysis was performed by using ConsensusClusterPlus package and k-means method with 80% item resampling and 1,000 times repetition (Wilkerson and Hayes, 2010). The cumulative distribution function and consensus heatmap were used to determine the optimal K-value. PPI network was constructed by using the STRING online tool⁴.

RNA-seq read counts were employed for DGE analysis by using edgeR package (Robinson et al., 2010). Immport database⁵ was used for immune function annotation of differentially expressed genes. A ranked list of genes created using −log10(P) × sign(log2(fold change)) (Reimand et al., 2019) from edgeR DGE output was conducted to generate enrichment scores by using the Preranked module of GSEA software (vision 4.1) (Subramanian et al., 2005) and C2: Canonical Pathway gene sets from the MSigDB database (vision 7.1) (Liberzon et al., 2011) with 1,000 permutations. The resulting enriched gene sets, with normalized enrichment score (NES) >2 and FDR < 0.01, were visualized in an enrichment map using Cytoscape (vision 3.7.2) software as described previously (Reimand et al., 2019). Fold change to median TPM value of each gene was calculated, then log2 transformed, and used for single-sample GSEA (ssGSEA), implemented in the GSVA package (Hänzelmann et al., 2013). ssGSEA enrichment score was calculated using gene sets composed of transcripts characteristic of an established T cell signature (Spranger et al., 2015), overall lymphocyte infiltration (Calabrò et al., 2009), immune cell populations (Charoentong et al., 2017), and signaling pathways (Liberzon et al., 2011). The enrichment scores generated by ssGSEA analysis were used to represent the relative abundance of immune cell or activity of signaling in each sample. CIBERSORT algorithm and signature matrix “LM22” were also used to deconvolve immune cell fractions (Newman et al., 2019). Cytotoxic activity (Davoli et al., 2017), MHC class-I (Lauss et al., 2017), and Batf3-dendritic cell (Batf3-DC) (Spranger et al., 2017) were calculated by taking mean expressions of genes in the gene sets. All gene sets we utilized are presented in Supplementary Table 2.

Weighted gene co-expression network analysis (WGCNA) was done by using WGCNA package to find modules of highly correlated genes (Langfelder and Horvath, 2008). Gene modules associated with both the CARD cluster and lymphocyte infiltrations were selected. Hub genes in the selected modules were employed to conducting LASSO penalized Cox proportional hazards regression by using glmnet package for identifing genes with best prognostic significance. The linear combination of gene expressions weighted by the estimated regression coefficients in the multivariate model was used to calculate the risk scores of each patient.

Tissue microarrays (TMAs) were constructed by using paraffin-embedded specimens from 103 IDH-wt glioma patients who were treated at the Sun Yat-sen University Cancer Center (SYSUCC) from 2010 to 2016. Patient informed consent was obtained, and the study was approved by the Ethics Committee of the Sun Yat-sen University Cancer Center. TMAs were stained with PTX3 (Abcam, Cambridge, United Kingdom) by immunohistochemistry (IHC) method as previously described (Hu et al., 2017). IHC scoring criteria were established by combining the positive proportion (1 for 0–25%, 2 for 26–50%, 3 for 51–75%, 4 for >75%) and staining intensity (0 for no staining,1 for light yellow, 2 for yellowish brown, 3 for brown) of the stained tumor cells. The final PTX3 expression score (0–12) was achieved by positive ratio × staining intensity. The best cutoff value of PTX3 staining was determined by using ROC curves with respect to OS, and applied for developing a two-level grade system of PTX3 expression. All the samples were scored separately by two independent neuropathologists, who were blinded to the patient data. The discrepancies were resolved by consensus under a microscope for multi-viewing.

Genomic, transcriptomic, and clinical data of GBM patients underwent ICI treatment was obtained from a recently published study by Zhao et al. (2019) IDH-wt gliomas with resected tissues obtained after recurrence were included for analysis. Median risk score, PD-L1 expression, and TMB were investigated for prognostic prediction.

Normalized gene expressions and ssGSEA scores between two groups were compared by Wilcoxon tests. Categorical variables were compared by chi-squared or Fisher Exact test. Kaplan–Meier analyses with log-rank tests were performed to assess survival differences. Univariate and multivariate Cox regression were conducted to calculate hazard ratio (HR). Time-dependent ROC curve was depicted for OS prediction with survivalROC package. R (vision 3.6.2) and GraphPad Prism (version 8.0.1) software were applied for the statistical analyses. P < 0.05 was considered statistically significant.

Previous studies from other cancer types suggested that the tumor-infiltrating lymphocytes, particularly activated T lymphocytes (T cell inflamed genotype), are linked to prognosis and response to checkpoint immunotherapy (Olson and Luke, 2019). To determine the association between CARD-containing genes and T lymphocyte infiltration in glioma, we analyzed the expression profiles of 31 known CARD genes and correlated them with an established T cell-inflamed score (Spranger et al., 2015) using the ssGSEA method in these 231 IDH-wt glioma patients from the TCGA dataset. Six CARD genes, namely, PYCARD, NLRC4, CASP1, CASP4, NOD2, and CARD16, were found to be associated with the T cell inflamed metrics (Figure 1A, Spearman R > 0.4, P < 0.001). K-means clustering based on the expressions of the six CARD genes identified two main clusters (CL1 and CL2) in the IDH-wt gliomas (Figures 1B,C). PCA revealed pronounced difference in the expression portraits between the two clusters (Figure 1D). ESTIMATE algorithm (Yoshihara et al., 2013) showed elevated immune cell infiltration in CL1 (Figure 1E). We then evaluated the infiltration of lymphocytes and T cells using five scoring systems, including ESTIMATE immune score, Lymphocyte infiltration score, T cell inflamed score, CYT score, and MHC score. We demonstrated that the CL1 cluster had an enrichment of lymphocytes and T cells with increased immunological activities including cytolysis and antigen presentation (Figure 1E).

To explore the difference of gene expression pattern and activation of signaling pathways between the two clusters, we compared the transcriptional profiles by using differential gene-expression (DGE) analysis. A total of 814 and 987 transcripts were identified to be up- and downregulated in CL1, respectively (| log₂ Fold change| > 1 and FDR < 0.01, Supplementary Table 3). We found that many of the top upregulated genes in CL1 had well-known immune-related functions in the ImmPort database annotation (Figure 2A and Supplementary Table 3), indicating differences in immune modulation between the two clusters. Meanwhile, Gene set enrichment analysis (GSEA) was performed, and it revealed a prominent enrichment of genes in key immune-related processes such as T cell receptor and signaling, innate immunity, and cytolysis activity in CL1. Activation of immune-stimulating pathways (e.g., IFN-γ and NOD-like receptor signalings) and those related with immune escape (e.g., PD-1 and IL-10 signalings) were observed in CL1 (Figures 2B,C and Supplementary Table 4).

To identify hub genes that are functionally linked to CARD clusters, WGCNA was performed to determine the co-expression modules. As shown in Figure 3A, transcripts in the green module displayed strong positive correlation with CARD clusters, as well as immune cell infiltration and cytotoxic T cell function (Figure 3A). Approximately half of the transcripts in the green module were differentially expressed between the two CARD clusters (Figure 3B). The 444 overlapped genes were thus used for least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazards regression to identify genes with best prognostic contribution in IDH-wt glioma (Figure 3C). A CARD-associated gene signature containing four genes (SPP1, PTX3, ABCC3, and BST1) with prognostic relevance was obtained (Figure 3D). Their expression was highly correlated with CARD-containing genes (Supplementary Figure 1A), and PPI network revealed underlying interactions (Figure 3E). To the survival significance of the CARD-associated genes, we investigated the protein level of one of the four genes, PTX3, in 103 IDH-wt glioma samples from the SYSUCC corhort by using TMA and IHC methods (Figure 3F). We found that the high expression of PTX3 was associated with a decreased OS (Figure 3G), which is consistent with the results of the tranciptional analysis from the TCGA dataset.

Then, a CARS was calculated for each patient using a formula based on the expression levels of these four CARD-associated genes weighted by their regression coefficients in multivariate analysis:

The distribution of risk score and the corresponding expression levels of these genes are shown in Supplementary Figures 1B–D. Higher CARS was observed in CL1 as compared with CL2 (Figure 3H). We segregated IDH-wt gliomas into high- and low-risk groups with the median CARS value. We found that GBM and mesenchymal subtypes were enriched in the high-risk group. According to the immune subtype classification proposed by Thorsson et al. (2019), we found that gliomas in high-risk group were dominantly attributed to the C4 subtype whereas the low-risk group had more tumors in the C5 subtype (Supplementary Figure 2). An alluvial diagram was used to visualize the attribute changes of individual samples (Figure 3I). More importantly, patients in the high-risk group exhibited significantly worse outcomes compared with those in the low-risk group (Figure 3J). The prognostic significance of CARS was validated in two additional independent datasets from the CGGA and GSE16011 cohorts (Figure 3H). Multivarite analyses also confirmed the independent prognostic value of the CARS (Table 1). Furthermore, CARS had a higher area under the curve (AUC) in receiver operating characteristic (ROC) analysis as compared with T cell inflamed score and clinical factors such as age, histological grade, and MGMT promoter status, indicating its superiority for prognosis–prediction (Figure 3K).

To further unveil the difference of tumor microenvironment (TME) in the high- and low-risk glioma groups, ssGSEA was employed to profile the level of immune cell infiltration and immune cell function in each group. We found a significant difference in major cell components associated with anti-tumor immunity in the two groups (Figure 4A and Supplementary Figure 3). Of note, TME of the high-risk group displayed a prominent T cell infiltration at baseline. We found that biomarkers for activated and memory T cells (Figure 4C, left and middle panels), as well as elevated cytolytic activity (Figure 4D), were dominantly enriched in the high-risk group, indicating the pre-existing adaptive immune response. It has been suggested that the recruitment and priming of effector T cells rely on chemotaxis and the presence of licensed professional antigen-presenting cell populations (Sanchez-Paulete et al., 2017). We showed that significant upregulation of chemokines, including CCL5, CXCL9, and CXCL10, that are pivotal for the intratumoral trafficking of effector T cells was observed in the high-risk group (Figure 4B). In addition, gliomas in the high-risk group had an elevated expression of genes involved in APC activation and function, including activated DC (Figure 4A), MHC-I class, and Batf3-DC signature, which represent cross-presentation of antigen (Figure 4D). Infiltrating cells that enhance innate immune response such as CD56^bright natural killer (NK) cells were also enriched (Figure 4A). Collectively, these results indicated that the high-risk glioma group had high CARS and exhibited T cell inflamed features with an increased pre-existing immunological activity.

As demonstrated above, patients in the high-risk group had a decreased survival despite the presence of cytotoxic T cells. The impaired function of inflamed T cells might partly account for this paradox. We found that multiple checkpoint receptors and their ligands such as PD-1, PD-L1, CTLA-4, LAG-3, TIM-3, IDO1, and CD39, as well as other immune inhibitory molecules including STAT1, were significantly upregulated in the high-risk group (Figure 4C, right panel) and thus might lead to inhibited function of T cells. The hyporesponsive T cells failed to incite effective anti-tumor immunity to maintain a durable tumor control. Moreover, the immunosuppressive microenvironment has been linked to the oncogenesis and tumor progression. Here, CIBERSORT deconvolution method revealed that macrophages dominantly infiltrated in gliomas of both groups (Supplementary Figure 4A and Figure 4E), whereas the proportion of overall and immunosuppressive M2-type macrophage (Figure 4E and Supplementary Figure 4B) as well as the macrophage/T cell ratio (Figure 4F) was higher in the high-risk group. Similarly, enrichment of regulatory T cells (Treg; Figure 4A) and an elevated Treg/T cell ratio (Figure 4G) were also observed in the high-risk group. In addition, we found an accumulation of other immune inhibitory cellular components (e.g., immature DC, CD56^dim NK cell, MDSC; Figure 4A) and cytokines (e.g., IL10, IL11, and TGFB1; Figure 4B) in the high-risk group. Consistent findings were also observed in the CGGA cohort (Supplementary Figure 3). In addition, we evaluated the immune homeostasis in these two groups by looking at the interaction between cells with anti-tumor or pro-tumor activities. Pearson’s correlation analysis revealed that the abundance of immune cells of these two populations was positively associated with each group (Figure 4H), indicating the activation of immunosuppressive components in the TME.

Accumulating evidence suggests that cancer-cell-intrinsic properties play a critical role in shaping the tumor immune microenvironment (Wellenstein and de Visser, 2018). In our study, we found that the major oncogenic pathways, such as PI3K/Akt/mTOR, JAK/STAT3, NF-κB, and TGF-β, were upregulated in the high-risk group (Figures 5A,B and Supplementary Figure 5), and activation of the WNT/β-catenin pathway was evident in the low-risk group (Figure 5C and Supplementary Figure 5). We also analyzed the genomic alterations in both groups. Significant enrichments of mutations of PTEN, PIK3R1, RB1, and MUC16 were observed in the high-risk group (Figure 6A and Supplementary Figure 6). There was no statistical difference in overall tumor mutation, aneuploidy, and neoantigen loads between the two groups (Figure 6B).

Although several multi-gene expression signatures were developed to delineate tumor immune microenvironment and implied to be candidate biomarkers for predicting ICI treatment response in gliomas, none of them have been validated (Zhu et al., 2019; Wu et al., 2020). Here, we investigated the predictive performance of CARS by using recently published data on transcriptome as well as clinical data of patients with recurrent IDH-wt gliomas under anti-PD-1 therapy (Zhao et al., 2019). In that study, responders, who were defined as radiologically stable for at least 6 months and/or having very few to no tumor cells in surgical tissues after anti-PD-1 therapy, had a favorable prognosis compared to non-responders (Figure 7A). We showed that recurrent gliomas with a high CARS were enriched in responders of anti-PD-1 therapy (Figure 7B). Moreover, we also found that an increased CARS was associated with a prolonged OS after anti-PD-1 therapy (Figure 7C). The median OS was not reached in patients with a high CARS and 8.1 months for those with a low CARS (HR: 0.12; 95% CI: 0.02–0.88; P = 0.036). In contrast, PD-L1 expression (Figure 7D) and tumor mutation burden (Figure 7E) failed to predict the prognosis of ICI in IDH-wt gliomas. Our data suggested that CARS had the capability to identify glioma patients who are more likely to benefit from anti-PD-1 therapy.