Identification of potential biomarkers associated with immune infiltration in the esophageal carcinoma tumor microenvironment

As a kind of gastrointestinal cancer, esophageal carcinoma (ESCA) ranks seventh according to incidence (572000 new cases) and sixth in mortality totally (509000 deaths) in 2018 [1]. The progression and chemoresistance of ESCA are important factors that impact patient survival. For those patients with locally advanced ESCA, the 5-year survival rate is approximately 20%. The standard treatment includes chemotherapy or chemoradiotherapy with surgery [2]. In recent years, tumor immunologic therapies have also shown promising early results in ESCA. And tumor immune cell infiltration significantly correlates with the progression and the effect of immunotherapy in cancers [3,4]. However, the specific molecular mechanism about immunotherapy for ESCA yet remains unclear and some genes associated with immune infiltration are required to be found.

With the development of gene chip technology and sequencing technology, bioinformatics is becoming more and more important [5–7]. Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) are two open databases, which provide gene chip data and sequencing data about the cancers. A lot of studies were published for analyzing the data from GEO and TCGA in ESCA, including the profiles of long noncoding RNA (lncRNA), microRNA (miRNA), and messenger RNA (mRNA). Interleukin 8 (IL8) and C–X–C chemokine receptor type 7 (CXCR-7) were predicted to differentially express in ESCA tissues compared with normal esophageal tissues and significantly correlate to the prognosis of patients with ESCA [8]. Liu et al. constructed the network of lncRNA, miRNA, and mRNA and found eight lncRNAs including WDFY3-AS2, CASC8, UGDH-AS1, RAP2C-AS1, AC007128.1, AC016205.1, AC092803.2, and AC079949.2 correlated with overall survival of the patients with ESCA [9]. Zeng et al. built the regulatory networks of miRNA-TF and miRNA-gene and identified the top five miRNAs (miRNA-93, miRNA-21, miRNA-4746, miRNA-196a-1, and miRNA-196a-2) with diagnostic value [10]. However, the immunotherapy-related genes were certainly not identified in these studies in ESCA.

In the present study, the common differentically expressed genes (DEGs) were screened between esophageal squamous cell carcinomas (ESCC) and esophageal adenocarcinoma (EAC) by analyzing the datasets from GEO. Then the gene ontology (GO) analysis, Reactome pathway analysis, and protein–protein interaction (PPI) analysis of common DEGs were performed. Subsequently, COL1A2 with its similar genes were selected in Gene Expression Profiling Interactive Analysis (GEPIA) and validated in UALCAN, meanwhile the prognostic significance analysis was performed. Finally, the immune infiltration analysis was carried out to explain the immunological significance of COL1A2 with its similar genes. The findings in this report revealed the important role of COL1A2 with its similar genes in regulating immune infiltration and suggested that these genes might be munotherapy-related genes in ESCA.

Two microarray datasets (GSE20347 and GSE26886) were downloaded from GEO. GSE20347 was based on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array) which contained 34 samples, including 17 pairs of tumors and their matching normal adjacent tissues from ESCC patients. GSE26886 was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) which contained 40 samples, including 21 EAC samples and 19 normal esophageal squamous epithelium samples.

We downloaded the raw CEL data and standardized them by the Affy package of R language [11]. Then the differential gene analysis was carried out using the limma R package and genes with an adjusted P-value <0.01 and |log2 fold change (FC)| > 2 were considered DEGs [12]. The heatmaps and volcano plots of DEGs were displayed. The common DEGs of two datasets were integrated using an online web tool, Venny 2.1, and the common DEGs including up-regulated and down-regulated genes were saved for subsequent analysis.

The GO PANTHER version 14.0 was used to perform GO functional and Reactome pathway analysis of the common DEGs [13]. And the GO functional analysis of integrated DEGs refers to three parts: biological process (BP), cell component (CC), and molecular function (MF). FDR < 0.05 was considered to have a statistically significant difference.

We constructed a PPI network of the common DEGs using the Search Tool for the Retrieval of Interacting Genes (STRING), an online web tool for predicting the correlation of protein function [14]. And the PPI network was visualized by Cytoscape software [15]. Then, the most significant module of up-regulated DEGs was selected by Molecular Complex Detection (MCODE), which was a plug-in application of Cytoscape [16].

UALCAN can analyze cancer OMICS data which is a comprehensive, user-friendly, and interactive web tool [17]. The expressions of COL1A2 and its similar genes were validated in primary ESCA samples compared with normal esophageal squamous epithelium samples. P<0.05 was considered to indicate a statistically significant difference.

GEPIA is a web-based tool used to deliver fast and customizable functionalities based on TCGA and GTEx data. Its primary functions included differential expression analysis, patient survival analysis, similar gene detection, and so on [18]. Similar genes and prognostic significance analysis were executed here to study our target genes and P<0.05 was considered to have a statistically significant difference.

Tumor Immune Estimation Resource (TIMER) is an interactive web tool to comprehensively investigate the molecular characterization of tumor-immune interactions. Six major functional modules allow users to interactively explore the correlations between immune infiltrates and a wide range of factors, including gene expression, somatic mutations, and somatic copy number alterations, and so on [19]. The correlations between the expression levels of COL1A2 with its similar genes and the infiltration levels of various immune cells were analyzed in ESCA, including B cells, CD8⁺ T cells, CD4⁺ T cells, neutrophils, macrophages, and dendritic cells (DCs). And we also analyzed the relationships between the expression levels of COL1A2 with its similar genes and the immune marker sets including CD8⁺ T cell, T cell (general), B cell, monocytes, tumor-associated macrophages (TAMs), M1 macrophages, M2 macrophages, neutrophils, natural killer (NK) cells, DCs, T-helper (Th) cells, T-helper 17 (Th17) cells, follicular helper T (Tfh) cells, regulatory T cells (Tregs) and exhausted T cells based on previous studies [20–22].

To find the common DEGs from the two subtypes of ESCA, we analyzed the ESCC chip expression dataset GSE20347 and EAC chip expression dataset GSE26886 (adjusted P<0.01 and |log2 FC| > 2). The GSE20347 dataset contained 225 DEGs, including 61 up-regulated genes and 164 down-regulated genes (Figure 1A,C). The GSE26886 dataset contained 691 DEGs, including 296 up-regulated genes and 395 down-regulated genes (Figure 1B,D). The common up-regulated DEGs and down-regulated DEGs from the two kinds of datasets were identified respectively, which were displayed by the Venn diagram, including 19 up-regulated genes and 109 down-regulated genes (Figure 1E,F).

DEGs GO and Reactome pathway analysis were performed with PANTHER. The data shown in the figure needed to meet the following criteria: (1) the FDR < 0.05; (2) the top five significant enrichment results. We can find the BP of up-regulated genes which were mainly enriched in the extracellular structure organization, extracellular matrix (ECM) organization, and collagen fibril organization. The down-regulated genes were mainly enriched in the cornification, epidermis development, keratinocyte differentiation, epidermal cell differentiation, and keratinization (Figure 2A). And the CC of up-regulated genes were mainly enriched in the complex of collagen trimers, ECM, collagen trimer, banded collagen fibril, and fibrillar collagen trimer. The down-regulated genes were mainly enriched in the cornified envelope, extracellular region, vesicle, extracellular organelle, and extracellular space (Figure 2B). The MF of up-regulated genes were mainly enriched in the ECM structural constituent, ECM structural constituent conferring tensile strength and structural molecule activity, but the results of down-regulated genes were not displayed because the FDR > 0.05 (Figure 2C). As for the pathway analysis, the up-regulated DEGs were significantly enriched in the ECM proteoglycans, collagen degradation, ECM organization, degradation of the ECM, and integrin cell surface interactions, while the down-regulated DEGs were significantly enriched in the formation of the cornified envelope, keratinization, and developmental biology (Figure 2D).

The PPI network of DEGs was constructed by using the Cytoscape, which contained 81 nodes and 174 edges (Figure 3A). The most significant module of up-regulated DEGs was identified by using MCODE, which contained COL1A2, COL5A2, COL4A2, LUM, and COL4A1 (Figure 3B).

COL1A2 was a subtype of Type I collagen, the pieces of evidence suggest its mRNA expression in ESCC tumors was higher than that in normal tissues. And it was associated with the survival of ESCC patients [23]. In the following study, we selected COL1A2 for our important target gene. In the first place, we explored the top ten similar genes of COL1A2 in ESCA by using GEPIA (Table 1). In order to validate our results, we analyzed the differential expression of COL1A2 and its similar genes between 11 normal esophageal squamous epithelium samples and 184 primary ESCA samples by using UALCAN (Figure 4A–K). Except for CTSK, SPARC, and CLEC11A, all of the genes expressed higher in primary ESCA samples compared with normal esophageal squamous epithelium samples. Meanwhile, we investigated whether these genes contributed to the prognosis of ESCA patients, GEPIA was used to analyze the disease-free survival (DFS) of COL1A2 and its similar genes in patients. As shown in Figure 5A–K, the higher expression of COL1A2 and its similar genes showed worse DFS in ESCA patients, except for COL5A1, LRRC15, ADAM12 and CLEC11A. These genes that have been validated and had prognostic value will be used for later study.

The infiltration of immune cells such as TAMs and Tregs in the microenvironment can promote esophageal carcinogenesis [24]. Therefore, we explored the relationship between the expression of COL1A2 with its similar genes and the infiltrating immune cells in ESCA using the TIMER database. Tumor purity plays an important role in the analysis of immune cell infiltration about clinical tumor samples by genomic approaches [25]. After adjusting the analysis by tumor purity, we found that the expression levels of COL1A2, COL1A1, COL3A1, ZNF469 and Periostin (POSTN) were significantly correlated with tumor purity in ESCA. And the expression levels of these genes were positively correlated with the infiltration levels of macrophages and DCs, as well as the expression levels of ZNF469 was also positively correlated with the infiltration levels of CD4⁺ T cells (Figure 6).

To investigate the relationship between COL1A2 with its similar genes and the different immune infiltrating cells, the correlations between COL1A2 with its similar genes and immune marker sets of various immune cells of ESCA were explored in the TIMER database. We analyzed the correlations between the expression of COL1A2 with its similar genes and immune marker genes of diverse immune cells, including CD8⁺ T cells, T cells (general), B cells, monocytes, TAMs, M1 and M2 macrophages, neutrophils, NK cells, and DCs in ESCA. Besides, the diverse T cells, for instance, Th1 cells, Th2 cells, Tfh cells, Th17 cells, Tregs, and exhausted T cells were also analyzed. Before the correlation analysis, we adjusted it for tumor purity because the tumor purity of clinical samples can impact the immune infiltration analysis. As shown in Table 2, the expression of COL1A2 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CCL2 and IL10), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), DCs (BDCA-4), T-helper (STAT4, GATA3, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). The expression of COL1A1 was significantly correlated with the expression of marker genes from tumor-infiltrating CD8⁺ T cell (CD8B), B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CCL2 and IL10), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), DCs (HLA-DPA1), T-helper (STAT4, STAT5A, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (TIM-3). The expression of COL3A1 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), M1 macrophages (COX2), neutrophils (CD66b and CD11b), DCs (HLA-DQB1 and HLA-DPA1), T-helper (STAT4, GATA3, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). The expression of ZNF469 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), neutrophils (CD66b and CD11b), T-helper (STAT4, STAT6, STAT5A, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). The expression of POSTN was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CD68), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), T-helper (STAT4, GATA3, STAT6, STAT5A, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ) and exhausted T cells (CTLA4 and TIM-3). All of the results demonstrated that COL1A2 with its similar genes were correlated with the infiltration of the immune cells in ESCA.

Over the past decade, there are many studies about bioinformatics analysis that have been published to explore the potential pathological mechanisms in ESCA. However, the mechanism of immune cell infiltration and the genes associated with immune infiltration in ESCA yet remain unclear. We firstly analyzed two subtypes of ESCA respectively, which contained the datasets of ESCC and EAC from GEO to find the common DEGs from the different histologic subtypes. Altogether 128 DEGs including 19 up-regulated genes and 109 down-regulated genes were screened. To demonstrate the functional mechanism of common DEGs in ESCA, we performed GO enrichment and Reactome pathway analysis. The results of GO enrichment and Reactome pathway analysis displayed that the up-regulated DEGs were mainly involved in the regulation of ECM such as extracellular structure organization and ECM organization, while the down-regulated DEGs were mainly involved in the regulation of cornification and keratinocyte differentiation. These results suggest the major variations in ESCA were involved in the remodeling of the tumor microenvironment (TME) and cell differentiation. The TME is an internal and complex environment in the tumor, which includes all kinds of ECM (for instance, collagen and laminin) and stromal cells including fibroblasts, immune cells, and vascular endothelial cells [26,27]. During tumor progression, the TME interacts with tumor cells depending on different kinds of immune cells (for example, DCs, macrophages, and myeloid-derived suppressor cells) to mediate their immunologic tolerance, subsequently impacting the clinical effect of immunotherapy [28,29].

Aiming to find the most significant module involved in immune infiltration, we constructed the PPI network of DEGs and identified the most significant module of up-regulated DEGs by using MCODE. The most significant module of up-regulated DEGs contained COL1A2, COL5A2, COL4A2, LUM, and COL4A1. One of five up-regulated genes, COL1A2, is a subtype of Type I collagen, which can be produced by stromal fibroblasts and cancer cells [30–32]. And COL1A2 plays an important role in the invasion and metastasis of ovarian cancer cells since it promotes the migration of cancer cells by interacting with α2β1-integrin [33]. However, what roles it plays in immune cell infiltration and immunotherapy about ESCA are still indistinct.

Subsequently, we explored the similar genes of COL1A2 in ESCA by using GEPIA. The top ten similar genes were screened including COL1A1, COL3A1, SPARC, ZNF469, POSTN, CTSK, COL5A1, LRRC15, ADAM12, and CLEC11A. We validated the expressions of COL1A2 and its similar genes in primary ESCA samples compared with normal esophageal squamous epithelium samples, while CTSK, SPARC, and CLEC11A were not significantly up-regulated. Among these genes, the high expression of COL1A2, COL1A1, COL3A1, SPARC, ZNF469, POSTN, and CTSK were negatively correlated with the DFS of ESCA patients. Then these genes (COL1A2, COL1A1, COL3A1, ZNF469, and POSTN) that have been validated and had prognostic value were used for immune-related research. The expression levels of COL1A2, COL1A1, COL3A1, ZNF469, and POSTN were positively correlated with the infiltration levels of macrophages and DCs, and the expression levels of ZNF469 was also positively correlated with the infiltration levels of CD4⁺ T cells. These results indicated these genes might be the candidate genes for assessing the immune infiltration levels in ESCA.

Finally, to further explore the relationship between COL1A2 with its similar genes and the infiltration levels of various immune cells comprehensively, we analyzed the correlation between these genes and immune marker sets. The expression of COL1A2 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CCL2 and IL10), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), DCs (BDCA-4), T-helper (STAT4, GATA3, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). The expression of COL1A1 was significantly correlated with the expression of marker genes from tumor-infiltrating CD8⁺ T cell (CD8B), B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CCL2 and IL10), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), DCs (HLA-DPA1), T-helper (STAT4, STAT5A, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (TIM-3). Type I collagen (COL1A1 and COL1A2) was correlated with the resembling of TAMs in non-small cell lung carcinoma (NSCLC) while it was not correlated with the paucity of T-cell accumulation in pancreatic ductal adenocarcinoma (PDAC), which were consistent with our results [34,35]. The expression of COL3A1 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), M1 macrophages (COX2), neutrophils (CD66b and CD11b), DCs(HLA-DQB1 and HLA-DPA1), T-helper (STAT4, GATA3, STAT6, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). COL3A1, a subtype of Type III collagen, which is similar to Type I collagen, and primarily refers to the growth and metastasis of tumors such as osteosarcoma and nasopharyngeal tumors [36,37]. While our results suggested that COL3A1 also can be a potential biomarker for assessing the infiltration of the immune cells in ESCA. The high expression of COL3A1 may promote the growth and metastasis of ESCA by recruiting various immune cells. The expression of ZNF469 was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), neutrophils (CD66b and CD11b), T-helper (STAT4, STAT6, STAT5A, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). Zinc finger protein (ZNF469) gene mainly refers to the pathological process of cornea-related disease and it is similar to three types of collagen including COL1A1, COL1A2, and COL4A1 [38–40]. Interestingly, we firstly found the important roles of ZNF469 in the assessment of prognosis and immune cell infiltration in ESCA. This gene may be an important biomarker associated with immune infiltration in the ESCA TME. The expression of POSTN was significantly correlated with the expression of marker genes from tumor-infiltrating B cells (CD19 and CD79A), monocytes (CD86 and CD115), TAMs (CD68), M1 macrophages (COX2), neutrophils (CD66b, CD11b, and CCR7), T-helper (STAT4, GATA3, STAT6, STAT5A, BCL6, and IL17A), Treg (FOXP3, CCR8, and TGFβ), and exhausted T cells (CTLA4 and TIM-3). POSTN is a soluble ECM protein and it is up-regulated in many cancers such as ESCA, head and neck squamous cell carcinoma, and non-small cell lung cancer [41–43]. Recent studies demonstrated that it was an important biomarker associated with the pathogenesis and prognosis of ESCA, and our results further underlined the important role of POSTN in immune cells infiltration [44,45].

Different from previous studies, we analyzed the common DEGs between the ESCC and EAC to find the most important variations in ESCA and identified the important roles of COL1A2 with its similar genes in influencing the immune infiltration in ESCA. However, some limitations still exist in the present study. Firstly, some technical defects of gene chip and sequencing technology existed here. Secondly, we only analyzed the potential biomarker in ESCA and the specific mechanism of immune cell infiltration induced by these genes are needed to be demonstrated in vivo and in vitro experiments.

In summary, we identified COL1A2 and its similar genes as the potential biomarkers for assessing the infiltration levels of immune cells using a comprehensive analysis of bioinformatics, which could help us better understand the mechanism of tumor immune cell infiltration in ESCA. However, more studies are needed to verify these results.

The gene expression data are available in public databases including GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA (https://portal.gdc.cancer.gov).

The authors declare that there are no competing interests associated with the manuscript.

This work was supported by the National Natural Science Foundation of China [grant number 81672688]; the Natural Science Foundation of Hunan Province, China [grant number 2019JJ50359]; the College Students’ Innovative Entrepreneurial Training Plan Program of China [grant number 2020105330254]; and the Graduate research innovation project of Central South University [grant number 2020zzts887].

Z.W., M.C. and Y.Q. designed the present study. Z.W., M.C., Y.Q., Y.Y., Y.H. and X.L. collected the raw data and prepared the analysis. Z.W., M.C. and W.Z. were involved in writing and revising the paper. All authors read and approved the final manuscript.

BP

biological process

CC

cell component

COL1A2

collagen type I alpha 2 chain

DC

dendritic cell

DEG

differentially expressed gene

DFS

disease-free survival

EAC

esophageal adenocarcinoma

ECM

extracellular matrix

ESCA

esophageal carcinoma

ESCC

esophageal squamous cell carcinoma

FC

fold change

FDR

false discovery rate

GEO

Gene Expression Omnibus

GEPIA

Gene Expression Profiling Interactive Analysis

GO

gene ontology

GTEx

Genotype-Tissue Expression

lncRNA

long noncoding RNA

MF

molecular function

miRNA

microRNA

mRNA

messenger RNA

NK

natural killer

POSTN

Periostin

PPI

protein–protein interaction

TAM

tumor-associated macrophage

TCGA

The Cancer Genome Atlas

Tfh

follicular helper T

Th

T-helper

Th17

T-helper 17

TIMER

Tumor Immune Estimation Resource

TME

tumor microenvironment

Treg

regulatory T cell