Network analysis of miRNA targeting m6A-related genes in patients with esophageal cancer

Data from TCGA and GEO datasets

We downloaded m6A-related gene expression profiling data for 162 esophageal cancer samples and 11 normal samples from the TCGA repository (https://portal.gdc.cancer.gov/). We predicted the upstream miRNAs for each m6A-related gene using the online tool starBase3.0 (http://starbase.sysu.edu.cn/). miRNA expression profiles for 185 esophageal cancer samples and 13 normal samples were also obtained from TCGA. Patient inclusion criteria included: (1) the availability of complete survival data and (2) an esophageal cancer histology type. Patients with a survival time of less than 100 days were excluded from our study. We enrolled 164 esophageal cancer patients in the survival analysis. We included 107 cases with full clinicopathological data for further Cox proportional hazard regression analysis. miRNA microarray expression profiling data and the clinical data for 119 patients with esophageal cancer were obtained from the GEO microarray dataset (GSE43732) as validation cohort data.

Screening of differentially expressed genes

We analyzed the differentially expressed m6A-related genes and miRNAs between esophageal cancer samples and matched healthy tissues using the “limma” package in R (Ritchie et al., 2015). The T-test in the limma package was used to calculate the differentially expressed genes (DEGs) in the senescent and quiescent. A P-value < 0.05 was considered to be significant.

Identification of prognosis-related genes

The univariate cox regression test and Kaplan–Meier (K–M) survival analysis were performed to screen out the potential overall survival (OS) relevant genes. A P-value less than 0.05 was considered to be statistically significant. We used multivariate cox analysis to evaluate the contribution of the selected genes. The analysis was conducted using the R package of survival.

Construction and evaluation of the miRNA-based prognostic signature

We used the “glmnet” package to perform LASSO analysis (https://CRAN.R-project.org/package=glmnet) on the selected miRNAs. The prognostic miRNA-based signature was designed using a risk scoring method with the following formula: $$\mathrm{R}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum _{\mathrm{i}=1}^{\mathrm{N}}(\mathrm{e}\mathrm{x}\mathrm{p}\ast \mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f})$$

in which exp was the expression value of the gene and coef was the coefficient of miRNA generated from the LASSO-Cox regression analyses. The median risk score was set as the cutoff point, according to which the esophageal cancer patients were divided into low-risk and high-risk groups. We then measured the survival difference between the high-risk and low-risk groups stratified by the OS-related miRNA prediction model. The area under the curve (AUC) of the receiver operating characteristic (ROC) curve was calculated to assess the efficacy and accuracy of the OS classifier. We conducted the univariate and multivariate Cox proportional hazards analyses to compare the relative prognostic value of the miRNA-based signature with that of different clinicopathological features.

Cell culture and transfection

Human normal esophageal epithelial cell line (HEEC) and esophageal carcinoma cell lines (KYSE30 and TE-1) were purchased from American Type Culture Collection (ATCC, Manassas, VA, USA). Both cell lines were cultured in RPMI 1640 medium (Gibco, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS, Belize, Gibco, Waltham, MA, USA), 100 U/ml penicillin and 100 mg/ml streptomycin (Invitrogen, Carlsbad, CA, USA). Mir-186 mimics, mimic-NC, sh-HNRNPC (sh-HNRNPC#1 and sh-HNRNPC#2) and sh-NC were constructed by RiboBio (Guangzhou, China). Esophageal carcinoma cell lines were placed into six-well plates until cell confluence was approximately 70%. The miR-186 mimics, mimic-NC, sh-HNRNPC, or sh-NC were transfected into esophageal carcinoma cell lines according to the Lipofectamine 2000 transfection reagent (Invitrogen, Carlsbad, CA, USA).

Quantitative real-time PCR (qRT-PCR)

Total RNA was extracted from cells using TRIzol reagent (Takara, Japan). Messenger RNA (mRNA) and miRNA were reverse-transcribed into complementary DNA with the cDNA Synthesis Kit (ABM, Vancouve, Canada). qRT-PCR was performed with a SYBR Green PCR kit (Takara, Kusatsu, Japan) and the Bio-Rad Real-Time PCR System (Bio-Rad, Hercules, CA, USA). GAPDH was treated as the internal control gene for HNRNPC, which U6 was the endogenous control gene for miRNAs. Primer information is presented in Table 1. The relative expression level was analyzed using the 2^−ΔΔCt. All assays were performed in triplicate.

MTT assay

Cell proliferation was evaluated by the 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-tetrazolium bromide (MTT) assay. Cells were plated in 96-well plates at a density of 4,000 cells per well in 100 μl of complete medium and cultured for 24 h. The MTT solution (5 mg/ml, 20 μl) was added to each well and incubated for 4 h at 37 °C. The media was then removed and 100 μl DMSO was added and thoroughly mixed for 15 min. The relative number of surviving cells was assessed by measuring the optical density (OD) of cell lysates at 490 nm. Each group contained five wells. All assays were performed in triplicate.

Cell migration and invasion assays

The ability of cell migration and invasion was detected with polycarbonate membrane Boyden chambers in a transwell apparatus (Costar, Irvine, CA, USA). We evenly covered the upper surface of the membrane with 50 μL of 1 mg/ml Matrigel (Sigma, St. Louis, MO, USA), incubated it for 1 h at 37 °C, and hydrated it in fetal bovine serum (FBS) for 2 h before use for the invasive assay. The cells were digested and added into the top chamber at the cell density of 3 × 10⁴ cells/ml. The lower chamber was filled with cell culture medium containing 20% FBS and was incubated at 37 °C for 36 h. After incubation, the invasive cells were fixed, stained, and counted under a microscope. In the cell migration assay, the same procedures were conducted without Matrigel on the membrane. Three independent experiments were performed.

Luciferase reporter assay

We inserted the 3′-UTR of wild-type HNRNPC (HNRNPC-WT) with mutations in the putative binding site downstream of the frefly luciferase reporter into the psiCHECK-2 vector (Promega, Madison, WI, USA). The corresponding mutant construct was created by mutating the seed regions of the miR-186 binding site. This site was named 3′-UTR HNRNPC-MUT. HEK293 cells were co-transfected with miR-186 mimics or the NC control, then combined with either HNRNPC-WT or HNRNPC-MUT reporter vectors for 48 h. Cells were harvested using the dual Glo™ Luciferase Assay System (Promega, Madison, WI, USA) to measure luciferase activity. Renilla activity was used as an internal control. Experiments were conducted independently at least three times.

Statistical analysis

Data were shown as mean ± standard deviation (SD). K–M survival analysis and univariate/multivariate Cox proportional hazards analyses were used to compare the two groups of patients. HR and 95% confidence intervals (CIs) were calculated. All statistical analyses were performed by SPSS 22.0 software (IBM Corp., Armonk, NY, USA). P < 0.05 was considered to be statistically significant.

The expression pattern of m6A-related genes in esophageal cancer

To explore the expression patterns of m6A-related genes in esophageal cancer, we first extracted and analyzed the expression data of 13 known m6A-related genes from the TCGA database. HNRNPC and YTHDF1 were significantly upregulated in tumor tissues compared with corresponding normal tissues, while ZC3H13, YTHDC2 and METTL14 notably decreased (Figs. 1A and 1B). These results suggested that m6A modification may affect esophageal cancer development.

Prognostic analysis of m6A-associated genes in esophageal cancer

The main clinicopathological characteristics of esophageal cancer patients are summarized in Table 2. We analyzed the correlation between m6A-related gene expression and clinical follow-up data using K–M survival analysis to investigate the prognostic role of m6A-associated genes in esophageal cancer. Patients were divided into the high- or low-expression groups based on the best cut-off value. The results showed that the low expression of ALKBH5 and the high expression of HNRNPC or WTAP were the best predictors for poor OS rates in esophageal cancer patients (Figs. 2A–2C).

We then performed univariate and multivariate Cox regression analyses to explore whether these genes were associated with the prognosis of esophageal cancer patients. The results of univariate cox regression revealed that ALKBH5 (P = 0.001) was a protective gene for esophageal cancer (Fig. 2D). Multivariate Cox regression revealed that both ALKBH5 (P = 0.006) and METTL14 (P = 0.015) were protective genes for OS. It also revealed that HNRNPC (HR = 4.171, P = 0.021, 95% CI [1.244–13.986]) may be an independent risk factor for OS (Fig. 2E). We focused on the m6A-related gene HNRNPC and its role in esophageal cancer in light of our analysis of m6A-related gene expression levels and the prognostic values above.

Constructing and validating the miRNA-based prognostic signature

We used starBase 3.0 to assess potential upstream miRNAs targeting m6A-related genes. A total of 522 miRNAs were obtained and the number of m6A-related genes targeted by each miRNA was calculated. The top seven miRNAs with the largest number of target genes were then determined. The relationships between the 13 target genes and their corresponding miRNAs are shown in Fig. 3A. We also explored these miRNAs’ expression patterns in the TCGA dataset. As shown in Fig. 3B, hsa-mir-320b-1 (P = 0.023), hsa-mir-548o (P = 0.036), and hsa-mir-320a (P = 0.042) were upregulated in esophageal cancer tissues compared with normal esophageal tissues.

We further evaluated the prognostic value of these seven miRNAs in esophageal cancer through univariate and multivariate cox regression analysis. The results of univariate cox regression analysis indicated that miR-186 (P = 0.003) had prognostic value (Fig. 3C). Multivariate cox regression demonstrated that miR-186 (P = 0.018) and miR-320c (P = 0.008) may be independent prognostic indicators (Fig. 3D). We used LASSO Cox regression analysis to determine the most prognostic miRNAs among esophageal cancer patients. Finally, four miRNAs (miR-186, miR-320c, miR-320d, miR-320b) were used for the construction of the prognostic signature according to the minimum criteria and the coefficients (Figs. 3E & 3F).

Risk scores for each sample in the TCGA data set were determined using the coefficients derived from the LASSO algorithm. The patients were classified into low- or high-risk groups on the basis of the median risk score. Kaplan–Meier survival curves indicated that patients with higher risk scores (n = 82) had significantly shorter OS than those with lower risk scores (n = 82, p < 0.001, Figu. 3G). A time-dependent ROC curve was performed to evaluate the sensitivity and specificity of the signature. The four-miRNA signature’s AUC values were 0.805, 0.681, and 0.645, respectively, for an OS of 1, 3, and 5 years (Fig. 3H). This indicates that our model could be used to predict esophageal cancer patient survival.

In addition, the accuracy of the risk scoring model was verified using the test set of esophageal cancer cohort in the GEO dataset. We calculated the risk scores for esophageal cancer patients according to the expression levels of miR-186, miR-320c, miR-320d, and miR-320b in the GEO microarray dataset. Patients were divided into high-risk or low-risk groups based on the median risk score of the cohort. Kaplan–Meier analysis showed that patients in the high-risk group had significantly poorer OS than those in the low-risk group (P < 0.001; Fig. 3I). ROC curve analysis also showed that our prognostic signature for 1 year had the best predictive performance, which was consistent with the results of TCGA database (AUC values of 1-year (0.669), 3-year (0.566), and 5-year (0.559); Fig. 3J).

Independence of the miRNA-based signature from other clinicopathological factors

The following clinicopathological characteristics were considered: age, gender, T staging, lymph node status, metastasis, and staging. Univariate and multivariate cox regression were used to detect the influence of the four-miRNA signature on OS. The lymph node status, stage, and risk score in esophageal cancer patients were associated with OS with univariate analysis (Fig. 4A). Multivariate analysis showed that the four-miRNA signature (HR = 1.023, P = 0.011) was an independent prognostic factor for esophageal cancer patients (Fig. 4B).

Construction of key miRNA-m6A related gene network in esophageal cancer

Our results indicated that HNRNPC may be the key m6A-related gene in esophageal cancer. Mir-186, miR-320c, miR-320d, and miR-320b, which were used for the construction of the prognostic signature, may be the key upstream miRNAs in m6A-related genes. We extracted a critical miRNA-m6A related gene network, found that HNRNPC was targeted by these four miRNAs, and then constructed the key miRNA-m6A related gene network (Fig. 5).

Silencing HNRNPC suppressed proliferation, migration, and invasion of esophageal cancer cells

We further evaluated the expression pattern of HNRNPC in esophageal cancer cells and found that the expression level of HNRNPC in the KYSE30 and TE-1 cell lines was higher than that in the normal esophageal carcinoma epithelial cell line HEEC (Fig. 6A). We transfected a short hairpin RNA (shRNA) targeting HNRNPC into KYSE30 and TE-1 cells to explore the function of HNRNPC. The qRT-PCR confirmed the knockdown efficiency of the two specific shRNAs in both cell lines (Fig. 6B). The MTT assay determined that the knockdown of HNRNPC significantly suppressed esophageal cancer cell proliferation (Figs. 6C & 6D). Furthermore, silencing HNRNPC may also inhibit the migration and invasion of esophageal cancer cells (Figs. 6E and 6F).

miR-186 targeted HNRNPC and suppressed HNRNPC expression

We hypothesized that the candidate miRNA would be downregulated in esophageal cancer with a potential binding sequence for the 3′-UTR of HNRNPC. We chose four miRNAs (miR-186, miR-320c, miR-320d, miR-320b) for further validation. miR-186 had low expression levels in both KYSE30 and TE-1 cells compared with HEEC cells (Fig. 7A). Therefore, we chose miR-186 for further studies. The luciferase reporter assay showed that miR-186 mimics significantly reduced the luciferase activity of the reporter gene in HNRNPC-WT compared with that of the negative control (P < 0.001) (Fig. 7B). The miR-186 mimics’ regulatory effect was inhibited when the predicted binding site in HNRNPC was mutated. KYSE30 or TE-1 cells treated with miR-186 mimics also decreased HNRNPC expression (Figs. 7C and 7D). These results together indicated that HNRNPC was a direct target of miR-186.