Time series expression pattern of key genes reveals the molecular process of esophageal cancer

Esophageal cancer is one of the world’s early cancers with poor diagnosis and high mortality. It has strong invasiveness and a fast growth rate [1]. Dysphagia and unconscious weight loss are the most common clinical symptoms [2]. Patients with esophageal cancer have further increased their difficulty in eating due to the staging and location of the tumor and the poor adjuvant treatment [3]. Studies have shown that there is a clear relationship between the development of esophageal cancer and Helicobacter pylori infection, gastroesophageal reflux disease, smoking and severe alcohol use, as well as diet and other genetic factors [4,5]. From a therapeutic point of view, esophageal cancer can be divided into early esophageal cancer, locally advanced resectable esophageal cancer, locally advanced unresectable esophageal cancer and metastatic esophageal cancer. Because of the anatomical features of esophageal cancer, esophageal cancer is usually detected in the late stage, which vitally affects the treatment and prognosis of patients [6,7]. Endoscopic therapy for early esophageal tumors is effective and safe. Optimal results can be obtained by using endoscopic mucosal resection, ablation therapy, and individualized methods combining both [8]. Treatment for advanced esophageal cancer is limited and may be hampered by the presence of micrometastatic disease [9].

In the development of esophageal cancer, the rs11473 polymorphism of the miR-483-5p binding site plays a vital role in the 3′-UTR of the basigin gene [10]. Single nucleotide polymorphisms (SNPs) in TERT may be associated with susceptibility to esophageal cancer and contribute to the development of esophageal cancer [11]. In terms of regulatory molecules, miR-502 regulates proliferation of esophageal cancer cells by promoting phosphorylation of AKT signaling [12]. B7-H1 can be used as a prognostic factor for human esophageal cancer and may be an important therapeutic target for immunotherapy against this malignant tumor [13]. miR-20b may play an essential role in the tumorigenesis of esophageal cancer by regulating PTEN expression, which may be a potential therapeutic target for the treatment of esophageal cancer [14]. MicroRNA-506 inhibits proliferation of esophageal cancer cells by targeting CREB1 [15]. MiR-21 targets critical proteins in the PTEN/PI3K/AKT signal transduction to promote proliferation, cell migration, cell invasion, and cell cycle, as well as inhibition of cell apoptosis in human esophageal cancer cells [16]. These findings have deepened our understanding of the pathogenesis of esophageal cancer and guided us in the direction of further research. Although the predecessors have reported a series of research results on esophageal cancer, the molecular pathogenesis of the disease is still elusive. To comprehensively and deeply explore the molecular processes of esophageal cancer progression and to explore potential therapeutic targets for the progress of esophageal cancer, we conducted a systematic module analysis. Overall, our work details the role of multifactorial mediated dysfunction modules in the overall growth of esophageal cancer, identifying essential genes and related biological processes (BPs), finding potential molecular mechanisms and therapeutic targets for esophageal cancer.

The Cancer Genome Atlas (TCGA) is a joint project of the National Cancer Institute and the American Human Genome Research Institute. High-throughput genomic analysis technology is a useful tool for people with better understanding of cancer, and it improves their abilities to prevent, diagnose, and treat disease. We first downloaded esophageal cancer RNA-Seq data from the TCGA database and screened ncRNA–mRNA interaction pairs with a score ≥0.5 from RNA associated interaction database (RAID) v2.0 [17] database, including 431937 interacting pairs, involving 5431 ncRNAs. All human transcription factor (TF) target data were downloaded and used in the general database-Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) v2 database for transcriptional studies, including 2492 TFs and 9396 interaction pairs.

In order to explore the molecular process of esophageal cancer staging, we selected four stages of esophageal cancer and normal samples for differential expression analysis, including healthy tissue samples vs stage 1 disease samples, stage 1 disease samples vs stage 2 disease samples, stage 2 disease samples vs stage 3 disease samples, stage 3 disease samples vs stage 4 disease samples. We used the limma package for analysis [18–20]. Using the Correct background function, we performed background correction and normalization on the data. The normalize Between Arrays function quantile normalization method can filter out the control probe and the low expression probe. The differentially expressed genes (DEGs) of the dataset were identified based on the lmFit and eBayes functions (^aP<0.01) using default parameters.

A protein–protein interaction (PPI) system was constructed based on STRING database data (score >500). The gene module of more than 30 nodes was screened throughout the network using the ClusterONE plug-in [21] of the Cytoscape software [22]. We use the Cytoscape plugin CytoHubba [23] to identify hub genes in the module subnet, while CytoHubba contains 12 methods for identifying hub genes. We obtained the top 10 genes and then screened the repeat genes in the 12 sets of genes for survival analysis.

The exploration of the functions and signal transduction involved in genes contributes to the study of the molecular mechanisms of disease. Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the R language Cluster profiler package [24]. The cluster Profiler is a Bioconductor software package that provides statistical analysis of functional clustering of gene sets.

The transcription and post-transcriptional regulation of genes are often dominated by the regulation of TFs and non-coding genes (ncRNA). Therefore, we have scientifically predicted its role in the esophageal cancer dysfunction module. If the regulatory effect between the regulator and the module exceeds 2, and the number of organizational relationships between the regulator and the module is essential (hypergeometric test, P-value <0.01), it can be confirmed that it is a regulator of the critical regulatory module.

For further validation of our findings, several dysregulated genes expression profiles (including FGF14, KCNH6, LOC100506136, RGS7, SH3GL2, THBS4, CPLX2, DPEP1, EPHA5, SCGB1A1, and ST18) were also determined in public GEO dataset of GSE13898, as well as cell (esophageal cancer cell line KYSE-150 and normal esophageal cell line Het-1A, both purchased from the Cell Bank of Chinese Academy of Sciences (Shanghai, China)).

Total RNA from cells were extracted using TRIzol reagent (Invitrogen). First-strand cDNA was generated using Reverse Transcription Reagents (Takara RR047) according to the manufacturer’s protocol. Real-time qPCR was performed in the Real-Time PCR Detection System (Bio-Rad) using SYBR Green (Takara RR820). Primer sequences were summarized in Table 1. Relative expression of genes normalized to GAPDH was calculated with the 2^–ΔΔC_t method. Moreover, the experimental strategy of the present study was shown in Supplementary Figure S1.

Biologists have conducted many experiments and studies on the pathogenesis of esophageal cancer and thus identified potential pathogenic genes for the deterioration of esophageal cancer. To observe molecular changes in the progression of esophageal cancer, we performed differential expression analysis based on RNA-Seq data from four stages of esophageal cancer in the TCGA database. Based on analysis of phase 1 disease samples of normal tissue samples and esophageal cancer, analysis of phase 1 disease samples and stage 2 disease samples, analysis of stage 2 disease samples and phase 3 disease samples, analysis of stage 3 disease samples and stage 4 disease samples, we obtained DEGs associated with each stage of four groups of esophageal cancer. A total of 7457 DEGs were received (Figure 1A). We believe that the presence of these DEGs is closely related to the development of various stages of esophageal cancer. Of the 7457 DEGs, there are 13 common genes (Figure 1B). The genes that continue to be down-regulated are CPLX2, DPEP1, EPHA5, SCGB1A1, and ST18. The genes that are continuously up-regulated are: FGF14, KCNH6, LOC100506136, RGS7, SH3GL2, and THBS4 (Figure 1C).

Gene module analysis helps us to study the complex collaborative relationships between multiple genes. Based on the protein interaction data of STRING database, the interaction network of DEGs was constructed, and 14 functional barrier modules were explored (Figure 2A). Using the 12 methods in Cyto-hubba, a total of 758 hub genes were identified in the interaction subnetwork of the module genes, including the gene SH3GL2, which is continuously up-regulated, in Module 8. Further, 23 hub genes shared by the top 10 gene set in 12 methods were screened for survival analysis. The results show that the gene in GNGT2, module 6, is the related gene (P=0.014) (Figure 2B). A decrease in survival rate accompanied the high expression of GNGT2 gene, and the expression level of GNGT2 gene was negatively correlated with survival rate. Function and pathway are essential mediators of the physiological response of the disease. We performed GC enrichment analysis on 14 module genes (Figure 3A) and KEGG (Figure 3B), and obtained 4844 BPs and 530 cell components (CCs), 773 molecular functions (MFs), and 173 KEGG pathways. The main BPs include positive regulation of protein transport, gastric acid secretion, and insulin-like growth factor receptor binding. The main signal transduction is the p53 signal transduction, the ERBB signal transduction, and the epidermal growth factor receptor (EGFR) signal transduction, which play a crucial role in the dysfunctional module for the functions and pathways involved in multiple genes.

From the perspective of systems biology and systems genetics, transcription and post-transcriptional regulation of genes have long been recognized as crucial regulators of disease development and development, while TFs and ncRNAs are universal regulators of expression and function. Although single or several TF and ncRNA regulation of esophageal cancer progression have been valued by many biologists, few studies have focused on their overall global effect on dysfunctional mechanisms and the role of bridges in development. Therefore, in the present study, based on the targeted regulation relationship between TF and ncRNA on the module gene, we performed a pivotal analysis of the conventional module to explore the crucial regulator that regulates the progression of esophageal cancer. The results showed that a total of 54 TFs involved 54 TF–module target pairs and 853 ncRNAs involved 944 ncRNA–module regulatory pairs. Statistical analysis revealed that TF HIF1A and ncRNA CRNDE regulate the most dysfunctional modules. These crucial TFs and ncRNAs may influence the development and progression of esophageal cancer by mediating dysfunctional modules. Thus, we identified these potential regulatory factors as regulators of dysfunction of esophageal cancer. Notably, has-miR-330-3p regulates the up-regulated DEG SH3GL2 throughout the esophageal cancer process (Figure 4), suggesting that hsa-miR-330-3p plays a crucial role in four stages of esophageal cancer.

DNA methylation often occurs during the development of disease, and methylation factors have multiple regulatory roles. Based on the test results, we identified 600 methylated DEGs in the first stage. The methylation level of these genes was mapped to each module, and methylation factors with the low degree of methylation and high expression were screened (Figure 5). It is speculated that these factors may change the expression level by methylation modification, and then play a role in regulating the development of esophageal cancer.

To confirm the reliability of the identified dysregulated genes, the FGF14, KCNH6, LOC100506136, RGS7, SH3GL2, THBS4, CPLX2, DPEP1, EPHA5, SCGB1A1, and ST18 expression profiles were verified in public GEO dataset of GSE13898. The results showed that EPHA5, LOC100506136, RGS7, THBS4, SCGB1A1, and DPEP1 were dysregulated between esophageal cancer and normal control in GSE13898 validation dataset (all, P<0.05). Moreover, SH3GL2, RGS7, SCGB1A1, ST18, and DPEP1 were found to be dysregulated between atypical hyperplasia of esophagus and normal control (all, P<0.05).

Similarly, as shown in Figure 6, compared with normal esophageal cell line, the expression levels of FGF14, KCNH6, LOC100506136, RGS7, SH3GL2, and THBS4 in esophageal cancer cells were significantly up-regulated in comparison with those in normal esophageal cells (all, P<0.05); however, the mRNA levels of CPLX2, DPEP1, EPHA5, SCGB1A1, and ST18 in esophageal cancer cells were significantly down-regulated (all, P<0.05), which is in accordance with the bioinformatics data above.

Esophageal cancer is one of the most deadly cancers, mainly because it is extremely aggressive and has a poor survival rate. Its 5-year survival rate is approximately 15–25% [1]. The underlying cause of this disappointing low survival rate is that most patients have reached the late stage of detection. For patients with metastatic and unresectable disease, their chances of survival are limited [25]. In the present study, we collected RNA-Seq data from TCGA esophageal cancer and selected four stages of esophageal cancer disease samples and normal samples for differential analysis, and finally obtained four sets of time series DEGs. After screening, we found 13 common genes in four groups of DEGs.

Komatsu et al. [26] studied clinical biomarkers of pulmonary neuroendocrine tumors (LNET) and found that CPLX2 was strongly positive in 16.3% of the examination groups. Importantly, positive CPLX2 expression is associated with lymphatic invasion, pathological staging, and adverse disease-specific survival in LNET patients. It was finally concluded that CPLX2 is a novel clinical biomarker for LNET [26]. In the study of breast cancer diagnostic markers, Fu et al. [27] found that changes in gene expression such as DPEP1 may lead to cancer progression. DPEP1 has been identified as a prognostic gene for colorectal cancer. We found that DPEP1 is overexpressed in CRC. After knocking out the DPEP1 gene, cells (SW480 and HCT116) essentially increased apoptosis and attenuated cell proliferation and cell invasion [28]. In the study of colorectal cancer, Eisenach [29] found that the expression of DPEP1 was essentially increased in colorectal cancer tissues compared with normal mucosa. Zhang et al. [30] also noted the DPEP1 gene in the study of pancreatic ductal adenocarcinoma and found that the gene expression was negatively correlated with histological grade and that lower expression of DPEP1 in tumors was associated with poor survival. Chen et al. [31] analyzed the gastric cancer-associated GEO data and found that THBS4 was up-regulated in patients with recurrent gastric cancer and was positively correlated with the pathological stage and poor prognosis of gastric cancer. THBS4 stimulates the proliferation of gastric cancer cells. The breast-related gene explored by Huang et al. [32] contains the gene THBS4, which is up-regulated in breast cancer. In the study of hepatocellular carcinoma, Su et al. [33] found that knockdown of ThBS4 inhibited migration and invasion of hepatocellular carcinoma cells, as well as hemangiocarcinoma-induced angiogenesis. THBS4 as a target is very promising for the treatment of advanced liver cancer. Both the above genes were present in the differential genes of the four stages of esophageal cancer in the present study and were continuously down-regulated. Moreover, it is identified as a clinical biomarker gene and a therapeutic target gene in various cancers. Therefore, we can reasonably speculate that this gene plays an important role in the occurrence and development of esophageal cancer, providing a reasonable direction for further study of esophageal cancer.

At the modular level, the modular gene is essentially involved in the p53 signal transduction, adrenergic signaling in cardiomyocytes, and insulin secretion. Liu et al. [34] found that lncRNA AK001796 was essentially up-regulated in esophageal squamous cell carcinoma. Cellular proliferation and cell cycle are regulated by regulating MDM2/P53 signaling in esophageal squamous cell carcinoma [34]. In the study of isoleucine anti-esophageal squamous cell carcinoma, it was found that iso-valine can inhibit tumor growth by activating the p53 pathway [35]. Zuev et al. [36] studied the insulin levels and islet β-cell morphology changes in 27 patients with esophageal cancer after one-stage esophagectomy and gastric pedicle graft repair. The secretory activity of insulin pancreatic organs is inhibited after surgery [36]. Insulin-like growth factor-II (IGF-II) secreted by Ld1 overexpressing esophageal cancer xenografts can promote the growth of distant esophageal tumors and promote the metastasis of circulating cancer cells [37]. We screened Hub genes that regulate dysfunctional dysfunction in 14 modules and used Hub genes for survival analysis. It was found that a decrease in survival rate accompanied the high expression of GNGT2 gene, and the expression level of GNGT2 gene was negatively correlated with survival rate, and GNGT2 was a key gene in module 6. Enrichment analysis showed that GNGT2 was involved in the detection of external stimuli and the external components of synthetic cell membranes. With cell differentiation, morphogenesis, and programmed cell death, precise control of organ size is undoubtedly one of the most critical processes in mammal development and regeneration. These processes are tightly regulated by complex, highly coordinated mechanisms to maintain a steady state of growth. Any level of distortion in these organ size regulation processes can lead to a variety of pathological conditions, of which cancer is the most terrifying [38]. Therefore, the completeness of the cell membrane and timely response to external stimuli are important manifestations of the cel l’s normal physiological state. We can speculate that the differential expression of GNGT2 may be involved in the development of cancer.

We predicted 853 ncRNAs and 54 TFs, possibly through a mediator module to participate in the mechanism of esophageal cancer progression. Comprehensive landscape analysis of the staging mechanism of esophageal cancer (Figure 4) revealed that the key gene SH3GL2 of ncRNA has-miR-330-3p regulates module 8, which is also a consensus gene in four stages of DEGs. Enrichment analysis showed that SH3GL2 is involved in the regulation of ERBB and EGFR signal transduction. The ERBB family of receptor tyrosine kinases plays a role in cell adhesion and migration, overexpressing in esophageal squamous cell carcinoma and esophageal adenocarcinoma [39], and is associated with the pathogenesis of esophageal adenocarcinoma [40]. EGFR, a transmembrane receptor kinase, is frequently overexpressed in various types of human cancers, including esophageal cancer. Recently, some studies have shown that EGFR inhibitors may exhibit anti-tumor effects, which are associated with the continued promotion of reactive oxygen species production and induction of apoptosis [41]. Gong et al. [42] studied the effect of pingyangmycin on the expression of EGFR in human esophageal cancer cells and the therapeutic effect of pingyangmycin and cetuximab on esophageal cancer xenografts. Pingyangmycin can down-regulate the expression of EGFR in esophageal cancer cells and enhance the effect of cetuximab on xenograft of esophageal cancer in nude mice [42]. Activation of the EGFR may trigger a series of intracellular signal transduction. It plays an important role in cell proliferation, apoptosis, angiogenesis, and metastasis [43]. It also involves suppressing the immune response by activating regulatory T cells and reducing the level of T cell chemoattractants [44]. The TF HIF1A regulates the key gene ACE of module 6, which is involved in BPs such as arachidonic acid secretion. He analyzed the difference between esophageal adenocarcinoma samples and normal samples and subsequent enrichment analysis and found that the KEGG pathway rich in DEG_S module is mainly in arachidonic acid metabolism, complement and coagulation cascade, and rheumatoid arthritis [45]. Zhi et al. [46] analyzed the difference between esophageal squamous cell carcinoma and normal esophageal squamous epithelium and identified nine of them related to arachidonic acid metabolism.

Till now, with the development of high-throughput sequencing technology, gene expression profiling has been widely used to identify genes or signal pathway related to the development of esophageal cancer [47–52]. However, most studies only identified DEGs or signaling pathways, which might be involved in esophageal cancer development based on different GEO datasets. In the present study, based on the test results above, we further identified 600 methylated DEGs. And the methylation level of these genes was mapped to each module, and methylation factors with the low degree of methylation and high expression were screened. The results of the methylation test showed that the SST gene was up-regulated a large number, which may be a key gene involved in methylation modification to regulate the progression of esophageal cancer. Jin et al. [53] found that hypermethylation of the SST promoter is common and is associated with early tumor progression in Barrett’s esophagus. The SH3GL2 gene is up-regulated. The gene is not only the common DEGs of the four-stage time series but also the Hub gene in module 6. It is also possible to play an important role in the regulation of esophageal cancer by methylation modification. Ghosh et al. [54] studied the effect of SH3GL2 methylation on the pathogenesis of head and neck squamous cell carcinoma, and the disorder of sh3gl2 is an independent pathway for early developmental abnormalities of the head and neck.

However, there are several limitations in the present study. First, among the 13 common DEGs in four groups, only DPEP1, CPLX2, and THBS4 have been mentioned in other literatures, and the remaining ten genes are yet to be studied. Second, we can only speculate that the methylated factors obtained by screening may change the expression level through methylation modification, thereby regulating the occurrence and development of esophageal cancer; however, the specific mechanism is still unclear, which deserves further study.

Based on the esophageal cancer-associated RNA-seq in TGCA, the present study studied the DEGs of esophageal cancer at various stages, constructed a PPI network, obtained 14 dysfunctional modules, and screened out Hub genes. We further performed enrichment analysis to predict ncRNA and TF, as well as methylation analysis of the genes in the module. A series of regulatory factors we predicted have a certain degree of regulation on the potential dysfunction mechanism of esophageal cancer, which provides a new idea for the subsequent study of esophageal cancer.

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

The authors declare that there are no sources of funding to be acknowledged.

Jiafu Wang performed the research. Xiang Xie acquired and analyzed the data. Yurong Sun designed the research and wrote the manuscript.

BP

Biological Process

CC

Cellular Component

CRNDE

Colorectal Neoplasia Differentially Expressed

DEG

Differentially Expressed Gene

EGFR

Epidermal Growth Factor Receptor

ERBB

Epidermal Growth Factor Receptor Family

GEO

Gene Expression Omnibus

HIF1A

Hypoxia-inducible Factor-1

KEGG

Kyoto Encyclopedia of Genes and Genomes

LNET

lung neuroendocrine tumors

PPI

Protein–Protein Interaction

qPCR

quantitative polymerase chain reaction

STRING

Protein-Protein Interaction Networks Functional Enrichment Analysis

TCGA

The Cancer Genome Atlas

TERT

Telomerase Reverse Transcriptase

TF

Transcription Factor