The gene expression profile, tumor mutation burden (TMB), and clinical follow-up information data of the 160 ESCA cancer samples were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov). The validation cohort (n = 70) was obtained from the GEO database (GSE54993). The clinical characteristics of TCGA-ESCA and GSE54993 are presented in Table 1 and Supplementary Table S1 , and TCGA training and validation sets are shown in Supplementary Table S2 . The IMvigor210 cohort containing transcriptome data was downloaded from the website (http://research‐pub.gene.com/IMvigor210CoreBiologies). The immune cell-related datasets were derived from the NCBI Gene Expression Omnibus (GEO) data portal (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi), including GSE13906, GSE23371, GSE27291, GSE27838, GSE28490, GSE28726, GSE37750, GSE39889, GSE42058, GSE49910, GSE59237, GSE6863, and GSE8059. These datasets included gene expression data for 14 immune cells ( Supplementary Table S3 ). The data were processed using the “RMA” function in the R package “affy” (12), followed by the application of the function “removeBatchEffect” in the R package “limma” (13). The flow chart of this article is shown in Figure 1 .

The WGCNA (14) package in R software was applied to construct an immune-related gene co-expression network to identify significant gene modules. The clusterProfiler (15) package in R software was used to implement the Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway and gene ontology (GO) analysis. The ConsensusClusterPlus (16) package in R software was used for consensus cluster analysis to determine the number of subtypes in ESCA samples. The rationality of clustering was verified through the resampling method. We finally determined the optimal number of clusters by considering the cumulative distribution function (CDF) and delta region graphs. The ssGSEA method was performed to quantify the infiltration levels of immune cell types, functions, and pathways in the cancer samples, using the GSVA package (17) in R software. Additionally, the Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm was applied to explore the infiltration degrees of 22 immune cell types in esophageal carcinoma (18). The microenvironment cell populations counter (MCPcounter) analysis method was used to determine the level of tumor-infiltrating immune cells with the MCPcount package (19) in R software (20). The correlation analysis used the corr.test() function in R, when adjust=“bonferroni”, output the validated P value.

The tumor immune dysfunction and exclusion (TIDE) score was computed online (http://tide.dfci.harvard.edu/). The two main mechanisms of tumor immune evasion were the induction of T cell dysfunction in tumors with high cytotoxic T lymphocyte (CTL) infiltration, and the prevention of T cell infiltration in tumors with low CTL levels. Identify genes that affect cytotoxic T cell function on patient survival outcomes based on immune escape mechanisms. The z-score for each gene was the interaction coefficient d divided by its standard error. The accessible data from patients in GSE78220 treated with immunotherapies was used to predict the clinical response using the subclass mapping method. Additionally, the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) was used to predict the efficacy of the three subtypes with chemotherapeutic drugs. The univariate Cox regression analysis was used to screen for prognostic differential genes with a significance level of p<0.05, employing the “coxph” function of the “survival” package in R. Subsequently, a multivariate Cox regression analysis, utilizing stepwise regression in our study, was applied to further select significant genes using the “stepAIC” function of MASS package (21) in R. This process started with the most complex model and then successively deleted variables to reduce AIC based on AIC Akaike information guidelines. The Kaplan–Meier survival curve was used to analyze the survival difference between the low- and high-risk groups using the “survminer” R package. The predictive accuracy of the risk model was evaluated using the “timeROC” package in R. The “timeROC” package had a built-in function, when adjusted=TRUE, output the validated P value (22). TIMER2.0 (http://timer.cistrome.org/) was used to explore the correlation of CHMP7 with immune invasion in esophageal cancer.

Normal esophageal epithelial (HEEC cells) and ESCA (KYSE150 and TE-1 cells) cell lines were purchased from the cell bank of the Chinese Academy of Sciences (Shanghai, China). They were cultured in RPMI1640 medium containing 10% FBS and 1% penicillin-streptomycin, maintained in a humidified incubator containing 5% CO₂ at 37°C. For each experiment, cells in the logarithmic phase were used.

We used Lipofectamine2000 (Invitrogen, USA) to transfect CHMP7-siRNA (Si-CHMP7#1: GGAGGTGTATCGTCTGTAT; Si-CHMP7#2: CAAGGTCTCTCCAGTCAAT; Si-CHMP7#3: GAGTGAACAGCTTCTCTCA); pcDNA3.1-CHMP7(RiboBio, China) and NC-siRNA (RiboBio, China) into cells. The transfection efficiency was assessed using RT-qPCR 48 h after transfection.

The cells were collected, and the total RNA was extracted using Trizol (Invitrogen, USA). Subsequently, 1 µg of total RNA was reverse transcribed into cDNA using the HiScript II Q RT SuperMix for qPCR kit (Vazyme, China). RT-qPCR reaction was conducted on different samples using the Taq Pro Universal SYBR qPCR Master Mix kit (Vazyme) for a total of 40 cycles. The reaction conditions included a 95°C initial denaturation for 10 s, followed by 30 s at 10°C, using the QuantStudio5 real-time fluorescence quantitative PCR detection system. The sequences of all primers are listed in Table 2 .

When the cell fusion degree reached 90%, a straight line was drawn in the cell petri dish with the 1000 µl pipette. It was rinsed with PBS, and 1640 medium containing 2% FBS was added. Images under the microscope were observed and captured at 0 and 24 h.

The cells were collected, and 300 µl serum-free medium was added to the upper chamber with matrix glue. Medium containing 10% FBS was added to the lower chamber, with 50,000 cells in each well. After 24 h of culture, the medium was discarded, and the cells were fixed with 4% paraformaldehyde for 30 min. They were stained with 2% crystal violet for 10 min, followed by rinsing with PBS several times, drying, observing under a microscope, and capturing the images.

After collecting cells, 1000 cells were added to each well of four 96-well plates. At 0, 24, 48, and 72 h after cell attachment, 10 μl of CCK8 reagent was added. After 4 h of incubation, the OD value was measured using an enzyme labeling instrument. The daily OD values were compared with the OD value obtained at 0 h, and the curve was obtained to evaluate the cell proliferation rate.

The total cell protein was extracted, and the protein concentration was determined. The protein was separated using 10% SDS-PAGE and transferred to a PVDF membrane. After 2 h of incubation at room temperature, diluted antibody was added and allowed to incubate overnight at 4°C. Following the removal of the primary antibody and washing the next day, the corresponding secondary antibody was added and incubated at room temperature for 1 h. The protein band images were captured using a gel imaging system with ECL, and ImageJ was used for quantitative analysis. The protein indicators tested were BCL2 (Proteintech, China), BAX (Proteintech), and β-actin (Proteintech).

The statistical analysis was primarily conducted using the R programming language. The chi-squared test was used to compare differences between groups for categorical variables. For non-normally distributed variables, the Wilcoxon rank-sum test was conducted. ANOVA and Kruskal–Wallis test were selected for comparing more than two groups. The log-rank test was performed to evaluate the statistical differences in overall survival among different groups in the Kaplan–Meier survival analysis. The statistical significance was set at p<0.05 or p<0.01.

Given the expression of CD8+ T cells may influence the survival of esophageal carcinoma patients, we investigated the CD8+ T cell-related prognosis in esophageal carcinoma. Thirteen immune cell datasets were merged to form a single dataset with batch effects eliminated. Data before and after normalization were explored using principal component analysis (PCA). The transition from scattered datasets to a mixed state indicated the successful application of the “removeBatchEffect” function ( Figure 2A, B ). Following the normalization of the gene expression profile, hierarchical clustering was performed on samples based on the 179 expression profile data of the immune cell datasets ( Figure 2C ). The Pearson correlation coefficient was then used to calculate the distance between each gene. Using the aforementioned WGCNA method, we constructed a scale-free network with a determined soft threshold of 9 ( Figure 2D ), resulting in the identification of 13 gene modules ( Figure 2E ). The grey module comprised genes that could not be aggregated into other modules. We further analyzed the correlations between each module and immune cells ( Figure 2F ). The results indicated that the purple module, containing 346 genes, was the most positively correlated with CD8+ T cells and had little correlation with other immune cells. The genes of all the modules are listed in Supplementary Table S2 .

Subsequently, we performed functional enrichment analysis on CD8+ T cell-related genes belonging to the purple module. For GO analysis, 103 biological processes (BP) showed significant differences (FDR<0.05; Supplementary Table S3 ), with the top ten BP displayed in Supplementary Figure S1A . Eight cellular components (CC) demonstrated significant differential enrichment (FDR<0.05; Supplementary Figure S1B , Supplementary Table S3 ), and 12 molecular function (MF) signatures were significantly differentially enriched (FDR<0.05; Supplementary Table S3 ), with the top ten presented in Supplementary Figure S1C . In KEGG analysis, eight pathways were significantly differentially enriched (FDR<0.05; Supplementary Figure S1D , Supplementary Table S3 ). The findings suggested that these genes were closely related to immune functions and pathways.

We conducted univariate analysis separately on genes related to CD8+ T cells in TCGA-ESCA and GSE54993 datasets. The results revealed 16 and 41 genes related to prognosis in TCGA ( Supplementary Table S4 ) and GSE54993 ( Supplementary Table S5 ), respectively. Only two genes were common to both cohorts ( Figure 3A ), indicating limited consistency in CD8+ T cell-related genes among datasets of different platforms. Therefore, we selected 55 CD8+ T cell-related genes, identified as prognostic genes in two datasets, for further analysis (p<0.05).

We applied consensus clustering to categorize 160 patients from TCGA-ESCA cohort. The clustering results remained relatively stable with three clusters ( Figure 3B ). We selected k = 3 as the optimal number, leading to three CD8+ T cell-related immune clusters (ICs; Figure 3C , Supplementary Table S6 ). Analyzing the prognostic signature of these three subtypes revealed significant differences ( Figure 3D ), with IC2 displaying the poorest outcomes and IC3 demonstrating generally favorable outcomes. Similar results were observed in the GSE54993 cohort ( Figure 3E , Supplementary Table S7 ), suggesting the transplantability of the three molecular subtypes based on CD8+ T cell-related genes across different cohorts.

We compared the distribution of various clinical features (survival events, TNM stages, stage, age, and gender) among the three immune subtypes in TCGA-ESCA cohort to investigate the differences among them during analysis. The results indicated a significant difference in the proportion of T stages among the three subtypes, with the IC2 group having the highest proportion of T3, indicating a poor prognosis. No other significant clinical differences were observed among the three immune subtypes ( Supplementary Figure S2A-G ).

No significant difference was observed in TMB and the number of mutated genes among the three subgroups ( Supplementary Figure S3A, B ). Furthermore, we screened 1224 genes with a mutation frequency greater than three ( Supplementary Table S8 ). Then, 124 genes with high mutation frequency in each subtype were identified (p<0.05; Supplementary Table S9 ). The mutation signatures of the top 15 genes are depicted in the mutation heatmap in Figure 4A .

Previous research has suggested that chemokines play a pivotal role in tumorigenesis and tumor development (23). Chemokines can attract various immune cells into the tumor microenvironment, providing T cells access to the tumor and influencing both tumor immunity and therapeutic effects. In this study, we analyzed whether there were differences in the expression distribution of chemokines among immune subtypes based on TCGA-ESCA cohort. Thirty-eight out of 41 chemokines significantly differed in expression among subtypes ( Figure 4B ), indicating variations in the level of immune cell infiltration among different immune subtypes. These differences could contribute to distinctions in tumor progression and immunotherapy efficacy. Additionally, the expression of chemokine receptor genes was compared, and 15 out of 18 (83.88%) significantly differed among immune subtypes ( Figure 4C ).

CD8+ T cells in the tumor microenvironment can produce interferon-γ (IFN-γ), which stimulates the upregulation of PD-1/PD-L1 and IDO1 (24). We obtained the Th1/IFN-γ gene signatures from previous research (25) and calculated the IFN-γ score of each patient using the ssGSEA method. We found significant differences among immune subtypes in IFN-γ scores. IC3 subtype exhibited a relatively high IFN-γ score, whereas the opposite trend was observed for IC1 and IC2 groups ( Figure 4D ). Moreover, intra-tumoral immune T cell lysis activity was assessed using the average of GZMA and PRF1 expression values, according to a previous study (26). Significant differences were observed among the three subgroups ( Figure 4E ). IC1 and IC2 exhibited relatively lower immune T cell lysis activity, whereas IC3 displayed the opposite trend.

Similarly, the angiogenesis score of each patient was evaluated based on an angiogenesis-related dataset derived from a previous study (27). The results suggested that different immune subtypes differed in angiogenesis score, with IC2 and IC3 exhibiting higher scores than IC1 ( Figure 4F ). Moreover, we obtained 47 immune checkpoint-related genes from prior research, and a significant difference was observed in 41 genes ( Figure 4G ) through further exploration. Immune checkpoint-related genes such as LAG3, CTLA4, PDCD1, PDCD1LG2, and IDO1 were highly expressed in IC3. These results indicated that different subgroups may exhibit varying responses to immunotherapy.

The CIBERSORT method was used to evaluate the infiltration scores of 22 immune cells in each TCGA-ESCA dataset sample. Overall, significant differences in immune signatures were observed among different subgroups. Additionally, CD8+ T cells, resting memory CD4+ T cells, and M0, M1, and M2 macrophages were highly expressed in each subtype, suggesting their potential role in ESCA ( Figure 5A, B ). By analyzing the differences in ten oncogenic pathways proposed in a previous study (28), six pathways were found to be significantly different between immune subtypes ( Figure 5C ). Subsequently, the immune infiltration analysis revealed that IC3 exhibited the highest immune microenvironment infiltration score, whereas IC1 had the lowest score ( Figure 5D ). Moreover, most of the immune checkpoint-related genes had the highest expression in IC3, possibly contributing to the favorable outcome of IC3.

To explore the relationship between our immune subtypes and the six pan-cancer immunotypes studied previously, we collected the data of molecular subtypes from previous research (29) for comparison. Our subgroups significantly differed from the six immunophenotype groups in the previous study ( Figure 5E, F ). The three subtypes we defined could serve as a supplement to the existing six subtypes.

Here, we applied the TIDE software to evaluate the potential clinical effects of immunotherapy. In TCGA dataset, the TIDE score of IC1 or IC3 was significantly higher than that of IC2 ( Figure 6A ), suggesting that immunotherapy had a greater impact on IC2 than on IC1 or IC3. The predictive T cell dysfunction scores of IC1 were relatively lower, whereas those of IC3 were higher ( Figure 6B ). However, for the predictive T cell rejection scores, the score of IC1 was significantly higher than that of IC3 ( Figure 6C ), potentially explaining why IC1 exhibited a poor prognosis while IC3 demonstrated a better outcome.

Subsequently, the subclass mapping method was applied to compare the similarity of the three defined subtypes with patients treated with immunotherapy from the available dataset, GSE78220. A lower p-value indicated a higher similarity. The IC3 subtype showed similarity to anti-PD-1 non-resistance (anti-PD-1 NR) in TCGA dataset ( Figure 6D ), indicating that the IC3 group patients were more likely to respond to anti-PD-1 agents.

We investigated the sensitivity of different subtypes to conventional chemotherapy drugs using the same approach. A lower IC50 value indicated a higher sensitivity. The results revealed that IC1 was more sensitive than other subtypes to these chemotherapeutic drugs ( Figure 6E-I ).

Ninety-six samples were eventually included in the training cohort, and 64 samples were included in the validation cohort, with the clinicopathological characteristics of patients listed in Table 3 . The grouped results were tested for rationality, and no significant difference existed between them (p>0.05). By applying univariate Cox regression on the training set, we identified ten genes with prognostic significance (p<0.05; Supplementary Table S10 ). Based on multivariate Cox regression, six genes were selected to construct a prognostic risk model to avoid overfitting ( Figure 7A ).

The following six genes were selected: CHMP7, DNAJB1, KLRB1, PIK3R1, RCAN3, and RNF157. The coefficients of these genes were -0.691, 0.550, 0.638, -0.314, -0.829, and 0.273, respectively. Ultimately, the complete risk score was calculated using the following formula: Risk score = (-0.691 * gene expression value of CHMP7) + (-0.550 * gene expression value of DNAJB1) + (-0.638 * gene expression value of KLRB1) + (-0.314 * gene expression value of PIK3R1) + (-0.829 * gene expression value of RCAN3) + (0.273 * gene expression value of RNF157).

We determined the risk score of each sample derived from TCGA training set based on their gene expression level and plotted the risk score distribution. The distribution showed that the risk scores and mortality of the sample in the low-risk cohort were lower than those in the high-risk cohort ( Figure 7B ). A time-dependent ROC analysis was conducted to assess the predictive power of the 6-gene-based model. The area under the curve (AUC) of 1-, 2-, and 3-year predictions were 0.83, 0.87, and 0.81, respectively, indicating the high prognostic diagnostic competence of the 6-gene signature ( Figure 7C ). Additionally, the Kaplan–Meier plots indicated that the overall survival probability in the low-risk group was significantly better than that in the high-risk group (p<0.0001; Figure 7D ).

Based on the risk score formula, similar results were obtained using the aforementioned validation methods on TCGA validation cohort ( Supplementary Figure S4A-C ). Furthermore, the conclusion drawn from the entire TCGA-ESCA cohort suggested these six genes were prognostic for ESCA and that the 6-gene-based model was effective ( Supplementary Figure S4D-F ). To further examine the extrapolation capability of the 6-gene-based model to external populations, we validated its predictive power using the independent validation cohort, GSE54993 ( Supplementary Figure S4G-I ).

When comparing the distribution of risk scores in TCGA-ESCA dataset among different groups based on clinical features (age, gender, TNM stage, and stage), the results indicated the following: (1) Significant differences were observed in terms of risk scores among the different N staging groups and immune subtypes (p<0.05). The risk score increased with a higher N stage. Particularly, the IC2 subtype, associated with the worst prognosis in various subtypes, exhibited the highest risk score, whereas the IC3 subtype, associated with the best prognosis, showed the opposite trend. (2) Additionally, no significant differences in risk scores were observed among the different groups based on the clinical characteristics ( Figure 8A-G ).

To identify the 6-gene signature as an independent prognostic factor for clinical features, both univariate and multivariate Cox regression were performed. The results from the univariable analysis indicated that the risk type was of prognostic significance, with a hazard ratio (95% CI) of 2.89 (1.71, 4.89) and p<0.0001 ( Figure 8H ). Similarly, the corresponding hazard ratio (95% CI) in the multivariable analysis was 2.36 (1.71, 4.89; p = 0.006; Figure 8I ), suggesting that the risk type was significantly related to survival. Overall, these results revealed the independent prognostic value of the 6-gene-based risk type for patients with ESCA.

Currently, constrained by the lack of effective biomarkers for predicting the clinical benefits of immunotherapy, the identification of new predictive biomarkers is essential for advancing precision immunotherapy. The immunotherapy dataset (Imvigor210) was used to explore whether the 6-gene model can predict the efficacy of immunotherapy. The Kaplan–Meier curve showed that in metastatic urothelial carcinoma (mUC) patients treated with immunotherapy, higher risk scores corresponded to worse survival rates ( Figure 9A ). The AUC curves indicated the AUC of our risk model was greater than that of previously published signatures, including TMB and immunogenic neoantigen (NEO; Figure 9B ). The high- and low-risk groups were assigned as previously mentioned. Significant differences between the two groups were observed in responders and non-responders to immunotherapy ( Figure 9C ). The MCPcounter analysis method was applied to calculate the immune-cell infiltrating level, followed by determining the correlation between the risk score and TMB, NEO, and immune cells. The results suggested a negative correlation between the risk score and immune cell score ( Figure 9D ).

Additionally, we compared the differences in the risk scores among different groups and found significant differences in risk scores among the effectiveness of immunotherapy groups. However, no significant difference in risk score was observed among other immune characteristics grouping, including immune cells, tumor cells, and immunophenotypic grouping ( Figure 9E-H ).

Based on the aforementioned results, we can conclude that the 6-gene prognostic risk model exhibits a strong ability to predict prognostic risk. Further analysis of the related genes in the 6-gene prognostic risk model revealed that PIK3R1, RCAN3, and CHMP7 were associated with poor prognosis in patients with ESCA. Subsequently, we conducted RT-qPCR analysis, revealing that only CHMP7 was overexpressed in both TE-1 and KYSE150 cells ( Figure 10A ). Given that upregulated genes are easier to manipulate than downregulated genes in biological and therapeutic systems, we focused on the function of CHMP7 in ESCA. After the transfection of Si-CHMP7 into TE-1 and KYSE150 cells to inhibit its expression ( Supplementary Figure S5A, B ), we observed a significant inhibition in terms of cell migration, invasion, and proliferation ( Figure 10B-F ). Western blot analysis indicated a downregulation of BCL2 and an upregulation of BAX ( Figure 10G ). In addition, we upregulated the expression of CHMP7 in esophageal cancer cells by transfecting overexpressing plasmids ( Supplementary Figure S5C, D ), and subsequently observed a significant increase in cell proliferation ( Figure 10H ). At the same time, in order to further explore the relevant mechanism of CHMP7 regulating tumor growth, we also obtained the signaling pathways that interact with CHMP7 through online analysis software (https://www.genecards.org/) ( Table 4 ). From these findings, we infer that CHMP7 possesses the ability to influence the phenotype and apoptosis of ESCA cells.

Further investigating the role of CHMP7 in immune invasion of esophageal cancer, we ultimately found a positive correlation with Treg cells (0.393;p=4.97e-08; Figure 11A ) and CD8+ T cells (0.147; p=4.97e-2; Figure 11B ), when CHMP7 was analyzed for correlation with different immune cells. Subsequent analysis showed that CHMP7 was negatively correlated with naive CD8+ Tcells (-0.124; p=4.01e-03; Figure 11C ) and positively correlated with central memory CD8+ Tcells (0.154; P=3.89e-02; Figure 11D ). In addition, CHMP7 was not significantly correlated with effector memory CD8+ Tcells (-0.003; P=9.73e-01). It indicating that CHMP7 may interact with immune cells to affect the occurrence and development of esophageal cancer.