RNA expression data for human gliomas were downloaded from The Cancer Genome Atlas (TCGA) (11) and the Chinese Glioma Genome Atlas (CGGA) (12). These datasets contain bulk mRNA expression profiles and corresponding clinical information of patients with glioma. A total of 702 glioma samples from TCGA and 693 glioma samples from CGGA were used in our study. At present, the WHO's classification on glioma utilized by TCGA and CGGA is the fourth WHO's classification of tumors of the central nervous system (13). RNA-seq data of 122 GBM samples with anatomic information from the Ivy Glioblastoma Atlas Project (IGAP) (14) were used in this study. RNA-seq data from 15 patients with recurrent GBM treated with adjuvant anti-PD-1 immunotherapy were accessed through GSE121810.

Survival curves were performed by Kaplan–Meier analysis between CD44 higher and lower group and were tested for significance using the Mantel–Cox log-rank test. Hazard ratio (HR) and confidence interval (CI) were also computed to confirm the prognostic value of CD44 in patients with glioma. Furthermore, Cox proportional hazards model, including patients' age and gender, was used to evaluate the predictor effect of CD44. A value of p < 0.05 was considered statistically significant.

Estimation of infiltrating stromal cells and immune cells in glioma samples was accomplished by the method ESTIMATE (15). Then, the xCell pipeline (16) was adopted to explore the relationship between CD44 expression level and glioma-infiltrating cell types. Because some cell types from xCell reference database do not distribute in the brain and glioma, we only used 41 cell states listed in Supplementary Figure 4 in this part analysis. Furthermore, we utilized CIBERSORTx portal (17) to infer immune cell compositions among TCGA and CGGA glioma using a set of 22 human immune cell reference profiles and predicting their absolute composition ratio within individual glioma sample.

Weighted gene coexpression network analysis (WGCNA), a bioinformatics algorithm (18), was used to generate unsigned coexpression networks in TCGA and CGGA glioma. WGCNA applies topological overlap measure, a robust measure of network interconnectedness and measures the connection strength between two adjacent transcripts and all other transcripts in a network, to cluster genes into network modules. The minimum size of modules was 30 transcripts and were randomly color labeled. Next, we related the gene modules to clinical traits and selected the module of interest which contained the CD44 and had the highest correlation coefficient with CD44 expression level. Intramodular connectivity of transcripts was used to identify hub genes in the module of interest.

Gene ontology (GO) (19) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (20) analyses were used to explore the biological function of gene sets of interest. Pearson's correlation coefficients between CD44 and all other genes were calculated, and CD44-associated genes were defined as genes with |r| >0.5 in TCGA glioma dataset, |r| > 0.25 in CGGA glioma dataset, and |r| > 0.5 in IGAP GBM dataset. Additionally, GSEA (Linux_4.1.0) software was used to determine whether the gene set of PD-1 signaling (http://www.gsea-msigdb.org/gsea/msigdb/cards/REACTOME_PD_1_SIGNALING.html) shows statistically significant difference between anti-PD-1 immunotherapy responder and nonresponder subgroups.

The immunomodulator genes were downloaded from the study introduced by Thorsson (21). These immunomodulators with CD44 were imported into the STRING website (https://string-db.org/) to construct protein–protein interaction network. All the parameters were default.

Raw count matrices of 15 patients with recurrent GBM treated with adjuvant anti-PD-1 immunotherapy were utilized in this part analysis. These 15 patients were divided into responder and nonresponder subgroups based on OS cutoff for 300 days. The patient with OS greater than 300 days was labeled as therapy responder, but the one with OS lower than 300 days was regarded as nonresponder. In addition, differentially expressed genes (DEGs) between nonresponder and responder groups were identified by the package DESeq2 (1.32.0) (22).

The human GBM signal cell RNA-seq data were downloaded from The Gene Expression Omnibus (GEO) (23) under accession number GSE103224. The expression matrices were analyzed using the R package Seurat (24). We removed all cells from the downstream analysis where >10% transcripts aligned to the mitochondrial genome, or which had either fewer than 400 or more than 2,500 RNA counts. Besides, gene was filtered out when it expressed in less than 10 cells among all the 7 GBM samples. The remaining 28,825 genes in 17,357 cells passed quality control and were used into further analysis. Malignant tumor cells were distinguished from stromal cells by copy number variants (CNVs), and nontumor cells were identified using known cell-type marker genes. FindMarkers in Seurat package was used to identify the DEGs between CD44+ and CD44- cell cluster. Additionally, the top 100 DEGs according to the avg_logFC were underwent the gene set enrichment analysis. Pairwise Pearson's correlation with CD44 was established in CD44+ cell cluster, and the top 100 correlated genes according to correlation coefficient were used in subsequent enrichment analysis. We took the standard workflow of pySCENIC (0.11.2) (25) to build the regulon network of different cell types. Regulons were ranked by the regulon specificity score from high to low, and the top 20 regulons of each cell type were used to plot the heatmap among different cell types. The regulon specificity scores were scaled by cell types. The Monocle 2 package (2.14.0) (26) was applied to construct single-cell pseudotime trajectories. Considering the small number of oligodendrocytes, endothelial cells, pericytes, and T cells, only malignant tumor cells and TAMs were used in trajectory analysis, respectively. The DEGs were identified by Monocle 2, and these genes filtered by qval <0.05 were further used in pseudotime order analysis.

The IHC data of unique GBM subtype samples from the IGAP and of normal brain tissues and glioma samples from the Human Protein Atlas (https://www.proteinatlas.org/) were downloaded and used in this study.

Formalin-fixed paraffin-embedded sections of glioma were collected and their pathological results were confirmed as glioma according to the 2016 WHO's classification of central nervous system tumors (13). Tumor sample use was approved by the Institutional Review Board at the Nanjing Brain Hospital Affiliated to Nanjing Medical University. Tumor sample was incubated with the first primary antibody against CD44 (1:5000 for IF, Servicebio) overnight at 4°C and then with the first matching secondary antibody at room temperature for 50 min in dark condition. Then, it was incubated with the second antibody CD163 (1:5000 for IF, Servicebio) overnight at 4°C and later with the second matching secondary antibody at room temperature for 50 min in dark condition. Slide was counterstained with DAPI for nuclei visualization. Finally, sample slide was imaged by Imaging System from Nikon. The CaseViewer software (3DHISTECH) can unmix and remove autofluorescence and analyze the multispectral images.

To answer the question whether CD44 is involved in malignant progress of glioma, we compared its expression level in different WHO grades, IDH types, 1p19q states, and recurrent status in TCGA and CGGA glioma mRNA-seq datasets. First, CD44 expression in both low-grade glioma (LGG) and glioblastoma (GBM) is higher than normal brain tissues, which means CD44 plays a role in glioma tumorigenesis (Figure 1A). It turns out that CD44 is significantly upregulated in GBM (WHO IV) than that in WHO grades II and III gliomas in TCGA glioma dataset (Figure 1B). Additionally, CGGA RNA-seq data validate well to this result (Figure 1C). Additionally, the IHC data of normal brain tissues and glioma samples reveal that high-grade glioma has the highest CD44 protein expression level, but that the normal brain tissue has the lowest CD44 expression (Supplementary Figure 2). These shows that a higher CD44 expression is associated with a higher malignancy in glioma.

Because IDH1 mutation frequency is >68% in LGG, but only 12% in GBM (27), and 1p19q codeletion is more common in LGG than GBM (28), we studied CD44 expression level in different IDH types and 1p19q states. In both the TCGA glioma and CGGA glioma, samples with strong CD44 transcription are primarily having wild-type IDH, but the ones with low CD44 expression have IDH mutation (Figures 1D,E), which can also be discovered among GBM groups (Supplementary Figures 3A,C). We also found a significant correlation between low CD44 transcription level and 1p19q codeletion status in both TCGA and CGGA glioma (Figures 1F,G, Supplementary Figures 3B,D). Glioma patients with IDH1 mutation have median OS three times longer than those without (27), codeletion of 1p and 19q is associated with better survival rates in glioma (29), and thus, low CD44 expression patients with glioma may have better clinical outcomes.

Furthermore, we uncovered up expression of CD44 after glioma recurrence in GBM, wild-type IDH, IDH mutation, and 1p19q non-codeletion state in CGGA glioma dataset (Figures 1H–J). Although multimodal treatments can prolong life, recurrence is inevitable in glioma, especially GBM, and it can obtain more malignant phenotype after recurrence (30). Herein, malignant glioma has relative higher CD44 expression level, and CD44 is involved in malignant progress of glioma.

CD44 transcription level is sufficient to predict OS of patients with glioma in TCGA and CGGA datasets (Figure 2). Compared with the lower expression of CD44 group in whole-grade glioma, patients with higher CD44 expression have significantly shorter survival. Given the obvious heterogeneity in glioma (31, 32), we further investigated the prognostic value of CD44 among different WHO grades, IDH types, 1p19q states, and recurrent status. Interestingly, we observed the same prognostic tendency of CD44 in these different compared groups (Figure 2). In addition, we did Cox-regression analysis of CD44 with patients' age and gender, namely multivariate Cox-regression analysis, in different glioma groups to evaluate their predictor values in glioma. In the Cox proportional hazards model, patients with higher CD44 expression level have higher mortality than those with lower CD44 expression in different compared groups (HR>1, p < 0.05), except in 1p19q codeletion status (p > 0.05) (Table 1). These findings suggest that CD44 could be a negative prognosticator in glioma.

First, we analyzed the distribution of CD44 in different GBM transcription subtypes to obtain an overview of the molecular expression pattern of CD44. The CD44 mRNA expression is quite different in the three different GBM characteristic subtypes (Figure 3A). Samples with high CD44 expression are mainly concentrated in the mesenchymal subtype, which is verified by IHC results (Figure 3B). Considering the mesenchymal GBM subtype is featured with more infiltrating tumor-associated macrophages (TAMs) and immunosuppression (33), we speculated that CD44 may play an important role in glioma immune environment. Next, we adopted the other glioma subtype model that portrays adult diffuse grades II, III, and IV gliomas comprehensively (34). They reveal that a subtype of IDH mutant glioma is associated with DNA hypermethylation and better outcome, and that a group of wild-type IDH diffuse glioma shows molecular similarity to pilocytic astrocytoma and relatively favorable survival. We explored the expression difference of CD44 among pan-glioma transcriptome subtype (Figure 4A), pan-glioma DNA methylation subtype (Figure 4C), and supervised pan-glioma subtype (Figure 4E). Through comparing the survival in different subtypes (Figures 4B,D,F), we discovered that the subtypes with high CD44 expression are related to bad clinical outcomes as well.

Because immune cells are the main participants in immune response, we characterized the relationship between CD44 and the infiltrated immune cells in glioma. The results indicate that CD44 is positively associated with the immune score and stromal score in TCGA and CGGA glioma databases (Figure 5A), which suggests that CD44 does work in glioma immunity. Additionally, the standard CIBERSORTx pipeline was applied to further evaluate the relationship between CD44 and infiltration of 22 immune cells in glioma microenvironment (17). In addition, the results show that CD44 has the most significantly positive correlation with M2-type macrophages (Figure 5B). The analysis results generated by xCell workflow also confirmed the positive correlation between CD44 expression level and macrophages M2 score (Supplementary Figure 4). It is well known that M2 TAMs accelerate tumor cell invasion and angiogenesis and suppress antitumor immunity (35), which explains why CD44 can be prognostic biomarker and confirms immunity function of CD44.

Then, to gain insight into the particular contribution of CD44 to glioma immunity, GO and KEGG analyses were used to screen relevant genes in TCGA and CGGA glioma databases. We adopt two different strategies, namely the WGCNA method and the Pearson's correlation method, to obtain the CD44 relevant genes. First, the general WGCNA analysis flow was taken to divide the genes in different modules and to correlate these module eigengenes with clinical traits including CD44 expression level. The modules were chosen as they had the highest correlation coefficient with CD44 transcription level and contained the CD44 as well. The brown, turquoise, and brown and blue module emerges as the most significant module for TCGA glioma, TCGA GBM, CGGA glioma, and CGGA GBM, respectively (Figures 6A,C,E,G, Supplementary Figure 5). Then, we picked out hub genes in these most significant modules and used them in subsequent enrichment analysis. These hub genes correlated with CD44 are involved in immune response (e.g., macrophage activation, cytokine production, microglial cell activation, and T cell activation) in both TCGA and CGGA glioma datasets (Figures 6B,D,F,H, Supplementary Table 1). Apart from the positively relevant genes, negatively relevant ones can also help us to understand the function of CD44. To acquire CD44 coexpressed genes in the second Pearson's correlation analysis part, we set different Pearson's |R| value cutoff in TCGA and CGGA glioma datasets as mentioned in the Methods section. Although the enrichment analysis results of these positively correlated genes are also most involved in immunity (Figure 7), we found that some of them are related to angiogenesis, cell adhesion, ECM-receptor interaction, and so on. (Supplementary Tables 2, 3), which is consistent with the former studies (36, 37). On the other hand, the negatively relevant genes of CD44 take parts in synaptic function (Supplementary Figure 6, Supplementary Tables 2, 3). To sum up, CD44 plays an important role in glioma immunity.

Glioblastoma is one of the most heterogeneous tumors, and multiple TCGA transcription subtypes can coexist in the close tumor regions (33). Spatial heterogeneity of glioma is the barrier for treatment. Thus, we compared the CD44 expression level in different GBM tumor regions, and genes were clustered according to their tumor region features to obtain the expression patterns whose spatial states are highly desirable as CD44. The result demonstrates that CD44 mainly distributes in the perinecrotic region within tumor (Figure 8A). Some of the top correlated genes are immune genes (e.g., TNC, CCL2, SOCS3, TNFRSF12A, PTX3, and VAT1) (38, 39), invasive genes (e.g., GBP2, EMP1, VIM, ANO6, RBM47, CHI3L1, HMOX1, and ICAM1) (40–42), proliferative genes (e.g., ANXA2, ZFP36L2, and CHI3L2) (43, 44), and apoptotic genes (e.g., SOD2) (45). From the spatial-based correlated gene pattern, CD44 may take part in GBM immunity, invasion, and proliferation process. Thus, we further did enrichment analysis of this spatial-based correlated gene pattern to explore underlying biological function of CD44. The results confirmed that CD44 is related to immunity, and invasion (Figures 8B,C, Supplementary Table 4). We also found that CD44 takes part in response to decreased oxygen levels, HIF-1 signaling pathway, and angiogenesis process (Figures 8B,C, Supplementary Table 4), which is consistent with the spatial distribution preference of perinecrotic zone (Figure 8A). Because HIF1A is a regulator of CD44 and increases CD44 expression, hypoxic region in breast cancer has also been reported to contain cells with a higher concentration of CD44 expression (46). In GBM, we uncovered that CD44+ cells mainly distribute in hypoxic region with high expression of immune genes.

Whereas, tumor is constituted by malignant cells and surrounding stromal cells, we know few about the heterogeneity of CD44 transcription in glioma at single-cell resolution until now. Global gene expression however is readily measurable, and the advent of single-cell RNA-seq technology enables us profile individual cell expression within tumors. Therefore, we applied single-cell RNA-seq dataset for GBM (47) to further explore CD44 expression at single-cell resolution. First, we annotated GBM cells into specific cell types according to CNVs and known cell markers (Supplementary Figures 7B,C). GBM contains malignant tumor cells and TAMs together with a small number of oligodendrocytes, endothelial cells, pericytes, and T cells (Figure 9A). Then, we compared the CD44 expression level among different cell types and found that only malignant tumor cells, TAMs, and T cells express CD44 (Figure 9B). TAMs are the major immune cells in GBM, which is in accord with our finding that CD44 acts in glioma immunity. Because only a part of TAMs express CD44 (Figure 9C), DEGs between CD44+ TAMs and CD44- TAMs were calculated and underwent further enrichment analysis. The functional analysis showed that CD44+ TAMs are enriched in leukocyte migration, leukocyte chemotaxis, positive regulation of cell-substrate adhesion, extracellular matrix organization, positive regulation of angiogenesis, and so on (Figure 9D, Supplementary Table 5). In addition, CD44+ TAMs are in M2-type polarization state presenting with high expression of CD163, CD206, and STAT3, and they can also activate angiogenesis by secreting VEGFA and promote cell invasion through secreting matrix metallopeptidase (e.g., MMP14, MMP19) (Figure 9E). Because CD44 expression level is positively correlated with M2-type macrophage infiltration (Figure 5B), we verified the infiltrating M2 TAMs among different CD44+ cell infiltration levels by IHC (Figure 10A) and found that the more CD44+ cells are detected, the more M2 TAMs infiltrate, which is consistent with the correlation of CD44 with CD163 in the TCGA and CGGA glioma datasets (Figure 10B). By checking expression of CD44 and CD163 in glioma samples through IF, CD44+CD163+ cells were uncovered (Supplementary Figure 8), which confirmed that the CD44+ TAMs are in M2 phenotype. We also uncovered that not all T cells express CD44 (Supplementary Figure 7D), and the enrichment analysis results of CD44-related genes in T cells are listed in Supplementary Table 6.

Neftel et al. declare that GBM malignant tumor cells can be grouped into a limited set of cellular states, namely astrocyte-like (AC-like), mesenchymal-1-like (MES1-like), mesenchymal-2-like (MES2-like), oligodendrocyte-progenitor-like (OPC-like), neural-progenitor-1-like (NPC1-like), and neural-progenitor-2-like (NPC2-like) states (32). Because CD44 was described as a marker of GBM CSCs, also known as glioma stem–progenitor cells or glioma-initiating cells (7), we also annotated the stem-like cells besides the six GBM tumor cellular states in malignant cells (Figure 9F). Furthermore, we revealed MES1-like state preference of CD44 (Figure 9G), which is similar to TCGA-mesenchymal subtype preference of CD44 (Figure 3). Surprisingly, no GBM stem-like cells express CD44, suggesting that it is not a suitable CSC marker (Figure 9G). Next, we compared the DEGs between CD44+ tumor cells and CD44- tumor cells (Figure 9H), and the enrichment analysis results verified the immunity function of CD44 (Figure 9I, Supplementary Table 7). In addition, we calculated the Pearson's correlated genes with CD44 in CD44+ tumor cells and uncovered that some positive related genes are immunosuppressive genes (e.g., C10orf54, TGFB1, VEGFA, HAVCR2, and IDO1) and that some negative related genes are immune stimulator genes (e.g., BTN3A2, BTN3A1, TNFSF9, TNF, TNFSF4, IL1A, IL12A, TNFRSF4, CD27, TNFRSF9, CD40, TLR4, HMGB1, and ENTPD1) (Figure 9J, Supplementary Table 8). The enrichment analysis results of CD44-positive related genes in CD44+ malignant cell showed that CD44 correlates with immune activities as well (Figure 9K, Supplementary Table 7).

CD44 is related to immunity in glioma, CD44+ malignant tumor cells are in MES-1-like state, and CD44+ TAMs are in M2 phenotype. Therefore, we hypothesized that glioma CD44+ cells take part in tumor immunosuppression. First, we used the immunomodulators (21) with CD44 to construct the protein–protein interaction network, and the immunomodulators with direct interaction to CD44 were identified (Figure 11A). Then, we found that some of the CD44-correlated genes in CD44+ cells are immunomodulators (Figure 9J, Supplementary Figures 7E,F). For example, in CD44+ T cells, some positive related genes are immunosuppressive genes (e.g., CD274, VEGFA, PDCD1, and CTLA4) and that some negative related genes are immune stimulator genes (e.g., CCL5, ITGB2, ICAM1, and PRF1). Thus, CD44+ malignant tumor cells, CD44+ TAMs, and CD44+ T cells may induce immunosuppression in glioma. Next, we explored the expression of the immunomodulators with direct interaction to CD44 among the glioma CD44+ cells (Figures 11B–D). Each CD44+ cell-type functions immunosuppression through different ways: CD44+ tumor cells express CD276, IL13, VEGFA, and IDO1; CD44+ TAMs express IL10, TGFB1, VEGFA, and HAVCR2; CD44+ T cells express CD274, TGFB1, CTLA4, LAG3, and PDCD1 (Figure 11C). To our surprise, we found that CD44+ T cells express both CD274 (namely PD-L1) and PDCD1 (namely PD-1) in glioma. In addition, we built the regulon network of glioma cells, and the regulon network of CD44+ cells is different from CD44- cells (Figure 12). DDIT3 regulon is upregulated in CD44+ tumor cells, which can inhibit type I interferon (IFN-I) and IFN-stimulated gene production (48). RUNX3 is upregulated in CD44+ T cells, which serves as immunosuppressive role in tumor (49). Therefore, CD44+ glioma cells participate in glioma immunosuppression.

It is an ongoing endeavor to conduct researches on the development of new drugs, and new strategies eradicating cancer target the immune cells with immune checkpoints (50). PD-1 blockade yielded promising results in many cancers, and the Ivy Consortium initiated a multiinstitution, randomized, open-label pilot study of pembrolizumab, an anti-PD-1 monoclonal antibody, in patients with recurrent GBM (51). To uncover the correlation between immunotherapy and CD44 in glioma, we extracted bulk RNA-seq data of patients with recurrent GBM receiving adjuvant, postsurgical PD-1 blockade therapy from this study. The responder and nonresponder patients are different in mRNA expression pattern based on principal component analysis (Figure 13A). Furthermore, the patients with nonresponder recurrent GBM are enriched in PD-1 signaling pathway, suggesting that they have immunosuppression microenvironment and would display therapy resistance (Figure 13B). CD44 and most of its directly interacted inhibitory immunomodulators express relatively higher in immunotherapy nonresponder subgroup than responder subgroup, although some genes are not statistically significant (Figure 13C). These results confirm immunosuppression role of CD44 in glioma.