Upregulated expression of MMP family genes is associated with poor survival in patients with esophageal squamous cell carcinoma via regulation of proliferation and epithelial‑mesenchymal transition.

Matrix metalloproteinases (MMPs) are involved in the cleavage of several components of the extracellular matrix and serve important roles in tumor growth, metastasis and invasion. Previous studies have focused on the expression of one or several MMPs in esophageal squamous cell carcinoma (ESCC); however, in the present study, the transcriptomics of all 23 MMPs were systematically investigated with a focus on the prognostic value of the combination of MMPs. In this study, 8 overlapping differentially expressed genes of the MMP family were identified based on data obtained from Gene Expression Omnibus and The Cancer Genome Atlas. The prognostic value of these MMPs were investigated; the receiver operating characteristic curves, survival curves and nomograms showed that the combination of 6 selected MMPs possessed a good predictive ability, which was more accurate than the prediction model based on Tumor‑Node‑Metastasis stage. Gene set enrichment analysis and gene co‑expression analysis were performed to investigate the potential mechanism of action of MMPs in ESCC. The MMP family was associated with several signaling pathways, such as epithelial‑mesenchymal transition (EMT), Notch, TGF‑β, mTOR and P53. Cell Counting Kit‑8, colony formation, wound healing assays and western blotting were used to determine the effect of BB‑94, a pan‑MMP inhibitor, on proliferation and migration of ESCC cells. BB‑94 treatment decreased ESCC cell growth, migration and EMT. Therefore, MMPs may serve both as diagnostic and prognostic biomarkers of ESCC, and MMP inhibition may be a promising preventive and therapeutic strategy for patients with ESCC.

Introduction

Esophageal cancer (EC) is the seventh most common type of cancer and the sixth highest cause of cancer-associated death (1). Based on the histopathological appearance, EC is primarily classified into two types, esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC) (2). ESCC is the most general type of EC, which accounts for ~90% of all EC patients in China, whereas EAC is more frequently observed in developed countries (3). Treatment of ESCC includes surgery, radiotherapy and chemotherapy, and have drastically improved survival. However, the five-year survival in patients with ESCC remains relatively low, primarily due to delayed diagnosis (4). Squamous cell carcinoma antigen and CYFRA 21-1 are tumor markers commonly used for managing patients with ESCC; however they have limited application for detection of early stage ESCC due to their lack of specificity and accuracy (5,6). Therefore, the identification of novel molecular biomarkers may assist in the development of novel diagnostic and therapeutic strategies for ESCC.

Matrix metalloproteinases (MMPs) are considered to serve vital roles in early carcinogenic events, tumor growth, invasion, metastasis and tumor-induced angiogenesis (7). The MMP family of proteins are zinc-dependent endopeptidases, and to date, 23 members have been identified in humans (8). MMPs participate in the cleavage of extracellular matrix (ECM) components, such as collagen and elastin (9). MMP dysregulation is associated with diagnostic and prognostic significance in breast, ovarian and colon cancer (10,11). Based on previous studies, the protein expression levels of several MMPs, such as MMP-1 and MMP-12, are associated with malignancy and metastasis in human ESCC (12,13). However, to date, there are no studies which have systematically assessed the impact of all MMP family members on the prognosis of ESCC to the best of our knowledge.

Bioinformatics analysis has been used to identify potential targets for cancer prevention and treatment (14,15), establishing a theoretical framework to efficiently guide experimental verification and research. Both Gene Expression Omnibus (GEO; ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA; cancer.gov/tcga.) databases contain clinical data on cancer patients and their gene expression profiles. In the present study, based on expression information from these two databases, whether MMPs could be used to predict the diagnosis and prognosis of patients with ESCC was assessed. The potential role of MMPs in the growth and migration of ESCC cells were assessed in vitro. The results showed that MMPs may serve as potential diagnostic markers of ESCC, and that inhibition of MMPs may be a potential therapeutic approach for treatment of ESCC.

Materials and methods

Patient data

Microarray data from GSE53625 (16), GSE23400 (17) and GSE38129 (18) datasets were obtained from the GEO database, whereas the RNA seq data from an ESCC cohort was downloaded from TCGA. Clinicopathological characteristics, including age, Tumor-Node-Metastasis (TNM) stage (19), tumor grade and sex, as well as survival data were included in the GSE53625 dataset (179 cases) and in TCGA (95 cases). Using SPSS version 23.0 (IBM, Corp.), the association between overall survival (OS) and clinicopathological characteristics of patients were investigated using univariate Cox regression analysis. In addition, characteristics with P<0.3 in the univariate analysis were further screened and used for multivariate Cox regression analysis.

Identification of differentially expressed genes (DEGs)

To identify important gene signatures within the MMP family of proteins, differential gene expression analysis was performed on the 4 datasets. The MMPs with a P value <0.05 were considered DEGs. Using R version 64 3.6.1 (20), an Empirical Bayes statistical test from the ‘limma’ package was used to analyze the GEO data, whereas an exact test in ‘edgeR’ package (which is analogous to the Fisher's exact test) was used to analyze TCGA data (21,22). The MMP expression levels in all 4 datasets were shown through construction of heat maps. Overlapping DEGs were shown in Venn diagrams, and their expression is presented as box diagrams.

Prognostic value analysis of the MMP family

Using the ‘survival’ package in R, the most suitable combination of DEGs, the one with the best predictive ability, was selected out through a stepwise multivariate Cox hazard regression analysis (23). Pearson correlation coefficients between all the MMP family members were also calculated. Subsequently, the screened MMPs were used to establish a prognosis scoring system. According to this analysis, the adjusted hazard ratio (HR) was calculated. The risk scores for patients with ESCC were calculated using the following formula:

where n indicates the gene number being used in the model, Coef is the coefficient of each gene, and Exp is the gene expression level. Based on the MMP signature, a nomogram which could predict the survival rate of patients was constructed (24). Also, the concordance index (C index), which can assess the accuracy of the prediction model, was calculated. Moreover, the receiver operating characteristic (ROC) curve was plotted and the area under the curve (AUC) was calculated to measure both the sensitivity and specificity of the prediction model which applied when the AUC value was >0.6 (25). Based on the median value of risk scores, the samples were stratified into low-risk and high-risk groups, and the prognostic difference between these two groups was investigated using Kaplan-Meier (K-M) survival curves. In addition, the predictive value of a survival prediction model based on the TNM stage using the ROC and the K-M curves was determined.

Nomogram integrating the signature of MMPs and TNM stage

Based on the coefficients from the multivariate Cox regression analysis, a nomogram integrating the MMP signature and TNM stage was constructed using the ‘rms’ package (26). The accuracy of the nomogram was evaluated based both on the C index and AUC values. Additionally, to graphically assess the performance of this nomogram, calibration curves were plotted. According to the median value of the risk scores calculated using Cox regression analysis, patients were separated into two groups; the distribution of the risk score and survival status of the patients were visualized to evaluate the prognostic difference between the two groups.

Gene set enrichment analysis (GSEA)

Using the expression profiles of tumor tissues, GSEA was used to identify significantly enriched pathways between the low-risk and high-risk groups defined by the MMP signature model (27). Oncogenic signatures gene sets (c6), Hallmark gene sets (h) and Kyoto Encyclopedia of Genes and Genomes gene sets (c2) in which the pathways are associated with the cancer process were used as references. Additional details can be found from the Molecular Signatures Database (28). Gene sets with P<0.05 were defined as indicators of significant differences. When gene sets had a normalized enrichment score (NES) >0, the pathway represented by this gene set was considered upregulated in the high-risk group, otherwise, it was considered downregulated. Subsequently, weighted gene co-expression network analysis (WGCNA) was performed using the ‘WGCNA’ package in R to identify co-expressed MMP genes within the model (29). The networks of co-expressed genes were drawn using Cytoscape (version 3.7.1) (30).

Cell lines and cell culture

KYSE30 and KYSE450 cell lines were purchased from The Cell Bank of Type Culture Collection of the Chinese Academy of Sciences. Cells were cultured in RPMI-1640 medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.), 100 U/ml penicillin (Gibco; Thermo Fisher Scientific, Inc.), and 100 µg/ml streptomycin (Gibco; Thermo Fisher Scientific, Inc.) and incubated at 37°C in a humidified incubator with 5% CO2.

Cell proliferation and colony formation assays

Cell proliferation was measured using a Cell Counting Kit-8 (CCK-8) assay (Beyotime Institute of Biotechnology). Briefly, ESCC cells were seeded in a 96-well plate with 2×103 cells/well. BB-94 (Selleck Chemicals, cat. no. S7155) was added to the wells at final concentrations of 0, 10, 20 or 40 µM. After 24, 48 or 72 h of treatment, 10 µl CCK-8 reagent was added to each well and incubated at 37°C for 2 h. The absorbance was measured at 450 nm using a microplate reader (Bio-Rad Laboratories Inc.).

To perform colony formation assays, 500 cells were seeded in a 6-well plate. After a 24 h incubation at 37°C, BB-94 was added to the wells to a final concentration of 0, 10, 20 or 40 µM. Cells were incubated for 7 days and the medium was removed and plates were washed with PBS. Cells were then fixed with methanol for 15 min at 25°C and stained with 0.1% crystal violet for 5 min at room temperature. Images of cell colonies were captured using a scanner (Canon).

Wound healing assay

ESCC cells were seeded in 6-well plates. Using the tip of a 20 µl micropipette, a scratch was made in the middle of the well, and cells were washed 3 times with PBS. ESCC cells were incubated in 2 ml serum-free medium containing different concentrations of BB-94 (0, 10, 20 or 40 µM). After 24 h, migration was observed and recorded under an inverted light microscope (magnification, ×40; Zeiss GmbH). The distances between the edges of the scratches were measured using ImageJ (FIJI distribution, version 1.52n, National Institutes of Health).

Western blotting

Proteins were extracted from ESCC cells using RIPA lysis buffer containing both a protease inhibitor cocktail and a phosphatase inhibitor cocktail (Beyotime Institute of Biotechnology). Protein concentrations were determined using a bicinchoninic acid assay kit (Beyotime Institute of Biotechnology). Proteins were resolved using SDS-PAGE with a 6–12% separation gel and 5% concentration gel using the Laemmli discontinuous buffer system and transferred to PVDF membranes (EMD Millipore). After blocking the membranes in 5% skimmed milk for 1 h at room temperature, they were incubated with primary antibodies diluted in PBS with 0.1% Tween 20 (PBST) overnight at 4°C. Primary antibodies used were anti-GAPDH (1:3,000; cat. no. 5174S) and an epithelial-mesenchymal transition (EMT) Antibody Sampler kit (cat. no. 9782T), which included E-Cadherin (1:2,000), Vimentin (1:3,000), β-Catenin (1:3,000), snail (1:3,000) and slug (1:3,000) antibodies, all of which were purchased from Cell Signaling Technology, Inc. After washing the membranes 3 times with PBST, they were incubated with secondary horseradish peroxidase (HRP)-conjugated antibodies (1:2,000; cat. no. 7074; Cell Signaling Technology, Inc.) for 1 h at room temperature. Signals were visualized using a Chemiluminescent HRP Substrate (EMD Millipore) and visualized with a UVP GelStudio PLUS Touch Imaging system (Analytik Jena).

Statistical analysis

R version 3.61, SPSS version 23.0 (IBM Corp) and GraphPad Prism version 8 (GraphPad Software, Inc.) were used to analyze obtained data. Two-tailed P values <0.05 were considered to indicate a statistically significant difference. The difference in MMP expression levels between normal and tumor tissues was compared using a paired Student's t-test. Univariate and multivariate Cox regression analyses were used to screen the independent prognostic variables of OS; and the variables considered significant (P<0.05) were used to establish the risk score formula and to construct the nomograms. A stepwise multivariate Cox hazard regression analysis was used to screen MMPs to construct the best prediction model. Samples were separated into two groups (low- and high-risk groups) based on the median value of the risk score. K-M survival curves were plotted, and survival was compared using a log rank test. The AUC was used to as a measure of the accuracy of the predictive models at the indicated times. The performance of the nomogram was validated by assessing the calibration curves as well. Experimental data are presented as the mean ± standard deviation of at least 3 independent experiments and results were analyzed using an ANOVA with a Bonferroni post hoc test.

Results

Cox regression analysis of ESCC patient clinical data

The GSE53625 dataset and TCGA database which included the complete clinical data of patients (Tables SI and SII, respectively) were used for Cox regression analysis. Univariate Cox regression analysis of the GSE53625 dataset showed that the TNM stage, age and N stage were significantly associated with OS (P<0.001, P=0.021 and P=0.030, respectively; Fig. 1A). Characteristics with P<0.3 in the univariate analysis were further screened and used for multivariate analysis. The TNM stage was an independent prognostic factor (P=0.001; Fig. 1B). Furthermore, univariate Cox regression analysis of the data from TCGA showed that sex, TNM stage and N stage were significantly correlated with OS (P=0.020, P=0.015 and P=0.012, respectively; Fig. 1C). Multivariate Cox regression analysis indicated that both sex and the N stage were independent prognostic factors (P=0.047 and P=0.012, respectively; Fig. 1D).

Figure 1.

Cox regression analysis of clinical data of patients with ESCC. Forest plots of univariate and multivariate Cox regression analysis of clinicopathological characteristics affecting OS in patients with ESCC in (A and B) the GSE53625 dataset and (C and D) TCGA. *P<0.05, **P<0.01, ***P<0.001. TNM stage, Tumor-Node-Metastasis stage; HR, hazard ratio; 95% CI, 95% confidence interval; T stage, stage of tumor invasion; N stage, stage of regional lymph node invasion; ESCC, esophageal squamous cell carcinoma; OS, overall survival; TCGA, The Cancer Genome Atlas.

Identification of differentially expressed MMPs in the 4 datasets

DEGs of the MMP family were identified; detailed results including the logarithm of fold change (log FC) and P values are shown in Tables SIII–SVI. Gene expression profiles of the MMP family members are displayed in heat maps and 20, 13, 12 and 17 DEGs were identified in GSE53625, GSE23400, GSE38129 and TCGA datasets, respectively (Fig. 2A-D). A Venn diagram showed 8 overlapping MMPs (MMP-1, −3, −9, −10, −11, −12, −13 and −14) were dysregulated in ESCC (Fig. 2E). In addition, box diagrams showing the expression levels of these overlapping DEGs from the 4 datasets showed that all these MMPs were upregulated in tumor tissues (Fig. 2F-I). These results suggest that these MMPs may serve as potential diagnostic markers for ESCC.

Figure 2.

Differential MMP expression analysis between tumor and normal tissues in ESCC. Heat maps, sorted in descending order of log FC, including every MMP in (A) GSE53625, (B) GSE23400 and (C) GSE38129 datasets, and in (D) TCGA database. Red and blue colors represent high and low expression, respectively. (E) Venn diagram showing overlapping dysregulated MMPs based on the 4 datasets mentioned above. (F-I) Box diagrams showing the expression levels of the 8 overlapping dysregulated MMPs in tumor tissues compared with normal tissues in all 4 datasets. *P<0.05, **P<0.01, ***P<0.001. MMP, matrix metalloproteinase; ESCC, esophageal squamous cell carcinoma; FC, fold change; TCGA, The Cancer Genome Atlas.

Prognostic value of the MMP signature and TNM stage models for ESCC

To investigate the prognostic values of MMPs, DEGs were screened to construct a prediction model based on GSE53625 and TCGA, both of which contained clinical data. Using the GSE53625 dataset, a heat map showing the log FCs and HRs of DEGs was constructed; from these, 6 differentially expressed MMPs (MMP-1, −8, −9, −13, −25 and −28) were screened out to construct a prediction model using a stepwise multivariate Cox hazard regression analysis (Fig. 3A). Subsequently, the adjusted HRs of the 6 selected MMPs were shown in a forest plot, which suggested that MMP-1 (HR=0.76; CI, 0.65–0.88; P<0.001), MMP-8 (HR=1.18; CI, 1.03–1.36; P=0.019) and MMP-25 (HR=0.72; CI, 0.53–0.98; P=0.036) were independent prognostic factors (Fig. 3B). According to the coefficients and gene expression levels (Exp) of these MMPs, the risk score of ESCC patients was calculated based on the following formula: Risk score=Exp MMP-13 ×0.82-Exp MMP-1 ×2.72-Exp MMP-25 ×2.23+Exp MMP-28 ×1.04+Exp MMP-8 ×0.71+Exp MMP-9 ×1.80. In addition, an MMP signature nomogram was constructed to predict the survival rates of ESCC patients (C index=0.617, Fig. 3C). To evaluate the accuracy of this model, ROC curves were drawn, and their relative AUC values were calculated. As shown in Fig. 3D and E, the AUC value of the MMP signature model reached 0.671, which was higher than that of the prediction model based on the TNM stage (AUC=0.637) (31). The latter is commonly used to predict the prognosis of cancer patients. Finally, based on the median risk score, ESCC patients were divided into low-risk and high-risk groups; a survival curve indicated that patients with lower risk scores had improved survival (P=0.009; Fig. 3F). Patients at TNM stage I/II also showed significantly higher survival rates compared with patients at TNM stage III (P=0.000; Fig. 3G). Survival and ROC curves of each selected MMP are shown in Fig. S1A-B, and the survival curves were generally consistent with the results from the website OSescc (bioinfo.henu.edu.cn/DBList.jsp), a tool to assess OS and relapse free survival based on the expression of given genes or probes (32). To analyze the correlation among each MMP gene, Pearson correlation coefficients of all the MMP family members were also calculated (Fig. S1C). Notably, similar results, including both ROC and survival curves of the prediction model, were obtained using TCGA (81 cases) which had a smaller sample size than the GSE53625 dataset (179 cases; Fig. S1D-G). Thus, MMPs were accurate prognostic predictive factors of ESCC.

Figure 3.

Prediction models based on MMP expression or TNM stage for ESCC. (A) Heat maps showing the log FCs and HRs of DEGs in the GSE53625 dataset. Based on six of the DEGs, a stepwise multivariate Cox hazard regression analysis was performed to screen the best MMP signature. (B) Forest plot of the adjusted HRs, 95% CIs of HR, and P values of the MMPs included in the model. (C) Nomogram to visualize the MMP signature prediction model. (D and E) ROC curves of both the MMP model and the TNM stage model. (F and G) Kaplan-Meier curves of the risk scores calculated using the prediction models. *P<0.05, ***P<0.001. MMP, matrix metalloproteinase; TNM stage, Tumor-Node-Metastasis stage; ESCC, esophageal squamous cell carcinoma; FC, fold change; HR, hazard ratio; DEG, differentially expressed gene; 95% CI, 95% confidence interval; ROC, receiver operating characteristic.

Prognostic value of a model combining the MMP signature model and the TNM stage model

To develop a more accurate prediction model, a nomogram was constructed integrating both the MMP signature and TNM stage (C index=0.652; Fig. 4A). The ROC curves indicated that the MMP-TNM-integrated nomogram had improved predictive ability compared with each model alone (AUC=0.721 for 3-year survival, AUC=0.717 for 5-year; Fig. 4B). The calibration curves indicated that the predicted outcome was broadly consistent with the actual outcome, suggesting that the new model accurately predicted the results (Fig. 4C). In addition, the risk scores of patients were calculated using the combined model, resulting in a reclassification of patients as belonging to either a low-risk or a high-risk group (Fig. 4D). Accordingly, in the distribution chart of survival status, the number of surviving patients from the low-risk group was 51, while the number of deaths was 39. By contrast, there were 22 patients alive and 67 dead in the high-risk group (Fig. 4E). The distribution map of survival status showed that the patients from the high-risk group had higher mortality rates. Therefore, the MMP-TNM-integrated nomogram may be a more effective tool for clinicians to predict prognosis of ESCC patients.

Figure 4.

Prediction model combining the MMP signature and TNM stage information. (A) Nomogram integrating MMP and TNM stage for prediction of the 3-year and 5-year survival rates. (B) ROC curves showing the predictive accuracy of the nomograms for 3-year and 5-year survival rates. (C) Calibration curves showing the calibration of the model either for a 3-year or a 5-year survival rate prediction in terms of the agreement between predicted and observed outcomes. (D) Risk score distribution of patients in ascending order, classified into low-risk (green) and high-risk (red) groups. (E) Survival time and status of patients plotted in order of increasing risk scores; red and green dots represent dead and alive patients, respectively. MMP, matrix metalloproteinase; TNM stage, Tumor-Node-Metastasis stage; ROC, receiver operating characteristic; AUC, area under the curve; OS, overall survival.

Biological function of MMP family members in ESCC

To determine how MMPs affect ESCC cells, GSEA was performed based on the risk score calculated by the MMP signature model. Tables SVII–IX show all the results which were significant (P<0.05), whereas Fig. 5 shows those the signaling pathways that were tightly associated with cancer. Based on hallmark gene sets, signaling pathways including EMT, TGF-β, P53 and Notch signaling were significantly enriched in the high-risk score group (Fig. 5A). Cancer-associated pathways from the KEGG gene sets included focal adhesion, ECM receptor interaction, Notch signaling pathway and mTOR signaling (Fig. 5B). In addition, results based on oncogenic signature gene sets indicated that the high-risk score group was significantly associated with TGF-β, Notch, KRAS and P53 pathways (Fig. 5C). Moreover, a co-expression gene network of the 6 MMPs included in the MMP signature was constructed (Fig. 5D). These 6 MMPs are marked as large red nodes, whereas smaller red nodes are used to represent the other MMP family members co-expressed with these 6 MMPs. Blue nodes represent the other co-expressed genes. The gene network clearly showed that vimentin and SNAI1 were associated with both MMP-13 and MMP-9. SNAI1 was also correlated with MMP-1, whereas TGFB1 was associated with MMP-13. In addition, COL12A1, COL1A1, COL5A2, COL6A1 and COL6A3, which belong to collagen family of proteins, were associated with MMP family members. These results suggest that MMPs may participate in ESCC cell growth and migration.

Figure 5.

Pathways and co-expressed genes associated with the MMP signature model. (A) Enriched pathways based on hallmark gene sets. (B) Enriched pathways based on KEGG gene sets. (C) Enriched pathways based on oncogenic signatures gene sets. (D) Gene co-expression network of MMPs in the prediction model. MMP, matrix metalloproteinase; KEGG, Kyoto Encyclopedia of Genes and Genomes; EMT-epithelial-mesenchymal transition; ECM, extracellular matrix.

Effect of MMP inhibition on ESCC cell growth

Based on the results of GSEA, the effect of inhibition of MMPs using BB-94, a pan MMP inhibitor on ESCC cell proliferation was assessed. BB-94 significantly reduced viability of ESCC cells after 72 h. (P<0.05; Fig. 6A and B). Colony formation assays also revealed that, compared with the control group, the number of colonies formed in the BB-94-treated groups decreased significantly (Fig. 6C-F). These results suggest that MMP inhibition reduced growth of ESCC cells.

Figure 6.

ESCC cell growth following inhibition of MMP activity using BB-94. Viability of (A) KYSE30 and (B) KYSE450 cells treated with different concentrations of BB-94 measured after 0, 24, 48 and 72 h. Representative images of colony formation assays in (C) KYSE30 and (D) KYSE450 cells treated with different concentrations of BB-94. Scale bar, 1 cm. Quantitative analysis of colony formation in (E) KYSE30 and (F) KYSE450 cells. **P<0.01, ***P<0.001 vs. Con. ESCC, esophageal squamous cell carcinoma; OD, optical density; MMP, matrix metalloproteinase.

MMPs affect migration and expression of EMT markers in ESCC cells

GSEA and gene co-expression analysis suggested that there was a tight association between MMPs and cell migration. Therefore, the effect of BB-94 on ESCC cell migration was examined. Wound healing rate was used to evaluate the migration of ESCC cells treated with different concentrations of BB-94. When compared with the control group, the wound healing rates in BB-94 treated groups were significantly reduced, both in KYSE30 and KYSE450 cell lines (Fig. 7A-D). Furthermore, inhibition of MMPs by BB-94 altered the protein expression levels of several EMT markers in both KYSE30 and KYSE450 cells. E-Cadherin expression was increased and expression of Vimentin, β-Catenin, Snail and Slug were decreased (Fig. 7E and F). Finally, the Pearson correlation coefficients between MMPs selected in the prediction model and EMT markers were calculated. Gene expression of MMP-9 and MMP-13 were positively associated with Vimentin. In addition, MMP-1, −9 and −13 were positively associated with SNAI1 expression (Fig. 7G). These results suggest that the MMP family members are involved in the regulation of ESCC migration.

Figure 7.

ESCC cell migration and association between expression of EMT markers and MMP expression. Representative images of wound healing assays in (A) KYSE30 and (B) KYSE450 cells treated with different concentrations of BB-94. Quantitative analysis of wound healing in the (C) KYSE30 and (D) KYSE450 cells. (E and F) Protein expression levels of the EMT markers, E-Cadherin, Vimentin, β-Catenin, Snail and Slug were detected by western blot. (G) Pearson correlation coefficient plots between expression of MMPs and EMT markers based on the GSE53625 dataset. ***P<0.001 vs. Con. ESCC, esophageal squamous cell carcinoma; EMT, epithelial-mesenchymal transition; MMP, matrix metalloproteinase; Con, control; Cor, correlation coefficient; VIM, vimentin.

Discussion

Poor prognosis for ESCC patients is largely due to delayed diagnosis. The 5-year survival rate for patients with ESCC is relatively higher when the cancer is diagnosed at an earlier stage (33), highlighting the importance of novel biomarkers for early detection of ESCC. MMP family members degrade ECM structural components, and this underlies tumor progression (34). Duffy et al (35) demonstrated that MMP family members participate in tumor initiation and progression of breast cancer. The MMP family has been extensively studied in different types of cancer, highlighting the importance of the involvement of MMPs in development and progression of various types of cancer, but previous studies were not systematic, instead focusing on one or a few MMPs (36,37). Furthermore, the prognostic role of specific MMPs remains controversial. For example, MMP-9, which is one of the best studied MMPs, was identified as a significant prognostic factor in certain studies (38,39), whereas in other studies no significant effects were observed (40,41). Unlike the previous studies that only focused on one or several MMP proteins detected by immunohistochemistry or ELISA, in the present study, the transcriptomics of all 23 MMPs were systematically investigated, with a focus on the prognostic value of the combination of MMPs in ESCC using bioinformatics analysis. The results of the bioinformatics analysis were subsequently confirmed in vitro.

Upregulated expression of MMP-1, −9, −13 has been previously observed in ESCC tumor tissues compared with normal tissues using immunohistochemistry (39). Han et al (13) also showed that MMP-12 expression was also upregulated in ESCC tumor tissues. Consistent with these studies, 8 MMPs, MMP-1, −3, −9, −10, −11, −12, −13 and −14) were differentially expressed in 4 different datasets; all of which were upregulated in tumor tissues and shown to be clinically significant for potential use as diagnostic markers for ESCC. The majority of previous studies reported to date suggest that increased gene expression levels of certain MMPs are significantly associated with a poor prognosis of ESCC (36,38). The use of multivariate Cox regression analysis to construct a prediction model incorporating several related genes provides a more robust predictive effect than models based on a single gene (42). In the present study, all the MMP genes deemed to be dysregulated were used to establish a survival prediction model. Using the GSE53625 dataset, after screening using a stepwise multivariate Cox hazard regression analysis, the best prognostic MMP signature was established. Of note, none of the genes in this signature alone was considered a significant prognostic factor (data not shown), but the combination of these genes accurately predicted the prognosis of ESCC patients. Compared with the survival prediction model based on the TNM stage that is commonly used (31), the model established in the present study predicted the survival rate of patients with ESCC with improved accuracy, highlighting its suitability for clinical use. The second model, comprised of screened MMPs from TCGA, also showed an improved predictive ability compared with the TNM stage model; the corresponding survival curves showed a similar tendency, that is, the high-risk group had a worse prognosis compared with low-risk group. Although the GSE53625 dataset was a considerably larger dataset than TCGA (179 samples vs. 81), similar results were obtained from both datasets, suggesting that the MMP-based prediction model was more powerful for predicting survival. In addition, a nomogram was constructed integrating both the MMP signature and TNM stage, and this exhibited improved accuracy compared with the MMP signature-based model and may thus be a more suitable tool for clinicians to predict the prognosis of patients with ESCC. This type of prediction model is in line with current trends toward personalized medicine and being more commonly used in other research areas (43,44). Therefore, the MMP-TNM nomogram may be a vital tool for predicting survival of patients with ESCC.

Several signaling pathways associated with cancer progression, such as EMT, Notch, TGF-β, mTOR and P53, are associated with the expression of different members of the MMP family of proteins. MMPs could stimulate processes associated with EMT (45). For instance, the expression of MMP-9 was associated with EMT genes in breast cancer samples (46). Pang et al (47) also showed that MMP-14 increases ESCC invasion and metastasis by reducing E-Cadherin expression and subsequently inducing EMT. Overexpression of MMP-13 was observed in ESCC clinical tissues, and its upregulated expression increased cancer cell aggressiveness (48). Knockdown of MMP-9 attenuated EMT induced by TGF-β1, and inhibited invasiveness and migration in ESCC (49). Notably, numerous studies have demonstrated the association between MMPs and Notch signaling in several types of cancer (50,51), the latter of which serves an important role in development and determination of cell fate (52). For example, a previous study identified Notch1 as an MMP-14 substrate in melanoma (53). Moreover, activation of the PI3K/PTEN/AKT/mTOR pathway upregulated MMP-9 expression in hepatocellular carcinoma (54). On the contrary, P53, as a metastasis suppressor, downregulated MMP-1 and MMP-9 expression (55).

Regarding co-expression, several MMP genes are co-expressed with PDGFRB, CREB3L1, COL1A1 and other collagen family members, such as COL12A1, COL5A2, COL6A1 and COL6A3, suggesting that MMPs participate in collagen-mediated metabolic processes in ESCC. MMP-13 expression is downregulated by CREB3L1, a metastasis suppressor which regulates the expression of a number of genes involved in angiogenesis (56). Co-expression analysis showed that certain collagen family members, which are involved in cell proliferation and migration, are co-expressed with specific MMP genes (58). In addition, co-expressed genes, including QKI, KDR, PDGFRB and COL1A1, are associated with vasculature development (GO:0001944). In a previous study, it was demonstrated that both MMP-13 and PDGFRB were upregulated in papillary thyroid carcinoma (59). It is widely known that MMPs are associated with EMT, in the gene co-expression network analysis, it was shown that MMPs were correlated with EMT markers such as VIM and SNAI1. In addition, MMP-13 was associated with TGFB1, which is important for tumor viability, migration and metastasis in multiple types of cancer (60). Taken together, these results highlight how MMP family members may influence the progression of ESCC and may be used to identify novels targets for the development of specific therapeutic strategies for ESCC treatment.

Cell proliferation and migration ability assays were performed to confirm the role of MMP family members in ESCC cell proliferation and migration using a pan-MMP inhibitor. It has been demonstrated that the pan-MMP inhibitor BB-94 inhibits growth and metastasis of human colon tumors in a patient-like orthotopic model in nude mice (61). However, to the best of our knowledge, there are no studies assessing the effect of BB-94 in ESCC. In the present study, the cytotoxic and phenotypic changes observed following treatment with BB-94 were consistent with previous studies (62,63). Inhibition of MMPs by BB-94 significantly reduced cell proliferation and suppressed EMT. Thus, BB-94 may be a potential therapeutic agent for treatment of ESCC. Further studies are required to determine the detailed mechanism and in vivo anti-ESCC activity of BB-94.

The present study has several limitations. There were 81 tumor samples, whereas only 11 normal samples in TCGA. In TCGA, the number of normal samples is frequently lower than the number of tumor samples in several types of cancer, such as bladder urothelial carcinoma, breast invasive carcinoma and colon adenocarcinoma (64). Thus, the GEO database GSE53625 dataset (including 179 normal samples and 179 tumor samples) was used to reduce the potential bias introduced by the large difference in the number of samples. The design of the inhibition experiment was not completely consistent with our prediction model. The optimal ESCC therapeutic regimen would promote the activation of MMP-1 and MMP-25, and inhibit MMP-8, −9, −13 and −28 concurrently. However, using a cocktail of agents able to exert these effects is difficult. BB-94 is one of the most widely used broad-spectrum inhibitors (63,65). However, using BB-94 may have suppressed several other MMPs upregulated in tumor tissues with different risk contributions. Thus, future studies should focus on developing specific inhibitors targeting certain MMPs with high risk scores. Another limitation was that gene expression profiles were used and the results showed that no single MMP had a significant predictive ability by itself, but instead the combination of specific MMPs exhibited good predictive ability. Thus, the expression of MMP protein expression levels or activity were not measured. Further testing regarding MMP protein expression or activity in ESCC cell lines treated with BB-94 and relevant clinical specimen are required to confirm their potential involvement.

In summary, members of the MMP family may be used as diagnostic and prognostic markers for ESCC. The MMP signature model was more accurate for predicting the survival of patients with ESCC compared with the TNM stage-based model. When integrating the MMP signature model and TNM stage to predict the survival rate, accuracy was further improved, suggesting that MMPs have considerable predictive value. Moreover, it was shown that MMP family members may influence the prognosis of ESCC through impacting signaling pathways involved in cancer, such as EMT, TGF-β, Notch, mTOR and P53. Finally, wound healing migration assays and western blotting showed that inhibition of MMPs using BB-94 reduced migration of ESCC cells by suppressing EMT. Thus, MMP family members may constitute potential therapeutic targets for prevention and treatment of ESCC.