MMP12 serves as an immune cell–related marker of disease status and prognosis in lung squamous cell carcinoma

Data acquisition and data processing

The high-throughput data related to LUSC for this study were obtained from multiple publicly available databases: the Gene Expression Omnibus (GEO), ArrayExpress, The Cancer Genome Atlas (TCGA), and Genotype-Tissue Expression (GTEx) databases. Among them, microarray and RNA-Seq data could be obtained directly from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Data for patients diagnosed with LUSC were downloaded from TCGA (http://xena.ucsc.edu/) and included RNA-Seq and clinical information (including clinical characteristics and survival data). A total of 33 criteria-compliant data sets with 2,738 samples, including 1,348 LUSC samples and 1,390 non-LUSC samples, were included in this study (Supplementary Material 1). In addition, we obtained four datasets (GSE37745, GSE29013, GSE30219, and GSE73403) involving a total of 221 LUSC cases to analyze whether MMP12 expression was associated with overall survival (OS) in LUSC patients.

The “limma” package (Ritchie et al., 2015) was used to process the data sets from public databases and to normalize them using log₂ (x + 1). Since the experimental batches differed for each data set, this study combined 33 LUSC cohorts into 10 new fellows based on the same platform. For example, GSE31552 and GSE44077 were merged into the reorganization cohort GPL6244 since both data sets were from the same platform: GPL6244. Finally, the “SVA” package (Leek & Storey, 2007) was used to eliminate batch effects.

Evaluation of MMP12 protein expression in LUSC tissues based on internal IHC staining

To clarify the differences in MMP12 protein expression between LUSC and non-LUSC patients, we collected 125 tissue samples (106 LUSC samples; 19 control samples) from the First Affiliated Hospital of Guangxi Medical University. A total of 125 internal samples were divided into three tissue microarray sections (Nos. LUC481, LUC1021, and LUC1502, constructed by Guilin Fanpu Biotech, Guangxi Zhuang Autonomous Region, China) for further IHC experiments. The relevant operations were performed according to the manufacturer’s standards. First, tissue sections were fixed by formalin and embedded in paraffin, then sections were dewaxed, and antigen repair was performed with ethylenediaminetetraacetic acid buffer (pH = 9.0) at boiling state. After inactivation of endogenous enzyme activity, rabbit antihuman MMP12 monoclonal antibody (EPR11944[B], ab170414, dilution ratio 1:100) was added and stored overnight at 4 °C ambient. Next, a second antibody-labeled horseradish peroxidase (D-3004-15, Changdao Biotechnology Co., Ltd., Shanghai, China) was added to the tissue sections and incubated at 25 °C for 30 min. Subsequently, the slides were stained with diaminobenzidine, taking care to rinse the slides with PBS between steps. After hematoxylin restaining, the slides were dehydrated with a gradient concentration of ethanol solution and finally sealed with neutral gum. The sections were scored by microscopy, with brown representing positive and blue representing negative. The staining scores were judged as follows: 0 (no staining), 1 (mild staining), 2 (moderate staining), and 3 (severe staining). Different numbers of positive cells represented different scores: 0 (number <5%), 1 (5%–25%), 2 (26%–50%), 3 (51%–75%), and 4 (number > 75%). The total IHC score was generated by the product of the sample staining score and the percentage of positive cells, and the above assessments were performed independently by two senior pathologists. The raw data of IHC staining were presented in Supplementary Material 2. This study was approved and supported by the Medical Ethics Review Committee of the First Affiliated Hospital of Guangxi Medical University (2021 [KY-E-246]).

MMP12 expression and its clinical relevance in LUSC

The ability of MMP12 expression to discriminate between tumor and control tissues could be determined using receiver operator characteristic (ROC) curves and summary ROC (sROC) curves based on sensitivity, specificity, and area under the curve (AUC) values. It was generally considered that larger sensitivity, specificity, and AUC values indicate a more remarkable ability of MMP12 to differentiate between cancer and control tissues. Based on Yoden Index, true positive rate, false positive rate, true negative rate, false negative rate and cut-off values were calculated for each data set , and Stata (v15.0) was used to plot ROC curves.

In addition, this study further analyzed the association between MMP12 expression and the survival outcomes of LUSC patients by establishing Kaplan–Meier curves and univariate and multivariate Cox regression analyses. The clinical indicators (gender, age, tumor stage, node stage, metastasis stage, clinical stage, and MMP12) and prognostic information utilized in the above analyses were available from the Xena database.

MMP12 potential mechanism in LUSC

Cistrome Data Browser (Zheng et al., 2019) contains a vast number of human transcription factors (TF), histone modifications, and chromatin accessibility samples. It is currently one of the most comprehensive chromatin immunoprecipitation sequence (ChIP-Seq)-related databases. To explore the potential molecular mechanism of MMP12 in LUSC, this study predicted the TFs of target genes based on the Cistrome Data Browser and the standardized mean difference (SMD). In addition, ChIP-Seq data in the Cistrome Data Browser were used to validate our results and IGV (2.16.0) was used to visualize them. ChIP-Seq data must pass all quality tests and target gene score > 1. JASPAR (Fornes et al., 2020) is a free and publicly available TF database that focuses on collecting relevant TFs with DNA-binding site motifs. Therefore, the relevant TF motif data and the potential promoter region of MMP12 in this study could be directly obtained from JASPAR and NCBI, respectively. Subsequently, the FIMO tool in the MEME suite was selected to mine the upstream transcription start site (TSS) of MMP12 and the potential binding sequences of these motifs.

Gene Set Enrichment Analysis (GSEA) was performed using the “clusterProfiler” package (Wu et al., 2021), hoping to further explore the potential signaling pathways that MMP12 may involve in LUSC. The signaling pathways shown in this analysis were obtained from the Kyoto Encyclopedia of Genes and Genomes database. To predict the relationship between the core genes of the enriched pathways, an interaction database platform, STRING (v11.0), was used to construct a protein–protein interaction (PPI) network (Crosara et al., 2018). This search was performed with “Homo sapiens” as the set species, the confidence score was set to 0.4, and other settings were set as default.

Correlation of MMP12 expression with the tumor microenvironment (TME)

Based on the TCGA-GTEx and GPL570 data sets, a single sample Gene Set Enrichment Analysis (ssGSEA) was used to analyze the differences in multiple immune cell scores between the high and low MMP12 expression groups. Subsequently, we analyzed the correlation between MMP12 and individual immune cell scores one by one and initially assessed the potential of MMP12 for immune prediction.

The Tumor Immune Estimation Resource (TIMER) (Li et al., 2020; Chen et al., 2022c; Liu et al., 2022a; Li et al., 2022b) is an online database that allows comprehensive detection of immune cell infiltration levels. The data for six immune-infiltrating cells were downloaded from the database. The six immune cells included B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells. In addition, to explore MMP12 expression levels in various cells (e.g., monocytes/macrophage), the Cancer Single-cell Expression Map (Zeng et al., 2022) was used to systematically investigate the heterogeneity of MMP12 in different cells.

Statistical analysis

The Wilcoxon rank-sum test and the SMD were used for comparisons between various groups. Among them, the Wilcoxon rank-sum test should be used when the combined data set was skewed and the variance of the totals to which the two samples belong was different. When the results show I² values < 50% or chi-square test p - values > 0.05, there was less heterogeneity between the data sets, and we need to use a fixed-effects model for SMD calculations. The results of the SMD were statistically significant when the corresponding 95% confidence interval (CI) did not contain 0. The following criteria were used to determine differentially expressed genes: high expression: SMD > 0 and 95% CI excluding 0; low expression: SMD < 0 and 95% CI excluding 0. Begg’s test and Egger’s test evaluated the publication bias of SMD results, and p < 0.1 indicated a significant publication bias. Sensitivity analysis was used to evaluate the reliability of the analysis results. For hazard ratio (HR), the presence of a 95% CI of 1 or p > 0.05 indicated no statistical significance. In addition to the above statistical methods, the Spearman correlation coefficient was used to detect the correlation between MMP12 and immune cells.

In this study, p < 0.05 represented a statistically significant difference. Stata (v15.0) was used for Egger’s tests, and sensitivity analysis, and R (v4.1.0) was used to complete all the remaining computational and visualization steps. The following was the research flow of this study (Fig. 1A).

MMP12 expression level differences between LUSC and its control tissues

In the 10 data sets included in this study, statistically significant differences in MMP12 expression between the LUSC group and control tissues were detected in nine cohorts (except for “GSE6044”). The results showed that MMP12 mRNA expression levels were increased in the LUSC group in nine of the combined data sets mentioned above (p < 0.05; Fig. 2A). Then, a random-effects model showed that MMP12 mRNA expression was elevated in the LUSC group but not in the control group (SMD = 3.13, 95% CI [2.51–3.75]; Fig. 2B), and no significant publication bias was detected (p > 0.1; Fig. 2C and Supplementary Material 3). Sensitivity analysis showed a small effect of low sample size on the overall MMP12 expression in LUSC and a large heterogeneity of the data (Supplementary Material 4). We believed that the greater heterogeneity was related to the data coming from different experimental platforms, and a random-effects model was used in this study.

Furthermore, the expression level of MMP12 in LUSC was further verified using internal immunohistochemical experiments. Compared to the non-LUSC group, MMP12 protein levels were significantly upregulated in the LUSC group samples (p < 0.001; Fig. 2D), consistent with the mRNA expression levels. Further, under the microscope, as seen in Figs. 3A–3D, MMP12 protein could be easily observed in LUSC tissues (Figs. 3B and 3D) but not in control tissues (Figs. 3A and 3C).

A potential molecular mechanism for the high expression of MMP12 in LUSC

As shown in the flow of Fig. 3E, based on the Cistrome Data Browser for calculation, we selected 11 TFs with a selection score greater than or equal to 0.5. Subsequently, the results of the SMD calculation and co-expression analysis determined that SNAI2 was a potential TF regulating MMP12 expression and positively correlated with MMP12 expression. The peak signal of SNAI2 ChIP-Seq was located within 2 kb upstream of MMP12, further validating the above conclusion (Supplementary Material 5).

The motif of SNAI2 was shown in Fig. 3F. Combined with JASPAR and FIMO tools, we obtained the binding sequence of SNAI2 upstream of MMP12 TSS:GCACCTGTC. This study also identified two potential binding sites between the SNAI2 motif and the TSS of MMP12 (Supplementary Material 6).

MMP12 might affect immune-related signaling pathways

In this study, GSEA was used to explore the potential mechanisms of MMP12 in multiple signaling pathways. The peaks of the curves appeared in the high-expression group, implying that these pathways were positively correlated with MMP12 expression and that these signaling pathways became more active when MMP12 expression was upregulated. Based on the LUSC samples from the data sets TCGA-GTEx (Fig. 4A) and GPL570 (Fig. 4B), GSEA suggested that MMP12 was associated with multiple signaling pathways, and the results of the two data sets overlapped to obtain seven common pathways (Fig. 4C). Notably, two pathways (the T-cell receptor signaling pathway and the natural killer cell–mediated cytotoxicity) showed a direct association of MMP12 with immune cells, suggesting that MMP12 was likely to be involved in essential mechanisms of cancer through immune-related signaling pathways (Supplementary Material 7).

Subsequently, for the core genes enriched to both pathways, differential expression of most genes was observed between the high and low MMP12 expression groups, and all of these genes were upregulated in the high MMP12 expression group (p < 0.05; Figs. 4D and 4E). Finally, this study performed a PPI analysis of the core genes enriched in both pathways and found that these genes closely interacted with each other (Fig. 4F). The black and green connecting lines between MMP12 and tumor necrosis factor (TNF) and CSF2 suggested their co-expression and previous studies between the genes, respectively. We further explored previous studies finding that MMP12 affects the activity of the TNF signaling pathway and monocyte-derived MMP12 is induced by CSF2 (Björnfot Holmström et al., 2017; Mao et al., 2022). Furthermore, we found that LCK proto-oncogene (LCK) had the highest number of PPI nodes with other genes, implying that LCK might be the central gene in both signaling pathways and had multiple different interactions with other genes, followed by PIK3R1 and TNF.

MMP12 positively correlated with immune cell levels

LUSC patients in the TCGA-GTEx cohort were divided into high and low MMP12 expression groups using median MMP12 expression. The correlation between MMP12 expression and 28 immune cells was investigated using ssGSEA. The results showed significant differences in immune cell scores between the high and low MMP12 expression groups, and higher levels of immune cells could be observed in the high MMP12 expression group (p < 0.01; Fig. 5A). Considering that the association between them could still not be intuitively determined by difference analysis alone, we further explored the correlation between MMP12 and these immune cell scores. As can be found in Fig. 5B, the expression level of MMP12 tended to be positively correlated with the immune cell score (ρ > 0.15, p < 0.01). Most of these results were validated by the GPL570 data set (Figs. 6A and 6B).

MMP12 might serve as an antigen-presenting cell–associated LUSC antigen

Considering that the immune response plays a crucial role in antitumor, this study further verified the potential association between MMP12 expression and immune cells in TME using different methods. Based on the TIMER algorithm (Fig. 7A), we found that MMP12 expression was significantly negatively correlated with tumor purity (ρ < 0, p < 0.001) and that upregulation of MMP12 expression was associated with increased infiltration levels of five immune cell types (except CD4+ T cells) (ρ > 0, p < 0.05). From the previous results, it was clear that MMP12 in LUSC was positively correlated not only with the level of specific immune cells (e.g., CD 8 + T cells) but also with the level of non-specific immune cell infiltration (e.g., neutrophils), suggesting that it might have the potential to activate the immune system.

Further analysis showed that multicenter single-cell sequencing data results also suggested that MMP12 expression was evident in immune cells in several classes of TME-related cells, such as dendritic cells (DC) and monocytes/macrophages (Figs. 7B–7D and Supplementary Material 8). This suggested that the high expression of MMP12 in LUSC might act as an antigen-presenting cell–associated tumor neoantigen and activate the body’s immune response.

MMP12 discriminated between LUSC tissues and normal lung tissues

The important clinical value of MMP12 in TME has already been observed, so this study expected that MMP12 could further demonstrate a valuable clinical role in LUSC by ROC and sROC curves. In the included 10 combined data sets, MMP12 mRNA expression possessed a very high precision to distinguish LUSC samples from control samples according to the ROC curves of nine of the cohorts (AUC > 0.9) (Fig. 8A). The sROC analysis showed that MMP12 expression made it feasible to distinguish LUSC samples from control tissue samples (sensitivity = 0.96, specificity = 0.94, AUC = 0.98) (Fig. 8B). These results suggested that MMP12 had significant potential in the identification of LUSC.

Correlation of MMP12 with the prognosis of LUSC patients

As shown above, MMP12 might be a potential marker for the identification of LUSC. Therefore, based on GEO and TCGA data, this study explored the relationship between MMP12 expression and the prognosis of LUSC patients. As shown in Figs. 9A and 9C, high MMP12 mRNA expression was associated with poor OS, and the OS of LUSC patients with MMP12 overexpression was poorer. A trend showed that overexpressed MMP12 could be a prognostic risk factor for LUSC patients, and only results based on the GSE37745 cohort were statistically significant (p = 0.026). Subsequently, the four GEO data sets were further combined and it was also demonstrated by prognostic analysis that overexpression of MMP12 leads to poorer prognosis in LUSC patients (Fig. 9B). The four GEO data sets covered 116 samples; compared to patients with reduced MMP12 expression, patients with elevated MMP12 expression tended to have a poorer OS (HR = 1.66, 95% CI [1.01–2.75]).

Further, the effect of different clinical indicators (gender, age, tumor stage, node stage, metastasis stage, pathologic stage, and MMP12 expression level/value on the OS of LUSC patients was investigated by establishing Cox regression analysis. The results showed that both increased age and upregulated MMP12 protein expression resulted in more adverse OS in LUSC patients (p < 0.05; Fig. 9D), and the multivariate Cox regression analysis showed similar results (p < 0.05; Fig. 9E).