Identification and validation of hub microRNAs dysregulated in esophageal squamous cell carcinoma

Screening and identification of differentially expressed miRNAs and target genes in ESCC

To identify differentially expressed miRNAs (DEMs) in ESCC, data from two independent miRNA expression arrays (GSE114110 and GSE43732) [12, 13] were downloaded from the GEO database and normalized using the limma software package of R language (Figure 1A–1D). According to the thresholds set (p < 0.05 and log2FC ≥ 1), 277 and 68 DEMs were found in GSE114110 and GSE43732, respectively. Subsequently, we obtained the top 20 DEMs in each dataset (Tables 1 and 2 and Figure 1E and 1F), and identified through Venn diagram analysis three upregulated and five downregulated DEMs held in common between the two datasets (Figure 1G and 1H). Based on these eight DEMs, miRNA-target gene interactions were evaluated using the experimentally validated miRTarBase database [14]. A total of 468 and 753 possible target genes, respectively, were thus identified for the three upregulated and five downregulated hub miRNAs.

Figure 1. Identification of differently expressed miRNAs (DEMs) in esophageal squamous cell carcinoma (ESCC). (A, B) GSE114110 data before and after normalization. (C, D) GSE43732 data before and after normalization. (E, F) Volcano plots of DEMs in GSE114110 and GSE43732, respectively. Black dots represent genes equally represented between ESCC and normal samples. Red and green dots represent upregulated and downregulated miRNAs, respectively. Volcano plots showing all DEMs. |log2FC| ≥ 1 and P < 0.05 were set as cut-off criteria. (G) Venn diagram analysis showing the top 10 upregulated and downregulated miRNAs in the two GEO datasets. (H) Identification of three upregulated and five downregulated miRNAs overlapping between both GEO datasets.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of hub miRNAs’ target genes

GO functional annotation analysis results indicated that for targets of the three upregulated miRNAs the most enriched terms were ‘regulation of transcription from RNA polymerase II promoter’ and ‘positive regulation of transcription’ in the biological process (BP) category, ‘focal adhesion’ and ‘ficolin-1-rich granule lumen’ in the cellular component (CC) category, and ‘RNA and DNA binding’ in the molecular function (MF) category (Figure 2A–2C). In turn, the target genes of the five downregulated miRNAs were mostly enriched in ‘negative regulation of apoptosis’ and ‘programmed cell death’ in the BP category, ‘focal adhesion’ and ‘actin cytoskeleton’ in the CC category, and ‘cadherin binding’ and ‘RNA binding’ in the MF category (Figure 2E–2G). KEGG pathway analysis showed that for targets of the upregulated miRNAs the most enriched pathways were ‘microRNAs in cancer’, ‘phosphoinositide 3-kinase (PI3K)-AKT signaling pathway’, and ‘cell cycle’ (Figure 2D). For targets of the downregulated miRNAs, the most enriched pathways were ‘cancer’, ‘focal adhesion’, and ‘proteoglycans in cancer’ (Figure 2H).

Figure 2. GO and KEGG analysis of genes targeted by DEMs commonly shared between the two GEO datasets. (A–C) Top 10 GO biological process (BP), cellular component (CC), and molecular function (MF) terms enriched in target genes of the three upregulated miRNAs. (D) Top 10 KEGG pathways enriched in target genes of the three upregulated miRNAs. (E–G) Top 10 GO BP, CC, and MF terms enriched in target genes of the five downregulated miRNAs. (H) Top 10 KEGG pathways enriched in target genes of the five downregulated miRNAs.

Construction of target gene-PPI and miRNA-hub gene networks

To identify hub genes among the targets of the eight DEMs defined above, PPI networks were constructed using the STRING database and Cytoscape software [15]. The top 10 hub genes thus identified are shown in Figure 3A and 3B. For the three upregulated miRNAs, the predicted hub genes were MYC, CCND1, HSP90AA1, PTEN, MAPK1, NOTCH1, CDH1, CASP3, HSPA4, and ACTB. For the five downregulated miRNAs, the predicted hub genes were AKT1, CCND1, CDC42, IL6, FN1, MAPK1, JUN, EGFR, ACTB, and HRAS. To verify potential importance in the development of ESCC, GO and KEGG analyses were also performed on these gene clusters. Results showed that these genes were significantly enriched in various GO functions, such as ‘positive regulation of cellular process’ and ‘positive regulation of cell proliferation’ in the BP category (Figure 3C), ‘cytoskeleton’ and ‘focal adhesion’ in the CC category (Figure 3D), and ‘purine ribonucleoside triphosphate binding’ and ‘protein kinase binding’ in the MF category (Figure 3E). These genes were also enriched in prominent cancer-related signaling pathways (Figure 3F). Subsequently, a miRNA-hub gene network was constructed using Cytoscape software. Interactional analysis indicated that among the 10 hub genes targeted by the upregulated miRNAs, six (MYC, ATCB, CCND1, MAPK1, PTEN, and CASP3) are potentially regulated by miR-196a-5p (Figure 4A), while another six (IL6, CCND1, ATCB, CDC42, EGFR, and FN1) are potentially regulated by miR-1-3p (Figure 4B). Utilizing the miRNACancerMAP database, we surveyed cellular signaling pathways involving miR-196a-5p and miR-1-3p. Results showed association of these miRNAs with signaling cascades closely related to tumor growth, including MAPK, adherens junction, mechanistic target of rapamycin (mTOR), and PI3K-AKT pathways (Figure 4C and 4D). Since these results overlapped with those obtained in our pathway analysis, we deemed that miR-196a-5p and miR-1-3p may be essential miRNAs in ESCC pathogenesis.

Figure 3. GO and KEGG analysis of hub-gene targets. (A) Top 10 hub-gene targets for the three upregulated miRNAs. (B) Top 10 hub-gene targets for the five downregulated miRNAs. (C–E) Top 10 GO BP, CC, and MF terms enriched in the top 20 hub-gene targets. (F) Top 10 KEGG pathways enriched in the 20 hub-gene targets.

Figure 4. Regulatory networks of miRNAs-hub-gene targets. (A) Interaction network of the three upregulated miRNAs and their hub-gene targets. (B) Regulatory network of the five downregulated miRNAs and their hub-gene targets. (C–D) Regulatory networks of miR-196a-5p and miR-1-3p and associated signaling pathways.

Subsequently, we accessed the UALCAN tool to assess within TCGA the expression of the most relevant target genes of miR-196a-5p and miR-1-3p in esophageal carcinoma [16]. Results indicated that the expression of four target genes (CCND1, CASP3, EGFR, and CDC42) was markedly upregulated in esophageal carcinoma (Figure 5B, 5E, 5G, and 5H). Since sampling size was relatively limited, further verification may still reveal differential expression of other hub genes regulated by miR-196a-5p and miR-1-3p in ESCC.

Figure 5. Relative expression of hub-gene targets of miR-196a-5p and miR-1-3p. (A–E) Relative expression of MYC, CCND1, PTEN, MAPK1, and CASP3 mRNAs in esophageal carcinoma, compared to normal esophageal tissue samples. (F–J) Relative expression of ACTB, EGFR, CDC42, IL6, and FN1 mRNAs in esophageal carcinoma, compared to normal esophageal tissue samples. Analysis of esophageal carcinoma RNA-seq datasets in TCGA was performed on the UALCAN database.

Prognostic value of miR-196a-5p and miR-1-3p expression in ESCC

To complement the above bioinformatics analyses, we verified the expression of miR-196a-5p and miR-1-3p in ESCC cell lines, the TCGA database, and clinical ESCC cases. Using qPCR, we determined that either upregulation or downregulation of these two hub DEMs occurred in different ESCC cell lines (KYSE30, KYSE140, KYSE410, KYSE150, KYSE510, Eca109, and TE-1), compared to normal esophageal epithelial (Het1A) cells (Figure 6A and 6D). Analysis of miRNA expression profiles from TCGA showed that the expression of miR-196a-5p was significantly higher, whereas that of miR-1-3p was significantly lower, in ESCC samples compared to normal esophageal tissue (Figure 6B and 6E). In addition, we collected 32 ESCC samples from our institution and summarized their clinicopathological data. Results from qPCR assays confirmed that both miR-196a-5p and miR-1-3p were differentially expressed in ESCC samples (Figure 6C and 6F). Analysis of the relationship between miR-196a-5p and miR-1-3p expression and clinicopathological features showed that the expression of these two miRNAs was associated with T classification in ESCC patients (Supplementary Tables 1–3). Meanwhile, survival analysis of TCGA data on the starBase database [17], showed that low miR-196a-5p expression was associated with better prognosis in patients with esophageal carcinoma. In contrast, no prognostic significance was found for miR-1-3p (Figure 6G and 6H).

Figure 6. Expression analysis and prognostic value of miR-196-5p and miR-1-3p in ESCC. (A) Analysis of miR-196a-5p expression by qPCR in ESCC and normal esophageal epithelial (Het1A) cells. (B) Expression of miR-196a-5p in ESCC and normal esophagus samples from TCGA. (C) Expression of miR-196a-5p in ESCC and paired normal tissue samples (n = 32) collected at our institution. (D) Expression of miR-1-3p in ESCC and Het1A cells. (E) Expression of miR-1-3p in ESCC and normal esophagus samples from TCGA. (F) Expression of miR-1-3p in ESCC samples and paired normal tissue samples (n = 32) collected at our institution. (G and H) Kaplan-Meier survival analysis based on miR-196a-5p and miR-1-3p expression in esophageal carcinoma samples from TCGA. *P < 0.05, **P < 0.01. ***P < 0.001, and ****P < 0.0001.

miR-196a-5p and miR-1-3p exert opposing roles in ESCC cell proliferation and migration

Next, we sought to elucidate the functional effects of miR-196a-5p and miR-1-3p dysregulation through studies on ESCC cells in vitro. To this end, we transfected ESCC cell lines with miR-196a-5p and miR-1-3p mimics or their negative controls (NC). Considering transfection efficiency and background expression of the two miRNAs, we upregulated miR-196a-5p in KYSE150 and KYSE410 cells and miR-1-3p in KYSE30 and KYSE410 cells. Subsequently, cell proliferation was measured using Cell Counting Kit-8 (CCK-8) and 5-ethynyl-2-deoxyuridine (EdU) assays. As shown in Figure 7A–7C (CCK-8 results), Figure 7D and 7F (EdU results), compared to the corresponding NCs, upregulation of miR-196a-5p promoted ESCC cell proliferation whereas upregulation of miR-1-3p showed inhibitory effects. We next assessed the impact of miR-196a-5p and miR-1-3p upregulation on the migrational potential of ESCC cells using Transwell migration assays. Results showed that upregulation of miR-196a-5p enhanced, whereas upregulation of miR-1-3p markedly suppressed, the migrational ability of the ESCC cell lines tested (Figure 7E and 7G).

Figure 7. Opposing effects of miR-196a-5p and miR-1-3p on proliferation and migration in cultured ESCC cells. (A–C) CCK-8 cell proliferation assay results from KYSE150, KYSE30, and KYSE410 cells transfected with mimics of miR-196a-5p, miR-1-3p, or negative control (NC). (D, F) EdU cell proliferation assay results from KYSE150, KYSE30, and KYSE410 cells transfected with mimics of miR-196a-5p, miR-1-3p, or NC. (E, G) Transwell migration assay results from KYSE150, KYSE30, and KYSE410 cells transfected with miRNA mimics. *P < 0.05, **P < 0.01, and ***P < 0.001.

Data collection and study design

Two ESCC miRNA expression microarray datasets (GSE114110 and GSE43732) were retrieved from the GEO repository (https://www.ncbi.nlm.nih.gov/geo). GSE114110 included 30 ESCC and 10 normal esophageal epithelial samples [12], while GSE43732 included 119 ESCC samples and paired adjacent normal tissues [13]. All samples were collected in China and originated from patients who had not received anticancer treatment. The above datasets were produced independently using the GPL24967 and GPL16543 platforms, respectively. Therefore, we normalized these data using R’s limma package.

DEM screening and Venn diagram analysis

After obtaining standardized chip data from the GSE114110 and GSE43732 datasets, differential miRNA expression in ESCC was investigated in relation to normal esophageal tissue (control) using the limma software package in R. We set P < 0.05 and |log2 fold change (FC)| ≥ 1 as the thresholds to identify DEMs. Overlapping miRNAs in the two datasets were screened using Venn diagram analysis with the VennDiagram R package [37].

Prediction of target genes and miRNA pathway analysis

MiRTarBase, an experimentally validated database of miRNA-target interactions, was used to predict potential target genes for the eight overlapping DEMs. MiRNACancerMAP database was used to perform miRNA pathway analysis.

GO and KEGG pathway analyses

Functional and pathway analyses of the predicted target genes of the eight overlapping DEMs were conducted using the GO and KEGG databases, respectively, by processing data with the DAVID online tool (https://david.ncifcrf.gov/) [38]. P < 0.05 was considered significant.

Construction of target gene-PPI and miRNA-gene networks

To assess the functional associations among the target genes of upregulated and downregulated DEMs, we uploaded target gene data to the STRING database. Interactions with a combined score > 0.4 were considered significant. Highly interconnected (hub) genes in the PPI network were analyzed using Cytoscape software (version 3.6.0). After hub DEMs and hub target genes were identified, Cytoscape was used to visualize the resulting miRNA-gene network.

Differential gene expression and survival analyses

After extracting expression data for miR-196a-5p and miR-1-3p from ESCC and normal esophageal epithelial tissues from the TCGA database, analysis of differential expression of target genes was performed on the University of Alabama Cancer Database (UALCAN; http://ualcan.path.uab.edu/) resource. The StarBase database was used to evaluate the prognostic values of hub miRNAs in esophageal carcinoma through Kaplan-Meier plots.

Cell lines and clinical samples

Human ESCC KYSE30, KYSE140, KYSE410, KYSE150, KYSE510, Eca109, and TE-1 and normal esophageal epithelial (Het1A) cell lines were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS; BI, Israel), 100 U/mL penicillin, and 100 μg/mL streptomycin at 37°C in an incubator with a humidified atmosphere containing 5% CO₂.

Thirty-two ESCC and paired normal adjacent tissue samples were acquired from patients undergoing surgical procedures at the Third Affiliated Hospital of Soochow University. Patient baseline characteristics were collected from medical records. All samples were obtained after written informed consent was provided, in accordance with the Code of Ethics of the World Medical Association (Declaration of Helsinki). All protocols for the use of patient samples in this study were approved by the Ethics Committee of the Third Affiliated Hospital of Soochow University.

Cell transfection

Negative control (NC), miR-196a-5p, and miR-1-3p mimics were purchased from GenePharma (Suzhou, China), and transfected into ESCC cell lines using siRNA-Mate (GenePharma) according to the manufacturer’s instructions.

RNA extraction and quantitative PCR

TRIzol reagent was used to extract RNA from clinical samples and ESCC cell lines. The miRNAs first-strand cDNA and q-PCR kits were purchased from TIANGEN (Beijing, China) and U6 snRNA was used as the internal reference. The forward primers used were: miR-196a-5p: 5'-CGCGTAGGTAGTTTCATGTTGTTGGG-3, miR-1-3p: 5'-GCGCGCTGGAATGTAAAGAAGTATGTAT-3', and U6 snRNA: 5'-CTCGCTTCGGCAGCACA-3'. The reverse primers were provided in the qPCR detection kit. Relative miRNA expression levels were calculated using the 2^-ΔΔCt method.

Cell proliferation assay

The CCK-8 assay (Dojindo, Kumamoto, Japan) was used to quantify cell proliferation. Twenty-four h after miRNA mimic transfection, cells were seeded into 96-well culture plates (1.5 × 10³ cells/well). At different time points, 10 μL CCK-8 solution was added to each well, incubated for 2 h, and samples’ optical density (OD) was measured immediately at 450 nm. The Cell-Light EdU Apollo567 In Vitro Kit (Cat.C10310-1, Ruibo Biotech., Guangzhou, China) was also used to assess cell proliferation. To this end, cells transfected with miRNA mimics for 24 h were seeded into 96-well plates (1.5 × 10⁴ cells/well). After 24 h, the cells were incubated in EdU working solution for 2 h, fixed with 4% paraformaldehyde, permeabilized, washed, and stained with 1x Apollo solution according to the manufacturer’s instructions. Results were analyzed from microphotographs taken using a fluorescence microscope.

Transwell migration assay

Transwell inserts were used to conduct migration assays. After transfection, cells were resuspended in serum-free culture solution and approximately 3 × 10⁴ cells were loaded into the upper Transwell chamber. The lower chamber was filled with 600 μL of culture medium containing 20% FBS. After 24 h, the cells that remained on the upper chamber were removed and those that migrated through the membrane were fixed in 4% paraformaldehyde for 30 min, stained with 0.1% crystal violet, and counted in five randomly selected fields under an inverted microscope.

Statistical analysis

GraphPad Prism 7 software was used to perform unpaired Student’s t-tests to analyze differences between two groups. P < 0.05 was regarded as significant.