Screening of genes related to programmed cell death in esophageal squamous cell carcinoma and construction of prognostic model based on transcriptome analysis

ABSTRACT

Objectives

To screen programmed cell death (PCD)-related genes in esophageal squamous cell carcinoma (ESCC) based on transcriptomic data and to explore its clinical value.

Methods

Differentially expressed PCD genes (DEPCDGs) were screened from ESCC transcriptome and clinical data in TCGA database. Univariate COX and LASSO COX were performed on prognostically DEPCDGs in ESCC to develop prognostic model. Differences in immune cell infiltration in different RiskScore groups were determined by ssGSEA and CIBERSORT. The role of RiskScore in immunotherapy response was explored using Tumor Immune Dysfunction and Exclusion (TIDE) and IMvigor210 cohorts.

Results

Fourteen DEPCDGs associated with prognosis were tapped in ESCC. These DEPCDGs form a RiskScore with good predictive performance for prognosis. RiskScore demonstrated excellent prediction accuracy in three data sets. The abundance of M2 macrophages and Tregs was higher in the high RiskScore group, and the abundance of M1 macrophages was higher in the low RiskScore group. The RiskScore also showed good immunotherapy sensitivity. RT-qPCR analysis showed that AUP1, BCAP31, DYRK2, TAF9 and UBQLN2 were higher expression in KYSE-150 cells. Knockdown BCAP31 inhibited migration and invasion.

Conclusion

A prognostic risk model can predict prognosis of ESCC and may be a useful biomarker for risk stratification and immunotherapy assessment.

KEYWORDS:

• Programmed cell death
• esophageal cancer
• prognosis
• immunotherapy
• immune cell infiltration

1. Introduction

Esophageal cancer is a common type of cancer, ranking eighth among all forms of cancer, and is the sixth leading cause of cancer-related deaths globally [1,2]. It is primarily made up of two main subtypes: esophageal squamous cell cancer (ESCC) and esophageal adenocarcinoma (EAC). ESCC is the most common type of esophageal cancer worldwide with a higher occurrence in East Asia and Africa. In contrast, EAC is more common in many developed nations [3]. Despite significant progress in managing ESCC through treatments such as surgery, endoscopic resection, chemoradiotherapy, and immunotherapy, this aggressive tumor poses a significant threat to patients due to its heterogeneity [4]. The development of molecular targeting and immunotherapy in recent years has improved the treatment of ESCC, but there are no reliable markers to predict the efficacy of treatment [5–7]. Therefore, the search for indicators that predict the efficacy of ESCC treatment and evaluate the prognosis has been a research hotspot.

Programmed cell death (PCD) is the spontaneous death of organismal cells in response to relevant signaling pathways and is essential for sustaining normal cell cycle and homeostasis in tissues. Malignant, stressed, damaged, or infected cells are lysed and effectively removed via PCD [8]. During tumorigenesis and progression, aberrant regulation of PCD often leads to uncontrolled proliferation and survival of tumor cells, making it an important aspect of tumor pathology. Earlier PCDs identified include apoptosis, autophagy, and necrosis [9]. In recent years, some novel PCD modalities have been gradually unearthed, including ferroptosis [10,11], cuproptosis [12,13]. Ferroptosis and cuproptosis are considered promising in immunotherapy and tumor radiation therapy [14]. It has been shown that tumor cell iron death is involved in the killing of tumor cells by CD8+ T cells. Activated CD8+ T cells secrete IFN-γ during treatment using immune checkpoint blockade with anti-PD-L1 antibody, lowering tumor cells’ intake of cystine by and increasing lipid peroxidation and ferroptosis [15]. Elemental copper can mediate cell death through mitochondria-dependent increased energy metabolism and reactive oxygen species accumulation-induced cytotoxicity [16]. Copper complexes are not only involved in pancreatic cancer cell death [17] but can also promote tumor growth and metastasis in breast cancer [18]. Cuproptosis was found to remodel the tumor microenvironment (TME) and to enhance the therapeutic response to anti-PD-L1 drugs [19]. It is very promising to explore novel therapeutic approaches for cancer species based on PCD mechanisms.

This study analyzed the prognostic features associated with PCD in ESCC based on multi-genomics and to construct a prognostic model. This study hopes to mine markers with clinical application value through systematic data analysis with the target of guiding the personalized treatment of ESCC patients.

2. Materials and methods

2.1. Data processing and analysis

Sequencing data (tumor: 162, normal: 11) of ESCC stored in TCGA source from the University of Cingifornia Sisha Cruz-Xena (UCSC-Xena) database (https://xenabrowser.net/) were downloaded. The expression values of each gene in the data are expressed as log2(FPKM +1). We also obtained survival information, clinical information (clinical stage, alcohol history, T stage, N stage, M stage, cigarettes per day, age, Barrett's esophagus) corresponding to each sample (Table S1). We also obtained single nucleotide variants (SNV) data for the tumor samples, containing the mutation site and mutation type for each gene within each sample.

2.1.1. Differential analysis and functional enrichment analysis

Differential analysis was performed using the limma package (version 3.32.3) [20] to identify differentially expressed genes (DEGs) in tumor samples in TCGA-ESCC, with normal samples as controls [21]. We extracted PCD genes from the Zou et al. [22] study and de-intersected with DEGs to identify differentially expressed PCD genes (DEPCDGs) in ESCC. To determine the function of DEPCDGs, up-regulated DEPCDGs and down-regulated DEPCDGs were uploaded to the DAVID database (https://david.ncifcrf.gov/) [23] for gene ontology enrichment analysis, respectively, and their involvement in the biological process. p value <0.05 for GO terms was considered to be statistically significant.

2.2. Prognostic model of PCD characteristics in ESCC

Univariate COX analysis of DEPCDGs in ESCC was performed to determine prognostic correlates (p < 0.05). Based on the function of glmnet [24], the effect of multicollinearity was removed by LASSO COX analysis to remove the DEPCDGs with high goodness-of-fit. Finally, according to the expression values of the screened DEPCDGs and the COX coefficients, we constructed the RiskScore for predicting the prognosis of ESCC. The formula for RiskScore was as follows.

$\mathit{RisksScore}=\sum \mathit{\beta i}\times \mathit{Expi}$

In TCGA-ESCC, the RiskScore of the tumor samples was evaluated. The samples were grouped into high RiskScore and low RiskScore groups by taking the median RiskScore value. Random sampling was divided into training and validation sets at 6:4 ratio. Kaplan–Meier (K-M) analysis was performed to show the survival region in TCGA-ESCC, training set, validation set. The survivalROC package [25] was used to plot receiver operating characteristic (ROC) curves for the three cohorts. The accuracy of RiskScore in predicting the survival of ESCC patients at 1, 3, and 5 years was assessed by area under the curve (AUC). Additionally, univariate COX and multivariate COX analyses were used to determine the prognostic independence of the RiskScore based on the RiskScore, clinical stage, alcohol history, T stage, N stage, M stage, cigarettes per day, age, and Barrett's esophagus of the samples in TCGA-ESCC. Based on the functionality of the rms package (version 3.4.0) [26], we integrated clinical information and RiskScore to draw nomograms predicting 1-, 3-, and 5-year survival in ESCC patients.

2.3. Correlation analysis between RiskScore and clinical information

To explore the relationship between ESCC clinical characteristics and RiskScore, we sequentially compared the distributions of clinical stage, alcohol history, T stage, N stage, M stage, cigarettes per day, age, and Barrett's esophagus among samples in the high/low RiskScore group.

2.4. Immune infiltration analysis

In the high/low RiskScore group, we calculated the difference in the degree of infiltration of each immune cell. From a previous study [27,28], marker genes of 28 immune cells were collected to perform ssGSEA. Meanwhile, by the CIBERSORT method [29], we assessed the infiltration abundance of 22 immune cells in the TME.

2.5. Immunotherapy sensitivity analysis

TIDE is an indicator that assesses tumor immune dysfunction and exclusion [30]. A higher TIDE score means a greater possibility of immune evasion, which indicates less immunotherapy benefit. We assessed the TIDE score of samples in TCGA-ESCC in the TIDE website (http://tide.dfci.harvard.edu/). Based on the functionality of the IMvigor210CoreBiologies package [31], we downloaded an immunotherapy treatment received by the sequencing data and treatment response information for the IMvigor210 cohort (Table S2). Based on treatment response, patients were defined as partial response (PR), progressive disease (PD), complete response (CR), and stable disease (SD). We also evaluated the survival trend of patients in the high/low risk group in IMvigor210.

2.6. Cell culture and transfection

Dulbecco’s Modified Eagle Medium (DMEM, Solarbio, Catalog no. 11995) was used to culture Het-1A and KYSE-150 cells (American Type Culture Collection, ATCC, Manassas, VA, USA) in a humidified atmosphere at 37°C with 5% CO₂. According to the manufacturer’s protocol of lipofectamine 2000 reagent (Invitrogen; Thermo Fisher Scientific), siBCAP31 (GeneChem Company) were used to knock down BCAP31 expression in KYSE-150 cells. siRNAs targeting BCAP31 (si-BCAP31, 5′- GUUGUUAGUUGGUGGCCGU-3′) and siRNA negative control (si-NC, 5′- UUCUCCGAACGUGUCACGU-3′) were synthesized by Shanghai GenePharma (China).

2.7. Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

TRIzol (Thermo Fisher, USA) reagent was used for the extraction of total RNA from quadriceps femoris of mice, then qRT-PCR was performed on a LightCycler 480 PCR System (Roche, USA) with the use of FastStart Universal SYBR Green Master. The reaction volume (20 μL) included an appropriate water volume,10 μL of PCR mixture, 0.5 μL of forward and reverse primers, and 2 μL of cDNA template. PCR reactions followed by these procedures: pre-degeneration at 95°C for 30 seconds (s), denaturation at 94°C for 15 s, annealing 56°C for 30 s, extension 72°C for 20 s under 45 cycles. Based on the 2^−ΔΔCT method, the expression level of genes was calculated and each sample was performed with three repetitions of technique. The primer pairs of target genes were listed in (Table 1), in which the GAPDH gene acted as a reference gene.

Table 1. Primers of genes.

Genes	Forward primer (5’-3’)	Reverse primer (5’-3’)
AUP1	AGACTGGCTGTGTAGACTTGAC	CTGTGGGTAGGGACTGAGTTC
BCAP31	CGGCTGGTGGAGTTGTTAGT	CGGGATTGTTCTGGAGGTT
DYRK2	GAAGCGCCAGGGCATGACAG	CGCAGGTGTTCCAGGATTCG
UBQLN2	GGAGGCGACATCAATGCAG	ACAGAAATGTGATTACGATGGCTGGGAG
TAF9	TCTTCCCTGAACGCTGGTTTC	CTTCCTTGTAGGATGCTGTGG
GAPDH	GGTGGTCTCCTCTGACTTCAACA	GTTGCTGTAGCCAAATTCGTTGT

2.8. Cell invasion

For the invasion assays, transwell plates were coated with Matrigel (BD Bioscience) and 5 × 10⁴ KYSE-150 cells were seeded onto the coated filter. Medium supplemented with 10% fetal bovine serum was added to the lower chamber. The cells remaining on the insert were removed after 24 h of incubation using cotton swabs. Next, 0.1% crystal violet was used to stain the cells, followed by photographing and cell number counting under a Qimaging Micropublisher 5.0 RTV microscope camera (Olympus).

2.9. Wound healing assay

KYSE-150 cells were seeded on a 6‐well plate. Cell debris was removed by washing the wounded monolayers using phosphate‐buffered solution. The morphological changes of the cells in each well were observed at 0 and 48 h under a microscope after the formation of wounds.

2.10. Statistical analysis

All the statistical analyses were performed in the R language (version 3.6.0) in the current study. Differences in continuous variables and categorical variables between the two groups were compared by the Wilcoxon rank-sum test and chi-square test, respectively. Log-rank test was used to compare differences in survival time between the two groups. Correlations between survival time and genes were assessed by univariate cox regression. In cellular experiments, we used Student’s t-tests to assess differences between treatment and control groups. Statistically significant difference was defined when p value <0.05.

3. Results

3.1. DEPCDGs in ESCC

From the tumor samples and normal samples in the TCGA-ESCC cohort, a total of 4497 DEGs (adjust p value <0.05 and |log2FC|>0.5) including3702 up-regulated DEGs and 795 down-regulated DEGs in the tumor samples were screened by performing of tumor samples and normal samples in the TCGA-ESCC cohort (Figure 1(a)). Taking its intersection with 1283 PCD genes, we extracted 340 DEPCDGs (Figure 1(b)). We found that the upwardly adjusted DEPCDGs are mainly concentrated in the apoptotic process, signal transduction, positive regulation of I-κB kinase/NF-κB signaling, cellular response to DNA damage stimulus, positive regulation of apoptotic process (Figure 1(c)). The downwardly adjusted DEPCDGs are mainly concentrated in the positive regulation of cell migration, positive regulation of autophagy, protein phosphorylation, positive regulation of macroautophagy, cellular response to glucose stimulus (Figure 1(d)). We also observed the mutational landscape of DEPCDGs in ESCC patients. We extracted the top 20 DEPCDGs in terms of mutation frequency and plotted a waterfall plot. We found that PRKDC, ACACA, and BRCA2 had mutation frequencies higher than 5% (Figure 1(e)).

Figure 1. Biological processes and mutational landscapes involved in DEPCDGs in ESCC. (a) Volcano plot of 4497 DEGs. Orange in the plot is up-regulated DEGs and green is down-regulated DEGs. (b) Wayne diagram of PCD genes and DEGs taking the intersection of 340 DEPCDGs obtained. (c) Upregulation of biological processes in which DEPCDGs are involved. (d) Down-regulation of biological processes involved in DEPCDGs. (e) Waterfall plot showing the mutational landscape of the top 20 DEPCDGs in terms of mutation frequency.

3.2. Prognostic modeling of PCD-related features in ESCC

Univariate COX analysis of 340 DEPCDGs was screened and 20 were associated with ESCC prognosis (p < 0.05). LASSO COX analysis identified 14 prognostically DEPCDGs (Figure 2(a,b)). BCAP31, DYRK2, UBQLN2, TAF9, AUP1, PPP1R15A, DEDD2, DDIT3, TBK1, PSMD10, HSPA1B, ZMPSTE24, CYCS, VAMP7 were screened (Figure 2(c)). In TCGA-ESCC, based on the median RiskScore value, the sample was clustered into a high RiskScore group and a low RiskScore group. The results of the K-M curves showed that the low RiskScore group exhibited a significant predominant survival profile (p < 0.001) (Figure 2(d)). Subsequently, according to the ratio of 6:4, the tumor samples in TCGA-ESCC were randomly sampled into training and validation sets. Tumor patients in the high RiskScore group exhibited significant inferior survival in both the two sets (p < 0.05) (Figure 2(e,f)). The ROC curves show that RiskScore has an impact on TCGA-ESCC (AUC at 1 year = 0.821, AUC at 3 year = 0.846, AUC at 5 year = 0.911) (Figure 2(g)), the training set (AUC at 1 year = 0.782, AUC at 3 year = 0.817, AUC at 5 year = 0.869) (Figure 2(h)), and the validation set (AUC at 1 year = 0.873, AUC at 3 year = 0.857, AUC at 5 year = 0.844) (Figure 2(i)) showed excellent robustness. In addition, to further validate the screened genes, we determined the mRNA expression levels of AUP1, BCAP31, DYRK2, TAF9, and UBQLN2 in KYSE-150 and normal cells, respectively, and found that all of these genes were significantly upregulated in expression (Figure S1(a–e)). Subsequently, BCAP31 as a key gene with high prognostic relevance to ESCC, therefore, we chose this gene for knockdown experiments and found that knockdown of BCAP31 inhibited the migration and invasion of KYSE-150 cells (Figure S1(f–i)).

Figure 2. Prognostic modeling of PCD-related features in ESCC. (a) Trends in the coefficient trajectories of the independent variables in the LASSO COX analysis. (b) Confidence interval for the penalty parameter λ based on 10-fold cross-validation. (c) LASSO COX regression coefficients for prognosis-related DEPCDGs. (d-f) K-M curves for patients in the high/low RiskScore group in TCGA-ESCC, training set, and validation set. (g-i) ROC curves for TCGA-ESCC, training set, and RiskScore in the validation set to predict 1-, 3-, and 5-year patient survival.

3.3. Correlation of RiskScore with clinical phenotypes

High-RiskScore patients tended to age older and smoked more cigarettes per day (Figure 3(a,b)). Compared with low RiskScore group, patients in the high RiskScore group had noticeably more long-term alcohol history and Barrett's esophagus pathology than in the (Figure 3(c,d)). Patients with advanced stages (T3+T4, N2+N3, M1, Stage III+Stage IV) were more concentrated in the high RiskScore group (Figure (e–h)).

Figure 3. Correlation of RiskScore with clinical phenotypes. (a) Age distribution of samples in the high/low RiskScore. (b) Distribution of the number of cigarettes smoked per day for samples in the high/low RiskScore. (c) Barretts esophagus pathology for samples in high/low RiskScore. (d) The alcohol history of the samples in the high/low RiskScore. (e-h) T stage, N stage, M stage, Clinical stage distribution of samples in high/low RiskScore.

3.4. Independence of nomogram and RiskScore

Whether RiskScore can serve as an independent indicator to assess the prognosis of ESCC patients was explored. RiskScore was confirmed as an independent indicator for ESCC prognosis (univariate COX: HR = 4.487, 95% CI = 3.058–6.585, p = 1.697e-14; multivariate COX: HR = 11.52, 95% CI=(3.098–42.8, p < 0.001)) (Figure 4(a,b)). Integrating RiskScore, clinical stage, alcohol history, T stage, N stage, M stage, cigarettes per day, age, and Barrett's esophagus, we plotted the predicted 1-, 3-, and 5-year survival of ESCC patients in the nomogram (C-index = 0.731) (Figure 4(c)). The difference between the overall survival time predicted by the nomogram at 1, 3, and 5 years, respectively, and the actual survival time was not significant (p>0.05) (Figure 4(d)).

Figure 4. Independence of Nomogram and RiskScore. (a,b) Univariate COX and multivariate COX analyses integrating RiskScore, T stage, N stage, M stage, clinical stage, alcohol history, cigarettes per day, age, and Barrett's esophagus. (c) Nomogram predicting 1-, 3-, and 5-year survival in ESCC patients. (d) Calibration curves for Nomogram prediction of 1-, 3-, and 5-year survival and actual survival in ESCC patients.

3.5. Relationship between immune cell infiltration microenvironment and RiskScore

According to the ssGSEA results of 28 immune cell activities, Immature dendritic cell, Type 2 T helper cell, CD56 dim natural killer cell within the high RiskScore group, effector memory CD4 T cell, CD56 bright natural killer cell, Type 17 T helper cell, Activated B cell had significantly lower scores than the low RiskScore group (p < 0.05) (Figure 5(a)). Based on the CIBERSORT results of the infiltration abundance of 22 immune cells, the infiltration abundance of Macrophages M2, T cells regulatory (Tregs) was noticeably higher in the high RiskScore group. On the contrary, compared with the low RiskScore group, high RiskScore group had significantly less infiltration (p < 0.05) abundance of Macrophages M1, Plasma cells, T cells CD8, and Dendritic cells activated (Figure 5(b)). The above data demonstrated that the immune microenvironment within the high RiskScore group was more suppressed.

Figure 5. Correlation between RiskScore and immune cell infiltration microenvironment. (a) ssGSEA results for 28 immune cell activity scores in the high/low RiskScore group. (b) CIBERSORT results for the abundance of 22 immune cell infiltrates in the high/low RiskScore group.

3.6. Relationship between immunotherapy and RiskScore

According to the results of TIDE score, patients in the low RiskScore group showed a higher TIDE score (Figure 6(a)). This suggests that ESCC patients with high RiskScore may be more suitable for immunotherapy. HAVCR2, LILRB4, and PDCD1 showed a trend toward higher expression in the high RiskScore group (Figure 6(b)). In the immunotherapy cohort IMvigor210, the proportion of patients with high RiskScore was higher in the CR/PR group (Figure 6(c)). Interestingly, in IMvigor210, patients with low RiskScore showed a significant survival advantage (p < 0.05) (Figure 6(d)). These results suggest that patients with a high RiskScore may be better suited for immunotherapy. RiskScore can also be used as a prognostic indicator for prognostic assessment inside the IMvigor210 dataset.

Figure 6. Correlation between RiskScore and immunotherapy. (a) TIDE score for the high/low RiskScore group. (b) Expression levels of immune checkpoint genes in the high/low RiskScore group. (c) Proportion of patients with high/low RiskScore in the CR/PR and SD/PD groups in the IMvigor210 cohort. (d) K-M survival curves for the high/low RiskScore group in the IMvigor210 cohort.

4. Discussion

ESCC is a common malignant tumor with a low survival rate, and chemotherapy as well as other methods fail to completely improve the prognosis of patients. Novel biomarkers are urgently required to enhance the survival of ESCC. In this study, DEPCDGs identified in ESCC were largely enriched with positive regulation of signal transduction, I-κB kinase/NF-κB signaling, apoptotic process, and positive regulation of apoptotic process, cellular response to DNA damage stimulus. DNA damage stimulus [32], NF-κB signaling [33] are important regulatory pathways in cancer progression.

In this study, 14 prognostically relevant DEPCDGs were included in the RiskScore, which were BCAP31, DYRK2, UBQLN2, TAF9, AUP1, PPP1R15A, DEDD2, DDIT3, TBK1, PSMD10, HSPA1B, ZMPSTE24, CYCS, and VAMP7. These genes are associated with cancer progression. Expression of BCAP31 is associated with cancer, and silencing BCAP31 expression inhibits emerin function thereby restoring β-catenin protein levels, a phenomenon that deregulates tumor suppression in colorectal cancer cells [34]. DYRK2 regulates the cell cycle of cancer cells and promotes tumor progression by modulating the G1/S transition. In addition, DYRK2 regulates epithelial mesenchymal transition and stemness in cancer cells [35]. UBQLN2 is associated with radiosensitivity in ESCC and is considered a reliable therapeutic target. It was found that silencing UBQLN2 in ESCC activated the p38 MAPK pathway and thus enhanced radiosensitivity [36]. TAF9 is also associated with radiosensitization, and TAF9 enhances the radiation tolerance of enhanced breast cancer cells [37]. AUP1 regulates lipid metabolic responses in renal clear cell carcinoma and induces lipid accumulation in cancer cells accelerating cancer progression [38]. A study by Wang et al. identified PPP1R15A as a promising therapeutic target in tumor immunotherapy by using single-cell data, whereby PPP1R15A counteracted the immune activity of antitumor drugs [39]. DDIT3 regulates gastric cancer cell stemness [40]. TBK1 regulates inflammatory cytokine production and immune system activity in the immune microenvironment by activating NF-κB signaling [41]. PSMD10 was found to promote tumor cell cycle progression [42]. Up-regulation of HSPA1B expression contributes to poor prognosis in colorectal cancer [43]. ZMPSTE24 is a DNA damage response-associated gene, which was found to exhibit concordance in colorectal and gastric cancers with the MSI-H phenotype [44]. CYCS is a pro-apoptotic factor, and CYCS inhibits cellular value-addition and promotes apoptotic responses in tumor cells to inhibit cancer progression [45]. VAMP7 is an important autophagy factor that regulates the immunogenicity of cancer cell death or inflammation [46]. These genes are associated with cancer progression or different PCD modalities, and ESCC prognostic models constructed with these genes are reasonable.

Among the RiskScore subgroups, the infiltration abundance of Macrophages M2, Tregs was high within the high RiskScore group. M2 type macrophages stimulates tumor growth development immunosuppression, promotion of tissue remodeling, metastasis, and neo-angiogenesis. Different stimulants, including IL-4, IL-13, IL-1α and other cytokines, as well as cellular metabolites such as spermine and spermidine, can all activate the polarization of M2 macrophages [47,48]. Tregs play a dual role in tumors, promoting tumor development or inhibiting tumor progression, which correlates with the tumor ecological niche. Tregs exhibit pro-tumorigenic effects when they limit the inflammatory response [49]. In contrast, the infiltration abundance of Macrophage M1 was higher in the low RiskScore group. Macrophage M1 showed anti-tumor activity by up-regulating the expression of pro-inflammatory molecules, such as inducible iNOS, tumor necrosis factor, and IL-1β, which in turn showed anti-tumor activity [50]. The results based on immune infiltration analysis suggest that a pro-tumorigenic environment was exhibited in the high RiskScore group, which could be a potential culprit for their poor prognosis. This phenomenon is also reflected in immunotherapy sensitivity. ESCC patients with high RiskScore may be more suitable for immunotherapy. And the expression levels of HAVCR2, LILRB4, and PDCD1 were higher in the high RiskScore group. The essence of the immune escape phenomenon is that the immune checkpoint genes on the surface of tumor cells inhibit the immune regulation of tumor cells, and the immune cells are restricted from immunotherapy effects [51]. Immunotherapy activates immune cell activity and relieves immune escape in several ways [52]. Based on this phenomenon, we hypothesized that excessive immune checkpoint proteins in the tumor microenvironment of patients in the high RiskScore group may have suppressed immune cell activity. It can be seen that the ESCC prognostic model may be able to assess the immune status of the tumor microenvironment of patients using RiskScore and provide more precise treatment guidance. However, this study still has some limitations. First, this study lacked cellular assays to validate the function of DEPCDGs in RiskScore in ESCC. Second, in order to be able to use the present model in the clinic, the accuracy should be further validated in a multicenter, large-sample cohort subsequently.

5. Conclusion

Overall, this study identifies a prognostic model for PCD characterization in ESCC through comprehensive and integrated data analysis and allows assessment of immune microenvironmental status and immunotherapeutic susceptibility in patients with different RiskScore subgroups. It demonstrated high robustness in the evaluation cohort and is a potentially reliable metric. Our study provides an important basis for future research on clinical treatment and immunotherapy in ESCC.

Abbreviations

ESCC	=	Esophageal squamous cancer
PCD	=	Programmed cell death
UCSC-Xena	=	University Of Cingifornia Sisha Cruz-Xena
TCGA	=	The cancer genome atlas
SNV	=	Single nucleotide variants
DEGs	=	Differential expression genes
DEPCDGs	=	Differential expression programmed cell death genes
LASSO	=	Least absolute shrinkage and selection operator
ROC	=	Receiver operating characteristic
AUC	=	Area under the curve
ssGSEA	=	Single-sample gene set enrichment analysis
TIDE	=	Tumor Immune Dysfunction and Exclusion
CR	=	Complete response
PR	=	Partial response
SD	=	Stable disease
PD	=	Progressive disease

Author contributions

All authors contributed to this present work: M Chen, YJ Qi, and HD Chen concepted and designed the research, SH Zhang, YB Du and SG Gao acquired the data, SG Gao, HD Chen analyzed and interpreted data. SG Gao, SH Zhang and M Chen drafted the manuscript, M Chen, HD Chen, YB Du and YJ Qi revised manuscript for important intellectual content. All authors read and approved the manuscript, and to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Declaration of interest

The authors have no relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript. This includes employment, consultancies, honoraria, stock ownership or options, expert testimony, grants or patents received or pending, or royalties.

Reviewer disclosures

Peer reviewers on this manuscript have no relevant financial or other relationships to disclose.

Data availability statement

The full analysis of raw data is available in the Github (https://github.com/1MinChen/Raw-data-for-article.git).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/14737140.2024.2377184.

Additional information

Funding

This manuscript was funded by National Natural Science Foundation of China “Porphyromonas gingivalis promotes epithelial mesenchymal transition in esophageal squamous cell carcinoma by activating AKT/β- Catenin signaling pathway” [No.81972571].