Abstract

   Neo-chemoradiotherapy (nCRT) before surgery is a standard treatment for
   locally advanced esophageal cancers. However, the treatment outcome of
   nCRT varied with different patients. This study aimed to identify
   potential biomarkers for prediction of nCRT-response in patients with
   esophageal squamous cell carcinoma (ESCC). Microarray datasets of nCRT
   responder and non-responder samples (access number [36]GSE45670 and
   [37]GSE59974) of patients with ESCC were downloaded from Gene
   Expression Omnibus (GEO) database. The mRNA expression profiles of
   cancer biopsies from four ESCC patients were analyzed before and after
   nCRT. Differentially expressed genes (DEGs) and miRNAs were screened
   between nCRT responder and non-responder ESCC samples. Functional
   enrichment analysis was conducted for these DEGs followed by
   construction of protein-protein interaction (PPI) network and
   miRNA-mRNA regulatory network. Finally, univariate survival analysis
   was performed to identify candidate biomarkers with prognostic values
   in ESCC. We identified numerous DEGs and differentially expressed
   miRNAs from nCRT responder group. GO and KEGG analysis showed that the
   dysregulated genes were mainly involved in biological processes and
   pathways, including “response to stimulus”, “cellular response to
   organic substance”, “regulation of signal transduction”, “AGE-RAGE
   signaling pathway in diabetic complications”, and “steroid hormone
   biosynthesis”. After integration of PPI network and miRNA-mRNA network
   analysis, we found eight genes, TNF, AKR1C1, AKR1C2, ICAM1, GPR68,
   GNB4, SERPINE1 and MMP12, could be candidate genes associated with
   disease progression. Univariate cox regression analysis showed that
   there was no significant correlation between dysregulated miRNAs (such
   as hsa-miR-34b-3p, hsa-miR-127-5p, hsa-miR-144-3p, and hsa-miR-486-5p,
   et al.) and overall survival of ESCC patients. Moreover, abnormal
   expression of MMP12 was significantly correlated with pathological
   degree, TNM stage, lymph nodes metastasis, and overall survival of ESCC
   patients (p < 0.05). Taken together, our study identified that MMP12
   might be a useful tumor biomarker and therapeutic target for ESCC.

   Keywords: esophageal neoplasms, differentially expressed genes,
   prognosis, MMP12, neo-chemoradiotherapy

Introduction

   Esophageal squamous cell carcinoma (ESCC) is a primary histological
   type of esophageal cancer worldwide. In the United States,
   approximately 18,440 new cases were diagnosed and 13,100 deaths
   occurred in 2020 ([38]Siegel et al., 2020). China is one of the areas
   with highest incidence of esophageal cancer. According to the cancer
   statistics in China, an estimate of 477,900 newly diagnosed esophagus
   carcinoma cases were found in China, and 375,000 individuals died of
   this disease ([39]Chen et al., 2016). Over the past decades, the
   incidence of esophageal cancer has been decreased slightly in the
   United States, whereas most cases of this disease occurred in
   developing countries. ESCCs in China bear more than half of global
   burden ([40]Bray et al., 2018; [41]He et al., 2019).

   Esophageal cancer has a poor prognosis when it is diagnosed at advanced
   stage. The 5 years survival rate would decrease to 4% when metastasis
   occurred ([42]Shaheen et al., 2017). As for locally advanced ESCC,
   neo-chemoradiotherapy (nCRT) before surgery is the primary management.
   However, patients always suffer from disease recurrence after
   resection. Also, the outcome of nCRT varied among patients. Patients
   who response to nCRT acquire superior survival time while those with no
   response suffer a poor prognosis ([43]Berger et al., 2005). Besides,
   nCRT may increase postoperative complications, and clinical parameters
   (TNM classification, tumor location) could not predict response to nCRT
   ([44]Rohatgi et al., 2005; [45]Wen et al., 2014). Therefore,
   identification of novel biomarkers that could predict response to nCRT
   would be beneficial for ESCC management. In addition, it would help
   clinical doctors to discontinue non-effective treatments, and thereby
   avoid overtreatment for non-responders.

   In recent years, the development of microarray analysis and RNA
   sequencing technology provided favorable information for exploring the
   association of gene expression and clinical outcomes ([46]Meyerson et
   al., 2010; [47]Ross and Cronin, 2011). In this study, based on RNA
   sequencing and bioinformatics analysis, we identified differentially
   expressed genes (DEGs) or differentially expressed miRNA between nCRT
   responders and non-responders in patients with ESCC. Additionally,
   protein-protein interaction (PPI) network integrated with miRNA-mRNA
   network were constructed to explore the candidate genes related to
   disease progression. Cox regression analysis was finally conducted to
   investigate the correlations between survival times and candidate genes
   or miRNAs in ESCC patients. A schematic diagram of bioinformatics
   analysis for ESCC datasets were shown in [48]Figure 1. Our study
   investigated the transcriptomic expression patterns of ESCC, and
   identified novel biomarkers with prognostic values that could provide a
   better understanding of ESCC progression.

FIGURE 1.

   [49]FIGURE 1
   [50]Open in a new tab

   A schematic diagram of bioinformatics analysis for ESCC datasets.

Materials and Methods

Data Source and Data Processing

   The gene expression profile [51]GSE45670 ([52]Wen et al., 2014) and
   miRNA expression profile [53]GSE59974 ([54]Fu et al., 2016) were
   downloaded from the Gene Expression Omnibus (GEO) database. The two
   datasets were derived from the same clinical samples that consist of 17
   treatment non-responders and 11 responders. [55]GSE45670 dataset were
   tested on platform of hg-u133_plus_2 Affymetrix Human Genome U133 Plus
   2.0 and [56]GSE45670 was analyzed on Agilent-038169 human miRNA novirus
   v18.0 platform.

   In addition, cancer biopsies were taken from four patients before and
   after nCRT in Changzhou No.2 People’s Hospital, the Affiliated Hospital
   of Nanjing Medical University. This study was approved by the Ethical
   Research Committee of hospital. Transcriptomics data (under access
   number [57]GSE137867,
   [58]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE137867) were
   obtained by RNA sequencing technology and mRNA expression was analyzed
   for these clinical specimens. [59]GSE137867 were considered as
   validation datasets. DEGs were screened from the three datasets and
   overlapped genes were considered as candidate genes related to ESCC.

   Data normalization was performed for these microarray data. We
   extracted the gene expression values from series matrix files, and
   converted corresponding probe ID into gene symbols. After removal of
   abnormal information, we selected the expressed values from ESCC
   responders or non-responders to conduct Principal Components Analysis
   (PCA).

PCA Analysis and DEGs Screening

   PCA is a multivariate technique that converted multiple correlated
   variables into cohort of values linearly, and uncorrelated variables
   represented principal components. It has been extensively used to
   explore high-dimensional data, such as genomic and transcriptomic
   expression data ([60]Ringnér, 2008; [61]Ellsworth et al., 2017). In
   this study, we used PCA method to investigate distribution of samples
   between the experimental group and control groups. After detection and
   removal of abnormal samples, we finally obtained a series of specimen
   with high similarity.

   The miRNA matrix and mRNA expression profiles were normalized based on
   quantile normalization method. Limma package ([62]Ritchie et al., 2015)
   was used to screen DEGs and miRNAs between experimental group and
   normal control group by setting thresholds as p < 0.1 and fold change
   >1.5 (|log2FC| > 0.585). Moreover, with the same cutoff criteria, we
   screened the overlapped DEGs of [63]GSE137867 and GEO datasets to
   conduct further analysis.

Functional Enrichment Analysis for DEGs Associated with ESCC

   Gene Ontology (GO) ([64]Ashburner et al., 2000) and Kyoto Encyclopedia
   of Genes and Genomes (KEGG) ([65]Ogata et al., 2000) signaling pathway
   enrichment analysis were conducted for DEGs. Using fisher's exact test,
   we screened a cohort of biological processes and pathway categories
   enriched by DEGs. p < 0.05 was considered as significant difference,
   and the column with smaller p value represented a closer association
   between DEGs and pathway categories.

Prediction of miRNA-mRNA Interactions

   We further predicted the miRNA-mRNA interactions by examining various
   databases [TargetScan ([66]Agarwal et al., 2015), miRTarBase
   ([67]Chih-Hung et al., 2017), miRDB ([68]Nathan and Wang, 2014) and
   miRanda ([69]Doron et al., 2008)]. As for the upregulated miRNAs, we
   explored their correlations with down-regulated target genes.
   Subsequently, we focused on the downregulated miRNAs, and investigated
   their correlations with upregulated target genes in nCRT responders.
   Cytoscape software was used to visualize the relationships of
   differentially expressed miRNAs and mRNAs.

Construction of Regulatory Network

   The online tool STRING ([70]https://string-db.org/) ([71]Christian et
   al., 2003) was used to analyze the interactions of proteins.
   Interaction score >0.7 (high confidence) was set as the cut-off
   criteria. After screening the correlation pairs, Cytoscape ([72]Shannon
   et al., 2003) was used to analyze the topological properties
   (connectivity degree, closeness degree and betweenness degree) of PPI
   network. The proteins with high scores were considered as hub factors
   in network and might be key candidate genes in disease progression.

   In addition, we integrated the PPI network and miRNA-mRNA regulatory
   network by using Cytoscape software to predict the candidate miRNAs
   involved in ESCC progression.

Survival Analysis

   Survival package ([73]Cox, 1972) in R software was used to construct
   cox regression model. Univariate cox regression analysis and survival
   analysis were performed to identify crucial genes and miRNAs related to
   ESCC survival based on the TCGA datasets.

Results

Screening DEGs From nCRT Responders of ESCCs

   PCA analysis and correlation analysis on samples were conducted for
   mRNA expression dataset [74]GSE45670. Under the threshold of p < 0.1
   and fold change >1.5 (|log2FC| > 0.585), we screened a total of 1311
   DEGs from nCRT responders and non-responders. We found 672 upregulated
   and 639 down-regulated genes. The same cutoff criteria were applied for
   RNA-Seq profile [75]GSE137867, and numerous DEGs were screened in
   specimens before and after nCRT.

   Moreover, we integrated the overlapped DEGs between RNA-Seq profile and
   GEO datasets. Of these gene overlaps, 17 genes were both upregulated in
   two datasets while eight genes were both downregulated in two datasets
   ([76]Table 1). Volcano plot and bidirectional clustering analysis were
   performed for these DEGs, and results were shown in [77]Figure 2.
   Furthermore, the differentially expressed miRNAs were identified by
   setting the thresholds of P < 0.1 and fold change >1.5 (|log2FC| >
   0.585). Scatter plot and volcano plot were used to visualize the
   results. Finally, we identified 44 upregulated miRNAs and 15
   downregulated miRNAs from the nCRT responders compared with
   non-responders. The differentially expressed miRNAs are listed in
   [78]Table 2 and heatmaps are visualized in [79]Figure 3.

TABLE 1.

   Twenty five DEGs were identified from two mRNA profiles ([80]GSE45670
   dataset and [81]GSE137867 dataset), including 17 upregulated genes and
   8 down-regulated genes in the ESCC responder groups compared with
   control non-responder groups.
   DEGs Gene names
   Up-regulated ANO4, BMP2, DTNB, GADD45A, GAS1, GNB4, GPR68, ICAM1, IL24,
   MAN1A1, MMP12, NRBF2, SERPINE1, SLC31A2, SNX10, TNF, TNNT1
   Down-regulated AKR1C1, AKR1C2, PCTP, PER2, RAB40B, TLE2, ZDHHC11,
   ZNF703
   [82]Open in a new tab

FIGURE 2.

   [83]FIGURE 2
   [84]Open in a new tab

   Screening the differential expression genes (DEGs) in ESCC-nCRT
   responder and non-responder samples. (A) Volcano plot visualized the
   distribution of DEGs. Red dots represented up-regulated genes and blue
   dots were down-regulated genes. The overlapped genes were named in
   figures represented the up-regulated or down-regulated genes in both
   datasets. (B) Bidirectional clustering analysis of DEG in nCRT
   responder and non-responder samples. These genes are overlapped genes
   dysregulated in both datasets. The color changed from blue to red
   represented the expression level from low to high.

TABLE 2.

   Fifety nine differential expressed miRNAs were identified from miRNA
   profiles ([85]GSE59974 dataset), including 44 upregulated genes and 15
   down-regulated genes in the ESCC responder groups compared with control
   groups.
   miRNAs Gene names
   Up-regulated hsa-miR-106b-3p, hsa-miR-122-5p, hsa-miR-1252,
   hsa-miR-127-5p, hsa-miR-1295a, hsa-miR-1343, hsa-miR-137,
   hsa-miR-155-3p, hsa-miR-18a-3p, hsa-miR-195-3p, hsa-miR-196a-5p,
   hsa-miR-206, hsa-miR-222-5p, hsa-miR-299-3p, hsa-miR-3144-5p,
   hsa-miR-323b-5p, hsa-miR-330-3p, hsa-miR-34b-3p, hsa-miR-3617-5p,
   hsa-miR-363-5p, hsa-miR-373-3p, hsa-miR-3975, hsa-miR-3978,
   hsa-miR-424-5p, hsa-miR-4448, hsa-miR-4528, hsa-miR-4645-3p,
   hsa-miR-4670-5p, hsa-miR-4704-5p, hsa-miR-4709-3p, hsa-miR-4715-5p,
   hsa-miR-4717-5p, hsa-miR-4773, hsa-miR-4774-5p, hsa-miR-4789-5p,
   hsa-miR-4798-5p, hsa-miR-503-5p, hsa-miR-504, hsa-miR-5189,
   hsa-miR-539-3p, hsa-miR-5579-3p, hsa-miR-637, hsa-miR-7-5p, hsa-miR-924
   Down-regulated hsa-miR-138-1-3p, hsa-miR-144-3p, hsa-miR-144-5p,
   hsa-miR-192-3p, hsa-miR-3145-5p, hsa-miR-3972, hsa-miR-486-5p,
   hsa-miR-489, hsa-miR-499a-5p, hsa-miR-550b-2-5p, hsa-miR-551b-3p,
   hsa-miR-5585-5p, hsa-miR-5690, hsa-miR-640, hsa-miR-675-3p
   [86]Open in a new tab

FIGURE 3.

   [87]FIGURE 3
   [88]Open in a new tab

   Volcano plot and heat map visualized the differentially expressed miRNA
   screening between nCRT responders and non-responders. (A, B) Scatter
   plot and Volcano plot represented the differentially expressed miRNAs
   found by correlation analysis of ESCC samples. (C) Clustering analysis
   of differentially expressed miRNAs between nCRT responder and
   non-responder samples. The color changed from blue to red represented
   the expression level from low to high.

GO and KEGG Analysis for These Candidate Genes

   The GO analysis results revealed that these dysregulated DEGs were
   mainly associated with several biological processes ([89]Table 3), such
   as “response to stimulus” (count = 18, FDR = 0.001617666), “regulation
   of signaling” (count = 12, FDR = 0.001658521), “response to organic
   substance” (count = 11, FDR = 0.001658521), “cellular response to
   organic substance” (count = 11, FDR = 0.001658521), and “regulation of
   signal transduction” (count = 11, FDR = 0.002048231). The KEGG pathways
   ([90]Table 4) enriched by these genes were “AGE-RAGE signaling pathway
   in diabetic complications” (count = 3, p = 5.23e−04), “Epstein-Barr
   virus infection” (count = 3, p = 3.91e−03), “steroid hormone
   biosynthesis” (count = 2, p = 3.91e−03), “African trypanosomiasis”
   (count = 2, p = 1.63e−03) and “malaria” (count = 2, p = 2.85e−03).

TABLE 3.

   Gene Ontology analysis of differentially expressed genes associated
   with ESCC (Top ten biological process terms).
   Category Description Count p value FDR
   BP GO:0071396∼cellular_response_to_lipid 6 2.80e−06 0.001617666
   BP GO:0050896∼response_to_stimulus 18 2.84e−06 0.001617666
   BP GO:0071395∼cellular_response_to_jasmonic_acid_stimulus 2 3.63e−06
   0.001617666
   BP GO:0023051∼regulation_of_signaling 12 6.06e−06 0.001658521
   BP GO:0010033∼response_to_organic_substance 11 7.33e−06 0.001658521
   BP GO:0071310∼cellular_response_to_organic_substance 10 7.44e−06
   0.001658521
   BP GO:0009966∼regulation_of_signal_transduction 11 1.32e−05 0.002048231
   BP GO:0071222∼cellular_response_to_lipopolysaccharide 4 1.46e−05
   0.002048231
   BP GO:0030155∼regulation_of_cell_adhesion 6 1.46−-05 0.002048231
   BP GO:2000351∼regulation_of_endothelial_cell_apoptotic_process 3
   1.54e−05 0.002048231
   [91]Open in a new tab

   Category stands for GO terms and BP refers to biological process.

TABLE 4.

   KEGG pathway analysis of differentially expressed genes in ESCC.
   Category Description Count p value FDR Genes
   KEGG hsa04933∼AGE-RAGE_signaling_pathway_in_diabetic_complications
   Homo_sapiens_(human) 3 5.23e−04 0.058021807 ICAM1/SERPINE1/TNF
   KEGG hsa05143∼African_trypanosomiasis Homo_sapiens_(human) 2 1.63e−03
   0.08638883 ICAM1/TNF
   KEGG hsa05144∼Malaria Homo_sapiens_(human) 2 2.85e−03 0.08638883
   ICAM1/TNF
   KEGG hsa05169∼Epstein-Barr_virus_infection Homo_sapiens_(human) 3
   3.91e−03 0.08638883 GADD45A/ICAM1/TNF
   KEGG hsa00140∼Steroid_hormone_biosynthesis Homo_sapiens_(human) 2
   4.24e−03 0.08638883 AKR1C1/AKR1C2
   KEGG hsa05217∼Basal_cell_carcinoma Homo_sapiens_(human) 2 4.67e−03
   0.08638883 BMP2/GADD45A
   KEGG hsa04115∼p53_signaling_pathway Homo_sapiens_(human) 2 6.06e−03
   0.095935394 GADD45A/SERPINE1
   KEGG hsa05323∼Rheumatoid_arthritis Homo_sapiens_(human) 2 9.54e−03
   0.095935394 ICAM1/TNF
   KEGG hsa04350∼TGF-beta_signaling_pathway Homo_sapiens_(human) 2 0.01
   0.095935394 BMP2/TNF
   KEGG hsa04713∼Circadian_entrainment Homo_sapiens_(human) 2 0.01
   0.095935394 GNB4/PER2
   [92]Open in a new tab

PPI Network Analysis and Prediction of miRNA-Gene Interactions

   By predicting the target genes of the upregulated miRNAs and
   downregulated miRNAs, we built two miRNA-target networks for
   upregulated miRNAs ([93]Figure 4A) and downregulated miRNAs ([94]Figure
   4B). In the regulatory networks, several molecules were identified as
   hub genes for high degree of connectivity, including PER2 (degree = 6),
   PCTP (degree = 6), and hsa-miR-486-5p (degree = 6).

FIGURE 4.

   [95]FIGURE 4
   [96]Open in a new tab

   Identification of candidate genes based on protein-protein interaction
   (PPI) network analysis and miRNA-mRNA regulatory network construction.
   (A, B) Regulatory network visualized the correlation of dysregulated
   miRNA and mRNAs in ESCC samples. The circles represent mRNA while
   diamonds refer to miRNAs. Red mean to up-regulated genes and blue refer
   to down-regulated gene. (C) PPI network analysis for these
   differentially expressed genes. A red dot represents an up-regulated
   gene and a blue dot is a down-regulated gene. (D) The miRNA-mRNA
   regulatory network analysis to identify crucial genes related to ESCC
   progression. Red colors represent up-regulated while blue colors refer
   to down-regulated genes. Diamond and circles represent miRNA and mRNA.

   Using the online tool STRING and Cytoscape software, we constructed a
   PPI network for DEGs ([97]Table 5; [98]Figure 4C) and obtained five
   interaction pairs among eight genes, including tumor necrosis factor
   (TNF), Aldo-keto reductase family one member C1 (AKR1C1), AKR1C2,
   intercellular adhesion molecule 1 (ICAM1), ovarian cancer G-protein
   coupled receptor 1 (GPR68), guanine nucleotide-binding protein subunit
   beta-4 (GNB4), plasminogen activator inhibitor-1 (SERPINE1) and matrix
   metalloproteinase-12 (MMP12).

TABLE 5.

   Protein-protein interaction network analysis for differential expressed
   genes in ESCC samples based on connectivity degrees evaluation.
   Gene     Degree
   TNF        3
   AKR1C1     1
   AKR1C2     1
   ICAM1      1
   GPR68      1
   GNB4       1
   SERPINE1   1
   MMP12      1
   [99]Open in a new tab

   Only AKR1C1 and AKR1C2 were down-regulated genes while other genes were
   upregulated.

   By integrating PPI network and miRNA-mRNA network, we finally
   constructed a regulatory network that consisted with miRNAs and genes.
   Among these genes, PER2, PCTP, TNF and hsa-miR-486-5p exhibited higher
   interactions than other genes. Thus, these genes were identified as hub
   genes related to ESCC progression.

Identifying Crucial Genes with Prognostic Values in nCRT Responders

   We conducted a survival analysis on independent cohort of patients (the
   data were derived from TCGA database) to analyze the prognostic
   significance of 8 genes in ESCC patients. The analyzed signatures
   included genes TNF, AKR1C1, AKR1C2, and ICAM1, et al. ([100]Figure 5)
   The patients were divided into high expression group and low expression
   group according to the median of gene expression. Survival probability
   was compared between two groups to predict candidate genes with
   prognostic value. The results revealed that patients with high level of
   MMP-12 exhibited poorer prognosis than patients in low expression
   groups (hazardous ratio for survival probability, 1.737; p < 0.05).

FIGURE 5.

   [101]FIGURE 5
   [102]Open in a new tab

   Kaplan–Meier survival curves of candidate genes in ESCC, including TNF,
   AKR1C1, ICAM1, MMP12, et al. Red line represents high expression of
   crucial genes while green line refers to low expression genes. The X
   axis represent overall survival time (day), and Y axis means survival
   probability.

   Moreover, univariate cox regression analysis and survival analysis were
   performed for these differentially expressed miRNAs, such as
   hsa-miR-34b-3p, hsa-miR-127-5p, hsa-miR-144-3p, and hsa-miR-486-5p et
   al. However, the results revealed that there was no significant
   correlation between any dysregulated miRNAs and prognosis of ESCC
   patients ([103]Figure 6).

FIGURE 6.

   [104]FIGURE 6
   [105]Open in a new tab

   Univariate cox regression analysis and survival analysis for
   differentially expressed miRNAs in ESCC responder samples, including
   hsa-miR-34b-3p, hsa-miR-127-5p, hsa-miR-144-3p, etc. The X axis
   represent overall survival time (month), and Y axis means survival
   probability. p < 0.05 represent a significant difference.

Discussion

   In this study, we attempted to identify potential tumor biomarkers for
   prediction of nCRT-response in patients with ESCC. After integrating
   overlapped genes between RNA sequencing dataset and GEO datasets, we
   identified 17 genes that were upregulated, and eight genes that were
   downregulated in these two datasets. After verification of the data set
   samples used in the project analysis and the actual number of samples
   included in the analysis, no samples were excluded from the project,
   and all samples were analyzed. Functional enrichment analysis showed
   that dysregulated genes were mainly involved in “response to stimulus”,
   “regulation of signal transduction”, “AGE-RAGE signaling pathway in
   diabetic complications”, “steroid hormone biosynthesis”. According to
   analysis of PPI network and miRNA-mRNA network, we identified eight
   candidate genes, TNF, AKR1C1, AKR1C2, ICAM1, GPR68, GNB4, SERPINE1 and
   MMP12, that were associated with progression of ESCC progression.
   Finally, univariate regression analysis in ESCC patients revealed that
   high expression of MMP12 was significantly correlated to poor
   prognosis.

   Previous studies have explored the difference of gene expression
   between nCRT responder and non-responder samples for prediction of nCRT
   response in esophageal cancer. [106]Maher et al. (2009) identified
   eight genes that were differentially expressed in patients responded to
   treatment. Based on analysis using predictive model, five genes were
   able to predict response to nCRT with high accuracy (95%) in a large
   proportion of esophageal cancer patients. In a recent study ([107]Chen
   et al., 2015), researchers identified nine dysregulated genes between
   treatment responders and non-responders. Functional enrichment analysis
   showed that four of the genes, miR-422, CDK4, Cyclin D2, and E2F3, were
   mainly related to G1/S checkpoint which could regulate tumor
   sensitivity to nCRT. Els [108]Visser et al. (2017) assessed the
   association of gene expression and clinical outcome in 22 studies of
   esophageal cancer, and they found a large heterogeneity in gene
   expression, response to nCRT, and lymph node metastasis.

   Our results demonstrated that eight genes, including TNF, AKR1C1,
   AKR1C2, ICAM1, GPR68, GNB4, SERPINE1 and MMP12, were candidate genes
   associated with ESCC progression. AKR1C1/C2 encodes enzymes belong to
   aldo/keto reductase superfamily. These enzymes are crucial in drug
   resistance of cancer cells and related to metabolism of polycyclic
   aromatic hydrocarbons ([109]Selga et al., 2008). AKR1C1 and AKR1C2
   share a high degree of homology in AKR1C subfamily and were different
   in seven amino acids. Previous study showed that AKR1C1/C2 were
   associated with EDHB-induced inhibition of esophageal cancer cell
   proliferation and this might promote treatment of esophageal cancer
   using EDHB, a substrate of the enzymes ([110]Li et al., 2016).

   In the miRNA-mRNA regulatory network, we identified PER2, PCTP, TNF and
   hsa-miR-486-5p as hub genes related to disease progression. MiRNA have
   been confirmed to have a major role in cancer development by regulating
   the expression of oncogenes or cancer suppressor genes. Aberrant
   expression of special miRNAs have been detected in ESCC patients such
   as upregulated of miR-10b ([111]Tian et al., 2010), miR-21 ([112]Mathe
   et al., 2009), miR-26a ([113]Shao et al., 2016) and downregulated
   miR-125b, miR-203, and miR-205 ([114]Matsushima et al., 2011; [115]Hu
   et al., 2016; [116]Fan et al., 2018). A recent study showed that
   downregulation of miR-486-5p was reported in esophageal carcinoma
   samples and it might function as a tumor suppressor gene in disease
   metastasis via regulating cellular migration ([117]Yi et al., 2016). By
   tissue microarrays analysis of 185 ESCCs samples, Ren et al. also found
   that decreased miR-486-5p was identified in 66.2% of cases and abnormal
   expression of miR-486-5p were related to prognosis of esophageal
   carcinoma ([118]Ren et al., 2016). Our results predicted down-regulated
   miR-486-5p, interact with target gene ICAM1 played critical roles in
   ESCC progression. ICAM1 or CD54 is a 90 kDa glycosylated transmembrane
   protein highly expressed on endothelial cells and mesenchymal stem
   cells. It plays major roles in cell metastasis, proliferation and
   multiple cellular immune response. Dysregulation of ICAM1 has been
   confirmed in liver cancer and ESCC stem cells ([119]Liu et al., 2013;
   [120]Tsai et al., 2015). ICAM1 promotes epithelial-to-mesenchymal
   transition by regulating metastasis-related genes in ESCC cells. ICAM1
   was identified as a target gene of lncRNA-ECM and involved in
   development and progression of ESCC ([121]Yao et al., 2018). However,
   the role of ICAM1 and miR-486-5p remains unclear, especially the
   interaction of the two molecules. A previous study reported that
   down-regulation of miR-486-5p could suppress tumor metastasis by
   regulating metastatic mediator of ICAM-1 in breast cancer
   ([122]Abdallah et al., 2017). It is speculated that ICAM-1 might be
   regulated by miR-486-5p and played a crucial role in development and
   metastasis of ESCC. Therefore, inhibition of ICAM-1 might be a
   potential strategy for ESCC treatment.

   In addition to hsa-miR-486-5p, survival analysis was performed for
   several other differentially expressed miRNAs, such as hsa-miR-34b-3p,
   hsa-miR-127-5p, and hsa-miR-144-3p, et al. in tumor samples. However,
   there was no significant correlation between any dysregulated miRNAs
   and prognosis of ESCC patients. According to literature, miR‐127 was
   downregulated in several types of cancers, including gastric cancer
   ([123]Guo et al., 2013), glioblastoma ([124]Jiang et al., 2014), and
   hepatocellular carcinoma ([125]Huan et al., 2016). In ESCC patients,
   Gao et al. revealed that miR-127 acted as a tumor suppressor in tumors
   by regulating oncogene Formin-like 3 (FMNL3) ([126]Gao et al., 2016).
   In another study, seven serum miRNAs were identified as ESCC
   biomarkers, including miR-127-3p ([127]Zhang et al., 2010). Also,
   tissue or serum miR-144 expression were evaluated in gastric cancer and
   low miR-144 expression was found to predict a poor prognosis in
   gastrointestinal cancer ([128]Liu et al., 2017). Moreover, miRNA-34b
   played an oncogenic role in ESCC development, the polymorphisms of
   rs4938723/pri-miR-34b/c were associated with ESCC susceptibility based
   an analysis on a large number of Chinese population ([129]Harata et
   al., 2010; [130]Zhang et al., 2014). Thus, the potential role of these
   miRNAs in response to nCRT in ESCC warrants further studied.

   Protein of MMP-12 is known as human macrophage metalloelastase (HME) or
   macrophage elastase (ME), and belongs to the MMPs family. MMP12 is
   associated with elastin degradation and macrophage migration in various
   diseases, such as chronic obstructive pulmonary disease, skin diseases
   and cancers ([131]Kerkela et al., 2000; [132]Hunninghake et al., 2009).
   However, the function of MMP12 in tumors is controversial. Aberrant
   expression of MMP12 was reported in several types of cancers, including
   hepatocellular carcinoma ([133]He et al., 2018), lung cancer
   ([134]Roman, 2017), colon cancer ([135]Klupp et al., 2016) and
   nasopharyngeal carcinoma ([136]Chung et al., 2014). Overexpression of
   HME/MMP12 mRNA in patients with colorectal carcinoma exhibited a
   significantly better survival outcome compared with patients with
   normal HME/MMP12 mRNA expression ([137]Yang et al., 2001). Cheng et al.
   observed a high level of MMP12/HME protein in patients with gastric
   carcinoma, and overexpression of MMP12 represented a better survival
   rate ([138]Cheng et al., 2010). The anti-tumorigenic effect of MMP12
   might be due to the generation of angiostatin, which is induced by
   MMP12 and could prevent tumor angiogenesis. However, in other types of
   cancers, upregulated MMP12 were reported to involve in short survival
   times. When grouping samples, the analysis is based on the average gene
   expression to distinguish high and low expression, so the sample size
   of high and low expression is different. The longest survival time of
   high expression group is around 2000 days, which is also consistent
   with the survival trend of grouping according to the median. To our
   knowledge, the potential role of MMP12 in ESCC remains uncertain. A
   recent study showed that overexpression of MMP12 was identified in
   cohorts of resectable tumor tissues compared with normal squamous
   epithelium, and overexpression of MMP12 was correlated with poor
   overall survival in ESCCs ([139]Han et al., 2017). Consistent with the
   previous study, our results also demonstrated that overexpression of
   MMP12 was correlated with clinical stage and poor survival outcomes by
   cox regression analysis. In miRNA-mRNA regulatory network, MMP12
   interacting with TNF was identified as hub genes related to disease
   progression. A previous study have reported that the expression and
   secretion of MMP-12 can be regulated by IL-1β and TNF-α in human airway
   smooth muscle cells, and thereby participated in diseases of the
   airway, such as chronic asthma or chronic obstructive pulmonary
   diseases (COPD) ([140]Xie et al., 2005). Yu et al. showed that
   TNFα-activated mesenchymal stromal cells can recruit CXCR2+ neutrophils
   to tumor microenvironment, resulting in upregulation of MMP12 and other
   metastasis-related genes in tumor cells, such as MMP13 and TGFβ, that
   promoted breast cancer metastasis ([141]Yu et al., 2017). However,
   there is few study on the interactions of MMP12 and TNFs in cancer
   development. Our findings suggested that MMP12 might interact with TNFs
   and play a role in ESCC progression.

   There were some limitations in our study. Experimental validation
   should be conducted to identify the exact biological behaviors of
   candidate DEGs or miRNAs in ESCC development. Meanwhile, the number of
   ESCC specimen was limited. Further validation in larger cohorts is
   necessary to investigate the disease predictive value of these genes.

Conclusion

   In conclusion, we identified eight genes, including TNF, AKR1C1,
   AKR1C2, ICAM1, GPR68, GNB4, SERPINE1, and MMP12, as candidate genes by
   performing integrative analysis on gene expression profiles of
   microarray datasets. We found that abnormal expression of MMP12 was
   significantly correlated with pathological degree, TNM stage, lymph
   nodes metastasis, survival time of ESCC patients. Further basic
   experiments and large-scale multi-center clinical research studies are
   required to validate our results since our study was conducted based on
   data analysis.

Data Availability Statement

   The datasets presented in this study can be found in online
   repositories. The names of the repository/repositories and accession
   number can be found below:
   [142]https://www.ncbi.nlm.nih.gov/,GSE137867.

Ethics Statement

   The studies involving human participants were reviewed and approved by
   the Ethics Committee of Changzhou Second Hospital Affiliated to Nanjing
   Medical University. The patients/participants provided their written
   informed consent to participate in this study.

Author Contributions

   Data curation: JW, Formal analysis: ZS, Funding acquisition: JW and JY.
   Resources: PY. Software: JW, JL, and ZS. Writing–original draft: JW.
   Writing–review and editing: JW and JY.

Funding

   We thanked the support from National Natural Science Foundation of
   China (11705095); Young Talent Development Plan of Changzhou Health
   Commission (CZQM2020074); Changzhou Sci and Tech Program (CJ20200053);
   The Key Research and Development Program of Jiangsu Province
   (BE2020637); Wuxi double hundred young and middle-aged medical and
   health top-notch talent project (No. 202014) and 3333 High-level
   personnel training Project (BRA2019025).

Conflict of Interest

   The authors declare that the research was conducted in the absence of
   any commercial or financial relationships that could be construed as a
   potential conflict of interest.

Abbreviations

   GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; OS,
   overall survival; BP, biological process; CC, cellular component; ESCC,
   esophageal squamous cell carcinomas; GEO, gene expression omnibus;
   DEGs, differentially expressed genes; nCRT, Neo-chemoradiotherapy; PPI,
   protein-protein interaction network; log2 FC, log2 fold
   change/logarithm of fold change.

References