This article is mentioned in:

Introduction

Esophageal cancer (EC), also known as oesophageal cancer, is one of the most aggressive tumors of the gastrointestinal tract, and ranks the sixth most common cause of cancer mortality (1). According to the statistical data of National Central Cancer Registry of China, an estimated 287,632 patients were diagnosed with and 208,473 patients succumbed to EC in 2010 in China (2). Although the EC-associated mortality rate has decreased over the past 30 years, it remains a major cancer burden in high risk areas (3). Thus, the early diagnosis and treatment of EC are vital.

Previous studies have indicated that tumor-specific microRNAs (miRNAs/miRs) are considered to be potential biomarkers for the early diagnosis and treatment of EC (4–6). miRNAs, a class of endogenous short non-protein-coding RNA molecules, bind to their target sites within the mRNA molecules via base-pairing with complementary sequences at a post-transcriptional level, and thus certain miRNAs are considered to be oncogenes or tumor suppressor genes (7–9). Previous studies have suggested that miRNAs are primarily expressed in EC tissue and plasma. For example, miR-145, miR-133a and miR-133b were demonstrated to function as tumor suppressors by targeting FSCN1 in esophageal squamous cell carcinoma (ESCC) (10). miR-375 may inhibit tumor growth and metastasis in ESCC by repressing insulin-like growth factor 1 receptor (11). miR-92a has been revealed to promote tumor metastasis of ESCC via epithelial cadherin (12). In addition, miR-21, miR-184 and miR-221 in the plasma were demonstrated to be correlated with the recurrence of ESCC, and serve as oncogenes in ESCC (13). However, miRNA expression in the saliva of patients with EC remains largely unknown.

Due to the extensive blood supply to the salivary glands, saliva is a promising bodily fluid for use in the early detection of diseases including cancer, infectious and cardiovascular diseases (14,15). Previous studies have also suggested that saliva molecules may be used to detect human organic and systemic diseases (16,17). In addition, Xie et al (18) analyzed the expression levels of miRNAs in saliva samples and revealed that miR-10b, miR-144 and miR-451 in whole saliva were significantly upregulated in patients with EC. These miRNAs demonstrate promise as biomarkers for the detection of EC due to the convenience and noninvasive nature of saliva collection (18). However, the target genes of theses miRNAs were not mentioned and the precise functions of these miRNAs remain to be fully understood. Using the same data as Xie et al (18), the present study aimed to identify the validated target genes of the miRNAs that were differentially expressed between whole saliva samples from patients with EC and healthy controls, and to analyze the potential functions of the targets by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.

Materials and methods

Microarray data and data preprocessing

The miRNA expression profile (GSE41268) was extracted from the Gene Expression Omnibus (GEO) database (19) (http://www.ncbi.nlm.nih.gov/geo/) based on the Agilent-021827 Human miRNA Microarray platform (V3; miRBase release 12.0 miRNA ID version; Agilent Technologies, Inc., Santa Clara, CA, USA). This dataset was deposited by Xie et al (18). The microarray contained 10 chips of whole saliva samples from 7 patients with EC and 3 healthy controls. The pathology of the patients with EC was squamous cell carcinoma: 1 patient was stage I, 1 patient was stage II, 3 patients were stage III and 2 patients were stage IV. The probe-level data in CEL files were converted into the miRNA expression profiles. For miRNAs corresponding to multiple probes that exhibited a plurality of expression values, the miRNA expression values of those probes were averaged. Log2 transformation was performed on probes that mapped with the miRNA names labeled in the annotation platform to normalize the data from a skewed distribution to the approximate normal distribution (20).

Screening of differentially expressed miRNAs

The Linear Models for Microarray Data package (21) was used to normalize the data and identify the differentially expressed miRNAs between patients with EC and healthy controls. miRNAs with the cutoff criteria of log2fold change (FC) >1 and P<0.05 were considered to be significantly differentially expressed. In order to explore whether the miRNAs were sample-specific, the Pheatmap package (22) in R was used to perform hierarchical clustering by comparing the value of each miRNA in 10 samples.

Target genes for differentially expressed miRNAs

miRecords (http://c1.accurascience.com/miRecords/download.php) is an integrated resource for animal miRNA-target interactions (23). The miRecords database currently contains 1,135 records of interactions between 644 miRNAs and 1,901 target genes in 9 animal species. In the present study, the validated target genes were extracted according to the miRecords database.

Interaction network analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database is a useful tool that provides multiple experimental and predicted data (24,25). In order to study the associations between miRNAs and their targets, the target genes were scanned by the STRING database, and protein-protein interaction (PPI) pairs with the cutoff criterion of a combined score >0.4 were selected. Then, the regulatory network that contained these miRNAs and the PPI pairs were constructed using Cytoscape software (version 3.2.1; http://cytoscape.org/) (26).

Function enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (27) was used to identify the enriched functions in DEGs. KEGG Orthology-Based Annotation System (KOBAS) is a software to identify significantly enriched pathways using hypergeometric tests (28). To explore the functions and pathways of target genes in the development of EC, the DAVID and KOBAS were used to identify the significantly enriched GO and KEGG categories, respectively. P<0.05 was used as a cutoff criterion.

Results

Screening of differentially expressed miRNAs

With the cutoff criteria of log2FC >1 and P<0.05, 18 differentially expressed miRNAs were identified in EC samples when compared with the healthy controls (Table I). The top 3 upregulated miRNAs were miR-144, miR-451 and miR-509-3p, while the top 3 downregulated miRNAs were miR-125b-1, miR-129 and miR-636. In addition, to investigate the differences among the 10 samples, a heat map was generated to compare their expression values. The hierarchical clustering analysis revealed a clearly distinct expression of all differentially expressed miRNAs between patients with EC and healthy controls (Fig. 1).

Figure 1. Hierarchical clustering heat map of 18 differentially expressed miRNAs. Each column corresponds to a single microarray, whereas each row indicates the expression profile of a single gene. Red and blue represent high and low miRNA expression, respectively. N-1, 2 and 3 indicate salivary samples in healthy controls, while EC-1, 2, 3, 4, 5, 6 and 7 indicate salivary samples in patients with esophageal cancer. miRNA, microRNA.

Table I. Differently expressed miRNAs in esophageal cancer tissue compared with normal tissue.

Table I.

Differently expressed miRNAs in esophageal cancer tissue compared with normal tissue.

miRNA ID	P-value	log2FC	miRNA ID	P-value	log2FC
hsa-miR-125b-1	0.0245	−4.5437	hsa-miR-144	0.0060	6.8525
hsa-miR-129	0.0024	−3.7921	hsa-miR-451	0.0107	6.8225
hsa-miR-636	0.0422	−2.1583	hsa-miR-509-3p	0.0408	5.0251
hsa-miR-1909	0.0487	−2.1353	hsa-miR-486-5p	0.0105	4.6848
hsa-miR-1274b	0.0008	−1.8183	hsa-miR-10b	0.0497	4.5171
hsa-miR-1274a	0.0119	−1.5339	hsa-miR-363	0.0396	4.2118
hsa-miR-720	0.0453	−1.4932	hsa-miR-139-3p	0.0206	3.7638
hsa-miR-371-5p	0.0206	−1.3268	hsa-miR-648	0.0454	3.7265
hsa-miR-933	0.0048	−1.3112	hsa-miR-98	0.0434	3.3743

[i] FC, fold change; miRNA/miR, microRNAs.

Screening of validated target mRNAs

According to the miRecords database, 43 validated target mRNAs corresponding to 7 upregulated miRNAs (including miR-144, miR-451 and miR-98) were acquired (Table II). No validated target gene for downregulated miRNAs was obtained.

Table II. Target genes of differentially expressed miRNAs.

Table II.

Target genes of differentially expressed miRNAs.

miRNA	Target	Test methods
hsa-miR-144	NOTCH1	Luciferase reporter assay
hsa-miR-144	FGA	Luciferase reporter assay
hsa-miR-144	FGB	Luciferase reporter assay
hsa-miR-144	FGG	ELISA; flow; luciferase reporter assay
hsa-miR-144	PLAG1	Luciferase reporter assay
hsa-miR-451	AKT1	RT-qPCR; western blotting
hsa-miR-451	MMP9	RT-qPCR; western blotting
hsa-miR-451	MMP2	RT-qPCR; western blotting
hsa-miR-451	ABCB1	RT-qPCR; luciferase reporter assay; western blotting
hsa-miR-451	BCL2	RT-qPCR; western blotting
hsa-miR-451	CAB39	RT-qPCR; luciferase reporter assay; western blotting; microarray
hsa-miR-451	MIF	ELISA; luciferase reporter assay; microarray; RT-qPCR; western blotting
hsa-miR-451	SSSCA1	RT-qPCR; western blotting
hsa-miR-509-3p	NTRK3	Luciferase reporter assay
hsa-miR-98	E2F1	Luciferase reporter assay
hsa-miR-98	E2F2	Northern blot
hsa-miR-98	MYC	Luciferase reporter assay
hsa-miR-98	ACADM	RT-qPCR
hsa-miR-98	AMMECR1	RT-qPCR
hsa-miR-98	CBX5	RT-qPCR
hsa-miR-98	CCL5	ELISA
hsa-miR-98	CHRNB1	RT-qPCR
hsa-miR-98	CNOT4	RT-qPCR
hsa-miR-98	EZH2	Northern blot
hsa-miR-98	HMGA2	Northern blot; RT-qPCR
hsa-miR-98	HNRPDL	RT-qPCR
hsa-miR-98	KLF13	RT-qPCR
hsa-miR-98	MEIS1	RT-qPCR
hsa-miR-98	NCOA3	Luciferase reporter assay
hsa-miR-98	SERP1	RT-qPCR
hsa-miR-98	SLC20A1	RT-qPCR
hsa-miR-98	SOCS4	Luciferase reporter assay; northern blot; RT-qPCR; western blotting
hsa-miR-98	THBS1	RT-qPCR
hsa-miR-98	TUSC2	Luciferase reporter assay; western blotting
hsa-miR-98	ZFP36L1	RT-qPCR
hsa-miR-98	ZMPSTE24	RT-qPCR
hsa-miR-10b	PPARA	Luciferase reporter assay; RT-qPCR; western blotting
hsa-miR-10b	KLF4	Luciferase reporter assay; western blotting
hsa-miR-10b	NCOR2	RT-qPCR; western blotting
hsa-miR-10b	HOXD10	Luciferase reporter assay
hsa-miR-10b	NF1	GFP reporter assay; luciferase reporter assay; microarray; RT-qPCR
hsa-miR-486-5p	CD40	Microarray; RT-qPCR
hsa-miR-363	CDKN1A	RT-qPCR; luciferase reporter assay; western blotting

[i] miRNA/miR, microRNA; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; GFP, green fluorescent protein.

Interaction network analysis

A total of 6 miRNAs and 27 validated target genes were mapped into the network using STRING and Cytoscape. In the network, 5 genes including Notch homolog 1 (NOTCH1), fibrinogen γ chain (FGG), fibrinogen α chain (FGA), and fibrinogen β chain (FGB) were modulated by miR-144; while 6 genes including AKT serine/threonine kinase 1 (AKT1), matrix metalloproteinase (MMP)9, and MMP2 were modulated by miR-451. A total of 11 genes, including v-myc avian myelocytomatosis viral oncogene homolog (c-MYC; MYC), E2F transcription factor (E2F)1, and E2F2 were modulated by miR-98. A total of 3 genes, including peroxisome proliferator-activated receptor α (PPARA) and Kruppel-like factor 4 (KLF4) were modulated by miR-10b. TNF receptor superfamily member 5 (CD40) and cyclin-dependent kinase inhibitor 1A (CDKN1A) were modulated by miR-485-5p and miR-363, respectively (Fig. 2).

Figure 2. miRNA-target regulatory network. The triangle nodes represent miRNAs, while oval nodes represent target genes. The different miRNAs and their target genes are marked with different colors. miRNA/miR, microRNA.

Function and pathway enrichment analysis

To gain additional insight into the underlying functions of all target genes and the target genes in the network, GO and KEGG enrichment analyses were performed. The results demonstrated that the top 3 enriched GO terms were cell cycle (E2F1, E2F2 and MYC), blood vessel development (AKT1, NOTCH1 and MMP2) and vasculature (AKT1, NOTCH1 and MMP2; Table III); and the top 3 enriched KEGG pathways were cell cycle (E2F1, E2F2 and MYC), complement and coagulation cascades (FGG, FGA and FGB) and epidermal growth factor receptor signaling pathway (AKT1, CDKN1A and MYC; Table IV).

Table III. Enriched GO terms of target genes in the microRNA-target regulatory network.

Table III.

Enriched GO terms of target genes in the microRNA-target regulatory network.

Term	Count	P-value	Genes
GO: 0007049~cell cycle	10	4.27×10−4	E2F1, E2F2, MYC, AKT1, SSSCA1, TUSC2, CDKN1A, BCL2, HMGA2, THBS1
GO: 0001568~blood vessel development	6	8.23×10−4	AKT1, NOTCH1, MMP2, ZFP36L1, NF1, THBS1
GO: 0001944~vasculature development	6	9.18×10−4	AKT1, NOTCH1, MMP2, ZFP36L1, NF1, THBS1
GO: 0022402~cell cycle process	8	1.39×10−3	E2F1, MYC, AKT1, SSSCA1, CDKN1A, BCL2, HMGA2, THBS1,
GO: 0043067~regulation of programmed cell death	9	2.65×10−3	MYC, AKT1, NOTCH1, MMP9, CDKN1A, BCL2, NF1, THBS1, MIF
GO: 0010941~regulation of cell death	9	2.71×10−3	MYC, AKT1, NOTCH1, MMP9, CDKN1A, BCL2, NF1, THBS1, MIF
GO: 0051726~regulation of cell cycle	6	3.09×10−3	E2F1, E2F2, MYC, AKT1, CDKN1A, BCL2
GO: 0022403~cell cycle phase	6	7.93×10−3	E2F1, AKT1, SSSCA1, CDKN1A, BCL2, HMGA2
GO: 0008283~cell proliferation	6	9.80×10−3	E2F1, MYC, TUSC2, BCL2, CD40, MIF

[i] GO, gene ontology.

Table IV. Enriched Kyoto Encyclopedia of Genes and Genomes pathways of target genes in the microRNA-target regulatory network.

Table IV.

Enriched Kyoto Encyclopedia of Genes and Genomes pathways of target genes in the microRNA-target regulatory network.

Term	Pathway	Count	P-value	Genes
hsa04110	Cell cycle	3	5.49×10−9	E2F1, E2F2, MYC
hsa04610	Complement and coagulation cascades	3	1.43×10−4	FGG, FGA, FGB
hsa04012	ErbB signaling pathway	3	1.52×10−3	AKT1, CDKN1A, MYC
hsa04620	Toll-like receptor signaling pathway	3	2.15×10−3	AKT1, CD40, CCL5

[i] ErbB, epidermal growth factor receptor.

Discussion

miRNA expression is known to be associated with various types of human cancer (29). In the present study, the salivary miRNA expression profile of EC deposited by Xie et al (18) in GEO database was analyzed, to screen for biomarkers for therapeutic intervention. In concordance with the study of Xie et al (18), the present study also demonstrated that miR-144, miR-451 and miR-10b were significantly upregulated in saliva samples from patients with EC and healthy controls. In addition, miR-144 was revealed to potentially be involved in EC progression by directly targeting 5 mRNAs, including NOTCH1, FGG, FGA and FGB. NOTCH 1 is a member of the NOTCH family. NOTCH signaling serves an important function in a variety of developmental processes by regulating cell differentiation, proliferation and apoptosis (30). NOTCH signaling has been demonstrated to be upregulated in EC, and may be used as a therapeutic target for EC (31,32). FGG, FGA and FGB belong to the fibrinogen family. Fibrinogen is a glycoprotein in vertebrates that assists with the formation of blood clots and functions as a common component of hemostasis (33). Takeuchi et al (34) identified that plasma fibrinogen concentration was significantly associated with chemoradiotherapy responsiveness in patients with ESCC. Matsuda et al (35) suggested that the plasma fibrinogen level may be used as a biomarker for postoperative recurrence in advanced EC. The fibrinogen level is correlated with tumor metastasis and recurrence in patients with ESCC (34,35). Thus, by targeting NOTCH1, FGG, FGA and FGB, miR-144 may be involved in cell invasion and metastasis of EC.

miR-451 was revealed to target 6 mRNAs in the present study, including AKT1, MMP9 and MMP2. The protein kinase B (AKT) family has been demonstrated to integrate extracellular signals with cell processes including proliferation, differentiation, migration and survival (36). Previous data have demonstrated that the phosphoinositide 3-kinase/protein kinase B/mechanistic target of rapamycin pathway is frequently activated in numerous types of human cancer, including colorectal cancer (37). Cancer cells are well known to exhibit a higher proliferation rate compared with normal cells, and cancer cells frequently fail to undergo apoptosis (38). Activated AKT may regulate its downstream targets to increase cell proliferation and decrease apoptosis (39). AKT activation was described as an early cellular response to carcinogen exposure, and may be a key step in environmental carcinogenesis (40). In addition, MMP2 and MMP9 were demonstrated to be overexpressed in ESCC and to be associated with esophageal tumorigenesis (41). Therefore, the increased expression of MMPs may promote tumor proliferation and inhibit apoptosis in the development of EC (42).

The present study revealed that miR-10b may regulate 3 mRNAs (including PPARA and KLF4). PPARA, also termed NR1C1, is a transcription factor that is a member of the peroxisome proliferator-activated receptor (PPAR) group. PPARs have been demonstrated to stimulate tumor cell proliferation and induce neoangiogenesis in different types of gastrointestinal cancer, including EC, liver and pancreatic cancer (43). KLF4, as a known tumor suppressor, may regulate cell proliferation, differentiation, development and apoptosis (44). KLF4 has been suggested to suppress EC cell migration and invasion, as a direct target of miR-10b (45). Therefore, the associated miRNAs, including miR-144 and miR-451 and miR-10b, may serve important functions in cellular processes including cell differentiation, proliferation, metastasis and apoptosis in EC by regulating NOTCH, FGG, FGA, FGB, AKT1, MMP9, MMP2, PPARA and KLF4.

Compared with the aforementioned previous study, among the upregulated miRNAs, miR-98 and miR-363 were novel (18). miR-98 belongs to the lethal-7 (let-7) family of miRNAs (let-7f, let-7d, miR-98 and let-7g). The miRNA let-7 may regulate key differentiation processes during development (46). The expression of let-7 has been suggested to predict the response to chemotherapy in EC (47). miR-98 may target almost half of the validated target genes in the network, including MYC, E2F1, and E2F2. MYC, a well-known transcription factor, exhibited a high degree of connectivity in the regulatory network of the present study. E2F1 and E2F2 are transcription factors and members of the E2F family. MYC and E2Fs have been identified as associated with the cell cycle (48). MYC overexpression enhances and MYC downregulation inhibits cell cycle progression (48). E2F1 expression is positively associated with cell proliferation in EC (49). E2F2 has been demonstrated to negatively regulate a subset of genes involved in the processes of DNA metabolism and cell cycle control (50). In addition, E2F1 and E2F2 are expressed in a cell cycle-regulated manner, exhibiting their highest levels in the late G1 and S phase (51). Cell cycle deregulation is a common feature of human tumors (52). The results of the present study also demonstrated that MYC, E2F1 and E2F2 were enriched in the cell cycle, which was the most significantly enriched GO function and KEGG pathway, suggesting that MYC, E2F1 and E2F2 serve critical functions in the cell cycle (Fig. 3). Therefore, by targeting MYC, E2F1 and E2F2, miR-98 may be involved in the cell cycle of EC cells, and act as novel biomarker for EC.

Figure 3. Detailed molecular network of pathways of the cell cycle (Kyoto Encyclopedia of Genes and Genomes ID, has04110). Red genes are the significant genes involved in the pathway. G1, G1 phase; S, S phase; G2, G2 phase; M, mitotic phase.

miR-363 was identified to be downregulated in head and neck squamous cell carcinoma (HNSCC) tissues, and to affect HNSCC cell invasion and metastasis (53). Hsu et al (54) indicated that miR-363 served a key function in the increment of gastric carcinogenesis via targeting c-Myc promotor binding protein 1, a negative regulator of MYC. In the present study, miR-363 was suggested to directly target CDKN1A, also termed p21, WAF1 or CIP1. The cyclin-dependent kinase inhibitor CDKN1A is a regulator of cell cycle progression at the G1 and S phases (55). Therefore, miR-363 may be involved in EC by regulating its target gene, and may be used as a novel biomarker for EC.

In addition, in contrast to the study of Xie et al (18), the present study identified that miR-486-5p was significantly upregulated in EC samples (P=0.01052). This may be caused by the different preprocessing methods used in these studies. The present study also identified that miR-486-5p may modulate EC by targeting CD40. CD40, a member of the TNF-receptor superfamily, was originally identified on the surface of B cells and was suggested to regulate B cell growth and differentiation (56). The interaction of CD40 with the CD40 ligand was revealed to modulate immune cells in the lymph nodes and to promote growth of local squamous cell cancer of the head and neck by protecting tumor cells (57,58). miR-486-5p was dysregulated in EC saliva samples and was demonstrated to regulate CD40 in the present study. However, additional studies are needed to verify the involvement of miR-486-5p in EC.

In conclusion, the present study demonstrated that the biomarkers (miR-144, and miR-451 and miR-10b) screened by Xie et al (18) may serve an important function in the regulation of EC by targeting NOTCH, fibrinogen, AKT1, MMPs, PPARA and KLF4. In addition, 2 novel biomarkers, miR-98 and miR-363, were identified. These are associated with the EC and modulate a serial of cell cycle-associated genes, including MYC, E2F1, E2F2 and CDKN1A. Additional experiments are needed to confirm the expression of miR-486-5p and its involvement in EC, due to the discordant conclusions between the present study and a previous study (18).

Acknowledgements

The present study was supported by Natural Science Foundation of Liaoning Province (grant no. 2015 020561).

References