A novel immune-related gene signature predicts survival in esophageal squamous cell carcinoma.

BACKGROUND: Immune-related genes (IRGs) are highly relevant to the progression and prognosis of esophageal squamous cell carcinoma (ESCC). A prognostic signature could be reliable in stratifying ESCC patients according to the risk score, which may help manage systematic treatments. In this study, a systematic and reliable immune signature was developed to estimate the prognostic stratification in ESCC.
METHODS: Ribonucleic acid (RNA) expression data of 79 ESCC samples from the Cancer Genome Atlas (TCGA) database and 269 normal esophageal mucosal samples from the Genotype-Tissue Expression (GTEx) project database were downloaded from the University of California, Santa Cruz (UCSC) website to form a TCGA-GTEx dataset. First, we screened differentially expressed genes (DEGs) and then filtered IRGs based on the Immunology Database and Analysis Portal (ImmPort) database to obtain immune-related DEGs (IRDEGs). Next, a novel prognostic signature based on IRDEGs was developed using multivariable Cox analysis. Immune infiltration status was evaluated via single-sample gene set enrichment analysis (ssGSEA). ESCC tissues were grouped into three clusters in terms of immune infiltration (Immunity-L, Immunity-M, and Immunity-H) by applying an unsupervised hierarchical clustering algorithm. Finally, the samples were divided into high- and low-risk groups using the median of the risk score scores for GSEA pathway enrichment analysis in the three clusters.
RESULTS: The prognostic signature based on IRDEGs (FCER1G, ISG20, and EGFR) performed moderately in prognostic predictions, with a concordance index (C-index) value of 0.73 [95% (confidence interval) CI: 0.63-0.84, P=2.02E-05] and an area under the curve (AUC) value of 0.817. The xenobiotic metabolism pathway was significantly enriched and up-regulated both in the high-risk group of the immunity-M and immunity-H clusters.
CONCLUSIONS: The novel immune-related prognostic signature we constructed has a good prognostic, predictive ability and can be used as an independent prognostic indicator. Our study provides clinicians with a quantitative tool to predict the probability of individual survival time and helps clinicians select targets for immunotherapies and individualized treatment strategies for ESCC patients. 2021 Translational Cancer Research. All rights reserved.

Introduction

Esophageal cancer (EC) is the sixth leading cause of cancer death worldwide (1) and has two common subtypes: esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC). ESCC is the most frequent subtype of EC (2). Although the emergence of new chemotherapeutic drugs and regimens has improved the overall survival (OS) rate of ESCC patients, drug resistance has significantly worsened patients’ prognostic outcome. At present, the pathogenesis of ESCC is unclear. Therefore, there is a pressing need for effective biomarkers and targeted therapeutics to improve patients’ prognostic outcome.

Cancer immunotherapy has been a major driver of personalized medicine, with aggressive efforts to leverage the immune system to fight tumors (3,4). It has shown encouraging results as an emerging treatment for some cancers, including melanoma (5), hepatocellular carcinoma (6), and non-small cell lung cancer (7). Tumor cells in the tumor microenvironment (TME) can mimic immune cells’ functions to affect the composition and function of host immune cells and induce immunosuppression by overexpressing IRGs, thereby promoting the proliferation and spread of tumor cells (8). The acquisition of immune functions maintains tumor cells that can survive in the immunosuppressive TME. Therefore, patients with different immune infiltration states may have different prognoses. Meanwhile, reliable biomarkers or stratification systems that can be used for more accurate prediction are essential in the TME.

Several immune-related parameters or signatures that predict ESCC patients’ prognosis have been reported (9-11), further suggesting that distinct immune infiltration status has a profound influence on ESCC patient outcomes. Herein, we established a predictive risk model based on IRGs to estimate patient outcomes in ESCC. Importantly, we evaluated this model’s predictive ability in assessing OS in different immune infiltration patterns and screened out pathways closely related to these patterns. We performed a systematic investigation of the immune phenotype within the ESCC microenvironment to better understand the complex antitumor response and guide effective immunotherapies in ESCC. We present the following article in accordance with the Materials Design Analysis Reporting (MDAR) checklist (available at http://dx.doi.org/10.21037/tcr-20-2665).

Methods

Data from the UCSC Xena website

Log2(FPKM+1) format ribonucleic acid (RNA)-sequencing data used in this study were downloaded from the University of California, Santa Cruz (UCSC) Xena website (https://xenabrowser.net/), which provides comprehensive information retrieval services based on the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) project databases. We collected and processed transcription data from 79 ESCC samples from the TCGA database and the expression data of 269 normal esophageal mucosa samples from the GTEx project database. For comparative analysis between ESCC samples and many normal samples, we combined the TCGA dataset with the GTEx dataset (TCGA-GTEx dataset) for subsequent analysis. Moreover, the corresponding survival and clinical information files of 79 patients with ESCC were also downloaded from the UCSC. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Screening immune-related differential genes

The limma R package was used to screen differentially expressed genes (DEGs) between normal and ESCC samples. |Log2FC| ≥2 and PFDR<0.05 were used as the interception criteria for screening DEGs. We also derived a list of 1,811 IRGs via the Immunology Database and Analysis Portal (ImmPort, https://immport.niaid.nih.gov) database. We then intersected DEGs with IRGs to obtain IRDEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were simultaneously conducted to explore potential molecular mechanisms of the IRDEGs, which were performed by the R bioconductor package clusterProfiler. PFDR<0.01 was set as the enrichment cut-off to screen meaningful enrichment results.

Construction of the immune-related prognostic signature

To provide clinicians with a quantitative tool to predict the individual probability of survival time, we developed an immune-related prognostic signature. A univariate Cox analysis was performed to select IRDEGs with a significant prognostic value (P<0.05) as candidates for the prognostic signature. Next, these selected IRDEGs were further refined with stepwise forward selection by employing the least absolute shrinkage and selection operator (LASSO) method. Finally, a multi-gene prognostic signature was constructed by multivariate Cox analysis, and an individual risk score was calculated based on a combination of coefficients and the expression status of these refined IRDEGs. The risk score was set as:

Ri, N, expression, and βgene i represented the patient’s risk score, the number of these refined IRDEGs, the gene expression level, and the coefficient value, respectively. The median value was used as the cutoff value for dividing high- and low-risk groups. Subsequently, Kaplan-Meier survival analysis, principal component analysis (PCA), and univariate and multivariate Cox analysis, were performed to evaluate the predictive value of this prognostic signature.

Gene set enrichment analysis (GSEA)

ESCC cases were divided into two risk groups according to the median value of the risk scores. GSEA was performed using Broad Institute GSEA software 4.1.0 (https://www.gsea-msigdb.org/gsea/downloads.jsp) in three immune clusters. The annotated gene set “h.all.v7.1.symbols.gmt” was selected and used for the enrichment analysis. Permutations were set to 1,000 to obtain a normalized enrichment score (NES). A normal P alue <0.05 and NES >1.5 were considered significantly enriched. We chose all hallmark gene sets with significant enrichment and displayed gene sets enrichment plots. The xenobiotic metabolism pathway was significantly enriched and up-regulated in both high-risk groups of the immunity-M and immunity-H clusters. Effective genes enriched by this pathway in two clusters were then intersected, and 34 significant correlation hub genes were imported to construct the protein-protein interaction (PPI) network.

Protein-protein interaction network (PPI)

STRING (https://string-db.org) was used for PPI network construction and hub gene screening. The PPI network was constructed by setting medium confidence at 0.400. After the irrelevant genes were excluded, 23 genes were used for constructing the PPI network.

Statistical analysis

All analyses were conducted using R software (version 3.5.3, https://www.r-project.org). PCA was also performed by using the ggplot2 package of R. Clinical correlation analyses were tested by the Student’s t-test or Wilcoxon test for two groups, and by one-way analysis of variance (ANOVA) or Kruskal-Wallis test for more than two groups. The Cox proportional hazards regression model was used for univariate and multivariate analyses. The performance of the multivariate Cox regression model was evaluated by the concordance index (C-index) using the survcompr software package. The area under the curve (AUC) of the multi-index receiver operating characteristic (multi-ROC) curve was calculated via the survivalROC R package to validate the prognostic signature’s performance. All tests were two sided, and P-values of less than 0.05 was considered statistically significant.

Results

Clinical characteristics

Seventy-nine patients diagnosed with ESCC were registered in the TCGA cohort, including 68 male (86.08%) and 11 female (13.92%) patients. Detailed clinical characteristics of all patients are presented in .

Table 1

Detailed clinical characteristics of all patients are listed below

Variables	Case, n (%)
Age at diagnosis (year)
≤55	34 (43.04)
>55	45 (56.96)
Gender
Male	68 (86.08)
Female	11 (13.92)
Histologic grade
G1	15 (18.99)
G2	38 (48.10)
G3	19 (24.05)
NA	7 (8.86)
Alcohol history
No	19 (24.05)
Yes	58 (73.42)
NA	2 (2.53)
Pathologic-M
M0	70 (88.61)
M1	3 (3.80)
NA	6 (7.59)
Pathologic-N
N0	45 (56.96)
N1	26 (32.91)
N2	5 (6.33)
N3	1 (1.27)
NA	2 (2.53)
Pathologic-T
T1	7 (8.86)
T2	27 (34.17)
T3	41 (51.90)
T4	3 (3.80)
NA	1 (1.27)
Pathologic-stage
I	6 (7.59)
II	47 (59.49)
III	22 (27.85)
IV	3 (3.80)
NA	1 (1.27)

NA, not available.

Screening of immune-related genes

As shown in , we initially identified a total of 1,575 DEGs in ESCC tissues compared with normal tissues (P<0.05, |log2FC| >2). Next, we intersected 1,811 IRGs from ImmPort with 1,575 DEGs from the TCGA-GTEx dataset. As a result, a total of 124 IRDEGs, among which 67 genes were up-regulated while 57 genes were down-regulated, were obtained and displayed using a Venn diagram ().

Figure 1

Differentially expressed genes. Heatmap (A) and volcano plot (B) demonstrating differentially expressed genes (DEGs) between ESCC and normal tissues. Brown dots represent differentially down-regulated genes, purple dots represent differentially up-regulated genes, and black dots represent no differentially expressed genes in the TCGA-GTEx dataset. (C) Venn diagrams showing the number of up-regulated or down-regulated IRDEGs. The top 10 significantly enriched GO-BP terms (D) and KEGG terms (E) of 124 IRDEGs. The red and blue dots in the graph represent up-regulated and down-regulated genes, respectively. The depth of the inner arc area shows the decrease or increase of the biological process.

The 124 IRDEGs examined in this study played various roles in immune functions, and many of them were involved in the TME. Through Gene Ontology (GO) enrichment analysis, we found that the top 10 GO biological process (BP) terms were mainly related to receptor signaling pathways (). The most enriched GO BP term correlated with those genes was the “immune response-activating cell surface receptor signaling pathway”, which contained 36 IRDEGs. The top 10 KEGG enrichment terms are presented in . Epstein-Barr virus infection and cytokine-cytokine receptor interaction were the most frequently identified KEGG pathway. The results indicate that these IRGs directly or indirectly affect tumor cells in the TME by acting on different signaling pathways.

Construction and evaluation of the prognostic signature based on IRDEGs

In our univariate Cox regression analysis, seven genes (FCER1G, PSMD6, ISG20, IRF9, CD14, EGFR, and PLXNA1) were selected as candidates for developing the predictive model (). We then chose the appropriate candidates using the LASSO method (), and subsequently constructed an immune-related prognostic signature with a C-index of 0.73 [95% confidence interval (CI): 0.63–0.84, P=2.02E-05] to predict the outcome of ESCC patients by multivariate Cox regression analysis (). The risk score distribution and survival status were ranked by the risk scores in the ESCC set (shown in ). The heatmap showing this signature’s expression profiles in the low- and high-risk groups is shown in . Through PCA analysis based on this immune signature, we found that high- and low-risk groups could be well distinguished (). Patients in the low-risk group had a longer OS compared to those in the high-risk group (). These results indicate that this prognostic signature exhibits good prognostic, predictive ability.

Figure 2

Development of the prognostic signature based on IRGs. (A) Univariate Cox analysis based on immune-related genes showed that seven genes (FCER1G, PSMD6, ISG20, CD14, IRF9, PLXNA1, and EGFR) have prognostic value. (B) LASSO coefficient profiles of the selected immune-related gene. (C) Partial likelihood deviance for LASSO coefficient profiles. (D) Multivariate Cox analysis was used to construct a prognostic model with three genes included (FCER1G, ISG20, and EGFR).

Figure 3

The predictive ability of this immune signature for the OS of ESCC patients was confirmed. (A) Risk score distribution. (B) Survival overview. (C) Heatmap showing the expression profiles of the signature in the low- and high-risk groups. (D) Principal Component Analysis (PCA) based on the immune signature. (E) Kaplan-Meier curves between the high- and low-risk score groups. Patients in high-risk group had shorter OS. (G) Multi-index ROC analysis. The AUC value of the prognostic signature was the highest (AUC =0.817).

Following this, we evaluated the prognostic, predictive value of five clinicopathological characteristics, including T stage, lymph node metastasis (N), distant metastasis of tumor (M), histologic grade (G), and stage, and determined the risk scores by Cox regression analysis. The results of the univariate Cox regression analysis suggested that the risk score (P<0.001) and N (P=0.036) were associated with prognosis, whereas the other clinicopathological factors were not (). Multivariate Cox regression analysis suggested that only the risk score could become an independent predictor after the other parameters were adjusted (). Furthermore, through multi-index ROC analysis, we affirmed that the diagnostic efficiency of this immune signature (AUC =0.817) was better than the above five clinicopathological characteristics in predicting survival, including T (AUC =0.515), N (AUC =0.635), M (AUC =0.541), G (AUC =0.582), and stage (AUC =0.581) (). Put simply, this signature is an independent prognostic indicator and is superior to the other five clinicopathological characteristics.

Table 2

Univariate and multiple regression analysis of esophageal squamous cell cancer

Variables	Univariate analysis			Multivariate analysis
	Hazard ratio (95% CI)	P value		Hazard ratio (95% CI)	P value
G	1.00 (0.55–1.82)	0.999		1.30 (0.59–2.86)	0.518
T	1.05 (0.57–1.93)	0.868		1.99 (0.70–5.68)	0.200
N (yes/no)	2.67 (1.07–6.69)	0.036*		4.70 (0.98–22.59)	0.053
M (yes/no)	3.20 (0.91–11.24)	0.070		6.36 (0.68–59.89)	0.106
stage	2.22 (0.96–5.15)	0.063		0.59 (0.11–3.05)	0.525
RiskScore	1.30 (1.13–1.49)	<0.001*		1.32 (1.11–1.56)	0.002*

*, significant results.

We also explored the clinical significance of the included genes and the risk score of the prognostic signature. The results showed that the expression of FCER1G was significantly higher in advanced G stage cases (), the expression of ISG20 was significantly lower in advanced T stage cases (), and the risk score was significantly higher in advanced G stage cases and in males () (all P<0.05).

Figure 4

Clinical correlation analysis. (A) Association between FCER1G and G stage. (B) Association between ISG20 and T stage. (C) Association between risk score and G stage. (D) Association between risk score and gender.

Evaluation of the prognostic, predictive ability of this prognostic signature in three immune clusters

The enrichment of immune cells and immune-related pathways were estimated using the single-sample GSEA (ssGSEA) algorithm. ESCC samples were clustered into low, intermediate, and high immune infiltration clusters (Immunity_L, Immunity_M, and Immunity_H) using the hierarchical clustering method (). The immunity_H cluster showed a higher risk score than both the immunity_M and immunity_L clusters (). Similarly, patients with a low risk score had a lower proportion of high immune infiltration than patients with a high risk score (). High immune infiltration in patients was linked to a high risk score and was associated with an unfavorable outcome (). The association between OS and the different clusters of ESCC was analyzed by log-rank test (). Due to the limited sample size, survival prediction was significant only in the immunity_M cluster. However, it can be seen from the results that the immunity_M and immunity_H clusters had an unfavorable survival probability compared with the immunity_L cluster.

Figure 5

Heterogeneous immune cell infiltration in the low- and high-risk score groups. (A) The heatmap showing the expression level of immune cells and immune-related pathways under different immune infiltration clusters. (B) The distribution of risk scores in three immune infiltration patterns. (C) The distribution of immune infiltration patterns in the low- and high-risk groups. (D) Sankey diagram of immune infiltration patterns in groups with different risk scores and survival outcomes. Kaplan-Meier curves between the high- and low-risk score groups in the immunity_L (E), immunity_M (F) and immunity_H (G) clusters, respectively.

Functional annotation and screening hub genes

To identify the constructed immune signature’s underlying biological characteristics, GSEA analysis was performed based on the high- and low-risk score groups in the three immune clusters, respectively. No significant pathway was enriched in the immunity_L cluster. The GSEA results showed that only xenobiotic metabolism was significantly up-regulated in the high-risk score group of the immunity_M cluster, which is the same result observed in the high-risk score group of immunity_H cluster. In contrast, the mitotic spindle, glycolysis, phosphatidylinositol 3-hydroxy kinase/protein kinase B/mammalian target of rapamycin (PI3K/AKT/MTOR) signaling, G2M checkpoint, and unfolded protein response pathways were significantly down-regulated in the high-risk score group of the immunity_H cluster ().

Figure 6

GSEA revealed that most significant hallmarks correlated with our immune signature in three immune infiltration clusters separately. The enriched pathways were up-regulated in the high-risk score group (A) and in the low-risk score group (B,C,D,E,F) in the immunity_H cluster. (G) The enriched pathway was up-regulated in high-risk tumors in the immunity_M cluster. A normal P value <0.05 and NES >1.5 were considered significantly enriched. (H) Effective genes enriched by the xenobiotic metabolism pathway in two clusters were intersected and depicted by Venn diagrams. (I) PPI analysis of intersection genes and depiction by STRING.

Next, effective genes enriched by the xenobiotic metabolism pathway in two clusters were intersected, and 34 significant correlation hub genes were imported to construct the PPI network (). Ultimately, the PPI network was constructed by 23 hub genes that were significantly related to this pathway (). We found that numerous genes, such as ADH7 (12), AKR1C2 (13), ALDH2 (14), NQO1 (15), ADH1C (12), EPHX1 (16), BLVRB (17), HMOX1 (18), ABCC3 (19), ALDH3A1 (20), and IDH1 (21,22), are related to ESCC tumor susceptibility, tumorigenesis, and progression.

Discussion

Cancer immunotherapy is a promising new avenue for ESCC treatment that acts primarily on the TME. Also, the heterogeneity of the TME is closely related to the efficiency of immunotherapy. In the present study, we attempted to identify an immune-related signature that contributes to the OS of ESCC. Using multiple indicators, we first evaluated its ability to predict prognosis in all samples. We discovered that this signature not only has good prognostic predictive ability but can also be used as an independent prognostic indicator. Importantly, we also assessed its ability to predict prognosis under different immune infiltration states, and used GSEA to determine which pathways were significantly enriched in the high- and low-risk groups. We found that this signature’s ability to predict prognosis was different under different infiltration patterns, and the pathways of significant enrichment in the high- and low-risk groups were also non-identical.

A multi-gene risk prediction model based on FCER1G, EGFR, and ISG20 was established in this study. FCER1G and ISG20 [also called cluster of differentiation (CD)25] are two risk genes, and EGFR is a protective gene (). Pauline et al. reported that FCER1G is essential for tumor development and squamous carcinogenesis by B cells, humoral immunity and activating FCER1G to establish chronic inflammatory programs (23). FCER1G is also engaged in many immune responses and plays a tumor-promoting role in numerous types of tumors, such as meningioma (24) and childhood leukemia (25). At present, there are no directly-related reports concerning the association between FCER1G and ESCC; however, it is generally conceded that chronic inflammation is closely related to oncogenesis. Therefore, we might assume that the overexpression of FCER1G regulates immune-related pathways of tumorigenesis. Furthermore, Guan et al. found that patients with ESCC had a significant increase in circulating anti-CD25 immunoglobulin G (IgG) antibodies compared to control subjects and that patients with early stage ESCC (stage I) had the highest expression of circulating IgG (26). In this study, we discovered that patients with T1 stage ESCC had the highest expression of ISG20 (). So, ISG20 may be considered as a potential biomarker for early diagnosis and prediction of ESCC recurrence or survival. EGFR expression was usually associated with a more advanced tumor stage and a reduced OS in patients with esophageal and esophagogastric junction adenocarcinomas (27). ESCC predominantly involves EGFR gene copy number alterations and protein overexpression (28-31), with little EGFR mutation (32-34). Furthermore, the literature regarding EGFR expression in ESCC contains conflicting data on the relationship between overexpression and survival (31,35,36). Although several clinical studies show the beneficial effects of EGFR inhibitors in the treatment of ESCC (37,38), evidence supporting the validity of EGFR-targeting therapies for esophageal cancer (39) is not robust. Yoshioka et al. discovered that one of the factors affecting the therapeutic effect of EGFR inhibitors in ESCC cells is cell phenotype. Mesenchymal-like ESCC cells are resistant to EGFR inhibitors because EGFR signaling is not blocked. On the contrary, EGFR signaling is affected by EGFR inhibitors, which exhibit antitumor effects on epithelial-like ESCC cells accompanied by the promotion of squamous cell differentiation. In the future, more potential mechanisms need to be studied to determine the effect of EGFR targeted therapy for ESCC (40).

We found that the OS of patients in the high-risk group was poor through model analysis in all samples compared with those in the low-risk group (). Patients in the high-risk group received more precise treatment options to improve their prognosis. We further evaluated the signature’s predictive ability for the OS of patients in different immune infiltration states. We discovered that the predictive ability and the predictive results of this signature could be different in different immune infiltration states (). Patients in the high-risk score group had a higher infiltration status and poorer survival than patients in the low-risk score group in the immunity_H and immunity_M clusters. However, further analysis of the prognostic evaluation for the OS of patients with low immune infiltration status in this paper is needed.

To understand which gene sets can be enriched in the high- and low-risk groups in different immune infiltration clusters, GSEA analysis was performed based on the high- and low-risk group in three immune clusters, respectively. Only xenobiotic metabolism was significantly up-regulated in the high-risk group of the immunity_M () and immunity_H () clusters. We found that some studies have been conducted on the association between xenobiotic metabolizing genes polymorphism and ESCC. Several studies have linked polymorphisms in genes encoding phase I and phase II xenobiotic metabolizing enzymes to the modulation of risk for ESCC (41-44). Based on the PPI network, we found numerous genes significantly related to the xenobiotic metabolism pathway, which are related to ESCC tumor susceptibility, tumorigenesis, and progression. For instance, Wang et al. have demonstrated that the ADH1B-ADH1C-ADH7 cluster single nucleotide polymorphisms (SNPs) confer susceptibility to ESCC in two case-control sets (12). Also, Ren et al. noted that ESCC tumors frequently expressed heme oxygenase 1 (HO-1), and knockdown of HO-1 promoted apoptosis through activation of a reactive oxygen species (ROS)-mediated caspase apoptosis pathway (18). Zhang et al. revealed that AKR1C2 could function as an oncogene (via activation of the PI3K/AKT pathway) and a novel prognostic biomarker and/or a potential therapeutic target in ESCC (13). In contrast, the mitotic spindle, glycolysis, PI3K/AKT/MTOR signaling, G2M checkpoint, and unfolded protein response pathways were significantly up-regulated in the low risk score group of the immunity_H cluster (). These pathways could be closely related to the good prognosis of the tumors. Shang et al. has confirmed that down-regulation of BIRC5 as a mitotic spindle checkpoint gene inhibits the migration and invasion of ESCC cells by interacting with the PI3K/AKT signaling pathway (45).

There are some limitations to our research that should be noted. This study is primarily based on database analyses, and lacks clinical practice confirmation. Thus, further verification through clinical trials is needed in the future.

Conclusions

In summary, this study established a stable multi-gene risk prediction and evaluated its predictive value for ESCC patients’ OS. We found several key pathways and hub genes closely related to ESCC patients’ prognosis under different immune infiltration states. Therefore, further verification with clinical trials in the future is needed.