The transcriptome expression profiles of 81 ESCC tissues were downloaded from TCGA and the corresponding clinical information (https://portal.gdc.cancer.gov/). The 340 normal esophageal mucosal epithelium samples were obtained from GTEx cohorts (https://www.gtexportal.org/home/index.html). The transcriptome expression datasets of 19 esophageal squamous carcinoma cell lines (TE1, TE4, TE5, TE6, TE8, TE9, TE10, TE11, TE14, TE15, KYSE30, KYSE70, KYSE140, KYSE150, KYSE180, KYSE270, KYSE410, KYSE510, and KYSE520) were downloaded from CCLE (http://www.broadinstitute.org/ccle/home). Using R software to integrate TCGA and GTEx data, 274 normal esophageal mucosa samples and 81 ESCC samples were selected for further analysis. The immune-related gene list was downloaded from the ImmPort dataset (https://immport.org/home).

The Limma package in R was used to screen the differentially expressed genes in ESCC and normal esophageal mucosa samples. A fold change >2 and p-value < 0.05 were the screening conditions (Law et al., 2014).

The ClusterProfiler, org. Hs.eg.db, enrichplot, ggplot2, and GOplot packages in R were used to perform functional enrichment analyses of the differentially expressed genes. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis results are shown in circular maps (Bindea et al., 2009). The Gene Set Enrichment Analysis (GSEA) was performed to analyze the distribution trend in gene table sorted by phenotypic correlation.

The WGCNA package in R was used for the cluster analysis of immune-related differentially expressed genes and screening of prognostic differentially expressed genes. Kaplan-Meier (KM) survival curve analysis was used to further verify the accuracy of the impact of these genes on prognosis (Langfelder and Horvath, 2008).

CIBERSORT is a useful tool for evaluating the abundance of various immune cell types in the tumor microenvironment. A p-value < 0.05 was a selective condition, and 73 ESCC samples were obtained for further analysis (Newman et al., 2015).

The area under the ROC curve (AUC) was used to calculate the predictive ability of the ESCC prognosis model and the COX regression analysis was used to calculate the risk coefficients of ESCC patients. The effect was compared with conventional indexes (Kamarudin et al., 2017).

Human ESCC cell lines (TE13, EC109) and normal human esophageal epithelial cells (HEEC) were provided by cardiovascular thoracic surgery department of Tianjin Medical University General Hospital. Cells were cultured with RPMI-1640 containing 10% Fetal Bovine Serum (FBS). All cells were maintained in a humidified chamber containing 5% CO₂ at 37°C.

Total RNA of cultured cells was isolated with RNAiso Plus (TaKaRa) and reverse-transcribed to cDNA with PrimeScript Strand cDNA Synthesis Kit (TaKaRa). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed with TB Green Premix Ex Taq II (TaKaRa). The primer sequences are shown in Supplementary Table S1. Relative gene expression was determined by the comparative 2^−ΔΔCT method.

The flow chart of the study is shown in Figure 1. We downloaded the transcriptome expression profiles of 95 ESCC tissues and 340 normal esophageal mucosal epithelium samples from TCGA and GTEx cohorts, then performed batch integration analysis on the data. Finally, 274 normal esophageal mucosa samples and 81 ESCC samples with 34,350 genes were selected for further analysis (Supplementary Data S1). Differentially expressed genes were screened by Limma. Overall, 6,231 genes were found to be differentially expressed between ESCC and normal tissues. Of these, 3,008 genes were downregulated and 3,223 genes were upregulated in ESCC tissues relative to normal tissues.

Next, we downloaded the list of 6,196 immune-related genes from the ImmPort dataset and screened for any differentially expressed genes among them. There were 647 differentially expressed immune-related genes, with 418 upregulated and 229 downregulated in ESCC tissues relative to normal tissues. The differentially expressed genes are shown in Supplementary Table S2. The heatmaps of differentially expressed immune-related genes and other differentially expressed genes are shown in Figure 2A; Supplementary Figure S1, respectively.

Immune reactions and viral infections may be the main factors that maintain the malignant state of ESCC. GO and KEGG analyses were performed on the differentially expressed immune-related genes, with the results showing that the first eight most significant GO terms enriched were all related to immune response. Additionally, three of the first eight KEGG terms enriched were related to immune response and three were related to viral infection. The enriched KEGG terms, B cell receptor signaling pathway, cytokine-cytokine receptor interaction, and chemokine signaling pathway, were consistent with the enriched GO terms, such as immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, and regulation of immune effector process (Figures 2B,C).

Pearson’s correlation coefficient was used to cluster the immune-related differentially expressed genes in ESCC and normal tissues. The outliers were all removed, and the clustering tree was built (Figures 3A,B). Next, the adjacency matrix and topological overlap matrix were constructed. The six most important modules reflecting the combined function of multiple genes were identified based on the clustering tree and the adjacency and topological overlap matrices (Figure 3C). The blue module with 202 genes and the turquoise module with 222 genes were highly related to the development of ESCC, so these modules were selected for further analysis. Then, we combined gene expression data and patient survival information to calculate the relationship between gene expression levels and days of survival. We identified 18 genes that were related to ESCC patient prognosis, 11 of which were related to good prognosis (RELB, HNRNPL, ITCH, ILF3, PSMC4, CBL, SKP2, DHX33, PRKDC, ZMYND11, and CCDC88A) and seven that were related to poor prognosis (PSME2, PSMD6, GBP2, RABEP2, ZC3HAV1, GRN, and STC2) (Figures 3D, 4A–R).

The 81 ESCC selected samples were randomly divided into two groups, training and test. Then, the risk coefficients of 18 prognostic genes were calculated by COX regression analysis. The sum of the product of risk coefficient and expression level for each gene was used to determine the risk degree of ESCC. The R survival package was then used to establish the prognostic risk model. Finally, we found that the prognostic risk model composed of seven of the 18 prognostic genes had high accuracy (PSMD6, RABEP2, GRN, STC2, ITCH, ILF3, and PSMC4). The expression levels of PSMD6, RABEP2, GRN, and STC2 were related to a poor ESCC prognosis, while the expression levels of ITCH, ILF3, and PSMC4 were related to a good ESCC prognosis. According to the prognostic risk model, we divided the training and test groups into high-risk and low-risk subgroups. For the training group, patients in the low-risk group had better survival rates, with similar results in the test group (Figures 5A,B). Independent prognostic analyses were then performed on the training and test groups. The results suggested that the risk score was related to ESCC prognosis in both groups (Figures 5C,D). We performed GSEA to verify the accuracy of the established prognostic risk model in predicting ESCC prognosis. There were significant differences in GSEA enrichment between the high-risk and low-risk groups. The GSEA enrichment in the high-risk group was mainly related to immune disease and immune response, while this enrichment in the low-risk group was mainly related to tumor disease and tumorigenesis (Figures 5E,F). The gene mutations between the high-risk and low-risk groups were also different (Supplementary Figure S2).

To examine the immune status of low-risk and high-risk group patients, CIBERSORT was used to evaluate the distribution of 23 immune cell types in the ESCC tumor microenvironment. The proportions of immune cells in the low-risk and high-risk groups are shown in Figure 6A. There was abundant immune cell infiltration in both the high-risk and low-risk groups. We then compared the immune cell composition between these groups, finding increased levels of M1 macrophages and resting dendritic cells in the high-risk groups (Figure 6B). Accordingly, we analyzed the relationship between the immune cell composition and survival time, observing that high levels of resting memory CD4 T cells, M0 macrophages, and regulatory T cells were related to a better survival time (Figures 6C–E). High infiltration of M1 macrophages, M2 macrophages, CD8 T cells, and resting dendritic cells showed opposite effects, suggesting that M1 macrophages and resting dendritic cells were increased in malignant ESCC patients and related to a poor prognosis (Figures 6F–I). Interestingly, high levels of M1 and M2 macrophages were related to poor prognosis, while high amounts of M0 macrophages showed the opposite. These data imply that the different states of macrophages are related to ESCC prognosis.

The presence of functional immune cells and functional molecules in the high-risk and low-risk groups were also compared. The results showed that a large number of functional immune cells and functional immune molecules were increased in the high-risk group (Figure 7A). Consistent with these data, high levels of functional immune cells and functional immune molecules were related to poor survival. In particular, patients with low levels of Chemokine receptors (CCR), immune checkpoint molecules, plasmacytoid DC (pDCs), dendritic cells (DCs), and T cell co-stimulation had better survival than those with high levels (Figure 7B).

We examined the Tumor Immune Dysfunction and Exclusion (TIDE) scores of the high-risk and low-risk groups to further analyze the difference of immune status between the two groups. The results showed that the TIDE scores were higher in the high-risk group than in the low-risk group (Figure 8A). The dysfunction score of high-risk patients was higher than that of low-risk patients, while the exclusion score showed the opposite trend (Figures 8B,D). These data indicate that the immune escape risk and poor prognosis were higher in the high-risk group. The high degree of microsatellite instability (MSI) in the high-risk group also confirmed this phenomenon (Figure 8C). After comparing the area under the ROC curve (AUC) for the risk score, TIDE score, and Tumor Inflammation Signature (TIS) score, we found the TIDE score to have the best ability of predicting survival time in ESCC patients (Figure 8E). The risk score had a good predictive effect on the 1, 2, and 3-years survival of ESCC patients, and the AUCs were 0.852, 0.842, and 0.931, respectively (Figure 8F). In conclusion, we found 18 immune-related prognostic genes related to the occurrence of ESCC and used several of them to establish a prognostic model for predicting disease severity.

We analyzed the relative mRNA expression levels of seven prognostic genes in the normal human esophageal epithelial cell (HEEC) and two ESCC cell lines (TE13, EC109) by qRT-PCR. As shown in Figure 9A, there were obvious differences in the expression of these genes. Meanwhile, we analyzed the expression of seven prognostic genes in 19 ESCC cell lines in the CCLE database. We found that the expression trends of these seven genes in poorly differentiated (TE5, TE9, KYSE70, KYSE150, KYSE410) and well differentiated (TE1, TE4, TE6, TE10, TE15, KYSE30, KYSE180, KYSE270, and KYSE510) ESCC cell lines were nearly consistent with our previous analysis (Figures 9B,C).