Identification of a PLCE1‑regulated competing endogenous RNA regulatory network for esophageal squamous cell carcinoma.

Phospholipase C epsilon 1 (PLCE1) and the competing endogenous RNA (ceRNA) network are crucial for tumorigenesis and the progression of esophageal squamous cell carcinoma (ESCC). However, whether PLCE1 can regulate the ceRNA network in ESCC has not been clarified. In the present study, we aimed to identify the PLCE1‑regulated ceRNA network and further elucidate the regulatory mechanisms by which ESCC is promoted. Microarray analysis was used to identify differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) from three pairs of samples of PLCE‑silenced Eca109 and control Eca109 cells. Next, the ceRNA regulatory network was established and visualized in Cytoscape, and functional enrichment analysis was performed to analyze DEGs from ceRNAs. Protein‑protein interaction (PPI) networks among the DEGs were established by the STRING database to screen hub genes. Kaplan‑Meier survival analysis was used to validate hub genes. Finally, PLCE1‑related hub gene/lncRNA/miRNA axes were also constructed based on the ceRNA network. A total of 105 DELs and 346 DEGs were found to be dysregulated in the microarray data (|log2FC| >1.5, adjusted P<0.05). We constructed a PLCE1‑regulated ceRNA network that incorporated 12 lncRNAs, 43 miRNAs, and 169 mRNAs. Functional enrichment analysis indicated that the DEGs might be associated with ESCC onset and development. A PPI network was established, and 9 hub genes [WD and tetratricopeptide repeats 1 (WDTC1), heat shock protein family A (Hsp70) member 5 (HSPA5), N‑ethylmaleimide sensitive factor, vesicle fusing ATPase (NSF), fibroblast growth factor 2 (FGF2), cyclin dependent kinase inhibitor 1A (CDKN1A or P21), bone morphogenetic protein 2 (BMP2), complement C3 (C3), GM2 ganglioside activator (GM2A) and discs large MAGUK scaffold protein 4 (DLG4)] were determined from the network. Kaplan‑Meier survival analysis validated four hub genes (BMP2, CDKN1A, GM2A, and DLG4) that were treated as prognostic factors. Ultimately, hub gene/lncRNA/miRNA subnetworks were obtained based on the 4 hub genes, 13 DEmiRNAs, and 10 DELs. In conclusion, the PLCE1‑regulated ceRNA contributes to the onset and progression of ESCC and the underlying molecular mechanisms may provide insights into personalized prognosis and new therapies for ESCC patients.

Introduction

Esophageal squamous cell carcinoma (ESCC), the most common subtype of esophageal cancer (ESCA), is the sixth most common cause of cancer-related death in the world (1). The treatment of ESCC has improved over the past few decades, with techniques that primarily include surgical resection, chemotherapy, and radiotherapy (2,3). Despite advances in treatment, the 5-year overall survival (OS) rate for ESCC patients remains unsatisfactory (1). Recent investigations have found that several prognostic factors, such as tumor stage, smoking, and alcohol, can be used to predict the OS of ESCC patients, yet the prognosis is still poor (4–7). Since we still do not fully understand how ESCC occurs and progresses, it is urgent to explore the underlying molecular mechanisms and identify novel prognostic biomarkers and potential therapeutic targets for ESCC.

Phospholipase C epsilon 1 (PLCE1), a new ESCC susceptibility locus at chromosome 10q236, has been identified by genome-wide association studies (GWAS) in Chinese Han populations (8). It has been reported that PLCE1 mediates different external signals that are relevant to tumor stage and survival in hepatocellular carcinoma, gastric carcinoma, bladder cancer and colorectal cancer (9–12). Our previous investigations also indicated that PLCE1 at the level of mRNA and protein is highly increased in ESCC from Han and Kazakh populations and that the overexpression of PLCE1 is correlated with metastasis and tumor aggressiveness (13–16). More importantly, PLCE1 is also associated with poor prognosis of patients with ESCC (13). Zhai and coworkers also confirmed that when PLCE1 is knocked out, the transcriptional activity of SNAIL is particularly repressed (17). Finally, the invasion and migration of ESCC are inhibited when PLCE1 is knocked out, suggesting that PLCE1 could be regarded as a marker for the molecular typing of ESCC (16), and we have also confirmed that PLCE1 can function as an oncogene by enhancing ESCC cell proliferation and inhibiting apoptosis. Nevertheless, while PLCE1 is overexpressed in ESCC, the specific mechanism by which it promotes tumorigenesis has not been clarified.

Competitive endogenous RNA (ceRNA) is a newly proposed regulatory mechanism by which lncRNAs can adsorb microRNAs (miRNAs), to participate in the regulation of target genes (18), thus playing an essential role in many cancers, such as ESCC (19), lung cancer (20) and gallbladder cancer (21). Therefore, it is clear that lncRNA-miRNA-mRNA networks participate in the initiation and development of not only ESCC but also other types of cancers. Nevertheless, no report has indicated that PLCE1-regulated ceRNA is involved in ESCC progression according to microarray or sequencing data.

In this study, microarray analysis was used to identify differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) from three pairs of samples of PLCE-silenced Eca109 and control Eca109 cells. The DELs and DEGs from microarray data were used to establish a ceRNA regulatory network. Subsequently, we predicted the cellular functions and pathways of DEGs involved in the ceRNA network by GO and KEGG Pathway analysis. Protein-protein interaction (PPI) analysis of DEGs was performed to screen for crucial hub genes. Finally, PLCE1-related hub gene/lncRNA/miRNA axes were obtained from the ceRNA, which provide insights into the regulatory mechanisms of ESCC.

Materials and methods

Cell culture and PLCE1-silenced microarray

The Eca109 cells purchased from Shanghai Institute of Biochemistry and Cell Biology of the Chinese Academy of Sciences were cultured at 37°C in DMEM or RPMI-1640 (Gibco; Thermo Fisher Scientific, Inc.), containing 10% FBS (Sigma-Aldrich; Merck KGaA) in moist air containing 5% CO2. Stable PLCE1 siRNA (siRNA-PLCE1, 5′-GGGUCUUGCCAGUCGACUATT-3′) transfected cells or control cells were previously constructed (16) by the transfection of PLCE1 siRNA or si-RNA control, respectively. Briefly, Eca109 cells were grown in 6-well plates and transfected with si-PLCE1 and siRNA control using HiPerFect Transfection Reagent (Qiagen) according to the manufacturer's instructions. At 72 h after transfection, cells were harvested for western blot analyses. In addition, three pairs of PLCE-silenced Eca109 cells and controls were selected for a full transcription chip, Affymetrix Clariom™ D Array (Affymetrix, Thermo Fisher Scientific, Inc.) which was performed by Compass Biotechnology Co., Beijing, China. Array signal intensities were analyzed with Expression Console Software (version 1.4.1, Affymetrix; Thermo Fisher Scientific, Inc.). Furthermore, we used another si-PLCE1 chip in EC9706 cells, the Affymetrix GeneChip Human genome 3 (Affymetrix, Thermo Fisher Scientific, Inc.), to validate results of the si-PLCE1 chip in Eca109 cells.

Western blot analysis

Total protein was extracted using radio immunoprecipitation assay (RIPA) buffer. Protein (50 µg) was loaded onto SDS-PAGE gels (8–10%) and electrophoresed (200 V, 2 h), transferred onto PVDF membranes (BD Biosciences) and hybridized with primary antibodies, against β-actin (diluted 1:1,000, Zhongshan), PLCE1 (diluted 1:200, Santa Cruz Biotechnology, Inc.) followed by incubation with the appropriate secondary antibodies diluted to 1:10,000 (Zhongshan).

Cell proliferation and apoptosis assay

Ethynyl-2′-deoxyuridine (EdU) assays were conducted to examine cell proliferation using a Beyotime Assay Kit (C0071S, Shihezi). Cells (1×105/well) were cultured using si-PLCE1 for 72 h, fixed by 4% formaldehyde for 25 min at 4°C and 1X Click Reaction Buffer (Beyotime Institute of Biotechnology) was used to process them for approximately 30 min in the dark, followed by DAPI staining. After being rinsed with PBS, images were obtained under a fluorescence microscope (×200 magnification).

Mitochondrial membrane potential was measured using JC-1. JC-1 can potential-dependently accumulate in mitochondria as indicated by a color change from green to red. Therefore, the decline in the red/green fluorescence intensity ratio indicates mitochondrial depolarization. After PLCE1 silencing treatment, the medium was replaced with 500 µl fresh medium, including 5 µg/ml JC-1 (C2005, Beyotime Institute of Biotechnology). One hour after labeling in the dark, the images were obtained under a fluorescence microscope (×200 magnification).

TUNEL assays applying the Beyotime Apoptosis Kit (C1086, Beyotime Institute of Biotechnology) were carried out on PLCE1-silenced cells. Cells (1×105/well) were cultured with si-PLCE1 for 72 h, fixed by 4% formaldehyde for 25 min at room temperature, and permeabilized with 0.3% Triton-X-100 for 5 min. Before being stained with DAPI for 5 min at room temperature and washed with PBS, cells were mixed with TUNEL reagent for 60 min at 37°C. Images were obtained under a fluorescence microscope (×200 magnification).

Identification of DELs and DEGs

We used Perl and R scripts (https://www.r-project.org/; http://www.perl.org/) for all data processing and normalized analysis and usage. Furthermore, the lncRNA and mRNA expression profiles were analyzed, respectively, to screen DELs and DEGs in the microarray. Dysregulated genes were identified by using the edgeR package of the R language (22). A |log2foldchange (log2FC)|>1.5 and adjusted P<0.05 were set as the thresholds to define differences.

ceRNA network

To better assess the relationships between PLCE1 and ceRNAs in ESCC, we constructed a ceRNA network based on the theory that miRNAs can lead to gene silencing by interacting with mRNAs, whereas lncRNAs can function as miRNA sponges to regulate gene expression by competitively binding to miRNAs. We used miRcode tools (23) (http://www.mircode.org/) to predict which DELs and DEmiRNAs could interact with each other. DEGs targeted by DEmiRNAs were retrieved based on several databases, including Targetscan (24), miRTarBase (25), and miRDB (26). Furthermore, we built the ceRNA network relationship through the above DELs, DEmiRNAs, and DEGs, and the established network was visualized by Cytoscape software (27).

Functional enrichment analysis

To identify the function of DEGs, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using Metascape (28). P-value <0.05 was chosen as the cutoff criterion. Significant GO terms with similar functions were visualized as interaction networks using the Enrichment Map plugin in Cytoscape (27).

Establishment of the PPI network and hub gene subnetworks

The DEGs were extracted by the STRING database (29) (https://string-db.org/) to assess the interactive relationships of DEGs from ceRNA, the network was visualized by Cytoscape, and hub genes were obtained by screening the degree of connectivity of each node in the network. Lastly, we established hub gene/lncRNA/miRNA networks to explore the molecular mechanisms.

Statistics and analysis procedure

Statistical analyses were performed with R software (https://www.r-project.org/; version 3.6.1). Results are expressed as mean ± standard deviation. P<0.05 was regarded as a statistically significant difference. Differences between means were investigated by Student's t-test within two groups. The survival value of hub genes were analyzed using Kaplan-Meier (KM) survival analysis, followed by a log-rank test. Fig. 1 demonstrates the analysis procedure of the data mining processes.

Figure 1.

Flow chart of the approach utilized in the present study. PLCE1, phospholipase C epsilon 1; DElncRNAs (DELs), differentially expressed long non-coding RNAs; DEmRNAs (DEGs), differentially expressed messenger RNAs; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ceRNA, competitive endogenous RNA; PPI, protein-protein interaction.

Results

Validation of PLCE1 knockdown

To confirm the oncogene function of PLCE1, si-PLCE1 was made to interfere with its expression in Eca109 cells. The expression of PLCE1 was determined by western blot analysis after transfection. The results revealed that PLCE1 protein expression levels were highly reduced in the PLCE1-silenced (si-PLCE1) cells compared with the controls (si-con) (Fig. S1), suggesting that the expression of PLCE1 was successfully silenced after siRNA transfection of Eca109 cells.

In this study, we used a microarray to study the expression of DELs and DEGs in three paired PLCE1-silenced ECa109 cells and controls (Fig. 2A-D). In total, 105 lncRNAs and 346 mRNAs were identified as dysregulated between two groups from the microarray (|log2FC|>1.5, adjusted P-value <0.05).

Figure 2.

lncRNA and mRNA expression profiles in PLCE1-silenced ESCC cells and control cells. (A) Hierarchical clustering shows a distinguishable lncRNA expression profile and (B) mRNA expression profile. (C) Volcano plot of the differentially expressed lncRNAs and (D) mRNAs in PLCE1-silenced ESCC cells and ESCC controls. The red and green points in the plot represent differentially expressed lncRNAs or mRNAs with statistical significance. PLCE1, phospholipase C epsilon 1; ESCC, esophageal squamous cell carcinoma; lncRNA, long non-coding RNA; mRNA, messenger RNA.

Functional enrichment analysis of PLCE1-related DEGs

Functional enrichment analysis was used to explore the role of DEGs (adjusted P-value <0.05). The GO analysis indicated that upregulated genes were related with negative regulation of cell growth (GO: 0030308), negative regulation of protein homodimerization activity (GO: 0090074), and cellular lipid catabolic processes (GO: 0044242) (Fig. 3A and B). In contrast, downregulated genes were related with positive regulation of protein serine/threonine kinase activity (GO: 0071902), positive regulation of neurogenesis (GO: 0050769) and I-kappaB kinase/NF-kappaB signaling (GO: 0007249) (Fig. 4A and B). Furthermore, the pathway analysis revealed that upregulated genes were involved in the metabolism of lipids, asparagine N-linked glycosylation, and constitutive signaling by AKT1 E17K in cancer (Fig. 3A and B); and downregulated genes were enriched in MAPK signaling pathways and the activation of anterior HOX genes in hindbrain development during early embryogenesis (Fig. 4A and B).

Figure 3.

GO and Pathway enrichment analysis for upregulated genes using Metascape analysis. (A) Network of enriched sets colored by ID. Threshold: 0.3 kappa score; similarity score >0.3. (B) Heatmap colored by P-values (P<0.05 cutoff was regarded as statistical importance). GO, Gene Ontology.

Figure 4.

GO and Pathway enrichment analysis for downregulated genes using Metascape analysis. (A) Network of enriched sets colored by ID. Threshold: 0.3 kappa score; similarity score >0.3. (B) Heatmap colored by P-values (P<0.05 cutoff was regarded as statistical importance). GO, Gene Ontology.

PLCE1 enhances proliferation and inhibits apoptosis in ESCC in vitro

To investigate the effect of PLCE1 on ESCC growth and apoptosis in vitro, several phenotype experiments were incorporated in PLCE1-silenced Eca109 cells. As shown in Fig. 5A-C, TUNEL and JC-1 staining experiments showed that the apoptosis rate of ESCC cells was markedly increased, whereas proliferation was significantly suppressed according to the EdU assays when PLCE1 was knocked down in the Eca109 ESCC cells. This indicated that PLCE1 in Eca109 cells can induce proliferation and attenuate apoptosis.

Figure 5.

(A) TUNEL cell apoptosis assay, (B) EdU cell proliferation and (C) JC-1 cell apoptosis assay were used to evaluate the apoptosis and proliferation ability of ESCC Eca109 cells after treatment with si-PLCE1. PLCE1, phospholipase C epsilon 1; ESCC, esophageal squamous cell carcinoma.

lncRNA/miRNA/mRNA ceRNA networks and functional assessment

To fully comprehend the role of DELs and DEGs by differentially expressed miRNAs (DEmiRNAs), a ceRNA network was constructed. A total of 566 interaction pairs were determined. After intersecting with the DEmiRNAs, 123 DEL/DEmiRNA pairs remained, containing 12 DELs and 43 DEmiRNAs. First, we searched for 105 DELs in the miRCode database and obtained 123 interacting pairs of lncRNAs and miRNAs using the Perl language, which contained 12 DELs and 43 DEmiRNAs. We further searched for mRNAs targeted by the 43 DEmiRNAs in TargetScan, miRDB, and miRanda databases, and selected those that overlapped with the identified DEGs. Finally, a total of 443 DEmiRNA/DEG interactions were added to the ceRNA network (Fig. 6). We also investigated the potential function of these DEGs in the ceRNA network using Metascape analysis (28) and found that they mainly participate in fatty acid metabolic processes, DNA damage response and negative regulation of signal transduction by p53 (Fig. S2; Table I). The pathway analysis indicated that DEGs were involved in metabolism of lipids and regulation of TP53 expression and degradation (Fig. S2; Table I).

Figure 6.

The PLCE1-regulated lncRNA-miRNA-mRNA network in ESCC. The triangle stands for the selected lncRNAs, the circle represents the mRNAs, and the rhombus represents the miRNAs. PLCE1, phospholipase C epsilon 1; ceRNA, competing endogenous RNA; lncRNA, long non-coding RNA; miRNA, microRNA; mRNA, messenger RNA; ESCC, esophageal squamous cell carcinoma.

Table I.

Functional enrichment analysis of DEGs involved in the ceRNA network.

GO	Description	Count	Log10(P)	Log10(q)
GO:0045664	Regulation of neuron differentiation	15	−4.65	−0.84
GO:0010575	Positive regulation of vascular endothelial growth factor production	4	−4.47	−0.84
GO:0006631	Fatty acid metabolic process	11	−4.39	−0.84
GO:0019228	Neuronal action potential	4	−4.19	−0.84
R-HSA-6798695	Neutrophil degranulation	12	−4.18	−0.84
GO:0046466	Membrane lipid catabolic process	4	−4.09	−0.84
R-HSA-556833	Metabolism of lipids	15	−4.07	−0.84
GO:0021675	Nerve development	5	−3.8	−0.73
GO:0006978	DNA damage response, signal transduction by p53 class mediator resulting in transcription of p21	3	−3.78	−0.73
GO:0051961	Negative regulation of nervous system development	9	−3.68	−0.73
GO:0007169	Transmembrane receptor protein tyrosine kinase signaling pathway	14	−3.62	−0.72
GO:0001763	Morphogenesis of a branching structure	7	−3.55	−0.7
GO:0098876	Vesicle-mediated transport to the plasma membrane	5	−3.37	−0.57
GO:0045807	Positive regulation of endocytosis	6	−3.31	−0.57
GO:1901797	Negative regulation of signal transduction by p53 class mediator	3	−2.91	−0.42
GO:0030308	Negative regulation of cell growth	6	−2.89	−0.42
GO:0051235	Maintenance of location	8	−2.87	−0.42
GO:0019216	Regulation of lipid metabolic process	9	−2.86	−0.42
GO:0001893	Maternal placenta development	3	−2.84	−0.41
R-HSA-6806003	Regulation of TP53 expression and degradation	3	−2.77	−0.37

DEGs, differentially expressed genes; ceRNA, competing endogenous RNA.

Establishment of the PPI network and hub gene subnetwork

PPI network was constructed that contained 93 nodes and 107 edges. According to the network, there were 9 major hub genes (WD and tetratricopeptide repeats 1 (WDTC1), heat shock protein family A (Hsp70) member 5 (HSPA5), N-ethylmaleimide sensitive factor, vesicle fusing ATPase (NSF), fibroblast growth factor 2 (FGF2), cyclin dependent kinase inhibitor 1A (CDKN1A or P21), bone morphogenetic protein 2 (BMP2), complement C3 (C3), GM2 ganglioside activator (GM2A) and discs large MAGUK scaffold protein 4 (DLG4)], and the gene with the highest degree (degree=7), can be observed in Fig. 7A and B, indicating that these genes may be the key genes that can eventually elicit a significant regulatory role in these networks. In other words, they belong to the lncRNA-miRNA-hub gene network.

Figure 7.

(A) Establishment of the PPI network for DEGs involved in the ceRNA network and (B) hub gene sub-network. PPI, protein-protein interaction; DEG, differentially expressed genes; ceRNA, competing endogenous RNA; WDTC1, WD and tetratricopeptide repeats 1; HSPA5, heat shock protein family A (Hsp70) member 5; NSF, N-ethylmaleimide sensitive factor, vesicle fusing ATPase; FGF2, fibroblast growth factor 2; CDKN1A or P21, cyclin dependent kinase inhibitor 1A; BMP2, bone morphogenetic protein 2; C3, complement C3; GM2A, GM2 ganglioside activator; DLG4, discs large MAGUK scaffold protein 4.

Kaplan-Meier survival analysis

There were 4 hub genes (BMP2, CDKN1A, GM2A, and DLG4) that were treated as prognostic factors by the Kaplan-Meier plotter analysis tool (P<0.05). However, the remaining genes were not related to prognosis. High expression levels of BMP2 (Fig. 8B) [hazard ratio (HR)=4.37, 95% CI: 1.43–7.95, P=0.0032], GM2A (Fig. 8D) (HR=4.39, 95% CI: 1.03–18.69, P=0.029), and CDKN1A (Fig. 8A) (HR=2.34, 95% CI: 1.03–5.34, P=0.038) were related to poor OS in ESCC patients, while low expression of DLG 4 (Fig. 8C) (HR=0.29, 95% CI: 0.09–0.98, P=0.035) was significantly associated with a favorable prognosis in ESCC patients.

Figure 8.

Kaplan-Meier survival analysis was used to validate hub genes. (A) CDKN1A, (B) BMP2, (C) GM2A and (D) DLG4. CDKN1A or P21, cyclin dependent kinase inhibitor 1A; BMP2, bone morphogenetic protein 2; GM2A, GM2 ganglioside activator; DLG4, discs large MAGUK scaffold protein 4; HR, hazard ratio.

Construction of a ceRNA network of 4 prognostic hub genes

To identify prognostic ceRNA networks, we constructed a prognostic hub gene/lncRNA/miRNA network (Fig. S3), including 4 hub genes, 13 miRNAs, and 10 lncRNAs.

Discussion

Esophageal squamous cell carcinoma (ESCC) is a malignant tumor that seriously threatens human health (30,31). To date, no biomarkers have been wholly proven to predict the survival of patients with ESCC. Therefore, to improve prognosis, it is imperative to look for reliable biomarkers that identify high risk patients with ESCC. Our previous research confirmed that phospholipase C epsilon 1 (PLCE1) is associated with poor prognosis in patients with ESCC (13) and can serve as an oncogene by improving cell proliferation and repressing apoptosis. Additionally, growing evidence has indicated that competitive endogenous RNA (ceRNA) mechanisms play a critical role in cancer initiation and progression. However, the underlying mechanisms by which PLCE1-regulated ceRNAs function in ESCC are still largely elusive.

Here, we constructed a PLCE1-related ceRNA network obtained from microarray. We also predicted the biological functions and pathways of DEGs involved in the ceRNA network by GO and KEGG pathway analysis. PPI networks among the DEGs from the ceRNA were constructed using the STRING database, and hub genes were screened from the PPI network. Finally, PLCE1-related lncRNA/miRNA/hub gene axes were obtained from the ceRNA network, which could potentially contribute to personalized prognosis and new therapies for ESCC patients.

Numerous studies have recently shown that the expression of lncRNAs is dysregulated in ESCC and represents opportunities for cancer prognosis, diagnosis, and therapeutics, suggesting that lncRNAs could act as tumor-associated biomarkers. For example, Liang et al showed that CASC9 is highly increased in ESCC tissues and cell lines and is associated with overall survival (OS) (32). More importantly, increasing evidence has uncovered that lncRNAs can serve as miRNA sponges during ESCC onset and development. For instance, Li et al (33) revealed that regulation of ZFAS1 by overexpression of miRNA-124 may mediate ESCC proliferation, invasion, migration, and apoptosis. In the present study, a total of 12 lncRNAs (SNORA9, SNORD3D, HCP5, LENG8-AS1, HCG15, CEACAM22P, SHC1P1, HOTAIRM1, TSPY25P, DNAJC3-AS1, ABCC6P2, and RAB43P1) were found to be involved in the ceRNA network. Among them, HCP5 serves as a ceRNA to negatively regulate miR-203 and enhance SNAI expression, thus it is a candidate therapeutic target for lung adenocarcinoma (34). Additionally, a previous report indicated that HOTAIRM1 functions as a ceRNA to regulate the expression of miR-129-5p and miR-495-3p in glioma progression and is positively correlated with OS in glioma patients (35). However, in our analysis, the PLCE1 expression level had a negative correlation with HOTAIRM1, but a positive relationship with HCG15 in Eca109 cells, which has been effectively validated by ESCA tissues based on the GEPIA (36) (Fig. S4A and B).

miRNAs, as a class of non-coding single-stranded RNA molecules encoded by endogenous genes with a length of about 22 nucleotides, are involved in regulating the degradation of targeted mRNAs and inhibiting mRNA translation (37), which allows them to further regulate cancer proliferation, differentiation, apoptosis, and migration (38). In our previous research, 14 DEmiRNAs in the ceRNA network were identified, among which hsa-miR-17-5p, hsa-miR-22-3p, and hsa-miR-1297 have been reported to be involved in tumorigenesis (39–41); hsa-miR-301b-3p and hsa-miR-455-5p have also been verified to be highly expressed in ESCA from the TCGA database (23) (Fig. S5).

To further identify the role of the ceRNA network, we established a PPI network that identified nine hub genes [WDTC1, HSPA5, NSF, FGF2, CDKN1A (P21), BMP2, C3, GM2A, and DLG4], indicating that PLCE1 could regulate these genes by a ceRNA mechanism. In this study, the PLCE1 expression level had a strongly negative relationship with HSPA5, NSF, CDKN1A, BMP2, GM2A, C3 and DLG4, while it had a positive correlation with WDTC1 and FGF2 in Eca109 cells, Our previous si-PLCE1 chip results in EC9706 also indicated that the PLCE1 expression level had an inverse relationship with HSPA5, NSF, FGF2, CDKN1A, BMP2, C3, and DLG4, while it correlated positively with WDTC1 and GM2A (Table II) (42), suggesting that their results are highly consistent with each other. Numerous reports have also demonstrated that CDKN1A (P21) and BMP2 (43) in many cancers can inhibit the cell cycle at the G1 and G2 phases to elicit tumor suppression (44–47). Intriguingly, our results also confirmed that PLCE1 inhibition may lead to an increase in CDKN1A (P21) and BMP2. And more importantly, CDKN1A-associated ceRNA network includes CDKN1A/hsa-miR-22-3p/HOTAIRMI and CDKN1A/hsa-miR-206-5p/HOTAIRMI; BMP2-associated ceRNA network harbors BMP2/hsa-miR-17-5p/HOTAIRMI, BMP2/hsa-miR-20b-5p /HOTAIRMI, BMP2/hsa-miR-17-5p/HCG15 and BMP2/hsa-miR-20b-5p/HCG15 in our research, revealing that PLCE1 triggers ESCC progression by decreasing CDKN1A (P21) and BMP2 expression via ceRNA. Functionally, DNA damage response, signal transduction by p53 class mediator resulting in transcription of p21 (GO: 0006978) and negative regulation of cortisol biosynthetic process (GO: 2000065) were enriched based on these two hub genes by means of Metascape (28). However, Kaplan-Meier plotter analysis showed that CDKN1A and BMP2 afforded poor OS in ESCC patients. One obvious reason for this is the difference between esophageal squamous cell carcinoma cells and tissue. In addition, too few samples were used for survival analysis. Collectively, these findings indicate that PLCE1 triggers ESCC progression by meditating these hub genes via ceRNA.

Table II.

Differentially expressed hub genes in PLCE1 siRNA-treated ESCC cells compared with control cells.

	Gene symbol	log2(Fold change)	Attribute		Gene symbol	Fold change	Attribute
Eca109-associated	WDTC1	−2.64	Down	EC9706-associated	WDTC1	0.78	Down
microarray	HSPA5	1.92	Up	microarray	HSPA5	2.31	Up
	NSF	2.44	Up		NSF	1.16	Up
	FGF2	−2.29	Down		FGF2	2.22	Up
	CDKN1A	3.86	Up		CDKN1A	3.27	Up
	BMP2	1.85	Up		BMP2	1.87	Up
	GM2A	2.25	Up		GM2A	0.73	Down
	C3	2.24	Up		C3	1.45	Up
	DLG4	2.91	Up		DLG4	1.14	Up

ESCC, esophageal squamous cell carcinoma; PLCE1, phospholipase C epsilon 1; WDTC1, WD and tetratricopeptide repeats 1; HSPA5, heat shock protein family A (Hsp70) member 5; NSF, N-ethylmaleimide sensitive factor, vesicle fusing ATPase; FGF2, fibroblast growth factor 2; CDKN1A or P21, cyclin dependent kinase inhibitor 1A; BMP2, bone morphogenetic protein 2; C3, complement C3; GM2A, GM2 ganglioside activator; DLG4, discs large MAGUK scaffold protein 4.

To date, most research involved in ceRNAs in cancer development has comprised lncRNAs, miRNAs and coding genes. There are few reports investigating whether other genes regulate ceRNAs in cancer. However, Sun's group (48) found that TGF-β activates lncRNA-ATB and upregulates ZEB1 or ZEB2 by means of competitive binding of miR-200; thereby inducing epithelial-mesenchymal transition in vivo and in vitro in hepatocellular carcinoma. In our study, PLCE1 was also found to be involved in the regulation of ceRNAs to promote the development of ESCC. Despite revealing the relationship between PLCE1 and ceRNAs in ESCC, for the first time, there are some limitations to this investigation. The results were obtained from microarray data and have not yet been verified by experimental means. In addition, a limited number of samples were used for the bioinformatics analysis.

In conclusion, for the first time, we comprehensively revealed the relationship between PLCE1 and ceRNAs in ESCC and constructed a PLCE1-regulated ceRNA network, which highlights the role of PLCE1 in promoting ESCC. Consequently, PLCE1 is a promising biomarker and potential novel therapeutic target for ESCC.