Abstract Gliomas, particularly glioblastomas, are highly malignant brain tumors with high recurrence rates and poor prognosis. Despite advances in treatment, recurrence remains a major challenge. Epithelial-mesenchymal transition (EMT) plays a key role in tumor invasion and recurrence. This study explores the transcriptional and regulatory mechanisms driving glioma recurrence, focusing on mesenchymal-like (MES-like) subpopulations. Single-nucleus RNA sequencing was performed on 52 IDH wild-type GBM specimens, including 26 primary and 26 recurrent tumors. Spatial transcriptomics data were also incorporated. Tumor subpopulations were identified through gene regulatory network analysis, copy number variation detection, and nonnegative matrix factorization. Functional validation was conducted using gene knockdown experiments, followed by xenograft studies. We discovered novel MES-like subpopulations in recurrent GBM enriched with EMT-related genes like EGR1 and SERPINE1. These subpopulations exhibited increased transcriptional activity and were associated with poor prognosis and invasiveness. Knockdown of SERPINE1 significantly reduced cell proliferation and migration. Spatial transcriptomics showed MES-like cells concentrated at the tumor margins, highlighting their role in invasion and recurrence. MES-like subpopulations, driven by EGR1 and SERPINE1, are critical in GBM. Targeting these regulators could offer new therapeutic strategies to reduce glioma recurrence and improve outcomes. Supplementary Information The online version contains supplementary material available at 10.1186/s40478-025-02036-6. Keywords: Glioblastoma, Epithelial-mesenchymal transition, Single-cell RNA sequencing, Recurrence, MES-like subpopulations Introduction Gliomas are the most common malignant tumors of the central nervous system. Despite receiving standard treatments, including maximal safe resection combined with radiotherapy and temozolomide chemotherapy, recurrence is inevitable [[46]1]. High-grade gliomas are particularly prone to recurrence, with a 5-year survival rate of less than 10% [[47]2, [48]3]. Additionally, the clinical significance of reoperation in recurrent glioblastoma patients remains unclear. Therefore, understanding the temporal evolution of tumors during treatment is crucial for addressing glioma recurrence. By analyzing gene expression data from primary and recurrent tumors, studying tumor evolution at the transcriptional level offers a novel perspective on the progression of glioblastoma. Recurrent gliomas are often associated with increased malignancy, making subsequent treatments more challenging and leading to poorer prognosis. The epithelial-mesenchymal transition (EMT) is regulated by inducers, transcription factor families, and a series of signaling pathway genes, and is implicated in the invasion and progression of gliomas [[49]4, [50]5, [51]6, [52]7, [53]8, [54]9]. The EMT process is frequently reactivated, enhancing the invasive and metastatic capabilities of cancer cells [[55]10, [56]11]. This suggests that recurrent gliomas, potentially driven by residual invasive cells, may undergo a multi-step reactivation process. This process facilitates adhesion to the extracellular matrix, secretion of enzymes that degrade surrounding matrices, and migration within the glioma microenvironment prior to tumor formation [[57]12, [58]13]. Therefore, investigating the mechanisms by which EMT influences GBM invasion and migration is crucial for developing new and effective therapeutic strategies. In recent years, the rapid development of single-cell sequencing technology has significantly advanced the study of glioma tumor plasticity, the characteristics of the tumor microenvironment, and the role of the microenvironment in tumor plasticity. This progress has also brought potential new opportunities for glioma treatment. Currently, most studies categorize glioma tumor cells into neural progenitor cell-like (NPC-like), mesenchymal-like (MES-like), oligodendrocyte precursor cell-like (OPC-like), and astrocyte-like (AC-like) states [[59]14, [60]15, [61]16]. Intratumoral heterogeneity and the dynamic plasticity of cell states are hallmarks of malignant brain tumors [[62]14, [63]16, [64]17, [65]18]. Previous studies, including those from our group and others, have demonstrated that the MES-like cell state is associated with glioma invasion [[66]16, [67]19, [68]20, [69]21]. However, research on tumor cell states remains limited, which may be a critical factor contributing to tumor recurrence and drug resistance. A comprehensive understanding of the transcriptional regulatory mechanisms in primary and recurrent gliomas may aid in developing new therapeutic strategies and improving patient outcomes. In this study, we utilized extensive single-cell sequencing data from matched primary and recurrent gliomas to identify genes with high transcriptional regulatory activity in different tumor cell states. Spatial transcriptomics analysis was conducted to determine the spatial localization of tumor cell subpopulations and the transcription factor EGR1. We validated the expression and differentiation of EGR1 and SERPINE in tumor microenvironment subpopulations using single-cell sequencing data from wild-type mice at different stages. Additionally, we performed cellular functional assays and animal experiments to verify our findings. Our research suggests that in recurrent gliomas, EGR1 may mediate EMT and that the SERPINE-MES subpopulation could be a key factor in glioma recurrence and drug resistance. Materials and methods Data collection This study performed single-nucleus RNA sequencing analysis on 52 human IDH wild-type GBM specimens, comprising 26 primary tumors and 26 matched recurrent tumors from the same patients. All patients underwent standard surgical treatment, radiotherapy, and chemotherapy. The collection and processing of tumor specimens followed protocols described in previous literature based on original data [[70]22]. Transcriptomic and clinical data was downloaded from the Chinese Glioma Genome Atlas (CGGA) database([71]http://www.cgga.org.cn/), including mRNA-693 cohort (n = 693), mRNA-325 cohort (n = 325), Rembrandt cohort (n = 475), and glioma cohort from The Cancer Genome Atlas (TCGA) (n = 702)([72]https://www.cancer.gov/ccg/research/genome-sequencing/tcg a). Spatial transcriptomics sequencing data was obtained from previous studies and is now publicly available ([73]GSE194329) [[74]19]. Patient sample collection and study design for spatial transcriptome data have received ethical approval from West China Hospital, Sichuan University (2020.837). The spatial transcriptomics data included three human IDH wild-type primary GBM specimens, two recurrent glioblastoma specimens, and one peritumoral tissue sample. Early (day 7) and late (day 28) IDH wild-type GBM mouse single-cell transcriptomic data was obtained from the publicly available GEO database ([75]GSE166218) [[76]23]. Single-cell quality control processing workflow The quality control parameters were set to allow a minimum of 500 genes and a maximum of 4000 genes per cell, with mitochondrial gene content limited to 20%. The FindVariableFeatures function was used to identify 2000 highly variable genes. Batch effects across different tumor samples were corrected using the R package harmony. Cell type identification and annotation were performed by referencing the CellMarker database and previous literature [[77]24, [78]25]. The R package SCP was used to create bar plots for cell proportion visualization and heatmaps for marker gene expression analysis. Single-cell copy number variation (CNV) analysis Single-cell CNV analysis was conducted based on single-cell transcriptome data to identify cells with copy number abnormalities, aiding in the inference of malignant tumor cells. To distinguish between tumor cells and non-malignant cells, we used the R package inferCNV [[79]26]. T cells and macrophages were set as the reference cell types. After quantifying the degree of copy number variation using the CNV score, cells were classified into tumor cells and normal cells. Subsequently, we annotated subpopulations of tumor cells using marker genes identified from previous literature [[80]25, [81]27, [82]28, [83]29, [84]30]. Gene regulatory network analysis The transcriptional states of microenvironmental cells in single-cell sequencing were derived from underlying gene regulatory networks [[85]31]. These networks included transcription factors, cofactors, and their downstream target genes. Using single-cell regulatory network inference and clustering (SCENIC) analysis, we effectively identified co-expression modules (regulons) between transcription factors (TFs) and potential target genes [[86]18, [87]31]. We analyzed the differences in transcriptional activity between IDH wild-type primary and recurrent glioma tumor cell subpopulations. Identifying MES subpopulation cells using non-negative matrix factorization (NMF) We constructed an expression matrix using genes from the EMT signaling pathway and excluded cells without EMT gene expression. EMT-related genes were obtained from The Molecular Signatures Database (MSigDB) (Supplement Table [88]1). The R package NMF was used to decompose the MES cell subpopulation [[89]32]. The function FindAllMarkers was employed to identify characteristic genes for each cluster. Subsequently, we redefined new MES cell subpopulations. Single-cell subpopulation stemness analysis We assessed the differentiation potential of glioma tumor cell subpopulations and MES subpopulations at single-cell level. The R package CytoTRACE was used to predict the differentiation capabilities of different subpopulation cells. The scores were visualized using bar plots and UMAP plots [[90]33]. Constructing cellular trajectories Cell subpopulation trajectory inference is crucial for studying the differentiation structure and sequence of cell lineages. We employed the R packages “monocle2” and “Slingshot” to infer the differentiation trajectories of MES cell subpopulations and displayed the enriched signaling pathways within these MES subpopulations. Pathway enrichment analysis Metabolic changes in cells can significantly influence tumorigenesis and development. We further analyzed pathway enrichment in MES stromal subpopulations using the R package scMetabolism for metabolic pathway enrichment analysis. Single-cell flux estimation analysis (scFEA) was used to evaluate metabolic changes in the epithelial-mesenchymal MES subpopulation [[91]34]. Heatmaps were generated to display the expression of immune-related signaling genes in the novel MES stromal subpopulation. Spatial transcriptomics analysis Spatial transcriptomics data were generated using the 10x Genomics Visium platform. Spatial cell type mapping was performed using single-cell transcriptome data as a reference [[92]35]. The reference data were preprocessed and modeled using the Cell2location Bayesian framework. We set the expected number of cells per location to 30 and applied empirical Bayes priors to avoid overfitting. The model was trained for 30,000 epochs with default hyperparameters, as the authors of Cell2location recommended. After training, we extracted the posterior cell abundance estimates and applied log-transformation for downstream visualization. To determine spatial enrichment, we computed Z-scores across tissue sections and considered a Z-score > 1.5 as indicative of significant enrichment for a given cell type in a spatial location. Computational gene perturbation Geneformer has been applied in various biological fields, including identifying potential disease drug targets [[93]36]. Geneformer effectively captures the relationships and importance of genes within single-cell transcriptomes. This focus mechanism allows Geneformer to prioritize the most relevant genes for specific learning objectives, optimizing prediction accuracy. In this study, we utilized Geneformer to investigate key genes associated with the transition from high EMT score MES cells to low EMT score MES cells in IDH wild-type recurrent patients. Survival analysis of transcriptomic cohorts We performed a prognostic analysis of MES stromal cell subpopulations. The R package GSVA was used to perform ssGSEA normalization scoring based on the gene sets of single-cell MES subpopulations within the expression matrix of tumor patients [[94]37]. Using the CGGA-693 cohort, survival analysis was conducted by grouping patients based on the median value of the ssGSEA scores. Immunostaining and analysis of human tissue microarrays A high-density tissue microarray containing 180 cases of gliomas of varying grades and control was purchased from Superchip Biotech (catalog number HBraG180Su01). An anti-EGR1 antibody (catalog number 55117-1-AP, Proteintech) was used for signal detection (1:500 dilution). Tissue spots with > 30% loss were excluded. The immunofluorescence (IF) score was calculated by multiplying the staining intensity by the percentage of stained cells Supplement Table [95]2. Cell culture, transfection and viral transduction The human glioma cell lines, U251 and T98G, were sourced from Procell (Wuhan, China) and cultured in Dulbecco’s Modified Eagle Medium (DMEM; Servicebio), supplemented with 10% fetal bovine serum (FBS, Excel) and 10 µL/mL penicillin-streptomycin solution (100×, Biosharp). The cells were incubated at 37 °C in a humidified atmosphere of 5% CO₂ (Hair). siRNA transfection For gene knockdown using siRNA, cells were seeded at 5 × 10⁴ cells per well in six-well plates and allowed to adhere for 24 h. siRNAs targeting EGR1 (siEGR1) and a non-targeting control (siControl), both obtained from Ribobio, were used. Transfections were performed using the Ribo FECT™ CP Transfection Kit following the manufacturer’s instructions when cells reached approximately 40% confluency. After six hours of incubation, the medium was replaced with fresh DMEM containing 10% FBS. Cells were collected 36 h post-transfection for subsequent experiments. Lentiviral transduction For stable SERPINE1 gene knockdown, cells were transduced with SERPINE1 shRNA lentivirus or control lentivirus (shControl) (Obio Technology). Polybrene (5 µg/ml, Yeasen) was used to enhance viral infection efficiency. After 24 h of infection, the medium was replaced with fresh DMEM containing 2.5 µg/ml puromycin (InvivoGen, USA) to select stably transduced cells for 10 days.siRNA and shRNA sequences used were listed in Supplement Table [96]3. RNA extraction, cDNA synthesis, and quantitative real-time PCR Total RNA was isolated from tissue samples stored at temperatures below − 80 °C, utilizing Takara RNAiso Plus (Takara Bio, Otsu, Shiga, Japan) according to the manufacturer’s instructions. For cDNA synthesis, 1 µg of RNA was reverse transcribed using the HiScript IV QRT SuperMix for qPCR kit (Vazyme Medical Technology). Quantitative real-time PCR (qPCR) data were analyzed using the 2⁻ΔΔCt method. Primer sequences are provided in Supplement Table [97]4. Cell counting Kit-8 proliferation assay CCK8 assay was utilized to gauge cell growth. U251 cells and T98G cells were respectively seeded in 96-well-plates at a density of 5,000 cells per 100 ul per well and treated for 24–72 h. Subsequently, 10ul of CCK8 reagent was added to each well, followed by a 30-minute incubation at 37 °C. Eventually, the Multiskan™ SkyHigh thermoscientific spectrophotometer was utilized to measure the absorbance at 450 nm. Wound healing assay The migration capacity of cells was gauged by performing wound healing assay. U251 cells and T98G cells were respectively seeded in 6-well-plates at a density of 500,000 cells per 1 ml per well. After 6 h, scratches were made using a 200-ul pipette tip. The cells were then washed three times with PBS and continued to be cultured in a medium containing 5% FBS. The healing progress of the wounds was monitored at 0 and 24 h. Flow cytometric cell cycle analysis To assess DNA fragmentation, U251 and T98G cells were plated in six-well plates at a density of 500,000 cells per 1 mL per well for 24 h. After incubation, the cells were collected by trypsinization, rinsed with a complete medium, washed with PBS, and fixed in 70% cold ethanol for 24 h. The cells were then gently resuspended in 500 µL of a hypotonic fluorochrome solution containing 450 µL propidium iodide (PI) and 50 µL RNase A to stain the nuclei. DNA content was measured using a CytoFLEX flow cytometer, and cell cycle distribution was analyzed with ModFit software, with 8,000 events recorded per sample. Animal experiments In a xenograft experiment, female athymic nude mice (n = 6/group) were acclimated for one week. U251 cells at 80% confluence were trypsinized after exposure to fresh medium for 24 h and inoculated 0.5 × 10^6 cells in 100 µl volumes into mice of each group. On day 25, tumor volume was measured by calipers and calculated as 1/2×(Length×Width×Height), and mice were euthanized for tumor fixation in formalin and paraffin embedding. Tissue sections were stained with H&E for proliferative analysis under a microscope. The animal experiments took place in the SPF Animal Laboratory at Wuhan University (ZN2023021). Statistical analysis Statistical analyses were performed using R (v4.3.1). For comparisons between two groups, unpaired two-tailed Student’s t-tests or Wilcoxon rank-sum tests were used as appropriate. For comparisons among more than two groups, one-way ANOVA with Tukey’s post hoc test was used for normally distributed data, and Kruskal–Wallis test for non-parametric data. Survival curves were analyzed using the Kaplan–Meier method and log-rank test. Multiple hypothesis testing was adjusted using the Benjamini–Hochberg method where applicable. A p-value < 0.05 was considered statistically significant. Results Distribution of cell subpopulations in IDH wild-type primary and recurrent gliomas The technical roadmap illustrates a multi-omics analysis of IDH wild-type primary and paired recurrent gliomas (Fig. [98]1A). We analyzed sequencing data from multiple cohorts to investigate the mechanisms of glioma recurrence. The findings were then validated through clinical tissue samples, cell-based experiments, and animal models. First, we analyzed transcriptional regulatory changes in cell subpopulations from single-cell data of primary and recurrent gliomas. Second, we identified key regulatory genes in specific MES-like subpopulations. Finally, we validated the roles of these key genes in gliomas through cellular and animal experiments. The UMAP plot shows the cell distribution of 52 IDH wild-type GBM samples, including 26 primary tumors and their 26 matched recurrent tumors (Supplement Fig. [99]1A). Cell subpopulations were annotated with markers for Oligodendrocytes, Astrocytes, Macrophages, Microglial cells, Neurons, Endothelial cells, and T cells (Fig. [100]1B). A total of 143,340 cells were included in the analysis for visualization (Fig. [101]1C). UMAP visualization depicts the distribution of cell subpopulations in primary and recurrent gliomas (Fig. [102]1D). The proportions of cell subpopulations in IDH wild-type primary gliomas were as follows: Astrocytes (60.5%), Oligodendrocytes (17.9%), Macrophages (7.3%), Microglial cells (5.1%), Neurons (7%), Endothelial cells (1.3%), and T cells (0.9%) (Fig. [103]1E). In recurrent gliomas, the proportions were higher for Oligodendrocytes (29.1%), Macrophages (7.9%), Microglial cells (6.4%), Neurons (7.6%), and Endothelial cells (4.1%) compared to primary gliomas. Enrichment analysis revealed that upregulated signaling pathways in recurrent gliomas include Salmonella infection, protein processing in the endoplasmic reticulum, endocytosis, Human T-cell leukemia virus 1 infection, and cell cycle pathways. Downregulated pathways include protein digestion and absorption, and cytoskeleton organization in muscle cells (Fig. [104]1F). Fig. 1. [105]Fig. 1 [106]Open in a new tab Comprehensive analysis of cell subpopulations and gene expression in gliomas using single-cell RNA sequencing (scRNA-seq). (A) Schematic workflow outlining the steps of single-cell RNA sequencing, data analysis, subpopulation classification, and experimental validation. (B) Heatmap displaying the expression levels of key genes across different glioma cell types, with rows representing genes and columns representing different cell subpopulations. (C) UMAP plot visualizes cell type distribution in the entire dataset, with each color representing a distinct cell type as annotated in the legend. (D) UMAP plots comparing cell populations between primary and recurrent gliomas, revealing differences in cell type composition between the two conditions. (E) Bar chart showing the relative proportions of different cell types in primary and recurrent gliomas, highlighting the higher abundance of certain cell types in recurrent gliomas. (F) Bar plot representing gene enrichment analysis, showing significantly upregulated and downregulated genes associated with glioma progression Distribution of tumor subpopulation expression and transcriptional activity To dissect the transcriptional heterogeneity of GBM at single-cell resolution, particularly across primary and recurrent states, we applied transcriptional signature analyses to identify distinct tumor subpopulations. We further identified tumor cell subpopulations using inferCNV. Violin plots illustrate the CNV scores of various cell clusters (Supplement Fig. [107]1B). Cell proportion plots revealed that primary gliomas comprised 86.4% tumor cells and 13.6% normal cells, while recurrent gliomas comprised 84.7% tumor cells and 15.3% normal cells (Supplement Fig. [108]1C). UMAP visualizations displayed the distribution of normal and tumor cells in both primary and recurrent gliomas (Supplement Fig. [109]1D). In the primary glioma, the number of malignant cells was 51,166, while the number of normal cells was 8,081. In recurrent glioma, the number of malignant cells increased to 71,221, and the number of normal cells was 12,872. Heatmaps showed the expression of marker genes in tumor cell subpopulations (Fig. [110]2A). UMAP visualizations revealed the spatial distribution of OC-like, OPC-like, AC-like, and MES-like tumor cell subpopulations (Fig. [111]2B). In the primary glioma, the numbers of the four types of malignant cells were as follows: OC-like (14,015), OPC-like (9,852), AC-like (26,709), and MES-like cells (590). In recurrent glioma, the numbers were: OC-like (28,222), OPC-like (11,632), AC-like (27,825), and MES-like cells (3,542). In primary gliomas, the proportions of OC-like, OPC-like, AC-like, and MES-like cells were 27.4%, 19.3%, 52.2%, and 1.2%, respectively (Supplement Fig. [112]1E). In recurrent gliomas, these proportions were 39.6%, 16.3%, 39.1%, and 5%, respectively. UMAP visualizations also showed gene expression and transcription factor activity in primary and recurrent tumor subpopulations (Fig. [113]2C). We conducted SCENIC analysis to reconstruct gene regulatory networks and identify cell states, generating nine modules in total (Fig. [114]2D). A total of 9 modules were formed through clustering. Modules 4 and 8 exhibited higher activity in recurrent glioma tissues (Fig. [115]2E). Our research primarily focused on mechanisms related to recurrent gliomas. Moreover, UMAP visualizations showed higher activity of MES subpopulations in Module 2 (Fig. [116]2F). CytoTRACE assessed the stemness and differentiation status of tumor cell subpopulations, with MES-like cells showing lower stemness scores (Fig. [117]2G, H). Fig. 2. [118]Fig. 2 [119]Open in a new tab Comparative analysis of primary and recurrent gliomas using single-cell RNA sequencing (scRNA-seq). (A) Heatmap showing the expression levels of key genes across different glioma subtypes in both primary and recurrent gliomas. Rows represent genes, and columns represent subtypes. (B) UMAP plots illustrate the distribution of 51,166 cells in primary gliomas and 71,221 cells in recurrent gliomas, with cells color-coded by subtype (OPC-like, AC-like, MES-like, OC-like). (C) UMAP plots comparing cell subtype distribution across different groups, showing changes in subtype composition between primary and recurrent gliomas. (D) Heatmap and transcription factor motif analysis highlighting key transcription factors associated with different glioma subtypes, including ETV5, JUN, and MYCN. (E-F) Expression patterns of various marker genes across glioma subtypes, visualized as UMAP feature plots, with intensity corresponding to gene expression levels. (G) CytoTRACE and phenotype plots showing differentiation states of cells across glioma subtypes, with predicted differentiation levels color-coded. (H) Boxplot representing the predicted ordering of cell differentiation levels by CytoTRACE across different glioma subtypes, indicating varying differentiation states between subtypes Enrichment pathway analysis and Spatial localization reveal characteristics of tumor subpopulations GSEA enrichment analysis of primary and recurrent glioma cell subpopulations indicated upregulated pathways in secondary MES-like subpopulation, including ECM-receptor interaction, WNT signaling, RAS signaling, and PI3K-AKT signaling pathways (Supplement Fig. [120]2A). This study supports previous findings on the significance of the interaction between the PI3K/AKT/mTOR and WNT/β-catenin pathways, suggesting that concurrently targeting these pathways may enhance the treatment of glioblastoma multiforme [[121]38, [122]39, [123]40]. In the AC-like subpopulation, upregulated pathways included inflammatory mediator regulation of TRP channels, oxytocin signaling, and cAMP signaling pathways. The OPC-like subpopulation showed upregulation in glioma, GABAergic synapse, transcriptional misregulation in cancer, MAPK signaling, and cGMP-PKG signaling pathway. The recurrent OC-like subpopulation had upregulated dopaminergic synapses and rap1 signaling pathways. Spatial transcriptomics analysis revealed the spatial localization of tumor cell subpopulations. The AC-like subpopulation occupied most of the spatial distribution in secondary tumors, consistent with our single-cell subpopulation analysis (Supplement Fig. [124]2B-F). Additionally, compared to primary gliomas, the spatial localization of MES-like cells was less prominent in secondary glioblastomas. Analysis of key signaling pathways, cell subpopulations, and EMT-related gene expression in recurrent gliomas We investigated the key signaling pathways and cell subpopulations involved in recurrent gliomas. Modules 4 and 8 were primarily expressed in the MES cell subpopulations of recurrent gliomas (Supplement Table [125]5). Additionally, transcriptional regulatory Module 2 was mainly expressed in the MES cell subpopulations of recurrent gliomas (Supplement Table [126]6). A Venn diagram of transcription factors from Modules 4, 8, and 2 in recurrent gliomas was constructed (Fig. [127]3A). Enrichment analysis revealed activation of the EMT pathway in the MES subpopulations, while IL2-STAT5 signaling, complement, and mitotic spindle pathways were inhibited (Fig. [128]3B). Key transcription factors activating the EMT pathway included EGR1, FLI1, and FOSL1 (Fig. [129]3C). An interaction map of EGR1 with its target genes was constructed (Fig. [130]3D). To further characterize intratumoral heterogeneity related to epithelial–mesenchymal transition (EMT), we performed unsupervised non-negative matrix factorization (NMF) clustering on malignant cells from both primary and recurrent glioblastoma samples, using the expression of EMT pathway-related genes. This approach enabled us to identify previously unrecognized EMT-like subpopulations, which may play distinct roles in tumor progression and treatment resistance. The heatmap displayed the expression of EMT-related genes in recurrent MES cell subpopulations (Fig. [131]3E). Compared to primary gliomas, the EGR1 regulon was highly expressed in the MES subpopulations of recurrent gliomas (p<0.001) (Fig. [132]3F). EGR1 gene, the EGR1 regulon, and the EMT pathway were spatially localized in secondary glioblastomas (Fig. [133]3G, H). Transcriptome sequencing validation showed high expression of EMT-related genes, including EGR1, FLI1, FOSL1, ERG, and EGR2, in recurrent gliomas (p<0.001) (Fig. [134]3I). Kaplan-Meier survival analysis indicated that high EGR1 expression was associated with poor prognosis in glioma patients (p<0.001) (Fig. [135]3J). The results also reveal a difference in tumor grading and the decision to undergo chemotherapy between patients with recurrent gliomas and those with primary gliomas (Supplement Table [136]7). Fig. 3. [137]Fig. 3 [138]Open in a new tab Gene expression analysis and survival prediction in glioblastoma subtypes (A) Venn diagram showing the overlap of differentially expressed genes between primary and recurrent glioblastomas, with the bar plot below illustrating the number of genes unique to each group. (B) Pathway enrichment analysis highlighting the activation and suppression of key signaling pathways, including IL2/STAT5 signaling, complement pathway, epithelial-mesenchymal transition (EMT), and mitotic spindle pathway, with dot size representing gene ratio and color indicating p-values. (C) Dot plot depicting the differential expression of transcription factors (EGR1, FLI1, FOXL1, etc.) between primary and recurrent glioblastomas. The dot size corresponds to the gene ratio, and the color indicates p-values. (D) Interaction network of transcription factors associated with glioblastoma subtypes, centered around EGR1 and related factors. (E) Heatmap showing the expression of key genes across glioblastoma subtypes (OC-like, OPC-like, AC-like, MES-like) in both primary and recurrent tumors. (F) Boxplot of EGR1 expression levels across glioblastoma subtypes, demonstrating elevated expression in MES-like cells, particularly in recurrent glioblastomas. (G) Spatial transcriptomic heatmaps showing the spatial distribution of EGR1 and EMT marker expression in secondary glioblastomas, with color intensity corresponding to expression levels. (H) Boxplot comparing the expression levels of various genes between primary and recurrent glioblastomas, showing significant differences in key regulators. (I) Kaplan-Meier survival curve showing the association between EGR1 expression levels and patient survival, with high EGR1 expression correlating with poorer prognosis (p < 0.001) Identification of novel EMT-like subpopulations and analysis of their signaling pathways and stemness differentiation We identified novel EMT-like subpopulations through unsupervised NMF clustering based on EMT pathway-related genes [[139]41]. The identified subpopulations included MGP + MES-like-C1 (15.7%), SLIT2 + MES-like-C2 (20.4%), SNTB1 + MES-like-C3 (7.9%), DPYSL3 + MES-like-C4 (7.1%), VEGFC + MES-like-C5 (6.8%), SERPINE1 + MES-like-C6 (10.4%), COL12A1 + MES-like-C7 (5.3%), unclear-MES-like-C8 (17.9%), and nongene-MES-like-C9 (8.4%) (Fig. [140]4A, B). Moreover, glioma patients with high expression of EMT-related genes, including MGP, SERPINE1, SNTB1, VEGFC, and COL12A1, had poor prognosis (Supplement Fig. [141]3A-G). Specifically, glioma patients with high expression of MGP + MES-like-C1, VEGFC + MES-like-C5, SERPINE1 + MES-like-C6, and COL12A1 + MES-like-C7 subpopulation cells had poor prognosis (Supplement Fig. [142]4A-G). The signaling pathways enriched in different EMT-like subpopulations included vascular smooth muscle contraction, ECM-receptor interaction, and focal adhesion pathways. Specifically, the SERPINE1 + MES-like-C6 and COL12A1 + MES-like-C7 subpopulations were enriched in the ECM-receptor interaction pathway (Fig. [143]4C). We assessed the stemness differentiation status of these tumor subpopulations using CytoTRACE (Fig. [144]4D). Box plots indicated that the MGP + MES-like-C1 subpopulation had the highest stemness score (Fig. [145]4E). Spatial transcriptomics analysis of recurrent gliomas revealed that most subpopulations were located at the tumor margins, with a broader distribution observed for the DPYSL3 + MES-like-C4 subpopulation (Fig. [146]4F, G). Using the Geneformer deep learning model, we inferred the key genes influencing the MES-like tumor subpopulations. The key genes identified were MGP, SLIT2, SNTB1, DPYSL3, VEGFC, SERPINE1, and COL12A1, which are associated with the transition from high EMT-score MES cells to low EMT-score MES cells in IDH wild-type recurrent gliomas (Fig. [147]4H, I). Fig. 4. [148]Fig. 4 [149]Open in a new tab Subtype analysis and pathway enrichment in secondary glioblastomas. (A) UMAP plot showing the distribution of secondary glioblastoma cells, color-coded by non-negative matrix factorization (NMF) subtypes. Each cluster represents a distinct subtype, annotated with marker genes. (B) Bar chart representing the relative proportions of different NMF-derived subtypes across secondary glioblastomas. Subtype MES-like C7 shows a higher prevalence compared to other subtypes. (C) Dot plot illustrating the pathway enrichment analysis for different subtypes. Pathways such as ECM-receptor interaction, focal adhesion, and protein digestion and absorption are highlighted, with dot size indicating the number of enriched genes and color representing significance (p-value). (D) CytoTRACE and phenotype UMAP plots showing differentiation states of the glioblastoma cells. Predicted differentiation states are color-coded, with less differentiated cells clustered in MES-like subtypes. (E) Boxplot comparing CytoTRACE-predicted ordering scores for glioblastoma subtypes, showing variation in predicted differentiation levels across subtypes. (F-G) Histological sections of secondary glioblastoma corresponding to each MES-like subtype. Color intensity reflects expression levels, highlighting the spatial heterogeneity of gene expression across tumor sections. (H-I) Volcano plots showing shifts in EMT (epithelial-mesenchymal transition) scores and differential gene expression in secondary glioblastomas. Key genes such as SERPINE1, MOAP1, and COLEC12 are highlighted based on their significance and EMT shift The trajectory analysis and changes in cell communication of MES-like tumor subpopulations were examined To explore the state changes of novel EMT-related MES-like tumor subpopulations, we conducted cell trajectory and cell communication analyses. The pseudotime analysis revealed that MES-like tumor cells could be classified into three states (Supplement Fig. [150]5A). Furthermore, the novel EMT-related MES-like tumor subpopulations formed the cell differentiation trajectory, with MGP and SERPINE1 located at the terminal stage of cell differentiation (Supplement Fig. [151]5B, C). Subsequently, the network diagram illustrates the number and strength of interactions between cell subpopulations (Supplement Fig. [152]5D, E), with MGP + MES-like-C1 showing a relatively high number and strength of interactions. Additionally, the heatmap indicated that both MGP + MES-like-C1 and COL12A1 + MES-like-C7 had prominent input and output signals (Supplement Fig. [153]5F). MGP + MES-like-C1 exhibited outgoing signaling pathways such as MK, BMP, IGF, and SEMA3, while incoming pathways included MK, BMP, PTN, EDN, and VISFATIN. COL12A1 + MES-like-C7 showed outgoing signals in the MK, BMP, IGF, SEMA3, and HGF pathways, with incoming signals from MK, BMP, PTN, IGF, and EDN. A hierarchy plot further demonstrated that the MGP + MES-like-C1 subpopulation interacted more frequently with other subpopulations through the MK and IGF signaling pathways (Supplement Fig. [154]5G, H). Lastly, the bubble plot displayed ligand-receptor relationships between subpopulations, highlighting that MDK, IGF, and BMP6 signaling pathways were involved in the regulation of MGP + MES-like-C1 and SERPINE1 + MES-like-C6 over the other subpopulations (Supplement Fig. [155]5I). Cell lineage differentiation and gene correlation analysis of MES-like tumor subpopulations We inferred the cell lineage differentiation of MES-like tumor subpopulations using Slingshot analysis, which identified two lineages (Fig. [156]5A, B). The MGP + MES-like-C1 and SERPINE1 + MES-like-C6 subpopulations were located at the end of the differentiation trajectories (Fig. [157]5C, D). Dynamic analysis revealed six cell clusters that were enriched to analyze biological processes, molecular functions, and cellular components (Fig. [158]5D). The correlation between EGR1 and MGP was high in both recurrent and primary gliomas (R[R]=0.45, R[P]=0.60) (Fig. [159]5E). Similarly, the correlation between EGR1 and SERPINE1 was high in both recurrent and primary gliomas (R[R]=0.54, R[P]=0.60) (Fig. [160]5F). We subsequently analyzed the immunofluorescence expression of EGR1 and SERPINE1 in primary and recurrent gliomas using tissue microarrays (Fig. [161]6A, B). Immunofluorescence experiments on tissue microarrays of primary and recurrent gliomas showed a high correlation between EGR1 and SERPINE1 (R = 0.47) (Fig. [162]6C). EGR1 exhibited high correlation in both primary and recurrent gliomas (Fig. [163]6D). Transcriptome sequencing analysis revealed a positive correlation between the expression of ERPINE1 and EGR1 (p<0.001) (Fig. [164]6E, F). Fig. 5. [165]Fig. 5 [166]Open in a new tab Lineage tracing and molecular characterization of glioblastoma subtypes (A) UMAP plot showing the distribution of secondary glioblastoma cells, color-coded by NMF-derived subtypes. Overlaid lines represent two distinct cell lineages identified within the dataset. (B) UMAP plots showing the two main lineages, Lineage 1 (left) and Lineage 2 (right), with pseudotime progression color-coded to indicate cellular differentiation along each lineage. (C) Expression trends of key genes (MGP, SLIT2, SNTB1, DYSF, VEGFC, SERPINE1, COL12A1) along the pseudotime axis in each lineage. The solid lines represent the smoothed gene expression curves for each gene, indicating dynamic expression changes during glioblastoma progression. (D) Heatmap displaying the gene expression profiles across subtypes in both lineages, with corresponding Gene Ontology (GO) enrichment analysis for biological processes (BP), cellular components (CC), and molecular functions (MF) shown on the right. The vertical lines divide the NMF subtypes into their respective lineages. (E-F) Scatter plots showing the correlation between SERPINE1 and MCP expression, highlighting significant correlations within the two main lineages. The color intensity reflects the strength of the correlation, with red lines representing the linear regression fits Fig. 6. [167]Fig. 6 [168]Open in a new tab Correlation between EGR1 and SERPINE1 expression in glioblastoma and functional validation. (A) Tissue microarray showing immunofluorescence staining for various markers across multiple glioblastoma samples. (B) Immunofluorescence images displaying EGR1 (green) and SERPINE1 (red) expression in primary and recurrent glioblastoma tissues. DAPI (blue) marks the nuclei, and the merged images show the co-localization of EGR1 and SERPINE1. (C) Correlation matrix showing relationships between EGR1, SERPINE1, and clinical factors, including gender, age, grade, laterality, progression-free survival (PFS), and overall survival (OS). A strong correlation between EGR1 and SERPINE1 (R = 0.47) is observed. (D) Bar plot showing Pearson’s correlation coefficient (R) between EGR1 and SERPINE1 in primary and recurrent glioblastomas, with a significant difference (*p < 0.05) between the two groups. (E) The relationship between SERPINE1 and EGR1 in the transcriptome sequencing of primary gliomas. (F) The relationship between SERPINE1 and EGR1 in the transcriptome sequencing of recurrent gliomas. The red regression line in the figure represents linear fitting, and the gray shaded area indicates the 95% confidence interval Expression and enrichment analysis of EMT-related genes, immune characteristics, and metabolic pathways in MES-like tumor subpopulations We examined the expression of EMT-related genes in various EMT-associated MES-like subpopulations (Supplement Fig. [169]6A). Next, we analyzed the expression of chemokines, chemokine receptors, MHC-related molecules, immune characteristics, and interferons in these subpopulations (Supplement Fig. [170]6B). Chemokines, chemokine receptors, and interferon-related genes were primarily expressed in the SERPINE1 + MES-like-C6 and VEGFC + MES-like-C5 subpopulations. The MGP + MES-like-C1 subpopulation was enriched in metabolic pathways such as sphingolipid metabolism, primary bile acid biosynthesis, and cysteine and methionine metabolism (Supplement Fig. [171]6C). The SERPINE1 + MES-like-C6 subpopulation was mainly enriched in steroid hormone biosynthesis, nicotinate and nicotinamide metabolism, nitrogen metabolism, and N-glycan biosynthesis. The DPYSL3 + MES-like-C4 and VEGFC + MES-like-C5 subpopulations enriched most metabolic pathways. Furthermore, we used scFEA to infer the metabolic flux of cell subpopulations through a graph neural network model (Supplement Fig. [172]6D). The fatty acid signaling pathway was upregulated in the MGP + MES-like-C1 and COL12A1 + MES-like-C7 subpopulations. The VEGFC + MES-like-C5 subpopulation showed downregulation of glucose and other metabolic signaling pathways. Expression and differentiation of epithelial-mesenchymal transition (EMT) related genes in IDH wild-type gliomas at different stages Given that our analysis of glioma single cells used nuclear sequencing, which contains fewer immune cells, we next examined the changes in the tumor immune microenvironment in IDH wild-type glioma mice at early and late stages. After annotating the cell subpopulation markers, the tumor microenvironment subpopulations were identified as microglial cells, monocytes, macrophages, T cells, NK cells, dendritic cells (DCs), B cells, and oligodendrocytes (Supplementary Fig. [173]7A). Compared to IDH-mutant glioma subgroups, IDH-wildtype gliomas exhibit higher proportions of monocytes (57.8%), macrophages (66.4%), T cells (65.5%), NK cells (74.2%), and dendritic cells (DCs) (65.1%) (Supplementary Fig. [174]7B, C). UMAP visualizations showed the distribution of the tumor microenvironment in IDH wild-type glioma mice at early and late stages (Supplementary Fig. [175]7D). Compared to early-stage IDH wild-type mice, the expression of Egr1 was significantly lower in immune cells within the tumor microenvironment of late-stage IDH wild-type mice (p < 0.001) (Supplementary Fig. [176]7E). Additionally, the CytoTRACE2 stemness analysis revealed that late-stage IDH wild-type cells exhibited a higher degree of stemness differentiation (Supplementary Fig. [177]7F-J). The pseudo-temporal analysis revealed that advanced-stage gliomas predominantly exist in the terminal differentiation phase (Supplementary Fig. [178]8A). The differential expression of Egr1 and Serpine1 also varied across different stages of glioma progression (Supplementary Fig. [179]8B, C). SERPINE1 Inhibition reduced glioma proliferation and migration To further explore the oncogenic role of SERPINE1 in gliomas, we designed siRNA sequences targeting SERPINE1 and assessed their knockdown efficiency in glioma cells for subsequent experiments (Fig. [180]7A). Previous studies have demonstrated that silencing EGR1 expression significantly inhibits glioma cell proliferation and induces G1 phase cell cycle arrest, indicating that EGR1 plays a crucial role in glioma cell proliferation [[181]42, [182]43]. The inhibition of SERPINE1 was found to diminish the proliferation and migration of glioma cells (Fig. [183]7B, C). The siRNA targeting SERPINE1 was designed and its efficacy in reducing SERPINE1 levels in U251 and T98G cells was confirmed. The CCK8 assay indicated a significant reduction in OD450 values in U251 cells and T98G cells 72 h after transfection, suggesting a decrease in cell proliferation (Fig. [184]7B). The wound healing assay showed a reduced healing rate in SERPINE1-depleted cells at 24 h, indicating impaired migration of glioma cells in vitro (Fig. [185]7C). Fig. 7. [186]Fig. 7 [187]Open in a new tab Functional validation of SERPINE1 knockdown in glioblastoma cells and in vivo tumor models. (A) Quantitative PCR results showing a significant reduction in SERPINE1 mRNA expression after si-SERPINE1 treatment in both T98G and U251 glioblastoma cells (****p < 0.0001). (B) Proliferation curves of T98G and U251 cells treated with si-SERPINE1, demonstrating reduced cell proliferation over 72 h (**p < 0.01). (C) Wound healing assay showing reduced migration distance in si-SERPINE1 treated T98G and U251 cells at 24 h. Quantification of migration distance shows a significant reduction compared to the negative control (NC) (****p < 0.0001). (D) Flow cytometry analysis of cell cycle distribution in si-SERPINE1 and NC-treated T98G and U251 cells. si-SERPINE1 treated cells show an increased percentage of cells in the G1 phase and a reduction in S phase (**p < 0.01, ***p < 0.001). (E) In vivo tumor growth in nude mice. Representative images of mice and excised tumors are shown for the NC and sh-SERPINE1 treated groups. (F) Tumor weight of excised tumors from sh-SERPINE1 and NC-treated mice. Tumor weight is significantly lower in the sh-SERPINE1 group (**p < 0.01). (G) Tumor volume measurement showing reduced tumor size in sh-SERPINE1 treated mice compared to NC-treated mice (*p < 0.05). (H) Histological analysis of tumor sections stained with H&E, Ki67, and SERPINE1. The sh-SERPINE1 treated group shows reduced Ki67 and SERPINE1 staining, indicating decreased proliferation and target knockdown efficacy Cell cycle analysis The cell cycle distribution of U251 cells and T98G cells treated with the siRNA sequence against SERPINE1 was measured by comparing the proportion of cells in the G0/G1, S, and G2/M phases. After treatment with the siRNA sequence against SERPINE1, the proportion of cells in G1 and S was significantly increased, while the proportion in the G2/M phase was decreased (Fig. [188]7D). These results indicate that the siRNA sequence against SERPINE1 facilitated the arrest of cells in the G1 and S phases of the cell cycle, thereby inhibiting U251 cell and T98G cell proliferation. SERPINE1 knockdown inhibits glioma cell proliferation in vivo Additionally, we examined the effects of SERPINE1 in vivo (Fig. [189]7E). The results from the tumor xenograft studies demonstrated that sh-SERPINE1 significantly reduced the size of the transplanted tumors, as evidenced by decreases in both tumor volume and weight (Fig. [190]7F, G). These findings are consistent with our in vitro results. In the xenograft tumor model, a comparison between the si-NC and si-SERPINE1 groups revealed that, while both groups formed tumor tissue based on HE staining, the tumor volume in the si-SERPINE1 group was significantly smaller (Fig. [191]7H). IHC staining for Ki67 indicated that the proportion of Ki67-positive cells, a marker of proliferation, was notably lower in the si-SERPINE1 group compared to the si-NC group, suggesting a marked reduction in tumor cell proliferation in the si-SERPINE1 group. Additionally, the expression level of SERPINE1 was significantly decreased in the si-SERPINE1 group, further supporting its potential role in inhibiting tumor growth. Overall, the results demonstrated that si-SERPINE1 exhibited a pronounced tumor-suppressive effect in the xenograft model, primarily reflected in reduced tumor proliferative activity and decreased SERPINE1 expression. Discussion In this study, we identified key regulatory genes and cellular mechanisms driving the recurrence of IDH wild-type gliomas. Through extensive single-cell RNA sequencing and spatial transcriptomics, we found that the transcription factor EGR1 significantly influences the EMT in glioma cells. Our analysis revealed that MES-like cell subpopulations exhibit high transcriptional regulatory activity, contributing to glioma invasion and drug resistance. Furthermore, we discovered that the SERPINE1-MES subpopulation is a critical factor in glioma recurrence, as evidenced by its high expression in recurrent gliomas and its association with poor patient prognosis. These findings were validated through functional assays and animal models, demonstrating the pivotal role of EGR1 and SERPINE1 in glioma recurrence and highlighting potential therapeutic targets to improve patient outcomes. The recurrence of gliomas, particularly IDH wild-type gliomas, is a significant clinical challenge. Our findings that EGR1 mediates EMT in MES-like subpopulations align with previous studies that have implicated EMT in tumor progression and metastasis [[192]4, [193]6, [194]44, [195]45, [196]46]. EGR1’s role in enhancing the invasive capabilities of glioma cells supports the notion that transcription factors are crucial in regulating the plasticity and adaptability of cancer cells under therapeutic pressure [[197]5, [198]42]. The upregulation of EMT-related genes in gliomas suggests that these factors could serve as biomarkers for aggressive glioma phenotypes and potential targets for novel therapeutic interventions [[199]7, [200]8, [201]42, [202]47, [203]48]. Our spatial transcriptomics analysis revealed that MES-like subpopulations are primarily located at the tumor margins in recurrent gliomas, indicating their role in tumor invasion and recurrence [[204]20]. This spatial distribution pattern is consistent with the enhanced invasive properties of these cells, as they interact with and degrade the extracellular matrix to facilitate tumor spread [[205]13, [206]25]. The identification of specific subpopulations, such as the SERPINE1-MES-like cells, which are enriched in ECM-receptor interaction pathways, further underscores the importance of the tumor microenvironment in glioma progression [[207]12]. These findings are corroborated by previous research demonstrating the critical role of the extracellular matrix in providing structural support and biochemical signals that drive tumorigenesis [[208]10]. The use of single-cell RNA sequencing allowed us to uncover the heterogeneity within glioma cell populations and identify subpopulations with distinct transcriptional profiles. This approach provides a more nuanced understanding of the tumor’s cellular landscape compared to bulk RNA sequencing [[209]17]. Our study revealed that recurrent gliomas exhibit higher proportions of oligodendrocyte-like and MES-like cells compared to primary gliomas, suggesting a shift in cellular composition that may contribute to therapy resistance and tumor aggressiveness [[210]16]. The differentiation potential of these subpopulations, as assessed by CytoTRACE, indicates that MES-like cells possess lower stemness scores, aligning with their specialized role in invasion rather than self-renewal [[211]33]. Functional assays further validated the roles of EGR1 and SERPINE1 in glioma recurrence. Knockdown of EGR1 in glioma cell lines led to a significant reduction in SERPINE1 expression, decreased cell proliferation, and impaired migratory capabilities, highlighting the interdependence of these factors in promoting tumor progression [[212]26]. In vivo experiments using xenograft models confirmed these findings, with SERPINE1 knockdown resulting in smaller tumor sizes and reduced invasiveness. These results suggest that targeting the EGR1-SERPINE1 axis could be a promising strategy for inhibiting glioma recurrence and improving patient outcomes [[213]21]. Despite these significant findings, our study has several limitations. First, while single-cell RNA sequencing provides detailed insights into cellular heterogeneity, it does not capture the full complexity of tumor-microenvironment interactions. Additional studies incorporating other omics approaches, such as proteomics and metabolomics, could provide a more comprehensive understanding of the mechanisms driving glioma recurrence. Second, our findings are based on samples from a relatively small cohort of patients, which may limit the generalizability of our results. Larger studies with diverse patient populations are needed to validate these findings. Third, while our in vivo models provided valuable insights into the roles of EGR1 and SERPINE1, they do not fully recapitulate the human tumor microenvironment. Future research should explore the use of patient-derived xenografts and organoid models to better mimic human glioma biology. While our functional validation experiments using T98 and U251 glioblastoma cell lines provided initial evidence supporting a regulatory relationship between EGR1 and SERPINE1, we acknowledge that these established cell lines do not fully recapitulate the molecular and phenotypic complexity of recurrent GBM. Specifically, they lack the therapy-induced evolutionary pressures and tumor microenvironmental context present in patient-derived recurrent tumors. As such, the functional implications of the EGR1–SERPINE1 axis in recurrence remain to be fully elucidated. Future studies utilizing patient-derived glioma stem-like cells, tumorsphere cultures, or recurrent GBM xenograft models will be essential to validate and extend these findings in a clinically relevant context. Conclusion Our study provides a comprehensive analysis of the regulatory mechanisms underlying the recurrence of IDH wild-type gliomas, emphasizing the critical role of the transcription factor EGR1 and the SERPINE1-MES subpopulation in mediating EMT and tumor invasion. This research highlights the complexity and heterogeneity of glioma cell subpopulations, revealing novel insights into their dynamic plasticity and interaction with the tumor microenvironment. These findings underscore the importance of targeting specific transcriptional regulators and signaling pathways to develop more effective therapeutic strategies for glioma treatment. Electronic supplementary material Below is the link to the electronic supplementary material. [214]40478_2025_2036_MOESM1_ESM.tif^ (8.6MB, tif) Supplementary Material 1: Supplement Figure 1. Copy number variation (CNV) analysis in primary and recurrent glioblastoma. (A) UMAP visualization of 143,340 cells from glioma samples, showing distinct cell populations, including oligodendrocytes, astrocytes, macrophages, microglia, endothelial cells, neurons, and T cells, color-coded based on cell type. (B) Violin plots showing CNV scores across different Seurat clusters. Each cluster represents a distinct cell population, and CNV scores are color-coded based on cluster identity. Higher CNV scores suggest a greater degree of malignancy. (C) Bar plot comparing the proportion of malignant and normal cells between primary and recurrent glioblastoma samples. Primary tumors consist of 86.4% malignant cells, while recurrent tumors contain 84.7% malignant cells. (D) UMAP plots visualizing the distribution of malignant (red) and normal (blue) cells in primary (left) and recurrent (right) glioblastomas. Malignant cells dominate both primary and recurrent samples, but recurrent tumors show a slightly higher proportion of normal cells. (E) Bar chart showing the relative proportions of cell subtypes in primary and recurrent gliomas. Recurrent gliomas display a higher percentage of MES-like and OC-like cells compared to primary gliomas. [215]40478_2025_2036_MOESM2_ESM.tif^ (1.5MB, tif) Supplementary Material 2: Supplement Figure 2. Spatial and pathway analysis of primary and secondary glioblastomas using spatial transcriptomics. (A) Pathway enrichment analysis showing significantly altered pathways across glioblastoma subtypes, including inflammatory mediator regulation of TRP channels, Wnt signaling pathway, ECM-receptor interaction, and MAPK signaling pathway. The x-axis represents -log10(q-value) and the dot size indicates the number of altered genes in each pathway. (B-C) Spatial transcriptomic heatmaps showing gene expression patterns for AC-like, MES-like, OC-like, and OPC-like cells in secondary glioblastoma samples. Each plot represents the spatial distribution of gene expression across the tumor sections, with color intensity corresponding to expression levels. (D-F) Spatial transcriptomic heatmaps illustrating gene expression patterns for AC-like, MES-like, OC-like, and OPC-like cells in primary glioblastoma samples. Similar to secondary glioblastoma, these heatmaps highlight the spatial heterogeneity of gene expression in the different subtypes. [216]40478_2025_2036_MOESM3_ESM.tif^ (5.7MB, tif) Supplementary Material 3: Supplement Figure 3. Kaplan-Meier survival curves for glioblastoma patients based on the expression levels of key genes. (A) Kaplan-Meier survival curve for MGP expression. Patients with high MGP expression show significantly worse survival (p < 0.001). (B) Kaplan-Meier survival curve for SERPINE1 expression. Higher SERPINE1 expression is associated with poorer survival (p < 0.001). (C) Kaplan-Meier survival curve for SLIT2 expression, showing no significant difference between high and low expression groups (p = 0.827). (D) Kaplan-Meier survival curve for SNTB1 expression. High SNTB1 expression correlates with worse survival (p = 0.010). (E) Kaplan-Meier survival curve for VEGFC expression. Patients with high VEGFC expression have significantly worse survival outcomes (p = 0.035). (F) Kaplan-Meier survival curve for COL12A1 expression. Higher COL12A1 expression is associated with significantly reduced survival (p < 0.001). (G) Kaplan-Meier survival curve for DPYSL3 expression. There is no significant difference in survival between high and low DPYSL3 expression groups (p = 0.155). [217]40478_2025_2036_MOESM4_ESM.tif^ (5.2MB, tif) Supplementary Material 4: Supplement Figure 4. Kaplan-Meier survival curves for glioblastoma patients based on gene expression in MES-like subtypes. (A) Kaplan-Meier survival curve for DPYSL3+MES-like C4 subtype. No significant difference in survival is observed between high and low expression groups (p = 0.722). (B) Kaplan-Meier survival curve for MGP+MES-like C1 subtype. High expression is significantly associated with worse survival (p < 0.001). (C) Kaplan-Meier survival curve for SERPINE1+MES-like C6 subtype. Patients with high expression have significantly poorer survival outcomes (p < 0.001). (D) Kaplan-Meier survival curve for SLIT2+MES-like C2 subtype. No significant difference in survival is observed between high and low expression groups (p = 0.231). (E) Kaplan-Meier survival curve for SNTB1+MES-like C3 subtype. No significant difference in survival is observed (p = 0.149). (F) Kaplan-Meier survival curve for VEGFC+MES-like C5 subtype. High expression is associated with significantly worse survival (p = 0.001). (G) Kaplan-Meier survival curve for COL12A1+MES-like C7 subtype. High expression is significantly associated with poorer survival (p < 0.001). [218]40478_2025_2036_MOESM5_ESM.tif^ (13.3MB, tif) Supplementary Material 5: Supplement Figure 5. Pseudotime analysis and signaling pathway interactions in glioblastoma subtypes. (A) Monocle2 trajectory plots showing the pseudotime progression of 3,536 cells across various NMF-derived subtypes. Each point represents a single cell, color-coded by subtype and pseudotime state. (B) Expression dynamics of MGP along the pseudotime trajectory. Cells with higher MGP expression cluster at the end of the trajectory, indicating association with later differentiation stages. (C) Expression dynamics of SERPINE1 along the pseudotime trajectory, showing similar clustering at later pseudotime points. (D) Cell-cell communication network analysis showing signaling interactions between different MES-like subtypes. The size of the nodes reflects the strength of outgoing signals, and the thickness of edges indicates interaction strength between subtypes. (E) Detailed signaling interaction network of MES-like subtypes, focusing on specific pathways, including those involving MGP, VEGFC, and SERPINE1. (F) Heatmap showing outgoing and incoming signaling patterns between glioblastoma subtypes. The intensity of green color indicates the strength of signaling interactions between subtypes. (G) MK signaling pathway network showing specific source-target relationships between subtypes, with key signals originating from MGP and targeting various MES-like subtypes. (H) IGF signaling pathway network highlighting interactions between MES-like subtypes, with MGP as a central node. (I) Dot plot showing pathway enrichment analysis for key transcription factors involved in glioblastoma subtype differentiation. Dot size represents the number of enriched pathways, and color reflects p-value significance. [219]40478_2025_2036_MOESM6_ESM.tif^ (8.6MB, tif) Supplementary Material 6: Supplement Figure 6. Metabolic pathway analysis and gene expression profiling across glioblastoma subtypes. (A) Heatmap showing the expression of key metabolic genes across NMF-derived glioblastoma subtypes. Rows represent genes, and columns represent subtypes, with color intensity reflecting normalized expression levels (z-scores). (B) Heatmap illustrating the expression of genes involved in the glycosylation pathway across glioblastoma subtypes, with similar color-coding to indicate gene expression differences between subtypes. (C) Dot plot depicting the enrichment of various metabolic pathways across NMF-derived subtypes. Dot size corresponds to the number of genes involved in each pathway, and color represents the significance of enrichment (p-value). (D) Pathway-specific heatmap visualizing key metabolic reactions (e.g., glucose to glucose-6-phosphate, fatty acid metabolism, and methionine to cysteine conversion) across glioblastoma subtypes, with each subtype displaying a distinct metabolic profile. Rows represent pathways, and columns represent subtypes, with color intensity reflecting pathway activity levels. [220]40478_2025_2036_MOESM7_ESM.tif^ (10.7MB, tif) Supplementary Material 7: Supplement Figure 7. Single-cell RNA sequencing analysis of immune cell populations and differentiation states in IDH wild-type (WT) and mutant (MUT) glioblastoma subtypes. (A) Heatmap showing differential gene expression across various immune cell subtypes in IDH WT and MUT glioblastomas. Rows represent genes, and columns represent cell subtypes. (B) Bar chart illustrating the relative proportions of immune cell types in WT and MUT glioblastomas. (C) Stacked bar chart showing the distribution of specific immune cell subtypes in WT and MUT glioblastomas, highlighting significant differences in monocyte and microglia populations. (D) UMAP plots visualizing the distribution of late-stage (left) and early-stage (right) glioblastoma cells across different immune cell types, including microglia, monocytes, macrophages, T cells, NK cells, dendritic cells (DCs), and oligodendrocytes. (E) Violin plots comparing Egr1 expression levels across immune cell subtypes in WT and RII glioblastomas, with significant differences marked (****p < 0.0001). (F) UMAP plot showing cell phenotypes across WT and RII glioblastoma samples, color-coded based on phenotype classifications. (G) UMAP plot showing CytoTRACE-predicted differentiation potential, with less differentiated cells clustered on the left and more differentiated cells on the right. (H) UMAP plot showing CytoTRACE-derived differentiation categories, dividing cells into undifferentiated and differentiated groups. (I) Boxplot comparing developmental potential between phenotypic groups, with stemness scores significantly higher in the WT group. (J) UMAP plot showing relative differentiation states as predicted by CytoTRACE, with color intensity indicating levels of differentiation. [221]40478_2025_2036_MOESM8_ESM.tif^ (1.4MB, tif) Supplementary Material 8: Supplement Figure 8. Pseudotime trajectory and gene expression dynamics in glioblastoma. (A) Monocle2 trajectory plots depicting the pseudotime progression of 38,931 cells across various states. Each point represents a cell, with colors indicating distinct pseudotime states (left), timepoints (middle), and overall pseudotime progression (right). (B) Pseudotime trajectory highlighting the expression of Egr1. Cells with higher Egr1 expression are predominantly located at specific branches of the trajectory, suggesting its association with certain differentiation stages. (C) Pseudotime trajectory showing the expression of Serpine1. Similar to Egr1, higher Serpine1 expression clusters at particular branches, indicating its role in specific cell states along the differentiation trajectory. [222]Supplementary Material 9^ (11.6KB, xlsx) [223]Supplementary Material 10^ (24.6KB, xlsx) [224]Supplementary Material 11^ (10.7KB, xlsx) [225]Supplementary Material 12^ (9.4KB, xlsx) [226]Supplementary Material 13^ (10KB, xlsx) [227]Supplementary Material 14^ (10.4KB, xlsx) [228]Supplementary Material 15^ (20.7KB, docx) Acknowledgements