- 1Department of Neurosurgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou University, Zhengzhou, Henan, China
- 2College of First Clinical Medicine, Shandong University of Traditional Chinese Medicine, Jinan, Shandong, China
- 3Department of Scientific Research, The First Affiliated Hospital of Zhengzhou University, Zhengzhou University, Zhengzhou, Henan, China
Background: Glioblastoma (GBM) was considered the most aggressive type of primary brain tumor, marked by poor clinical outcomes and a high tendency to relapse. The therapeutic efficacy of GBM was significantly compromised by tumor heterogeneity, dysregulated metabolic pathways, the formation of an immunosuppressive microenvironment, and treatment resistance. Therefore, multi-dimensional therapeutic strategies targeting GBM-specific molecular features, its intrinsic properties, and microenvironmental regulatory networks were considered to potentially provide new breakthroughs for overcoming treatment resistance in GBM.
Methods: We analyzed single-cell RNA sequencing (scRNA-seq) data processed with the Seurat package to accurately identify cell types. Spatial transcriptomics integrated Multimodal Intersection Analysis, TransferData, and Robust Cell Type Decomposition techniques to characterize the spatial distribution patterns of key cell subtypes. CellChat was employed to assess intercellular communication networks. Furthermore, in vitro experiments confirmed the main regulatory role of YEATS4 (key transcription factor of C2 PCLAF+ subtype) in GBM malignant progression.
Results: Through scRNA-seq, we identified the C2 PCLAF+ subtype in GBM and analyzed its molecular characteristics and functional role in tumor progression. This subtype exhibited a unique malignant phenotype, marked by significant proliferative activity, characteristic metabolic reprogramming, and dysregulated cell death regulation mechanisms. Spatial transcriptomics revealed its preferential localization within specific tumor niches. Furthermore, the C2 PCLAF+ subtype established a specific interaction with fibroblasts through the MDK-LRP1 ligand–receptor pair. Critically, silencing YEATS4 in vitro significantly inhibited GBM malignancy. Additionally, the prognostic risk score model based on the C2 PCLAF+ subtype demonstrated significant clinical translational value.
Conclusion: Our study systematically elucidated the malignant characteristics of the C2 PCLAF+ subtype and its molecular mechanisms driving GBM progression. This subtype promoted therapeutic resistance through unique metabolic reprogramming, MDK-LRP1-mediated microenvironmental interactions, and immunosuppressive properties. YEATS4 knockdown effectively suppressed malignant tumor behaviors, highlighting its therapeutic potential. These findings provided novel targeted intervention strategies to address GBM heterogeneity and treatment resistance, offering promising avenues for overcoming current therapeutic limitations.
Introduction
Glioblastoma (GBM) was recognized as the most frequently diagnosed and aggressive primary brain malignancy (1), responsible for more than 48% of cases (2) and occurring at a rate of approximately 3.2 per 100,000 people each year (3). GBM predominantly affected individuals aged 50 years and older (4), with incidence figures being notably elevated in male patients (2). GBM was associated with high mortality and poor prognosis, exhibiting a 5-year survival rate of approximately 5% (5) and a median survival of only about 15 months (6, 7). Based on clinical status and disease progression, GBM was categorized into two types: newly diagnosed Glioblastoma (ndGBM) and recurrent Glioblastoma (rGBM) (8). NdGBM exhibited highly malignant and rapidly proliferative biological characteristics (1), while rGBM not only demonstrated reduced sensitivity to chemoradiotherapy but also typically displayed more advanced malignancy (9, 10).
For GBM patients, the Stupp treatment paradigm typically comprised initial surgical intervention, followed by concurrent radiotherapy and chemotherapy, and subsequent adjuvant chemotherapy (1, 11). However, it was associated with a high recurrence rate, treatment resistance, and poor prognosis (12, 13). In recent years, emerging strategies such as immunotherapy (e.g., PD-1 inhibitors) (14), targeted metabolic pathways (e.g., IDH mutations) (15), and nanoparticle drug delivery systems (16, 17) had made some progress. However, these approaches were still limited by the blood-brain barrier (BBB), immunosuppressive tumor microenvironment (TME), and resistance (12, 18, 19). Additionally, the incidence of GBM had significantly increased in elderly patients, and their prognosis was worse, further highlighting the challenges in treating GBM. Therefore, it was of great significance to explore the molecular mechanisms of GBM, resistance, and the regulation of the microenvironment to develop novel therapeutic strategies.
The breakthrough advancement of single-cell RNA sequencing (scRNA-seq) technology granted unprecedented high-definition analytical capabilities to GBM research (20–22). This technology enabled the precise identification of key cell subtypes and immune microenvironment characteristics, offering important insights into the molecular mechanisms underlying tumor recurrence and resistance. Meanwhile, scRNA-seq could dynamically monitor the evolution of tumor cell heterogeneity and microenvironment remodeling during treatment, providing molecular clues for adjusting personalized treatment strategies. The introduction of spatial transcriptomics (ST) further expanded the analytical dimensions, through its unique capacity to resolve the spatial architecture of diverse cell types in GBM tissues, especially the spatial enrichment patterns of key cell subtypes in specific tumor niches.
This study systematically analyzed the cellular heterogeneity characteristics of GBM through integrated scRNA-seq and ST technologies. We identified the C2 PCLAF+ as a critical driver of tumor progression, characterized by pronounced proliferative activity, metabolic reprogramming (such as enhanced oxidative phosphorylation and glycolysis), and distinct cell death regulation, collectively contributed to its resistance to treatment. ST analysis further showed that this subtype had distinct distribution patterns in specific tumor niches, significantly co-localizing with regions of high heterogeneity and proliferative activity. We also elucidated the interaction mechanism between the C2 PCLAF+ subtype and fibroblasts through the MDK-LRP1 ligand–receptor pair and identified key transcription factors (TFs), such as YEATS4, that regulated its malignant phenotype. Functional validation demonstrated that silencing YEATS4 significantly inhibited GBM cell proliferation, migration, and invasion. Guided by these findings, we formulated a prognostic risk scoring model that effectively distinguished patient outcomes and correlated with specific immune microenvironment characteristics, genomic variation patterns, and drug sensitivity. Our findings not only deepened the understanding of the mechanisms underlying GBM heterogeneity but also provided important theoretical foundations and potential targets for developing precision therapies targeting specific subtypes.
Materials and methods
Collection of GBM data
We retrieved data from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE182109. Following data acquisition, 16 GBM samples were selected for subsequent analysis. To complement our findings, bulk RNA-seq datasets were obtained from The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/). Meanwhile, we conducted an external validation using the dataset with GSE43378. As all data were sourced from publicly available repositories, no additional ethical approval was required for this study.
Quality control and computational analysis pipeline for scRNA-seq data
First, the initial processing of data was conducted using R (v4.3.3) and Seurat (v4.3.0) (23–25). Then, low-quality cells were eliminated through purification, filtration, and doublet exclusion in accordance with these standards: (1) nFeature detection threshold (300–5,000); (2) nCount detection threshold (500–50,000); (3) mitochondrial gene expression ratio (≤10%); (4) erythrocyte gene expression ratio (≤5%). This stringent quality control process resulted in the retention of 105,873 high-quality cells for subsequent analysis.
We first normalized the raw expression data using the “NormalizeData” function (26–28) to eliminate technical biases such as sequencing depth. Subsequently, we identified the top 2,000 highly variable genes with “FindVariableFeatures” (29, 30) to capture the most biologically significant expression differences between cells (31, 32). Following this, we applied “ScaleData” to standardize all genes, ensuring comparability across different expression levels (33). Next, we performed Principal Component Analysis (PCA) (34, 35) and corrected for batch effects using the Harmony R package (v0.1.1), improving data integration quality (36–38). Based on the top 30 principal components, we further conducted clustering analysis with “FindNeighbors” (39–41) and “FindClusters” (42–44), and finally achieved nonlinear dimensionality reduction and visualization using Uniform Manifold Approximation and Projection (UMAP) (45–47).
Cell subtype annotation and characterization
Based on the cell clustering results, we systematically annotated and characterized cell subtypes by calculating the mean expression levels of cell-type-specific marker genes obtained from the CellMarker database (http://xteam.xbio.top/CellMarker/).
Functional enrichment and AUCell analysis
We first systematically identified differentially expressed genes (DEGs) among cell subtypes using “FindAllMarkers” (48, 49), followed by comprehensive functional characterization of these DEGs with the “ClusterProfiler” package. Specifically, we simultaneously analyzed Gene Ontology (GO) terms (50–52) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (53–56) to elucidate key biological processes and signaling pathways significantly associated with the DEGs. To gain deeper functional insights, we employed Gene Set Enrichment Analysis (GSEA) to evaluate coordinated expression patterns of predefined gene sets (57–59). Furthermore, we implemented the AUCell algorithm to calculate gene activity and enrichment scores for each cell subtype, providing quantitative functional characterization at single-cell resolution.
Analysis of ST
We selected two samples from the ST dataset (GSE194329). For the ST 1 slide, we identified spatial enrichment patterns of key cell subtypes using Multimodal Intersection Analysis (MIA) and performed spatial localization with Seurat’s TransferData function. The ST 2 slide was analyzed using the Robust Cell Type Decomposition (RCTD) technology to map the key cell subtypes identified in the scRNA-seq data to the ST data, thereby further characterizing their spatial distribution features.
Differentiation trajectory analysis
We systematically investigated the developmental dynamics of glioma subtypes through multi-dimensional analyses. First, CytoTRACE was employed to evaluate the stemness levels across distinct glioma subtypes. To validate these findings, multilineage trajectory inference was performed with Slingshot (v2.6.0), where developmental paths were determined using the “getLineages” and smooth curves were generated via “getCurves”, ultimately visualizing the trajectories in UMAP. The identification of trajectory initiation and termination points was based on established biological characteristics of distinct glioma subtypes and their associated molecular signatures.
Intercellular communication network
We systematically deciphered the cell-cell interaction network within the GBM communication microenvironment using CellChat (v1.6.1) to reconstruct intercellular communication networks (60–62). The “netVisual_diffInteraction” was employed to visualize differential communication strength patterns among distinct cell types. Based on the CellChatDB database, we identified significantly activated signaling pathways and specific ligand-receptor regulatory networks.
Transcriptome regulatory landscape profiling
We employed pySCENIC (v0.10.0) to infer both the clustering of glioma cell subtypes in GBM and their specific transcriptional regulatory modules. Regulatory activity was quantified using AUCell, with significance determined through 100 permutation tests.
Cell culture
LN229 and A1207 GBM cell lines were propagated in DMEM containing 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin (P/S) at 37°C in a humidified 5% CO2 atmosphere.
Cell transfection
We employed siRNA synthesized by GenePharma (Suzhou, China) to target and silence YEATS4. For the transfection experiments, cells were incubated in 6-well plates until reaching 50% confluence, followed by transfection via the Lipofectamine 3000 RNAiMAX system (Invitrogen, USA). Two specific siRNA sequences were designed: siYEATS4-1 (GUGUAAGAAUGGAUGCTAU) and siYEATS4-2 (AAUCCUUUAAGAGUUGUUA), with non-targeting siRNA (si-NC) serving as the negative control.
Cell viability monitoring
Cells were seeded in 96-well plates at a density of 5×10³ cells per well and cultured for 24 h. Subsequently, 10 μL of CCK-8 reagent (A311-01, Vazyme) was added to each well, followed by incubation for 2 h at 37°C in the dark. Absorbance was measured at 450 nm using a microplate reader (A33978, Thermo). The experiment was conducted daily for 4 consecutive days, with the results presented as time-viability curves to illustrate dynamic changes.
Quantitative real-time polymerase chain reaction
RNA was extracted from cells using TRIzol reagent, followed by reverse transcription with the PrimeScript™ Kit. Real-time monitoring was carried out using a fluorescent dye. The primer sequences were as follows: F: TGTTCAAGAGAATGGCCGAA, R: CCACTGATGAGTGTGCCCAT.
5-ethynyl-2’-deoxyuridine proliferation assay
Transfected LN229 and A1207 cells were plated in 6-well plates at 5×10³ cells/well and cultured for 24 h. After incubation with 2× EdU working solution (2 h, 37°C), cells were washed twice with PBS, fixed with 4% paraformaldehyde (30 min), and permeabilized with 2 mg/mL glycine and 0.5% Triton X-100. For detection, cells were stained with a 1:1 mixture of 1X Apollo and Hoechst 33342 (30 min). EdU-positive cells were quantified using fluorescence microscopy.
Wound healing assay
The cells were transfected, then transferred to 6-well plates and allowed to grow to near-confluence. Linear scratches were generated in the cell monolayer via a sterile 200 μL pipette tip. After washing the wells with PBS, serum-free medium was added. Photographs of the wound areas were taken at 0 h and 48 h, and the variations in gap width were quantified with Image-J software.
Transwell assay
Serum-starved cells (24 h) were suspended in Matrigel (BD Biosciences, USA) and plated in transwell upper chambers, with serum-containing medium in the lower chambers as chemoattractant. After 48 h, cells were fixed with 4% paraformaldehyde and stained with crystal violet for quantification.
Prognostic risk scoring system for patient stratification
Based on clinical data obtained from the TCGA database, we initially identified 6 risk genes through univariate Cox regression analysis (63–66). To optimize model performance, LASSO (67–69) and multivariate Cox regression analyses were performed, ultimately selecting 5 key genes with independent prognostic value. A risk score model was constructed using the following formula: , where X represented the coefficient and Y denoted gene expression. Patients were stratified into the high and low risk groups based on the optimal cutoff value, and further evaluated through Kaplan-Meier and Receiver Operating Characteristic (ROC) curves (51, 70–72).
Assessment of immunological profiles and pharmacological vulnerability
By applying the CIBERSORT algorithm, we first systematically analyzed the immune cell infiltration profiles of both groups. Based on these findings, we further calculated immune signature scores using the ESTIMATE. Additionally, to investigate potential differences in immune evasion between the two groups, we compared the Tumor Immune Dysfunction and Exclusion (TIDE) scores. Furthermore, to evaluate the differences in therapeutic drug sensitivity, we predicted the half-maximal inhibitory concentration (IC50) values of commonly used anticancer drugs using the “pRRophetic” package (v0.5) (73).
Statistical analysis
For statistical analysis of the data, we employed R and Python. Differences between groups were assessed using the Wilcoxon rank-sum test, while Spearman’s correlation analysis was applied to examine associations between variables. All statistical tests were conducted as two-tailed tests, with a P-value <0.05 considered statistically significant, P-value < 0.01 considered highly significant, P-value < 0.001 considered very highly significant, and P-value < 0.0001 considered extremely significant.
Results
Single-cell transcriptomic profiling of GBM uncovered glioma cell heterogeneity
We performed a comprehensive bioinformatics analysis of the obtained GBM dataset, with the complete analytical workflow spanning from raw data quality control to final visualization as illustrated in Figure 1. First, we systematically analyzed 16 samples from this dataset. Through rigorous cell quality filtering and standardized procedures, we ultimately obtained 105,873 high-quality cells for subsequent investigations. As illustrated in Figure 2A, the UMAP plots visualized the distribution characteristics of 16 samples and 8 cell types (T cells and NK cells, endothelial cells (ECs), fibroblasts, oligodendrocytes, mixed, glioma cells, B-plasma cells, and myeloid cells), while revealing the distribution characteristics between 2 groups (ndGBM and rGBM) and among 3 cell cycle phases (G1, G2/M, and S). Notably, through in-depth analysis of cell cycle and tissue origin across the 8 cell types, we observed a prominent feature: glioma cells showed significantly higher distribution proportion in S phase compared to G1 and G2/M phases (Figure 2B), with these cells predominantly originating from the ndGBM (Figure 2C). Further characterization revealed that glioma cells exhibited relatively high levels in multiple biological indicators, including pMT, Cell-Stemness-AUC, nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score (Figures 2D, E). Figure 2F displayed the mean expression levels of the top 5 DEGs in each cell type. Of particular note, the signature genes of glioma cells included CLU, CHI3L1, PTPRZ1, MT3, and PTN, which may be closely associated with their malignant phenotype. Functional enrichment analysis demonstrated that biological processes related to glioma cells were primarily concentrated in energy metabolism pathways, including aerobic electron transport chain, mitochondrial ATP synthesis coupled electron transport, ATP synthesis coupled electron transport, and proton motive force-driven ATP synthesis (Figure 2G). Finally, we quantified metabolic pathway activities using the AUCell scores. The results identified oxidative phosphorylation, glycolysis/gluconeogenesis, pyruvate metabolism, cysteine and methionine metabolism, and propanoate metabolism as the top 5 most active metabolic pathways among the 8 cell types (Figure 2H). These findings demonstrated that glioma cells possessed high heterogeneity, proliferative activity, and metabolic characteristics, which was consistent with previously reported mechanisms whereby glioma cells drive tumor heterogeneity through metabolic abnormalities and reactive oxygen species accumulation (74). Based on these findings, we further systematically characterized the subtypes of glioma cells to elucidate the intrinsic relationship between their molecular heterogeneity and clinical phenotypes.

Figure 1. The workflow diagram for GBM scRNA-seq. The analysis steps for GBM scRNA-seq encompassed data preprocessing, characterization of key cell subtypes, and exploration of their biological significance.

Figure 2. Comprehensive Single-Cell Characterization of GBM. (A) UMAP plots illustrated the distribution of 16 different samples (upper left), 8 cell types (upper right), various groups (lower left), and different cell cycle phases (lower right) across the entire cell population. (B) The stacked bar chart displayed the proportional distribution of each cell type across different cell cycle phases. (C) Distribution of different cell types across the two groups in proportion. (D, E) Bar plots and UMAP plots respectively illustrated the expression and distribution of pMT, Cell-Stemness-AUC, nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score across different cell types. (F) The bubble plot demonstrated the mean expression levels of the top 5 DEGs in each cell type. (G) GSEA described the enrichment scores of biological processes associated with glioma cells. ****P < 0.0001. (H) Heatmap assessed the AUCellscores of the top metabolism-related pathways across 8 cell types.
Identification and multi-omics characterization of the C2 PCLAF+ subtype
Given the critical role of glioma cellular heterogeneity in GBM pathogenesis and treatment resistance, we systematically characterized their molecular features. Based on marker gene expression patterns, we identified 9 distinct glioma cell subtypes (Figure 3A), including C0 CHI3L1+ glioma cells, C1 MMP16+ glioma cells, C2 PCLAF+ glioma cells, C3 BIRC5+ glioma cells, C4 NDRG1+ glioma cells, C5 TYROBP+ glioma cells, C6 STMN2+ glioma cells, C7 CXCR4+ glioma cells, and C8 MX1+ glioma cells. UMAP visualization revealed significant transcriptomic differences and unique molecular characteristics among these subtypes. To further characterize the biological properties of these subtypes, we systematically analyzed their transcriptional activity and cell cycle features. The results demonstrated that the C2 PCLAF+ subtype exhibited significantly higher expression levels of nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score compared to other subtypes, particularly in S.Score (Figure 3B). These findings suggested that this subtype likely possessed stronger proliferative activity. Through in-depth analysis of cell cycle distribution and tissue origin across different subtypes (Figures 3C, D), we observed that the C2 PCLAF+ subtype showed remarkably higher proportion in S phase and was predominantly derived from rGBM, characteristics that were highly consistent with its hyperproliferative phenotype. To elucidate the molecular characteristics of each subtype, we presented the top 5 marker genes for each subtype in Figure 3E. The C2 PCLAF+ subtype was characterized by PCLAF, TYMS, PCNA, DUT, and HIST1H4C as its top marker genes.

Figure 3. Molecular characterization and spatial distribution of the C2 PCLAF+ subtype. (A) The central part of the circular plot outlined the clustering of 9 glioma cell subtypes using contour lines, while the outer, middle, and inner axes displayed the logarithmic scale of clusters for each glioma cell subtype, as well as the distribution of groups and phases. The peripheral UMAP plots illustrated the distribution of nCount-RNA (upper left), nFeature-RNA (upper right), G2/M.Score (lower left), and S.Score (lower right) within glioma cells. (B) Bar plots visualized the expression levels of nCount-RNA, nFeature-RNA, G2/M.Score, and S.Score in each glioma cell subtype. (C) The stacked bar chart displayed the proportional distribution of each glioma cell subtype across different cell cycle phases. (D) Box plots illustrated the frequency of 9 glioma cell subtypes within each group. (E) The bubble plot displayed the top 5 marker genes in each glioma cell subtype. The pie charts showed the proportions of different phases, while the violin plots illustrated the expression levels of G2/M.Score, S.Score, and nCount-RNA. (F) Heatmap showed the enrichment analysis results of GB-BP and KEGG for 9 glioma cell subtypes. (G) The bubble plot illustrated the metabolic pathways in each glioma cell subtype. (H) The heatmaps assessed the AUCell scores of the top 20 metabolism-related pathways across the C2 PCLAF+ subtype and S phase. (I, J) Bar graphs and UMAP plots, based on the AUCell algorithm, respectively evaluated the activity and distribution of oxidative phosphorylation, glycolysis/gluconeogenesis, and pyruvate metabolism across different glioma cell subtypes. (K) Heatmap displayed the mean expression levels of programmed cell death pathways in each glioma cell subtype. (L) The AUCell algorithm was used to calculate the activities of oxeiptosis, cuproptosis, and parthanatos in each glioma cell subtype. (M) Heatmap displayed the enrichment degree of 9 spatial clusters associated with various cell types on ST 1 slide, as assessed by MIA analysis. (N) The visualization results of the spatial spots corresponding to the 9 spatial clusters were presented. (O) ST feature map illustrated the specific spatial expression pattern of the C2 PCLAF+ subtype.
To further investigate the biological characteristics of different glioma cell subtypes, we conducted functional enrichment analysis (Figure 3F). The results indicated that the C0 CHI3L1+ subtype was mainly enriched in biological processes such as response to oxygen levels, response to decreased oxygen levels, response to hypoxia, cell-substrate adhesion, and cell growth, and was associated with pathways including proteoglycans in cancer, ECM-receptor interaction, cytoskeleton in muscle cells, focal adhesion, and HIF-1 signaling pathway. The C1 MMP16+ subtype was primarily enriched in biological processes such as cytoplasmic translation, oxidative phosphorylation, aerobic electron transport chain, proton transmembrane transport, and aerobic respiration, and was related to amyotrophic lateral sclerosis and pathways of neurodegeneration-multiple diseases. The C2 PCLAF+ subtype was predominantly enriched in biological processes such as DNA replication, DNA-templated DNA replication, DNA recombination, double-strand break repair, and chromosome segregation, and was related to pathways such as DNA replication and cell cycle. The C3 BIRC5+ subtype was mainly enriched in biological processes such as chromosome segregation, sister chromatid segregation, mitotic sister chromatid segregation, mitotic nuclear division, and nuclear chromosome segregation, and was related to pathways such as amyotrophic lateral sclerosis and Alzheimer disease. The C4 NDRG1+ subtype was mainly enriched in biological processes such as cytoplasmic translation, pyruvate metabolic process, glycolytic process, NADH regeneration, and canonical glycolysis, and was related to pathways such as ribosome, biosynthesis of amino acids, glycolysis/gluconeogenesis, and carbon metabolism. The C5 TYROBP+ subtype was mainly enriched in biological processes such as positive regulation of leukocyte activation, positive regulation of cell activation, leukocyte cell-cell adhesion, positive regulation of lymphocyte activation, and chemotaxis. The C6 STMN2+ subtype was mainly enriched in biological processes such as axonogenesis, forebrain development, regulation of neuron projection development, telencephalon development, and axon guidance, and was related to pathways such as motor proteins, ATP-dependent chromatin remodeling, salmonella infection, and synaptic vesicle cycle. The C7 CXCR4+ subtype was mainly enriched in biological processes such as leukocyte cell-cell adhesion, T cell differentiation, regulation of T cell activation, lymphocyte differentiation, and antigen receptor-mediated signaling pathway, and was related to pathways such as ribosome, antigen processing and presentation, Th17 cell differentiation, and Th1 and Th2 cell differentiation. The C8 MX1+ subtype was mainly enriched in biological processes such as response to virus, defense response to virus, defense response to symbiont, regulation of viral process, and negative regulation of viral genome replication, and was related to pathways such as antigen processing and presentation, and proteasome. Regarding metabolic features, we noticed that the metabolic pathways of 3 glioma cell subtypes are quite distinct (Figure 3G). The C2 PCLAF+ subtype was primarily active in pathways such as pyruvate metabolism, propanoate metabolism, oxidative phosphorylation, glutathione metabolism, and cysteine and methionine metabolism. The C3 BIRC5+ subtype was predominantly enriched in pyruvate metabolism and propanoate metabolism, whereas the C4 NDRG1+ subtype showed significant enrichment in pyruvate metabolism, propanoate metabolism and glycolysis/gluconeogenesis. Quantitative comparison of metabolic pathway activities using AUCell scores (Figure 3H) revealed that the top 20 metabolic pathways of the C2 PCLAF+ subtype matched those of S phase. Further visualization analysis (Figures 3I, J) revealed that compared with other subtypes, the C2 PCLAF+ subtype showed significantly enhanced activity in 3 critical metabolic pathways: oxidative phosphorylation, glycolysis/gluconeogenesis, and pyruvate metabolism. Additionally, analysis of programmed cell death pathways (Figures 3K, L) demonstrated that this subtype showed relatively higher AUC scores in oxeiptosis, cuproptosis, and parthanatos pathways, but lower AUC scores in alkaliptosis and netotic cell death pathways. This distinctive expression pattern implied the existence of unique cell death regulation mechanisms in this subtype.
STs revealed the spatial expression patterns of proliferative C2 PCLAF+ subtype
To further investigate the spatial heterogeneity of GBM, we selected two representative GBM sections for ST analysis. In ST 1 slide, MIA and TransferData analyses revealed that the C2 PCLAF+ subtype was predominantly enriched in spatial cluster C1 (Figure 3M). Comparative analysis of spatial distribution patterns demonstrated a high degree of similarity between the spatial localization pattern of cluster C1 (Figure 3N) and the distribution characteristics of the C2 PCLAF+ subtype (Figure 3O). To validate this finding, we employed RCTD analysis to map the spatial regions with high expression of the C2 PCLAF+ subtype (Supplementary Figure 1A), which showed significant spatial co-localization with regions exhibiting high PCLAF gene expression (Supplementary Figure 1B). In ST 2 slide, the spatial distribution patterns of nCount-Spatial, nFeature-Spatial, G2/M.Score, and S.Score were highly consistent with the characteristics of the C2 PCLAF+ subtype (Supplementary Figure 1C), further confirming its unique spatial organization within the GBM architecture.
C2 PCLAF+ subtype-driven differentiation promoted GBM progression
To probe the developmental origin of GBM and elucidate its heterogeneity, we systematically investigated the differentiation characteristics of glioma cell subtypes. We used CytoTRACE to predict the stemness ranking of 9 glioma cell subtypes and found that the C2 PCLAF+ subtype exhibited high stemness (Figures 4A, B). This result suggested that the C2 PCLAF+ subtype might possess stronger stem cell-like properties, laying the foundation for further investigation into its differentiation potential. We evaluated the mean expression of top stemness genes in each cell cluster and visualized the distribution of top stemness genes (EZH2, LGR5, NOTCH1, and ABCG2) in the C2 PCLAF+ subtype (Figures 4C, D). To comprehensively delineate the differentiation trajectories of glioma cells, we employed Slingshot to assess the differentiation lineages of the 9 glioma cell subtypes (Figure 4E). Notably, the analysis demonstrated that all 4 lineages originated from the C2 PCLAF+ subtype, further supporting its pivotal role as a differentiation initiator. Subsequently, we visualized the progression of lineage 1, lineage 2, lineage 3, and lineage 4 along the inferred pseudotemporal trajectories (Figure 4F), which provided more intuitive evidence for the central position of the C2 PCLAF+ subtype in the differentiation process and revealed distinct differentiation pathways among the lineages. Finally, to explore the biological characteristics of the different differentiation lineages, we further examined the major biological processes associated with the 4 lineages during glioma cell differentiation (Figure 4G). Functional enrichment analysis revealed that cluster C1 was primarily enriched in healing wound and chemotaxis processes, cluster C2 was mainly associated with mitotic segregation and spindle-related functions, cluster C3 showed strong correlations with leukocyte activation and immune responses, while cluster C4 was significantly linked to immune-mediated and adaptive mechanisms.

Figure 4. Differentiation dynamics and lineage characteristics of glioma cell subtypes. (A) The two-dimensional graphs utilized CytoTRACE to assess the predicted order (left) and distribution characteristics (right) of the 9 glioma cell subtypes. (B) Box plot demonstrated the ordering of the 9 glioma cell subtypes as predicted by the CytoTRACE analysis. (C) Heatmap displayed the mean expression of top stemness genes in the 9 cell clusters. (D) UMAP plots visualized the distribution of top stemness genes in the C2 PCLAF+ subtype. (E) The Slingshot analysis depicted 4 differentiation lineages of the 9 glioma cell subtypes, with the C2 PCLAF+ subtype appearing in the early stages of differentiation. (F) UMAP plots represented the progression of 4 lineages along the inferred pseudotime trajectories. (G) The GO-BP enrichment analysis results revealed key biological processes associated with the 4 lineages during glioma cell differentiation.
Interplay between the C2 PCLAF+ subtype and fibroblasts through the MDK-LRP1 ligand–receptor pair
Building upon the pivotal role of the C2 PCLAF+ subtype in GBM differentiation, we deeply analyzed its molecular mechanisms in regulating microenvironmental communication. First, we quantitatively assessed the number and intensity of interactions between the C2 PCLAF+ subtype and all other cell types, revealing enhanced communication between this subtype and fibroblasts (Figure 5A). To comprehensively characterize the intercellular communication network, we further analyzed the outgoing and incoming communication patterns among 9 glioma cell subtypes and 7 other cell types (Figure 5B). Notably, the analysis demonstrated that both outgoing and incoming signals of the C2 PCLAF+ subtype were predominantly regulated by Pattern 1, with significant contributions from key molecules including MPZ, CDH, NRXN, NCAM, PTN, MK, and CD99. These molecules likely collectively formed the molecular basis of C2 PCLAF+ subtype communication. Of particular note, our analysis of communication probabilities within the MK signaling network (Figure 5C) showed increased communication likelihood between the C2 PCLAF+ subtype and fibroblasts. To elucidate the mechanism underlying this specific communication, we performed in-depth analysis of functional roles for various cell types within the MK signaling pathway network (Figure 5D). The results indicated that the C2 PCLAF+ subtype exhibited multifaceted communication functions, primarily serving as sender, receiver, mediator, and influencer, while fibroblasts mainly participated as receivers and influencers in the communication network. Additionally, through visual analysis, we found that the key ligand MDK and its receptor LRP1 in the MK signaling pathway were highly expressed in both the C2 PCLAF+ subtype and fibroblasts (Figures 5E, F). Based on this, we further focused on the MDK-LRP1 signaling axis and confirmed significant cellular interactions between the C2 PCLAF+ subtype and fibroblasts (Figures 5G, H).

Figure 5. Single-cell communication atlas dominated by the C2 PCLAF+ subtype. (A) Circle diagrams illustrated the interactions of the C2 PCLAF+ subtype, presenting it as a source (left) and a target (right) in relation to other cell types, highlighting both the number (upper) and the intensity (lower) of these interactions. (B) Heatmaps displayed the outgoing (left) and incoming (right) communication patterns of 9 glioma cell subtypes and 7 other cell types under 3 cell communication patterns, as well as the key proteins involved in these 3 communication patterns. (C) Heatmap depicted the communication probabilities of 9 glioma cell subtypes and 7 other cell types under the MK signaling network. (D) Heatmap illustrated the importance of 9 glioma cell subtypes and 7 other cell types within the MK signaling pathway network, including sender, receiver, mediator, and influencer. (E, F) The bubble plot and violin plots visualized the expression patterns of essential ligands and receptors within the MK signaling pathway, encompassing 9 glioma cell subtypes and 7 additional cell types. (G, H) The hierarchical graph and circle diagram depicted the interactions among 9 glioma cell subtypes and 7 other cell types within the MDK-LRP1 ligand–receptor pair.
TF activity mapping identified proliferation-driving characteristics in C2 PCLAF+ subtype
To elucidate the molecular characteristics of different glioma cell subtypes, we performed a re-clustering analysis of glioma cells based on regulon activity and integrated the cell cycle phases. The results revealed that the C2 PCLAF+ subtype represented a large proportion in the S phase (Figure 6A), suggesting its high DNA replication activity and potential association with tumor proliferation. To further explore the TF network regulating the C2 PCLAF+ subtype, we divided functionally and expression-similar TFs into two regulatory modules, M1 and M2, based on AUCell similarity scores (Figure 6B). UMAP visualization demonstrated distinct distribution patterns of TFs in M1 and M2 across glioma cells (Figure 6C). Notably, the C2 PCLAF+ subtype exhibited higher expression in M2 (Figure 6D), and S-phase cells within M2 showed elevated TF expression (Figure 6E). Further analysis revealed that both the C2 PCLAF+ subtype and S-phase cells displayed relatively higher regulon activity scores in the M2 (Figures 6F, G), further supporting the link between this subtype and proliferative activity. To decipher the core regulatory mechanisms of the C2 PCLAF+ subtype, we identified the top 5 TFs with the highest specificity scores in each subtype (Figure 6H). In the C2 PCLAF+ subtype, the top 5 TFs were YEATS4, E2F1, TFDP1, MYBL1, and CREM, indicating their potential cooperative role in driving the proliferative phenotype of this subtype. Ultimately, we delineated the phase-specific distribution characteristics of TFs during G1, G2/M and S phases using UMAP (Figure 6I) and specifically highlighted the distribution patterns of the top 5 TFs in the C2 PCLAF+ subtype within the M2 module (Figure 6J). We observed significantly higher expression of YEATS4 in the C2 PCLAF+ subtype, prompting further in vitro experiments to validate the functional role of YEATS4 in GBM progression.

Figure 6. Decoding the transcriptional regulatory network of the C2 PCLAF+ subtype. (A) UMAP plots displayed the distribution of clustering for 9 glioma cell subtypes based on regulon activity level, and the pie charts illustrated the proportions of their phases. (B) Heatmap revealed two regulatory modules (M1 and M2) identified based on TF recognition in glioma cell subtypes. (C) UMAP plots showed the distribution of TFs expression in M1 and M2. (D, E) Bar plots visualized the AUC values of different glioma cell subtypes and phases within M1 and M2. **P < 0.01, ****P < 0.0001. (F, G) Scatter plots exhibited the ranking of regulon activity scores for various glioma cell subtypes and phases within M1 and M2. (H) Scatter plots presented the ranking of specificity scores for the top 5 TFs in each subtype of glioma cells. (I) UMAP plots visualized the expression distribution of TFs in the G1, G2/M, and S phases. (J) UMAP plots visualized the expression distribution of the top 5 TFs in the C2 PCLAF+ subtype within M2.
In vitro experimental validation of the functional characterization of YEATS4 in GBM progression
We conducted in vitro experiments using LN229 and A1207 GBM cell lines to investigate the functional consequences of YEATS4 knockdown. These two cell lines were divided into 3 treatment groups: si-NC (negative control), siYEATS4-1, and siYEATS4-2. As shown in Figure 7A, the mRNA and protein expression levels of YEATS4 were significantly reduced following knockdown. We monitored the average optical density values of the 2 GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4–2 over time (0 to 4 days) and found that YEATS4 knockdown significantly inhibited cell viability (Figure 7B). Subsequently, colony formation assays revealed a marked reduction in the number of colonies after YEATS4 knockdown (Figures 7C, D). Aligned with these discoveries, EDU staining confirmed that cell proliferation was impaired following YEATS4 knockdown (Figures 7D, E). Moreover, wound healing assays demonstrated significantly delayed wound closure at 48 h post-knockdown (Figures 7F, G). Most strikingly, transwell assays showed that the migratory and invasive capacities of both GBM cell lines were simultaneously impaired following YEATS4 knockdown (Figures 7G, H).

Figure 7. In vitro experimental verification. (A) Bar plots presented the relative expression levels of YEATS4 mRNA and protein in LN229 and A1207 GBM cell lines. si-NC, siYEATS4-1, and siYEATS4–2 represented different siRNA treatment groups, which were used to compare the effects of YEATS4 knockdown on expression levels. (B) The line graphs illustrated the cell viability over time (0 to 4 days) in two GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4-2. (C) The colony formation assay was conducted to compare the effects of si-NC, siYEATS4-1, and siYEATS4–2 on colony formation in two GBM cell lines. (D) Bar plots illustrated the colony numbers and cell proliferation in two GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4-2. (E) The EDU staining assay was performed to assess cell proliferation in two GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4-2. (F) The cell wound healing assay evaluated the wound healing ability of two GBM cell lines at 0 h and 48 h following treatment with si-NC, siYEATS4-1, and siYEATS4-2. (G) Bar plots illustrated the wound healing percentage, cell migration, and cell invasion in two GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4-2. (H) Transwell assay evaluated the cell migration (upper) and invasion capabilities (lower) of two GBM cell lines treated with si-NC, siYEATS4-1, and siYEATS4-2. *P < 0.05, **P < 0.01, ***P < 0.001.
The construction and analytical validation of a prognostic risk model based on C2 PCLAF+ glioma cells
The prognosis assessment of GBM urgently required reliable molecular markers. To address this, we established a risk scoring model founded on C2 PCLAF+ glioma cells. In the process of constructing the model, we first preliminarily screened 6 risk genes closely tied to patient survival (Figure 8A). To further optimize the specificity of the model, LASSO regression analysis was employed to select candidate genes (Figure 8B), followed by multivariate Cox regression analysis, which ultimately identified 5 genes with independent prognostic value (Figure 8C). Notably, coefficient value analysis revealed that the positive coefficients of IRF7, FOSL2, IKZF3, and CEBPD indicated their high expression portended poorer patient outcomes (Figure 8D). To validate the clinical significance of these genes, patients were dichotomized by the PCLAF+ glioma cells risk score (PGRS) cutoff (Figure 8E). Comparative analysis demonstrated that significant upregulation of IRF7, FOSL2, IKZF3, and CEBPD, along with downregulation of SOX6, in the high PGRS group compared to controls (Figure 8F). This result was highly consistent with survival analysis: Kaplan-Meier curves showed that the OS rate was lower in the high PGRS group (Figure 8G), further supporting the prognostic value of PGRS stratification. Meanwhile, external datasets were used for validation, and the results showed that the two sets of data had good consistency (Figure 8H). The predictive performance of the model was validated using time-dependent ROC curves. As shown in Figure 8I, the model demonstrated high sensitivity and specificity in predicting 1-year, 3-year, and 5-year survival, with AUC values indicating reliable discriminative ability. More importantly, the successful validation of the external validation dataset has significantly enhanced the clinical applicability and generalization ability of the model (75). Follow-up studies assessed connections between prognostic genes and clinical outcomes: IRF7, FOSL2, IKZF3, and CEBPD were not only significantly negatively correlated with OS but also positively correlated with risk scores, whereas SOX6 exhibited the opposite pattern—its high expression was associated with prolonged OS and reduced risk scores (Figure 8J). This association was further confirmed by gene expression stratification analysis: as shown in Figure 8K, patients with high expression of IRF7, FOSL2, and CEBPD had significantly reduced OS rates, while those with high SOX6 expression showed survival advantages. Finally, we observed significant differences in the expression patterns of these 5 prognostic genes between the two groups (Figure 8L): IRF7, FOSL2, IKZF3, and CEBPD were generally highly expressed in the high PGRS group, whereas SOX6 was specifically enriched in the low PGRS group.

Figure 8. Development of a risk scoring model. (A) The forest plot demonstrated the results of univariate Cox regression analysis for key genes. Hazard ratios (HR) and their confidence intervals (lines) were used to evaluate prognostic associations: an HR < 1 indicated a protective factor, while an HR > 1 indicated a risk factor. (B) LASSO regression results. Optimal lambda selection via cross-validation (upper); Coefficient trajectories (upward for significant genes) at optimal lambda (lower). (C) The forest plot displayed the results of multivariate Cox regression analysis for key genes. (D) Bar plot compared the distribution of coefficients for the selected genes. (E) The survival curve combined with scatter plot revealed the survival rate trends over time in the high PGRS group (orange) and low PGRS group (green), with event endpoints including death and alive. (F) Heatmap demonstrated the comparative results of prognostic gene expression levels between high PGRS group and low PGRS group. (G, H) The Kaplan-Meier survival curves compared the OS rates between the high and low PGRS groups in the research (upper) and external cohorts (lower). The table at the bottom displayed the number at risk at each annual time point. (I) The ROC curves depicted the area AUC for 1-year, 3-year, and 5-year survival predictions in the research (left) and external cohorts (right). (J) The matrix heatmap combined with the scatter plot depicted the correlations between prognostic genes and both OS and risk scores. (K) The Kaplan-Meier survival curves compared the OS rates among patients stratified by expression levels of the 5 prognostic genes. (L) Box plots showed differential expression of the 5 prognostic genes between high and low PGRS groups. ***P < 0.001.
The PGRS risk scoring system revealed GBM immune microenvironment characteristics and predicted immunotherapy sensitivity
Building upon our previously established PGRS scoring system, we proceeded to explore its immunomodulatory implications. We noted significant heterogeneity in the infiltration proportions of various immune cell types in Figure 9A, which prompted us to explore their relationship with the PGRS score. In-depth analysis revealed that B cells naive, Neutrophils, and Mast cells activated were significantly enriched in the high PGRS group, while B cells memory and T cells follicular helper were predominant in the low PGRS group (Figure 9B). Notably, this distribution pattern was consistent with the correlation analysis of risk scores: T cells CD4 memory activated, Monocytes, and Neutrophils showed significant positive correlations, whereas T cells follicular helper, B cells memory, and Eosinophils exhibited clear negative trends (Figure 9C). Collective evidence collectively suggested that the PGRS score was intimately related to the infiltration characteristics of specific immune cell types. Further examination of the immune microenvironment features revealed a robust interaction network among different immune cell types, with generally high pairwise correlation strengths (Figure 9D). More importantly, the 5 key prognostic genes demonstrated significant correlation patterns with multiple immune cell types (Figure 9E), among which the strong positive correlations between IRF7 and Macrophages M1, as well as IKZF3 and B cells naive, were particularly prominent. These molecular-immune associations provided new insights into the regulatory mechanisms of prognostic genes. In the genomic feature analysis, we systematically evaluated the mutation profiles of 461 samples and found somatic mutations in 405 (87.85%) cases, with missense mutations being the predominant type. Notably, PTEN, TP53, and EGFR exhibited relatively high mutation frequencies (Figure 9F). The CNV analysis indicated the CNV gains in the 5 prognostic genes, particularly prominent for IRF7 and IKZF3 (Figure 9G). To comprehensively assess the potential for immunotherapy response, we compared TIDE and signature scores between high and low PGRS groups. The results showed that the high PGRS group had higher TIDE, stromal, immune, and ESTIMATE scores (Figures 9H, I). This difference was further corroborated in the immune checkpoint analysis: most immune checkpoint-related genes not only showed positive correlations with risk scores but also exhibited significant co-upregulation with IRF7, FOSL2, IKZF3, and CEBPD, while only SOX6 displayed a negative correlation pattern (Figure 9J). Gene expression profiling confirmed that most immune checkpoint-related genes, such as CD276, CD44, and C10orf54, were significantly upregulated in the high PGRS group (Figure 9K). These results jointly pointed to increased immune checkpoint inhibitor susceptibility in the high PGRS group.

Figure 9. The pivotal role of PGRS scoring in the regulation of the GBM immune network. (A) Box plot demonstrated the immune cell infiltration landscape. (B) Box plot revealed the proportions of immune cell infiltration between high and low PGRS groups. (C) The lollipop plot depicted the association analysis between different immune cell types and risk scores. (D) The matrix heatmap displayed pairwise correlation strengths between different immune cell types. (E) Heatmap showcased correlation strengths between 5 prognostic genes, OS, risk scores, and various immune cell types. (F) The waterfall plot integrated with bar charts displayed the mutation frequencies and variant types of the top 20 mutated genes. (G) The stacked bar plot presented the frequency analysis of CNV for 5 prognostic genes. (H) Box plot combined with violin plot compared TIDE scores between high and low PGRS groups. (I) Box plot illustrated the stromal, immune, and ESTIMATE scores between high and low PGRS groups. (J) The bubble plot revealed the correlation analysis between 5 prognostic genes, OS, risk scores, and immune checkpoint-related genes. (K) Box plot highlighted differential expression of immune checkpoint-related genes between high and low PGRS groups. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.
Functional enrichment characteristics of DEGs and drug sensitivity analysis
To delineate the molecular characteristic distinctions across PGRS stratification groups and their potential clinical significance, we comprehensively analyzed the DEGs, functional pathway characteristics, and drug sensitivity differences between the two groups. In Figure 10A, we displayed the significantly upregulated or downregulated genes between the high and low PGRS groups and further visualized their expression patterns (Figure 10B). Subsequently, we performed systematic functional annotation analysis on these DEGs. Biological processes analysis revealed that the DEGs were highly abundant in humoral immune response, antimicrobial humoral response, and antimicrobial humoral immune response mediated by antimicrobial peptide (Figure 10C). Regarding cellular components, the DEGs were primarily associated with structures such as specific granule lumen, specific granule and tertiary granule (Figure 10D). Molecular functions analysis indicated that these genes were significantly enriched in chemokine activity, cytokine activity and chemokine receptor binding (Figure 10E). KEGG enrichment analysis further supported these findings, showing that the DEGs were significantly enriched in pathways including viral protein interaction with cytokine and cytokine receptor, and IL-17 signaling pathway (Figure 10F). To further elucidate the biological significance of the DEGs, we conducted GSEA enrichment analysis, which demonstrated that these genes were significantly involved in key biological processes such as adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, adaptive immune response, neutrophil chemotaxis, and granulocyte chemotaxis (Figure 10G). Finally, the drug sensitivity results disclosed that the high PGRS group exhibited lower IC50 values for Cyclopamine, Dasatinib, Sunitinib, Erlotinib, Lapatinib, PF.02341066, and Tipifarnib, while showing higher IC50 values for ATRA, Axitinib, and AZD.2281, revealing a significant differential drug sensitivity between the high and low PGRS groups (Figure 10H). The study outcomes provided a potential basis for the selection of clinical treatment strategies.

Figure 10. Functional enrichment and drug response in PGRS-stratified groups. (A) The volcano plot highlighted the significantly upregulated or downregulated genes between high and low PGRS groups. (B) Heatmap illustrated the DEGs expression levels between high and low PGRS group. (C–E) GO enrichment analysis of DEGs across biological processes, cellular components, and molecular functions. (F) KEGG enrichment analysis of DEGs revealed significantly enriched biological pathways. (G) GSEA enrichment analysis of DEGs showed significantly enriched biological processes, with positive and negative normalized enrichment scores (NES) indicating directional activation or suppression. (H) Box plots combined with violin plots evaluated the IC50 values of different drugs between high and low PGRS groups. ***P < 0.001.
Discussion
GBM was widely recognized as one of the most challenging malignant tumors due to its high heterogeneity, aggressive invasiveness, high recurrence rate, and treatment difficulty (16). The treatment difficulties stemmed from its complex pathogenesis and multiple biological characteristics, including high heterogeneity, the presence of stem-like cells, along with the generation of an immunosuppressive milieu. These features intertwined and collectively created significant obstacles to treatment. The tumor’s high heterogeneity not only led to poor efficacy of targeted therapies but also caused GBM to develop intrinsic resistance to most existing anticancer drugs (12). At the same time, the stem-like cell population within the TME further exacerbated multidrug resistance (16, 76), while the BBB severely limited the effective delivery of therapeutic drugs (6, 12). Even more troublesome, tumor cells escaped immune surveillance by shaping an immunosuppressive microenvironment (77, 78), making it highly prone to recurrence even after aggressive treatment. These factors worked together, causing the current comprehensive treatment regimen, primarily based on surgical resection combined with radiotherapy and chemotherapy, to only extend survival to a certain extent, but the patient’s prognosis remained very poor. In the face of this dire situation, developing new therapeutic strategies capable of overcoming tumor heterogeneity, breaching the BBB, and regulating the immune microenvironment became an urgent direction of research (79, 80).
Since GBM is fundamentally a type of glioma (16), we specifically focused on the biological characteristics of glioma cells. The findings revealed that glioma cells exhibited significantly higher heterogeneity, proliferative activity, and metabolic properties compared to other cell types, prompting us to conduct a more in-depth analysis of their cellular subtypes. Utilizing scRNA-seq, we identified 9 distinct glioma cell subtypes. Among these, the C2 PCLAF+ subtype displayed the most unique proliferative and metabolic features. From a metabolic perspective, the C2 PCLAF+ subtype showed significantly markedly enhanced activity in key metabolic pathways such as oxidative phosphorylation, glycolysis/gluconeogenesis, and pyruvate metabolism. This not only revealed its characteristic of maintaining energy supply through the synergy of aerobic glycolysis and mitochondrial metabolism but also directly suggested that this subtype might drive treatment resistance through metabolic reprogramming (81). Specifically, the elevated glycolytic activity helped GBM cells survive in hypoxic microenvironments by rapidly generating ATP and NADH, and provided substrates for DNA repair after chemotherapy. Meanwhile, the enhanced oxidative phosphorylation inhibited apoptotic signals by maintaining mitochondrial membrane potential, thereby resisting cell death induced by radiotherapy or targeted therapy. Existing studies confirmed that the enhancement of glycolysis in GBM was closely related to temozolomide resistance (82). In ST 1 slide, MIA and TransferData analysis demonstrated the unique spatial distribution pattern of C2 PCLAF+ subtype in GBM tissue, suggesting that this subtype might establish specific tumor niches to promote recurrence and drug resistance. In ST 2 slide, RCTD analysis showed that high-expression regions of C2 PCLAF+ subtype exhibited significant co-localization with spatial expression patterns of PCLAF gene, and its distribution characteristics closely matched areas of high heterogeneity and proliferative activity. These observations were consistent with single-cell sequencing data, further supporting the crucial role of C2 PCLAF+ subtype in rapid tumor proliferation and invasive growth. Moreover, the spatial distribution of this subtype within tumor tissue was not random but rather preferentially enriched in specific microenvironmental regions, which might be associated with its metabolic adaptability and treatment resistance capabilities.
The C2 PCLAF+ subtype played a central regulatory role in the differentiation of GBM, driving tumor heterogeneity evolution through its higher stemness levels and early differentiation potential. This subtype exhibited significant stem cell-like characteristics, marked by the high expression of key stemness genes such as EZH2, LGR5, NOTCH1, and ABCG2, which together formed a core regulatory network maintaining tumor stem cell properties. Slingshot analysis further confirmed that all 4 differentiation lineages originated from the C2 PCLAF+ subtype, suggesting it might be a key origin subtype driving the multi-directional differentiation of GBM. Furthermore, the functional characteristics of the 4 lineages revealed major biological processes during GBM differentiation. These findings systematically delineated the dynamic evolutionary landscape of GBM, providing new theoretical foundations for understanding the mechanisms underlying tumor heterogeneity formation.
By analyzing intercellular communication patterns, we found that in the MK signaling network, the communication likelihood between the C2 PCLAF+ subtype and fibroblasts was significantly increased. The MK signaling pathway exhibited significant functional heterogeneity across different tumor types: In cervical cancer, this pathway synergistically promoted HPV-related progression through ECM remodeling mediated by cancer-associated fibroblasts and immunosuppression (83); in hepatocellular carcinoma, it promoted angiogenesis and metastasis, and coupled lipid metabolism to drive immunosuppression (84); and in lung adenocarcinoma, it collaboratively drove tumor progression with EGFR signaling and induced T cell rejection (85). Previous studies have demonstrated that the MK signaling pathway played a pivotal role in the malignant progression of GBM through multiple mechanisms, including mediating TME remodeling, inducing the formation of an immunosuppressive microenvironment, and driving malignant tumor behaviors (86). Of particular importance, we identified that both the key ligand MDK and its receptor LRP1 in the MK signaling pathway exhibited high expression characteristics in both C2 PCLAF+ subtypes and fibroblasts. This pivotal discovery prompted us to focus our investigation on the MDK-LRP1 signaling axis. Ultimately, we demonstrated that the C2 PCLAF+ subtype established specific interaction mechanisms with fibroblasts in the TME through the MDK-LRP1 ligand–receptor pair.
To elucidate the transcriptional regulatory mechanisms of the C2 PCLAF+ subtype in GBM, we systematically analyzed its specific transcriptional regulatory network. When further dissecting the core regulatory mechanisms of the C2 PCLAF+ subtype, we identified the top 5 TFs with the highest specificity scores in this subtype. Notably, the significant upregulation of YEATS4 in the C2 PCLAF+ subtype suggested its potential critical role in GBM progression. Current studies demonstrated that YEATS4 (GAS41), as a crucial TF (87, 88), played a pivotal role in the occurrence and development of various tumors. In colorectal cancer, YEATS4 drove tumor cell proliferation by promoting cell cycle transition and inhibiting apoptosis (89); in pancreatic cancer, YEATS4 enhanced the stemness of tumor cells and mediated resistance to gemcitabine by binding to H2A.Z.2 and activating the Notch1 signaling pathway (90); in non-small cell lung cancer, YEATS4 amplification antagonized cellular senescence by inhibiting the p53-p21 pathway, promoted proliferation, and led to cisplatin resistance (91). It was worth noting that YEATS4 also played a crucial role in the malignant progression of GBM. For instance, studies showed that miR-203 could affect the proliferation and migration of GBM cells by regulating the GAS41/miR-10b axis (87). Therefore, it was necessary for us to further explore the specific function of YEATS4 in GBM and verify its specific role in regulating the proliferation and malignant progression of GBM cells through in vitro experiments.
The experimental results demonstrated that knockdown of YEATS4 in two GBM cell lines, LN229 and A1207, led to significant suppression of multiple malignant phenotypes. Specifically, YEATS4 knockdown resulted in decreased cell viability, markedly reduced colony-forming capacity, impaired cell proliferation, as well as inhibited migratory and invasive abilities in GBM cells. Collectively, these experimental findings indicated that YEATS4 played an essential contributor in regulating the proliferation, migration and invasion of GBM cells, potentially promoting tumor malignancy through modulation of these key biological processes. Therefore, YEATS4 might represent a potential therapeutic target for GBM treatment.
To establish a reliable molecular marker for predicting the prognosis of GBM, we formulated a prognostic risk score model founded on C2 PCLAF+ glioma cells, which effectively distinguished different patient groups with different prognosis risks. Based on this, we further investigated the immune microenvironment characteristics of GBM, and the differences in the composition of immune cells reflected the heterogeneity of the TME between the high and low PGRS groups. The important point is that a strong positive correlation was observed between the key prognostic genes (IRF7 and Macrophages M1, and IKZF3 and B cells naive), suggesting that these genes may regulate the functions of specific immune cells (such as macrophage polarization and B cell activation status) to shape the immunosuppressive or activating phenotype of the TME (92–94). Genomic analysis indicated that the high PGRS group exhibited a more active somatic mutation landscape and more extensive CNVs, which might influence their regulatory functions through dosage effects. This not only drives the malignant progression of tumors but may also indirectly shape the immunosuppressive TME by influencing cytokine secretion or interferon signaling pathways. Evaluation of immunotherapy response further supported the potential sensitivity of the high PGRS group to immunotherapy, reflecting a stronger immune response capacity. Analysis of immune checkpoint genes reinforced these findings, with several immune checkpoint-related genes (e.g., CD276, CD44, and C10orf54) significantly upregulated in the high PGRS group, potentially promoting immune evasion and influencing immunotherapy efficacy.
To comprehensively characterize the molecular differences across PGRS stratification groups, we analyzed DEGs and their functional enrichment profiles. Of particular note was the marked activation of the IL-17 signaling pathway (95, 96), in agreement with the pro-inflammatory microenvironment observed in the high PGRS group. Furthermore, the drug sensitivity analysis provided direct clues for individualized treatment based on the PGRS. In summary, the PGRS scoring system was not only a prognostic tool, but also its underlying immune cell map and the strong correlation between key genes and immune cells profoundly revealed the inherent immunosuppressive characteristics of GBM. This immunosuppressive microenvironment was a key barrier that drove tumor progression and treatment resistance. Our findings emphasized the importance of combined treatment strategies: that is, while directly killing tumor cells, it was necessary to actively intervene in the TME, reverse immunosuppression (such as targeting immunosuppressive cells or their secreted factors), release the suppressed anti-tumor immune response, so as to effectively improve treatment efficacy and overcome drug resistance.
This study thoroughly explored the unique molecular characteristics and spatial distribution patterns of the C2 PCLAF+ subtype in GBM and clarified its critical role in tumor progression, metabolic reprogramming, and cell death regulation. More importantly, based on the potential mechanisms of this subtype in GBM heterogeneity and treatment resistance, we further discussed potential clinical translational directions, providing important insights for the development of novel therapeutic strategies. Specifically, we proposed subtype-specific targeted therapies for C2 PCLAF+, particularly targeting its metabolic features (such as the development of glycolysis inhibitors) and key regulatory factor YEATS4. Meanwhile, the predictive models constructed by integrating ST data have shown potential clinical application value. They can not only be used to predict patients’ sensitivity to specific treatment regimens but also provide new technical means for real-time monitoring of the dynamic evolution of the TME during treatment. Based on this, we suggested a combined therapeutic strategy to block the MDK-LRP1 signaling axis, aimed at disrupting malignant cell-niche interplay. Finally, we proposed a personalized treatment strategy based on PGRS stratification. For instance, in clinical practice, the model can be integrated into the preoperative assessment and follow-up procedures, and combined with traditional imaging and histopathological information. Additionally, efforts can be made to enhance the coupling development of drug combination strategies and immune microenvironment regulation methods. For example, high PGRS patients can be included in the combined treatment plans of immunostimulants or metabolic intervention drugs to break through the current bottleneck of treatment resistance.
However, this research had certain limitations that merit consideration. Although the study highlighted the critical role of the C2 PCLAF+ subtype, other subtypes might also play important roles at different stages of GBM or in different microenvironments. Future research will need to further explore the functional characteristics of other subtypes to furnish a broader perspective on GBM heterogeneity. Additionally, the ST analysis was conducted using only two tissue sections, which limited the universality of the spatial localization results of the C2 PCLAF+ subtype. Future studies should expand the spatial dataset. Some of the findings in this study were validated through in vitro experiments; however, in vitro models might not fully replicate the complex in vivo TME, and key rescue experiments (such as overexpression of YEATS4 after siRNA treatment) are absent to confirm the causal relationship of the PCLAF+/YEATS4 axis. Moreover, the downstream targets of YEATS4 (such as PCLAF or cell cycle regulators) still required further validation at the protein level. The tumor immune microenvironment was analyzed using only computational deconvolution methods, and while these methods were feasible, the reliability of the immune feature analysis would be significantly enhanced if they could be verified through experimental means. The complexity of the tumor immune microenvironment was another area that needed further exploration (97–100), particularly the interactions between immune cells, immune evasion mechanisms, and immune tolerance, all of which remained to be investigated.
Conclusion
This study employed scRNA-seq and ST technologies to systematically analyze the molecular characteristics and spatial heterogeneity distribution patterns of the C2 PCLAF+ subtype in GBM. Through integrated multi-omics analysis, we found that the C2 PCLAF+ subtype exhibited significantly enhanced proliferative activity, unique metabolic reprogramming features, and aberrant cell death regulation mechanisms, all of which contributed to rapid tumor growth and the development of a drug-resistant phenotype. Notably, in vitro functional experiments confirmed the critical role of YEATS4 in promoting GBM progression, providing a potential new therapeutic target for targeted treatment strategies. Based on the molecular characteristics of the C2 PCLAF+ subtype, we successfully constructed an effective prognostic risk score model and integrated genomic, immune microenvironment, and drug sensitivity data. These findings offered important theoretical insights and practical guidance for a fuller comprehension of the molecular pathogenesis of GBM and the development of personalized treatment strategies.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author/s.
Author contributions
SM: Writing – original draft, Software, Formal analysis, Project administration, Methodology, Conceptualization, Validation, Investigation, Writing – review & editing. YS: Writing – review & editing, Methodology, Writing – original draft, Validation, Investigation, Conceptualization. SZ: Writing – review & editing, Methodology, Writing – original draft, Supervision, Investigation, Software, Data curation, Conceptualization. YF: Formal analysis, Data curation, Methodology, Visualization, Project administration, Software, Investigation, Validation, Supervision, Writing – review & editing, Conceptualization. LW: Supervision, Methodology, Investigation, Conceptualization, Software, Writing – review & editing, Project administration, Visualization, Formal analysis, Validation, Data curation. DL: Formal analysis, Validation, Methodology, Writing – review & editing, Supervision, Project administration, Conceptualization, Investigation, Software, Visualization. HJ: Supervision, Methodology, Conceptualization, Software, Project administration, Investigation, Visualization, Formal analysis, Validation, Writing – review & editing, Data curation. XZ: Investigation, Writing – review & editing, Visualization, Formal analysis, Validation, Conceptualization, Supervision, Project administration, Methodology. XL: Formal analysis, Writing – review & editing, Project administration, Data curation, Methodology, Visualization, Investigation, Supervision, Conceptualization, Validation. DY: Funding acquisition, Resources, Supervision, Project administration, Validation, Visualization, Writing – review & editing. DC: Writing – review & editing, Supervision, Funding acquisition, Visualization, Validation, Resources, Project administration. ZY: Project administration, Validation, Visualization, Funding acquisition, Supervision, Writing – review & editing, Resources.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. Thanks for the funding of the National Natural Science Foundation of China (82401623), China Postdoctoral Science Foundation (2024M752952), Henan Natural Science Foundation (242300421481), Henan Provincial Key Research and Development Special Project (251111313300) and the Henan Medical Science and Technique Research Project (LHGJ20230239, SBGJ202402043).
Acknowledgments
We sincerely appreciated the assistance provided by Figdraw (Image ID: WTRRW3ccd8) for Figure 1 in our study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1614549/full#supplementary-material
Supplementary Figure 1 | Spatial transcriptome atlas. (A, B) ST feature maps showed the spatial expression pattern of the C2 PCLAF+ subtype and PCLAF gene on ST 2 slide. (C) ST feature maps revealed spatial distributions of nCount-Spatial, nFeature-Spatial, G2/M.Score, and S.Score across ST 2 slide.
References
1. Lee IY, Hanft S, Schulder M, Judy KD, Wong ET, Elder JB, et al. Autologous cell immunotherapy (IGV-001) with IGF-1R antisense oligonucleotide in newly diagnosed glioblastoma patients. Future Oncol. (2024) 20:579–91. doi: 10.2217/fon-2023-0702
2. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2012-2016. Neuro Oncol. (2019) 21:v1–100. doi: 10.1093/neuonc/noz150
3. Ostrom QT, Patil N, Cioffi G, Waite K, Kruchko C, and Barnholtz-Sloan JS. Corrigendum to: CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2013-2017. Neuro Oncol. (2022) 24:1214. doi: 10.1093/neuonc/noaa269
4. Olar A and Aldape KD. Using the molecular classification of glioblastoma to inform personalized treatment. J Pathol. (2014) 232:165–77. doi: 10.1002/path.4282
5. Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The epidemiology of glioma in adults: a “state of the science” review. Neuro Oncol. (2014) 16:896–913. doi: 10.1093/neuonc/nou087
6. Ahmady F, Sharma A, Achuthan AA, Kannourakis G, and Luwor RB. The role of TIM-3 in glioblastoma progression. Cells. (2025) 14:346. doi: 10.3390/cells14050346
7. Anghileri E, Di Ianni N, Paterra R, Langella T, Zhao J, Eoli M, et al. High tumor mutational burden and T-cell activation are associated with long-term response to anti-PD1 therapy in Lynch syndrome recurrent glioblastoma patient. Cancer Immunol Immunother. (2021) 70:831–42. doi: 10.1007/s00262-020-02769-4
8. Liau LM, Ashkan K, Brem S, Campian JL, Trusheim JE, Iwamoto FM, et al. Association of autologous tumor lysate-loaded dendritic cell vaccination with extension of survival among patients with newly diagnosed and recurrent glioblastoma: A phase 3 prospective externally controlled cohort trial. JAMA Oncol. (2023) 9:112–21. doi: 10.1001/jamaoncol.2022.5370
9. Youssef G, Rahman R, Bay C, Wang W, Lim-Fat MJ, Arnaout O, et al. Evaluation of standard response assessment in neuro-oncology, modified response assessment in neuro-oncology, and immunotherapy response assessment in neuro-oncology in newly diagnosed and recurrent glioblastoma. J Clin Oncol. (2023) 41:3160–71. doi: 10.1200/JCO.22.01579
10. Campos B, Olsen LR, Urup T, and Poulsen HS. A comprehensive profile of recurrent glioblastoma. Oncogene. (2016) 35:5819–25. doi: 10.1038/onc.2016.85
11. Krex D, Bartmann P, Lachmann D, Hagstotz A, Jugel W, Schneiderman RS, et al. Aurora B kinase inhibition by AZD1152 concomitant with tumor treating fields is effective in the treatment of cultures from primary and recurrent glioblastomas. Int J Mol Sci. (2023) 24:5016. doi: 10.3390/ijms24055016
12. Ou A, Yung W, and Majd N. Molecular mechanisms of treatment resistance in glioblastoma. Int J Mol Sci. (2020) 22:351. doi: 10.3390/ijms22010351
13. Ma C, Qiu O, Huang C, Huang J, and Qu S. Roles of tenascin-XB in the glioma immune microenvironment. Bio Integration. (2023) 4:18–29. doi: 10.15212/bioi-2022-0014
14. Chen RQ, Liu F, Qiu XY, and Chen XQ. The prognostic and therapeutic value of PD-L1 in glioma. Front Pharmacol. (2018) 9:1503. doi: 10.3389/fphar.2018.01503
15. Hernandez A, Domenech M, Munoz-Marmol AM, Carrato C, and Balana C. Glioblastoma: relationship between metabolism and immunosuppressive microenvironment. Cells. (2021) 10:3529. doi: 10.3390/cells10123529
16. Thiruvengadam R, Dareowolabi BO, Moon EY, and Kim JH. Nanotherapeutic strategy against glioblastoma using enzyme inhibitors. BioMed Pharmacother. (2024) 181:117713. doi: 10.1016/j.biopha.2024.117713
17. Lakshmi BA and Kim YJ. Modernistic and emerging developments of nanotechnology in glioblastoma-targeted theranostic applications. Int J Mol Sci. (2022) 23:1641. doi: 10.3390/ijms23031641
18. Rosen J, Blau T, Grau SJ, Barbe MT, Fink GR, and Galldiks N. Extracranial metastases of a cerebral glioblastoma: A case report and review of the literature. Case Rep Oncol. (2018) 11:591–600. doi: 10.1159/000492111
19. Niu X, Li G, Kahlert UD, Ding L, Zheng J, Li C, et al. Integrative disulfidptosis-based risk assessment for prognostic stratification and immune profiling in glioma. J Cell Mol Med. (2025) 29:e70429. doi: 10.1111/jcmm.70429
20. Darmanis S, Sloan SA, Croote D, Mignardi M, Chernikova S, Samghababi P, et al. Single-cell RNA-Seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Rep. (2017) 21:1399–410. doi: 10.1016/j.celrep.2017.10.030
21. Tang W, Sun G, Ji GW, Feng T, Zhang Q, Cao H, et al. Single-cell RNA-sequencing atlas reveals an FABP1-dependent immunosuppressive environment in hepatocellular carcinoma. J Immunother Cancer. (2023) 11:e007030. doi: 10.1136/jitc-2023-007030
22. Traynor S, Jakobsen MK, Green TM, Komic H, Palarasah Y, Pedersen CB, et al. Single-cell sequencing unveils extensive intratumoral heterogeneity of cancer/testis antigen expression in melanoma and lung cancer. J Immunother Cancer. (2024) 12:e008759. doi: 10.1136/jitc-2023-008759
23. Gu X, Cai L, Luo Z, Shi L, Peng Z, Sun Y, et al. Identification and validation of a muscle failure index to predict prognosis and immunotherapy in lung adenocarcinoma through integrated analysis of bulk and single-cell RNA sequencing data. Front Immunol. (2022) 13:1057088. doi: 10.3389/fimmu.2022.1057088
24. Lin Z, Wang F, Yin R, Li S, Bai Y, Zhang B, et al. Single-cell RNA sequencing and immune microenvironment analysis reveal PLOD2-driven Malignant transformation in cervical cancer. Front Immunol. (2024) 15:1522655. doi: 10.3389/fimmu.2024.1522655
25. Lin Z, Sui X, Jiao W, Chen C, Zhang X, and Zhao J. Mechanism investigation and experiment validation of capsaicin on uterine corpus endometrial carcinoma. Front Pharmacol. (2022) 13:953874. doi: 10.3389/fphar.2022.953874
26. Zhao Z, Dong Y, Zhao Z, Xiahou Z, and Sun C. Single-cell atlas of endothelial cells in atherosclerosis: identifying C1 CXCL12+ ECs as key proliferative drivers for immunological precision therapeutics in atherosclerosis. Front Immunol. (2025) 16:1569988. doi: 10.3389/fimmu.2025.1569988
27. Zhao F, Jiang X, Li Y, Huang T, Xiahou Z, Nie W, et al. Characterizing tumor biology and immune microenvironment in high-grade serous ovarian cancer via single-cell RNA sequencing: insights for targeted and personalized immunotherapy strategies. Front Immunol. (2024) 15:1500153. doi: 10.3389/fimmu.2024.1500153
28. Hou M, Zhao Z, Li S, Zhang Z, Li X, Zhang Y, et al. Single-cell analysis unveils cell subtypes of acral melanoma cells at the early and late differentiation stages. J Cancer. (2025) 16:898–916. doi: 10.7150/jca.102045
29. Wang J, Zhao F, Zhang Q, Sun Z, Xiahou Z, Wang C, et al. Unveiling the NEFH+ Malignant cell subtype: Insights from single-cell RNA sequencing in prostate cancer progression and tumor microenvironment interactions. Front Immunol. (2024) 15:1517679. doi: 10.3389/fimmu.2024.1517679
30. Zhao Z, Ding Y, Tran LJ, Chai G, and Lin L. Innovative breakthroughs facilitated by single-cell multi-omics: manipulating natural killer cell functionality correlates with a novel subcategory of melanoma cells. Front Immunol. (2023) 14:1196892. doi: 10.3389/fimmu.2023.1196892
31. Zhao Z, Zhao Z, Lin Z, Fan L, Xiahou Z, Dong Y, et al. Decoding multiple myeloma: single-cell insights into tumor heterogeneity, immune dynamics, and disease progression. Front Immunol. (2025) 16:1584350. doi: 10.3389/fimmu.2025.1584350
32. Li H, Bian Y, Xiahou Z, Zhao Z, Zhao F, and Zhang Q. The cellular signaling crosstalk between memory B cells and tumor cells in nasopharyngeal carcinoma cannot be overlooked: Their involvement in tumor progression and treatment strategy is significant. J Cancer. (2025) 16:288–314. doi: 10.7150/jca.101420
33. Feng X, Luo Z, Zhang W, Wan R, Chen Y, Li F, et al. Zn-DHM nanozymes enhance muscle regeneration through ROS scavenging and macrophage polarization in volumetric muscle loss revealed by single-cell profiling. Adv Funct Mater. (2025), 2506476. doi: 10.1002/adfm.202506476
34. Jiang A, Bao Y, Wang A, Gong W, Gan X, Wang J, et al. Establishment of a prognostic prediction and drug selection model for patients with clear cell renal cell carcinoma by multiomics data analysis. Oxid Med Cell Longev. (2022) 2022:3617775. doi: 10.1155/2022/3617775
35. Zhang P, Wang D, Zhou G, Jiang S, Zhang G, Zhang L, et al. Novel post-translational modification learning signature reveals B4GALT2 as an immune exclusion regulator in lung adenocarcinoma. J Immunother Cancer. (2025) 13:e010787. doi: 10.1136/jitc-2024-010787
36. Li X, Lin Z, Zhao F, Huang T, Fan W, Cen L, et al. Unveiling the cellular landscape: insights from single-cell RNA sequencing in multiple myeloma. Front Immunol. (2024) 15:1458638. doi: 10.3389/fimmu.2024.1458638
37. Sun L, Shao W, Lin Z, Lin J, Zhao F, and Yu J. Single-cell RNA sequencing explored potential therapeutic targets by revealing the tumor microenvironment of neuroblastoma and its expression in cell death. Discov Oncol. (2024) 15:409. doi: 10.1007/s12672-024-01286-5
38. Zhao F, Hong J, Zhou G, Huang T, Lin Z, Zhang Y, et al. Elucidating the role of tumor-associated ALOX5+ mast cells with transformative function in cervical cancer progression via single-cell RNA sequencing. Front Immunol. (2024) 15:1434450. doi: 10.3389/fimmu.2024.1434450
39. Ni G, Sun Y, Jia H, Xiahou Z, Li Y, Zhao F, et al. MAZ-mediated tumor progression and immune evasion in hormone receptor-positive breast cancer: Targeting tumor microenvironment and PCLAF+ subtype-specific therapy. Transl Oncol. (2025) 52:102280. doi: 10.1016/j.tranon.2025.102280
40. Lin L, Zou J, Pei S, Huang W, Zhang Y, Zhao Z, et al. Germinal center B-cell subgroups in the tumor microenvironment cannot be overlooked: Their involvement in prognosis, immunotherapy response, and treatment resistance in head and neck squamous carcinoma. Heliyon. (2024) 10:e37726. doi: 10.1016/j.heliyon.2024.e37726
41. Zhang Y, Zhao Z, Huang W, Kim BS, Lin L, Li X, et al. Pan-cancer single-cell analysis revealing the heterogeneity of cancer-associated fibroblasts in skin tumors. Curr Gene Ther. (2024) 25:e15665232331353. doi: 10.2174/0115665232331353240911080642
42. Jin W, Zhang Y, Zhao Z, and Gao M. Developing targeted therapies for neuroblastoma by dissecting the effects of metabolic reprogramming on tumor microenvironments and progression. Theranostics. (2024) 14:3439–69. doi: 10.7150/thno.93962
43. Nie W, Zhao Z, Liu Y, Wang Y, Zhang J, Hu Y, et al. Integrative single-cell analysis of cardiomyopathy identifies differences in cell stemness and transcriptional regulatory networks among fibroblast subpopulations. Cardiol Res Pract. (2024) 2024:3131633. doi: 10.1155/2024/3131633
44. Huang W, Kim BS, Zhang Y, Lin L, Chai G, and Zhao Z. Regulatory T cells subgroups in the tumor microenvironment cannot be overlooked: Their involvement in prognosis and treatment strategy in melanoma. Environ Toxicol. (2024) 39:4512–30. doi: 10.1002/tox.24247
45. Shao W, Lin Z, Xiahou Z, Zhao F, Xu J, Liu X, et al. Single-cell RNA sequencing reveals that MYBL2 in Malignant epithelial cells is involved in the development and progression of ovarian cancer. Front Immunol. (2024) 15:1438198. doi: 10.3389/fimmu.2024.1438198
46. Zhou W, Lin Z, and Tan W. Deciphering the molecular landscape: integrating single-cell transcriptomics to unravel myofibroblast dynamics and therapeutic targets in clear cell renal cell carcinomas. Front Immunol. (2024) 15:1374931. doi: 10.3389/fimmu.2024.1374931
47. Lin Z, Li X, Shi H, Cao R, Zhu L, Dang C, et al. Decoding the tumor microenvironment and molecular mechanism: unraveling cervical cancer subpopulations and prognostic signatures through scRNA-Seq and bulk RNA-seq analyses. Front Immunol. (2024) 15:1351287. doi: 10.3389/fimmu.2024.1351287
48. Ge Q, Zhao Z, Li X, Yang F, Zhang M, Hao Z, et al. Deciphering the suppressive immune microenvironment of prostate cancer based on CD4+ regulatory T cells: Implications for prognosis and therapy prediction. Clin Transl Med. (2024) 14:e1552. doi: 10.1002/ctm2.1552
49. Ding Y, Zhao Z, Cai H, Zhou Y, Chen H, Bai Y, et al. Single-cell sequencing analysis related to sphingolipid metabolism guides immunotherapy and prognosis of skin cutaneous melanoma. Front Immunol. (2023) 14:1304466. doi: 10.3389/fimmu.2023.1304466
50. He Y, Luo Z, Nie X, Du Y, Sun R, Sun J, et al. An injectable multi-functional composite bioactive hydrogel for bone regeneration via immunoregulatory and osteogenesis effects. Adv Compos Hybrid Mater. (2025) 8:128. doi: 10.1007/s42114-025-01213-4
51. Zhao ZJ, Wei DP, Zheng RZ, Peng T, Xiao X, and Li FS. The gene coexpression analysis identifies functional modules dynamically changed after traumatic brain injury. Comput Math Methods Med. (2021) 2021:5511598. doi: 10.1155/2021/5511598
52. Chen F, Wang Y, He J, Chen L, Xue G, Zhao Y, et al. Molecular mechanisms of spawning habits for the adaptive radiation of endemic East Asian cyprinid fishes. Res (Wash D C). (2022) 2022:9827986. doi: 10.34133/2022/9827986
53. Zhao Z, Li T, Dong X, Wang X, Zhang Z, Zhao C, et al. Untargeted metabolomic profiling of cuprizone-induced demyelination in mouse corpus callosum by UPLC-orbitrap/MS reveals potential metabolic biomarkers of CNS demyelination disorders. Oxid Med Cell Longev. (2021) 2021:7093844. doi: 10.1155/2021/7093844
54. Zhao ZJ, Zheng RZ, Wang XJ, Li TQ, Dong XH, Zhao CY, et al. Integrating lipidomics and transcriptomics reveals the crosstalk between oxidative stress and neuroinflammation in central nervous system demyelination. Front Aging Neurosci. (2022) 14:870957. doi: 10.3389/fnagi.2022.870957
55. Wang Y, Zhao ZJ, Kang XR, Bian T, Shen ZM, Jiang Y, et al. lncRNA DLEU2 acts as a miR-181a sponge to regulate SEPP1 and inhibit skeletal muscle differentiation and regeneration. Aging (Albany NY). (2020) 12:24033–56. doi: 10.18632/aging.104095
56. Li XY, Zhao ZJ, Wang JB, Shao YH, Hui-Liu, You JX, et al. m7G methylation-related genes as biomarkers for predicting overall survival outcomes for hepatocellular carcinoma. Front Bioeng Biotechnol. (2022) 10:849756. doi: 10.3389/fbioe.2022.849756
57. Chen L, He Y, Duan M, Yang T, Chen Y, Wang B, et al. Exploring NUP62’s role in cancer progression, tumor immunity, and treatment response: insights from multi-omics analysis. Front Immunol. (2025) 16:1559396. doi: 10.3389/fimmu.2025.1559396
58. Du H, Li S, Lu J, Tang L, Jiang X, He X, et al. Single-cell RNA-seq and bulk-seq identify RAB17 as a potential regulator of angiogenesis by human dermal microvascular endothelial cells in diabetic foot ulcers. Burns Trauma. (2023) 11:tkad020. doi: 10.1093/burnst/tkad020
59. Huo Y, Shao S, Liu E, Li J, Tian Z, Wu X, et al. Subpathway analysis of transcriptome profiles reveals new molecular mechanisms of acquired chemotherapy resistance in breast cancer. Cancers (Basel). (2022) 14:4878. doi: 10.3390/cancers14194878
60. Xing J, Cai H, Lin Z, Zhao L, Xu H, Song Y, et al. Examining the function of macrophage oxidative stress response and immune system in glioblastoma multiforme through analysis of single-cell transcriptomics. Front Immunol. (2023) 14:1288137. doi: 10.3389/fimmu.2023.1288137
61. Liu P, Xing N, Xiahou Z, Yan J, Lin Z, and Zhang J. Unraveling the intricacies of glioblastoma progression and recurrence: insights into the role of NFYB and oxidative phosphorylation at the single-cell level. Front Immunol. (2024) 15:1368685. doi: 10.3389/fimmu.2024.1368685
62. Zhang P, Yang Z, Liu Z, Zhang G, Zhang L, Zhang Z, et al. Deciphering lung adenocarcinoma evolution: Integrative single-cell genomics identifies the prognostic lung progression associated signature. J Cell Mol Med. (2024) 28:e18408. doi: 10.1111/jcmm.18408
63. Lin Z, Fan W, Sui X, Wang J, and Zhao J. Necroptosis-related lncRNA signatures for prognostic prediction in uterine corpora endometrial cancer. Reprod Sci. (2023) 30:576–89. doi: 10.1007/s43032-022-01023-9
64. Lin Z, Sui X, Jiao W, Wang Y, and Zhao J. Exploring the mechanism and experimental verification of puerarin in the treatment of endometrial carcinoma based on network pharmacology and bioinformatics analysis. BMC Complement Med Ther. (2022) 22:150. doi: 10.1186/s12906-022-03623-z
65. Zou J, Lin Z, Jiao W, Chen J, Lin L, Zhang F, et al. A multi-omics-based investigation of the prognostic and immunological impact of necroptosis-related mRNA in patients with cervical squamous carcinoma and adenocarcinoma. Sci Rep. (2022) 12:16773. doi: 10.1038/s41598-022-20566-0
66. Tang JY, Feng XQ, Huang XX, Zhang YP, Guo ZT, Chen L, et al. Development and validation of a predictive model for patients with post-extubation dysphagia. World J Emerg Med. (2023) 14:49–55. doi: 10.5847/wjem.j.1920-8642.2023.021
67. Zheng RZ, Zhao ZJ, Yang XT, Jiang SW, Li YD, Li WJ, et al. Initial CT-based radiomics nomogram for predicting in-hospital mortality in patients with traumatic brain injury: a multicenter development and validation study. Neurol Sci. (2022) 43:4363–72. doi: 10.1007/s10072-022-05954-8
68. Lin Z, Zou J, Sui X, Yao S, Lin L, Wang J, et al. Necroptosis-related lncRNA signature predicts prognosis and immune response for cervical squamous cell carcinoma and endocervical adenocarcinomas. Sci Rep. (2022) 12:16285. doi: 10.1038/s41598-022-20858-5
69. Zhao J, Jiao W, Sui X, Zou J, Wang J, and Lin Z. Construction of a prognostic model of luteolin for endometrial carcinoma. Am J Transl Res. (2023) 15:2122–39.
70. Zhao ZJ, Chen D, Zhou LY, Sun ZL, Wang BC, and Feng DF. Prognostic value of different computed tomography scoring systems in patients with severe traumatic brain injury undergoing decompressive craniectomy. J Comput Assist Tomogr. (2022) 46:800–07. doi: 10.1097/RCT.0000000000001343
71. Jiang A, Ye J, Zhou Y, Zhu B, Lu J, Ge S, et al. Copper death inducer, FDX1, as a prognostic biomarker reshaping tumor immunity in clear cell renal cell carcinoma. Cells. (2023) 12:349. doi: 10.3390/cells12030349
72. Ding J, Chen J, Zhou J, Jiang Z, Xiang D, and Xing W. Association between renal surface nodularity and increased adverse vascular event risk in patients with arterial hypertension. Clin Exp Hypertens. (2023) 45:2228518. doi: 10.1080/10641963.2023.2228518
73. Zhao J, Zou J, Jiao W, Lin L, Wang J, and Lin Z. Construction of N-7 methylguanine-related mRNA prognostic model in uterine corpus endometrial carcinoma based on multi-omics data and immune-related analysis. Sci Rep. (2022) 12:18813. doi: 10.1038/s41598-022-22879-6
74. He Z, Liu Z, Wang Q, Sima X, Zhao W, He C, et al. Single-cell and spatial transcriptome assays reveal heterogeneity in gliomas through stress responses and pathway alterations. Front Immunol. (2024) 15:1452172. doi: 10.3389/fimmu.2024.1452172
75. Qiu C, Wang W, Xu S, Li Y, Zhu J, Zhang Y, et al. Construction and validation of a hypoxia-related gene signature to predict the prognosis of breast cancer. BMC Cancer. (2024) 24:402. doi: 10.1186/s12885-024-12182-0
76. Kang DW, Hwang WC, Noh YN, Park KS, and Min DS. Phospholipase D1 inhibition sensitizes glioblastoma to temozolomide and suppresses its tumorigenicity. J Pathol. (2020) 252:304–16. doi: 10.1002/path.5519
77. DeCordova S, Shastri A, Tsolaki AG, Yasmin H, Klein L, Singh SK, et al. Molecular heterogeneity and immunosuppressive microenvironment in glioblastoma. Front Immunol. (2020) 11:1402. doi: 10.3389/fimmu.2020.01402
78. Wang Z, Dai Z, Zhang H, Zhang N, Liang X, Peng L, et al. Comprehensive analysis of pyroptosis-related gene signatures for glioblastoma immune microenvironment and target therapy. Cell Prolif. (2023) 56:e13376. doi: 10.1111/cpr.13376
79. Xu S, Liang J, Shen T, Zhang D, and Lu Z. Causal links between immune cells and asthma: Insights from a Mendelian Randomization analysis. J Asthma. (2025) 62:346–53. doi: 10.1080/02770903.2024.2403740
80. Duan Y, Wu Y, Tian J, Yin Y, Yuan Z, Zhu W, et al. Elucidation of the mechanism Underlying the promotion of ferroptosis and enhanced antitumor immunity by citrus polymethoxyflavones in CRC cells. Front Pharmacol. (2025) 16:1571178. doi: 10.3389/fphar.2025.1571178
81. Luo ZW, Sun YY, Xia W, Xu JY, Xie DJ, Jiao CM, et al. Physical exercise reverses immuno-cold tumor microenvironment via inhibiting SQLE in non-small cell lung cancer. Mil Med Res. (2023) 10:39. doi: 10.1186/s40779-023-00474-8
82. Dong J, Peng Y, Zhong M, Xie Z, Jiang Z, Wang K, et al. Implication of lncRNA ZBED3-AS1 downregulation in acquired resistance to Temozolomide and glycolysis in glioblastoma. Eur J Pharmacol. (2023) 938:175444. doi: 10.1016/j.ejphar.2022.175444
83. Wen S, Lv X, Li P, Li J, and Qin D. Analysis of cancer-associated fibroblasts in cervical cancer by single-cell RNA sequencing. Aging (Albany NY). (2023) 15:15340–59. doi: 10.18632/aging.205353
84. Dai H, Tao X, Shu Y, Liu F, Cheng X, Li X, et al. Integrating single-cell RNA-Seq and bulk RNA-Seq data to explore the key role of fatty acid metabolism in hepatocellular carcinoma. Sci Rep. (2025) 15:2077. doi: 10.1038/s41598-025-85506-0
85. Fu Y, Li S, Zhao Y, Zhang X, Mao X, and Xu R. Integrative single-cell and spatial transcriptomics analysis reveals MDK-NCL pathway’s role in shaping the immunosuppressive environment of lung adenocarcinoma. Front Immunol. (2025) 16:1546382. doi: 10.3389/fimmu.2025.1546382
86. Yuan F, Cai X, Cong Z, Wang Y, Geng Y, Aili Y, et al. Roles of the m(6)A modification of RNA in the glioblastoma microenvironment as revealed by single-cell analyses. Front Immunol. (2022) 13:798583. doi: 10.3389/fimmu.2022.798583
87. Pal D, Mukhopadhyay D, Ramaiah MJ, Sarma P, Bhadra U, and Bhadra MP. Regulation of Cell Proliferation and Migration by miR-203 via GAS41/miR-10b Axis in Human Glioblastoma Cells. PloS One. (2016) 11:e0159092. doi: 10.1371/journal.pone.0159092
88. Wang Z, Zhao N, Zhang S, Wang D, Wang S, and Liu N. YEATS domain-containing protein GAS41 regulates nuclear shape by working in concert with BRD2 and the mediator complex in colorectal cancer. Pharmacol Res. (2024) 206:107283. doi: 10.1016/j.phrs.2024.107283
89. Tao K, Yang J, Hu Y, and Deng A. Knockdown of YEATS4 inhibits colorectal cancer cell proliferation and induces apoptosis. Am J Transl Res. (2015) 7:616–23.
90. Han S, Cao C, Liu R, Yuan Y, Pan L, Xu M, et al. GAS41 mediates proliferation and GEM chemoresistance via H2A.Z.2 and Notch1 in pancreatic cancer. Cell Oncol (Dordr). (2022) 45:429–46. doi: 10.1007/s13402-022-00675-8
91. Pikor LA, Lockwood WW, Thu KL, Vucic EA, Chari R, Gazdar AF, et al. YEATS4 is a novel oncogene amplified in non-small cell lung cancer that regulates the p53 pathway. Cancer Res. (2013) 73:7301–12. doi: 10.1158/0008-5472.CAN-13-1897
92. Yang S, Liu H, Zheng Y, Chu H, Lu Z, Yuan J, et al. The role of PLIN3 in prognosis and tumor-associated macrophage infiltration: A pan-cancer analysis. J Inflammation Res. (2025) 18:3757–77. doi: 10.2147/JIR.S509245
93. Tang J, He J, Guo H, Lin H, Li M, Yang T, et al. PTBP2-mediated alternative splicing of IRF9 controls tumor-associated monocyte/macrophage chemotaxis and repolarization in neuroblastoma progression. Res (Wash D C). (2023) 6:33. doi: 10.34133/research.0033
94. Xu Y, Zeng Y, Xiao X, Liu H, Zhou B, Luo B, et al. Targeted imaging of tumor associated macrophages in breast cancer. Bio Integration. (2023) 4:114–24. doi: 10.15212/bioi-2022-0010
95. Amatya N, Garg AV, and Gaffen SL. IL-17 signaling: the yin and the yang. Trends Immunol. (2017) 38:310–22. doi: 10.1016/j.it.2017.01.006
96. Li X, Bechara R, Zhao J, McGeachy MJ, and Gaffen SL. IL-17 receptor-based signaling and implications for disease. Nat Immunol. (2019) 20:1594–602. doi: 10.1038/s41590-019-0514-y
97. Shen J, Hu R, Lin A, Jiang A, Tang B, Liu Z, et al. Characterization of second primary Malignancies post CAR T-cell therapy: real-world insights from the two global pharmacovigilance databases of FAERS and VigiBase. EClinicalMedicine. (2024) 73:102684. doi: 10.1016/j.eclinm.2024.102684
98. Zhang P, Wen B, Gong J, Liu Z, Zhang M, Zhou G, et al. Clinical prognostication and immunotherapy response prediction in esophageal squamous cell carcinoma using the DNA damage repair-associated signature. Environ Toxicol. (2024) 39:2803–16. doi: 10.1002/tox.24155
99. Liang L, Zhu Y, Li J, Zeng J, Yuan G, and Wu L. Immune subtypes and immune landscape analysis of endometrial carcinoma. J Immunol. (2022) 209:1606–14. doi: 10.4049/jimmunol.2200329
Keywords: glioblastoma, scRNA-seq, spatial transcriptomics, YEATS4, tumor heterogeneity
Citation: Ma S, Sun Y, Zheng S, Fu Y, Wang L, Liu D, Jiao H, Zhu X, Li X, Yan D, Chen D and Ye Z (2025) Single-cell and spatial atlas of glioblastoma heterogeneity: characterizing the PCLAF+ subtype and YEATS4’s oncogenic role. Front. Immunol. 16:1614549. doi: 10.3389/fimmu.2025.1614549
Received: 19 April 2025; Accepted: 11 July 2025;
Published: 25 July 2025.
Edited by:
Jing Zhang, University of South Dakota, United StatesReviewed by:
Chen Li, Free University of Berlin, GermanyAimin Jiang, Fudan University, China
Shengshan Xu, Jiangmen Central Hospital, China
Kang Qin, University of South Dakota, United States
Copyright © 2025 Ma, Sun, Zheng, Fu, Wang, Liu, Jiao, Zhu, Li, Yan, Chen and Ye. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dongming Yan, bXJkbXlhbkAxNjMuY29t; Di Chen, bXJkaWNoZW5AMTYzLmNvbQ==; Zi Ye, eWV6aTAzMjZhQDE2My5jb20=
†These authors have contributed equally to this work