Identification of disulfidptosis in esophageal squamous cell carcinoma based on single-cell and bulk RNA-seq data to predict prognosis and treatment response

Purpose: Our study aims to identify the molecular subtypes of genes associated with disulfidptosis in esophageal squamous cell carcinoma(ESCC), develop a prognostic model, and identify potential therapeutic targets.

Methods: Based on the GSE53625 expression profile data, we identified molecular subtypes with significant survival differences through consensus cluster analysis. Subsequently, univariate Cox, multivariate Cox, and LASSO-Cox regression analysis were used to establish risk stratification models. The transcriptome data of the TCGA-ESCC cohort and the GSE160269 single-cell sequencing dataset were integrated to verify the biological significance of the model, and further analyze the heterogeneity of the tumor immune microenvironment, explore the differences in the intercellular communication network, and screen potential targeted drugs, providing a theoretical basis for subsequent translational research.

Results: We identified two distinct patterns of disulfidptosis expression with significant differences in overall survival. Then, we constructed the prognostic signature of disulfidptosis, and results showed patients with high score had worse prognosis. Univariate and multivariate Cox analysis demonstrated that the constructed prognostic signature was an independent prognostic factor and was validated in an independent validation set. The two subgroups differed in the proportion of immune cell infiltration and related signaling pathways in ESCC. The exploration of immunotherapy data confirmed our prognostic signature also had certain predictive power for immunotherapy. Drug screening suggested AZD8186 and JQ1 as potential therapies for high-score patients.

Conclusion: This study provides a new prognostic signature for ESCC, explores new therapeutic targets, and provides new theoretical support for personalized treatment.

Introduction

Esophageal cancer (EC) is one of the world’s deadliest cancers, consisting primarily of esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (1). Among them, ESCC is a highly lethal cancer in the esophagus, primarily concentrated in Asia and Africa (2, 3). Currently, the primary clinical treatments for ESCC are surgical resection, radiotherapy, and chemotherapy. However, most patients are already in the middle to late stage when detected because the early symptoms of ESCC are not obvious, leading to unsatisfactory treatment outcomes (4). Additionally, one of the core features of cancer is the ability to escape cell death, posing a great challenge for oncology treatment (5). Therefore, developing reliable and effective biomarkers and guiding individualized and optimal treatment modalities for ESCC patients is important.

Previous studies discovered that cancer cells can autonomously alter the flux of various metabolic pathways to increase bioenergy and biosynthetic demands and to mitigate the oxidative stress required for cancer cell proliferation and survival (6). When the metabolite accumulation exceeds the metabolic load of tumor cells, excess oxidative stress is generated, resulting in cell death (7). A recent study has identified a previously unexplained cell death caused by the rapid accumulation of excess cystine intracellular due to disulfide stress (8). This study discovered that in glucose-deprived SLC7A11-high cancer cells, intracellular NADPH is rapidly depleted (8, 9). The massive accumulation of disulfide molecules leads to abnormal disulfide bond cross-linking between actin cytoskeletal proteins and cytoskeletal contraction, disrupting their organization and resulting in actin network collapse and cell death, possibly leading to a new form of cell death, namely disulfidptosis (8, 10). Additionally, glucose transporter protein inhibitors could effectively inhibit cellular glucose uptake, causing NADPH depletion, actin cytoskeleton cross-linking, and disulfidptosis in SLC7A11-high cancer cells (10, 11). In vivo experiments on mice also demonstrated that GLUT inhibitors significantly inhibited tumor growth and induced aberrant disulfide cross-linking of actin cytoskeletal proteins in SLC7A11-high (10). This finding promises to be a new area of tumor therapy, but further research and exploration are needed to understand its specific mechanisms and therapeutic applications.

Currently, the study of disulfidptosis is still in its infancy, but relevant data are available to present its huge relationship with tumor development, such as lung and bladder cancers (12, 13). Similarly, we performed an analysis based on ESCC to explore the disulfidptosis-related genes expressions in specific carcinomas, resolve possible molecular linkages, and investigate potential prognostic targets. Therefore, this study aims to identify the molecular isoforms of disulfidptosis-related genes in ESCC in public databases, resolve the regulation pattern of disulfidptosis in ESCC, and construct a scoring model for disulfidptosis-related genes. Single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq data were combined to explore the different tumor growth patterns, clinical outcomes, immune microenvironment, and cellular communication and predict new potential compounds to provide new theoretical support for clinical decision-making.

Materials and methods

Data collection and preprocessing

We obtained the GSE53625 dataset (gene expression profiles and clinical data) from the Gene Expression Omnibus (GEO). After filtering, 179 tumor samples with complete expression and survival (overall survival [OS]) data were included as the training cohort. The R package TCGAbiolinks was used to download FPKM expression profiles (for log transformation), survival data, and clinical data from the TCGA-ESCA EC and kept 80 squamous cancer samples with both expression and survival data for validation. The expression profiles (UMI-count) of 60 squamous cancer samples from the GSE160269 single-cell dataset were selected for analysis in this project. A dataset of cutaneous melanomas treated with anti-CTLA-4 was downloaded and used to assess the predictive efficacy of signature in the immunotherapy cohort (14). Interaction relationships of disulfidptosis-related genes were obtained based on the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://www.string-db.org/), and Protein-Protein Interaction Networks (PPI) were constructed based on the interactions (10).

Consensus unsupervised clustering

Based on the expression profile data of 10 disulfidptosis-related genes, unsupervised cluster analysis was applied to identify different expression patterns of disulfidptosis. The ConsensusClusterPlus package was used for the operation, the distance used for clustering was Pearson, the clustering method was Pam, and 1000 replications were performed to ensure the stability of the classification. The R packages survminer and survival were used to generate survival curves for prognostic analysis using the Kaplan-Meier method. The log-rank test was used to determine the significance of differences and resolve the correlation between samples with different expression patterns and OS. Additionally, the R package limma was used to identify differentially expressed genes between different subgroups. We identified differentially expressed genes with |log2FC| ≥ 0.585 (fold change ≥1.5) and FDR < 0.05.

Constructing a prognostic signature

The hazard ratio (HR) and prognostic significance of differential genes and screened genes with p < 0.05 as prognostic correlates were determined using univariate Cox regression analysis.

We applied LASSO regression (glmnet R package) to select key prognostic genes. The risk score was calculated by summing the product of each gene’s expression level (exp) and its LASSO coefficient (coef): Score = Σexp*coef. The samples were divided into high and low groups according to the score. We generated Kaplan-Meier survival curves to compare patient outcomes. The significance of the differences was determined using the log-rank test to further resolve the correlation between these two types of samples and OS. The predictions of the scoring system on scoring were evaluated using the receiver operating characteristic curve (ROC), and the area under the curve (AUC) was visualized using the R package timeROC. Moreover, univariate and multivariate Cox analysis were performed to explore the independent prognostic value of the score.

Gene set variation analysis and functional enrichment

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses using the clusterProfiler R package (parameters pvalueCutoff = 0.05, pAdjustMethod = “BH”), and the R package GSVA was used to perform genes differentially expressed in the two groups of samples were screened using the difference multiplicity |log2FC| ≥ 0.585 (difference multiplicity great ne set variation analysis to annotate the potential functions of key genes and predict their potential molecular functions. Gene sets were downloaded from the Molecular Signatures Database (MSigDB) for the HALLMARK, KEGG, and Gene Ontology Biological Processes (GOBP) sublibraries and used to conduct GSVA analysis.

Tumor immune microenvironment assessment

The immune score, stromal score, and tumor purity were calculated for each tumor sample using the ESTIMATE algorithm. The Spearman correlation was calculated between score and immune score, stromal score, and tumor purity. According to previous studies, the single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to estimate the relative abundance of each cell infiltrate in TME. The distribution of immune cell infiltrates were compared between different grouped samples using the Wilcoxon test.

Predicting drug sensitivity

Based on the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) and Chemical and Systems Biology Program (CTRP) (https://portals.broadinstitute.org/ctrp/) cancer genomics drug susceptibility databases, R package oncoPredict’s calcPhenotype algorithm was used to evaluate the drug IC₅₀ values for each sample in the training set. The spearman correlation was calculated between score and drug IC₅₀ to assess the correlation between drug sensitivity and signature. The difference between drugs IC₅₀ in high and low groups was compared.

Single-cell transcriptome data quality control and identification of malignant cells

After the original authors’ quality control, 60 ESCC samples were analyzed from the GSE160269 dataset using the R package Seurat (v4.1.0) and normalized using the NormalizeData function. The FindVariableFeatures function was used to identify highly variable genes based on their average expression values (greater than 0.1 and less than 8) and dispersion (greater than 1) to identify highly variable genes. Batch correction between samples was performed using the R package Harmony to avoid batch effects interfering with downstream analysis. Then, the data were scale-transformed, principal component analysis was used for dimensionality reduction, the top 50 principal components were selected for downstream analysis, and the RunUMAP function was used for visualization.

The cell type annotation information provided by the original authors of the dataset was extracted. Among them, the malignant epithelial cells were identified using the copykat method. The CellScore was calculated using the AddModuleScore function of the Seurat package based on the signature-containing genes. Malignant epithelial cells were classified into high and low groups based on the median CellScore of malignant epithelial cells. Cell stemness scores were calculated using the AddModuleScore function and stemness gene sets from previously published studies.

Trajectory analysis and cellular communication

The R package monocle2 was used for trajectory analysis, time series analysis of malignant epithelial cells was proposed, and the R package CellChat was used for cell-to-cell communication analysis.

Statistical analysis

All analyses were performed using R version 4.1.2. For significance analysis between various values (expression, infiltration ratio, and various eigenvalues), the Wilcoxon rank sum test was used to compare differences between two groups of samples, whereas the Kruskal-Wallis was used to compare differences between multiple groups of samples. For plot presentation, ns indicates p > 0.05, * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001, and **** indicates p < 0.0001. Survival curves were generated for the prognostic analysis by the Kaplan-Meier method, and determined the significance of differences using the log-rank test.

Results

Consensus clustering identifies sample subgroups

We compared their expression differences between normal and tumor samples to assess whether the expression of disulfidptosis-related genes affects tumor progression in ESCC. The results revealed that six genes were significantly highly expressed in tumor samples, whereas two genes were significantly low (Figure 1A). NDUFS1 showed the highest connectivity in the protein interaction network (Figure 1B), suggesting its potential as a hub gene regulating disulfidptosis.

Figure 1

Figure 1. Consensus clustering of disulfidptosis-related genes in ESCC samples (A) expression of disulfidptosis-related genes in normal and tumor samples; (B) PPI of disulfidptosis-related genes; (C) consensus clustering of disulfidptosis-related genes, 1 and 2 represent two subgroups; (D) consensus clustering CDF plot; (E) consensus clustering cumulative distribution function Delta area; (F) OS curves of two subgroups; (G) clinicopathological characteristics and expression level heat plot of disulfidptosis-related genes in two subgroups. (ns, p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

We performed unsupervised cluster analysis using the R package ConsensusClusterPlus based on 10 disulfidptosis-related genes expression profiles of ESCC samples from the GSE53625 dataset. We identified two subgroups, named Cluster1 and Cluster2 (N = 85/94, Figures 1C–E). The Kaplan-Meier curves revealed that the two subgroups had significantly different prognoses, with Cluster1 having a worse OS (Figure 1F). We also compared the clinicopathological characteristics of different subgroups of ESCC using a heat plot. We discovered significant differences in the distribution of esophageal carcinoma sites (locations) between subgroups of patients (p < 0.05, Figure 1G).

Construction and validation of the prognostic signature

We constructed a disulfidptosis signature for predicting the prognosis of ESCC patients based on differential genes between the disulfidptosis expression patterns. Initially, we screened 493 differentially expressed genes between the two disulfidptosis expression patterns using the limma package (Figures 2A, B). GO and KEGG enrichment analysis of these differentially expressed genes using the clusterProfiler package. The KEGG results indicated that these genes were primarily enriched in pathways related to cytochrome P450, DNA adducts, reactive oxygen species (ROS), suggesting potential associations with drug toxicity, drug resistance, or metabolism of environmental carcinogens, involving mechanisms of carcinogenesis driven by DNA damage or oxidative stress. For the GO analysis: Biological Process (BP) terms were predominantly linked to xenobiotic substance metabolism and stress response, lipid metabolism and inflammatory regulation. Molecular Function (MF) terms were mainly associated with oxidation-reduction and detoxification functions, monooxygenase and lipid metabolism, and heme binding and oxygen metabolism. Cellular Component (CC) terms were enriched in processes related to cornified envelope formation, neural synapses, and regulation of cell polarity. (Figure 2C).

Figure 2

Figure 2. Molecular differences in disulfidptosis expression patterns (A) differential volcano plot; (B) differential gene expression heat plot; (C) functional enrichment analysis of differential genes.

Next, the univariate Cox regression analysis exposed that 47 of these differential genes were significantly associated with patient OS in the GSE53625 esophageal squamous cancer cohort (Supplementary Table 1), and the forest plot visualizes the Cox regression analysis results of the 20 genes with the smallest p-values (Figure 3A). Supplementary Figure 1 displays the KM curves of the six genes with the smallest to the greatest P-values.

Figure 3

Figure 3. signature construction and validation (A) forest plot of prognostic efficacy of top20 prognostic genes; (B) confidence interval of each Lambda of LASSO regression; (C) change trajectory of LASSO regression independent variables, the horizontal coordinate indicates the logarithm of the independent variable Lambda, and the vertical coordinate indicates the coefficient of the independent variable; (D) LASSO regression coefficient of key prognostic genes; (E) Prognostic survival curves for high- and low- score groups in the training set; (F) ROC curves in the training set; (G) Prognostic survival curves for high- and low- score groups in the validation set; (H) ROC curves in the validation set; (I) Independent prognosis for high and low groups of risk scores in the training set; (J) Independent prognosis for high and low groups of risk scores in the validation set.

Furthermore, LASSO-Cox regression analysis was performed based on these 47 genes. We performed 10-fold cross-validation under optimal conditions to determine the model’s penalty parameter (λ), eliminating 19 key prognostic factors affecting patient OS (Figures 3B–D). Based on the key prognostic factor expression levels and the linear combination of the corresponding weights, we construct the signature that can assess each patient’s prognosis. Table 1 illustrates the coefficients of each factor.

Table 1

Table 1. Key factors and the corresponding coefficients.

We calculate the risk score for each patient in the training set based on the constructed prognostic signature and divide them into high- and low-risk groups based on the median value. KM curve analysis and the log-rank test indicated that high-risk group patients had significantly shorter OS (p-value < 0.05, Figures 3E). The AUC for the predicted outcome of the sample was 0.712, 0.726, and 0.756 at one, two, and three years respectively (Figures 3F), indicating that the Score can provide a good characterization of sample OS. Then, we explored the independence of prognostic signature in the training set. Univariate and multivariate Cox regression models were constructed based on prognostic signature and clinical characteristics, and the results revealed that prognostic signature was an independent prognostic factor (HR = 2.31, p-value < 0.05, Figure 3I).

We used TCGA-ESCC as an independent validation set to assess the reliability of the prognostic signature. Patients were divided into high- and low-risk groups according to the prognostic signature risk score. Patient OS was also significantly lower in the high-risk group than in the low-risk group (Figure 3G), with predicted outcome AUCs of 0.649, 0.803, and 0.625 for the validation set samples at one, two, and three years, respectively (Figure 3H). Univariate and multivariate Cox regression models were constructed based on prognostic signature and clinical characteristics, and the results were consistent with the test set, again supporting prognostic signature as an independent prognostic factor (HR = 3.22, p-value < 0.05, Figure 3J).

Single-cell transcriptome analysis of prognostic signature

The original authors of the dataset quality-controlled scRNA-seq data of 192,078 cells (containing 20,335 B cells, 10,346 Endothelial, 44,547 Epithelial, 27,881 Fibroblast, 1,136 Fibroblastic reticular cells (FRC), 18,514 Myeloid, 3,023 Pericytes, and 66,296 T cell). Moreover, 17,986 genes were detected. The PCA results presented a significant batch effect between samples (Figure 4A), and the batch effect between samples was removed using Harmony (Figure 4B). UMAP demonstrated the distribution of different cell types (Figure 4C), and the proportion of cell distribution in each sample was heterogeneous (Figure 4D).

Figure 4

Figure 4. CellScore of single cell sequencing data and trajectory analysis (A) PCA analysis; (B) cell distribution before and after batch effect; (C) cell type distribution; (D) cell proportional distribution of samples. (E, G, I, K) trajectory distribution of State, Stemness, Pseudotime, CellScore; (F, J, M) UMAP plots of Stemness, CellScore, CellGroup of UMAP plot; (H, L) violin distribution of Stemness, CellScore in different States; (N) proportional distribution of CellGroup in different States.

We extracted epithelial cells identified by copykat as malignant and with at least one model gene detected (N = 19159, such as at least one model gene with expression greater than 0), calculated Stemness and CellScore for each cell, and classified CellGroup into high and low groups according to CellScore, using UMAP for presentation (Figures 4F, J, M). Trajectory analysis of malignant epithelial cells revealed three differentiation states (Figure 4E). State1 cells, characterized by low stemness (Figures 4G, H), represented the least aggressive phenotype. Tumor progression was associated with a gradual shift from State1 to State2/3, marked by increased stemness and loss of differentiation (Figure 4I). State1 → State2 and State1 → State3 trajectory routes have an increased CellScore (Figures 4K, L) and an increased proportion of high in the CellGroup (Figure 4N), indicating an increase in tumor malignancy.

Prognostic signature and immune microenvironment

Based on the bulk sequencing data, we calculated the pathway/biological process activities of KEGG and GOBP using the GSVA algorithm. We compared the differences between the activities of different groups using the rank sum test. The results demonstrated that the immune-related biological processes activities, such as T cell activation and T cell-mediated tumor cell immune response, were significantly higher in the high-score group than in the low-score group. High-score tumors exhibited stronger activation of pro-tumor pathways, including p53 signaling, chemokine production, and immune cell recruitment (Figure 5A).

Figure 5

Figure 5. Functional enrichment analysis of different groups (A) GOBP and KEGG pathway enrichment analysis of bulk sequencing data; (B,C) GOBP and KEGG pathway enrichment analysis of scRNA sequencing data; (D, E) GSEA analysis of bulk and scRNA sequencing data. (*, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

Based on scRNA-seq data, we calculated differentially characterized genes between high and low groups of CellGroup and performed GO and KEGG enrichment analysis using clusterProfiler. The results presented that genes were significantly enriched in immune-related biological processes, such as regulation of T-cell activation and proliferation, and tumor progression-related pathways, such as cell cycle, P53 signaling pathway, and antigen processing and presentation (Figures 5B, C), consistent with the results of bulk sequencing data.

Next, we analyzed the enrichment differences between high and low groups using GSEA. The results disclosed that the high-risk group was significantly enriched in immune-related pathways, such as B-cell activation and lymphocytes, in both bulk and scRNA-seq data (Figures 5D, E).

We calculated the percentage of immune cell infiltration in each tumor sample using the bulk sequencing data. The results revealed that the score of prognostic signature was significantly positively correlated with a stromal score, immune score, and ESTIMATE score, and significantly negatively correlated with tumor purity (Figure 6A). The percentage of infiltrated cells, such as Natural killer cells, Type1 T helper cells, and regulatory T cells (Tregs), was significantly lower in the Low group samples than in the High group (Figure 6B).

Figure 6

Figure 6. Immune landscape of different groups (A) correlation of immune stroma score and tumor purity with the score of prognostic signature; (B) distribution of the proportion of immune cell infiltration in different groups. (ns, p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

Differences in specific cellular communication between high-and low-prognostic signature score groups

Next, we performed intercellular communication analysis using the CellChat package, with extensive cellular communication among the cell populations (Figures 7A, B). When distinguishing between incoming and outgoing signaling, Malignant Epithelial Cells, Non-Malignant Epithelial Cells, Pericytes, Low neoplastic, Endothelial, and FRC are outgoing signaling dominant senders, while B cell, T cell, Fibroblast, High neoplastic, and Myeloid are incoming signaling dominant receivers (Figure 7C). Low and high neoplastic are in multiple tumor malignant progression-related pathways such as MIF, ITGB2, and FN1 with different cells for messaging (Figures 7D–F).

Figure 7

Figure 7. Cellular communication (A) all network diagram; (B) signaling dominant statistical heat plot; (C) signaling dominant statistical point diagram; (D, E, F) network of signaling pathways.

Potential treatment strategies for prognostic signature

We explored the predictive efficacy of prognostic signature on sample prognosis in the immunotherapy cohort (14). We discovered that patients with the same low score in the immunotherapy cohort had a better prognosis (Figure 8A). However, there was no significant difference in risk scores among patients in the immunotherapy with/without response group (Figure 8B), nor was there a significant difference in the proportion of samples responding to immunotherapy with/without immunotherapy in the different risk groups (Figure 8C).

Figure 8

Figure 8. Prognostic signature predicts drug sensitivity (A) KM survival curves in different score groups; (B) Score distribution in different response groups; (C) Distribution of immunotherapy response in different score groups; (D, E) Sensitivity analysis of Top 6 drugs with positive correlation with prognostic scores in the GDSC database; (F, G) Sensitivity analysis of Top 6 drugs with negative correlation with prognostic scores in the GDSC database; (H, I) Sensitivity analysis of Top 6 drugs with positive correlation with prognostic scores in the CTRP database; (J, K) Sensitivity analysis of Top 6 drugs with negative correlation with prognostic scores in the CTRP database. (HR, Hazard Ratio; CI, Confidence Interval; ***, p < 0.001; ****, p < 0.0001).

We predicted the IC₅₀ values of the drugs in the training set sample using the R package oncoPredict, GDSC, and CTRP database drug information combined with the expression profile of the training set. We compared the spearman correlation between the score and the IC₅₀ of each drug, ranked the drugs according to the absolute value of the correlation coefficient from the largest to the smallest, and selected the top six drugs with significant positive and negative correlations, respectively (Figures 8D, F, H, J, p < 0.05), with significant differences in drug IC₅₀ across Score groupings (Figures 8E, G, I, K).

Discussion

Recent studies have demonstrated that cells with high SLC7A11 expression exhibit increased susceptibility to cell death under glucose deprivation. This vulnerability primarily stems from reduced glycolysis under glucose-deficient conditions, leading to insufficient NADPH production. Without adequate NADPH, cystine cannot be reduced to cysteine, resulting in disulfide accumulation and subsequent disulfide stress. The accumulated disulfides interact with the actin cytoskeleton, where cross-linking between actin and cytoskeletal proteins disrupts actin architecture, ultimately triggering cell death (15, 16). Termed “disulfidptosis,” this process has been progressively recognized as a form of metabolic cell death, joining ferroptosis and cuprotosis as novel paradigms targeting tumor metabolic vulnerabilities (17–20). While these findings establish fundamental mechanisms of disulfidptosis in solid tumors, its regulatory networks and clinical translation potential in specific malignancies such as ESCC remain underexplored. In this study, we analyzed disulfidptosis-related gene expression differences in ESCC samples from public databases. Unsupervised clustering analysis revealed that higher risk scores were significantly associated with poorer prognosis. To investigate prognostic disparities between high- and low-risk patients, we performed differential gene expression analysis, followed by enrichment analysis, immune infiltration profiling, and drug sensitivithe mechanistic framework of disulfidptosisy evaluation of identified differentially expressed gsurvival curves for prognostic analysis were generated using the Kaplan-Meier methodenes.

ESCC is highly heterogeneous, leading to different clinical outcomes and treatment sensitivity (21). To address this issue, we enriched and analyzed the differential genes between the expression patterns of disulfidptosis and wanted to develop a signature to enable risk stratification and personalized treatment prediction. The findings suggest that tumor progression-related pathways such as DNA adducts, ROS, and lipid metabolism may contribute to prognostic differences in ESCC. The rapid proliferation and metabolic dysregulation of ESCC cells lead to excessive ROS generation, posing a significant threat to cellular survival. To counteract this damage, cells produce glutathione to scavenge excess ROS (22, 23). Furthermore, ROS can oxidize DNA to form adducts, activate DNA damage repair mechanisms, deplete NAD+ reserves, and exacerbate energy crises, indirectly promoting disulfide accumulation (24, 25). Cytochrome P450 enzyme-mediated metabolism of procarcinogens generates DNA adducts while potentially exacerbating reductive stress through NADPH consumption, thereby accelerating disulfide accumulation (26, 27). These signaling pathways are closely associated with the mechanism of disulfidptosis, also indirectly demonstrating that disulfidptosis is a form of cell death induced by disulfide accumulation resulting from cellular metabolic abnormalities.

Tumor microenvironment heterogeneity represents a critical factor in ESCC chemoradiotherapy resistance (18). Our study identified significantly enhanced activity of immune-related BP and tumor progression signaling pathways in the high-risk group, along with substantial differences in immune cell infiltration proportions, including natural killer cells, Th1 cells, and Tregs. Notably, IFN-γ secreted by CD8+ T cells was found to suppress the expression of disulfidptosis-related genes SLC3A2 and SLC7A11—two subunits of the cystine transport system (28, 29). This suggests that the anti-tumor effects of disulfidptosis may involve immunomodulatory mechanisms, potentially through CD8+ T cell functional inhibition or SLC3A2/SLC7A11-mediated cystine metabolism regulation impacting immune escape in ESCC. Importantly, analysis of immunotherapy cohorts confirmed enhanced treatment responsiveness in low-risk patients, indicating the prognostic model’s utility in identifying potential immunotherapy beneficiaries. However, the lack of significant correlation between risk scores and objective response rates to immune checkpoint inhibitors implies that this signature primarily reflects global TME characteristics (e.g., immunosuppressive status or metabolic stress levels) rather than directly influencing antigen presentation efficiency. This observation aligns with features of SLC7A11-high tumors, where previous studies have shown that SLC7A11 overexpression enhances tumor antioxidant capacity through GSH synthesis, fostering an immunosuppressive microenvironment (30, 31).

Targeted therapy, as a new therapeutic method, plays an important role in the treatment of ESCC (32), such as targeting the metabolic vulnerability of SLC7A11-high cancer cells, glucose transporter type 1 biosynthesis, and glutathione induces NADPH dissipation, significant disulfide molecules accumulation, such as cysteine, and ROS accumulation, thereby inhibiting tumor growth and spread (33, 34). Using the oncoPredict algorithm, we identified significantly correlated IC₅₀ values for therapeutic agents. Notably, AZD8186 and JQ1 exhibited lower IC₅₀ values and stronger cytotoxic effects in high-score tumors, suggesting synergistic interactions between their targeted pathways (PI3K, BET) and score-characterized metabolic vulnerabilities (35, 36). Conversely, high-score ESCC cells demonstrated significant resistance to EGFR inhibitors (e.g., Gefitinib) and BTK inhibitors (e.g., Ibrutinib), potentially mediated through compensatory EGFR pathway activation or drug efflux pump upregulation. While previous studies in bladder, lung, and breast cancers have validated associations between disulfidptosis-related models and chemotherapeutic sensitivity (12, 13, 37), this work establishes the first disulfidptosis-associated drug prediction framework for ESCC, providing a basis for personalized therapeutic strategies.

While mechanistic understanding of disulfideptosis continues to advance, the disease-specific regulatory networks in ESCC remain to be fully elucidated. Real-time dynamic monitoring of metabolic biomarkers including glutathione and ROS is crucial for addressing adaptive therapeutic resistance. Furthermore, the observed dissociation between our risk stratification model and immunotherapy outcomes implies the existence of additional regulatory factors modulating the disulfideptosis-immune crosstalk, which necessitates multi-omics integration for accurate predictive modeling. Finally, computationally identified drug candidates such as AZD8186 and JQ1 require rigorous experimental validation to confirm their therapeutic potential. Future investigations should incorporate systematic approaches combining in vivo and in vitro experiments to decipher the precise regulatory architecture of disulfideptosis in ESCC. Prospective clinical trials should be implemented to systematically evaluate the therapeutic efficacy of targeted strategies such as SLC7A11/GLUT1 inhibitors, while concurrently advancing the development of enhanced multi-omics-integrated prognostic frameworks.

Conclusions

This study provides a new prognostic signature based on disulfidptosis for ESCC, which can be used as an effective tool for predicting prognosis. We analyzed the prognostic indicators of ESCC and explored different tumor growth patterns, immune microenvironments, cell communication, and drug therapy, providing new theoretical support for clinical decision-making.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

XZ: Conceptualization, Data curation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. JD: Conceptualization, Methodology, Visualization, Writing – original draft, Writing – review & editing. XL: Investigation, Software, Visualization, Writing – original draft, Writing – review & editing. SZ: Data curation, Investigation, Methodology, Software, Visualization, Writing – review & editing. TZ: Conceptualization, Investigation, Methodology, Supervision, Writing – review & editing. MC: Conceptualization, Investigation, Resources, Writing – review & editing. GH: Methodology, Software, Supervision, Writing – review & editing. CC: Conceptualization, Investigation, Project administration, Resources, Supervision, Writing – review & editing. BZ: Conceptualization, Funding acquisition, Project administration, Resources, Software, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Acknowledgments

We thank Dr. Jerry for revising the English of this manuscript.

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.1567793/full#supplementary-material

Abbreviations

AUC, Area under the curve; BP, Biological Process; CC, Cellular Component; CTRP, Chemical and Systems Biology Program; CI, Confidence Interval; EC, Esophageal cancer; ESCC, Esophageal squamous cell carcinoma; FRC, Fibroblastic reticular cells; GEO, The Gene Expression Omnibus; GSVA, Gene set variation analysis; GO, Gene Ontology; GOBP, Gene Ontology Biological Processes; GDSC, Genomics of Drug Sensitivity in Cancer; HR, Hazard Ratio; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, Molecular Function; MSigDB, Molecular Signatures Database; OS, Overall survival; PPI, Protein-Protein Interaction Networks; ROC, Receiver operating characteristic curve; ROS, Reactive Oxygen Species; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins database; ssGSEA, Single-sample gene-set enrichment analysis; scRNA-seq, Single-cell RNA sequencing; TME, Tumor immune microenvironment; Tregs, regulatory T cells.

References

1. Smyth EC, Lagergren J, Fitzgerald RC, Lordick F, Shah MA, Lagergren P, et al. Oesophageal cancer. Nat Rev Dis Primers. (2017) 3:17048. doi: 10.1038/nrdp.2017.48

PubMed Abstract | Crossref Full Text | Google Scholar

2. Morgan E, Soerjomataram I, Rumgay H, Coleman HG, Thrift AP, Vignat J, et al. The global landscape of esophageal squamous cell carcinoma and esophageal adenocarcinoma incidence and mortality in 2020 and projections to 2040: new estimates from GLOBOCAN 2020. Gastroenterology. (2022) 163:649–658.e2. doi: 10.1053/j.gastro.2022.05.054

PubMed Abstract | Crossref Full Text | Google Scholar

3. Abnet CC, Arnold M, Wei WQ. Epidemiology of esophageal squamous cell carcinoma. Gastroenterology. (2018) 154:360–73. doi: 10.1053/j.gastro.2017.08.023

PubMed Abstract | Crossref Full Text | Google Scholar

4. Wang Y, Yang W, Wang Q, Zhou Y. Mechanisms of esophageal cancer metastasis and treatment progress. Front Immunol. (2023) 14:1206504. doi: 10.3389/fimmu.2023.1206504

PubMed Abstract | Crossref Full Text | Google Scholar

5. Pavlova NN, Zhu J, Thompson CB. The hallmarks of cancer metabolism: still emerging. Cell Metab. (2022) 34:355–77. doi: 10.1016/j.cmet.2022.01.007

PubMed Abstract | Crossref Full Text | Google Scholar

6. Lei G, Zhang Y, Koppula P, Liu X, Zhang J, Lin SH, et al. The role of ferroptosis in ionizing radiation-induced cell death and tumor suppression. Cell Res. (2020) 30:146–62. doi: 10.1038/s41422-019-0263-3

PubMed Abstract | Crossref Full Text | Google Scholar

7. Zheng P, Zhou C, Ding Y, Duan S. Disulfidptosis: a new target for metabolic cancer therapy. J Exp Clin Cancer Research: Cr. (2023) 42:103. doi: 10.1186/s13046-023-02675-4

PubMed Abstract | Crossref Full Text | Google Scholar

8. Liu X, Olszewski K, Zhang Y, Lim EW, Shi J, Zhang X, et al. Cystine transporter regulation of pentose phosphate pathway dependency and disulfide stress exposes a targetable metabolic vulnerability in cancer. Nat Cell Biol. (2020) 22:476–86. doi: 10.1038/s41556-020-0496-x

PubMed Abstract | Crossref Full Text | Google Scholar

9. Yan Y, Teng H, Hang Q, Kondiparthi L, Lei G, Horbath A, et al. SLC7A11 expression level dictates differential responses to oxidative stress in cancer cells. Nat Commun. (2023) 14:3673. doi: 10.1038/s41467-023-39401-9

PubMed Abstract | Crossref Full Text | Google Scholar

10. Liu X, Nie L, Zhang Y, Yan Y, Wang C, Colic M, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol. (2023) 25:404–14. doi: 10.1038/s41556-023-01091-2

PubMed Abstract | Crossref Full Text | Google Scholar

11. Joly JH, Delfarah A, Phung PS, Parrish S, Graham NA. A synthetic lethal drug combination mimics glucose deprivation-induced cancer cell death in the presence of glucose. J Biol Chem. (2020) 295:1350–65. doi: 10.1016/S0021-9258(17)49891-7

PubMed Abstract | Crossref Full Text | Google Scholar

12. Zhao S, Wang L, Ding W, Ye B, Cheng C, Shao J, et al. Crosstalk of disulfidptosis-related subtypes, establishment of a prognostic signature and immune infiltration characteristics in bladder cancer based on a machine learning survival framework. Front Endocrinol. (2023) 14:1180404. doi: 10.3389/fendo.2023.1180404

PubMed Abstract | Crossref Full Text | Google Scholar

13. Qi C, Ma J, Sun J, Wu X, Ding J. The role of molecular subtypes and immune infiltration characteristics based on disulfidptosis-associated genes in lung adenocarcinoma. Aging. (2023) 15:5075–95. doi: 10.18632/aging.204782

PubMed Abstract | Crossref Full Text | Google Scholar

14. AlDubayan SH, Conway JR, Camp SY, Witkowski L, Kofman E, Reardon B, et al. Detection of pathogenic variants with germline genetic testing using deep learning vs standard methods in patients with prostate cancer and melanoma. JAMA. (2020) 324:1957–69. doi: 10.1001/jama.2020.20457

PubMed Abstract | Crossref Full Text | Google Scholar

15. Chen Y, Zhang D, Yang H, Wu J, He W. Advances in the study of disulfidptosis in digestive tract tumors. Discover Oncol. (2025) 16:186. doi: 10.1007/s12672-025-01875-y

PubMed Abstract | Crossref Full Text | Google Scholar

16. Xie J, Deng X, Xie Y, Zhu H, Liu P, Deng W, et al. Multi-omics analysis of disulfidptosis regulators and therapeutic potential reveals glycogen synthase 1 as a disulfidptosis triggering target for triple-negative breast cancer. MedComm. (2024) 5:e502. doi: 10.1002/mco2.v5.3

PubMed Abstract | Crossref Full Text | Google Scholar

17. Mao C, Wang M, Zhuang L, Gan B. Metabolic cell death in cancer: ferroptosis, cuproptosis, disulfidptosis, and beyond. Protein Cell. (2024) 15:642–60. doi: 10.1093/procel/pwae003

PubMed Abstract | Crossref Full Text | Google Scholar

18. Huang TX, Fu L. The immune landscape of esophageal cancer. Cancer Commun (london England). (2019) 39:79. doi: 10.1186/s40880-019-0427-z

PubMed Abstract | Crossref Full Text | Google Scholar

19. Jiang X, Stockwell BR, Conrad M. Ferroptosis: mechanisms, biology and role in disease. Nat Rev Mol Cell Biol. (2021) 22:266–82. doi: 10.1038/s41580-020-00324-8

PubMed Abstract | Crossref Full Text | Google Scholar

20. Xiong C, Ling H, Hao Q, Zhou X. Cuproptosis: p53-regulated metabolic cell death? Cell Death differentiation. (2023) 30:876–84. doi: 10.1038/s41418-023-01125-0

PubMed Abstract | Crossref Full Text | Google Scholar

21. Fang P, Zhou J, Liang Z, Yang Y, Luan S, Xiao X, et al. Immunotherapy resistance in esophageal cancer: possible mechanisms and clinical implications. Front Immunol. (2022) 13:975986. doi: 10.3389/fimmu.2022.975986

PubMed Abstract | Crossref Full Text | Google Scholar

22. Li T, Song Y, Wei L, Song X, Duan R. Disulfidptosis: a novel cell death modality induced by actin cytoskeleton collapse and a promising target for cancer therapeutics. Cell Communication Signaling: CCS. (2024) 22:491. doi: 10.1186/s12964-024-01871-9

PubMed Abstract | Crossref Full Text | Google Scholar

23. Marini HR, Facchini BA, di Francia R, Freni J, Puzzolo D, Montella L, et al. Glutathione: lights and shadows in cancer patients. Biomedicines. (2023) 11:2226. doi: 10.3390/biomedicines11082226

PubMed Abstract | Crossref Full Text | Google Scholar

24. Zhao X, Ma Y, Li J, Sun X, Sun Y, Qu F, et al. The AEG-1-USP10-PARP1 axis confers radioresistance in esophageal squamous cell carcinoma via facilitating homologous recombination-dependent DNA damage repair. Cancer Lett. (2023) 577:216440. doi: 10.1016/j.canlet.2023.216440

PubMed Abstract | Crossref Full Text | Google Scholar

25. Mou Y, Wang J, Wu J, He D, Zhang C, Duan C, et al. Ferroptosis, a new form of cell death: opportunities and challenges in cancer. J Hematol Oncol. (2019) 12:34. doi: 10.1186/s13045-019-0720-y

PubMed Abstract | Crossref Full Text | Google Scholar

26. Huang J, Zhang J, Zhang F, Lu S, Guo S, Shi R, et al. Identification of a disulfidptosis-related genes signature for prognostic implication in lung adenocarcinoma. Comput Biol Med. (2023) 165:107402. doi: 10.1016/j.compbiomed.2023.107402

PubMed Abstract | Crossref Full Text | Google Scholar

27. Stipp MC, Acco A. Involvement of cytochrome P450 enzymes in inflammation and cancer: a review. Cancer Chemotherapy Pharmacol. (2021) 87:295–309. doi: 10.1007/s00280-020-04181-2

PubMed Abstract | Crossref Full Text | Google Scholar

28. Lee J, Roh JL. SLC7A11 as a gateway of metabolic perturbation and ferroptosis vulnerability in cancer. Antioxidants (basel Switzerland). (2022) 11:2444. doi: 10.3390/antiox11122444

PubMed Abstract | Crossref Full Text | Google Scholar

29. Wang W, Green M, Choi JE, Gijón M, Kennedy PD, Johnson JK, et al. CD8+ T cells regulate tumour ferroptosis during cancer immunotherapy. Nature. (2019) 569:270–4. doi: 10.1038/s41586-019-1170-y

PubMed Abstract | Crossref Full Text | Google Scholar

30. Wang T, Guo K, Zhang D, Wang H, Yin J, Cui H, et al. Disulfidptosis classification of hepatocellular carcinoma reveals correlation with clinical prognosis and immune profile. Int Immunopharmacol. (2023) 120:110368. doi: 10.1016/j.intimp.2023.110368

PubMed Abstract | Crossref Full Text | Google Scholar

31. Zhou S, Liu J, Wan A, Zhang Y, Qi X. Epigenetic regulation of diverse cell death modalities in cancer: a focus on pyroptosis, ferroptosis, cuproptosis, and disulfidptosis. J Hematol Oncol. (2024) 17:22. doi: 10.1186/s13045-024-01545-6

PubMed Abstract | Crossref Full Text | Google Scholar

32. Yang YM, Hong P, Xu WW, He QY, Li B. Advances in targeted therapy for esophageal cancer. Signal Transduction Targeted Ther. (2020) 5:229. doi: 10.1038/s41392-020-00323-3

PubMed Abstract | Crossref Full Text | Google Scholar

33. Liu X, Zhang Y, Zhuang L, Olszewski K, Gan B. NADPH debt drives redox bankruptcy: SLC7A11/xCT-mediated cystine uptake as a double-edged sword in cellular redox regulation. Genes Dis. (2021) 8:731–45. doi: 10.1016/j.gendis.2020.11.010

PubMed Abstract | Crossref Full Text | Google Scholar

34. Koppula P, Zhuang L, Gan B. Cystine transporter SLC7A11/xCT in cancer: ferroptosis, nutrient dependency, and cancer therapy. Protein Cell. (2021) 12:599–620. doi: 10.1007/s13238-020-00789-5

PubMed Abstract | Crossref Full Text | Google Scholar

35. Zhu W, Zhang Y, Yang L, Chen L, Chen C, Shi Q, et al. Construction of a lung adenocarcinoma prognostic model based on KEAP1/NRF2/HO−1 mutation−mediated upregulated genes and bioinformatic analysis. Oncol Lett. (2025) 29:155. doi: 10.3892/ol.2025.14902

PubMed Abstract | Crossref Full Text | Google Scholar

36. Ruiz de Porras V, Bernat-Peguera A, Alcon C, Laguia F, Fernández-Saorin M, Jiménez N, et al. Dual inhibition of MEK and PI3Kβ/δ-a potential therapeutic strategy in PTEN-wild-type docetaxel-resistant metastatic prostate cancer. Front Pharmacol. (2024) 15:1331648. doi: 10.3389/fphar.2024.1331648

PubMed Abstract | Crossref Full Text | Google Scholar

37. Liang J, Wang X, Yang J, Sun P, Sun J, Cheng S, et al. Identification of disulfidptosis-related subtypes, characterization of tumor microenvironment infiltration, and development of a prognosis model in breast cancer. Front Immunol. (2023) 14:1198826. doi: 10.3389/fimmu.2023.1198826

PubMed Abstract | Crossref Full Text | Google Scholar