Bioinformatic analyses of microRNA-targeted genes and microarray-identified genes correlated with Barrett's esophagus

ABSTRACT

Barrett's esophagus (BE) is defined as a metaplasia condition in the distal esophagus, in which the native squamous epithelium lining is replaced by a columnar epithelium with or without intestinal metaplasia. It is commonly accepted that BE is a precancerous lesion for esophageal adenocarcinoma. The aim of this study was to investigate the aberrant microRNAs (miRNAs) and differentially expressed genes (DEGs) associated with BE based on online microarray datasets. One miRNA and five gene expression profiling datasets were retrieved from the Gene Expression Omnibus Database. Aberrant microRNAs and DEGs were obtained using R/Bioconductor statistical analysis language and software. 23 dysregulated miRNAs and 632 DEGs demonstrating consistent expression tendencies in the five gene microarrays were identified in BE. Moreover, 1962 target genes of aberrant miRNAs were predicted using three bioinformatic tools, namely TargetScan, RNA22-HSA and miRDB. Ultimately, 93 target DEGs were obtained, after which functional annotation was performed on DAVID Bioinformatics Resources. Among Gene Ontology (GO) biological processes, digestive tract development and epithelial cell differentiation have demonstrated significant associations with BE pathogenesis. In addition, analysis of the KEGG pathways has revealed associations with cancer. To enable further study, one miRNA-target DEGs regulatory network was constructed using Cytoscape. 6 target DEGs demonstrated higher-degree distributions in the network, and ROC analysis indicated that FNDC3B may be the best potential biomarker for BE diagnosis. The data presented herein may provide new perspectives for exploring BE pathogenesis and may offer hits with regard to potential biomarkers in BE diagnosis, prediction and therapeutic evaluation.

KEYWORDS:

• Barrett's esophagus
• bioinformatic analyses
• differentially expressed genes
• microarray
• microRNA

Introduction

Barrett's esophagus (BE) is defined as a metaplasia condition in the distal esophagus, in which the native squamous epithelium lining is replaced by a columnar epithelium with or without intestinal metaplasia [1]. It is a commonly accepted risk factor and important precancerous lesion for esophageal adenocarcinoma (EAC). In recent years, the incidences of BE and EAC have risen consistently, and the risk of BE progressing to EAC has increased by 30–125 times in the general population [2]. It is noteworthy that most studies regarding BE pertain to its carcinogenesis, the underlying pathogenesis of which remains unclear. Therefore, further study is required to clarify BE pathogenesis [3].

MicroRNAs (miRNAs) are a class of endogenous small (18-25 nucleotides) non-coding RNAs derived from long primary transcripts and are found in eukaryotic organisms [4]. Mature miRNAs can identify and bind to target mRNA 3'-untranslated regions (3'-UTR) by means of complementary base pairing, by which the degradation or translational repression of mRNA targets at the post-transcriptional level are regulated [5,6]. Known as post-transcriptional regulators, miRNAs have received substantial attention from researches, and many pertinent studies have already been performed. Recent studies have indicated that miRNAs may be involved in regulating various biological activities, such as development, organ formation, cell proliferation, differentiation, migration and apoptosis [7–9]. Furthermore, some miRNAs have been reported to function as oncogenes or tumor suppressor genes [10–12]. Interestingly, miRNAs present several advantages. First, the number of human miRNAs is significantly smaller than that of mRNAs. In addition, as small non-coding RNA fragments, miRNAs exhibit greater stabilities than mRNAs. Finally, one miRNA could regulate translation of hundreds of mRNAs, and one mRNA could be regulated by a number of miRNAs. Therefore, expression analyses of miRNAs may improve risk stratification analyses for different types of disorders. Microarrays have been widely used to explore expression genes, including miRNAs expression profiling. A series of studies [11,13–17] have indicated that microarrays provide a more effective means of identifying and verifying potential biomarkers in BE at the molecular level.

A number of studies have used microarrays for expression profiling of miRNAs in BE patients [15,18–21]. Compared with healthy controls, these differentially expressed miRNAs may regulate different mRNA targets and influence the pathogenesis of BE. Hence, miRNAs have the potential to become new biomarkers for BE diagnosis, prediction and therapeutic evaluation [10,11,15,18,20,22]. However, the functional roles of these aberrant miRNAs in BE pathogenesis remain unclear and their mRNA targets are poorly studied.

Focusing on BE pathogenesis, the present study investigates the potential miRNAs that are expressed differently in miRNAs expression profiling and differentially expressed genes (DEGs) in five mRNA expression profiling datasets that correlate with BE as determined using online microarrays. Bioinformatic analyses were performed to explore the enrichment in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms. For further study, the interaction network between aberrant miRNAs and target DEGs was constructed, and biomarkers were investigated.

Results

Expression profiling datasets retrieved from GEO Datasets of NCBI

All the microarray samples considered in this study were fresh frozen tissues from human esophagi. First, we applied quality control to the expression microarray chips with QC function in R/Bioconductor to guarantee the quality of our study. In this study, the miRNAs profiling dataset, which included results from 10 BE and 10 NE samples, was GSE24839 using the GPL8871 Platform, which was generated by the University of Padova. The mRNA expression profiling datasets were GSE13083, GSE13898, GSE26886, GSE34619 and GSE39491, respectively. GSE13083 originated from the University of Pennsylvania, which used the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array Platform, and 7 BE and 7 NE samples. GSE26886, generated from Charite-Berlin on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array Platform, included 20 BE and 19 NE samples originally. However, the sample numbers were reduced to 8 BE and 10 NE after quality control. All pertinent details are listed in Table 1.

Table 1. Information of expression profiling microarray datasets.

						Numbers of raw samples		Final Numbers of samples
Accession number	Platform	Organism	Experiment type	Organization	Status	BE	NE	BE	NE
GSE24839	GPL8871	Homo sapiens	Non-coding RNA profiling by array	University of Padova	Public on Jan 01, 2011	10	10	10	10
GSE13083	GPL96	Homo sapiens	Expression profiling by array	University of Pennsylvania	Public on Oct 08, 2008	7	7	7	7
GSE13898	GPL6102	Homo sapiens	Expression profiling by array	University of Texas	Public on Jan 04, 2011	2	28	2	28
GSE26886	GPL570	Homo sapiens	Expression profiling by array	Charite-Berlin	Public on Oct 18, 2013	20	19	8	10
GSE34619	GPL6244	Homo sapiens	Expression profiling by array	MRC Cancer Cell Unit	Public on Mar 30, 2012	10	8	10	8
GSE39491	GPL571	Homo sapiens	Expression profiling by array	National Cancer Institute, National Institutes of Health	Public on Apr 15, 2014	40	40	28	30

Identification of aberrant miRNAs and DEGs in Barrett's esophagus

A total of 23 miRNAs were identified as dysregulated in BE tissues compared to those in NE mucosa (|FC|> 2, P values < 0.05, see Supplementary Table 1 for details). Among these miRNAs, 7 miRNAs were significantly up-regulated (hsa-miR-215, hsa-miR-192, hsa-miR-194, hsa-miR-10a, hsa-miR-145*, hsa-miR-200a and hsa-miR-217), while the expression levels of the other 16 miRNAs were decreased conversely in BE. The top 6 dysregulated miRNAs were hsa-miR-18a*, hsa-miR-205, hsa-miR-215, hsa-miR-203, hsa-miR-192 and hsa-miR-194, for which the values of Log₂|FC| were -5.06, -3.92, 3.83, -3.39, 2.21 and 2.13, respectively. Bi-cluster analysis was conducted, and the data showed that most of the BE tissues investigated could be separated well from NE tissues (Figure 1).

Figure 1. Bi-cluster analysis of differentially expressed miRNAs in Barrett's esophagus tissue Each row represents one of the dysregulated miRNA, and each column represents a tissue sample from Barrett's esophagus and normal esophagus mucosa. Column “BE” represents Barrett's esophagus tissue, while column “NE” represents tissue from normal esophagus. The heat map was constructed based on miRNAs expression value in each sample.

For analyses of the mRNA expression profiling datasets, the limma algorithm was utilized to obtain DEGs in BE, and five lists of DEGs were generated from the above five datasets (P values < 0.05). Then, we applied statistical methods to acquire a final list of common DEGs, and 632 DEGs were identified in all five microarrays (Supplementary Table 2). The gene numbers of each dataset and the intersections among these datasets are shown in Figure 2(A). Of the 632 commonly expressed DEGs, 286 genes were up-regulated and 346 genes were down-regulated in BE tissues (Figure 2(B)). Moreover, each gene of the 632 common DEGs displayed the same expression tendency in all five datasets.

Figure 2. Exhibition of differentially expressed genes among the five mRNA expression profiling datasets in Barrett's esophagus; (A) Venn diagram of individual dataset comparisons among the five mRNA expression profiling datasets. The venn diagram was drawn based on each gene expression profiling data from the five datasets. Each circle depicts the number of expressed genes. Overlapping expressed genes shared more than one dataset are represented in the areas of intersection among corresponding circles. The 632 genes in the center of the venn diagram represents the genes commonly expressed in all the five datasets; (B) Scatter diagram of differentially expressed genes among the five mRNA expression profiling datasets. The scatter diagram was constructed based on the Log₂(FC) values of the 632 commonly expressed DEGs derived from five datasets. These five datasets were GSE39491, GSE26886, GSE12083, GSE13898 and GSE34619 as the diagram shows.

Bioinformatic analyses of interactions between the potential targets of aberrant miRNAs and DEGs

MiRNAs could interact with their mRNA targets, resulting in post-transcriptional suppression of their target genes. Three main bioinformatic tools (TargetScan, RNA22-HSA and miRDB) were used to predict the potential candidate targets of the 23 dysregulated miRNAs. The target genes for each aberrant miRNA were considered common genes identified using all three bioinformatic tools. Then, we determined the intersection between the target genes and common DEGs in the five mRNA expression profiling datasets (Figure 3(A)). There were 93 common target DEGs in this intersection.

Figure 3. Bioinformatic analysis of predicted targets of aberrant miRNAs and DEGs in Barrett's esophagus; (A) Venn diagram of comparison between predicted targets of aberrant miRNAs and DEGs in the five mRNA expression profiling datasets. Each circle depicts the number of expressed genes. Overlapping expressed genes shared more than one dataset is represented in the area of intersection between 2 circles. The 93 genes in the center of the venn diagram represent the target genes of aberrant miRNAs, which are also commonly expressed in all the five mRNA expression profiling datasets; (B) Regulatory network of aberrant miRNAs and target common DEGs in Barrett's esophagus. Triangles represent aberrant miRNAs (yellow for up-regulated miRNAs and blue for down-regulated miRNAs), and circles represent target common DEGs (red for up-regulated DEGs and green for down-regulated DEGs). The direction of arrows is from dysregulated miRNAs to target DEGs.

To explore the biological functions of the 93 common target DEGs with regard to the regulation of BE pathogenesis, DAVID functional annotation was performed to obtain GO and KEGG pathway enrichment data. First, we uploaded the list containing these 93 DEGs for functional annotation. Then, data pertaining to the GO and KEGG pathway was collected as statistically significant enrichment (P values < 0.05) for further analyses (Tables 2 and 3). The information we obtained regarding biological function and pathway annotation indicated that the 93 DEGs were enriched in 5 Gene Ontology_Biological Process (GO_BP), 4 Gene Ontology_Cellular Component (GO_CC) and 1 Gene Ontology_Molecular Funtion (GO_MF) terms as well as in 5 KEGG pathway terms. Among these GO biological processes, digestive tract development and epithelial cell differentiation were found to be significantly related to BE pathogenesis. Genes included in the digestive tract development process were CLDN18, TRPS1 and COL3A1, with a fold enrichment value of 18.24. Genes enriched in the epithelial cell differentiation process were TAGLN, CNN3 and ANXA4, with a fold enrichment value of 10.09. Moreover, analysis of the enriched KEGG pathway revealed that pathways in cancer were involved with a fold enrichment of 2.90. The related genes were GNAI3, FGFR3, ROCK2, RASGRP1, TGFA, PLCB1 and DAPK1. The predominant, related pathways were MAPK and VEGF signaling pathways, which could regulate sustained angiogenesis and cell proliferation [23–26].

Table 2. GO terms of the commonly expressed 93 DEGs in five mRNA expression profiling datasets.

	GO terms	P	Fold enrichment	Gene symbols
Biological process	digestive tract development	0.01	18.24	CLDN18, TRPS1, COL3A1
	cellular response to starvation	0.02	12.82	CADPS2, DSC2, RRAGD
	epithelial cell differentiation	0.03	10.09	TAGLN, CNN3, ANXA4
	negative regulation of inflammatory response	0.04	8.95	TYRO3, PTGER4, RORA
	intracellular signal transduction	0.02	3.69	ARHGEF4, TUFT1, WNK1, PRKCH, PLCB1, DAPK1
Cellular component	mast cell granule	0.04	43.49	RASGRP1, CXCR2
	intracellular membrane-bounded organelle	0.02	6.50	HSD17B11, AP2B1, VCAN, AHCYL2
	focal adhesion	0.03	3.57	RND3, DAB2, CNN3, LASP1, EPHA2, SLC9A1
	intracellular	0.03	2.61	ARHGEF4, RND3, TYRO3, PHACTR4, TUFT1, MAP3K9, PRKCH, EPHA2
Molecular function	RNA polymerase II regulatory region sequence-specific DNA binding	0.00	7.65	HLF, KAT2B, PAX9, JARID2, TRPS1, RORA, TFAP2C

Table 3. KEGG pathway terms of the selected 93 DEGs.

KEGG terms	P	Fold enrichment	Gene symbols
Platelet activation	0.01	6.56	GNAI3, ROCK2, RASGRP1, COL3A1, PLCB1
Glutamatergic synapse	0.03	5.97	PLD1, GNAI3, PLCB1, SLC1A1
Pathways in cancer	0.03	2.90	GNAI3, FGFR3, ROCK2, RASGRP1, TGFA, PLCB1, DAPK1
Sphingolipid signaling pathway	0.04	5.33	PLD1, GNAI3, ROCK2, PLCB1
Axon guidance	0.04	5.16	GNAI3, SEMA4G, ROCK2, EPHA2

For further elucidation of interactions between the aberrant miRNAs considered herein and their potential common target DEGs, one miRNA-target gene network was constructed using Cytoscape (Figure 3(B)). The network was composed of two basic elements (miRNAs and genes) and their connections, which represented miRNA-target gene interactions. Our results showed that one miRNA could target multiple genes, and one gene could be affected by multiple miRNAs. The miRNAs with the highest numbers of interactions were hsa-miR-20a, let-7c, hsa-miR-24, hsa-miR-106a, hsa-miR-206, hsa-miR-200c, hsa-miR-221, hsa-miR-222, hsa-miR-18* and hsa-miR-203, each one possessing no less than 5 connections. Moreover, 6 target DEGs (ABHD2, FNDC3B, LASP1, ARHGAP12, TRPS1 and PHACTR4) were found to have higher-degree distributions in the network. We assumed that genes with higher-degree distributions in the regulatory network played an important role in BE pathogenesis. ABHD2, FNDC3B and LASP1 were predicted to be regulated by 4 miRNAs, while ARHGAP12, TRPS1 and PHACTR4 were each predicted to be the targets of 3 miRNAs. Next, we performed ROC analysis on the 6 genes based on GSE39491 dataset because of its sample number. The ROC curves indicated that the expression of these 6 genes provided a means of discriminating between BE and NE tissues, as shown in Figure 4. Among these, FNDC3B expression provided better differentiation conditions with a sensitivity of 92.86, a specificity of 93.33 and an AUC value of 0.930 (P values < 0.05) (Figure 4(B)).

Figure 4. Receiver operating characteristic curve analysis for discriminating BE from NE (BE vs. NE); (A) ROC curve of ABHD2; (B) ROC curve of FNDC3B; (C) ROC curve of LASP1; (D) ROC curve of ARHGAP12; (E) ROC curve of TRPS1; (F) ROC curve of PHACTR4.

Discussion

Gastroesophageal reflux disease (GERD) [27] is a common chronic disorder that occurs predominantly when the contents of the stomach are refluxed into the esophagus. The clinical symptoms not only include gastrointestinal symptoms, such as acid reflux and heartburn, but also extraesophageal symptoms. In 2006, the first GERD global consensus was determined using evidence-based medicine, referred to as the Montreal definition [27]. Experts in the Asia Pacific Region made several revisions, after which the Asia Pacific Regional Consensus of GERD was established and renewed [1]. Using upper gastrointestinal endoscopy, GERD patients could be divided into cases of non-erosive reflux disease (NERD) and reflux esophagitis (RE), and were found to involve Barrett's esophagus and esophageal adenocarcinoma (EAC) among other complications [27]. In addition, patients who presented no GERD symptoms may also have been diagnosed with GERD upon exhibiting GERD-related complications or pathological reflux, as determined by examinations. As one of the complications of GERD, BE is characterized by pathological metaplasia in the distal esophagus, in which the squamous epithelial lining is replaced by columnar epithelium with or without intestinal goblet cells [1]. Because intestinal metaplasia with goblet cells is more relevant to EAC, an aggressive malignant cancer, we focused on Barrett's esophagus with intestinal metaplasia [1,28]. MiRNAs are a group of non-coding small RNAs that play an important role in regulating a spectrum of basic cellular processes [7,8]. Endogenous miRNAs interact with the 3'-UTR of target mRNA to suppress post-transcription of genes [5,6]. Previous studies have underscored the relevance of miRNAs as potential biomarkers for disease diagnosis, cancer prediction and prognosis assessment, including BE and EAC [10,11,13,15,17–20,22,29–33]. These studies mostly focused on BE carcinogenesis and were less concerned with its pathogenesis. Therefore, the pathological mechanism of intestinal metaplasia associated with BE is not fully understood. However, studies regarding BE pathogenesis could assist with preventing its occurrence and progression, which would constitute a significant development. To date, only a few studies regarding BE pathogenesis have been reported, so the resulting lack of data entices more researches in this field [3]. Microarrays provide a new methodology for exploring changes in BE pathogenesis on the molecular level. A series of studies have indicated that BE pathogenesis may be related to various factors, including aberrant miRNAs, and that dysregulated miRNAs could be relevant to BE development [15,18–21]. However, studies involving the interactions between aberrant miRNAs and target DEGs in BE are still rare. To our knowledge, the present study is the first to reveal interactions between dysregulated miRNAs and differentially expressed genes in miRNAs and mRNA expression profiling datasets. Identification and analyses of aberrant miRNAs and target DEGs in BE may reveal the underlying pathogenesis of BE on the molecular level and provide improvements in BE diagnosis, treatment and prevention.

After sophisticated arithmetic based on R [34], the microarray results will be visualized and presented as tables corresponding to differentially expressed genes [35]. In the present study, GSE24839, an miRNAs expression profiling dataset containing 10 BE and 10 NE tissue samples, was retrieved and analyzed. Then, 23 dysregulated miRNAs were identified, and bi-cluster analysis was conducted to evaluate data reliability. We found that most of the BE samples could be well differentiated from NE tissue samples, so these dysregulated miRNAs were credible and could be used for further analysis. Among these 23 miRNAs, 7 were significantly up-regulated in BE and the other 16 miRNAs were down-regulated in BE. The 7 up-regulated miRNAs were hsa-miR-215, hsa-miR-192, hsa-miR-194, hsa-miR-10a, hsa-miR-145*, hsa-miR-200a and hsa-miR-217, while the top 6 dysregulated miRNAs were hsa-miR-18a*, hsa-miR-205, hsa-miR-215, hsa-miR-203, hsa-miR-192 and hsa-miR-194, with high values of Log₂|FC|. A series of other studies [10,11,15,18,20,22] has already revealed that miRNAs expression in BE is quite different from that in NE. Luzna P et al. [13] proposed that miR-192 and miR-196a were up-regulated with an increase during BE progression. Cabibi D et al. [17] found increased expression levels of miR-143, miR-145, miR-194 and miR-215 in BE, while miR-203 and miR-205 were down-regulated. Wijnhoven B P et al [19]. indicated that miR-21, miR-143, miR-145, miR-193 and miR-215 were expressed in higher levels, and miR-203 and miR-205 were down-regulated in BE compared to NE. Moreover, miR-205, which was down-regulated in BE, was expressed at levels equivalent to NE after endoscopic ablation [33]. Leidner R S et al. [10] also investigated miRNAs expression by next-generation sequencing and found that, among 26 miRNAs related to EAC, as many as 23 miRNAs demonstrated early-stage (BE) dysregulation. Furthermore, Slaby O et al. [11] found that miR-21, miR-25, miR-194 and miR-196a expression levels increased gradually during BE progression to EAC, while miRNA-203, miRNA-205, miRNA-210 and miRNA-378 decreased during the same period. Previous studies have demonstrated the down-regulation of miR-203 in BE as well in other squamous epithelia such as the skin and upper respiratory tract [19,36,37]. This action may target p63 and Zfp281 to regulate stem cell renewal in stratified epithelia [38]. In contrast, expression levels of miR-215 were high in BE. Braun CJ et al. [39] demonstrated that miR-215 may be up-regulated by genotoxic stress, inducing cell cycle arrest, and that p53 played an important role in induction. Therefore, it is possible that miR-215 may act as a tumor suppressor during BE development. Although miRNAs expression profiling in BE has received much attention of late, the regulating mechanisms involved remain unclear, such that the functions of most dysregulated miRNAs remain to be clarified.

In addition to the miRNAs microarray, we applied R/Bioconductor statistics to the mRNA expression data processing [34]. Finally, 632 common DEGs were identified among the five mRNA expression profiling datasets. Interestingly, all 632 DEGs in these datasets displayed consistent expression tendencies. That is, the five datasets selected are of high quality and our results are reliable. To investigate the interactions between the aberrant miRNAs and DEGs derived from microarrays, we first predicted target genes for these 23 dysregulated miRNAs utilizing three main bioinformatic tools (TargetScan, RNA22-HSA and miRDB) [40–42]. In this way, 1962 target genes were obtained. To improve precision, we obtained the intersection of these 1962 target genes and the 632 DEGs and identified 93 potential DEGs.

Subsequently, we utilized DAVID Bioinformatics Resources to analyze these 93 potential DEGs. Functional annotation was performed to reveal the enriched GO and KEGG terms. We found that 5 GO_BP, 4 GO_CC and 1 GO_MF terms were involved. Among these GO biological processes, digestive tract development and epithelial cell differentiation were found to be significantly related to BE pathogenesis. The genes involved in digestive tract development process were CLDN18, TRPS1 and COL3A1 with a fold enrichment value of 18.24. The genes involved in epithelial cell differentiation process were TAGLN, CNN3 and ANXA4 with a fold enrichment value of 10.09. BE is defined as a metaplasia condition in the distal esophagus, in which the squamous epithelial lining is replaced by columnar epithelium [1]. Therefore, the digestive tract development and the epithelial cell differentiation processes are particularly important. Moreover, pathways in cancer were also involved in the 5 enriched KEGG pathway terms with a fold enrichment value of 2.90. The predominant related pathways included the MAPK and VEGF signaling pathways, which have been shown to regulate sustained angiogenesis and cell proliferation [23–26]. BE is the only well-recognized precancerous lesion of EAC, so it was no surprise to find that pathways in cancer were involved. The MAPK (mitogen-activated protein kinase, MAPK) signaling pathway involves an important group of proteins that carry extracellular signals from cell membrane receptors to the nucleus of the cell. This pathway is comprised of the following three principal components: ERK (extracellular signal-regulated kinase, ERK), JNK (c-Jun NH2-terminal kinase, JNK) and p38 kinase [43]. These proteins regulate cell proliferation, differentiation, and many other important biological progresses depending on signaling mediators [23,24,43]. Souza RF et al. [44] have shown that acid can activate MAPK signaling in BAR-T and EAC cells, mediating cell proliferation and apoptosis suppression. Angiogenesis is quite important for the formation of new vascular supplies and is mainly regulated by VEGF (vascular endothelial growth factor, VEGF) signaling pathway [25]. VEGF not only affects angiogenesis in endothelial cells but also plays an important role in tumor cell proliferation and apoptosis inhibition [26]. Many different tumor cells have been reported to co-express VEGF and its receptors, including Barrett's metaplastic epithelium [26,45]. However, the underlying mechanisms remain to be revealed, in particular, the regulatory role of miRNAs.

For further elucidation, the miRNA-target gene network was constructed using Cytoscape with miRNAs and genes as the two basic elements. The miRNAs demonstrating the highest numbers of interactions were hsa-miR-20a, let-7c, hsa-miR-24, hsa-miR-106a, hsa-miR-206, hsa-miR-200c, hsa-miR-221, hsa-miR-222, hsa-miR-18* and hsa-miR-203, all of which possessed no fewer than 5 connections. Moreover, ABHD2, FNDC3B and LASP1 were predicted to be targets of 4 miRNAs, while ARHGAP12, PHACTR4 and TRPS1 were predicted to be targets of 3 miRNAs. It has been shown that one miRNA can target multiple genes and one gene can be affected by multiple miRNAs, consistent with previous researches [31,46]. The aforementioned 6 genes demonstrating higher-degree network distributions were selected for ROC analysis. The ROC curves showed that FNDC3B expression provided superior differentiation, exhibiting a sensitivity of 92.86, a specificity of 93.33 and an AUC value of 0.930 (Figure 4(B)). FNDC3B, which is located at 3q26 and may amplified in several cancers [47], was identified as an important oncogenic driver of 3q amplicons, such as hepatocytes, and mammary and kidney epithelial cells. Recently, Xu H et al. [48] speculated that FNDC3B may be related to cellular process inhibition, including viability, proliferation, migration and invasiveness. However, FNDC3B studies are somewhat rare, and we could find no published data regarding FNDC3B and BE. This facet indicates the importance of the present study with regard to understanding BE pathogenesis; indeed, this work may provide insight into new biomarkers for BE diagnosis, prediction and therapeutic evaluation. However, more studies are required to elucidate the exact mechanisms of dysregulated miRNAs and miRNA-targeted DEGs.

Although microarray analysis has some limitations, this methodology can provide additional information without selection bias [49]. The sample size of the microarray data utilized herein was not large enough for most of the datasets, so we selected five mRNA expression profiling datasets to reduce the bias caused by insufficient sample size. In the same way, the influence of tissue heterogeneity could also be minimized. However, the quantity of the miRNA profiling datasets was quite small online, so further investigation with larger samples is needed urgently. A study that correlates miRNAs profiling with the clinical outcomes would be a very valuable addition to this field.

Altogether, our results provide evidence that biological processes like digestive tract development and epithelial cell differentiation may be involved in BE pathogenesis and that MAPK and VEGF pathways may also participate in the processes. Furthermore, our results indicate the potential of FNDC3B as a new biomarker in BE diagnosis and prediction.

Materials and methods

Retrieval of microarray datasets from an online database

Microarray datasets of miRNAs and gene expression profiling from BE and normal esophageal mucosa (NE) were obtained and downloaded from Gene Expression Omnibus (GEO) Datasets on the National Center for Biotechnology Information (NCBI) Database (http://www.ncbi.nlm.nih.gov/gds/). The search keywords were Barrett's esophagus. The selection was restricted to Homo sapiens and expression profiling. Finally, the non-coding miRNA profiling dataset (GSE24839) and the mRNA expression profiling datasets (GSE13083, GSE13898, GSE26886, GSE34619 and GSE39491) were downloaded for further analyses.

Bioinformatic analyses of aberrant miRNAs and DEGs

R/Bioconductor statistical analysis language and software [34,50] were utilized for miRNAs and mRNA expression microarray processing, and limma algorithm package was used for identification of differentially expressed miRNAs and mRNA genes. NE was used as control in all these datasets, such that comparisons between BE and NE were conducted. Differentially expressed miRNAs with cut-off-values of at least 2-FC and P values < 0.05 were considered as aberrant miRNAs. DEGs in BE compared to those in NE were defined as commonly expressed and exhibiting consistent tendencies in all five mRNA expression profiling datasets with cut-off-values of at least 2-FC and P values < 0.05.

The target genes of aberrant miRNAs were predicted using three main bioinformatic tools, namely TargetScan (http://www.targetscan.org/), RNA22-HSA (https://cm.jefferson.edu/) and miRDB (http://mirdb.org/cgi-bin/search.cgi). The target genes of each aberrant miRNA were defined as common in all the lists generated from the above three bioinformatic tools. Subsequently, the Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources (https://david.ncifcrf.gov/) was used for functional annotation [35], and the miRNA-target genes interaction network was constructed using Cytoscape [51].

Statistical analysis

All statistical analyses were performed using SPSS software (version 18.0) and GraphPad Prism software (version 6). Furthermore, receiver operating characteristic (ROC) curve analysis was performed using Med-Calc statistical software (version 15.2.2.2), and P values < 0.05 were considered statistically significant when sensitivity, specificity and the area under the curve (AUC) were calculated.

Disclosure of potential conflicts of interest

The authors report no conflict of interest.

Acknowledgements

This work was supported by grants from the National Natural Science Foundation of China (No. 81702067), China Postdoctoral Science Foundation (#2017M613176) and Shaanxi Natural Science Foundation of China (#2017JQ8041).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (NSFC) [grant number 81702067]; Shaanxi Natural Science Foundation of China [grant number 2017JQ8041]; China Postdoctoral Science Foundation [grant number 2017M613176].