Construction and Validation of an Immune-Related Gene Prognostic Index for Esophageal Squamous Cell Carcinoma

Ji, Qinghua; Cai, Yingying; Shrestha, Sachin Mulmi; Shen, Duo; Zhao, Wei; Shi, Ruihua

doi:https://doi.org/10.1155/2021/7430315

BioMed Research International

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusions Data Availability Disclosure Conflicts of Interest Authors’ Contributions Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2021 | Article ID 7430315 | https://doi.org/10.1155/2021/7430315

Construction and Validation of an Immune-Related Gene Prognostic Index for Esophageal Squamous Cell Carcinoma

Qinghua Ji,¹Yingying Cai,¹Sachin Mulmi Shrestha,^1,2Duo Shen,¹Wei Zhao,^1,3and Ruihua Shi^1,2

Academic Editor: Guoxian Yu

Received27 Jul 2021

Accepted20 Sept 2021

Published21 Oct 2021

Abstract

Immune checkpoint inhibitor (ICI) therapy may benefit patients with advanced esophageal squamous cell carcinoma (ESCC); however, novel biomarkers are needed to help predict the response of patients to treatment. Differentially expressed immune-related genes within The Cancer Genome Atlas ESCC dataset were selected using the weighted gene coexpression network and lasso Cox regression analyses. Based on these data, an immune-related gene prognostic index (IRGPI) was constructed. The molecular characteristics of the different IRGPI subgroups were assessed using mutation information and gene set enrichment analysis. Differences in immune cell infiltration and the response to ICI therapy and other drugs were also analyzed. Additionally, tumor and adjacent control tissues were collected from six patients with ESCC and the expression of these genes was verified using real-time quantitative polymerase chain reaction. IRGPI was designed based on CLDN1, HCAR3, FNBP1L, and BRCA2, the expression of which was confirmed in ESCC samples. The prognosis of patients in the high-IRGPI group was poor, as verified using publicly available expression data. KMT2D mutations were more common in the high-IRGPI group. Enrichment analysis revealed an active immune response, and immune infiltration assessment showed that the high-IRGPI group had an increased infiltration degree of CD8 T cells, which contributed to the improved response to ICI treatment. Collectively, these data demonstrate that IRGPI is a robust biomarker for predicting the prognosis and response to therapy of patients with ESCC.

1. Introduction

Esophageal carcinoma is one of the leading causes of cancer-related death worldwide [1]. The main types of esophageal carcinoma are squamous cell carcinoma (ESCC) and adenocarcinoma, which exhibit different pathomorphological, epidemiological, and molecular characteristics [2]. Moreover, ESCC accounts for approximately 90% of the morbidity and mortality of all esophageal cancers [3]. The promotion of large-scale population screening has effectively improved the survival rate of patients with ESCC; nonetheless, a considerable number of patients with ESCC are in an advanced stage at diagnoses [4]. Only 27% of patients with ESCC benefit from traditional chemotherapy and radiotherapy [5, 6]. Therefore, novel methods to predict the prognosis of patients with ESCC and their corresponding treatment are urgently needed.

An increasing number of researchers have recently focused on immune checkpoint inhibitor (ICI) therapy, which prevents tumor cell immune escape and induces an immune response by inhibiting immune checkpoints, such as programmed death-1 (PD-1), programmed death-ligand 1 (PD-L1), and CTL-associated protein-4 (CTLA-4) pathways [7–9]. Phase III ICI clinical trials, including those on nivolumab and pembrolizumab for ESCC, showed that ICIs can significantly improve the patient survival rate and reduce the incidence of treatment-related adverse events compared with traditional chemotherapy drugs [7–9]. Most patients achieve cancer control upon ICI therapy, although some develop recurrence or drug resistance [10]. Moreover, approximately 10% patients treated with ICIs may have unconventional response patterns (pseudoprogression), which challenges the evaluation of treatment efficacy [11]. Thus, methods for evaluating whether patients are suitable for ICI therapy and assessing ICI therapy efficacy are urgently needed.

The response of patients to ICI therapy is mainly affected by tumor cell-intrinsic factors (such as tumor mutational burden and microsatellite instability) and the tumor microenvironment [12, 13]. Hence, evaluation of the immune status of the tumor microenvironment by immune gene signatures may effectively help predict the benefits of ICIs [14]. The expression level of immune-related genes can predict the response of patients to ipilimumab treatment [15]. Additionally, recent studies showed that immune- or tumor microenvironment-related gene expression scores can predict the survival and response to immunotherapy in hepatocellular carcinoma and lung cancer [16, 17]. Therefore, it is helpful to evaluate the immune microenvironment, prognosis, and response to immunotherapy by examining gene expression in patients before treatment.

In this study, we used immune gene signatures to develop prognostic and ICI therapy indicators for patients with ESCC. We also performed weighted gene coexpression network analysis (WGCNA) and lasso regression analysis to construct an immune-related gene prognostic index (IRGPI). The molecular and immune characteristics of the IRGPI subgroups were evaluated, and the potential of IRGPI for assessing immunotherapy efficacy in patients with ESCC was determined. The study design is shown in Figure 1.

2. Materials and Methods

2.1. Collection of Patient Information and Databases

ESCC transcriptome data, clinical information, and gene mutation data were downloaded from The Cancer Genome Atlas (TCGA) database, which included 81 tumor and 11 adjacent noncancerous samples. Transcriptome data and clinical information of the validation cohort GSE53625 with 179 ESCC tumor samples and 179 adjacent normal tissues were downloaded from the NCBI Gene Expression Omnibus database [18]. Immune-related gene lists were downloaded from InnateDB (https://www.innatedb.com), ImmPort, and IRIS [19–21].

2.2. Identification of Immune Genes Correlated with Prognosis

Differentially expressed ESCC genes were identified using the R package limma (version 3.44) based on TCGA transcriptome data, with a false discovery and fold [22]. The differentially expressed immune ESCC genes were determined after the intersection of immune-related gene sets (InnateDB) with differentially expressed genes in ESCC. Gene Ontology (GO) analysis of these genes was performed using the R package clusterProfiler (version 3.17.5) [23].

The WGCNA (version 1.46) method was used to identify hub genes that were significantly associated with ESCC [24]. We used the scale-free topology criterion to determine the soft threshold of . Under this selection, the scale-free topology fitting index . Based on the gene expression matrix, the similarity of gene expression was calculated to obtain the adjacency matrix, which was then transformed into a topological overlap matrix. Genes were grouped by hierarchical clustering and then divided into different expression modules according to the coexpression pattern. We then calculated the correlation between these gene modules and ESCC occurrence. Modules with an absolute value of the correlation were selected for further analysis.

Univariate prognostic analysis was performed for genes in the selected modules with the R package survival (version 3.2). The R package glmnet (version 4.0) was used for lasso analysis of survival-related genes in univariate analysis [25]. The immune-related genes selected using lasso analysis were used to construct the prognosis index model.

2.3. IRGPI Construction and Reliability Evaluation

For immune-related genes included in the IRGPI, the R package survminer (version 0.4.8) was used to determine the optimal cutoff expression value for prognosis and the logrank test and Kaplan–Meier survival curve were used to determine its relationship with overall survival [26, 27]. Based on the immune-related genes selected, multivariate regression analysis using the R package rms (version 6.1) was performed to construct the IRGPI model. The IRGPI score of each sample was assessed by calculating the sum of each immune-related gene expression value multiplied by its corresponding regression coefficient based on the Cox proportional hazard model. The risk score plot of IRGPI was obtained using the R package ggrisk (version 1.2). The survivalROC (version 1.0.3) package was used to test the diagnostic significance of IRGPI in patients with ESCC in different years. TCGA patient data were divided into two groups according to the IRGPI median value. The Kaplan–Meier survival curve and logrank test were used to evaluate the prognosis value of IRGPI in TCGA ESCC cohort, which was further confirmed in the GSE53625 dataset.

2.4. Identification of Mutation and Enrichment Characteristics in Different IRGPI Subgroups

Based on the median value of IRGPI (0.13), TCGA ESCC patients were divided into high- and low-IRGPI groups. The R package limma (version 3.44) was used to analyze the differentially expressed genes between groups. The R package maftools (version 2.6.05) was used to summarize and visualize mutation information between subgroups [28]. The R package Pi (version 2.2.1) was used for gene set enrichment analysis (GSEA) between IRGPI subgroups based on the MSigDB hallmark gene set [29].

2.5. Evaluation of Immune Cell Infiltration

The ssGSEA function of the R package GSVA (version 1.36) was used to calculate the enrichment score of 28 types of immune cells for each sample [30]. The R packages ggpubr (version 0.4) and ggplot2 (version 3.3.0) were used to compare and visualize the immune cell enrichment score between groups. The R package ComplexHeatmap (version 2.7.9) was used to display the landscape map of the relationship between the infiltration of immune cells and clinical information of the samples [31].

2.6. Comparison of Immunotherapy Effectiveness and Chemotherapeutic Response

The SubMap module in GenePattern (https://cloud.genepattern.org/gp) was applied to predict the response effectiveness of the IRGPI subgroups to immunotherapy [32]. Expression data from patients with melanoma who responded to CTLA-4 and PD-1 therapy were compared with those of TCGA ESCC samples upon Bonferroni correction [33].

According to the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/), we predicted the different reactions of the two IRGPI subgroups to chemotherapy. Based on the gene expression of samples and GDSC training set, the half-maximal inhibitory concentration (IC₅₀) of each chemotherapy was evaluated by ridge regression using the R package pRRophetic (version 0.5) [34]. Tenfold crossvalidation was used to ensure prediction accuracy.

2.7. ESCC Tissue Collection

After obtaining ethics approval (no. 2019ZDSYLL023-Y01) from Zhongda Hospital, ESCC and adjacent-control tissues were collected from six patients with ESCC. The adjacent tissues were collected from the esophageal tissues that were more than 2 cm and less than 5 cm away from the ESCC tissues. Fresh ESCC and adjacent normal tissues were collected during surgery. After rapid freezing in liquid nitrogen, all tissue samples were stored at −80°C. The obtained tumor tissues were pathologically verified.

2.8. RNA Extraction and Quantitative Polymerase Chain Reaction (qPCR)

The FastPure Cell/Tissue Total RNA Isolation Kit (Vazyme, Nanjing, China) was used according to the manufacturer’s instructions to extract total RNA from patient tissue samples. cDNA was synthesized using HiScript III RT SuperMix for qPCR (Vazyme). Using ChamQ SYBR qPCR Master Mix (Vazyme), qPCR was performed. The gDNA filter columns in the RNA extraction kit and subsequent gDNA wiper mix before reverse transcription ensured that there was little or no gDNA residue in the qPCR system. Relative gene expression was calculated using the method, with GAPDH as an internal reference. All samples were evaluated three times. The corresponding primers used are listed in Supplementary Table S1.

2.9. Statistical Analyses

All statistical analyses were implemented using R (3.6.1 version). The Wilcoxon rank test was used to verify the statistical significance between continuous variables, and the chi-squared test was used to compare classified variables. A value < 0.05 was considered to indicate statistically significant results.

3. Results

3.1. Immune-Related Genes in Esophageal Carcinoma

In total, 81 tumor and 11 healthy tissue samples from TCGA ESCC were used for differential expression analysis and 4838 differentially expressed genes were screened (false discovery , as a standard). Cross-analysis with immune-related gene databases (Immport, IRIS, and InnateDB) revelated 1271 differentially expressed immune-related genes in the TCGA ESCC dataset. Overall, 713 and 558 genes were up- and downregulated, respectively. GO enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses further showed that these differentially expressed immune-related genes were significantly correlated with 479 GO terms and 56 KEGG pathways (Supplementary Figure S1 and Table S2).

Further assessment of the differentially expressed immune-related ESCC genes by WGCNA according to the coexpression patterns revealed eight modules. Based on the correlation coefficient between each gene module and ESCC, genes in the brown and turquoise modules (correlation coefficient with ) were selected for further analysis (Supplementary Figure S2).

The 547 genes in the brown and turquoise modules were analyzed by univariate Cox regression, which revealed that 25 genes were associated with the prognosis of ESCC. By lasso Cox regression, four genes (CLDN1, HCAR3, FNBP1L, and BRCA2) were selected as the best model (Supplementary Figure S3). Using the maximally selected test statistics from R package survminer, we determined the best cutoff points for survival analysis of the above four genes. The cutoff points of CLDN1, HCAR3, FNBP1L, and BRCA2 for survival analysis were 9.39, 2.18, 5.17, and 4.98, respectively. The results showed that high expression of CLDN1 and HCAR3 was associated with poor prognosis, whereas high expression of FNBP1L and BRCA2 was associated with good prognosis in the TCGA ESCC cohort (Figure 2).

(a)

(b)

(c)

(d)

(e)

3.2. Prognostic Value of IRGPI in TCGA and the Validation ESCC Cohort

We used the aforementioned four genes to construct an immune-related gene prognosis model of ESCC by multivariate Cox regression. The formula for IRGPI was as follows:in which the correlation coefficient from the Cox proportional hazard model of each gene was multiplied by the expression of the corresponding gene. The TCGA ESCC cohort was divided into high- and low-IRGPI groups according to the risk score, with the median IRGPI defined as the cutoff point (Figure 3(a)). Using the median IRGPI of 0.13 as the cutoff value, the high- and low-IRGPI groups showed different prognosis outcomes, immune infiltration, mutation characteristics, and responses to chemotherapy and immunotherapy. The expression of CLDN1 in the high-IRGPI group was higher than that in the low-IRGPI group (Student’s -test, ), and the median expression levels were 8.56 and 7.4, respectively. The expression of HCAR3 in the high-IRGPI group was higher than that in the low-IRGPI group (Student’s -test, ), and the median expression levels were 1.89 and 0.51, respectively. The expression of FNBP1L in the high-IRGPI group was lower than that in the low-IRGPI group (Student’s -test, ), and the median expression levels were 4.90 and 5.48, respectively. The expression of BRCA2 in the high-IRGPI group was lower than that in the low-IRGPI group (Student’s -test, ), and the median expression levels were 4.13 and 4.65, respectively.

(a)

(b)

(c)

(d)

(e)

The prognostic value of IRGPI in the TCGA ESCC cohort was evaluated using a time-dependent receiver operating characteristic (ROC) curve. The area under the curve (AUC) values of IRGPI at 1, 3, and 5 years were 0.791, 0.807, and 0.845, respectively (Figure 3(b)). Moreover, Kaplan–Meier analysis showed that the prognosis of the high-IRGPI group in the TCGA ESCC cohort was significantly worse than that of the low-IRGPI group (logrank test, ) (Figure 3(c)). The prognostic value of IRGPI was confirmed in the ESCC cohort (GSE53625), for which the AUC values of 1, 3, and 5 years in the ROC curve were 0.837, 0.841, and 0.893, respectively, and the prognosis of the high-IRGPI group was worse than that of the low-IRGPI group (logrank test, ) (Figures 3(d) and 3(e)).

3.3. Different Molecular Characteristics of the IRGPI Subgroups

First, we analyzed the differences in the mutational status of the different IRGPI subgroups. In ESCC, the most common mutation type was a missense mutation, followed by nonsense mutation, and was the most common single-nucleotide variant mutation. In both IRGPI groups, TP53, TTN, CSMD3, DNAH5, MUC16, NFE2L2, and PIK3CA were the most commonly mutated genes, showing mutation rates of over 10% in ESCC samples. KMT2D, MUC17, and TGFBR2 mutations were more common in the high-IRGPI group, whereas FLG, ZNF750, and NOTCH1 mutations were more common in the low-IRGPI group (Figure 4).

(a)

(b)

Next, GSEA was performed on different IRGPI subgroups. The top five gene sets enriched in the high-IRGPI group were an interferon-α response, KRAS signaling DN, interferon-gamma response, P53 pathway, and TNFα signaling via NF-κB (Figure 5(a) and Supplementary Table S3). The top five gene sets enriched in the low-IRGPI group were epithelial-mesenchymal transition, hedgehog signaling, angiogenesis, E2F targets, and G2M checkpoint (Figure 5(b)).

(a)

(b)

3.4. Different Immune Characteristics of IRGPI Subgroups

We used the immune cell gene set to analyze immune cell infiltration within the tumor in the TCGA ESCC cohort using the ssGSEA method. The Wilcox rank test was used to distinguish differences in immune cell infiltration between the IRGPI groups (Figure 6(a)).

(a)

(b)

Figure 6

Immune cell infiltration of ESCC samples in the IRGPI subgroups. (a) Wilcoxon rank test was used to inspect the statistical differences of normalized enrichment score (NES) in 28 immune cells between IRGPI subgroups (ns: , , , , and ). The plot shows median (thick lines), quartiles (bottom and top of the boxes), and kernel density estimation (outlines) for each NES distribution. (b) Landscape of the tumor microenvironment in TCGA ESCC subgroups. Clinical features and group information are used as column annotations.

In the high-IRGPI group, activated CD8 T cells, monocytes, neutrophils, and type 17 T helper cells were highly infiltrated. In the low-IRGPI group, memory B cell infiltration was higher than that in the high-IRGPI group. Additionally, we used the cellular landscape to assess the infiltration of immune cells according to the different clinical features and groups of the samples (Figure 6(b)), including the tumor stage and grade, age, gender, and other tumor classification models.

According to the information on the microenvironment and mutation data in TCGA database, tumors were divided into six immune subtypes: C1 (wound healing), C2 (interferon-γ-dominant), C3 (inflammatory), C4 (lymphocyte-depleted), C5 (immunologically quiet), and C6 (TGF-β-dominant) [35]. Compared with the IRGPI subgroups of the TCGA ESCC cohort, there were more C2 types in the high-IRGPI group and more C1 types in the low-IRGPI group (, chi-squared test) (Figure 7(a)). According to the characteristics of the immune cell microenvironment of TCGA squamous cell carcinoma, these cancers could be divided into six immune subtypes, identified as IS 1–6 [36]. IS4 and IS6 subtypes were more prevalent in the high-IRGPI group than in the low-IRGPI group, and there were more IS1 in the low-IRGPI group (, chi-squared test) (Figure 7(b)). According to the molecular characteristics of ESCC, these samples were divided into three subtypes: ESCC1, ESCC2, and ESCC3 [2]. Overall, the high-IRGPI group contained more ESCC2 samples (, chi-squared test) (Figure 7(c)).

(a)

(b)

(c)

3.5. Treatment Strategies Vary between IRGPI Subgroups

We used SubMap analysis in GenePattern to predict the possible IRGPI subgroup response to immunotherapy by comparing with immunotherapy data of melanoma samples. Overall, the results suggested that anti-PD-1 therapy was more likely to be effective in patients in the high-IRGPI group (, Bonferroni corrected ) (Figure 8(a)).

(a)

(b)

Figure 8

Different responses of IRGPI subgroups to immunotherapy and chemotherapy. (a) Compared with published transcriptome data of melanoma immunotherapy using SubMap, the response of IRGPI subgroups to immunotherapy was predicted. The TCGA ESCC high-IRGPI group may be more sensitive to the PD-1 inhibitor. (b) Based on the transcriptome data of GDSC cell lines and IRGPI subgroups, a ridge regression model was constructed to predict the different estimated IC₅₀ values of the drugs. Drugs with a significant difference () between the IRGPI subgroups are shown.

Based on the expression profile of GDSC and TCGA ESCC cohorts, we constructed a ridge regression model to predict the IC₅₀ of drugs from IRGPI subgroups. Drugs such as gefitinib showed a significant difference in IC₅₀ between the IRGPI groups () (Figure 8(b)).

3.6. Verification of the Expression of Immune-Related Genes in ESCC

We verified the expression of genes in the IRGPI in the collected ESCC and adjacent paracancerous tissues. The trend in the qPCR results was the same as that in TCGA ESCC cohorts, confirming the role of these genes in ESCC (Figure 9). The median relative expression of CLDN1 was 1.11 in paracancerous tissues and 1.91 in ESCC tissues, which was significantly increased in tumor tissues (Wilcox test, ). The median relative expression of HCAR3 was 1.87 in paracancerous tissues and 3.04 in ESCC tissues, which was increased in tumor tissues (Wilcox test, ). The median relative expression of FNBP1L was 1.31 in paracancerous tissues and 0.57 in ESCC tissues, which was significantly decreased in tumor tissues (Wilcox test, ). The median relative expression of BRCA2 was 1.04 in paracancerous tissues and 0.46 in ESCC tissues, which was significantly decreased in tumor tissues (Wilcox test, ).

(a)

(b)

(c)

(d)

4. Discussion

The therapeutic effect of ICI therapy results from the interactions of the tumor cells, tumor microenvironment, and immune system. Previous clinical trials showed that ICI therapy can effectively prolong the survival time and reduce treatment-related adverse reactions in patients with advanced ESCC compared with traditional chemotherapy and radiotherapy [5, 6, 37]; however, there are uncertainties, such as drug resistance and pseudoprogression, in some populations [11]. Therefore, it is vital to establish a new and robust method for evaluating and predicting the clinical efficacy of ICI therapy in patients with ESCC.

Expression of immune-related genes in tumor samples impacts the tumor immune microenvironment. Thus, an immune gene signature can effectively help predict the clinical benefit of patients receiving immunotherapy [14]. In this study, differentially expressed immune ESCC genes were identified from an immune gene set. The gene modules closely related to the occurrence of ESCC were screened by WGCNA, and genes related to the prognosis of patients with ESCC were further screened with univariate Cox regression and lasso regression. Taken together, we constructed a prognostic index of immune-related genes in ESCC. The IRGPI was demonstrated to be a stable and robust prognostic indicator using publicly available ESCC data, with a high-IRGPI score associated with a poor prognosis. The predictive potential of the IRGPI was confirmed in a second, validation publicly available ESCC dataset.

The IRGPI is based on the expression of four genes: CLDN1, HCAR3, FNBP1L, and BRCA2. CLDN1 is a membrane protein involved in the formation of tight junctions between cells and regulates the proliferation and metastasis of various tumors [38], including ESCC by inducing autophagy through the AMPK/STAT1/ULK1 signaling pathway [39], and is closely related to lymphocyte reactions in colorectal cancer [40]. HCAR3 is a member of the G protein-coupled receptor superfamily. Previous studies showed that HCAR3 is a potential target for regulating cellular immunity and metabolism [41], with its activation by uracil acid exerting an immunosuppressive effect [42, 43]. Moreover, HCAR3 is essential for the metabolism and proliferation of breast cancer cells [44]. FNBP1L is involved in connecting the cell surface signal to the actin cytoskeleton by interacting with CDC42 and N-WASP. FNBP1L promotes epidermal growth factor-induced cell migration and invasion in epidermal and breast cancer [45, 46]. BRCA2 is a common tumor suppressor gene, and its mutation increases the risk of ESCC [47–49]. In lung and breast cancer cells, long-term induced inactivation of BRCA2 leads to upregulation of interferon-stimulated genes and activation of the cGAS/STING/STAT pathway, confirming that inactivation of BRCA2 triggers cellular innate immune responses [50]. In the calculation of IRGPI, the coefficients of CLDN1 and HCAR3 were positive, whereas those of FNBP1L and BRCA2 were negative. Therefore, IRGPI is positively correlated with the expression of CLDN1 and HCAR3 and negatively correlated with that of FNBP1L and BRCA2.

To explore the differences in the molecular characteristics between IRGPI subgroups, we analyzed their mutational status. transitions are the most common single-nucleotide variant type of ESCC. A high frequency of substitution may be associated with CpG methylation, and the change of germ line methylation can lead to substitution rate variation at the CpG region [51]. The most common mutation type in the two subgroups was missense mutations, followed by nonsense mutations. KMT2D mutation was more common in the high-IRGPI group. It has been reported that KMT2D mutation is the main modulator of ICI in several tumors [52], as KMT2D mutation enhances the immune infiltration and immunogenicity of tumors, thereby making tumors more sensitive to ICI therapy. In turn, ZNF750 mutation was more common in the low-IRGPI group. Studies have shown that ZNF750 is a commonly mutated gene in ESCC, mainly with nonsense mutations. ZNF750 can inhibit epithelial-mesenchymal transition by directly depressing the SNAI1 promoter [53]. In addition, a decrease in ZNF750 levels promotes angiogenesis in ESCC by activating the DANCR/miR-4707-3p/FOXC2 axis [54]. Therefore, analysis of the mutation information of the IRGPI subgroups suggested that patients with high-IRGPI scores are more sensitive to ICI treatment and that those with low-IRGPI values are more prone to the epithelial-mesenchymal transition and angiogenesis phenotype.

GSEA revealed that the enriched gene sets differed between the IRGPI subgroups. The interferon-α response gene set was found to be enriched in the high-IRGPI group. Interferon-α can effectively activate the immune response and reverse the immunosuppressive effect of mesenchymal stem cells [55, 56]. Additionally, the interferon-γ response gene set was enriched in the high-IRGPI group. Interferon-γ maintains immune homeostasis in the tumor microenvironment, limits adaptive and innate immune killing, and, thus, limits the response of patients to ICI treatment [57, 58]. The high-IRGPI group was enriched in the KRAS signaling pathway gene set, which can affect the immune escape of tumor cells [59]. Furthermore, the P53 and NF-κB pathways, which play a role in the tumor immune microenvironment and ICI therapy, were enriched in the high-IRGPI group [60, 61]. The epithelial-mesenchymal transition, angiogenesis, and E2F gene sets were enriched in the low-IRGPI group, further suggesting that this group is prone to tumor invasion, metastasis, and cell cycle. Therefore, the gene enrichment results suggest that high-IRGPI values are associated with tumor immunity, in contrast to the low-IRGPI group.

Next, we analyzed the difference in immune cell infiltration between the IRGPI subgroups. In the high-IRGPI group, activated CD8 T cells, monocytes, neutrophils, and type 17 T helper cells showed higher infiltration, whereas in the low-IRGPI group, memory B cell infiltration was more common. Studies have shown that CD8 T cells are closely related to the expression of PD-L1, suggesting their value for predicting the prognosis of patients and response to ICI treatment [62–64]. The increase in the neutrophil count in tumors is often closely related to poor prognosis, explaining the poor prognosis of patients in the high-IRGPI group [65]. Interleukin-17, produced by type 17 T helper cells, stimulates tumor and stroma cells to produce tumor-promoting factors, whereas interleukin-8, produced by type 17 T helper cells, recruits neutrophils [66]. Memory B cell infiltration was higher in the low-IRGPI group. A previous study showed that the new subgroup of memory B cells can promote angiogenesis [67]. Therefore, analysis of immune cell infiltration in the IRGPI subgroups showed that under the effect of immune-infiltrating cells, high-IRGPI values were indicative of adverse prognosis and predicted the ICI response, whereas low-IRGPI indicated angiogenesis.

Combined with other immune subtypes, we can determine the immune status of IRGPI subgroups. Compared with pan-cancer immune subtypes, there were more C2-type patients in the high-IRGPI group and more C1-type patients in the low-IRGPI group [35]. The C2 type is characterized by an IFN-γ response and high CD8 T cell markers and lymphocyte infiltration rate, indicating a better immune response but poor prognosis. The C1 type indicates more angiogenic gene expression. In squamous cell carcinoma immune subtypes, IS4 and IS6 were more prevalent in the high-IRGPI group than in the low-IRGPI group [36]. IS4 showed the highest T cell expression and IFN-γ response, with a good immune activation phenotype, whereas IS1 showed an immunosuppressive phenotype. The high-IRGPI group tended to undergo ESCC2 classification, which involves higher leukocyte infiltration [2]. Overall, the immunophenotype of the high-IRGPI group was more active than that of the low-IRGPI group.

5. Conclusions

Traditional chemotherapy and radiotherapy are effective in few patients with ESCC. Hence, for these patients, ICI therapy may be beneficial for prolonging their survival time and reducing the incidence of treatment-related adverse reactions. Taken together, this study may fill the gap in the need for new biomarkers and the proposed IRGPI may be used as a biomarker to evaluate ESCC prognosis and the response to ICI therapy.

Data Availability

The datasets generated and analyzed during the current study are available in The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov) and the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo).

Disclosure

The previous version of this research was published as a conference abstract in Annals of Oncology and was further published after adding bioinformatics analyses and verification of clinical samples.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

Qinghua Ji and Yingying Cai contributed equally to this work.

Supplementary Materials

Supplementary Figure S1: Gene ontology (GO) enrichment analysis of differentially expressed immune-related genes in the TCGA ESCC cohort. Supplementary Figure S2: WGCNA method used to identify the hub genes. (A) Determination of soft threshold in WGCNA. (B) WGCNA analysis of differentially expressed immune-related genes in ESCC and the gene modules. (C) Gene expression clustering tree and coexpression topological heat map. (D) Relationship between the gene modules and features of ESCC. Supplementary Figure S3: differentially expressed immune-related genes analyzed using lasso Cox regression to determine genes included in the prognostic model. (A) Lasso Cox regression fitting process. (B) Tenfold crossvalidation was used to determine the λ values. The lambda value in lasso Cox reflects the degree of regularization. The larger the lambda value, the fewer effective variables are screened. Tenfold crossvalidation was used to determine the lambda values; the minimum crossvalidation error is set to the best lambda value. Supplementary Table S1: primers used for qPCR. Supplementary Table S2: GO terms and KEGG pathway enrichment analyses of differentially expressed immune-related genes of ESCC. Supplementary Table S3: GSEA enrichment analysis between IRGPI subgroups. (Supplementary Materials)

References

E. C. Smyth, J. Lagergren, R. C. Fitzgerald et al., “Oesophageal cancer,” Nature Reviews Disease Primers, vol. 3, no. 1, article 17048, 2017.
View at: Publisher Site | Google Scholar
A. Ally, M. Balasundaram, R. Carlsen et al., “Integrated genomic characterization of oesophageal carcinoma,” Nature, vol. 541, no. 7636, pp. 169–175, 2017.
View at: Publisher Site | Google Scholar
C. C. Abnet, M. Arnold, and W. Q. Wei, “Epidemiology of esophageal squamous cell carcinoma,” Gastroenterology, vol. 154, no. 2, pp. 360–373, 2018.
View at: Publisher Site | Google Scholar
G O C Collaborators, “The global, regional, and national burden of oesophageal cancer and its attributable risk factors in 195 countries and territories, 1990-2017: a systematic analysis for the Global Burden of Disease Study 2017,” The Lancet Gastroenterology & Hepatology, vol. 5, no. 6, pp. 582–597, 2020.
View at: Google Scholar
T. Kojima, M. A. Shah, K. Muro et al., “Randomized phase III KEYNOTE-181 study of pembrolizumab versus chemotherapy in advanced esophageal cancer,” Journal of Clinical Oncology, vol. 38, no. 35, pp. 4138–4148, 2020.
View at: Publisher Site | Google Scholar
K. Kato, B. C. Cho, M. Takahashi et al., “Nivolumab versus chemotherapy in patients with advanced oesophageal squamous cell carcinoma refractory or intolerant to previous chemotherapy (ATTRACTION-3): a multicentre, randomised, open-label, phase 3 trial,” Lancet Oncology, vol. 20, no. 11, pp. 1506–1517, 2019.
View at: Publisher Site | Google Scholar
S. Bagchi, R. Yuan, and E. G. Engleman, “Immune checkpoint inhibitors for the treatment of cancer: clinical impact and mechanisms of response and resistance,” Annual Review of Pathology: Mechanisms of Disease, vol. 16, no. 1, pp. 223–249, 2021.
View at: Publisher Site | Google Scholar
X. Liu, G. D. Hogg, and D. G. DeNardo, “Rethinking immune checkpoint blockade: 'beyond the T cell',” Journal for Immuno Therapy of Cancer, vol. 9, no. 1, 2021.
View at: Publisher Site | Google Scholar
J. A. Marin-Acevedo, E. O. Kimbrough, and Y. Lou, “Next generation of immune checkpoint inhibitors and beyond,” Journal of Hematology & Oncology, vol. 14, no. 1, p. 45, 2021.
View at: Publisher Site | Google Scholar
A. Ribas and J. D. Wolchok, “Cancer immunotherapy using checkpoint blockade,” Science, vol. 359, no. 6382, pp. 1350–1355, 2018.
View at: Publisher Site | Google Scholar
M. Nishino, N. H. Ramaiya, H. Hatabu, and F. S. Hodi, “Monitoring immune-checkpoint blockade: response evaluation and biomarker development,” Nature Reviews Clinical Oncology, vol. 14, no. 11, pp. 655–668, 2017.
View at: Publisher Site | Google Scholar
G. Manson, J. Norwood, A. Marabelle, H. Kohrt, and R. Houot, “Biomarkers associated with checkpoint inhibitors,” Annals of Oncology, vol. 27, no. 7, pp. 1199–1206, 2016.
View at: Publisher Site | Google Scholar
M. Yarchoan, A. Hopkins, and E. M. Jaffee, “Tumor mutational burden and response rate to PD-1 inhibition,” New England Journal of Medicine, vol. 377, no. 25, pp. 2500-2501, 2017.
View at: Publisher Site | Google Scholar
G. T. Gibney, L. M. Weiner, and M. B. Atkins, “Predictive biomarkers for checkpoint inhibitor-based immunotherapy,” Lancet Oncology, vol. 17, no. 12, pp. e542–e551, 2016.
View at: Publisher Site | Google Scholar
R. R. Ji, S. D. Chasalow, L. Wang et al., “An immune-active tumor microenvironment favors clinical response to ipilimumab,” Cancer Immunology, Immunotherapy, vol. 61, no. 7, pp. 1019–1031, 2012.
View at: Publisher Site | Google Scholar
Y. Dai, W. Qiang, K. Lin, Y. Gui, X. Lan, and D. Wang, “An immune-related gene signature for predicting survival and immunotherapy efficacy in hepatocellular carcinoma,” Cancer Immunology, Immunotherapy, vol. 70, no. 4, pp. 967–979, 2021.
View at: Publisher Site | Google Scholar
J. Wu, L. Li, H. Zhang et al., “A risk model developed based on tumor microenvironment predicts overall survival and associates with tumor immunity of patients with lung adenocarcinoma,” Oncogene, vol. 40, no. 26, pp. 4413–4424, 2021.
View at: Publisher Site | Google Scholar
J. Li, Z. Chen, L. Tian et al., “LncRNA profile study reveals a three-lncRNA signature associated with the survival of patients with oesophageal squamous cell carcinoma,” Gut, vol. 63, no. 11, pp. 1700–1710, 2014.
View at: Publisher Site | Google Scholar
S. Bhattacharya, P. Dunn, C. G. Thomas et al., “ImmPort, toward repurposing of open access immunological assay data for translational and clinical research,” Scientific Data, vol. 5, no. 1, article 180015, 2018.
View at: Publisher Site | Google Scholar
J. Kelley, B. de Bono, and J. Trowsdale, “IRIS: a database surveying known human immune system genes,” Genomics, vol. 85, no. 4, pp. 503–511, 2005.
View at: Publisher Site | Google Scholar
K. Breuer, A. K. Foroushani, M. R. Laird et al., “InnateDB: systems biology of innate immunity and beyond—recent updates and continuing curation,” Nucleic Acids Research, vol. 41, no. D1, pp. D1228–D1233, 2013.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, p. e47, 2015.
View at: Publisher Site | Google Scholar
G. Yu, L. G. Wang, Y. Han, and Q. Y. He, “Cluster Profiler: an R package for comparing biological themes among gene clusters,” Omics-A Journal of Integrative Biology, vol. 16, no. 5, pp. 284–287, 2012.
View at: Google Scholar
P. Langfelder and S. Horvath, “FastRFunctions for robust correlations and hierarchical clustering,” Journal of Statistical Software, vol. 46, no. 11, 2012.
View at: Publisher Site | Google Scholar
J. Friedman, T. Hastie, and R. Tibshirani, “Regularization paths for generalized linear models via coordinate descent,” Journal of Statistical Software, vol. 33, no. 1, pp. 1–22, 2010.
View at: Publisher Site | Google Scholar
E. L. Kaplan and P. Meier, “Nonparametric estimation from incomplete observations,” Journal of The American Statistical Association, vol. 53, no. 282, pp. 457–481, 1958.
View at: Publisher Site | Google Scholar
J. M. Bland and D. G. Altman, “The logrank test,” BMJ, vol. 328, no. 7447, p. 1073, 2004.
View at: Publisher Site | Google Scholar
A. Mayakonda, D. C. Lin, Y. Assenov, C. Plass, and H. P. Koeffler, “Maftools: efficient and comprehensive analysis of somatic variants in cancer,” Genome Research, vol. 28, no. 11, pp. 1747–1756, 2018.
View at: Publisher Site | Google Scholar
The ULTRA-DD Consortium, H. Fang, H. de Wolf et al., “A genetics-led approach defines the drug target landscape of 30 immune- related traits,” Nature Genetics, vol. 51, no. 7, pp. 1082–1091, 2019.
View at: Publisher Site | Google Scholar
S. Hanzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, 2013.
View at: Publisher Site | Google Scholar
Z. Gu, R. Eils, and M. Schlesner, “Complex heatmaps reveal patterns and correlations in multidimensional genomic data,” Bioinformatics, vol. 32, no. 18, pp. 2847–2849, 2016.
View at: Publisher Site | Google Scholar
Y. Hoshida, J. P. Brunet, P. Tamayo, T. R. Golub, and J. P. Mesirov, “Subclass mapping: identifying common subtypes in independent disease data sets,” PLoS One, vol. 2, no. 11, p. e1195, 2007.
View at: Publisher Site | Google Scholar
W. Roh, P. L. Chen, A. Reuben et al., “Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance,” Science Translational Medicine, vol. 9, no. 379, article eaah3560, 2017.
View at: Publisher Site | Google Scholar
P. Geeleher, N. Cox, and R. S. Huang, “pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels,” PLoS One, vol. 9, no. 9, article e107468, 2014.
View at: Publisher Site | Google Scholar
V. Thorsson, D. L. Gibbs, S. D. Brown et al., “The immune landscape of cancer,” Immunity, vol. 48, no. 4, pp. 812–830.e14, 2018.
View at: Publisher Site | Google Scholar
B. Li, Y. Cui, D. K. Nambiar, J. B. Sunwoo, and R. Li, “The immune subtypes and landscape of squamous cell carcinoma,” Clinical Cancer Research, vol. 25, no. 12, pp. 3528–3537, 2019.
View at: Google Scholar
K. Kato, M. A. Shah, P. Enzinger et al., “KEYNOTE-590: phase III study of first-line chemotherapy with or without pembrolizumab for advanced esophageal cancer,” Future Oncology, vol. 15, no. 10, pp. 1057–1066, 2019.
View at: Publisher Site | Google Scholar
A. A. Bhat, N. Syed, L. Therachiyil et al., “Claudin-1, a double-edged sword in cancer,” International Journal of Molecular Sciences, vol. 21, no. 2, 2020.
View at: Google Scholar
J. Wu, F. Gao, T. Xu et al., “CLDN1 induces autophagy to promote proliferation and metastasis of esophageal squamous carcinoma through AMPK/STAT1/ULK1 signaling,” Journal of Cellular Physiology, vol. 235, no. 3, pp. 2245–2259, 2020.
View at: Publisher Site | Google Scholar
D. Süren, M. Yıldırım, V. Kaya et al., “Loss of tight junction proteins (Claudin 1, 4, and 7) correlates with aggressive behavior in colorectal carcinoma,” Medical Science Monitor: international medical journal of experimental and clinical research, vol. 20, pp. 1255–1262, 2014.
View at: Google Scholar
N. J. Kapolka and D. G. Isom, “HCAR3: an underexplored metabolite sensor,” Nature Reviews Drug Discovery, vol. 19, no. 11, p. 745, 2020.
View at: Publisher Site | Google Scholar
N. J. Kapolka, G. J. Taghon, J. B. Rowe et al., “DCyFIR: a high-throughput CRISPR platform for multiplexed G protein-coupled receptor profiling and ligand discovery,” Proceedings of the National Academy of Sciences, vol. 117, no. 23, pp. 13117–13126, 2020.
View at: Google Scholar
I. Cervenka, L. Z. Agudelo, and J. L. Ruas, “Kynurenines: Tryptophan’s metabolites in exercise, inflammation, and mental health,” Science, vol. 357, no. 6349, 2017.
View at: Google Scholar
C. Staubert, O. J. Broom, and A. Nordstrom, “Hydroxycarboxylic acid receptors are essential for breast cancer cells to control their lipid/fatty acid metabolism,” Oncotarget, vol. 6, no. 23, pp. 19706–19720, 2015.
View at: Publisher Site | Google Scholar
J. Hu, A. Mukhopadhyay, and A. W. Craig, “Transducer of Cdc42-dependent Actin Assembly Promotes Epidermal Growth Factor- induced Cell Motility and Invasiveness,” Journal of Biological Chemistry, vol. 286, no. 3, pp. 2261–2272, 2011.
View at: Publisher Site | Google Scholar
H. Chander, P. Truesdell, J. Meens, and A. W. Craig, “Transducer of Cdc 42-dependent actin assembly promotes breast cancer invasion and metastasis,” Oncogene, vol. 32, no. 25, pp. 3080–3090, 2013.
View at: Google Scholar
J. M. Ko, L. Ning, X. K. Zhao et al., “BRCA2loss‐of‐function germline mutations are associated with esophageal squamous cell carcinoma risk in Chinese,” International Journal of Cancer, vol. 146, no. 4, pp. 1042–1051, 2020.
View at: Publisher Site | Google Scholar
M. R. Akbari, R. Malekzadeh, D. Nasrollahzadeh et al., “Germline BRCA2 mutations and the risk of esophageal squamous cell carcinoma,” Oncogene, vol. 27, no. 9, pp. 1290–1296, 2008.
View at: Publisher Site | Google Scholar
N. Hu, C. Wang, X. Y. Han et al., “Evaluation of BRCA2 in the genetic susceptibility of familial esophageal cancer,” Oncogene, vol. 23, no. 3, pp. 852–858, 2004.
View at: Publisher Site | Google Scholar
T. Reisländer, E. P. Lombardi, F. J. Groelly et al., “BRCA2 abrogation triggers innate immune responses potentiated by treatment with PARP inhibitors,” Nature Communications, vol. 10, no. 1, p. 3143, 2019.
View at: Publisher Site | Google Scholar
C. F. Mugal and H. Ellegren, “Substitution rate variation at human CpG sites correlates with non-CpG divergence, methylation level and GC content,” Genome Biology, vol. 12, no. 6, p. R58, 2011.
View at: Publisher Site | Google Scholar
G. Wang, R. D. Chow, L. Zhu et al., “CRISPR-GEMM pooled mutagenic screening identifies KMT2D as a major modulator of immune checkpoint blockade,” Cancer Discovery, vol. 10, no. 12, pp. 1912–1933, 2020.
View at: Publisher Site | Google Scholar
P. Kong, E. Xu, Y. Bi et al., “Novel ESCC-related gene ZNF750 as potential prognostic biomarker and inhibits epithelial-mesenchymal transition through directly depressing SNAI1 promoter in ESCC,” Theranostics, vol. 10, no. 4, pp. 1798–1813, 2020.
View at: Publisher Site | Google Scholar
Y. Bi, S. Guo, X. Xu et al., “Decreased _ZNF750_ promotes angiogenesis in a paracrine manner via activating DANCR/miR-4707-3p/FOXC2 axis in esophageal squamous cell carcinoma,” Cell Death & Disease, vol. 11, no. 4, p. 296, 2020.
View at: Publisher Site | Google Scholar
E. C. Borden, “Interferons α and β in cancer: therapeutic opportunities from new insights,” Nature Reviews Drug Discovery, vol. 18, no. 3, pp. 219–234, 2019.
View at: Publisher Site | Google Scholar
P. Shou, Q. Chen, J. Jiang et al., “Type I interferons exert anti-tumor effect via reversing immunosuppression mediated by mesenchymal stromal cells,” Oncogene, vol. 35, no. 46, pp. 5953–5962, 2016.
View at: Publisher Site | Google Scholar
J. L. Benci, L. R. Johnson, R. Choa et al., “Opposing functions of interferon coordinate adaptive and innate immune responses to cancer immune checkpoint blockade,” Cell, vol. 178, no. 4, pp. 933–948.e14, 2019.
View at: Publisher Site | Google Scholar
C. J. Nirschl, M. Suárez-Fariñas, B. Izar et al., “IFNγ-Dependent Tissue-Immune Homeostasis Is Co-opted in the Tumor Microenvironment,” Cell, vol. 170, no. 1, pp. 127–141.e15, 2017.
View at: Publisher Site | Google Scholar
S. Hamarsheh, O. Groß, T. Brummer, and R. Zeiser, “Immune modulatory effects of oncogenic KRAS in cancer,” Nature Communications, vol. 11, no. 1, p. 5439, 2020.
View at: Publisher Site | Google Scholar
C. Munoz-Fontela, A. Mandinova, S. A. Aaronson, and S. W. Lee, “Emerging roles of p53 and other tumour-suppressor genes in immune regulation,” Nature Reviews Immunology, vol. 16, no. 12, pp. 741–750, 2016.
View at: Publisher Site | Google Scholar
C. S. Garris, S. P. Arlauckas, R. H. Kohler et al., “Successful anti-PD-1 cancer immunotherapy requires T cell-dendritic cell crosstalk involving the cytokines IFN-γ and IL-12,” Immunity, vol. 49, no. 6, pp. 1148–1161.e7, 2018.
View at: Publisher Site | Google Scholar
W. H. Fridman, L. Zitvogel, C. Sautès–Fridman, and G. Kroemer, “The immune contexture in cancer prognosis and treatment,” Nature Reviews Clinical Oncology, vol. 14, no. 12, pp. 717–734, 2017.
View at: Publisher Site | Google Scholar
D. Bruni, H. K. Angell, and J. Galon, “The immune contexture and Immunoscore in cancer prognosis and therapeutic efficacy,” Nature Reviews Cancer, vol. 20, no. 11, pp. 662–680, 2020.
View at: Publisher Site | Google Scholar
F. Petitprez, N. Fossati, Y. Vano et al., “PD-L1 expression and CD8(+) T-cell infiltrate are associated with clinical progression in patients with node-positive prostate cancer,” European Urology Focus, vol. 5, no. 2, pp. 192–196, 2019.
View at: Google Scholar
A. J. Templeton, M. G. McNamara, B. Seruga et al., “Prognostic role of neutrophil-to-lymphocyte ratio in solid tumors: a systematic review and meta-analysis,” JNCI: Journal of the National Cancer Institute, vol. 106, no. 6, p. u124, 2014.
View at: Google Scholar
D. R. Powell and A. Huttenlocher, “Neutrophils in the tumor microenvironment,” Trends in Immunology, vol. 37, no. 1, pp. 41–52, 2016.
View at: Publisher Site | Google Scholar
W. van de Veen, A. Globinska, K. Jansen et al., “A novel proangiogenic B cell subset is increased in cancer and chronic inflammation,” Science Advances, vol. 6, no. 20, article z3559, 2020.
View at: Google Scholar

Copyright

Copyright © 2021 Qinghua Ji et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

408

Downloads

994

Citations