Identification and validation of disulfidptosis-related clusters in glioblastoma

The status of seven programmed cell death (PCD) processes was evaluated in the pancancer dataset via ssGSEA based on well-established gene sets obtained from the literature and MSigDB database (Supplementary Table S1) (Liu et al. 2023; Zhang, et al. 2024). We found that the disulfidptosis score in GBM was the highest among all cancers, whereas the scores of other types of PCD were relatively low, suggesting disulfidptosis may be distinctive to GBM (Fig. 1A). Thus, we conducted NMF analysis based on the disulfidptosis related genes (DRGs) in the TCGA cohort (Supplementary Table S1). According to the consensus heatmap and silhouette value, GBM patients were finely classified into two clusters, which were referred to as C1 and C2 (Supplementary Fig. S1A). PCA confirmed that the disulfidptosis-related clusters were distinct from each other (Supplementary Fig. S1B). Nearest template prediction (NTP) analysis was applied to identifiy disulfidptosis clusters of GBM patients in the CGGA and Gravendeel datasets (Supplementary Fig. S1C). Survival analysis suggested that C1 patients had a worse prognosis than C2 patients did in all three datasets (Fig. 1B, Supplementary Fig. S1D, E). We then compared the clinical and molecular pathological features between the two clusters. In C1, higher proportions of classical and mesenchymal subtypes, and a greater tendency toward IDH mutation were observed (Fig. 1C and Supplementary Fig. S1F, G). We also noted that the proportions of secondary GBM and proneural subtypes were greater in C2 (Fig. 1C and Supplementary Fig. S1F). These findings indicate that the transcript levels of DRGs are correlated with the malignancy of GBM, which prompted the question of whether the survival of GBM patients could be predicted through DRG expression.

Fig. 1

Identification of disulfidptosis-related clusters. A Scaled averaged ssGSEA scores of seven PCD in pan-cancer. B Kaplan–Meier analysis (log-rank test) of the disulfidptosis-related clusters in TCGA cohort. C Clinical pathological features of the disulfidptosis-related clusters in pie charts (Fisher’s test)

Construction and validation of the disulfidptosis-related gene signature

We first utilized WGCNA to acquire gene modules associated with the disulfidptosis clusters. During WGCNA, we acquired 12 coexpression modules by setting the beta at 6 when the nonscale R2 was 0.8. The module–trait relationships revealed that the blue module was correlated with C1, and the red module was related to C2 (Fig. 2A). We also obtained differentially expressed genes (|log2FC|≥1 and P value < 0.05) from the two clusters (Fig. 2B). By intersecting cluster markers and module genes, we obtained characteristic genes, which were narrowed down to 30 on the basis of the betweenness scores calculated via protein-protein interaction (PPI) network analysis (Fig. 2C, D and Supplementary Table S2). Next, the Lasso-Cox regression was conducted with characteristic genes to construct a prognostic model (disulfidptosis-related gene signature, DRS) (Supplementary Fig. S2A and Supplementary Table S3). GBM patients were then separated into a high-DRS group and a low-DRS group for further analyses according to the median score. The Sankey plots showed that the majority of C1 had a relatively high DRS (Fig. 2E and Supplementary Fig. S2B).

Fig. 2

Construction of the disulfidptosis-related gene signature. A Identification of modular genes associated with disulfidptosis related clusters by WGCNA. B Volcano plot of differential expressed genes of the two clusters. C Schematic illustration of acquisition of the characteristic genes of the clusters. D PPI network of the top 30 characteristic genes ordered by their betweenness. E Sankey plot demonstrated the correspondence between the disulfidptosis-related clusters and scores. F Kaplan–Meier analysis of high-DRS and low-DRS groups in the TCGA cohort. G, H ROC curve and DCA curve depicting the performance of the DRS and other published gene signatures predicting 1-year OS. I Univariate and multivariate Cox regression of OS in the TCGA cohort. J The nomogram to predict 0.5-, 1.0-, and 1.5-year OS. K Calibration curve eliciting the performance of the nomogram in the TCGA cohort. L, M ROC curve and DCA curve depicting the performance of the DRS and other well-established biomarkers of GBM predicting 1-year OS

The high-DRS group had a better prognosis than the low-DRS group did (P < 0.001) (Fig. 2F and Supplementary Fig. S2C, H). The receiver operating characteristic (ROC) curve and the decision curve analysis (DCA) showed that DRS outperformed other published signatures in predictive performance (Fig. 2G, H and Supplementary Fig. S2D, E, I, J) (Zhang et al. 2023; Wang et al. 2023b; Sun and Liu 2024; Fan et al. 2022; Zhang et al. 2023). Additionally, the univariate and multivariate Cox regression analyses revealed that the DRS, alone with age, was an independent prognostic indicator (Figure 2I and Supplementary Fig. S2F, K).

However, we recognize that there is still room for improvement in the predictive power of DRS. Thus, we constructed a nomogram to further enhance its performance (Fig. 2J). Subsequent calibration analysis, ROC curve analysis, and DCA demonstrated a further improvement in the ability of the nomogram to predict 1-year survival compared to DRS alone, suggesting that integrating DRS with other clinical variables enhances predictive performance and provides a potential strategy for its clinical application. (Fig. 2K–M and Supplementary Fig. S2G, L).

Evaluation of the relationship between DRS and therapeutic response in GBM

Treatment resistance is believed to account for the poor prognosis of GBM patients (Johannessen et al. 2013). We compared the therapeutic response of the DRS groups in the Ducray cohort. The response rate to radiotherapy was lower in the high-DRS group (55.6% vs. 63.2%), while there was no difference in the chemotherapeutic response rate between the two groups (50.0% vs. 50.0%) (Fig. 3A). In addition, compared with that in untreated GBM cells, DRS was increased in radioresistant GBM cells but not in temozolomide-resistant GBM cells (Fig. 3B). For patients who underwent radiotherapy or even radiotherapy alone in all three cohorts, patients with high DRSs had shorter survival (Fig. 3C, D). However, for the patients who did not receive radiotherapy, no prognostic difference was observed between the two groups in any of the three cohorts (Fig. 3E). These findings indicated that the DRS was positively correlated with the radioresistance of GBM patients.

Fig. 3

The correlation between the DRS and radioresistance. A The therapeutic response rate of GBM patients with different DRS score in the Ducray cohort. B The DRS score among therapeutic response at cell level. C–E Kaplan-Meier survival analysis of DRS in subgroups in three cohorts (log-rank test). RT radiotherapy, RR radioresistant, CR chemoresistant

Single-cell analysis and validation in vitro identified EFEMP2 as the key gene of DRS associating with GBM cellular radioresistance

The causes of radioresistance can be attributed to in two factors: intrinsic resistance in tumor cells and microenvironment-related resistance (Xia et al. 2018; Yin et al. 2020). To reveal the mechanism of radiosensitivity heterogeneity in the DRS groups, we acquired a single cell dataset, GSE131928, and reannotated immune cells with canonical cell markers (Fig. 4A). Nonmalignant cells were annotated as macrophages (marked with CD163, F13A1, IL-10, MRC1, IL-1B, PTGS2, CCL4, and GPNMB), microglials (marked with TMEM119, and P2RY12), neutrophils (marked with S100A9 and S100A8), CD4+ T-cells (marked with CD3E, CD3D, CD4, and CTLA4), CD8+ T-cells (marked with CD3E, CD3D, CD8A, and CD8B) and oligodendrocyte (marked with OLIG1 and OLIG2) (Szulzewsky et al. 2015; Zeiner et al. 2019; You et al. 2023; Hu et al. 2022; Rajendran et al. 2023; Zhang et al. 2018). Higher proportions of nonmalignant cells were found in the high DRS group (Fig. 4B).

Fig. 4

The DRS was correlated with cellular radioresistance. A UMAP plot of tumor-infiltrating cells in GSE131928 (left) and dot plot showing the signature gene expression of cell types (right). B The proportion of cells that contributed to each cell type by the DRS groups. Color annotations match those of the UMAP plot. C Difference in ssGSEA scores of ‘HP_INCREASED_SENSITIVITY_TO_IONIZING_RADIATION’ pathway between malignant cells from each DRS group. D Pseudotime (left), the DRS group (middle) and radiosensitivity index (right) visualized on the cell trajectory plot of malignant cells. E Violin plots showing expression of EFEMP2 among cell types. F The transcriptional differences of EFEMP2 between GBM cells and normal brain cells. G The evaluation of radioresistance of EFEMP2-OE GBM cells by CCK8 assay. H The evaluation of radioresistance of EFEMP2-OE GBM cells utilizing clonogenic survival assay

The malignant cells from the high-DRS group had lower ssGSEA scores for radiosensitivity (Fig. 4C). Subsequently, we carried out cell trajectory analysis and found that the malignant cells from the high-DRS group were mainly located mainly at the beginning of the cell trajectory (Fig. 4D). The radiosensitivity index is a gene signature revealing the radioresistance of tumor cells (Eschrich et al. 2009). Consequently, we calculated the RSI for each malignant cell and discovered that high-RSI cells also resided at the start of the trajectory, which was in accordance with the notion of the low radiosensitivity of glioma stem cells (GSCs) (Fig. 4D). These results demonstrated that the DRS was negatively correlated with radiosensitivity. Additionally, the transcriptome levels of mesenchymal glioma stem cell (MES-GSC) markers (e.g., CD44 and YKL40) were greater in the high-DRS group than in the low-DRS group and decreased in pseudotime. Conversely, the expression levels of proneural GSC markers (e.g., CD133 and OLIG2) elevated in the low DRS group and increased by the pseudotime (Supplementary Fig. S3A).

Among the five DRS genes, EFEMP2 and TIMP1 were the prominent genes expressed in malignant cells (Fig. 4E and Supplementary Fig. S3B). However, qPCR showed that only EFEPM2 was highly expressed in all glioma cell lines, suggesting that EFEMP2 may be the key gene regulating the radiosensitivity of GBM cells (Fig. 4F; Supplementary Fig. S3C). On the other hand, EFEMP2 had the highest weight in the DRS (Supplementary Table S3). Based on these findings, we chose EFEMP2 for subsequent analysis. The overexpression of EFEMP2 significantly promoted the tolerance of GBM cells to ionizing radiation and their ability to proliferate (Fig. 4G, H; Supplementary Fig. S3D; Supplementary Fig. S4A, B). These results indicate that DRS reflects the cellular radiosensitivity of GBM.

Exploration of the correlation of DRS with the immunosuppressive microenvironment of GBM

GSEA was applied to investigate the mechanisms underlying the differences in prognosis and radiosensitivity between the DRS groups. Several immune related biological processes were enriched in the high DRS group (Fig. 5A; Supplementary Fig. S5A, B). The TIMER algorithm was used to evaluate the infiltration of immune cells in the GBM cohort, which revealed that the fractions of macrophages, dendritic cells, neutrophils, CD4+ T cells and CD8+ T cells increased in the high-DRS group (Fig. 5B; Supplementary Fig. S5C, D). Similar results were observed in the GSE131928 (Fig. 4B). These findings illustrated that the high DRS group might be characterized by a relatively active immune microenvironment and good prognosis, which contradicts our previous results. To elucidate the possible mechanisms of this paradoxical phenomenon, we conducted tracking tumor immunophenotype (TIP) analysis. TIP analysis revealed that the release and presentation of cancer cell antigens (Steps 1 and 2) and the trafficking and infiltration of immune cells (Steps 4 and 5) were significantly activated in the high DRS group (Fig. 5C; Supplementary Fig. S5E). Interestingly, highly activated early steps failed to enhance the cytotoxicity of T cells (Step 7) in the high DRS group (Fig. 5C; Supplementary Fig. S5E). Downregulated expression of cytotoxic makers was observed in CD8+ T-cells in the high DRS group (Fig. 5D). Immune escape is orchestrated by immune checkpoints, attenuating of cytotoxicity (Aktar et al. 2022). In the three GBM cohorts, high expression of immune checkpoints was detected in the high DRS group (Fig. 5E; Supplementary Fig. S5F). Further analysis revealed that increased expression of PD-1 on CD8+ T-cells and PD-L1 on malignant cells in the high DRS group (Fig. 5F). These results indicated that the antitumor immunity of the high DRS group was hindered by immune checkpoints at the execution stage.

Fig. 5

The DRS was associated with immunosuppression. A GSEA results of the DRS groups. B Differential analysis of infiltration of immune cell estimated by TIMER algorithm between the DRS groups. C Tracking tumor immunophenotype analysis of the DRS groups in the CGGA cohort. D Differential expression analysis showing genes related to T cell activation, cytotoxic and inhibition in CD8+ T cell. E Differential expression of immune check points in the TCGA cohort. F Differential expression analysis of genes associated with tumor immune escape in CD8+ T cells and malignant cells

Immunotherapy reversed the radioresistance phenotype of GBM cells overexpressing EFEMP2 in vivo

Given that the antitumor immunity of radioresistant GBM with a high DRS is diminished by immune checkpoints, immunotherapy might be an optimal treatment. By using the TIDE and submap algorithms, the immunotherapy response of patients from the TCGA and CGGA cohorts was estimated, suggesting that the high DRS group was sensitive to immune checkpoint blockade (P < 0.05) (Fig. 6A, C). On the other hand, there were more immunotherapy responders (PR/CR) in the high DRS group than in the low DRS group in the GSE91061 and GSE100797 cohorts (Fig. 6B). To determine whether high DRS GBM could benefit from immunotherapy, GL261/GL261-OE-EFEMP2 cells were implanted into C57BL/6J mice, which were then treated with an anti-PD-1 antibody alone, radiation alone, or their combination (n = 6) (Fig. 6D). The upregulated expression of EFEMP2 triggered radioresistance in the GL261-bearing mouse model (Fig. 6E, F). Additionally, high expression of EFEMP2 sensitized GBM to immune checkpoint blockade (Fig. 6E, F). Notably, the combination of radiation and an anti-PD-1 antibody exhibited a synergistic efficacy in mice overexpressing EFEMPE2, which was not observed in a mouse model without this combination (Fig. 6F). These findings revealed that increased expression of EFEMP2 was responsible for GBM resistance to radiation but sensitivity to ICB, which could inform the development of a synergistic radioimmunotherapeutic strategy.

Fig. 6

The association between DRS and immunotherapy response. A Differential proportions of the immunotherapy responders of the DRS groups in TCGA and CGGA cohorts predicted by the TIDE algorithm. B The Immunotherapy benefit rate of different DRS groups in GSE91061 and GSE100797 cohorts. C The potential response to immune checkpoint inhibitors of each DRS group estimated by submap algorithm. D Schematic illustration of combination irradiation with immune checkpoint blockade in EFEMP2OE and untreated C57BL6 mice. E Images of bioluminescence in mice during 21 days after intracranial injection of indicated GL261 cells. F Survival of C57BL6 mice bearing EFEMP2OE or untreated GL261 under different treatments (log-rank test, n = 6). PD progressive disease, SD stable disease, PR partial response, CR complete response, IR irradiation, ICB immune checkpoint blockade