Identification of a transcription factor-microRNA network in esophageal adenocarcinoma through bioinformatics analysis and validation through qRT-PCR

Introduction

Esophageal carcinoma has a five-year survival rate of 12–20% in Western populations^1,2 and induces the deaths of more than 400,000 people worldwide annually.³ The incidence of esophageal adenocarcinoma (EAC), one major histologic subtype of esophageal carcinoma, has increased by approximately 6-fold in the past 40 years.⁴ EAC is a form of cancer with a rapidly increasing incidence in Western countries.^5,6 In the last few decades, the prognosis of patients with EAC has mildly improved, but it remains worse than that of most other types of cancer.^7–9 However, the molecular mechanisms giving rise to EAC remain unclear. Therefore, a thorough exploration of EAC pathogenesis is urgently needed to improve EAC diagnosis, treatment and prognosis.

Transcription factors (TFs) are DNA-binding proteins that can function as tumor suppressors or oncogenes.¹⁰ TFs play significant roles in the regulation of gene expression and can induce avoidance of apoptosis and uncontrolled growth.¹¹ microRNAs (miRNAs) are short, noncoding RNAs consisting of 18–25 nucleotides that regulate the translation of mRNAs.¹² Mature miRNAs can recognize and bind to the 3ʹ untranslated region of mRNAs and regulate translation at the posttranscriptional level through repression or degradation of the targeted mRNAs.¹³ Numerous studies have revealed that TFs regulate gene expression by interacting with miRNAs. For instance, Chandra et al found that two TFs, ZEB1 and ZEB2, could contribute to epithelial-mesenchymal transition (EMT) through interaction with miR-101–1.¹⁴ In addition, Hua et al demonstrated that ETS1 could inhibit cell proliferation, invasion and migration by regulating miR-139-5p in hepatocellular carcinoma (HCC).¹⁵ In addition to miRNAs, lncRNAs and circRNAs were also proven to interact with TFs. Sun et al reported that NRG-1-ICD could be regulated by circACTA2 via the NRG-1-ICD/circACTA2/miR-548f-5p axis in vascular smooth muscle cells (VSMCs).¹⁶ Additionally, Guo et al found that KIF25-AS1, LINC01355 and AC092171.2 were regulated by POLR2A in ovarian cancer tissues.¹⁷

Recently, bioinformatics analysis has gradually been utilized to explore the mechanisms of cancer initiation and progression to identify important cancer-related genes, noncoding RNAs and TFs for further experimental verification. In this study, dysregulated genes and miRNAs in EAC were identified from the Gene Expression Omnibus (GEO) database. Subsequently, to understand the underlying molecular mechanisms in EAC, we used a public database to find the TFs that regulate the dysregulated genes. miRNAs were further predicted according to the TFs that targeted the DEGs, and interactions between miRNAs and TFs were predicted with miRTarBase. Based on our findings, TF-miRNA-mRNA networks were constructed to reveal a transcriptional regulation model and identify key genes and miRNAs associated with EAC. Functional enrichment analysis and protein-protein interaction (PPI) network analysis were also performed, and we used the Cancer Genome Atlas (TCGA) database to identify correlations between the expression of key genes and miRNAs in the TF-miRNA-mRNA network and overall survival. Finally, quantitative real-time PCR (qRT-PCR) was utilized to validate the results.

Methods

Collection of datasets

A flowchart of the study design is shown in Figure 1. Both mRNA and miRNA expression profiles were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/gds/). The GSE92396, GSE13898, GSE26886 and GSE1420 datasets were used for mRNA, while the GSE16456 dataset was used for miRNA. All of the datasets were selected based on the following criteria: (1) there were no fewer than six samples, (2) the samples were from human esophageal tissue, and (3) the datasets included case-control groups. Overall, 111 EAC and 80 normal esophagus (NE) samples were analyzed.

Figure 1 Flowchart of the bioinformatics analysis. Abbreviations: DEGs, differentially expressed genes; DEMs, differentially expressed miRNAs; miRNA, microRNA; PPI, protein-protein interaction; TCGA, the Cancer Genome Atlas; TF, transcription factor; TRED, Transcriptional Regulatory Element Database; TRRUST, Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining.

Aberrantly expressed miRNA and DEG analysis

The limma algorithm package¹⁸ in R software was utilized to identify DEGs and aberrantly expressed miRNAs in EAC tissues compared with NE tissues. DEGs were identified based on the following criteria: (1) they had a |log₂FC|≥1 and a P-value <0.05, (2) they were shared among the four datasets, and (3) they showed consistent trends in all four datasets. The Benjamini and Hochberg false discovery rate (FDR) and adjusted P-values (adj.p) were applied to provide a balance between the probability of obtaining false positives and the probability of discovering statistically significant genes. The cutoff criteria for aberrantly expressed miRNAs were a |log₂FC|≥1 and a P-value <0.05. Finally, an UpSet plot was constructed to visualize the DEGs in the different mRNA datasets.

Prediction of TFs and miRNAs and construction of a TF-miRNA-mRNA network

The Transcriptional Regulatory Relationships Unraveled (TRRUST, http://www.grnpedia.org/trrust/) database presents the relationships between TFs and genes based on the existing literature,¹⁹ and the Transcriptional Regulatory Element Database (TRED, https://cb.utdallas.edu/cgibin/TRED/tred.cgi?process=home) includes 36 TFs related to cancer and TF-target gene regulatory networks.²⁰ TFs that not only were cancer-related but also targeted DEGs and the TF-targeted DEGs (key genes) were identified by a comprehensive analysis of the above two databases. Subsequently, miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php) was used to predict the miRNAs that regulate the key genes and to predict the interactions between the miRNAs and TFs.²¹ We identified overlapping miRNAs by intersecting the dysregulated miRNAs from the GSE16456 dataset and the miRNAs obtained from miRTarBase. Finally, a TF-miRNA-mRNA network was constructed to show the potential molecular mechanisms of EAC initiation and progression.

Functional enrichment and pathway analysis

Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were conducted to analyze the functions of the key genes. Then, the Database for Annotation Visualization and Integrated Discovery (DAVID: http://david.abcc.ncifcrf.gov/) was used to perform GO and KEGG pathway analysis of key genes.²² The cutoff criterion was a P-value <0.05. The GOplot R package was used to visualize the results of GO analysis, while the results of KEGG pathway analysis were visualized with Cytoscape software.

PPI network construction

To identify relevant PPIs, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://www.string-db.org/) was utilized to construct a PPI network for the key genes.²³ The PPIs were analyzed using the combined score from the STRING database, and a combined score >0.4 was used as the cutoff criterion. Additionally, Cytoscape software was used to visualize the PPI network. The modules in the PPI network were extracted using the Cluster Finding algorithm in MCODE with a max depth of 100, a node score cutoff of 0.2, a degree cutoff of 2, and a k-core of 2.²⁴

Survival analysis validation

TCGA, which collates information from the National Human Genome Research Institute (NHGRI) and the National Cancer Institute (NCI), has generated multidimensional and comprehensive maps for multiple types of cancers. To verify the results obtained thus far, data for 89 EAC tissue samples and the clinical data of the patients were downloaded from TCGA, and the prognostic value of the key genes and miRNAs was identified by Kaplan-Meier analysis. The survival package in R software was used for survival analysis. A P-value <0.05 was set as the cutoff criterion.

Total RNA extraction and qRT-PCR validation

A total of 20 EAC tissues and adjacent noncancerous tissues were collected from Qingdao University Hospital. This study was approved by the Ethical Committee of Qingdao University Hospital and was conducted in accordance with the Declaration of Helsinki, and each patient offered written informed consent. All tissues were immediately frozen in liquid nitrogen after resection and stored at −80 °C.

Tissue sample RNA was extracted using TRIzol (Invitrogen, CA, USA) according to the manufacturer’s protocol. Total RNA including miRNAs was reverse transcribed into cDNA with a miRNA First Strand cDNA Synthesis Kit (Sangon, China) and a 1st strand cDNA Synthesis Kit (Takara, China). miRNA and mRNA expression was then examined by qRT-PCR with a MicroRNAs qPCR Kit (Sangon, China), a Biosystems 7500 Fast Real-Time PCR System, and TransStart® Top Green qPCR SuperMix (TransGen Biotech, China). miRNA and mRNA expression was normalized to U6/GAPDH expression. The 2^−∆∆Ct method was used to calculate the levels of miRNA and mRNA expression. The primers were designed using Primer 5.0 software, and the primer sequences are shown in Table 1. GraphPad Prism version 7.0 was used to analyze the qRT-PCR data, and comparisons between two groups were performed using Student’s t-test, with a P-value <0.05 considered significant.

Table 1 Primer sequences of qRT-PCR

Results

Microarray datasets from the GEO

miRNA profiles were obtained from the GSE16456 dataset, while mRNA expression profiles were obtained from the GSE13898, GSE26886, GSE92396 and GSE1420 datasets. Twenty Barrett’s esophagus, 19 normal and 21 EAC tissues were included in the GSE26886 dataset, while 15 Barrett’s esophagus, 28 normal and 64 EAC tissues were included in the GSE13898 dataset. The GSE92396 dataset included a total of 9 normal tissues and 12 EAC tissues, while the GSE1420 dataset included 8 Barrett’s esophagus, 8 normal and 8 EAC tissues. The miRNA dataset GSE16456 included 10 Barrett’s esophagus, 16 normal and 6 EAC tissues. All relevant details are shown in Table 2.

Table 2 Information of expression profiling microarray datasets

Analysis of DEGs and differentially expressed miRNAS

To identify DEGs, the GEO mRNA expression profile datasets were analyzed. A total of 294 DEGs were obtained through examining genes that showed the same expression trends among the mRNA datasets (Figure 2). Among the DEGs, 113 genes were upregulated in EAC, while the remaining 181 genes were downregulated (|log₂FC|≥1, adjusted P-value <0.05)

Figure 2 Identification of overlapping DEGs in the mRNA expression profile datasets GSE1420, GSE26886, GSE13898, and GSE92396.Abbreviation: DEG, differentially expressed gene.

A total of 61 DEMs were identified between the EAC and NE tissues in the GSE16456 dataset. Among them, 37 miRNAs were significantly downregulated in EAC, while the other 24 miRNAs were upregulated (|log₂FC|≥1, adjusted P-value <0.05).

Construction of the TF-miRNA-mRNA network

Upon analysis of 294 DEGs for potential TFs in the TRRUST database, which is a database of transcriptional regulatory networks, 46 TFs and 70 TF-targeted genes were identified. Then, we cross-referenced the above results with the 36 cancer-related TF families and TF-target gene regulatory networks in the TRED database and ultimately identified 16 TFs and 41 key genes in EAC. Subsequently, miRTarBase was used to identify miRNAs predicted to be linked to the 41 key genes and to predict interactions between the TFs and miRNAs. Furthermore, 12 overlapping miRNAs were identified from among the DEMs in the GSE16456 dataset and the miRNAs predicted from miRTarBase. A TF-miRNA-mRNA interaction network was constructed based on the TRED, TRRUST and miRTarBase analyses. Finally, TP53-(miR-125b)-ID2 and JUN-(miR-30a)-IL1A were identified as possibly playing significant roles in EAC. Cytoscape software was used to visualize the TF-miRNA-mRNA network in EAC (Figure 3 and Table 3).

Table 3 Regulatory relationships among TFs, key genes, and miRNAs

Figure 3 TF-miRNA-mRNA network. The blue rhombuses indicate TFs. The blue circles indicate downregulated genes, and the red circles indicate upregulated genes. The blue squares represent downregulated miRNAs, and the red squares represent upregulated miRNAs. Abbreviations: DEG, differentially expressed gene; TF, transcription factor; miRNA, microRNA.

Functional and pathway analysis of key genes

To further study the functional roles of key genes in the pathogenesis of EAC, GO and KEGG pathway enrichment analyses were performed. In the GO biological process (BP) group, the key genes were mainly enriched for terms related to the inflammatory response, extracellular matrix disassembly, positive regulation of I-kappaB kinase/NF-kappaB signaling, and epidermis development, among others. In the cellular component (CC) group, the key genes were mainly enriched for the extracellular space term, while in the molecular function (MF) group, the key genes were mainly enriched for the cytokine activity term (Figure 4 and Table 4).

Table 4 GO terms of key genes

Figure 4 GO analysis results for the key genes. Abbreviation: GO, Gene ontology.

The results of KEGG pathway analysis suggested that the key genes were significantly enriched in the TNF signaling pathway, rheumatoid arthritis, cytokine-cytokine receptor interaction, transcriptional misregulation in cancer and pathways in cancer (Figure 5 and Table 5). In particular, the KEGG results further verified that the regulation of the key genes was associated with TFs.

Table 5 KEGG pathways of key genes

Figure 5 KEGG pathway analysis results for the key genes. Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI analysis

PPI networks, which are used to analyze the interaction of proteins, were constructed in the STRING database with a combined score >0.4 as the cutoff criterion. Then, we used the MCODE plug-in to identify functional modules. One module consisting of 10 nodes and 31 edges was found and included IGFBP3, MMP3, SPP1, MMP1, PLAUR, VCAN, IL1A, MMP9, ECM1, and CXCL1 (Figure 6).

Figure 6 PPI network of the key genes. The yellow rectangles represent the genes in the functional modules of the PPI network. Abbreviation: PPI, protein-protein interaction.

Survival analysis validation in TCGA

The clinical data and RNA expression data of patients from TCGA were used to validate the key genes and miRNAs. The survival R package was used to analyze patient survival, and three key genes (CCL20, IL1A, and ABCC3) and two miRNAs (hsa-miR-23b and hsa-miR-191) were identified as being significantly associated with prognosis in EAC (Figure 7). More importantly, one miRNA-gene pair (hsa-miR-30a-5p and IL1A) was identified as being related to prognosis in EAC (Figure 8).

Figure 7 Prognostic value of ABCC3, CCL20, IL1A, hsa-miR-23b, and hsa-miR-191 in EAC patients from TCGA datasets (A–E). Abbreviations: EAC, esophageal adenocarcinoma; TCGA, the Cancer Genome Atlas.

Figure 8 Prognostic value of hsa-miR-30a-5p and IL1A in EAC patients from TCGA datasets. Abbreviations: EAC, esophageal adenocarcinoma; TCGA, the Cancer Genome Atlas.

Validation of miRNAs and key genes by qRT-PCR

qRT-PCR was performed to validate the expression of the key genes and miRNAs. We found that three key genes (CCL20, IL1A, and ABCC3) and two miRNAs (hsa-miR-23b and hsa-miR-191) were significantly dysregulated in EAC. Of these genes and miRNAs, IL1A and miRNA-23b were significantly downregulated, while CCL20, ABCC3 and hsa-miR-191 were significantly upregulated in EAC (P<0.001). The outcome revealed that the expression of the miRNAs and key genes was consistent with the results of our survival analysis, indicating that the results were reproducible and reliable (Figure 9).

Figure 9 qRT-PCR validation of the expression of ABCC3, CCL20, IL1A, hsa-miR-191 and hsa-miR-23b in EAC patients (A–E). Notes: **P<0.01, ***P<0.001, ****P<0.0001. Abbreviations: EAC, esophageal adenocarcinoma; qRT-PCR, quantitative real-time PCR.

Discussion

EAC is one of the most common types of esophageal cancer, and the five-year survival rate is 12–20% in Western populations.^1,2 However, the precise pathological mechanism of EAC remains unclear. Therefore, exploring the mechanisms of EAC progression is essential to improve survival and prevent progression. In recent years, microarray technology has been widely used to investigate gene expression in multiple tumor types, providing a novel approach for identifying important genes and leading to further study on tumor initiation and progression.

In this study, 294 DEGs were identified in four mRNA datasets, including 181 downregulated and 113 upregulated genes. Then, public databases were used to predict the TFs that targeted the DEGs. Generally, not all DEGs correspond to TFs. Through a series of comprehensive analyses, 16 TFs (TP53, NFKB1, JUND, RELA, POU5F1, CEBPA, ETS1, JUN, FOS, ETS2, TFAP2A, FLI1, ETV4, TFAP2C, CEBPD and RARA) and 41 key genes were identified to be related to cancer progression. In recent years, many studies have stated that TFs regulate gene expression by interacting with miRNAs. Thus, we further studied the correlations between the 41 key genes and miRNAs and the interactions between the TFs and miRNAs with the miRTarBase database. Subsequently, a TF-miRNA-target gene network was constructed to study the potential molecular mechanisms, and we found that TP53-(miR-125b)-ID2 and JUN-(miR-30a)-IL1A might play significant roles in EAC. Once the network was constructed, enrichment analysis and PPI analysis were performed to understand the underlying functions of these key genes. Intriguingly, the key genes were enriched in the TNF signaling pathway, transcriptional misregulation in cancer and pathways in cancer, which further verified that the key genes were associated with misregulation of TFs leading to tumorigenesis. Furthermore, clinical data and RNA-Seq data from TCGA were used to analyze the correlations of the key genes and miRNAs with patient survival. The analysis predicted ABCC3, IL1A, and CCL20 as independent factors associated with EAC prognosis. ABCC3 is a member of the ATP-binding cassette (ABC) family that is associated with the prognosis of several cancers.²⁵ Zhao et al revealed that overexpression of ABCC3 is related to TNM stage and poor prognosis in non-small cell lung cancer.²⁶ Liu et al showed that increased expression of ABCC3 is correlated with worse overall survival in urinary bladder cancer patients.²⁷ CCL20 belongs to the chemotactic cytokine family, which is correlated with the occurrence, invasion, and metastasis of cancer.^28–30 A recent study indicated that the expression of CCL20 is increased in esophageal squamous cell carcinoma and can help classify patients into low-risk, medium-risk and high-risk groups.³¹ In addition, the prognostic value of CCL20 in multiple tumors has been evaluated. For example, it has been reported that elevated expression of CCL20 is inversely associated with prognosis in smoker patients with lung cancer, providing a potential therapeutic target for this deadly disease.³² IL1A is a member of the interleukin-1 family of cytokines that have multiple functions in cancer and inflammation; when expressed on the cytomembrane, IL1A heightens the immunity of antitumor cells, as indicated by tumor regression.³³ However, IL1A has received less attention as a potential tumor suppressor gene than ABCC3 and CCL20. Although the prognostic value of ABCC3, CCL20 and IL1A in EAC has not been reported in previous studies, the importance of these three genes should not be underestimated. In addition, we discovered that hsa-miR-23b and hsa-miR-191 were significantly associated with patient survival in EAC. Many previous studies have reported miR-23b as a tumor suppressor that combats the occurrence and progression of different cancers. The expression of miR-23b is reduced in clear cell renal cell carcinoma, and decreased miR-23b expression predicts poor survival and increased risk of disease progression.³⁴ Downregulation of miR-23b is significantly correlated with poor prognosis in epithelial ovarian cancer (EOC) patients and may be a potential prognostic marker for EOC.³⁵ However, the prognostic value of miR-23b and miR-191 in EAC has not been validated in previous studies. miR-191 has been reported to participate in the initiation and progression of multiple cancers, including pancreatic cancer, colon cancer, lung cancer, and gastric cancer. Overexpression of miR-191 has been discovered in four gastric cancer cell lines and in gastric tumor tissues, indicating miR-191 that may be a novel biomarker for gastric cancer.³⁶ Additionally, inhibition of miR-191 decreases HCC cell proliferation and tumor growth,³⁷ and miR-191 has been found to be closely correlated with TNM stage, metastasis and poor prognosis in pancreatic cancer.³⁸ In our study, a miRNA-mRNA pair (hsa-miR-30a and IL1A) was found to be a prognostic marker for EAC. However, the prognostic value of hsa-miR-23b, hsa-miR-191, and the hsa-miR-30/IL1A pair in EAC has not previously been evaluated. Thus, qRT-PCR was performed to verify the results of bioinformatics analysis and revealed that the expression of the miRNAs and key genes was consistent with our analysis.

Conclusion

In this study, we constructed a TF-miRNA-mRNA network to analyze the potential molecular mechanisms of EAC and identified 16 TFs, 41 key genes, 12 miRNAs and the TF-miRNA-mRNA sets TP53-(miR-125b)-ID2 and JUN-(miR-30a)-IL1A. Among the identified molecules, CCL20, IL1A, ABCC3, miR-23b and miR-191 were determined to be potential independent prognostic factors via survival analysis. qRT-PCR was used to validate the results. In addition, a miRNA-mRNA pair (hsa-miR-30a-5p and IL1A) was determined to be significantly associated with patient survival through integrative survival analysis; this finding might provide a novel direction for further experiments, and the pair could be a potential target for precision treatment. The dysregulated genes identified in this study may be promising biomarkers for EAC diagnosis, treatment and prognosis and can provide a basis for follow-up experiments.