Identification of key pathways and genes in Barrett's esophagus using integrated bioinformatics methods.

Barrett's esophagus (BE) is a premalignant lesion of esophageal adenocarcinoma. The aim of the present study was to investigate the possible mechanisms and biomarkers of BE. To identify the differentially expressed microRNAs (DEmiRNAs) and genes (DEGs) in BE, the miRNA expression profile GSE20099 and the gene expression profiles GSE26886, GSE13083 and GSE34619 were obtained from the Gene Expression Omnibus (GEO) database. DEGs and DEmiRNAs were screened for using the GEO2R tool. Using DAVID, functional and pathway enrichment analysis was performed to explore the biological function of identified DEGs. The protein‑protein interaction (PPI) network was detected using STRING and constructed by Cytoscape software. Furthermore, targets of identified DEmiRNAs were predicted by the miRecords database, then integrated with the identified DEGs to obtain key genes involved in BE. In total, 311 DEGs were identified. These genes were significantly enriched in the pancreatic secretion, metabolic pathways and cytochrome P450 drug metabolism pathways. In the PPI network, 16 hub genes, including keratin 16, cystic fibrosis transmembrane conductance regulator, involucrin, protein kinase C α and cadherin 17 were identified. Following integration of the predicted target genes of DEmiRNAs with DEGs, three key BE genes were identified: PRKCA, CDH17 and epiregulin. In conclusion, a comprehensive bioinformatics analysis of identified DEGs and DEmiRNAs was performed to elucidate potential pathways and biomarkers involved in the development of BE.

Introduction

Barrett's esophagus (BE) is defined as a metaplasia of the lower esophagus, where the normal squamous epithelium is replaced by intestinalized columnar epithelium in response to gastroesophageal reflux disease (GERD). Population based studies estimate that the prevalence of BE is increasing, with the current incidence at 1.3–1.6% in general population of Sweden and Italy (1). Notably, BE is a premalignant lesion of esophageal adenocarcinoma (EAC) and significantly enhances the risk of EAC development (2). However, the molecular mechanism of BE development is still largely unknown.

A number of studies focused on the development of BE have been carried out in the previous decade. Bile acid is a pathogenic factor of BE (3) and has been reported to promote its development through the upregulation of caudal type homeobox 2 expression (4). A recent study also demonstrated that the hedgehog pathway may contribute to the development of BE via regulation of the forkhead box A2 gene (5). Furthermore, expression levels of the intercellular adhesion molecule 1 and Kruppel like factor 4 genes were demonstrated to be significantly higher in BE tissues compared with control samples (6,7), whereas Rho-kinase gene expression was almost unchanged (8). Although these studies have made progress, most have only focused on disparate genes, which is not sufficient for providing an overall mechanism of BE development. Thus, the aim of the current study was to explore possible molecular mechanisms and potential biomarkers of BE using bioinformatic methods. Datasets were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) and identified differentially expressed genes (DEGs) and differentially expressed microRNAs (DEmiRNAs) using the GEO2R program. Subsequently, DEGs were analyzed using functional and pathway enrichment analysis, and then their regulatory networks were constructed. Following prediction of the target genes of the DEmiRNAs, the target genes were integrated with the relevant DEGs to further identify the key genes involved in BE.

Materials and methods

Data resources

mRNA expression profiles (GSE26886, GSE13083 and GSE34619) and the miRNA expression profile of GSE20099 were downloaded from the GEO data repository. The dataset of GSE26886 contains 20 BE samples and 19 normal esophageal samples. GSE13083 includes 7 BE samples and 7 normal samples. GSE34619 consists of 10 BE samples and 8 normal samples. The miRNA expression dataset of GSE20099 contains 14 BE samples and 14 normal samples.

Differential expression analysis

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is an online tool that performs comparisons on GEO datasets based on the GEOquery and Limma R packages (9). The BE group and normal group were selected, and the GEO2R program was subsequently applied for differential expression analysis. The genes and miRNAs that met the cut-off criteria of the adjusted P-value (adj. P) <0.01 and |log fold change| >2 were considered as DEGs and DEmiRNAs.

Functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) is a gene functional enrichment program, providing a large series of functional annotation tools for researchers to decipher the biological implications behind huge amounts of genes. To understand the deeper biological meaning of DEGs, Gene Ontology (GO) and KEGG pathway enrichment analysis of identified DEGs was performed using DAVID (10). GO functional analysis consists of three categories: Biological process (BP), cellular component (CC) and molecular function (MF). P<0.05 was set as the threshold value.

Protein-protein interaction (PPI) network

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://string.embl.de/) is an online database containing known and predicted PPI networks. In this research, a PPI network of identified DEGs in all three datasets was identified using the STRING database (combined score >0.4) and subsequently visualized using Cytoscape (http://www.cytoscape.org/) software (version 3.4.0; The Cytoscape Consortium, San Diego, CA, USA) (11,12). The Molecular Complex Detection (MCODE) program within Cytoscape was used to detect modules of the PPI network (13). The cut-off criteria were set as follows: Degree, ≥2; node score, ≥0.2; k-score, 0.2; maximum depth, 100. The function and pathway enrichment analysis of the identified modules was then performed using the DAVID database.

Prediction of miRNA targets

miRecords (http://c1.accurascience.com/miRecords/) is an online program which integrates miRNA target predictions made by 11 established prediction databases (NBmiRTar, RNA22, MicroInspector, PicTar, PITA, RNAhybrid, and TargetScan, DIANA-microT, miRanda, MirTarget2 and miTarget) (14). In this study, potential targets of DEmiRNAs were identified by at least four databases.

Results

Gene expression profiles GSE26886, GSE13083 and GSE34619 identified 2,112, 941 and 604 DEGs respectively, and 311 DEGs were identified across all three datasets (Fig. 1). Genes present in upregulated and downregulated groups in different datasets were excluded, and 311 DEGs still presented stable trends in at least two datasets, including 163 upregulated genes and 148 downregulated genes in BE samples compared with normal esophageal samples.

Figure 1.

Identification of differentially expressed genes in gene expression profile datasets GSE26886, GSE13083 and GSE43619. Percentages displayed below the number of genes in each section indicates the proportion of genes in that section in regards to all three data sets.

Functional and KEGG pathway enrichment analysis

The upregulated genes enriched in BP terms were mainly associated with digestion, tissue homeostasis and microvillus organization, upregulated genes enriched in CC terms were mainly associated with extracellular vesicles and upregulated genes enriched in MF terms were mainly associated with actin binding. Furthermore, KEGG pathway enrichment analysis demonstrated that the upregulated genes were mainly enriched in pancreatic secretion, metabolic pathways and salivary secretion pathways (Table I). The downregulated genes enriched in BP terms were mainly associated with epidermis development, keratinocyte differentiation and skin development, in CC terms the genes were mainly associated with extracellular regions and in MF terms the genes were mainly associated with receptor antagonist activity. In addition, KEGG pathway enrichment analysis revealed that the downregulated genes were mainly associated with cytochrome P450 drug metabolism and the amoebiasis pathway (Table II).

Table I.

Functional and KEGG pathway enrichment analysis of upregulated genes in BE.

A, Biological process
Term	Description	Count	P-value
GO:0007586	Digestion	13	3.91×10−8
GO:0001894	Tissue homeostasis	13	3.50×10−7
GO:0032528	Microvillus organization	6	9.49×10−6
GO:0030277	Maintenance of gastrointestinal epithelium	5	1.09×10−5
GO:0016266	O-glycan processing	7	1.68×10−5
B, Cellular component
Term	Description	Count	P-value
GO:1903561	Extracellular vesicle	73	1.24×10−15
GO:0043230	Extracellular organelle	73	1.26×10−15
GO:0070062	Extracellular exosome	72	3.61×10−15
GO:0031988	Membrane-bounded vesicle	78	1.55×10−12
GO:0044421	Extracellular region part	80	7.07×10−12
C, Molecular function
Term	Description	Count	P-value
GO:0003779	Actin binding	14	6.49×10−5
GO:0022853	Active ion transmembrane transporter activity	7	4.72×10−4
GO:0050839	Cell adhesion molecule binding	13	0.001010
GO:0008092	Cytoskeletal protein binding	18	0.001705
GO:0022804	Active transmembrane transporter activity	11	0.001831
D, KEGG
Term	Description	Count	P-value
hsa04972	Pancreatic secretion	8	4.20×10−5
hsa01100	Metabolic pathways	27	1.16×10−4
hsa04970	Salivary secretion	6	0.001815
hsa05110	Vibrio cholerae infection	5	0.002050
hsa04976	Bile secretion	5	0.005352

Count, the number of enriched genes in each term; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table II.

Functional and KEGG pathway enrichment analysis of downregulated genes in BE.

A, Biological process
Term	Description	Count	P-value
GO:0008544	Epidermis development	21	3.50×10−13
GO:0030216	Keratinocyte differentiation	14	1.24×10−11
GO:0043588	Skin development	17	7.84×10−11
GO:0009913	Epidermal cell differentiation	15	1.71×10−10
GO:0031424	Keratinization	8	1.40×10−7
B, Cellular component
Term	Description	Count	P-value
GO:0005576	Extracellular region	73	5.14×10−9
GO:0044421	Extracellular region part	63	7.59×10−8
GO:0070062	Extracellular exosome	50	3.22×10−7
GO:1903561	Extracellular vesicle	50	3.76×10−7
GO:0043230	Extracellular organelle	50	3.81×10−7
C, Molecular function
Term	Description	Count	P-value
GO:0048019	Receptor antagonist activity	3	0.002808
GO:0030547	Receptor inhibitor activity	3	0.004791
GO:0005149	Interleukin-1 receptor binding	3	0.006381
GO:0005198	Structural molecule activity	14	0.010380
GO:0005509	Calcium ion binding	13	0.013760
D, KEGG
Term	Description	Count	P-value
hsa00982	Drug metabolism-cytochrome P450	4	0.013555
hsa05146	Amoebiasis	4	0.042967

Count, the number of enriched genes in each term; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network and modules analysis

A total of 204 nodes and 382 edges were obtained from the PPI network program. A degree of >10 was set as the cut-off criterion for hub gene identification. A total of 16 hub genes were identified: Keratin 16 (KRT16), cystic fibrosis transmembrane conductance regulator (CFTR), involucrin (IVL), protein kinase C α (PRKCA), mucin 2, oligomeric mucus/gel-forming (MUC2), amyloid beta precursor protein (APP), cadherin 17 (CDH17), mucin 6, oligomeric mucus/gel-forming (MUC6), MET proto-oncogene, receptor tyrosine kinase (MET), envoplakin (EVPL), desmoglein 1 (DSG1), keratin 8 (KRT8), peroxisome proliferator activated receptor g (PPARG), prominin 1 (PROM1), cytochrome P450 family 2 subfamily E member 1 (CYP2E1) and filaggrin (FLG; Table III and Fig. 2). The three most significant modules were then selected from the PPI network using MCODE (Fig. 3). Notably, not all of the hub genes were present in these modules. Functional enrichment analysis of these modules was then performed. Genes in module A were mainly enriched in skin development, keratinocyte differentiation and keratinization. Genes in module B were predominantly enriched in O-glycan processing, mucin type O-glycan biosynthesis and metabolic pathways. Finally, genes in module C were mainly enriched in homeostatic processes, pancreatic secretion and the bile secretion pathway (Table IV).

Table III.

Identified hub genes in the protein-protein interaction network.

	GSE26686		GSE13083		GSE34619
Gene	logFC	adj. P	logFC	adj. P	logFC	adj. P
KRT16	−4.62	3.33×10−4	−5.07	5.63×10−3	−3.52	2.61×10−5
CFTR	2.35	8.6×10−5	5.33	1.08×10−6	4.31	1.61×10−11
IVL	−6.29	2.10×10−6	−4.35	2.53×10−4	−4.27	1.79×10−11
PRKCA	2.79	6.77×10−7	3.31	2.07×10−5	2.36	2.63×10−11
MUC2	4.21	4.34×10−5	6.72	3.16×10−6	2.20	2.77×10−3
APP	2.41	4.04×10−6	3.56	1.08×10−5	2.09	2.31×10−8
CDH17	5.58	8.00×10−5	7.36	5.66×10−8	4.31	2.77×10−4
MUC6	6.11	1.02×10−6	3.10	7.45×10−3	5.94	1.46×10−11
MET	2.60	7.40×10−9	4.21	1.27×10−6	2.32	9.04×10−9
EVPL	−3.66	2.23×10−7	−3.55	6.81×10−4	−2.10	7.23×10−12
DSG1	−3.85	2.13×10−3	−4.84	4.51×10−5	−5.48	1.06×10−7
KRT8	5.73	5.13×10−15	6.42	2.18×10−9	4.01	3.50×10−11
PPARG	4.56	2.02×10−12	3.37	1.48×10−6	3.02	1.66×10−11
PROM1	8.47	1.02×10−17	7.37	4.24×10−11	5.44	2.03×10−14
CYP2E1	−4.34	1.66×10−7	−4.88	1.43×10−6	−2.06	1.59×10−8
FLG	−4.40	5.87×10−4	−6.16	3.09×10−5	−4.96	1.14×10−10

FC, fold change; adj. P, adjusted P-value.

Figure 2.

The protein-protein network of identified differentially expressed genes. Red nodes and green nodes indicate upregulated genes and downregulated genes, respectively; V-shape nodes represent the identified hub genes.

Figure 3.

Top three modules identified from the protein-protein interaction network. (A) Module A, (B) module B and (C) module C. Red nodes and green nodes indicate upregulated genes and downregulated genes, respectively.

Table IV.

Functional and KEGG pathway enrichment analysis of gene modules.

A, Module A
Term	Description	Count	P-value
GO:0043588	Skin development	7	2.13×10−10
GO:0030216	Keratinocyte differentiation	6	9.43×10−10
GO:0031424	Keratinization	5	4.99×10−9
GO:0009913	Epidermal cell differentiation	6	7.20×10−9
GO:0008544	Epidermis development	6	9.84×10−8
B, Module B
Term	Description	Count	P-value
GO:0016266	O-glycan processing	7	1.41×10−13
GO:0006493	Protein O-linked glycosylation	7	4.24×10−12
GO:0005975	Carbohydrate metabolic process	9	1.91×10−10
GO:0043413	Macromolecule glycosylation	7	1.65×10−9
GO:0006486	Protein glycosylation	7	1.65×10−9
hsa00512	Mucin type O-Glycan biosynthesis	3	1.16×10−4
hsa01100	Metabolic pathways	4	0.019425
C, Module C
Term	Description	Count	P-value
GO:0042592	Homeostatic process	10	5.77×10−6
GO:0019725	Cellular homeostasis	8	6.83×10−6
GO:0030003	Cellular cation homeostasis	7	1.09×10−5
GO:0006873	Cellular ion homeostasis	7	1.24×10−5
GO:0055080	Cation homeostasis	7	2.09×10−5
hsa04972	Pancreatic secretion	6	1.07×10−6
hsa04976	Bile secretion	4	3.98×10−4
hsa04971	Gastric acid secretion	4	4.70×10−4
hsa04970	Salivary secretion	4	7.60×10−4
hsa04070	Phosphatidylinositol signaling system	4	0.001112

Count, the number of enriched genes in a certain term; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Integrated analysis of DEmiRNAs and paired DEGs

The miRNA expression profiles identified five DEmiRNAs, including two upregulated miRNAs [miRNA-215 (miR-215) and miR-192] and three downregulated miRNAs (miR-205, miR-203 and miR-486-5p) in BE samples compared with normal esophageal samples. Among them, miR-215 and miR-205 were the most significantly upregulated and downregulated miRNAs, respectively. In addition, predicted targets of DEmiRNAs were obtained from miRecords database. Considering the fact that an inverse association exists between the expression of miRNA and its target mRNA, DEmiRNAs with target genes that were identified as DEGs were selected. Remarkably, 33 pairs of DEmiRNAs and DEGs with an inverse relationship of expression met this criterion. Among these target genes, PRKCA was the target of two downregulated miRNAs (miR-203 and miR-205), while epiregulin (EREG) was the target of two upregulated miRNAs (miR-215 and miR-192). Furthermore, two hub genes (PRKCA and CDH17) were predicted as targets of miR-203 (Table V).

Table V.

Differentially expressed miRs and their paired DEGs in Barret's esophagus.

miR	Adjusted P-value	Log fold change	Paired DEGs
hsa-miR-215	1.05×10−11	5.972098	EREG
hsa-miR-192	1.89×10−6	2.665345	EREG
hsa-miR-205	1.76×10−3	−3.301436	PRKCA, ABHD2,
			KDELR3, DPP4,
			HSD17B11, MLPH
hsa-miR-203	3.04×10−3	−2.784608	PRKCA, CDH17,
			CPS1, TOX3,
			CPD, STS, PIGR,
			FUT8, KRT20,
			PBLD
hsa-miR-486-5p	2.77×10−4	−2.506280	TM4SF20

miR, microRNA; DEGs, differentially expressed genes.

Discussion

BE is the precursor lesion of EAC and is still not completely understood (2). Therefore, the aim of the current study was to provide an overall view of the molecular mechanism and biomarkers in BE. In the present research, a total of 311 DEGs and five DEmiRNAs were identified from GEO databases using the GEO2R program. The functional analysis of 163 upregulated DEGs demonstrated these genes were mainly associated with digestion, tissue homeostasis and microvillus organization, while the functional analysis of 148 downregulated DEGs demonstrated that they were associated with epidermis development, keratinocyte differentiation and skin development. Following construction of the PPI network, 16 hub genes were identified. The target genes of DEmiRNAs were subsequently integrated with the DEGs, leading to the identification of three key genes that may be regulated by miRNAs.

KEGG pathway enrichment analysis revealed that the upregulated DEGs were mainly associated with pancreatic secretion, metabolic pathways, salivary secretion, vibrio cholera infection and bile secretion. Many clinical studies have focused on these pathways. Greer et al (15) reported that the insulin/insulin-like growth factor (IGF) pathway, which is associated with pancreatic secretion, may have a role in the development of BE. Regarding the importance of metabolic pathways, accumulating evidence suggests that obesity is a central risk factor for BE (16–18). Furthermore, high serum leptin and adiponectin have been shown to be associated with the development of BE (19–21).

As mentioned previously, BE is a columnar metaplasia in response to GERD. GERD can be a result of the impairment of the anti-reflux barrier, which in turn is associated with salivary secretion function (22,23). In addition, bile exposure has confirmed effects in the development of BE (24,25). Reveiller et al (26) demonstrated that the expression of squamous differentiation genes such as IVL and DSG1 was inhibited by bile exposure in human esophageal epithelial cells. Notably, IVL and DSG1 are among our identified hub genes.

The downregulated DEGs identified were mainly associated with the cytochrome P450 drug metabolism pathway. Supporting this finding, a study previously demonstrated that cytochrome P450 family 1 subfamily A member 2 (CYP1A2), CYP3A4 and CYP2E1 were highly expressed in regions of active cell proliferation in BE (27).

The results of the functional and pathway analysis performed on the gene modules are largely similar to the DEGs mentioned above. The module enrichment analysis also indicated that skin development, mucin type O-glycan biosynthesis and homeostatic processes were associated with the development of BE. The normal esophageal epithelium is made up of squamous epithelium, which also constitutes the skin. The genes associated skin development were all downregulated. Clinical manifestation of this change in gene expression is observed in the replacement of squamous epithelium with columnar epithelium that is typical of BE (1). Notably, genes associated with mucin type O-glycan in module B were all upregulated in BE samples, which was correlated with the degree of intestinal metaplasia in BE. Mucin type O-glycan has been reported to have important roles in intestinal homeostasis (28). A previous study indicated that dysfunction of homeostatic processes may be vital in the development of BE (29). Therefore, a focus on these pathways may provide us with greater knowledge of BE development.

miRNAs are a group of endogenous non-coding RNA molecules that act as negative regulators in post-transcriptional gene regulation. Recent studies reported that aberrant miRNA expression is associated with the development of BE (30,31). The current study found five differentially expressed miRNAs, including two upregulated DEmiRNAs and three downregulated DEmiRNAs. Among them, miR-215 and miR-205 were the most significantly upregulated and downregulated DEmiRNAs, respectively. Other studies have also confirmed the aberrant expression of these two miRNAs in BE and EAC (32,33). miR-215 is involved in many diseases, including diabetic nephropathy, ovarian cancer, non-small cell lung cancer and multiple carcinomas. miR-215 commonly acts as a tumor suppressor by promoting apoptosis (34–36). miR-205 also has an effect in various diseases and may act as a prognostic marker for many malignant neoplasms (37,38). This integrated study identified that two hub genes (PRKCA and CDH17) were target genes of miR-203, and EREG was a target of two upregulated miRNAs (miR-215 and miR-192).

As a member of the serine/threonine-specific protein kinase family, PRKCA is involved in various biological processes, including the innate immune response, angiogenesis, insulin secretion, cell proliferation, adhesion and migration. In addition, PRKCA was also identified to be involved in several pathways, including the mitogen-activated protein kinase (MAPK) pathway, phosphoinositide 3-kinase (PI3K)-Akt pathway and Wnt pathway. Previous studies have demonstrated that the MAPK and PI3K pathways are involved in the proliferation of BE-associated EAC (39), while the insulin/IGF and Wnt pathways have important roles in BE development (15,40).

CDH17, which is located on chromosome 8, encodes liver-intestinal cadherin (LI-cadherin), which is involved in cell adhesion and oligopeptide transport. LI-cadherin is localized at cell junctions and is selectively expressed in enterocytes and intestinal goblet cells, but not in the esophagus (41). Multiple studies have reported that aberrant LI-cadherin expression is a sensitive marker for detection of various metaplastic diseases such as gastric intestinal metaplasia and BE (42,43).

EREG, located on chromosome 4, encodes the epiregulin protein, which is involved in various biological processes, including the innate immune response, cytokine production, wound healing and epithelial cell proliferation. A previous study also demonstrated that epiregulin regulates the differentiation of airway epithelial cells through the epidermal growth factor receptor pathway (44). Furthermore, epiregulin has roles in a series of other malignant diseases as well as esophageal cancer, and can be used as biomarker in certain diseases (45–47). Taken together, this data indicates that the genes PRCKA, CDH17 and EREG may be involved in the development of BE through various pathways.

There are certain limitations to the present study. Due to the fact that there is limited research available on the pathways involved in BE, the identified KEGG pathways could not be discussed in detail. Additionally, excluding miRNA data, the datasets together provided only 37 BE samples and 34 control samples. Samples obtained from different platforms may also suffer from some bias. The high level of redundancy in the GO database was also a limiting factor in the study. Further research should focus on methods to evaluate possible bias that may arise in this type of research.

In conclusion, the present study performed a comprehensive bioinformatics analysis of identified DEGs to identify potential pathways and biomarkers involved in the development of BE. In the current research, 311 DEGs and 5 DEmiRNAs were identified. Following integration of the DEGs with any corresponding DEmiRNAs, three key genes (PRKCA, CDH17 and EREG) were identified that may be associated with BE. This research provides a novel insight into molecular mechanisms that may underlie the development of BE. However, biochemical analysis and further research are still necessary to validate these results.