Abstract

   Glioblastoma (GBM) is classified into subtypes according to the
   molecular expression profile; the proneural subtype has a relatively
   good prognosis, and the mesenchymal type is the most aggressive form
   with the worst prognosis. GBM undergoes proneural–mesenchymal
   transition (PMT) during its evolution or in response to changes in the
   microenvironment or therapeutic interventions. PMT is accompanied by
   infiltration of non-tumor cells, decreased tumor purity, and immune
   evasion. However, the cellular and molecular mechanisms underlying PMT
   remain unclear. Differentially expressed genes (DEGs) were identified
   using GBM transcriptome datasets, and prognostic analysis was performed
   to screen for PMT-related genes (PMTRGs). Consensus cluster analysis
   was followed by Gene Set Enrichment Analysis, Gene Ontology, and Kyoto
   Encyclopedia of Genes and Genomes analyses of DEGs to determine the
   biological functions and pathways regulated by PMTRGs. CIBERSORT,
   TIMER, MCPCOUNTER, and XCELL algorithms were used to analyze immune
   cell infiltration patterns. The TIDE algorithm was used to examine
   immunotherapy scores. The Lasso, Cox, and Step machine learning
   algorithms were used to construct and screen the optimal risk
   assessment prognostic model. PMTRG expression patterns in patient
   tissues and different cell subsets were examined by proteomics and
   single-cell transcriptome data analysis. Seventeen DEGs and prognostic
   PMTRGs were identified in proneural and mesenchymal subtypes.
   PMTRG-related mRNA interactions and protein–protein interaction
   networks were associated with the immune activity of GBM. Consensus
   cluster analysis based on PMTRGs divided GBM into three independent
   subclusters. Functional and pathway analyses showed that PMTRGs were
   highly expressed in the C1 subcluster, which was associated with GBM
   mesenchymal isoforms, pathways, and poor prognosis, and showed stronger
   immune responses. Four immune evaluation algorithms and TIDE analysis
   showed that the C1 cluster had high levels of immune cell infiltration
   and immune molecule scores. The prognostic risk assessment model based
   on PMTRGs can effectively predict the prognosis of GBM patients.
   Proteomic data from immunohistochemistry and single-cell transcriptome
   data suggested that PMTRGs are predominantly expressed in monocytes,
   macrophages, and blood vessels rather than in tumor cells. This study
   identified 17 key genes associated with PMT in GBM. These PMTRGs are
   mainly expressed on immune cells and blood vessels in the GBM
   microenvironment and are associated with poor prognosis, suggesting
   that PMT events mainly arise from the infiltration and activation of
   immune cells derived from the bone marrow and blood vessels. These
   findings provide new evidence and targets for the treatment of GBM.

   Keywords: Glioblastoma, Proneural–mesenchymal transition, Machine
   learning, Multi-omics analysis, Immune microenvironment

Introduction

   Glioma is a type of brain tumor that originates from glial cells, and
   glioblastoma (GBM) is a highly malignant form of glioma classified as
   grade 4 by the World Health Organization. GBM accounts for
   approximately 50% of all central nervous system tumors; the survival of
   GBM patients remains at an average of 15 months, and the 5-year
   survival rate is < 10% despite aggressive treatment such as surgical
   resection combined with postoperative chemoradiotherapy^[52]1,[53]2.
   GBM is characterized by a high degree of intratumoral and intertumoral
   heterogeneity and complex molecular mechanisms, and the development of
   effective drugs is thus challenging. In addition, infiltrating
   monocytes, macrophages, neutrophils, fibroblasts, and aberrant
   vascularization in the GBM microenvironment constitute a niche that
   favors tumor growth and promotes treatment resistance^[54]3–[55]6. The
   infiltration of non-tumor cells and the abnormal vascularization in GBM
   result in the formation of an immunosuppressive microenvironment, which
   leads to a higher degree of malignancy and poor prognosis of tumors;
   the immune cell infiltration profiles in GBM vary significantly among
   different molecular subtypes^[56]7,[57]8. To precisely target GBM, it
   is necessary to elucidate the subtype-specific molecular mechanisms and
   immune infiltration patterns that lead to resistance to GBM
   therapy^[58]9,[59]10.

   In 2010 and 2017, Verhaak and Wang et al. from The Cancer Genome Atlas
   (TCGA) team proposed three or four typing schemes for GBM in which the
   gene expression profiles of the proneural subtype and the mesenchymal
   subtype were significantly different; the prognosis of the proneural
   subtype was relatively good, whereas that of the mesenchymal subtype
   was poor^[60]11–[61]13. The mesenchymal subtype may have evolved from
   the proneural subtype, and proneural-mesenchymal transition (PMT)
   occurs when GBM adapts to changes in the microenvironment or
   chemoradiotherapy^[62]14–[63]16. Single-cell lineage tracking suggests
   that cell subtype transitions under therapeutic conditions may be an
   escape mechanism that can result in the failure of precision
   therapy^[64]17. The proneural subtype of GBM can change into a
   mesenchymal subtype upon recurrence and is accompanied by changes in
   tumor purity and immune cell infiltration profiles; PMT in GBM promotes
   tumor growth, invasion, and treatment resistance, resulting in
   decreased patient survival^[65]18,[66]19. Bone marrow and
   vascular-derived monocytes, as well as resident microglia in the brain
   are transformed into tumor-associated macrophages and undergo M2-type
   polarization after infiltration into the GBM microenvironment, and
   M2-type macrophages are directly involved in PMT events in
   GBM^[67]20,[68]21. Spatial transcriptome data indicate that tumor
   purity decreases over time during the natural evolution of glioma, in
   parallel with an increase in neuronal and oligodendrocyte marker genes
   and tumor-associated macrophages, which is also characteristic of the
   evolution of GBM toward a mesenchymal subtype^[69]22. However, the
   cellular and molecular mechanisms underlying PMT remain to be
   elucidated.

   The development of bioinformatics and single-cell sequencing
   technologies has provided oncologists with new methods to study the
   molecular expression profile of tumors from a macroscopic perspective,
   which enables a deeper understanding of the cellular and molecular
   heterogeneity of GBM and helps guide the development of personalized
   and precision therapies^[70]23–[71]25. Combined analysis of machine
   learning, transcriptome, proteomics, and single-cell sequencing data
   identified 17 differentially expressed genes (DEGs) in proneural and
   mesenchymal subtypes of GBM that were associated with patient
   prognosis. Four different algorithms were used to evaluate the immune
   and prognostic characteristics of these PMT-related genes (PMTRGs).

Methods

Data acquisition and differential gene screening

   RNA sequencing (RNAseq) data and matched clinical information from GBM
   patient tissues for differential gene analysis were downloaded from
   TCGA database ([72]https://portal.gdc.cancer.gov/) and the CGGA
   database ([73]http://www.cgga.org.cn/) ^[74]26,[75]27. The results of
   585 Agilent microarrays in the TCGA dataset, 325 mRNA sequencing
   results, and 301 RNA microarrays in the CGGA database were used for
   analysis of DEGs. RNAseq data from 5 normal brain tissue, 671 low-grade
   glioma (LGG) tissue, and 153 GBM tissue samples were used for
   expression analysis in TCGA database. Analysis of DEGs was performed
   using the R package DESeq2 (v1.32.0). Genes with a false discovery rate
   FDR < 0.05, | Fold change|≥ 2 were considered differentially expressed.
   The RNAseq data of 153 GBM patients in the TCGA database were subjected
   to univariate Cox regression analysis, and the statistical significance
   of differences was evaluated using the log-rank test. Kaplan–Meier
   survival curves were plotted. Spearman correlation coefficient analysis
   was used to identify correlations between the 17 PMTRGs. Protein
   interactions were analyzed and functional annotations were performed
   using the GeneMANIA website ([76]https://genemania.org/)^[77]28. The
   expression of different isoforms of PMTRGs and the levels of GBM were
   analyzed using GraphPad (v10.2.3).

Consensus cluster analysis

   Consensus ClusterPlus (v1.54.0) was used to perform consistency
   analysis; the maximum number of clusters was 6, and 80% of the total
   samples were extracted from 100 replicates, clusterAlg = “hc”,
   innerLinkage = 'ward. D2’. The clustering heatmaps were analyzed using
   the R package Pheatmap (v1.0.12). Survival differences between
   subclusters were analyzed using Survival (v3.5–5), and the statistical
   significance of differences was analyzed using the log-rank test;
   Kaplan–Meier survival curves were plotted. The proportional histogram
   tool from the Bioinformatics website
   ([78]http://bioinformatics.com.cn/) was used to plot the distribution
   of different subclusters in the four molecular isotypes of GBM^[79]29.
   The expression of PMTRGs in different subclusters was analyzed using
   GraphPad (v10.2.3).

Gene set enrichment analysis (GSEA), gene ontology (GO), and kyoto
encyclopedia of genes and genomes (KEGG) analysis

   The marker gene sets were downloaded from the MSjgDB website
   ([80]https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). The
   permutation was set to 1000, and GSEA enrichment analysis was performed
   on patients with different clusters using GSEA software
   (v4.3.3)^[81]30,[82]31. The differential analysis tool of the Assistant
   for Clinical Bioinformatics website
   ([83]https://www.aclbi.com/static/index.html#/) was used to analyze
   DEGs between the C1 and C3 clusters. GO functional enrichment and KEGG
   pathway analysis of the DEGs were performed using the R package
   ClusterProfiler (v3.18.0).

Immune landscape analysis

   The infiltration profiles and estimated proportions of 22 immune cells
   per sample were calculated using the CIBERSORT method in the R package
   IOBR^[84]32. Immune cell infiltration profiles were examined and immune
   scores, matrix scores, and immune microenvironment scores were
   calculated using the XCELL algorithm on the XCELL website
   ([85]https://aran-lab.com/software/xcell/)^[86]33. TIMER and MCPCOUNTER
   algorithms in The Immune Scoring Tool
   ([87]https://www.aclbi.com/static/index.html#/) were used to assess the
   immune cell infiltration profile and analyze the correlation between
   PMTRGs and immune cell infiltration^[88]34,[89]35. Two different tumor
   immune escape mechanisms and PD-L1, CD8, CAF, IFN-γ, and Merck18 scores
   were evaluated using the TIDE algorithm
   ([90]http://tide.dfci.harvard.edu/)^[91]36. The expression of immune
   checkpoint molecules in different clusters was mapped using GraphPad
   (v10.2.3).

Construction of a prognostic model

   A prognostic model based on the 17 PMTRGs was constructed using the
   Lasso, Cox, and Step algorithms in the prognostic model tool of the
   Assistant for Clinical Bioinformatics
   ([92]https://www.aclbi.com/static/index.html#/) website. Decision curve
   analysis (DCA) was performed using the ggDCA R package to examine the
   predictive performance of different models. The risk score of the best
   Cox multivariate regression model = (0.1303)*SAA1 + (0.0795)*NPY2R +
   (0.1204)*CCL2 + (0.709)*SAA2 + (0.113)*FCGR2B + (0.0952)*RARRES1 +
   (0.1077)*GJB2 + (− 0.0691)*LTF + (0.1798)*DDIT4L + (− 0.1022)*C1S +
   (− 0.0358)*CHI3L2 + (− 0.0947)*MARCO + (0.124)*AQP9 + (− 0.0394)*LY96 +
   (0.0245)*CCL18 + (− 0.0713)*FAM81B + (− 0.0578)*LOX. Receiver Operating
   Characteristic (ROC) curve analysis was performed at 1, 2, and 3 years
   using R package Proc (v1.17.0.1). The CI function of Proc was used to
   evaluate the area under the curve (AUC) and confidence intervals (CIs).
   The prognostic models were analyzed using the R package Survival.
   Sankey plots of prognostic risk according to patient subclusters, GBM
   typing, age, and gender were generated using the Sankey plot tool on
   the OKX website ([93]https://cloud.oebiotech.com/#/home).

Immunohistochemical staining analysis

   Ten PMTRGs with available immunohistochemical data of protein
   expression in glioma tissues from the HPA website
   ([94]https://www.proteinatlas.org/) were analyzed. Representative
   staining patterns of tumor tissues from two high grade glioma patients
   and two low grade glioma patients were selected, and the staining
   intensity, specificity, and positivity rate of the target protein in 11
   or 12 patients were determined.

Single-cell sequencing analysis

   Single-cell RNA sequencing (scRNA-seq) data from 3589 cells from the
   tumor core and peritumoral core of four GBM patients were obtained from
   the TISCH database ([95]http://tisch.comp-genomics.org/gallery/)
   [96]GSE84465 dataset^[97]37, including tumor cells and representative
   cells from each of the major central nervous system cell types
   (vascular, immune, neuronal, and glial cells). Full-length scRNA-seq
   data of 8252 T cells in 31 GBM patients were acquired from the
   GSE163108_10x dataset^[98]38. The single-cell data were processed and
   analyzed using the R package MAESTRO and the Seurat package, and cell
   clustering was confirmed using the t-SNE method to analyze the
   expression levels of different PMTRGs in the two datasets.

Statistical analysis

   All R package analyses were performed using R software (v4.0.3). SPSS
   (v26.0) and GraphPad (v10.2.3) were used to perform statistical
   analyses, and comparisons between the two groups were performed using
   unpaired t-tests. The Wilcox test was used for comparisons between more
   than two groups. * P < 0.05, ** P < 0.01, *** P < 0.001, ns: not
   significant.

Results

Identification of PMTRGs in GBM

   DEGs were identified in patients with the mesenchymal and proneural
   subtypes of GBM (control group) in TCGA and CGGA datasets with the
   threshold set to |Fold Change|> 2 and FDR < 0.05. A total of 141
   upregulated genes and 148 downregulated genes were identified from TCGA
   dataset, 3156 upregulated genes and 1621 downregulated genes were
   identified in the CGGA325 dataset, and 366 upregulated genes and 121
   downregulated genes were identified in the CGGA301 dataset
   (Fig. [99]1A). A Venn plot was generated to identify differentially
   upregulated genes that intersected in the three datasets, and 83
   commonly upregulated genes were identified (Fig. [100]1B). Batch
   survival analysis of the 83 genes in 153 patients with GBM from TCGA
   identified 17 target genes with prognostic significance; these genes
   were categorized as PMTRGs. Spearman correlation analysis of the mRNA
   expression of PMTRGs showed a broad correlation among the 17 genes
   (Fig. [101]1C). Protein–protein interaction analysis and functional
   annotation were performed. The first functional pathways identified
   showed that the proteins encoded by these genes were involved in immune
   response, chemokine response, and chemokine receptor binding,
   suggesting that the main functions of PMTRGs were related to tumor
   immune regulation (Fig. [102]1D). The expression of PMTRGs in each
   subtype of GBM was characterized, and the results showed that 16 genes
   except DDIT4L were significantly overexpressed in the mesenchymal
   subtype (Fig. [103]1E).

Fig. 1.

   [104]Fig. 1
   [105]Open in a new tab

   Screening and expression analysis of PMTRGs. (A) Volcanic plot of
   proneural and mesenchymal GBM transcriptome differential genes from the
   TCGA, CGGA325, and CGGA301 platforms. (B) Venn diagram of
   differentially upregulated genes in the three platforms. (C) mRNA
   correlation analysis of 17 PMTRGs. (D) Protein–protein interaction
   network and functional annotation of PMTRGs. (E) Expression profiles of
   PMTRGs in four subtypes of GBM.

Prognostic analysis of PMTRGs in GBM

   The forest plot of univariate Cox regression analysis showed that the
   expression of PMTRGs in GBM significantly affected the survival of
   patients (P < 0.05, hazard ratio > 1) (Fig. [106]2A). The histogram
   included the expression of PMTRGs in normal tissues, LGGs, and GBM.
   Except for FAM81B, the expression of 16 genes was higher in GBM than in
   LGG and normal brain tissues (Fig. [107]2B). Kaplan–Meier survival
   analysis of 153 GBM patients showed that survival time was
   significantly higher in patients with high PMTRG expression level than
   in those with low expression (Fig. [108]2C). These results suggest that
   PMTRGs are significantly differentially expressed in glioma and are
   associated with shorter survival in GBM patients.

Fig. 2.

   [109]Fig. 2
   [110]Open in a new tab

   Prognostic and expression analysis of PMTRGs. (A) Forest plot of PMTRGs
   by univariate Cox analysis. (B) Expression profile of PMTRGs in normal
   tissues, LGG, and GBM. (C) Kaplan–Meier survival curves of 17 PMTRGs in
   GBM patients.

GBM consensus cluster analysis based on PMTRGs

   Consensus cluster analysis of 17 PMTRGs divided the 153 GBM patients
   into three subtypes according to the decrease in the Delta curve as
   follows: C1 (n = 57), C2 (n = 55), and C3 (n = 41) (Fig. [111]3A–C).
   Principal component analysis (PCA) showed that patients were clustered
   into three distinct subclusters (Fig. [112]3D), and the clustering
   heatmap showed that PMTRGs were significantly differentially expressed
   in the three subclusters (Fig. [113]3E). Kaplan–Meier analysis showed
   that patients in the C1 and C2 clusters had a significantly shorter
   survival than those in the C3 cluster (Fig. [114]3F). The histogram
   suggested a potential correlation between the three subclusters and the
   four GBM subtypes, with the C1 subcluster predominantly resembling the
   mesenchymal subtype (61.4% of patients have a mesenchymal origin), the
   C3 subtype showing high similarity to the proneural subtype (61% of
   patients have a proneural origin), and the C2 subtype showing a
   transitional subtype between the other two (The proportions of
   mesenchymal and proneural patients are 25.9% and 16.7%) (Fig. [115]3G).
   The expression profiles of PMTRGs in the three clusters showed that the
   expression of 16 genes except DDIT4L was significantly higher in the C1
   cluster than in the other two clusters (Fig. [116]3H).

Fig. 3.

   [117]Fig. 3
   [118]Open in a new tab

   Consensus cluster analysis. (A–C) Consensus cluster analysis of PMTRGs
   divided GBM patients into three subclusters. (D) PCA showed a
   relatively independent distribution trend of the three subclusters. (E)
   Expression heatmap of PMTRGs in the three subclusters. (F) Kaplan–Meier
   survival curves in different clusters of patients. (G) Histogram of
   three subcluster patients and the four subtypes of GBM. (H) Expression
   of PMTRGs in patients in the three clusters.

Subcluster functional enrichment and DEG analysis

   The association of subclusters with regulatory pathways was examined by
   GSEA of the C1 and C3 subclusters, which had the largest differences in
   PMTRG expression and survival rates. The results showed that the C1
   cluster was significantly correlated with the mesenchymal subtype of
   GBM, epithelial-mesenchymal transition, angiogenesis, hypoxia, fatty
   acid metabolism, inflammatory response, and IL6_JAK_STAT3 pathway, and
   negatively correlated with the proneural subtype of GBM. The C3 cluster
   showed the opposite pattern of enrichment than the C1 cluster; the C3
   isoform was significantly positively correlated with the proneural
   subtype of GBM and negatively correlated with the remaining pathways
   (Fig. [119]4A). DEGs were analyzed in the C1 and C3 clusters, and the
   results of the differential analysis are shown in the volcano map and
   cluster heat map in Fig. [120]4B,C. Next, we performed GO function
   enrichment and KEGG pathway analysis based on the DEGs. The results of
   the GO analysis showed that the molecular functions of immune
   receptors, chemokines, and cytokines were significantly upregulated.
   Cytokine receptor interactions in the KEGG pathway and a significant
   activation of chemokine signaling pathways suggested that cell-to-cell
   communication and immune activity are stronger in the C1 cluster than
   in the C3 cluster (Fig. [121]4D).

Fig. 4.

   [122]Fig. 4
   [123]Open in a new tab

   Functional enrichment and differential analysis of subclusters. (A)
   Enrichment analysis of the C1 and C3 clusters in eight GSEA datasets.
   (B,C) Volcano map and cluster heat map of DEGs in the C1 and C3
   clusters. (D) GO and KEGG analyses of DEGs in the C1 and C3 clusters.

Analysis of subcluster immune patterns

   Previous protein–protein interaction and functional enrichment analyses
   suggested a potential link between PMTRGs and GBM immune activity. We
   therefore examined immune cell infiltration in the three clusters of
   GBM. The results of the CIBERSORT algorithm showed that M2 macrophages
   were the most abundant immune cells, and they were present in
   significantly higher numbers in the C1 cluster than in the C2 and C3
   clusters (Fig. [124]5A,B). Consistently, the MCPCOUNTER algorithm
   suggested that the main infiltrating immune cells in each subcluster
   were monocytes and macrophages, and the infiltration abundance was
   higher in the C1 cluster than in the C2 and C3 clusters (Fig. [125]5C).
   The TIMER algorithm showed that neutrophils, macrophages, and bone
   marrow dendritic cells were significantly more abundant in the C1
   cluster than in the C2 and C3 clusters (Fig. [126]5D). The XCELL
   algorithm showed that the infiltration abundance of monocytes and
   macrophages was significantly higher in the C1 cluster than in the C2
   and C3 clusters, and the immune score, matrix score, and immune
   microenvironment score were significantly higher in the C1 cluster than
   in the C2 and C3 clusters (Fig. [127]5E,F). The TIDE score did not
   differ significantly between the three clusters; however, immune cell
   dysfunction was higher in the C1 subcluster than in the C2 and C3
   groups, whereas the exclusion score indicated predominant immune cell
   dysfunction in the C1 subcluster and insufficient immune cell
   infiltration in the C3 subcluster (Fig. [128]5G). Further analysis
   showed that the immune scores for PD-L1, CD8, CAF, IFN-γ, and Merck18
   were significantly higher in the C1 cluster than in the C2 and C3
   subclusters (Fig. [129]5H). Assessment of 24 immune checkpoint
   molecules in the three subclusters showed that the expression of
   molecules such as TGFB1, VSIR, PD-L1, IDO1, IL10, HAVCR2, CTLA4, TIGIT,
   and CD276 was significantly higher in the C1 subcluster than in the C2
   and C3 subgroups (F[130]ig. [131]5I).

Fig. 5.

   [132]Fig. 5
   [133]Open in a new tab

   Immune landscape analysis in different subclusters. (A,B) The CIBERSORT
   algorithm was used to evaluate the infiltration of different immune
   cells in the three clusters. (C,D) The MCPUNTER and TIMER algorithms
   were used to evaluate immune cell infiltration in the three clusters.
   (E,F) The XCELL algorithm evaluates immune cell infiltration abundance
   and immune score, matrix score, and immune microenvironment score.
   (G,H) TIDE immunoassay and immune molecule scores. I: Gene expression
   analysis of 24 immune checkpoints in the three clusters.

GBM risk prognosis model based on PMTRGs

   The prognostic value of PMTRGs in GBM was further examined by
   constructing a prognostic model using three methods: the Lasso
   algorithm, multivariate Cox regression, and Lasso modeling followed by
   Step function calculation. ROC curves were generated to evaluate the
   predictive performance of the model (Fig. [134]6A–C), and the optimal
   predictive performance of the Cox multivariate regression model was
   selected by DCA (Fig. [135]6D). Kaplan–Meier survival analysis showed
   that the survival time was significantly shorter in the high-risk group
   (Fig. [136]6E). The relationship between the high- and low-risk groups
   identified by the prognostic model and the different subclusters of
   patients, GBM subtypes, age, and gender is shown in the Sankey diagram
   in Fig. [137]6F. Through a set of three consecutive figures, we further
   demonstrated the correlations among the prognostic model scores, the
   patients’ survival, and the expressions of the 17 PMTRGs within the
   model (Fig. [138]6G).Finally, we evaluated the correlation between the
   17 PMTRGs and the infiltration of different immune cells. The TIMER and
   MCPCOUNTER algorithms identified a strong correlation between PMTRGs
   and monocytes, macrophages, and bone marrow dendritic cells
   (F[139]ig. [140]6H,I).

Fig. 6.

   [141]Fig. 6
   [142]Open in a new tab

   GBM prognostic model construction and analysis. (A–C) ROC curve of the
   prognostic risk model constructed by three machine learning algorithms:
   Lasso, multivariate Cox analysis, and Step function. (D) The optimal
   prognostic model was screened by decision curve analysis. (E)
   Kaplan–Meier survival curve analysis of the high-risk and low-risk
   groups identified by the prognostic model; prognostic risk assessment
   using a Sankey diagram of the relationship between risk groups and
   patient subclusters, GBM subtype, age, and gender. (G) Triplet diagram
   established by the prognostic model. (H) Analysis of the correlation
   between PMTRGs and various immune cells.

Immunohistochemical analysis of PMTRGs in glioma tissues

   To examine the expression of proteins encoded by PMTRGs in glioma
   tissues, 10 PMTRGs with available protein expression information in all
   glioma tissues from the HPA website were identified. A total of 11 or
   12 stained tissues for each protein were used for analysis, and
   immunohistochemical staining results were obtained from two high-grade
   and two low-grade patients as representative samples. The intensity,
   sensitivity, and positivity rate of immunohistochemical staining was
   examined in all 11 or 12 patients (Fig. [143]7A,B). The results showed
   that the retinoic acid receptor responder 1 (RARRES1) gene was
   expressed at high levels in glioma; expression data on gap junction
   beta-2 (GJB2) and chitinase 3-like 2 (CHI3L2) was available for few
   patients, and C–C motif chemokine ligand 2 (CCL2) was expressed in the
   tumor stroma and perivascular areas, whereas no positive staining was
   observed in glioma cells. Fc Gamma Receptor IIb (FCGR2B), LTF,
   complement component 1S (C1S), macrophage receptor with collagenous
   structure (MARCO), aquaporin 9 (AQP9), and FAM81B were not detected in
   the tumor cells of these patients. These results indicate that most
   PMTRGs are not expressed in GBM tumor cells or are expressed in only a
   few patients, suggesting that the mechanism by which PMTRGs promote
   mesenchymal transition and result in poor clinical outcomes in glioma
   patients may not occur within tumor cells.

Fig. 7.

   [144]Fig. 7
   [145]Open in a new tab

   Immunohistochemical staining of PMTRGs in the tissues of glioma
   patients. (A) Schematic diagram of the expression of 10 PMTRGs in
   glioma tissues detected by immunohistochemical staining. (B) Staining
   intensity, sensitivity, and positivity rate of 10 PMTRGs in all 11 or
   12 glioma patient tissue samples.

Single-cell sequencing analysis of the expression of PMTRGs in different
subsets of GBM

   Single-cell sequencing analysis of 3589 cells from four GBM patients in
   the GSE86645 dataset classified cells into eight subsets. Analysis of
   the expression of PMTRGs in these cell subsets showed that CCL2,
   FCGR2B, MARCO, AQP9, and lymphocyte antigen 96 (LY96) were mainly
   expressed in monocytes and macrophages, whereas C1S, CCL2, and LOX were
   highly expressed in blood vessels. Moreover, C1S, LOX, CCL2, RARRES1,
   and GJB2 are also expressed at low levels in cells such as Astrocyte,
   Neuron, and OPC. (Fig. [146]8A,B). Full-length scRNA-seq of 8252 cells
   from 31 GBM patients in the GSE163108_10x dataset classified cells into
   six subsets including T cells, oligodendrocytes, monocytes, and
   macrophages. The results showed that FCGR2B, LY96, CCL2, AQP9, MARCO,
   and RARRES1 were mainly expressed in monocytes and macrophages
   (Fig. [147]8C,D). These results indicated that the effect of PMTRGs on
   PMT and the poor prognosis of GBM patients is related to their
   expression in immune cells, especially monocytes and macrophages.

Fig. 8.

   [148]Fig. 8
   [149]Open in a new tab

   Expression of PMTRGs in different cell subsets. (A,B) Schematic and
   statistical analysis of the expression of 16 PMTRGs in cell subsets
   from the [150]GSE84465 dataset. (C,D) Schematic and statistical
   analysis of the expression of 10 PMTRGs in cell subsets from the
   GSE163108_10x dataset.

Discussion

   GBM is a refractory central nervous system tumor that presents several
   treatment challenges such as resistance and recurrence, and its
   heterogeneity and invasive growth pattern complicate surgery and drug
   treatment. The molecular expression pattern and immune microenvironment
   differ among GBM subtypes, underscoring the importance of considering
   the different subtypes of GBM for research purposes and in the
   development of treatments to meet the requirements of precision
   therapy^[151]4,[152]39,[153]40. PMT in GBM occurs in response to
   various treatments and changes in the microenvironment, and results in
   the transformation of tumors from the proneural subtype with a
   relatively good prognosis into the most aggressive mesenchymal subtype,
   thereby promoting immune evasion^[154]41. Elucidating the cellular and
   molecular biological mechanisms underlying PMT is critical for the
   development of targeted therapies^[155]42. We identified 17 PMTRGs that
   were commonly differentially upregulated and had significant prognostic
   value by screening DEGs in patients with proneural and mesenchymal GBM
   in different databases. There is a wide range of correlations and
   interactions between the RNA and protein expression of these PMTRGs in
   GBM, and the protein function annotation suggests that PMTRGs are
   mainly involved in immune response, chemokine receptor binding, and
   chemokine response. Consensus cluster analysis according to PMTRG
   expression divided the 153 GBM patients into three distinct
   subclusters, with the C1 cluster showing the highest PMTRG expression
   levels and the worst prognosis. GSEA, GO, and KEGG analysis showed that
   the C1 cluster had mesenchymal and epithelial-mesenchymal transition
   characteristics and was highly correlated with tumor hypoxia,
   angiogenesis, fatty acid metabolism, inflammatory response, and the
   STAT3 pathway, providing evidence that high expression levels of PMTRGs
   are associated with the malignant biological behavior and poor
   prognosis of GBM. Targeting macrophages in the hypoxic environment has
   shown therapeutic potential regarding vasculature normalization in
   GBM^[156]43, and monocytes in the GBM microenvironment have also been
   implicated as important targets for glioma mesenchymal
   transition^[157]44.

   Immunotherapy is highly successful in tumors other than GBM, but has
   had limited efficacy in GBM, an "immune-cold" tumor^[158]45. In
   addition to tumor cells, GBM tumor tissues show infiltration of various
   immune cells and immune matrices, which provide a niche for tumor
   growth and contribute to the immunosuppressive microenvironment. High
   levels of immune cell infiltration are associated with lower tumor
   purity, a high degree of malignancy, an invasive growth pattern, and
   treatment resistance in GBM^[159]46–[160]48. The recruitment and
   activation of bone marrow-derived monocytes and macrophages plays a key
   role in the tumor immunosuppressive microenvironment^[161]3,[162]49.
   Time-resolved single-cell transcriptomic evidence suggests that bone
   marrow-derived monocytes differentiate into M2 macrophages within 24 h
   of entry into GBM tissue and exert immunosuppressive functions^[163]50.
   The development of drugs targeting immune cells in the tumor
   microenvironment is considered one of the promising therapeutic options
   for GBM in the future^[164]51,[165]52. Toosendanin, a small molecule
   compound targeting tumor-associated macrophages in GBM, can reverse the
   sensitivity of immunosuppressive microenvironment-sensitizing CAR-T
   therapies for the treatment of GBM^[166]53. The CSF1R-targeting drug
   GW2580 can support antitumor T cell responses and directly exert
   antitumor effects by reprogramming M2-type macrophages to the
   M1-type^[167]54.

   In this study, analysis of four different immune scoring algorithms
   indicated that the C1 cluster had high levels of monocyte and
   macrophage infiltration, indicating that the expression of PMTRGs may
   be involved in immune cell infiltration and the formation of a tumor
   immunosuppressive microenvironment. PMTRGs provide a new understanding
   of the complex immune-related expression profiles in proneural and
   mesenchymal GBMs. The Cox multivariate regression risk model based on
   PMTRGs constructed in this study showed a good predictive performance
   for the prognosis of GBM patients. A strong correlation between PMTRGs
   and monocytes, macrophages, and bone marrow-derived dendritic cells
   supported their involvement in the recruitment and activation of bone
   marrow-derived immune cells in addition to monocyte infiltration.
   Evidence from immunohistochemistry and single-cell sequencing in
   clinical GBM patients suggests that the expression of PMTRGs in GBM is
   not localized to tumor cells, but mainly to monocytes, macrophages, and
   blood vessels. These results support that PMTRGs are involved in the
   recruitment of bone marrow and vascular-derived monocytes and
   macrophages through their expression in immune cells and blood vessels
   and contribute to the immunosuppressive microenvironment of mesenchymal
   GBM. Because of their high levels of expression in monocytes,
   macrophages, and blood vessels, PMTRGs are expected to be novel targets
   in the immune microenvironment for the treatment of GBM.

   The expression and function of 17 PMTRGs were intrinsically different.
   The results of immunohistochemistry and single-cell transcriptome
   analysis showed that CCL2, FCGR2B, AQP9, LY96, and MARCO were the main
   immune-related genes expressed on monocytes and macrophages and played
   a central role in PMT in GBM. CCL2 is a major chemokine involved in the
   recruitment of bone marrow and vascular-derived monocytes in
   GBM^[168]55. FCGR2B is an important regulator and therapeutic target of
   the GBM immune microenvironment^[169]56. The membrane channel protein
   AQP9 is associated with the GBM tumor microenvironment and glioma
   grade^[170]57. Myeloid differentiation factor LY96 is a marker gene for
   metastasis and poor prognosis in prostate tumors^[171]58. The
   macrophage receptor MARCO was identified as a pro-tumor marker of
   mesenchymal GBM in a previous study, which is consistent with our
   findings^[172]59. RARRES1 was the only PMTRG confirmed in GBM cells and
   had a strong correlation with all five PMTRGs expressed on monocytes
   and macrophages (Spearman R > 0.5), suggesting that it is involved in
   the interaction of GBM cells and macrophages in the tumor
   microenvironment. In addition, the high vascular expression of
   complement C1S, CCL2, and L-amino acid oxidase LOX suggests that they
   interact with vascular-derived immune cells^[173]60,[174]61. However,
   immunohistochemistry and single-cell transcriptome data of the 17
   PMTRGs were incomplete, and expression data were missing for some of
   these genes, limiting a comprehensive and in-depth analysis. Further
   studies are needed to elucidate their role in PMT events.

Conclusion

   We screened 17 PMTRGs with prognostic significance based on GBM
   transcriptome and clinical information from multiple sequencing
   platforms. Consensus cluster analysis and function and pathway
   enrichment analysis of DEGs determined the correlation between PMTRGs
   and GBM malignant behavior and immune cell infiltration. Four immune
   scoring algorithms identified different immune infiltration profiles
   based on PMTRGs, and the Cox risk prognostic model evaluated the
   predictive value of PMTRGs. Protein expression and single-cell
   sequencing demonstrated that PMTRGs mediate the recruitment and
   activation of monocyte-macrophages and the immunosuppressive
   microenvironment of GBM through expression in immune cells and blood
   vessels. This study provides new insight into the cellular and
   molecular mechanisms underlying PMT events in GBM and suggests
   molecular targets in monocyte-macrophages in the tumor microenvironment
   to reverse the immunosuppressive microenvironment in GBM.

Author contributions

   CX, ZZ, and RX conceptualized the manuscript. HX, XC, and YZ collected
   transcriptome data and performed immune infiltration analysis. MG, LH,
   CH and QF performed the acquisition and analysis of single-cell
   transcriptome and immunohistochemical staining data. WL, YW, JZ, and YY
   performed consensus cluster analysis and prognostic model construction.
   CX, JY and ZZ prepared the figures and wrote the paper. The research
   work was funded by RX. All authors read and approved the final
   manuscript.

Funding

   This research was supported by grants from the National Natural Science
   Foundation of China (No.82171355). National Natural Science Foundation
   of China (No: 82403762). China Postdoctoral Science Foundation (No:
   2023M733164).

Data availability

   All data generated or analysed during this study are included in this
   published article.

Declarations

Competing interests

   The authors declare no competing interests.

Footnotes

   Publisher’s note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

   Can Xu and Jin Yang contributed equally to this work and share the
   first authorship.

Contributor Information

   Zhaomu Zeng, Email: zzmhemisphere@163.com.

   Ruxiang Xu, Email: xuruxiang1123@163.com.

References