Single-cell and spatial transcriptome analyses reveal MAZ(+) NPC-like clusters as key role contributing to glioma recurrence and drug resistance

Abstract

Background

Temozolomide (TMZ) is an important chemotherapeutic agent for glioma treatment. However, the emergence of drug resistance impedes its application. Traditional population-level studies are limited in elucidating resistance mechanisms. Advances in single-cell and spatial transcriptomics technologies provide feasible resolution for studying the cellular composition and dynamics of tumors. In this study, we investigated the heterogeneity of gliomas associated with TMZ resistance at the single-cell and spatial transcriptome levels to identify the resistance mechanisms and potential therapeutic strategies.

Methods

Single cell sequencing technology was utilized to identify the cellular clusters of gliomas. Drug perturbation analysis and cellular propensity analysis revealed key cluster responding to TMZ. Enrichment analysis was preformed to explore the function of clusters. Transcription factor activity analysis revealed key transcription factors contributing to tumor resistance. Spatial transcriptome data and bulk RNA-seq data validates the role of key transcription factors. Downstream targets of key transcription factors were predicted and validated using gene regulation assays. Drug sensitivity analyses were used to seek viable strategies to overcome drug resistance.

Results

Glioma cells from before and after temozolomide treatment samples were classified into six clusters: NPC-like cluster, OPC-like cluster, MES-like cluster, AC-like cluster, OC, and Neuron. NPC-like clusters exhibited strong stemness and DNA repair capacity. The activity of MAZ in NPC-like cluster was significantly enhanced after TMZ treatment. The proportion of MAZ( +)_NPC-like cluster was higher in TMZ treated samples. Patients with high proportion of MAZ( +)_NPC-like cluster had poorer survival. Upregulation of MAZ is able to enhance drug resistance in glioma cells, but this phenomenon disappeared when FoxM1 expression was further silenced. The combination of paclitaxel and Trametinib is a promising strategy to overcome resistance.

Conclusions

NPC-like cluster is prevalent in recurrent and drug-resistant gliomas. MAZ transcription factors are critical regulators that promote the development of drug resistance in NPC-like clusters by enhancing the capacity of DNA repair and stemness. Patients with high proportions of MAZ( +)_NPC-like clusters have poor TMZ sensitivity and prognosis. MAZ enhances stemness and drug resistance in glioma cells by upregulating FOXM1 expression. The combination of paclitaxel and paclitaxel is a promising therapeutic strategy for treating gliomas and overcoming drug resistance.

Introduction

Glioma is the most prevalent malignant tumor in the central nervous system [1], characterized by the heterogeneity, invasiveness, and drug resistance, which poses great challenges to treatment [2]. Despite advancements made in recent years with surgery, radiotherapy and chemotherapy, the overall prognosis for glioma patients remains poor. Temozolomide (TMZ), as an alkylating agent, is capable of penetrating the blood–brain barrier and significantly improving patient survival, has become the standard chemotherapeutic agent for treating glioma in clinical practice [3]. However, as treatment proceeds, tumors inevitably develop resistance to temozolomide [4], which limits its clinical efficacy.

The development of drug resistance is closely associated with tumor heterogeneity [5]. Gliomas display intricate cellular architectures, with distinct cell clusters exhibiting obvious differences in gene expression patterns, signaling pathway activities, and drug responsiveness [6]. This heterogeneity not only increases the complexity of treatment strategies but also complicates the study of drug resistance mechanisms. Traditional population-level research often ignores differences at the single-cell and spatial levels [7], which makes researchers difficult to elucidate the underlying molecular mechanisms of drug resistance. In response, the methods based on single-cell and spatial transcriptomics sequencing technology provide researchers with effective tools that allow them to delve deeper into the cellular composition of tumors and dynamic changes [8].

In this study, we combine single-cell and spatial transcriptome sequencing technologies to comprehensively and systematically analyze the heterogeneity associated with temozolomide resistance in glioma. We explore the mechanisms of drug resistance and search for potential therapeutic targets. Finally, we provide reliable theoretical support and experimental basis for overcoming temozolomide resistance.

Materials and methods

Data sources

Single-cell and spatial transcriptome data for gliomas were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo). The GSE190129 dataset provides single cell sequencing data for gliomas, which includes the TMZ untreatment sample (GSM5714943) and corresponding treatment sample (GSM5714951 and GSM5714952). The spatial transcriptome sequencing data for gliomas were obtained from the GSE235672 dataset, where GSM7507316 and GSM7507321 contain both before and after temozolomide treatment tissue slices. The bulk RNA-sequencing data for glioma patients were accessed from the TCGA database [9] (https://portal.gdc.cancer.gov) and the CGGA database [10] (http://www.cgga.org.cn), which contains the CGGA_325, TCGA_LGG and TCGA_GBM datasets. The patients treated with temozolomide were screened for further analysis.

Processing of glioma single-cell data

Single-cell RNA sequencing data were processed and analyzed using the Seurat [11] package. First, UMI matrices were read to construct Seurat objects. Next, the data were quality controlled, normalized, batch differences removed, and clustered. Finally, the cell clusters were annotated with the help of literature and the CellMarker website (http://bio-bigdata.hrbmu.edu.cn/CellMarker).

Inference of cellular origin

Tumor cells are inferred using the InferCNV [12] package. Neuron cluster served as the normal reference cells. Clusters with significant changes in chromosome copy number (e.g., amplifications or deletions) were classified as tumor clusters.

Single-cell perturbation and drug sensitivity analysis

The Augur [13] package was used to assess the effects of TMZ perturbation on clusters, with AUC values indicating the degree of the effects. The beyondcell [14] package further predicts the directionality of perturbations. This package contains two feature sets: the drug perturbation feature set (PSc) and the drug sensitivity feature set (SSc), which can evaluate cell-drug response patterns. The BCscore and Switch Points are calculated by the Beyondcell algorithm. When the BCscore is greater than the Switch Points, it indicates that the cell x gene matrix is consistent with drug-induced transcriptional changes. When BCscore is less than Switch Points, it indicates that the cell x gene matrix is opposite to the transcriptional changes.

Single-cell stemness analysis and pseudotime analysis

The CytoTRACE [15] package was used to evaluate the developmental potential of single cells. The higher CytoTRACE score, the stronger the developmental potential and stemness in the cells. The Monocle [16] package was used for pseudotime analysis. The process of cellular differentiation over (pseudo) time was tracker by constructing trajectories of cellular change.

Single-cell transcription factor activity analysis

The pySCENIC [17] package was used to infer transcription factor (TF) activity in single-cell data. First, gene co-expression networks were constructed. Second, regulatory modules for transcription factors and their target genes were generated. Finally, the regulatory modules were scored using the AUCell algorithm to calculate transcription factor activity.

Single-cell clusters propensity analysis

The STARTRAC [18] package was used to evaluate single cell clusters propensity. The Ratio of observed to expected (Ro/e) values was calculated to observe clusters preference. Higher Ro/e values indicated a greater likelihood of cluster enrichment in a particular group.

Virtual knockout experiments

The scTenifoldKnk [19] package was used to perform virtual knockout experiments. The single-cell data were used for analysis. Genes that were significantly changed after target gene knockdown were acquired.

Enrichment analysis

GSEA analysis was performed using the clusterprofiler [20] package, selecting the c2.all.Hs.symbols gene set as the reference. Single-cell enrichment analysis was performed for each cell using the AddModuleScores function in the Seurat, followed by visualization via UMAP plots.

Spatial transcriptome analysis

Spatial transcriptome data were processed using the Seurat package. The UMI matrix and corresponding image slice data were read to construct Seurat objects. Data were normalized using SCTransform. Enrichment scoring was performed for each splot using the AddModuleScores function. The annotation results of cell type from single-cell data were integrated into the spatial data through the RCTD [21] package. The t-test was used to compare the difference in the proportion of cell types in the tissues before and after TMZ treatment. Correlation analysis was used to analyse the relationship between cluster abundance and function.

Inference of relative proportions of cell types in bulk RNA-seq for glioma patients

The relative proportion of cell types in bulk RNA-seq data was calculated through the BayesPrism [22] package. The annotation results of single-cell data were employed as a priori information to infer the distribution of different cell types in bulk RNA-seq data. The inference results were further validated by pathology information.

Gene correlation analysis

For single-cell data, the data were first normalized and then analyzed for correlation between the MAZ genes and all other genes using the ‘Pearson’ method. In the results, genes with p-values less than 0.05 and correlation coefficients greater than 0.35 were considered as MAZ high correlation genes. For the bulk RNA-seq data from TCGA and CGGA, TPM data were selected. The expression correlation between the MAZ gene and the FoxM1 gene was analyzed using the ‘Pearson’ method, and the results were visualized through scatter plots.

Survival analysis and COX regression analysis

Survival analysis and COX regression analysis were performed by the survival package. Patients treated with temozolomide were screened from the TCGA_LGG, TCGA_GBM, and CGGA_325 databases. Survival rates were calculated for patients in the MAZ(+)_NPC-like high and low groups. The results were shown as survival graphs.

Cell culture

Human glioblastoma cell line LN229 were obtained from the Pricella company. Cells were cultured in DMEM containing 10% fetal bovine serum and 1% penicillin/streptomycin at the condition of 37 °C with 5% CO2. When the cell density reached 90%, the cells were detached using 0.25% trypsin-EDTA and passaged at a ratio of 1:2 or 1:3.

Regulation of gene expression in glioma cells

The MAZ-overexpressing lentiviral was generated by co-transfecting 293 T cells with the lentiviral packaging plasmids and the MAZ gene transfer plasmid. After 8 hours, the medium was replaced with fresh medium. and cells were incubated for an additional 48 hours. The cell medium was collected, centrifuged, and filtered to obtain lentiviral particles. Lentivirus were used to infect glioma cells. After 12 hours of transfection, the medium was changed to fresh medium, and cells were cultured for another 48 hours. Uninfected cells were eliminated by treatment with 2 μg/mL puromycin. Finally, a steady overexpression cell line of MAZ was established. The FoxM1 gene of glioma cells was knocked down by treatment with FoxM1 siRNA (GGAAAUGCUUGUGAUUCAACA). For the rescue experiment, an siRNA-resistant FoxM1 overexpression plasmid was constructed. Synonymous mutations were designed within the siRNA-targeted region to disrupt complementarity with the siRNA (Mutation of GGAAATGCTTGTGATTCAAC sequence into GGAGATGCTCGTTATCCAGCA). For identification and detection of mutant proteins, FLAG tag sequences were inserted in the plasmid. Glioma cells were transfected with the mutated Flag-FoxM1 overexpression plasmid for 12 hours. After 12 hours of transfection, the culture was changed with fresh medium. The efficiency of overexpression and knockdown was detected using PCR and WB.

Reverse Transcription quantitative PCR (RT-qPCR)

The total RNA from glioma cells was extracted using the TRIzol reagent. The RNA was then reverse transcribed into cDNA using PrimeScript™ RT reagent Kit (Takara). Real-time quantitative PCR was performed with the SYBR Green Premix Ex Taq kit (Takara) under the following conditions: pre-denaturation at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 15 s and extension at 55 °C for 30 s. Relative gene expression level was calculated using the ΔΔCt method. primer sequence: MAZ: Forward: 5'-CTACAACTGCTCCCACT-3'; Reverse: 5'-TGGTGAAGCCTTTGTTG-3'. FoxM1: Forward: 5'-GCGCACGGCGGAAGATGAA-3'; Reverse: 5'- CCACTCTTCCAAGGGAGGGCTC −3'.

Western blot (WB)

Glioma Cells were lysed using RIPA buffer and then centrifuged at 12,000 × g for 15 mins to get the protein. The BCA assay was used to calculate protein concentrations. The protein samples were added to the gel tank for electrophoresis, and then the proteins transferred onto PVDF membranes at 300 mA for 1.5 hours using a wet transfer system. The membranes were blocked with 5% non-fat milk for 2 hours at room temperature. The membranes were incubated with primary antibodies for 12 hours at 4 °C, and secondary antibodies for 1 hour at room temperature. Finally, protein bands were visualized using enhanced chemiluminescence (ECL) reagents.

Temozolomide toxicity assay

Cells from different experimental groups were spread into 96-well plates with 5000 cells per well. After 24 hours of incubation, different concentrations of TMZ (0 μM, 400 μM, 1000 μM, 1600 μM, 2200 μM, 2800 μM) were used to treat the cells. After 48 hours of drug treatment, 10 μL of CCK8 was added to each well. After 1 hour, the absorbance at 450 nm was measured to calculate the cell viability. Drug-response curve was constructed and the IC50 of temozolomide was extrapolated.

Plate cloning assay

Glioma cells were seeded into 6-well plates with 80,000 cells per well. There are two groups in the experiment. Control group: Cells were cultured for 14 days without any drug treatment. TMZ group: After 24 hours of seeding, cells were treated with medium containing TMZ (50 μM). After 48 hours of TMZ exposure, the medium was replaced with fresh medium, and cells were incubated for an additional 11 days. After 14 days of culture, cells were fixed with methanol and stained with 0.1% crystalline violet.

Statistical analysis

All statistical analyses were performed based on R software (version: 4.4.1). Correlation analysis was performed through the cor function. The comparisons of two groups were performed using the t-test method, while comparisons of three or more groups were performed using the ANOVA method. The level of difference between groups was quantified using effect size. P-values less than 0.05 were considered statistically significant.

Results

NPC-like clusters are the dominant clusters for glioma temozolomide resistance and recurrence

Single-cell data of gliomas before and after temozolomide treatment were obtained from the GSE190129 dataset and were clustered into 11 subsets (Fig. 1A). These subsets were then further annotated into six clusters: AC-like cluster (Astrocyte like cluster), OPC-like cluster (Oligodendrocyte Progenitor like cluster), NPC-like cluster (Neural Progenitor like cluster), MES-like cluster (Mesenchymal like cluster), OC cluster (Oligodendrocyte cluster) and Neuron cluster (Fig. 1B)) based on their marker gene expression (Fig. 1C). Copy number variant (CNV) analysis results showed that the neuron and OC cluster exhibited no significant copy number alterations, whereas the AC-like, OPC-like, NPC-like, and MES-like clusters displayed chromosomal amplifications and deletions (Fig. 1D). Consequently, AC-like, OPC-like, NPC-like, and MES-like were classified as tumor cell types, while OC and Neuron were considered as normal cell types. The'Augur'package was employed to evaluate the perturbation efficacy of TMZ. The calculated results showed that the AUC values of AC-like, OPC-like, NPC-like, and MES-like clusters exceeded 0.9, thus indicating that TMZ exerted a significant perturbation effect on these clusters (Fig. 1E). Further analysis used the'beyondcell'package to explore the specific perturbations produced by TMZ. Resistance perturbation signatures of temozolomide were extracted from the Drug Perturbation Signatures collection (PSc). The analysis results showed that, in the control group, the BCScore of most cells was ranged between 0 and 0.5, which was opposite to the change of the resistance feature, suggesting drug sensitivity. In contrast, in the TMZ group, the BCScore of most cells was located between 0.5 and 1, which was consistent with the change of the resistance feature, suggesting drug resistance (Fig. 1F). The t-test comparing the difference in BCScore of tumor clusters before and after TMZ treatment showed that the BCScore of tumor clusters were all increased, with the most significant increasement was in the NPC-like cluster (Fig. 1G). Furthermore, cell propensity analysis revealed that the Ro/e index for NPC-like cluster was the higher in the resistant group and recurrent group in the GSE190129 and GBmap datasets, suggesting that NPC-like cluster are prevalent in the resistant and recurrent tumor (Fig. 1H). These findings demonstrate that all clusters develop drug resistance in TMZ group, with the strongest responses observed in the NPC-like cluster.

Fig. 1

A Glioma single cell were classified into 11 cell subsets. B These 11 subsets were further annotated into 6 clusters: AC-like cluster, OPC-like cluster, NPC-like cluster, MES-like cluster. C The marker genes of the 6 clusters were demonstrated through the bubble diagram. D Neuron cluster were as the normal reference cells for the inferCNV analysis. The OC cluster exhibited no significant copy number alterations, whereas the AC-like, OPC-like, NPC-like, and MES-like clusters displayed chromosomal amplifications and deletions. E The AUC values of OPC-like, AC-like, MES-like and NPC-like clusters were above 0.9. F In the control group, the BCScore of most cells rangedd between 0 and 0.5, indicating that this sample was sensitive to temozolomide, whereas in the TMZ group, the BCScore between 0.5 and 1, indicating resistance to temozolomide. G The BCScore of tumor clusters increased after temozolomide treatment, with the most significant increasement was in the NPC-like cluster. H The Ro/e index of NPC-like was highest in the resistant and recurrent group, indicating that NPC-like cluster tend to be enriched in the drug resistant and recurrent samples

The changes in the function of NPC-like cluster before and after TMZ treatment

The analysis using ‘CytoTRACE’ and ‘monocle’ packages showed the NPC-like cluster possessed a powerful degree of stemness (Fig. 2A) and were capable of differentiating into OPC-like, MES-like, and AC-like clusters along pseudotemporal trajectory (Fig. 2B). The differential genes for each cluster before and after TMZ treatment were calculated (Fig. 2C). GSEA analysis of the differential genes for the NPC-like cluster revealed that the PI3K pathway and Erbb4 pathway were activated, adipocytokine and glycosaminoglycan metabolism was enhanced, and cell stemness and intercellular communication were elevated after temozolomide treatment (Fig. 2D). Further functional exploration for the NPC-like cluster found that stemness features, such as Embryonic Stem Cell and Pediatric Cancer Markers, were significantly enriched, telomere extension and maintenance were reinforced, and multiple DNA repair functions, such as Mismatch Repair, Base Excision Repair, Base Excision Repair, Nucleotide Excision Repair, and Homology Directed Repair appeared, and various metabolic functions, such as Adipogenesis, Metabolism Of Nucleotides, Onecarbon Metabolism, and Folate Biosynthesis were enhanced, and the above functions were further strengthened in the TMZ group (Fig. 2E). The above results demonstrate that the NPC-like cluster possess a robust capacity for DNA repair and stemness, enabling resistance to TMZ-induced toxicity and facilitating tumor regeneration with enhanced drug resistance.

Fig. 2

A The NPC-like cluster exhibited a high degree of stemness. B NPC-like cluster is able to differentiate into MES-like cluster, OPC-like cluster, and AC-like cluster. C Volcano plots displayed differential genes of each cluster before and after temozolomide treatment. D GSEA revealed that after temozolomide treatment, the NPC-like activated the PI3K pathway and Erbb4 pathway, enhanced adipocytokine and glycosaminoglycan metabolism, and increased cell stemness and cell-to-cell communication. E The function enrichment of NPC-like cluster showed that stemness features, such as Embryonic Stem Cell and Pediatric Cancer Markers, were significantly enriched, telomere extension and maintenance were reinforced, and multiple DNA repair functions, such as Mismatch Repair, Base Excision Repair, Base Excision Repair, Nucleotide Excision Repair, and Homology Directed Repair appeared, and various metabolic functions, such as Adipogenesis, Metabolism Of Nucleotides, Onecarbon Metabolism, and Folate Biosynthesis were enhanced, and the above functions were further strengthened in the TMZ group

MAZ transcription factors are important regulators that promote the development of drug resistance in NPC-like clusters

The ‘pySCENIC’ package was employed for transcription factor (TF) analysis, and the results showed that there were distinct transcriptional activity profiles in the various tumor clusters (Fig. 3A). Specifically, a comparison of the TF activities of NPC-like clusters before and after TMZ treatment revealed that several transcription factors, including JUND, MZA, MXD4, and KLF9, were enhanced after treatment (Fig. 3B). Among the transcription factors, MAZ was identified as a pivotal regulator in the development of drug resistance in NPC-like clusters. In the control group, MAZ transcription factors exhibited high activity in a small number of NPC-like cells, whereas in the TMZ group, MAZ transcription factors showed high activity in the majority of NPC-like cells (Fig. 3C). Based on MAZ activity, NPC-like cluster was categorized into MAZ(+)_NPC-like and MAZ(-)_NPC-like cluster. The comparisons of functional characteristics showed that the MAZ(+)_NPC-like cluster exhibited enhanced DNA repair, stemness, metastasis, cell cycle, and adipogenesis (Fig. 3D). The ‘scTenifoldKnk’ package was used to perform virtual knockouts to predict the function of MAZ. Analysis results showed that 365 genes were found to be obviously altered (FC > 0.01) after MAZ virtual knockdown. These genes were significantly enriched in several key biological processes, such as cell division, DNA repair, DNA replication, DNA damage response, regulation of cell cycle, as well as several important pathways, such as Cell cycle, DNA replication, Mismatch repair, Nucleotide excision repair, and Platinum drug resistance in GO and KEGG databases (Fig. 3E).

Fig. 3

A Transcription factor (TF) analysis showed that there were distinct transcriptional activity profiles across the four tumor clusters. B Comparison of the TF activities of NPC-like clusters before and after TMZ treatment revealed that several transcription factors, including JUND, MZA, MXD4, and KLF9, were enhanced after treatment. C In the control group, MAZ transcription factors exhibited high activity in a small number of NPC-like cells, whereas in the TMZ group, MAZ transcription factors showed high activity in the majority of NPC-like cells. DThe MAZ(+)_NPC-like cluster exhibited enhanced DNA repair, stemness, metastasis, cell cycle, and adipogenesis. E MAZ gene virtual knockouts experiments identified 365 genes obviously altered, these genes were significantly enriched in several key biological processes, such as cell division, DNA repair, DNA replication, DNA damage response, regulation of cell cycle, as well as several important pathways, such as Cell cycle, DNA replication, Mismatch repair, Nucleotide excision repair, and Platinum drug resistance in GO and KEGG databases

Spatial transcriptome data confirm that MAZ(+)_NPC-like has high stemness and DNA repair characteristics

Spatial transcriptomic data for glioblastoma was obtained from the GSE235672 dataset[23], with the pathology slices containing samples before and after temozolomide treatment (Fig. 4A). MAZ gene expression levels between the before and after temozolomide-treated samples were compared, and the results showed that MAZ gene expression was significantly higher in the temozolomide-treated samples (Supplementary Material 1). Tissue stemness score and DNA repair score were also compared between the before and after temozolomide-treated samples, and results showed that tissue stemness score and DNA repair score were significantly higher in temozolomide-treated samples (Supplementary Material 2). The annotation of single-cell data was deconvolved into the spatial transcriptome data using the “RTCD” algorithm. The results showed that the content of MAZ(+)_NPC-like clusters was significantly increased in the TMZ treatment tissues (Fig. 4B). The enrichment scores of stemness and DNA repair gene sets were calculated for each spatial point. It was found that the content of MAZ(+)_NPC-like cluster showed weak associations with both stemness and DNA repair scores in untreated samples, while there was a high positive correlation in treated samples (Fig. 4C, D), suggesting that TMZ can stimulate MAZ(+)_NCP-like clusters to develop strong stemness and DNA repair capacity. Additionally, MAZ(+)_NPC-like cells were distributed diffusely across the tissue, suggesting that MAZ(+)_NPC-like can rebuild tumors multicentrally, which may be one of the reasons why glioma recurrence and progression are so rapid.

Fig. 4

A Glioma tissue slices contain specimens before and after temozolomide treatment (From the GSE235672 dataset in the GEO database). The left section is temozolomide untreated tissues, while the right section is temozolomide treated tissues. B The content of MAZ(+)_NPC-like clusters in each spot was calculated through “RCTD” deconvolution method. Differences in the content between tissues before and after TMZ treatment were compared. The results showed that the content of MAZ(+)_NPC-like clusters was higher in temozolomide-treated tissues. C In the untreated sample, there was a weak correlation between the Stem score and the content of MAZ(+)_NPC cluster for each spatial point. In the treated sample, there was a significant positive correlation between the Stem score and the content of MAZ(+)_NPC cluster for each spatial point. D In the untreated sample, there was a weak correlation between the DNA repair score and the content of MAZ(+)_NPC cluster for each spatial point. In the treated sample, there was a significant positive correlation between the DNA repair score and the content of MAZ(+)_NPC cluster for each spatial point

Relationship between the MAZ(+)_NPC-like cluster and clinical characteristics of glioma patients

The ‘BayesPrism’ package was adopted to deconvolve single-cell data into bulk RNA sequencing data of glioma patients in the TCGA and CGGA databases. Protein-coding RNAs were selected as features for inference (Fig. 5A). The cellular clusters composition of each patient was calculated by deconvolution and visualized by proportional charts (Fig. 5B). The distribution patterns of distinct cellular clusters across different pathological tissue types were analyzed, and the results are presented by heatmap (Fig. 5C). According to the heatmap, the samples with a high proportion of AC-like cluster predominantly corresponded to astrocytomas, while those with a high proportion of OPC-like cluster were mainly oligodendrogliomas (Fig. 5C). This pathological correspondence validates the accuracy of our single-cell annotations and the reliability of the ‘BayesPrism’ algorithms. Notably, both MAZ(+)_NPC-like cluster and MAZ(-)_NPC-like cluster were enriched in GBM, and immunohistochemical data from The Human Protein Atlas showed significantly higher expression levels of NPC-like cluster marker genes in high-grade gliomas (Fig. 5D, Supplemental Materials 3). These results suggest that the NPC-like clusters may be a characteristic subtype associated with GBM. For the characterization of MAZ(+)_NPC-like clusters, we found that that the prognosis of glioma patients with high proportion of MAZ(+)_NPC-like cluster is poorer (Fig. 5E). Moreover, the proportion of MAZ(+)_NPC-like clusters increased significantly in patients after TMZ treatment, while MAZ(-)_NPC-like clusters did not change significantly (Fig. 5F). Additionally, the proportion of MAZ(+)_NPC-like cluster increased with increasing WHO grade, and the proportion of MAZ(+)_NPC-like cluster was higher in older patients, IDH wild type patients, and recurrent gliomas patients (Fig. 5G). The scoring of DNA repair and tumor stemness revealed a significant positive correlation with the proportion of MAZ(+)_NPC-like cluster (Fig. 5H). The higher proportion of MAZ(+)_NPC-like cluster, stronger tumor stemness and DNA repair capacity. These findings further support that MAZ(+)_NPC-like clusters play an important role in tumor drug resistance.

Fig. 5

A The protein-coding RNA has a high correlation between single-cell data and bulk data. Thus, protein-coding RNA was selected as feature data for the deconvolution analysis. B Relative proportions of cell types (AC-like cluster, OPC-like cluster, NPC-like cluster, MES-like cluster, OC, and neuron) in glioma patients were assessed using BayesPrism algorithms. C Heatmap illustrating the relationship between clusters occupancy and pathological diagnosis. D The immunohistochemical images from The Human Protein Atlas showed significantly higher expression of NPC-like cluster marker genes (CENPF, GMNN, HMGB2, TOP2A) in high-grade gliomas. E Survival analysis revealed that glioma patients with high proportion of MAZ(+)_NPC-like cluster had poorer prognoses. F The proportion of MAZ(+)_NPC-like clusters increased significantly after TMZ treatment, whereas MAZ(-)_NPC-like clusters did not change significantly. G the proportion of MAZ(+)_NPC-like cluster increased with increasing WHO grade, and the proportion of MAZ(+)_NPC-like cluster was higher in older patients, IDH wild type patients, and recurrent gliomas patients. H The scores for DNA repair and tumor stemness were significantly positively correlation with the proportion of MAZ(+)_NPC-like cluster

MAZ transcription factors alter drug resistance of tumor cells by regulating the FoxM1 gene

In order to explore the downstream targets of MAZ, five databases (ENCODE, hTFtarget, Chip-Atlas, GTRD, and JASPAR) were employed to predict target genes of MAZ. The predicted results were then overlapped to obtain 1094 potential MAZ target genes. The correlation between MAZ and all genes in the single cell dataset was calculated, and 17 genes with high correlation with MAZ were screened out. The intersection of MAZ high-correlation genes and MAZ predicted target genes highlighted FOXM1 as the most promising downstream target of MAZ (Fig. 6A). The further calculation of the correlation of glioma samples in TCGA_GBM_LGG and CGGA_325 showed a significant positive correlation between MAZ and FOXM1 (Fig. 6B). These results suggest that MAZ is likely to affect the drug resistance of glioma by regulating the FOXM1 gene. In order to validate this result, cellular experiments were conducted. The results of PCR and WB experiment showed that when the MAZ gene was overexpressed, the mRNA and protein expression levels of the FOXM1 gene increased accordingly (Fig. 6C, D). Furthermore, CCK8 experiments showed that MAZ gene upregulation enhanced the survival of glioma cells in different concentrations of TMZ, with the IC50 increasing from 914.6 μM to 1626 μM. However, when the FOXM1 gene was further silenced, the phenomenon of MAZ gene promoting tumor cell resistance disappeared, and the IC50 decreased from 1626 μM to 1050 μM (Fig. 6E). Additionally, plate cloning assays showed that the clone formation ability of MAZ-overexpressing glioma cells against TMZ was significantly enhanced. However, when silencing FOXM1 in MAZ-overexpressing glioma cells, the cells'clone formation ability against TMZ was diminished (Fig. 6F). To enhance the specificity of the results, rescue experiments were conducted. The siRNA-resistant flag-FoxM1 plasmids (mutated in the siRNA target region) were used to restore FoxM1 protein expression levels (Fig. 6G). Restoration of FoxM1 expression successfully reversed si_FoxM1-induced alterations in TMZ sensitivity. Following the restoration of FoxM1 expression, the IC50 of glioma cells to TMZ increased from 924.3 μM to 1426 μM (Fig. 6H), and the clonal survival of glioma cells under TMZ treatment was significantly increased (Fig. 6I).

Fig. 6

A Integration of five databases (ENCODE, hTFtarget, Chip-Atlas, GTRD, and JASPAR) identified 1094 potential MAZ target genes. The intersection of MAZ high-correlation genes and MAZ predicted target genes highlighted FOXM1 as a key downstream target. B The relationship between MAZ and FOXM1 showed a significant positive correlation in glioma patients in TCGA_GBM_LGG and CGGA_325. C, D when the MAZ gene was overexpressed, the mRNA and protein expression levels of the FOXM1 gene increased accordingly. E TMZ toxicity assays showed that MAZ overexpression enhanced glioma cell survival in different concentrations of TMZ, with the IC50 increasing from 914.6 μM to 1626 μM. Silencing FOXM1 in MAZ overexpressing cells removed this phenomenon, with the IC50 decreased from 1626 μM to 1050 μM. F Plate cloning assays showed that MAZ overexpression enhanced the ability of glioma cells against TMZ. silencing FOXM1 in MAZ overexpressing glioma cells diminished this effect. G The flag-FoxM1 plasmids (mutated in the siRNA target region) were used to restore FoxM1 protein expression levels. H TMZ toxicity assays showed that, following the restoration of FoxM1 expression, glioma cell survival against TMZ enhanced significantly, with the IC50 increasing from 924.3 μM to 1426 μM. I Plate cloning assays showed that, following the restoration of FoxM1 expression, the clone capability of glioma cells under TMZ treatment was significantly increased

Drug combinations as a promising strategy

The oncoPredict package was used to evaluate the impact of MAZ(+)_NPC-like cluster on drug sensitivity. The results showed that the MAZ(+)_NPC-like cluster high occupancy group exhibited a higher drug score, indicating a higher IC50 for the drug (Fig. 7A). This result suggests that patients with high occupancy of MAZ(+) NPC-like clusters responded poorly to temozolomide, which is consistent with the previous results. Additionally, several drugs were identified as potentially effective for MAZ(+) NPC-like cluster high occupancy patients (Fig. 7B). However, single-cell level analysis revealed a different phenomenon. For example, dactinomycin and docetaxel were predicted to effectively kill glioma cells, but at the single-cell level, they were ineffective against the AC-like cluster and exhibited high neuronal toxicity. Dasatinib and trametinib were effective against AC-like, MES-like, and OPC-like significantly, with little neurotoxicity, but were weak on NPC-like (Fig. 7C). Considering the highly heterogeneity of gliomas, it is difficult to target all tumor clusters with a single drug. However, analysis of the combination of the two drugs revealed that the combination of Paclitaxel and Trametinib is a promising strategy. This combination can cover all cluster and have low neuronal toxicity (Fig. 7D). To verify this hypothesis, in vitro experiments were performed. We evaluated the survival of LN229 glioma cells at different drug concentrations of paclitaxel and trametinib (Fig. 7E). The results showed that monotherapy with paclitaxel (0.2 μM) or trametinib (0.1 μM) had limited cytotoxicity, but the combination of paclitaxel (0.2 μM) and trametinib (0.1 μM) produced significant toxic effects on glioma cells (Fig. 7F). Thus, the combination of paclitaxel and trametinib may be an effective way for glioma treatment.

Fig. 7

A the MAZ(+)_NPC-like high occupancy group exhibited higher drug scores, suggesting a higher IC50 for TMZ. B Drugs ranking by score revealed several potential effective drugs for MAZ(+) NPC-like high occupancy patients. C dactinomycin and docetaxel were predicted to effectively kill glioma cells, but at the single-cell level, they were ineffective against the AC-like cluster and exhibited high neuronal toxicity. Dasatinib and trametinib were effective against AC-like, MES-like, and OPC-like significantly, with little neurotoxicity, but were weak on NPC-like. D The combination of paclitaxel and Trametinib can cover all cluster and have low neuronal toxicity. E Drug toxicity assays showed that the IC50 for LN229 glioma cells was 0.866 μM in paclitaxel and 0.317 μM in trametinib. F Monotherapy with paclitaxel (0.2 μM) or trametinib (0.1 μM) showed limited cytotoxicity, while combination of paclitaxel (0.2 μM) and trametinib (0.1 μM) produced significant toxic effects on glioma cell viability

Discussion

Gliomas are highly heterogeneous tumors characterized by complex internal environments [24]. Nevertheless, Advances in single-cell technologies have progressively elucidated this complexity. Daan [25], Toshiro [26] et al. classified the tumor cells into four clusters (AC-like cluster, OPC-like cluster, NPC-like cluster, and MES-like cluster) based on the similarity of internal characteristics. Our study got the same results, we further identified the MAZ(+)_NPC-like cluster is the key cluster leading to the drug resistance in glioma. In the results of the functional enrichment analysis, MAZ(+)_NPC-like cluster shows highly stemness and activation of numerous DNA repair functions (such as MMR, BER, FA, and HR), which makes it capable to resist drug damage and regenerate tumor. However, considering the limited sample size in the single-cell and spatial transcriptomic datasets, we integrated the single-cell data into the TCGA and CGGA databases, which contain a large number of glioma patients. In the results, the samples with high proportions of AC-like are pathologically predominantly corresponding to astrocytoma, while the samples with high proportions of OC-like and OPC-like mainly corresponded to oligodendroglioma. This pathological matching validates both the accuracy of the single-cell annotation and the reliability of the BayesPrism algorithm. Further clinical correlation analysis revealed significant enrichment of the NPC-like cluster in GBM, suggesting that it may be a potential characteristic cluster of high-grade gliomas. Considering that this cluster exhibits high stemness and DNA repair capacities, which is also consistent with the biological behavior of high-grade gliomas. Notably, in the NPC-like cluster, the proportion of MAZ(+)_NPC-like cluster was closely associated with the prognosis of glioma patients treated with temozolomide. And compared to MAZ(-)_NPC-like clusters, MAZ(+)_NPC-like cluster had stronger stemness and DNA repair capacity. Combining the above results, we confirmed that the MAZ(+)_NPC-like cluster plays an important role, which is able to resist the efficacy of temozolomide during treatment and differentiate into various tumor clusters to rebuild the tumors after treatment. Furthermore, we also found some important functions in the MAZ(+)_NPC-like cluster. Such as, Lipid metabolism is enhanced in the MAZ(+)_NPC-like cluster, and several studies have demonstrated that lipid metabolism can provide energy for tumor cell proliferation, raw materials for cell membrane synthesis, maintain the stemness phenotype of tumor cells, and enhance drug resistance of tumor cells. In the study of Hao[27], researchers found that temozolomide combined with atorvastatin showed more significant antitumor activity than temozolomide alone. Thus, lipid metabolism may be an important mechanism for stemness characterization and drug resistance of MAZ(+)_NPC-like cluster, which needs more experiments to explore in the future.

In subsequent analyses, FoxM1 was identified to be a downstream target of MAZ. FoxM1 has been confirmed by a large number of studies to be closely related to multiple oncogenic pathways in glioma [28, 29]. FOXM1 has been shown to enhance the proliferation [30], migration [31], and drug resistance [32] of glioma cells. When upregulating MAZ in glioma cells, the expression of FoxM1 increased accordingly, and the tumors showed higher stemness and enhanced drug resistance. However, after silencing of FoxM1 in MAZ-overexpressing cells, the stemness and drug resistance of tumor cells were reduced. These results indicate that MAZ regulates the stemness and drug resistance of tumor cells by modulating FoxM1. However, the specific regulatory mechanism of MAZ on FoxM1 is not thoroughly elucidated, and in-depth exploration and research will be conducted in the future.

In the drug sensitivity analysis, we compared temozolomide sensitivity between high and low MAZ(+)_NPC-like proportion groups. As expected, patients with the high MAZ(+)_NPC-like proportion group had a higher IC50 for temozolomide, indicating that they were fewer sensitive to the drug, which was also consistent with the results of the single-cell and spatial transcriptome analysis. We further calculated all drugs from the GDSC2 database and found several candidates are sensitive to glioma patients. However, single-cell analysis revealed that these drugs were effective only against the few clusters. We attribute to the extremely heterogeneity of gliomas, where each cluster is characterized differently and therefore has different sensitivity to the drugs. To address this problem, we considered drug combinations as a feasible strategy. After exploration, we found that the combination of Paclitaxel and Trametinib can cover all clusters and has low toxicity to neurons. Paclitaxel and Trametinib have been shown in many studies to exert good anti-tumor effects[33, 34]. Thus, the combination of Paclitaxel and Trametinib is a promising drug combination for treating gliomas and overcoming drug resistance.

Conclusion

NPC-like cluster is prevalent in recurrent and drug-resistant gliomas.

MAZ transcription factors are important regulators that promote the development of drug resistance in NPC-like clusters.

The MAZ(+)_NPC-like cluster possesses a robust capacity for DNA repair and stemness, enabling resistance to TMZ-induced toxicity and facilitating tumour regeneration with enhanced drug resistance.

Patients with a high proportion of MAZ(+)_NPC-like clusters have a low sensitivity to temozolomide and a poor prognosis.

MAZ enhances stemness and drug resistance in glioma cells by upregulating FOXM1 expression.

The combination of paclitaxel and paclitaxel is a promising drug combination for treating gliomas and overcoming drug resistance.