High Matrix Metalloproteinase 28 Expression is Associated with Poor Prognosis in Pancreatic Adenocarcinoma

Abstract

Purpose

Pancreatic adenocarcinoma (PAAD) is a devastating disease with high mortality and morbidity. Matrix metalloproteinase 28 (MMP28) has been associated with carcinogenesis of many human cancers. However, little is known about the potential prognostic value and underlying regulatory mechanisms of MMP28 in PAAD.

Methods

The relationship between MMP28 expression level and various clinicopathological parameters was analyzed in TCGA-PAAD cohorts. MMP28-correlated genes in the TCGA-PAAD cohort were identified and enrichment analysis according to the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes was conducted using LinkedOmics. Protein–protein interaction and transcription factors-miRNA co-regulatory networks were constructed with the use of NetworkAnalyst. Then, the distribution of immune cells related to MMP28 expression in blood was analyzed using the Human Protein Atlas, and the tumor microenvironment of PAAD was analyzed by the TIMER 2.0 database. To investigate the biological function of MMP28 in PAAD, siRNA was constructed to knock down the MMP28 gene in vitro.

Results

High MMP28 expression is associated with poor overall survival and disease-free survival in PAAD patients. The expression of MMP28 in PAAD is most significantly correlated with KRT19, IL1RN, and ANXA2 genes. Network analysis revealed that MIR-181 family, TAFs, and CDC6 are potential regulators of MMP28. Furthermore, naive CD4⁺ T cell, naive CD8⁺ T cell, and mucosal-associated invariant T cell enrichment in blood were correlated with MMP28 expression. Furthermore, high MMP28 expression was correlated with a decrease in B cell, naive CD4⁺ T cell, naive CD8⁺ T cell, and endothelial cell presence in the tumor microenvironment in PAAD. Finally, genetic knockdown of MMP28 could restrain the proliferation, migration, and invasion of PAAD cells.

Conclusion

Our findings indicate that high MMP28 expression in PAAD is associated with cancer progression, invasion, and metastasis. Hence, MMP28 might serve as an independent prognostic biomarker and a prospective therapeutic target for PAAD.

Keywords:

• matrix metalloproteinase 28
• pancreatic adenocarcinoma
• metastasis
• prognosis

Introduction

Pancreatic adenocarcinoma (PAAD) is one of the most lethal malignancies globally, characterized by poor prognosis and high morbidity.¹ Compared with other types of cancer, PAAD has the lowest 5-year relative survival rate (ie, 9%). In line, PAAD is the fourth leading cause of cancer death in the United States of America.² Although many research efforts have been made to improve treatment strategies for PAAD, the prognosis of this cancer remains poor, as it is commonly diagnosed at a late stage. More than 90% of patients with PAAD have metastasis at diagnosis.^3,4 Therefore, there is an urgent need to identify effective prognostic factors for this cancer.

Matrix metalloproteinases (MMPs) are a family of zinc-dependent endopeptidases, which plays an important role in the degradation of proteins in the extracellular matrix (ECM). ECM remodeling is regulated by MMPs, thus MMPs can modify tumor microenvironment (TME) and contribute to the formation of a premetastatic niche and stimulate angiogenesis.⁵ On top of that, MMPs are also involved in various pathological processes related to carcinogenesis, including regulation of cancer progression and apoptosis, promotion of cancer invasion, metastasis, and epithelial-to-mesenchymal transition (EMT).^6,7 Several previously identified members of the MMP family such as MMP1, MMP2, MMP7, MMP11 and MMP14 have been found to be highly expressed in PAAD and their roles in PAAD have been studied.^8–13 However, the effect of MMP28 overexpression on the prognosis of PAAD and the possible mechanism remain unclear.

MMP28, also known as epilysin, belongs to the MMP-19 subfamily. MMP28 gene is located on the chromosome 17q11.2 and contains eight exons and seven introns. It is expressed in the adult pancreas and its expression is elevated in PAAD.¹⁴ Previous studies have suggested that MMP28 plays an important role at multiple stages of tumor growth. MMP28 is involved in fibrosis formation and in inflammatory cell recruitment in pulmonary emphysema.¹⁵ MMP28 is correlated with atrial fibrillation and indicates poor prognosis of a heart failure.¹⁶ In hepatocellular and glioblastoma carcinoma, upregulated MMP28 accelerates the migration and invasion and is linked to poor prognosis.^17,18 Inhibiting the transcription of MMP28 could suppress the invasion and metastasis of gastric cancer.¹⁹ In PAAD, MMP28 shows a cancer-promoting effect and is involved in tumor formation.^20,21 However, the exact prognostic relevance of MMP28 and an underlying mechanism remain poorly understood in the context of PAAD.

In this study, we investigated the association of MMP28 expression with clinical and genomic features in PAAD patient cohorts using public databases. Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed and functional networks related to MMP28 in PAAD were constructed. Distribution of immune cells associated with MMP28 in the blood was also analyzed to shed light on possible effects of MMP28 in the TME. To investigate the impact of MMP28 on PAAD cell biology, small interfering RNA (siRNA) was constructed to knockdown the MMP28 gene in vitro. Our results indicate that MMP28 might be a new prognostic biomarker in PAAD and support targeting it in therapy.

Materials and Methods

Oncomine Database Analysis

Oncomine (https://www.oncomine.org/) is a cancer microarray database and a web-based genetic data-mining platform. The mRNA expression of MMP28 in various types of cancer was examined in the Oncomine database. The data in the Pei’s Pancreas dataset were also used to compare the MMP28 mRNA expression levels between PAAD and pericarcinomatous tissue, using a Student’s t-test for P value analysis. The threshold was defined according to the following criteria: P value of 0.01, fold change of 1.5, and gene ranking of the top 10%.

GEPIA Database Analysis

The Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/) is a newly developed interactive web for evaluating the RNA sequencing expression data, which includes 9,736 tumors and 8,587 normal samples from TCGA and the Genotype Tissue Expression (GTEx) projects.²² GEPIA was used to explore the expression of MMP28 mRNA in PAAD survival analysis.

Human Protein Atlas Database Analysis

Human Protein Atlas (HPA) (https://www.proteinatlas.org/), a public protein database, is a useful diagnostic tool based on compound clinical information and pathological images of various types of cell lines, cancer tissues, and normal tissues.²³ In this study, the distribution of MMP28 RNA in immune cells present in the blood was investigated.

LinkedOmics Database Analysis

The LinkedOmics database (http://www.linkedomics.org/login.php) is a web-based platform, enabling the analysis of 32 TCGA cancer-associated multi-angle datasets.²⁴ Genes correlated with MMP28 in PAAD were statistically analyzed using Pearson’s correlation coefficient, and visualized by volcano plots and heat maps. Enrichment analysis of GO term, KEGG pathways, miRNA, kinases, and transcription factors was performed using LinkedOmics functional tools. False discovery rate (FDR) <0.05 and 1000 simulations were set as the thresholds.

cBioPortal Database Analysis

The cBio Cancer Genomics Portal (http://www.cbioportal.org/) contains cancer genomics data sets and is used to study genetic changes in various cancer types.²⁵ In our study, mutations, copy number alterations, and genes correlated with MMP28 in PAAD were analyzed with cBioPortal tools.

NetworkAnalyst Database Analysis

NerworkAnalyst 3.0 tool (https://www.networkanalyst.ca/) is a useful tool used for creating generic or tissue-specific protein-protein interaction (PPI) networks, gene co-expression networks, and transcription factor (TF)-miRNA coregulatory networks.²⁶ In our study, MMP28-correlated genes in PAAD, as listed by cBioPortal, were used to construct PPI and TF-miRNA coregulatory networks. PPI networks were executed based on the method of Search Tool for the Retrieval of Interacting Genes (STRING).

Cytoscape Software

Cytoscape (version 3.8.0) is an open-source software project, which provides basic tools useful in the analysis of biological networks, for instance to illustrate the relation of a given network with different gene expression profiles.²⁷ The plug-in Molecular Complex Detection (MCODE) (version 1.6.1) of Cytoscape is an APP for clustering a given network according to its topology to find intensively connected regions.²⁸ Following, PPI and TF-miRNA coregulatory networks were drawn using Cytoscape, and the most significant modules and related proteins or genes were identified using MCODE. The criteria for selection were as follows: degree cutoff =2, node score cutoff =0.2, K-core =2 and Max depth =100.

TIMER and TIMER 2.0 Database Analysis

TIMER (https://cistrome.shinyapps.io/timer/) is a comprehensive resource for immune infiltrate analysis in various cancer types based on the TCGA.²⁹ It utilizes a deconvolution method to estimate the abundance of tumor-infiltrating immune cells (TIICs).³⁰ In our study, we used the somatic copy number alteration (SCNA) module to compare the infiltration levels of TIICs in PAAD in relation to the somatic copy number alteration of MMP28. We also analyzed the correlation between MMP28 expression and the infiltration levels of various TIICs, including B cells, CD4⁺ T cells, CD8⁺ T cells, dendritic cells, cancer-associated fibroblasts, and endothelial cells by TIMER 2.0. The associations were assessed by partial Spearman correlation adjusted for tumor purity.

Cell Culture and PAAD Cell Lines

Cell lines PANC-1 and SW1990 were purchased from Shanghai Enzyme Research Biotechnology Co., Ltd. (Shanghai, China) and were donated to us by the Science Park of West China Hospital, Sichuan University (Chengdu, China).²¹ The use of these cell lines were approved by the Ethics Committees of West China Hospital of Stomatology, Sichuan University. All cell lines were maintained in Dulbecco’s Modified Eagle Medium (DMEM) and supplemented with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100μg/mL streptomycin (GIBCO, Grand Island, NY, USA) at 37 °C with 5% CO₂.

SiRNA and Transfections

Small interfering RNAs (siRNAs) targeting hMMP28 were purchased from Sangon Biotech (Shanghai, China); their sequences are shown in Table S1.

Transient transfections into PANC-1 and SW1990 cell lines were performed with 10 nM Lipofectamine RNAiMAX (Thermo Fisher Scientific, Waltham, MA, USA).

In the blank group, only DMEM medium was used, without siRNA. The negative control vector contained an insert with no significant homology to mammalian genes present in humans, and it was used as a negative control for MMP28 siRNA vectors. In negative controls, Lipofectamine Max (Invitrogen) and negative control siRNA were used.

Quantitative Real-Time PCR Analysis

Total RNA of PANC-1 and SW1990 cells were isolated using TRIzol reagent (Invitrogen, USA). cDNA samples were prepared using the RT reagent Kit gDNA Eraser (TaKaRa). cDNA expression levels were detected by SYBR-Green (TaKaRa) with qRT-PCR analysis. GAPDH was used as an internal reference gene. Primers with the following sequences were used: MMP28, Forward (F):5ʹ-TCCCACCTCCACTCGATTCAG-3ʹ, Reverse (R):5ʹ-GCCGCATAACTGTTGGTATCT-3ʹ; glyceraldehyde 3-phosphate dehydrogenase (GAPDH), Forward (F): 5ʹ-ACAACTTTGGTATCGTGGAAGG-3ʹ, Reverse (R):5ʹ- GCCATCACGCCACAGTTTC-3ʹ.

Western Blot Analysis

Total protein was extracted with RIPA lysis buffer (Solarbio, Beijing, China) and PMSF (Solarbio, Beijing, China). 20μg of cell protein lysates was loaded on 10% tris-polyacrylamide gels (SDS–PAGE, Bio-Rad, USA) and transferred into polyvinylidene difluoride (PVDF) membranes (Solarbio, China). The primary antibodies were used at the following dilutions: hMMP28 (Abcam Corporation, Cambridge, MA, USA), 1:2000; GAPDH (Abcam Corporation, Cambridge, MA, USA), 1:3000 (used as an internal loading control). Secondary antibodies were purchased from ZSGB-BIO (Beijing, China), and used in the 1:5000 concentration.

Cell Counting Kit-8 (CCK8) Assay

The cells were seeded at a density 8.0×10³ cells/well in the 96 well plate. Then cells were incubated at 37°C in 5% CO₂ for 0 h, 12 h, 24 h, or 48 h. To each well, 10 µL of cell counting Kit-8 (Dojindo, Japan) solution was added and incubated for 1.5 h at 37°C. The cell viability was detected by a microplate reader (Thermo Scientific, Varioskan Lux). The absorbance of samples was measured at 450 nm.

Migration and Invasion Assays

PANC-1 and SW1990 cells were transiently transfected with the siRNA for 48 h prior to migration and invasion assays. A 24 well plate with 8-μm pore filters (Corning, New York, NY, USA) was used in these assays. The transwell chamber was coated with Matrigel (BD Bioscience, San Diego, CA, USA) for 1 h at 37 °C. 8×10⁴ cells were seeded in the upper chamber in 200 μL serum-free basic DMEM medium, and the lower chamber contained 500 μL of complete culture medium. After 24 h, cells that had moved across the transwell membrane were stained with 0.5% crystal violet for 20 min and were counted under a microscope in five randomly selected fields.

Statistical Analysis

Various R packages were applied for statistical analysis and visualization, such as ggsci, ggplot2 etc, PAAD data of TCGA database. Detailed statistical methods related to other public databases are described in the corresponding paragraphs.

In in vitro studies, the Prism 7.0 software (GraphPad, La Jolla, CA) was used for all data analyses. Differences between groups were determined with a Student’s t-test, one-way ANOVA and two-way ANOVA. All quantitative data were displayed as means ±SD. P <0.05 was considered as statistically significant.

Results

Elevated Expression of MMP28 in PAAD

Initially, we analyzed the mRNA expression levels of the MMP family in the pericarcinomatous tissue and PAAD in the TCGA database. We found that MMP28 was highly expressed in PAAD. In addition, some other members of the MMP family, such as MMP1, MMP2, MMP7, MMP11 and MMP14, were also highly expressed in PAAD, but their roles in PAAD have been studied. (Figure 1A).^8–13 Next, we explored the mRNA expression level of MMP28 across multiple cancer cohorts from TCGA using the Oncomine (Figure 1B) and the TIMER 2.0 (Figure 1C). The PAAD cohort was characterized by a relatively high expression of MMP28. We also found that MMP28 was overexpressed in PAAD patients compared with the pericarcinomatous tissue (Figure 1D and E).

Figure 1 The mRNA expression of MMP28 in PAAD. (A) The MMP family expression in pericarcinomatous tissue and PAAD. The red squares represent higher expression and the blue squares represent lower expression. The darker the color is, the greater the difference is. Normal group means the pericarcinomatous tissue. (B) The MMP28 mRNA expression in different types of cancer (Oncomine). The red squares represent higher expression and the blue squares represent lower expression compared with normal tissue. The number in the legend below the figure means the percentage of MMP28 expressed differently in all cancers. The number inside the colored box represents the number of samples with statistically significant analyses for this type of cancer. (C) The MMP28 mRNA expression profile across all cancer types and normal tissues (TIMER 2.0). The red and blue dots represent cancer and normal tissues, respectively. TPM, transcripts per million. TPM is a method of standardizing the read counts of genes or transcripts in RNA-seq analysis. PAAD. Normal means pericarcinomatous tissue. (D) The mRNA expression of MMP28 in Pei Pancreas dataset (Oncomine). Pancreas means pericarcinomatous tissue. P <0.01, Fold change=1.5, gene ranking=top 10%. (E) The mRNA expression of MMP28 between PAAD and TCGA normal and GTEx data (GEPIA). The red square represents PAAD and the grey square represents TCGA normal and GTEx data. The red asterisk represents statistical significance. The cutoffs of P value and |log2 fold change (Fc)| were determined as 0.01 and 1, respectively. And log2 (TPM+1) is used for log scale. *P <0.05. Pancreas means the pericarcinomatous tissue.

Relationship Between the mRNA Expression of MMP28 and Clinicopathological Parameters of PAAD Patients

We used R to analyze PAAD data of the TCGA database. MMP28 expression was higher in PAAD patients aged 81–100 years compared with those aged 61–80 years (P<0.05, Figure 2D). Patients with chronic pancreatitis also had higher expression of MMP28 compared with those without chronic pancreatitis (P<0.05, Figure 2G). However, no significant difference of MMP28 expression was found between PAAD patients with different gender, tumor stages, smoking or drinking history, history of diabetes and tumor metastasis (Figure 2A–C, E, F and H).

Figure 2 The mRNA expression of MMP28 in PAAD based on individual cancer stages and in subgroups of patients with PAAD classified based on smoking, age, gender, and other criteria. (A) The mRNA expression of MMP28 in PAAD based on different cancer stages. (B) Boxplot showing relative mRNA expression of MMP28 in different smoking years of PAAD patients. (C) Boxplot showing relative mRNA expression of MMP28 in male or female PAAD patients, respectively. (D) Boxplot showing relative mRNA expression of MMP28 in PAAD patients at different ages. (E) Boxplot showing relative mRNA expression of MMP28 in PAAD patients with different drinking habits. (F) Boxplot showing relative mRNA expression of MMP28 in PAAD patients with or without diabetes. (G) Boxplot showing relative mRNA expression of MMP28 in PAAD patients with or without chronic pancreatitis. (H) Boxplot showing relative mRNA expression of MMP28 in PAAD patients with different nodal metastasis status. Normal group means the pericarcinomatous tissue. The central mark is the median; the edges of the box are the 25th and 75th percentiles. *P<0.05.

Abbreviation: NS, no significance.

High mRNA Expression of MMP28 is Associated with Poor Prognosis in PAAD Patients

We further explored the prognostic value of MMP28 in patients with PAAD. The data from the GEPIA were used to analyze the correlation between mRNA expression of MMP28 and prognostic efficiency of patients with PAAD. The Kaplan-Meier curves revealed that high MMP28 mRNA expression was significantly associated with low overall survival (OS) and disease-free survival (DFS) rates (P<0.001) (Figure 3A and B). Subsequently, we analyzed the genes with the highest prognostic potential related to OS (Table 1) and DFS (Table S2) in PAAD based on the GEPIA database. Disease free survival indicated the period during which the patient lived after the treatment of PAAD and no tumor recurrence occurred. Median value was used to define the low and high expression level of MMP28 for survival analysis. These results showed that MMP28 was a potential prognostic molecular marker of PAAD.

Table 1 Most Differential Survival Gene (Overall Survival) Top 10 in PAAD

Gene Symbol	Gene ID	P-value
ANKRD19P	ENSG00000187984.12	1.06e-6
MYEOV	ENSG00000172927.7	1.42e-6
MMP28	ENSG00000271447.5	2.76e-6
KRT16P3	ENSG00000214822.8	4.58e-6
PPEF2	ENSG00000156194.17	5.67e-6
ATP6V0E2-AS1	ENSG00000204934.10	6.02e-6
EFR3B	ENSG00000084710.13	7.59e-6
CH17-360D5.2	ENSG00000276850.4	8.01e-6
KCNC1	ENSG00000129159.6	8.15e-6
PLA2G16	ENSG00000176485.10	8.20e-6

Figure 3 The relationship between high mRNA expression of MMP28 and prognostic survival in PAAD (GEPIA). (A) Overall survival (OS). (B) Disease free survival (DFS).

GO Term and KEGG Pathway Enrichment Analysis of MMP28-Correlated Genes in PAAD

To define the possible biological significance of MMP28 in PAAD, we examined the MMP28-correlated mode in the TCGA-PAAD cohort using LinkedOmics database. 3402 genes positively-correlated with MMP28 and 3465 negatively-correlated genes were identified based on the P value for Pearson’s correlation with false discovery rate adjustment (FDR < 0.01) (Figure 4A). The top 50 significant genes positively and negatively correlated with MMP28 were then illustrated in heat maps (Figure 4B and C). A detailed description of these correlated genes is listed in (Table S3). The three most positively correlated genes were keratin 19 (KRT19) (r=7.487e-01, p=1.127e-27), Interleukin 1 Receptor Antagonist (IL1RN) (r=7.434e-01, p=4.117e-27), and Annexin A2 (ANXA2) (r=6.736e-01, p=8.857e-21), all of which are cancer-related genes.^31–33 These results suggested that MMP28 was linked to carcinogenesis as well.

Figure 4 MMP28 correlated genes in PAAD (LinkedOmics). (A) The overall MMP28 highly correlated genes identified by the Pearson test in the PAAD cohort. (B) Heat maps showing top 50 genes positively correlated with MMP28 in PAAD. (C) Heat maps showing the top 50 genes negatively correlated with MMP28 in PAAD. The red squares present positively correlated genes and the blue squares present negatively correlated genes. (D–G) GO term annotation and KEGG pathway enrichment analysis of MMP28 correlated genes. (D) cellular component (CC), (E) biological process (BP), (F) molecular function (MF), and (G) KEGG pathway.

Significant GO term annotation identified the functional roles of MMP28-correlated genes based on three aspects, including cellular components (CC), biological process (BP) and molecular functions (MF). The main CC were the cornified envelope, cell-substrate junction, and ECM (Figure 4D). These genes primarily participate in BP such as skin development, NADH dehydrogenase complex assembly, ephrin receptor signaling pathway, response to type I interferon, and isoprenoid metabolic processes (Figure 4E). The MF terms such as isoprenoid binding, peptidase regulator activity, extracellular matrix binding, steroid dehydrogenase activity, and calcium-dependent protein binding were also significantly regulated by MMP28-correlated genes (Figure 4F). KEGG pathway enrichment analysis revealed that proteasome-related pathways, base excision repair, oxidative phosphorylation, glycolysis, metabolism of xenobiotics by cytochrome P450, ribosome and necroptosis were closely related pathways (Figure 4G). Therefore, we can conclude that MMP28 is involved in a series of complex biological processes in PAAD.

Underlying Regulators of MMP28-Correlated Genes in PAAD

To further explore the underlying regulators of MMP28-correlated genes in PAAD, we estimated the enrichment of miRNA, TF, and kinases of these genes in the LinkedOmics database. The data were derived from the Broad Institute of MIT and Harvard and the University of North Carolina. The Pearson correlation test was used for statistical analysis. The enrichment of miRNA was mainly related to the microRNA (MIR)-181 family, including MIR-181A, MIR-181B, MIR-181C, MIR-181D (Figure 5A). In a recent study, MIR-181 family members expression was shown to be elevated in primary breast cancer, where these genes act as oncogenes.³⁴ The following three transcription factors were found enriched: activator protein 1 (AP1), ecotropic viral integration site 1 (EVI1) and organic cation uptake transporter 1 (OCT1) (Figure 5B). Previous studies have shown that these genes are all related to prognosis in different cancer type.^35–37 No significant kinase was found to be enriched. In essence, these findings suggested that MMP28 is tightly connected to prognosis and survival of PAAD patients.

Figure 5 The underlying regulators in miRNA, transcription factor, and kinase of MMP28 correlated genes in PAAD (LinkedOmics) and MMP28 genomic alterations in PAAD (cBioPortal). (A and B) came from LinkedOmics database, (C and D) came from cBioPortal database. (A) miRNA. miRNA means microRNA. (B) Transcription factor. AP1,activator protein 1; SRF,serum response factor; STAT, signal transduction and activator of transcription; TEF,transcript elongation factor; ELK1, Ets-like transcription factor-1; DR1, down-regulator of transcription 1; EVI1,ecotropic viral integration site 1; OCT1,organic cation uptake transporter 1; POU6F1, POU Class 6 Homeobox 1. (C) Oncoprint of MMP28 alterations in the PAAD cohort. (D) MMP28 expression frequency in different copy-number alteration types.

Genomic Alterations of MMP28 in PAAD

Subsequently, we used the cBioPortal tool to determine genomic alterations of MMP28 in PAAD based on the TCGA data. MMP28 was altered in approximate 10 of 185 (5%) patients with PAAD (Figure 5C). These alterations included mRNA high expression in about 7 cases (3.5%) and amplification (AMP) in 3 cases (1.5%). Thus, high mRNA expression was the most common type of MMP28 gene alteration in PAAD. As shown in Figure 5D, the box-and-whisker diagram illustrated the relationship between MMP28 mRNA expression and its copy number alteration type. Compared with the amplification group, the diploid group had higher MMP28 mRNA expression frequency. Finally, we defined genes correlated to MMP28 in PAAD (Table S4), which were used for further network analysis. In summed, MMP28 mRNA was found highly expressed in PAAD.

Construction of PPI and TF-miRNA Coregulatory Networks Based on MMP28-Correlated Genes in PAAD

To better understand the interplay between MMP28-correlated genes, we analyzed generic PPI networks of MMP28-correlated genes (both positively and negatively) based on STRING using the NetworkAnalyst database. Revealing functional interactions between proteins may provide insights into the mechanism behind the prognostic potential of MMP28 in PAAD. As shown in Figure 6A, the most significant module of the PPI network of positively-correlated genes was formed by an MCODE plugin in Cytoscape. The results indicated that the ribosomal protein L8 (RPL8), ribosomal protein L30 (RPL30), and ribosomal protein S16 (RPS16) were the top three positively correlated proteins. Then, we predicted the most significant module of the PPI network of negatively correlated genes (Figure 6B). Ribosomal protein L15 (RPL15), ribosomal protein S4 X-Linked (RPS4X), and ribosomal protein S3A (RPS3A) were selected as the top three negatively correlated proteins. Finally, we further studied the TF-miRNA coregulatory network of positively correlated genes, and the most significant module was presented in Figure 6C. The top five regulators were TATA-Box Binding Protein Associated Factor 1 (TAF1), TAF2, TAF5, TAF8, and cell division cycle 6 (CDC6). Generally speaking, the PPI networks demonstrated that the most strongly associated protein was a ribosomal protein contributing to cancer development, progression and metastasis.^38–40 In line, previous studies have shown that TAFs and CDC6 could promote cancer development and progression.^41–43 These results suggested that MMP28 might be involved in promoting pancreatic cancer development, progression, and metastasis.

Figure 6 PPI and TF-miRNA coregulatory networks based on MMP28 correlated genes in PAAD (NetworkAnalyst). (A) The most significant module of PPI network of positively correlated genes. (B) The most significant module of PPI network of negatively correlated genes. (C) The most significant module of TF-miRNA coregulatory network of positively correlated genes. Genes are represented as nodes and their interactions were denoted by lines. The color and size of the nodes represent degree values. The gene of darker color and greater circles show the higher degree values, whereas the lighter color and the smaller circles or triangles show the smaller degree values in these networks.

The Relationship Between MMP28 RNA Expression and Immune Cells in Blood

Furthermore, we investigated whether there was any link between MMP28 RNA expression and the presence of immune cells in blood, based on the HPA database (Figure 7A). The Blood Atlas in HPA, which included single cell type information on genome-wide RNA expression profiles of human protein-coding genes expressed by various B and T cells, monocytes, granulocytes and dendritic cells. The results showed that MMP28 RNA was specific for the T-cell lineage of blood cells, especially for naive CD4⁺ T-cells and naive CD8⁺ T-cells. The results suggested that high MMP28 expression could activate T cell-mediated cellular immunity. Interestingly, the mucosal-associated invariant T (MAIT) cell type was enhanced compared with other immune cell, and this cell type is a new hotspot in cancer immunotherapy, discovered only recently.⁴⁴ Therefore, these results showed a close correlation of MMP28 with immune cells in the blood.

Figure 7 Relationship between MMP28 and immune cells in blood and in tumor microenvironment. (A) MMP28 RNA specifically correlated immune cell types in blood (HPA). NX is the normalized expression of transcript expression values. (B) MMP28 copy number variety (CNV) affects the infiltrating levels of CD4+ T cell, B cell, CD8+ T cell, and dendritic cell in PAAD (TIMER). *P<0.05; **P<0.01; ***P<0.001.

MMP28 is Correlated with Immune Cells in the Tumor Microenvironment in PAAD

To better understand the relationship between MMP28 and TME in PAAD, we utilized TIMER and TIMER 2.0 to shed light on the influence of MMP28 expression on immune cells. The SCNA of MMP28 was associated with the infiltration level of CD4⁺ T cells, B cells, CD8⁺ T cells and dendritic cells in PAAD. Arm-level deletion was the most relevant factor related to CD4⁺ T cell infiltration level (Figure 7B). We analyzed the correlation between MMP28 expression with immune infiltration level in TME of diverse cancer types using the TIMER 2.0 database. The immune infiltration levels of B cells, CD4⁺ T cells, and CD8⁺ T cells in PAAD were almost negatively correlated with MMP28 expression in PAAD (Figure 8A–C). In addition, the number of endothelial cells was negatively correlated with MMP28 expression in PAAD, despite a positive correlation in other cancer types (Figure 8D). Figure 8E–H shows the correlation between MMP28 and naive B cells (CIBERSORT), naive CD4⁺ T cells (XCELL), naive CD8⁺ T cell (QUANTISEQ), and endothelial cells (XCELL) in PAAD. In short, high MMP28 expression might reduce the infiltration level of immune cells in the TME of PAAD, indicating that MMP28 could be a potential immune therapy target.

Figure 8 The correlation of MMP28 expression with immune infiltration level in TME in diverse cancer types (TIMER 2.0). (A) B cell immune infiltration level. (B) CD4+ T cell immune infiltration level. (C) CD8+ T cell immune infiltration level. (D) Endothelial cell immune infiltration level. The red square represents positive correlation, and the blue square represents negative correlation. The darker the color is, the stronger the correlation is. (E) Correlation between MMP28 and B cell naive CIBERSORT in PAAD. (F) Correlation between MMP28 and T cell CD4+ naive XCELL in PAAD. (G) Correlation between MMP28 and T cell CD8+ naive QUANTISEQ in PAAD. (H) Correlation between MMP28 and Endothelial cell XCELL in PAAD.

Knockdown of MMP28 Suppresses the Proliferation, Migration and Invasion of PAAD Cells

To investigate the influence of MMP28 on biological function of PAAD cells, siRNA was constructed to knockdown the MMP28 gene. The knockdown efficiency of MMP28-siRNA was defined by RT-qPCR and siRNA-1 showed the highest knockdown efficiency (Figure 9A). The knockdown efficiency of MMP28 protein level was also verified by Western blot (Figure 9B). MMP28 knockdown inhibited the proliferation, migration and invasion of PANC-1 and SW1990 cells in vitro (Figure 9C–F). These results suggested that MMP28 might promote proliferation, migration and invasion of cancer cells in PAAD.

Figure 9 Knock-down of the MMP28 gene reduced proliferation, migratory and invasive abilities of PAAD cells. (A and B) The knockdown efficiencies of MMP28 in PAAD cells were detected by RT-qPCR and Western blot. (C) CCK8 assays were performed to determine the influence of MMP28 knock-down on the proliferation abilities in PANC-1 and SW1990 cells. (D–F) Transwell assays were performed to determine the influence of MMP28 down-regulation on the migratory and invasive abilities in PANC-1 and SW1990 cells. ***P<0.001.

Abbreviation: NC, negative control.

Discussion

Analysis of MMP28 mRNA expression from several different public web resources and TCGA PAAD data confirmed that MMP28 mRNA expression was higher in PAAD than in other cancer types and pericarcinomatous tissue. The expression of MMP28 in PAAD was not be affected by most common risk factors. These observations indicate that MMP28 might constitute an advanced diagnostic marker for PAAD. Moreover, high mRNA expression of MMP28 was strongly linked to poor survival and disease-free state in many cohorts of patients. Overall, our results suggest that MMP28 needs further clinical validation as a potential prognostic biomarker for PAAD.

Furthermore, the three genes most-highly correlated with MMP28, including KRT19, IL1RN, ANXA2, are strongly expressed in many cancers and are related to cancer progression and prognosis.^31–33 GO term annotations like cornified envelope, cell-substrate junction, ECM, skin development, and extracellular matrix binding focused on the following biological behaviors: degradation and binding of ECM, skin formation, or keratin formation. Dysregulation of the proteasome contributes to several diseases, including cancer.⁴⁵ NADH dehydrogenase complex I is overexpressed in metastatic cells and not in non-metastatic cells. Prior to cancer metastasis, the increased NADH dehydrogenase complex can satisfy high energy demands.⁴⁶ Hence, the BP of the NADH dehydrogenase complex assembly might indicate that MMP28 could promote cancer metastasis. These findings are in accord with previously reported functions of MMP28, related to promotion of tumor invasion and metastasis. Finally, KEGG pathway analysis indicated that MMP28 might be involved in regulating the protein translation process by the ribosome, promoting protein degradation and inducing necroptosis, which is linked to poor prognosis. Furthermore, a recent study found that tumor cells themselves can adapt to changes in the metabolic environment by switching between oxidative phosphorylation and glycolysis.⁴⁷ Three intermediate filament proteins (keratins, vimentin, and nestin) are highly present in cancer.⁴⁸ Further research regarding related biological processes is still needed. In brief, biological functions of MMP28-correlated genes indicate that MMP28 might play a role in accelerating cancer.

Next, we analyzed the regulators of MMP28-correlated genes, including miRNAs, transcription factors, and kinases. While no significantly enriched kinase was detected, many miRNAs and TFs had been considered as regulators of MMP28-correlated genes. This might be explained by genomic alterations of MMP28 in PAAD, which we investigated subsequently. The main genomic alteration of MMP28 in PAAD was mRNA overexpression, therefore, the regulators were mainly related to miRNA and transcription factors. The function of protein kinases is to transfer phosphate groups from adenosine triphosphate (ATP) to a substrate, leading to its activation. This process often follows mRNA transcription and protein translation. MIR-181 family has been shown to impact tumor progression, invasion and metastasis, prognosis and survival in many cancers.^49–52 Our results suggest that the MIR-181 family is an important tumor regulator and that MMP28 might be related to the prognosis of PAAD via its relation with the MIR-181 family. Further research is needed to test this hypothesis. Finally, all of the three mainly enriched TF regulators, ie, AP1, EVI1 and OCT1, are closely related to cancer prognosis.^35–37 Thus, these regulators indicate that MMP28 is a prognostic biomarker in PAAD, and its potentially regulated pathways might include MIR-181 family, AP1, EVI1, and OCT1.

In our research, the main genomic alteration of MMP28 in PAAD was high mRNA expression, and the most common type of MMP28 mRNA copy number was a diploid. Previous studies revealed that, compared with normal diploid cells, aneuploid cancer cells are characterized by heterogeneous genomic landscapes, including sub-diploid, diploid and supra-diploid regions and higher gene copy number abnormalities. These increased mRNA levels of some single transcribed cancer genes (induced by the genetic instability) may restore functional haplo-insufficient mRNA in cancer.⁵³ Hence, MMP28 mRNA copy number was increased mainly through the formation of mRNA diploids, which might play a role in recovering the functionality of cancer cells in PAAD. Detailed related mechanisms remain to be further explored. Ribosomal proteins take part in composing the ribosome, which is responsible for protein synthesis in all living cells. The inherited mutations of the ribosome contribute to cancer development and progression.^39,40 Malfunctions of the ribosomal protein expression and translation mechanisms can promote breast cancer metastasis.³⁸ Recent studies reported that TAF1 plays a crucial role in leukemogenesis and high expression of the TATA-box binding protein (TBP) induces oncogene DNA damage and genomic instability in cancer.^42,43 Furthermore, CDC6 reduces the mitotic rate, which is related to drug resistance in cancer cells.⁴¹ In our research, PPI networks revealed that ribosomal proteins were strongly associated with MMP28. It indirectly indicated that MMP28 may produce numerous mRNA transcriptions, which is in accordance with the genomic alteration of MMP28. Then, the two main TF-miRNA coregulators of genes positively correlated with MMP28 were TAFs and CDC6. In essence, these related proteins and regulators suggest that MMP28 contributes to cancer progression and metastasis. TAFs and CDC6 could serve as potential MMP28 regulators on TF and miRNA levels.

MAIT cell is one type of the unconventional T cell, which is only presented in a particular molecular context. MAIT cell is characterized by a limited range of non-peptide antigens that are recognized via T cell receptors.⁴⁴ The presence of MAIT was correlated with the expression of pro-inflammatory cytokines and cancer initiation, growth, and metastases.^54,55 Furthermore, both MMP28 and MAIT cells are associated with promoting fiber formation and repair.^15,56 Here, we found that high MMP28 expression led to high expression of MAIT cell in blood. The mechanism and pathways related to the interaction between MMP28 and MAIT cells need to be further studied.

Immune cells in the TME play crucial roles in tumorigenesis. Anticancer immune cells infiltrating the TME could target and kill cancer cells at an early stage of tumorigenesis. However, cancer cells can escape from immune surveillance and even suppress the cytotoxic function of anticancer immune cells by a series of complex mechanisms. Immune evasion provides a new strategy of cancer immunotherapy.⁵⁷ High MMP28 expression was correlated with reduced anticancer immune cells infiltration, including B cells, naive CD4⁺ T cells, and naive CD8⁺ T cells in PAAD. Therefore, MMP28 can be regarded as a potential immune therapy target in PAAD. High MMP28 expression was also linked to a decrease in the number of endothelial cell in PAAD, compared with other cancers. This may be related to the strong invasiveness of PAAD.

The knockdown of MMP28 inhibited proliferation, migration and invasion of PAAD cells in vitro, which in turn proves that MMP28 is closely correlated with PAAD progression.

Conclusion

Our findings indicate that MMP28 is highly expressed in PAAD where it is associated with cancer progression, invasion, and metastasis. Hence, MMP28 might serve as an independent prognosis biomarker and a prospective therapeutic target for PAAD.

Abbreviations

PAAD, pancreatic adenocarcinoma; PDAC, pancreatic ductal adenocarcinoma; MMP28, matrix metalloproteinase 28; ECM, extracellular matrix; TME, tumor microenvironment; EMT, epithelial mesenchymal transition; TPM, transcripts per million; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; siRNA, small interfering RNA; GEPIA, gene expression profiling interactive analysis; GTEx, genotype tissue expression; HPA, human protein atlas; FDR, false discovery rate; PPI, protein-protein interaction; STRING, search tool for the retrieval of interacting genes; TF, transcription factor; MCODE, molecular complex detection; SCNA, copy number alteration; TIICs, tumor-infiltrating immune cells; DMEM, dulbecco’s modified eagle medium; FBS, fetal bovine serum; PVDF, polyvinylidene difluoride; OS, overall survival; DFS, disease-free survival; CC, cellular components; BP, biological process; MF, molecular functions; MIR, microRNA; KRT19, keratin 19; IL1RN, interleukin 1 receptor antagonist; ANXA2, annexin A2; AP1,activator protein 1; SRF, serum response factor; STAT, signal transduction and activator of transcription; TEF, transcript elongation factor; ELK1, Ets-like transcription factor-1; DR1, down-regulator of transcription 1; EVI1, ecotropic viral integration site 1; OCT1, organic cation uptake transporter 1; POU6F1, POU Class 6 Homeobox 1; AMP, amplification; RP, ribosomal protein; TAF, TATA-Box Binding Protein Associated Factor; CDC, cell division cycle; MAIT, mucosal-associated invariant T; ATP, adenosine triphosphate; TBP, TATA-box binding protein; CCK8, cell counting kit-8.

Data Sharing Statement

The datasets used during the present study are available from different public databases.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Acknowledgments

We would like to thank the Science Park of West China Hospital, Sichuan University (Chengdu, China) for the help in acquisition and culturing of the PAAD cell lines.

Disclosure

These authors declare that they have no competing interests.

Additional information

Funding

This work was supported by the following Grants: National Natural Science Foundation of China (81572663, 82071125, 81972551, 81702701, 82071137, 81602375), the CAMS Innovation Fund for Medical Sciences (CIFMS, 2019-I2M-5-004), 111 Project of MOE (B14038), China and the Sichuan Science and Technology Program (2020YJ0286).