Comprehensive Analysis of Tumor-Infiltrating Immune Cells and Relevant Therapeutic Strategy in Esophageal Cancer

A growing body of evidence has indicated that behaviors of cancers are defined by not only intrinsic activities of tumor cells but also tumor-infiltrating immune cells (TIICs) in the tumor microenvironment. However, it still lacks a well-structured and comprehensive analysis of TIICs and its therapeutic value in esophageal cancer (EC). The proportions of 22 TIICs were evaluated between 150 normal tissues and 141 tumor tissues of EC by the CIBERSORT algorithm. Besides, correlation analyses between proportions of TIICs and clinicopathological characters, including age, gender, histologic grade, tumor location, histologic type, LRP1B mutation, TP53 mutation, tumor stage, lymph node stage, and TNM stage, were conducted. We constructed a risk score model to improve prognostic capacity with 5 TIICs by least absolute shrinkage and selection operator (lasso) regression analysis. The risk score = −1.86∗plasma + 2.56∗T cell follicular helper − 1.37∗monocytes − 3.64∗activated dendritic cells − 2.24∗resting mast cells (immune cells in the risk model mean the proportions of immune cell infiltration in EC). Patients in the high-risk group had significantly worse overall survival than these in the low-risk group (HR: 2.146, 95% CI: 1.243-3.705, p = 0.0061). Finally, we identified Semustine and Sirolimus as two candidate compounds for the treatment of EC based on CMap analysis. In conclusion, the proportions of TIICs may be important to the progression, prognosis, and treatment of EC.


Introduction
Esophageal cancer (EC) is one of the most common digestive tract cancers; the morbidity and mortality of which is ranked 9th and 6th in all malignant tumors. It is estimated that, in 2018, 508,585 people worldwide would die of EC [1]. Thus, it becomes extremely urgent to find effective treatments. In most cases, for cancers in early stages, the general excision of cancers is chosen as the surgical treatment. However, due to the lack of early symptoms, most affected EC patients lose the optimal opportunity of surgery [2]. Although in recent years, new adjuvant therapy, chemotherapy, and precise radiation therapy for EC have made some progresses, unfortunately, the overall efficacy is still not ideal, and the 5-year survival rate is around 30% [3]. Against this backdrop, identifying other effective treatments is significant for improving the survival rate of EC patients.
Increasing evidence indicates that behaviors of cancers are defined by not only the intrinsic activities of tumor cells but also the tumor-infiltrating immune cells (TIICs) in the tumor microenvironment. TIICs are the heterogeneous immune populations existing in the tumor tissues and which play a key role in host antigen-specific tumor immune response. The cells including T cells, B cells, natural killer (NK) cells, macrophages, dendritic cells (DC), and polymorphonuclear leukocytes all belong to TIICs [4]. In this regard, it is reported in hepatocellular carcinoma (HCC) that the functional interaction of tumor-infiltrating T cells and B cells was contributed to the prognosis of HCC patients through immune activation [5]. In breast cancer, the findings demonstrated that the increased fractions of regulatory T cells and M0 macrophages were linked to a lower pathological complete response rate, shorter disease-free survival (DFS), and worse overall survival (OS) [6]. For EC, the reports similarly have linked the presence of TIICs with its treatment response and outcome. In these stage II+III patients, the densities of NK cells and macrophages found a significant relation with patients' postoperative prognoses [7]. For esophageal squamous cell carcinoma (ESCC), the expression of PD-L1 was positively associated with TIIC density and that was also correlated with worse prognosis [8,9]. Another evidence indicated that regulatory T cell (Treg) infiltrate presented in the tumor had an association with the pathological response and exhibited a favorable value in predicting cancer-specific survival [10]. Thus, the full analysis of the types and the range of TIICs is a promising strategy to make a huge change for the treatment of EC.
In colorectal cancer, it has more recently been established that the type, density, and location of TIICs within tumor samples are more powerful to predict patient survival than the histopathological methods currently used for its tumor stage [11,12]. However, for EC, it still lacks a wellstructured and comprehensive analysis of TIICs and its therapeutic value. Due to the limitations of methods and techniques, previous studies thus were focused merely on the finite areas of immune response. Recently, with the development of a novel metagene approach of CIBERSORT [13], it is possible to computationally dissect the density of TIICs and then to predict the clinical value in EC. In this study, 20 differential fractions of TIICs were identified in normal tissues and tumor tissues by CIBERSORT. Based on the least absolute shrinkage and selection operator (lasso) regression model, it was revealed that TIICs exhibited potential effect in EC prognosis. Additionally, following the gene ontology (GO) analysis for illuminating the differential enrichment signals between the low-and high-risk groups, we also found candidate compounds for the treatment of EC.
2.2. Assessment of Immune Cell Infiltration. The CIBERSORT (RRID:SCR_016955) algorithm was applied to evaluate the proportions of TIICs in tissues. CIBERSORT was a method that was designed and robustly validated to identify 22 human immune cell phenotypes, outperforming other methods in the matter of noise and unknown mixed content [13]. 22 human immune cell phenotypes were analyzed in the study, including seven T cell types (T cell CD8, naïve T cell CD4, resting T cell CD4 memory, activated T cell CD4 memory, T cell follicular helper, regulatory T cells (Tregs), and T cell gamma delta); naïve and memory B cells; plasma cells; resting and activated NK cells; monocytes; macrophages M0, M1, and M2; resting and activated dendritic cells; resting and activated mast cells; eosinophils; and neutrophils. The CIBERSORT p value and root mean squared error (RMSE) were calculated for each tissues. Only samples with a CIBER-SORT p value < 0.05 were enrolled. Eventually, there were 144 normal samples from GTEx and 6 normal samples and 141 tumor samples from TCGA were eligible in the study (Supplementary Figure 1).

Connectivity Map (CMap)
Analysis. CMap (https:// portals.broadinstitute.org/cmap/, RRID:SCR_015674) was used to find connections between drugs and genes [18]. Up tag file and down tag file were uploaded into the quick query of the tool web. The value of connectivity score fluctuated between -1 and 1, and a high negative connectivity score manifested that the drug reversed the expression of the query signature.

Statistical Analyses.
For each immune cell phenotypes, we calculated the quartile, median, and third quartile of the normal and tumor groups. The Wilcoxon signed-rank test was used to compare the different TIICs. Correlation analysis was performed by package "corrplot"of R. Correlations between clinicopathological characters and proportions of TIICs were realized by the Wilcoxon signed-rank test and visualized by GraphPad Prism 8.0.1 (RRID:SCR_002798). The Lasso regression model was built by package "glmnet" and "survival" of R. Univariate and multivariate cox analyses and forestplot were completed by "survival" and "forestplot" of R. GO analysis between the high-and low-risk groups was conducted by "clusterProfiler," "http://org.Hs.eg.db/," "enrichplot," and "ggplot2." Gene set enrichment analysis (GSEA) was performed to get the core genes of top sets (RRID:SCR_003199) [19,20]. All tests were two-tailed p value, and p value < 0.05 was considered significant. Survival was evaluated with the hazard ratio (95% confidence intervals).

Correlations between Clinicopathological Characters and
Proportions of TIICs in EC. For basic information of EC patients, the proportions of naïve B cell and regulatory T cells (Tregs) were higher in older patients (age > 65), while the proportion of activated dendritic cells was significantly reduced in older male patients (Figures 2(a) and 2(b)). The percentage of monocytes was slightly increased in male patients ( Figure 2(b)). The higher the histologic grade was, the higher proportions of naïve B cell and regulatory T cells (Tregs) were found in EC (Figure 2 (Figure 3(c)). Then, we divided patients into two groups: the high-risk group (risk score > median value of all patients) and the low-risk group (risk score ≤ median value). The survivorship curve (Figure 3(d)) indicated that patients in the high-risk group had significantly worse overall survival than that in low-risk group (HR: 2.146, 95% CI: 1.243-3.705, p = 0:0061).  great diagnostic performance of the risk model. As a validation of GSE19417 in an indirect manner, the risk score of poorly differentiated group was higher than that of the well-differentiated group (p = 0:025). Besides, the risk score of patients with more than 5 positive lymph nodes was higher than that in patients with only 1 positive lymph node (p = 0:033) (Supplementary Figure 3). It is reported that the increasing number of positive nodes is related to poor outcome of EC [21]. Furthermore, from the results of univariate (HR: 3.150, 95% CI: 1.435-6.915, p = 0:004, Figure 4(a)) and multivariate (HR: 2.342, 95% CI: 1.012-5.416, p = 0:047, Figure 4(b)) cox analyses, the risk score was the independent risk factor for EC.

Identification of Two Candidate
Compounds for the Treatment of EC Based on CMap Analysis. GO analysis was conducted to illuminate the differential enrichment signals between the low-and high-risk groups stratified by immune cell infiltration of EC. The top ten enriched signals were shown in the matter of biological process (BP), cellular    Tumor stage   T3-T4  T1-T2  T3-T4  T1-T2  T3-T4  T1-T2  T3-T4  T1-  . And the terms "mRNA processing," "chromosomal region," and "chromatin binding" were the most enriched signals of BP, CC, and MF, respectively ( Figure 5(a)). Then, GSEA was performed to identify the core genes. Figure 5(b) showed the most enriched sets mRNA binding, mRNA metabolic process, and M phase of mitotic cell cycle in the high-risk group, while Figure 5(c) displayed the most enriched gene sets receptor-mediated endocytosis, cell substrate adhesion, and integrin binding in the low-risk group. The core genes of high-risk group obtained from mRNA binding, mRNA metabolic process, and M phase of the mitotic cell cycle were visualized in Figure 6(a), and the core genes of low-risk group from gene sets, receptor mediated endocytosis, cell substrate adhesion, and integrin binding, were revealed in Figure 6(b). Then, 11 upregulation core genes (LSM3, LSM5, PABPC1, PABPC3, PABPC4, EIF4A3, SMC3, ESPL1, CDC25B, CCNA2, and AURKA) and 1 downregulation core gene (THBS3) were uploaded into CMap as upregulated tags and downregulated tags, respectively. Two compounds (Semustine and Sirolimus) might be the candidate compounds for the treatment of EC.

Discussion
To date, malignant tumors remain the main killer of human health [1]. The focus of the research in the past with regard to the tumor development is often on the tumor cells themselves. In 1989, Paget first proposed the theory of "seed and soil," which has been widely recognized and extended since its introduction [22]. The theory holds that the occurrence and development of tumors not only depend on the change in tumor cytogenetics and epigenetics but also link to the tumor microenvironment as a "fertile soil" for the growth and breeding of malignant seeds [22,23]. In normal tissues, the microenvironment is an important barrier for defense against tumors. However, the inhibitory function in tumor cells would be reduced by recruiting integrated fibroblasts,   Disease Markers regulating antitumor immune cells, secreting immunosuppressive molecules, and thus transforming into tumor microenvironment suitable for tumor growth [24]. It is just the interaction and coevolution of tumor cells and their microenvironment that promote the production of tumors. In recent years, increasing attention has been given to the tumor microenvironment in terms of the potential response to therapy. In EC, it is also reported that tumor microenvironment has a close relationship with tumor growth. The cells, including NK cells, T cells, B cells, plasma cells, monocytes, macrophages, dendritic cells, mast cells, eosinophils, and neutrophils, have all been reported to have exhibited significant association with EC treatment [7,10,[25][26][27][28][29][30][31][32]. The elevation or reduction of these different types of TIICs in tumor microenvironment is commonly associated with poor prognosis of EC. Nevertheless, prior studies only focused on the density and clinical value of one or two TIICs in EC, lacking systematic analysis of all TIICs with a different density. In this study, following the operation of CIBERSORT algorithm, we found that 20 immune-infiltrating cells had statistic difference simultaneously between normal and tumor tissues. Of note, the proportions of B cell memory, monocytes, and resting mast cells were significantly decreased in EC, while the proportions of activated T cell CD4 memory, Tregs, macrophage M0, macrophage M1, resting dendritic cells, and activated dendritic cells were significantly increased.
Tregs are a subpopulation of CD4+ T lymphocytes. It could inhibit the antitumor immunity and assist in tumor immune escape through secreting immunosuppressive fac-tors, directly killing or inhibiting effector cell proliferation, as well as affecting T cell activation [33,34]. The clinical research also points that the high-density infiltration of Tregs often indicates the poor clinical outcomes. In severely immunodeficient mice, the repeated reinfusion of patient-derived CD3+CD25 T cells could prevent tumor growth. For the negative regulation of Tregs against tumor immune function, in line with this, it was found that the reinfusion of Tregs could reverse such protection [35,36]. Similarly, Nishikawa et al. found that leukocytes isolated from patients with melanoma and ovarian cancer had responded to selective tumor antigens only after removal of Tregs [37]. In this study, by analyzing the correlations between clinicopathological characters and TIICs, we found that Treg density infiltration in EC patients had a close association with age, histologic grade, tumor location, histologic type, LRP1B mutation, tumor stage, lymph node stage, and TNM stage. In fact, it has been acknowledged that Tregs increased according to the EC progression [38]. In ESCC, IL32 expression and Treg infiltration were found to play an important synergistic role in tumor growth and invasion; the combination of which was one of poor independent factors [39]. In operative specimens, generally, the frequency of local Tregs negatively correlated with the pathological response and overall patient survival [40]. However, up to now, there are still no relevant reports on Tregs and LRP1B, and it remains to be studied whether the poor prognosis of patients in the high age group can be considered from the perspective of Treg infiltrating density.

Disease Markers
Additionally, from the results in Figure 2, it was shown that macrophages M0 and M1 differed greatly in different tumor stages and histologic types, respectively. Unlike Tregs, infiltrated macrophages are developed from a precursor of bone marrow monocytes that are activated into different subtypes by stimulation signals in different microenvironments [41,42]. Studies have shown that macrophages induced immune disability through promoting angiogenesis, inducing tumor metastasis, provoking chemotherapy tolerance, and interacting with the other immune cells or raising other immune inhibitory cells [43]. Macrophage is one of the important factors in tumor microenvironments that cause the poor prognosis of patients. But for the moment, it is relatively frequent on the research about the relationship between macrophages and ESCC. Compared with adenocarcinoma, there is no basic data to support whether the infiltration of macrophages in ESCC is relatively high. In addition to macrophages, dendritic cells are also associated with a variety of clinicopathological features in EC. Dendritic cell is one of the most potent antigen-presenting cells [44]. It was realized that the impaired immune function and the decreased number of dendritic cells were the significant causes for the pathogenesis and progression of EC [45]. In EC, a research once has pointed that LRP1B possessed the recurrent copynumber variant character, significantly promoting cancer cell proliferation, migration, and invasion [46]. A deep research of the relationship between dendritic cells and LRP1B gene is a great necessity since there are no such reports until now.
Considering the low survival rate of EC patients, subsequently, univariate cox regression analysis was performed in exploring the prognostic value of these immune cell infiltration. It was found that one proportion of immune cell infiltration cannot predict the outcome. According to the lasso regression model, it was shown that the joint detection of plasma, T cell follicular helper, monocytes, activated dendritic cells, and resting mast cell infiltration exhibited great performance for its diagnostic. In EC, as indicated above, for instance, Svensson et al. guessed that the antitumoral effects on time to recurrence (TTR) and OS were largely dependent on a functional interplay between T and B lymphocytes or plasma cells [47]. In the study by Lu et al., the higher number of CD1a dendritic cells had a correlation with significantly improved OS of patients with ESCC [48]. Nonetheless, due to the more complexity of TIICs, the full consideration of the link between TIICs and tumor prognosis is still in the blank stage. In this study, we have not only addressed the issue but also found that the risk score was an independent risk factor for EC. The risk factor presented a great link with EC prognosis including age, sex, location of tumor, histologic type, tumor stage, and nodal invasion. Overall, in assessing the EC prognosis, it is necessary to fully evaluate the above factors. In this regard, while multiple reports in   [49], many others have found no exact assessment for considering all involved factors. Thus, in such circumstances, the acquisition of an independent risk factor would favor tumor prognosis assessment. In such later cases, sometimes, for simple and convenient assessment, the risk score could be chosen as an independent risk factor. Additionally, the study found that mRNA binding, the mRNA metabolic process, and the M phase of mitotic cell cycle were enriched in the high-risk group, while receptormediated endocytosis, cell substrate adhesion, and integrin binding were enriched in the low-risk group. Through the assessment of these differential enrichment signals, 11 upregulation core genes (LSM3, LSM5, PABPC1, PABPC3, PABPC4, EIF4A3, SMC3, ESPL1, CDC25B, CCNA2, and AURKA) and 1 downregulation core gene (THBS3) were identified. LSM3 was found downregulated in cervical cancer, correlating with progression free survival [50,51]. In colorectal cancer, LSM3 also was significantly associated with lymphatic metastasis [52], but in EC, there is still no report. Similarly, there is just a report in breast cancer that the high expression of PABPC1 would promote cancer tumorigenesis and resistance [53] and in lung adenocarcinoma, it might be involved in tumor development [54]. In follicular thyroid cancer, it was identified as a mutated cancer driver gene. Except for PABPC1, PABPC3 was also found as a cancer driver gene in follicular thyroid cancer [55]. For PABPC4, it exhibited significant differences in expression in immune cell-infiltrating breast tumors, hopefully regarding as a candidate for its diagnostics and therapy [56]. Besides, Liu et al. pointed that PABPC4 likely played a role in the pathogenesis of colorectal cancer [57]. In various types of cancers, numerous studies have now documented a link between EIF4A3, SMC3, ESPL1, CDC25B, CCNA2, or AURKA and their response to therapy [58][59][60][61][62][63]. For the downregulation core gene THBS3, conversely, it was expressed at significantly 9 Disease Markers high levels in osteosarcoma, which was a predictor of worse OS at diagnosis [64]. All in all, there is few report of the above genes in EC.
In order to find candidate drugs that may provide novel insights into tumorigenesis therapeutic advancements for EC based on immune cell infiltration, we performed CMap analysis. As a result, we found two compounds Semustine and Sirolimus that were promising candidates for the treatment of EC. Semustine is a nitrosourea antitumor drug, which belongs to nonspecific cell cycle drug. Summarizing the indications, Semustine has a good curative effect on malignant melanoma, lymphoma, brain tumor, and lung cancer. However, up to now, there are few clinical studies of Semustine in patients with EC. The finding thus is expected to expand the application range of Semustine. In addition to such chemotherapeutic, the use of immunosuppressants is an important change in the way the cancer is treated. For EC, at present, most immunosuppressants are still in the clinical research stage. For instance, pembrolizumab is currently in Phase Ib study with PD-L1-positive esophageal cancer patients, exhibiting tolerable toxicity and effective antitumor effect [65]. In clinical II phase, nivolumab monotherapy has achieved good results in patients with advanced EC of which the objective response rate is 17.2% and the disease control rate is 42% [66]. Although in this study, for Sirolimus, we have identified it as a potential candidate inhibitor for EC; a large number of trials are needed to evaluate its antitumor activity and safety.
Finally, we have to admit some limitations of the study. Owing to the insufficiency of clinical data and independent validation cohort, we cannot further prove our finding. We have endeavoured to collect as many samples as possible from public database. Regretfully, we just found a validation to preliminarily verify the value of the risk score. Besides, we set about collecting sample of EC in our own institute. Our findings would be further validated with basic biology experiment.

Conclusion
In EC, 20 TIICs were identified as having a statistical difference between normal and tumor tissues. These TIICs had great links with clinicopathological characters. A five-TIIC risk score model was built to guide prognosis of EC. Based on the enrichment analysis of stratification by the risk score, two candidate compounds Semustine and Sirolimus were found as therapeutic strategy for EC. In general, the proportions of TIICs are likely to be important to the progression, prognosis, and treatment of EC.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no competing interests.