Common gene signatures and key pathways in hypopharyngeal and esophageal squamous cell carcinoma: Evidence from bioinformatic analysis.

BACKGROUND: Hypopharyngeal and esophageal squamous cell carcinoma (ESCC) are the most common double primary tumors with poor prognosis. Intensive work has been made to illuminate the etiology, but the common carcinogenic mechanism remains unclear. Thus, we conducted the study to seek to find the common gene signatures and key functional pathways associated with oncogenesis and treatment in hypopharyngeal squamous cell carcinoma (HSCC) and ESCC by bioinformatic analysis.
METHODS: Three independent datasets (GSE2379, GSE20347, and GSE75241) were screened out from the Gene Expression Omnibus (GEO) database and the overlapping differentially expressed genes (DEGs) were identified using GEO2R online platform. Subsequently, the Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis of DEGs were conducted using database for annotation, visualization and integrated discovery (DAVID). Meanwhile, the protein-protein interaction network (PPI) constructed by search tool for the retrieval of interacting genes (STRING) was visualized using Cytoscape. Afterwards, the most key module and hub genes were extracted from the PPI network using the Molecular Complex Detection plugin. Moreover, the gene expression profiling interactive analysis (GEPIA) was applied to verify the expression differences and conduct the survival analyses of hub genes. Finally, the interaction network of miRNAs and hub genes constructed by encyclopedia of RNA interactomes (ENCORI) was visualized using Cytoscape.
RESULTS: A total of 43 DEGs were identified, comprising 25 upregulated genes and 18 downregulated genes, which were mainly involved in the extracellular matrix-receptor interaction, collagen metabolic, epidermis development, cell adhesion, and PI3K/Akt signaling pathways. Subsequently, 12 hub genes were obtained and survival analysis demonstrated SERPINE1 and SPP1 were closely related to poor prognosis of patients with HSCC and ESCC. Finally, hsa-miR-29c-3p, hsa-miR-29a-3p, and hsa-miR-29b-3p were confirmed as the top 3 interactive miRNAs that target the most hub genes according to the interaction network of miRNAs and hub genes.
CONCLUSION: The common gene signatures and functional pathways identified in the study may contribute to understanding the molecular mechanisms involved in the carcinogenesis and progression of HSCC and ESCC, and provide potential diagnostic and therapeutic targets.

Introduction

Hypopharyngeal squamous cell carcinoma (HSCC) is one of the most common head and neck malignancies worldwide. Due to its frequent regional lymphatic metastasis and delayed diagnosis, the prognosis of HSCC is the worst among head and neck malignancies.[ Besides, compared with other head and neck cancers, second primary esophageal squamous cell carcinoma (ESCC) is the most common in HSCC patients, with an incidence of 10% to 50%.[ Even with aggressive treatment, the prognosis of patients with double primary tumors remains poor, with a 5-year overall survival of only 9% to 11%.[ Thus, it is crucial to develop more efficacious diagnostic and therapeutic strategies.

In recent years, intensive work has been made to illuminate the etiology, the common risk factors containing cigarette, alcohol, and betel nut have been proven to may trigger field cancerization in the hypopharynx and esophagus.[ However, the common molecular mechanisms involved in the carcinogenesis of HSCC and ESCC remains unclear. With the advent of microarray and high throughput sequencing technology, bioinformatic analysis has been widely used to identify the differentially expressed genes (DEGs) and functional pathways involved in the carcinogenesis and progression of tumors, which may contribute to developing effective diagnostic and therapeutic strategies. Therefore, we conducted the study to seek to find the common gene signatures and key functional pathways associated with oncogenesis and treatment in HSCC and ESCC by bioinformatic analysis.

Materials and methods

Microarray datasets search

Datasets containing gene expression differences between HSCC, ESCC, and normal tissues were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo),[ which offers massive public available gene expression data to conduct comprehensive genes analysis. To control the heterogeneity, all included datasets must meet the following inclusion criteria: search term: hypopharyngeal squamous cell carcinoma or esophageal squamous cell carcinoma; sample source: homo sapiens; sample size ≥30; study type: expression profiling by array; publication date: January 1, 2005 to January 1, 2020.

Identification of DEGs

The DEGs between HSCC, ESCC, and normal tissues were extracted and analyzed by GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). GEO2R is an integrative web analysis tool that contributes to identifying DEGs between different groups of samples.[ The cut-off values of DEGs were defined as adjusted P value <.05 and |log fold-change| (|logFC|) values >2. Afterwards, a Venn diagram was drawn to obtain the common DEGs.

GO annotation and KEGG pathways enrichment analyses of DEGs

Gene ontology (GO) analysis is an important part of functional genomics research, which is applied to annotate the biological process (BP), cellular composition (CC), and molecular function (MF) of all genes in the genome.[ The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was applied to clarify relevant signaling pathways of the DEGs involved.[ The Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.7; http://david.ncifcrf.gov) is an online bioinformatic database that provides integrated biometric annotation information of genes.[ To analyze the biological function of DEGs, GO annotation and KEGG pathway enrichment analyses were conducted by DAVID, and P-value < .05 was considered statistically significant.

PPI network construction and key module analysis

Analysis of the functional interactions between proteins could provide important insights for the carcinogenesis and progression of tumors. Thus, the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) was applied to construct the protein–protein interaction (PPI) network of DEGs with a confidence score ≥0.4.[ Afterwards, the PPI network was visualized using Cytoscape software version 3.7.0 (Cytoscape Consortium, San Diego, CA, USA), which is a public platform for visualizing molecular interaction networks from attributing data.[ And the key module in the PPI networks was identified using the Molecular Complex Detection (MCODE) plugin of Cytoscape with degree cut-off = 10, Max depth = 100, node score cut-off = 0.2, and k-score = 2.[

Hub genes selection and analysis

The genes with degree ≥10 in the PPI network were defined as hub genes. Subsequently, the cBioPortal online platform (http://www.cbioportal.org) was utilized to define the coexpressed genes of hub genes according to spearman correlation coefficient >0.8, and the coexpression network was constructed using Cytoscape.[ Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn) is an interactive web server for analyzing gene expression profiling of cancer and normal tissues.[ The GEPIA was applied to verify the expression differences of hub genes between tumor samples with normal tissues, and conduct the survival analyses of hub genes using the Kaplan–Meier method. Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/) is a public web server for comprehensive analysis of immunological, clinical, and genomic features of diverse malignancies.[ The expression levels of hub genes in multiple malignancies were explored using TIMER.

miRNA-hub gene network construction

The aberrant expression of miRNAs plays a crucial role in the coordinate regulation of target gene expression. To further study the interactions between miRNAs and hub genes, the Encyclopedia of RNA Interactomes (ENCORI) (http://starbase.sysu.edu.cn/; version 3.0) was utilized to predict the targeted miRNAs of hub genes. In the study, the targeted miRNAs of hub genes were defined according to the positive results of ≥3 miRNA-target predicting databases, including TargetScan, miRanda, PicTar, and PITA.[ Finally, the interaction network of miRNAs and hub genes was constructed using Cytoscape.

Results

Datasets for research

After the rigorous screening of all relevant datasets, 3 gene expression datasets including GSE2379,[ GSE20347,[ and GSE75241[ met the inclusion criteria. Among them, GSE2379 used 38 samples to detect the gene expression difference between hypopharyngeal squamous cell carcinoma and normal tissues. While the others concentrate on the gene expression differences between esophageal squamous cell carcinoma and normal esophageal tissues. The baseline characteristics of the 3 included datasets are summarized in Table 1.

Table 1

Basic characteristics of the GEO datasets included in the study.

The original data of 3 datasets (GSE2379, GSE20347, GSE75241) were obtained from GEO and then subjected to differential expression analysis by the GEO2R online platform. Based on the predefined cut-off values, DEGs (254 in GSE2379, 292 in GSE20347, and 257 in GSE75241) were identified (Fig. 1A–C). Afterwards, the overlapping DEGs among 3 datasets were identified by drawing a Venn diagram (Fig. 1D), comprising 25 upregulated genes and 18 downregulated genes (Table 2).

Figure 1

Identification of differentially expressed genes (DEGs) in 3 gene expression omnibus datasets. (A) Volcano plot of DEGs in GSE2379. (B) Volcano plot of DEGs in GSE20347. (C) Volcano plot of DEGs in GSE75241. Red, green, and gray color represents the relatively high, low, and equal expression of genes in the corresponding dataset, respectively. (D) Venn diagram of overlapping DEGs from the intersection of the 3 independent datasets.

Table 2

A total of 43 DEGs were identified from the 3 independent gene expression datasets.

GO and KEGG enrichment analyses of DEGs

In order to analyze the biological function of DEGs, GO annotation and KEGG pathways enrichment analysis were conducted by DAVID. The biological processes of DEGs were mainly involved in the extracellular matrix organization, collagen metabolic process, multicellular organismal metabolic process, epidermis development, and cell adhesion (Fig. 2A). The cellular compositions of DEGs mainly include extracellular matrix, extracellular space, extracellular region, collagen, secretory granule, basement membrane (Fig. 2B). The changes in molecular function (MF) of DEGs were mainly concentrated on peptidase activity, polysaccharide binding, pattern binding, extracellular matrix structural constituent, and platelet-derived growth factor binding (Fig. 2C). As for the KEGG functional pathways, DEGs mainly involved in the extracellular matrix-receptor interaction and focal adhesion (Fig. 2D).

Figure 2

GO annotation and KEGG pathways enrichment analyses of the differentially expressed genes (DEGs). (A) Top 20 of the biological process of the DEGs. (B) The cellular compositions of the DEGs. (C) The molecular function of the DEGs. (D) KEGG signaling pathways of the DEGs. GO = Gene Ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes.

To analyzing the functional interactions between DEGs, the PPI network constructed by STRING was visualized using Cytoscape (Fig. 3A), consisting of 33 nodes and 129 edges. Afterwards, the most key module of PPI network was extracted using MCODE (Fig. 3B), and 12 upregulated genes (SERPINE1, SPP1, LUM, POSTN, COL1A1, COL1A2, COL3A1, MMP1, MMP13, FN1, SPARC, VCAN) involved in the module were identified as hub genes based on the predefined criteria. Subsequently, the GO annotation demonstrated the biological processes of hub genes were mainly concentrated on extracellular matrix organization, collagen fibril organization, cellular response to growth factor stimulus, response to cytokine, regulation of cell-substrate adhesion, and regulation of cell migration (Fig. 3C). Finally, the coexpression network of hub genes obtained from cbioportal was constructed using Cytoscape (Fig. 4).

Figure 3

The PPI network construction and the key module analysis of the differentially expressed genes (DEGs). (A) The PPI network of DEGs was constructed using Cytoscape. (B) The most key module of the PPI network was extracted using MCODE. Upregulated genes are marked in red; downregulated genes are marked in blue. (C) The biological process analysis of hub genes. PPI = protein–protein interaction.

Figure 4

The coexpression network of hub genes. Red nodes represent hub genes, blue nodes represent the coexpression genes.

Survival analysis and verification of hub genes

Thorny difficulties were encountered in verifying hub genes due to limited data associated with HSCC in public databases. Referring to previous studies, the Head-Neck cancer dataset was used to evaluate hub genes in HSCC. In the present study, the expression differences of hub genes between tumors and normal tissues were verified using GEPIA. As shown in Fig. 5, the expression levels of hub genes in tumor samples were significantly elevated, which was consistent with our results. Moreover, survival analysis conducted by GEPIA demonstrated that SERPINE1 and SPP1 were closely related to poor prognosis of patients with HSCC and ESCC (Fig. 6), and the expression levels of SERPINE1 and SPP1 were found to be upregulated in multiple malignancies including breast carcinoma, esophageal carcinoma, kidney clear cell carcinoma, hepatocellular carcinoma, stomach adenocarcinoma, colorectal adenocarcinoma, thyroid carcinoma, and head and neck cancer (Fig. 7).

Figure 5

The expression differences of hub genes between tumors and normal tissues were verified using GEPIA online platform. ESCA = esophageal carcinoma; HNSC = head and neck squamous carcinoma.

Figure 6

Survival analysis of hub genes in patients with hypopharyngeal carcinoma and esophageal carcinoma. P-value < .05 was considered statistically significant.

Figure 7

The expression levels of SERPINE1 and SPP1 in multiple tumors. P-value significant codes: 0 ≤ ∗∗∗ < 0.001 ≤ ∗∗ < 0.01 ≤ ∗ < 0.05.

To illustrate the regulatory relationships between miRNA and hub genes, the interaction network of miRNAs and hub genes constructed by ENCORI was visualized using Cytoscape. As illustrated in Fig. 8, the interaction network consists of 11 hub genes and 116 miRNAs. After analyzing the network, all miRNAs and hub genes were ranked by degree score. Among them, COL1A2 (degree score = 32), SERPINE1 (degree score = 30), COL3A1 (degree score = 22), SPARC (degree score = 22), and MMP1 (degree score = 20) were confirmed as the top 5 interactive hub genes. Meanwhile, hsa-miR-29c-3p (degree score = 7), hsa-miR-29a-3p (degree score = 6), and hsa-miR-29b-3p (degree score = 6) were confirmed as the top 3 interactive miRNAs that target the most hub genes. Previously, Qiu et al[ reported that miR-29a/b could promote cell invasion and migration by inducing SPARC and COL3A1 gene expression in nasopharyngeal carcinoma. Similarly, our results also confirmed that miR-29a/b can regulate the expression of the SPARC and COL3A1. Therefore, the interaction network may contribute to understanding the molecular mechanisms involved in the carcinogenesis and progression of HSCC and ESCC.

Figure 8

Interaction network of hub genes and targeted miRNAs. Hub genes are presented in red circles, whereas targeted miRNAs are shown in blue circle.

Discussion

Hypopharyngeal and esophageal squamous cell carcinoma are the most common double primary tumors worldwide. Even with aggressive treatment, the prognosis remains poor. Recently, genomic studies have demonstrated that the genomic changes in ESCC are similar to those in HSCC based on the identified common risk factors.[ However, the common carcinogenic mechanisms remain unclear, as no relative study has been carried out. With the development of pan-cancer research and high throughput sequencing technology, the bioinformatic analysis may contribute to identifying the common DEGs and functional pathways involved in the carcinogenesis and progression of HSCC and ESCC.

In the study, 3 datasets were included to identify the common DEGs between HSCC, ESCC, and normal tissues to offset the false-positive rates in independent datasets analysis. As a result, a total of 43 DEGs were obtained, consisting of 25 upregulated genes and 18 downregulated genes. Subsequently, GO annotation, KEGG pathways enrichment analysis, and PPI network construction were conducted to explore interactions between the DEGs. The results indicate that the DEGs are mainly involved in the extracellular matrix organization, collagen metabolic process, extracellular matrix–receptor interaction, focal adhesion, and epidermis development. As reported, the extracellular matrix as the medium of cell communication exerts a significant impact on the carcinogenesis and progression of tumors, and activation of the extracellular matrix is also regarded as a landmark event for the formation of tumors.[ In addition, aberrant adhesion of tumor cells to the extracellular matrix also plays a critical role in tumor invasion and metastasis.

To conduct deeper research, 12 hub genes involved in the key module were identified based on the predefined criteria. After verification, the expression levels of hub genes were significantly elevated in tumor samples. Subsequently, we conduct survival analysis and found SERPINE1 and SPP1 were closely related to poor prognosis. SERPINE1 is known as the plasminogen activator inhibitor, which plays a crucial role in enhancing tumor cell migration and invasion through the PI3K-Akt pathway, promoting angiogenesis, protecting tumor cells from Fas/Fas-L mediated apoptosis.[ Meanwhile, SERPINE1 overexpression was proven to be strongly associated with poor prognosis in multiple cancers including breast cancer, fibrosarcoma, esophageal cancer, and head and neck cancer.[ At present, SERPINE1 has been established as a prognostic marker in patients with early lymph-node negative breast cancer.[ In addition, previous studies have demonstrated SERPINE1 was closely related to chemoradiotherapy resistance in head and neck cancer. One hand, the expression level of SERPINE1 is significantly upregulated after reactive oxygen exposure or irradiation, which activates hypoxia-related factors, thereby contributing to radiation resistance.[ On the other hand, the overexpression of SERPINE1 could protect tumor cells from inducing apoptosis after cisplatin treatment, which is mediated by PI3K-Akt pathway activation.[ In the context, the antitumor activity of several small molecules inhibitors targeting SERPINE1 is currently being evaluated. Among them, Tiplaxtinin is proven to block tumor cells’ growth and induce apoptosis in head and neck cancer.[ However, more preclinical and clinical trials are necessary to explore the application of the specific SERPINE1 inhibitors in the treatment of patients with HSCC and ESCC. SPP1, also known as osteopontin, is a secreted glycophosphoprotein. Kim et al[ found that the ectopic overexpression of SPP1 could activate ITGB1/FAK/AKT pathway, thereby enhancing the metastatic ability of head and neck cancer cells. Meanwhile, SPP1 overexpression was also reported to be involved in tumor proliferation, invasion, and angiogenesis in multiple malignancies including lung, breast, colorectal, and head and neck cancer. Moreover, our results demonstrated SPP1 was closely related to poor prognosis in patients with HSCC and ESCC.[ Overall, SERPINE1 and SPP1 play vital roles in the carcinogenesis and progression of HSCC and ESCC, which may present therapeutic targets and prognostic markers for patients with HSCC and ESCC in the future.

Although the other hub genes were not found to be directly associated with prognosis in our study, they are still proved to be involved in the carcinogenesis and progression of tumors. COL1A1 and COL1A2 encode the pro-α1 chain and pro-α2 chain of type I collagen, respectively, which is a key structural component of the extracellular matrix.[ Previous studies have demonstrated ESCC cells could produce COL1A1 endogenously,[ and miR-133a-3p could inhibit the ESCC cell proliferation, invasion, and migration by targeting COL1A1.[ However, the specific roles of COL1A2 in various tumors remain controversial. In bladder and colorectal cancer, COL1A2 was significantly downregulated and mainly plays an anticarcinogenic role in tumor development.[ However, in other malignancies such as ovarian cancer, pancreatic cancer, and head and neck cancer, overexpressed COL1A2 was found to promote tumor cell invasion and migration.[ Thus, the unique roles of COL1A2 in various cancers may be partially attributed to specific genetic characteristics of the different malignancies. Matrix metalloproteinases (MMPs) are zinc-dependent endopeptidase involved in the degradation of extracellular matrix and basement membrane, which were regarded as strong predictors of tumor metastasis.[ In the study, MMP1 and MMP13 were screened out as DEGs of HSCC and ESCC, implying that they may play important roles in tumor progression. Previous studies have demonstrated that MMP1 could facilitate tumor cells metastasize into the blood or lymphatic circulation by degrading interstitial collagen and vascular endothelium.[ In addition, it also supports tumor angiogenesis by activating protease-activated receptor-1.[ However, distinct from the function of MMP1 in angiogenesis, MMP13 not only promoted capillary tube formation but also induced vascular endothelial growth factor-A (VEGF-A) secretion from endothelial cells, which can indirectly stimulate tumor angiogenesis.[

FN1, an extracellular matrix glycoprotein, mainly mediates the interaction between tumor cells and extracellular matrix.[ Previous studies have reported that FN1 could activate the PI3K/Akt pathway to stimulate tumor cell proliferation and invasion through binding to α5β1 integrin receptors.[ Meanwhile, POSTN has also been proven to be able to bind to multiple integrin receptors including αvβ3, αvβ5, and α6β4, thereby regulating the intracellular PI3K/Akt signaling pathway.[ Abnormal activation of PI3K/Akt signaling pathway could directly mediate the epithelial-mesenchymal transition, which is crucial for tumors to obtain malignant properties.[ Thus, targeting PI3K/Akt signaling pathway may be a promising strategy for anti-cancer therapy. VCAN belongs to the large chondroitin sulfate proteoglycans family and plays a vital role in the formation of tumor-specific extracellular matrices. As reported, VCAN could induce tumorigenesis by inhibiting TNF signaling-mediated apoptosis and promote tumor invasion and metastasis through inducing MMPs expression.[ SPARC, also termed osteonectin, mainly regulates cell–matrix interactions and cell adhesion.[ Che et al[ have demonstrated that overexpressed SPARC is closely related to lymph node metastasis and distant metastasis in ESCC patients. Meanwhile, as a vital downstream target of SPARC, the expression level of COL3A1 has been found to be positively regulated by SPARC. Moreover, consistent with previous studies,[ our results also confirmed that miR-29a/b can regulate the expression of the SPARC and COL3A1. As for LUM, conflicting data have been reported in the previous studies with regard to its exact role in tumor progression.[ For instance, LUM could modulate the expression of MMP-9 and MMP14 to inhibit tumor cell migration in breast cancer, which mainly acts as the anticancer effector.[ However, in pancreatic carcinoma, LUM overexpression plays a carcinogenic effect by promoting cell invasion and metastasis.[ Thus, the exact role of LUM in ESCC and HSCC deserves further exploration.

In conclusion, the present study found 43 common DEGs between ESCC, HSCC, and normal samples using bioinformatics analysis, which were mainly involved in the extracellular matrix–receptor interaction, collagen metabolic, epidermis development, cell adhesion, and PI3K/Akt signaling pathways. Among DEGs, 12 genes (SERPINE1, SPP1, LUM, POSTN, COL1A1, COL1A2, COL3A1, MMP1, MMP13, FN1, SPARC, VCAN) were identified as hub genes, which may serve as diagnostic and therapeutic targets in ESCC and HSCC. In addition, SERPINE1 and SPP1 were proven to be closely related to poor prognosis, which may be potential prognostic biomarkers. Meanwhile, the interaction network of miRNAs and hub genes illustrates the regulatory relationships of the hub genes and miRNA, and hsa-miR-29c-3p, hsa-miR-29a-3p, and hsa-miR-29b-3p were confirmed as the top 3 interactive miRNAs that target the most hub genes. These results provide important ideas for a comprehensive understanding of cancer characteristics, however, further studies are needed to validate the current findings and elucidate the specific molecular mechanisms of these genes in ESCC and HSCC.

Author contributions

Conceptualization: Tao Zhang.

Data curation: Rui Zhou, Denghua Liu.

Formal analysis: Rui Zhou, Denghua Liu.

Software: Rui Zhou, Denghua Liu.

Supervision: Jing Zhu, Tao Zhang.

Validation: Rui Zhou, Denghua Liu. Visualization: Rui Zhou, Denghua Liu.

Writing – original draft: Rui Zhou.

Writing – review & editing: Tao Zhang.