MMP dysregulation is associated with relapse-free and overall survival

Given the relevance of MMPs in the first step of the metastatic process, we explored the differences in mRNA expression levels between non-metastatic TNBC and normal tissue samples in The Cancer Genome Atlas (TCGA) cohort. Seven MMPs were downregulated and ten MMPs were upregulated in TNBC tumors (n = 90) compared to normal tissues (n = 99, Figure S1a). Interestingly, six consecutive MMP genes with increased expression in TNBC tumors are encoded on the same genomic locus (Fig. 1a), in a region that harbors nine different MMP genes, located at Chr11q22.2 (Figure S1b). Furthermore, coexpression patterns were observed in TNBC tumors between genes located on each side of the MMP region, revealing two different expression clusters (Fig. 1b), which were defined as 5’MMP region (containing MMP7, MMP8, MMP20, and MMP27; in addition to neighbor non-MMP genes TMEM123 and BIRC3) and 3’MMP region (containing MMP1, MMP3, MMP10, MMP12, and MMP13). The ratio of the expression levels between MMP genes located at Chr11q22.2 (3’MMPs/5’MMPs ratio) was significantly higher in TNBC (Figure S1c).

Fig. 1

MMP expression and its clinical relevance in TNBC. a. Gene expression of MMP genes located on MMP region (chr11q22.2) comparing normal (N = 99) and TNBC (N = 90) patient-derived samples from TCGA consortium. Mann-Whitney test. ns: no significant, *P < 0.05, **P < 0.01. b. Correlation matrix for the expression genes located at the MMP locus located genes (P < 0.01). c. Gene distribution at MMP region (top). Genes were classified according to their location as 5’ or 3’. Gene location was used to define MMP signatures to perform Kaplan-Meier curves of relapse-free survival (left) and overall survival (right)

We wondered whether MMP expression differences might be a consequence of chromosomal or regional copy number alterations (CNAs). Despite TNBC and HER2 + exhibiting significant changes in copy number compared to the HR + subtype (Figure S2a), this alteration was not correlated with the expression of genes located at the MMP locus (Figure S2b). We also explored the potential influence of tumor heterogeneity on the MMP expression signatures. Using transcriptome deconvolution to estimate the cellular composition of each tumor tissue [34], we found that none of the MMP signatures displayed a correlation with stromal content. 5’MMPs and 3’MMPs signatures showed a poor correlation with immune infiltration (r = 0.23 and 0.27, respectively) and TME content (r = 0.26 and 0.33, respectively). Importantly, the 3’MMPs/5’MMPs ratio signature was independent of stroma content, immune cell score, and microenvironment content (Figure S2c).

Thus, we hypothesized that the alternative expression of these two clusters might be associated with clinical outcomes. In the survival analyses of patients with TNBC (n = 417), we observed that tumors with a higher 3’MMP/5’MMP ratio have a significantly shorter relapse-free survival (RFS; Fig. 1c, [log rank P = 0.023; hazard ratio (HR) = 1.57, 95% CI = 1.06 − 2.33]) and overall survival (OS, Fig. 1c, [log rank P = 0.005; HR = 2.65, 95% CI = 1.31 − 5.37]). We also conducted survival analyses considering the expression of 5’MMPs alone or 3’MMPs alone instead of the ratio. We observed that while the expression of 5’MMP genes was associated with better RFS and OS, the expression of 3’MMP genes did not show an association with survival (Figure S3a). Importantly, the evaluation of patients affected by other breast cancer subtypes showed that the 3’MMPs/5’MMPs ratio was also associated with a worse prognosis in hormone receptor-positive breast cancer, but not in patients with HER2-positive disease. Interestingly, when considering non-breast solid tumors, this signature was associated with a worse prognosis in liver cancer, lung adenocarcinoma, and sarcoma, whereas it correlated with better survival in gastric cancer (Figure S3b).

An insulator element near the MMP8 gene promoter region is involved in regional MMP gene regulation

Considering the expression pattern exhibited by the MMP genes at Chr11q22.2, we explored the distribution of gene regulatory elements in this region (Figure S4a). We found 29 potential insulator elements (IEs) and 13 potential enhancer elements (EEs) at Chr11q22.2 (Fig. 2a). The IEs play a pivotal role in the topologically associating domain (TAD) formation, being DNA elements recognized by the CCCTC-binding factor (CTCF) and contributing to chromatin loop formation [35]. We explored potential TADs at Chr11q22.2 using data from CTCF ChIA-PET (Figure S4b). This data suggested that the IE located between the MMP8 promoter region and the MMP10 gene body (Chr11:102,732,800 − 102,733,900; hg38), hereinafter IE8, demarcates a boundary for a TAD contributing to the expression signature observed by MMP genes at Chr11q22.2. This is supported by data from eQTLs since we observed that SNPs tend to modulate the expression of genes located on the same side of TADs at Chr11q22.2 (Figure S4c). It is worth noting that the CTCF binding site of the IE8 was depleted in mutations according to data from the Pan-Cancer Analysis of Whole Genomes (PCAWAG) database (N = 2,658) and whole genome sequencing of TNBC samples (n = 237) (Figure S4d), indicating that its normal function is required in cancer.

Fig. 2

Identification and edition of gene regulatory elements at chr11q22.2. a. Scheme of identified enhancer and insulator elements. CRISPR/Cas9-disrupted Insulator element 8 (IE8) is highlighted. b. Expression profile of 5’ and 3‘MMP genes in TNBC cell lines. c. Hi-C contact frequency matrix for the 3 Mb genomic region surrounding IE8 binned at 10-kb resolution. d. CTCF binding is represented as fold enrichment (top) relative to isotype control in MDA-MB-231 and MDA-MB-436 and (bottom) relative fold enrichment after IE8 disruption. Student’s T-test. *P < 0.05, **P < 0.01

Construction of cell models to explore the impact of IE8 activation on MMP gene expression

To select cell models that resemble the expression patterns detected in clinical specimens, we classified TNBC patients based on the expression of MMP genes into high 3’MMP/5’MMP ratio, intermediate 3’MMP/5’MMP ratio, and normal-like MMP expression profile (Figure S5a). Following, we profiled the expression levels of MMP genes in a panel of well-characterized TNBC cell lines (n = 17; Fig. 2b) and selected two cell lines that recapitulate the most predominant MMP abnormal expression patterns in TNBC clinical specimens: high 3’MMP/5’MMP ratio (MDA-MB-231) and intermediate 3’MMP/5’MMP ratio (MDA-MB-436). Chromatin interactions were established using the Hi-C method in the two cell models. Interestingly, suggesting an active insulator role, we observed that IE8 colocalizes with TAD boundaries in both models (Fig. 2c), in concordance with data from ChIA-PET.

Both TNBC cell lines were transiently transfected with a Cas9-containing plasmid and a sgRNA to stably disrupt the CTCF binding motif on IE8 using CRISPR/Cas9 technology. A single clone for each condition was selected to perform further experiments (Figure S5b). The CTCF binding ability to IE8 was tested in both cell lines through Cleavage Under Targets and Release Using Nuclease (CUT&RUN) followed by qPCR [36]. Basal levels of CTCF occupancy were higher in MDA-MB-231 than in MDA-MB-436, showing the biological variability reflected in clinical specimens. After IE8 disruption, a significant decrease in CTCF binding was observed in both cell lines (Fig. 2d).

In addition, we performed Hi-C on the cell models after IE8 disruption. We found that disruption of IE8 did not change the higher-order TAD chromosomal organization (Figure S6a). However, the distribution of Hi-C interaction signals indicates that IE8 interacts with 3’MMP distal regulatory regions in both models (Figure S6b). Interestingly, some of these interactions contain additional IEs disposed of in a convergent orientation to IE8, which is a requirement for IE interactions [37]. Thus, we identified that insulation induced by IE8 activation preferably involves IEs located at the 3’MMP region of the MMP locus Chr11q22.2 (Figure S6b). Apart from regional interactions, we also explored high-confidence interchromosomal interactions in both models (Figure S6c). Interestingly, we found interchromosomal interactions between the two ends of the MMP locus and super-enhancer elements located in different chromosomes (Figure S6d).

Disruption of IE8 leads to local chromatin accessibility changes

To address the implications of the IE8 disruption on chromatin accessibility, we performed ATAC-seq on our TNBC cell line models. 34,047 common peaks between all replicates were identified in MDA-MB-231, whereas 35,347 common peaks were found in MDA-MB-436. Accessibility analysis was performed by assessing differential accessibility peaks, which were defined by either their presence in only one condition (WT or IE8 dis) or by the significant change in intensity of shared peaks between conditions (Figure S7a). Thus, we identified 3,083 and 2,232 regions that gained and lost accessibility upon IE8 disruption in MDA-MB-231, respectively. Regarding MDA-MB-436, 8,673 and 370 regions were more and less accessible after IE8 disruption, respectively (Fig. 3a). Importantly, a significant overlap between the differentially accessible regions (n = 1,033) of both cell line models was observed.

Fig. 3

Local consequences of IE8 disruption. a. (top) Description of differentially accessible regions in MDA-MB-231 and MDA-MB-436 upon IE8 disruption and (bottom) representation of commonly modulated regions in both cell lines. b. (top) ATAC-seq peak intensity signals of all the TSS present in the whole genome, chr11 and 10 MB around IE8 (middle) Heatmap of active TSS in TNBC models and (bottom) variations on the accessibility of promoter regions of 5’ MMPs (purple) and 3’MMPs (green) regions in MDA-MB-231 and MDA-MB-436 before and after IE8 disruption. c. Illustrative examples of variation in accessibility upon IE8 disruption (IE8dis minus WT) in the IE8 region (grey), the 5’ MMP region (purple), and the 3’MMP region(green)

Despite changes in chromatin accessibility being detected across the genome, we examined whether they were also specifically enriched around IE8. We explored the differentially accessible regions located on chromosome 11 – where our region of interest is located – as well as in different width windows around IE8 (± 10 MB to ± 0.5 MB). Significant enrichment in the number of differentially accessible regions was observed in the IE8 genome vicinity whereas no differences were observed across chromosome 11 (Figure S7b).

Moreover, we aimed to identify potential gene regulatory elements associated with differential chromatin accessibility reprogramming after IE8 disruption. The relative abundance of promoter, enhancer, and insulator elements that exhibited differential accessibility after IE8 disruption was similar between MDA-MB-231 and MDA-MB-436 (Figure S7c). We focused on changes in chromatin accessibility around (± 2 kb) the gene transcription start sites (TSS). As expected, the resulting heatmap displayed a similar accessibility profile on all conditions with increased peak density on the TSS (Fig. 3b). Similar profiles were also observed between all conditions when peaks were centered in insulators and enhancers (Figure S7d). However, we found interesting differences when we focused on the TSS of the genes located at the Chr11q22.2. Promoters in the 5’MMP region were more accessible upon IE8 disruption, but no changes were observed in promoter regions located toward the 3’MMP region (Fig. 3b and c). In this regard, we also found changes in the accessibility of enhancer elements in both 5’ and 3’ MMP regions (Figure S7e). We additionally identified a decrease in chromatin accessibility at the IE8 CTCF binding site disrupted by CRISPR/Cas9 (Fig. 3c). In agreement with a higher CTCF occupancy detected by CUT&RUN (Fig. 2d), the decrease in chromatin accessibility of this site was more evident in the MDA-MB-231 cells than in the MDA-MB-436 cells.

IE8 disruption modulates regional MMP expression patterns

We determined whether the observed alterations in chromatin accessibility resulted in differential expression of the MMP genes near IE8. We, therefore, assessed the mRNA expression levels of MMP genes in the wild-type and IE8-disrupted cells through RNA-seq. We observed 237 mRNAs with a significant differential expression upon IE8 impairment (Fig. 4a). We performed a gene ontology (GO) analysis to classify the 166 significantly upregulated and 71 downregulated genes by biological process. Importantly, we found significant enrichment of extracellular matrix organization or Ca2+-dependent cell-cell adhesion (Figure S8a). In addition, pathways associated with fatty acid import were also upregulated, which has been confirmed as a feature of TNBC tumors with better prognosis [38]. Therefore, we assessed the functional impact of these gene expression changes by evaluating fatty acid uptake ability. Importantly we found that after IE8 disruption, both cell lines increased fatty acid uptake (Figure S8b). Since the CRISPR/Cas9-mediated IE8 disruption occurred in chromosome 11, we assessed whether significant changes were enriched on this chromosome. We did not observe significant variation in gene expression when considering all the genes, but the genes located 1 Mb around IE8 were significantly modulated after IE8 disruption in MDA-MB-231 (Fig. 4b).

Fig. 4

mRNA expression changes after IE8 disruption. a. Volcano plot summarizes the RNA-seq results in MDA-MB-231. 166 genes were significantly upregulated (Padj < 0.05) (blue dots), whereas only 71 were downregulated (red dots). b. Representation factor of differentially expressed genes in MDA-MB-231. Hypergeometric test. ns P > 0.05. c. Variation of MMP8 and MMP1 RNA expression in MDA-MB-231 in RNA-seq. Student’s T-test. *P < 0.05, **P < 0.01. d. MMP1 (top) and MMP8 (bottom) mRNA expression levels determined by qPCR in MDA-MB-231 and MDA-MB-436 models. Student’s T-test. ns P > 0.05, *P < 0.05, **P < 0.01 e. Variation of accessibility in enhancers close to MMP2 (top) and MMP9 (bottom) promoters. f. MMP2 (top) and MMP9 (bottom) mRNA expression levels determined by qPCR in MDA-MB-231 and MDA-MB-436 models. Student’s T-test. ns P > 0.05, *P < 0.05, **P < 0.01

Regarding the Chr11q22.2 encoded genes; the modification of regional gene expression patterns supports the promoter accessibility changes (Fig. 3c). RNA-seq revealed two interesting changes at the mRNA level, an increase in MMP8 and a decrease in MMP1 (Fig. 4c) after IE8 disruption. These alterations were confirmed by qPCR (Fig. 4d). MDA-MB-231 showed a shift between MMP1 and MMP8 after the IE8 disruption, exhibiting a decrease in the pro-invasive enzyme MMP1 and an increase in MMP8, associated with antitumor activity. Although variations were not statistically significant, RNA expression levels displayed a similar tendency in MDA-MB-436 (Fig. 4d). Thus, the ratio between MMP1 and MMP8 is decreased after IE8 disruption, resembling the profile observed in healthy breast samples as opposed to tumor samples (Figure S8c). Importantly, changes in MMP expression after the switch after IE8 disruption are not associated with the reprogramming of interchromosomal interactions (Figure S6d).

We further explored whether reprogramming of the MMPs located in the genomic vicinity of IE8 triggered any compensatory mechanism either modulating the expression of tissue inhibitor of metalloproteinases (TIMP) genes, a four-member family that balances the MMP activity or altering the RNA levels of two other relevant metalloproteinases, MMP2 and MMP9. TIMP1, TIMP2, TIMP3, and TIMP4 did not show significant changes at the mRNA level upon IE8 disruption (Figure S8d). IE8 depletion was not associated with changes in chromatin accessibility at MMP2 or MMP9 promoter regions. However, an important increase in accessibility was reported on enhancer elements located 8 kb and 15 kb upstream of TSS of the MMP2 and MMP9 genes, respectively (Fig. 4e). These changes were translated into a cell line-dependent increase in MMP mRNA levels. MMP2 was upregulated only in MDA-MB-231 whereas MMP9 was increased in MDA-MB-436 after IE8-disruption (Fig. 4f).

We considered whether changes in chromatin accessibility and the concomitant modulation of gene expression may also be observed in TNBC patient samples. We explored the six TNBC samples with ATAC-seq and mRNA-expression data available at TCGA. We found different levels of chromatin accessibility at IE8 (Figure S8e). Interestingly, we could correlate these changes to variations in 3’MMP/5’MMP expression. We observed that those patient-derived samples with higher levels of accessibility at IE8 showed a higher ratio of 3’MMPs/5’MMPs (r = 0.87, p = 0.02, Figure S8f).

IE8 disruption interferes with MMP1 release and decreases invasive potential in Breast cancer

We observed that IE8 disruption triggered a significant increase in MMP8 protein levels both in MDA-MB-231 and MDA-MB-436 (Fig. 5a). Conversely, a decrease in both MMP1 abundance and MMP1 activity was observed upon IE8 disruption (Fig. 5b and c). We checked the functional consequences of IE8 disruption on relevant features of cancer cells. We did not observe differences either in cell proliferation (Figure S9a) or in clonogenic ability (Figure S9b). Migration capacity assessed by wound healing assay showed no changes between wild-type and IE8 disruption conditions (Figure S9c). However, the ability to degrade collagen-type I – which is a major component of ECM and the breast basement membrane [39] – was decreased after IE8 disruption in MDA-MB-231, whereas we did not observe changes in MDA-MB-436 (Figure S9d). Similarly, anchorage-independent growth assay was significantly reduced in IE8-disrupted MDA-MB-231 cells, but no differences were observed in the IE8-disrupted MDA-MB-436 cells (Fig. 5d). In addition, when we employed collagen I, we observed a significant decrease in the number of invasive cells after IE8 disruption in both the MDA-MB-231 and MDA-MB-436 cells (Fig. 5e).

Fig. 5

Functional impact of changes in MMP protein expression and activity after IE8 disruption. Evaluation of MMP8 (a) and MMP1 (b) levels released to the extracellular space after IE8 disruption in MDA-MB-231 and MDA-MB-436. c. Time-course of MMP-1 activity in MDA-MB-231 and MDA-MB-436 upon IE8 disruption. d. (left) An illustrative example of anchorage-independent spheroid growth of MDA-MB-231 and MDA-MB-436 cell models and (right) volume quantification. e. Colorimetric quantification of cell invasion MDA-MB-231 and MDA-MB-436 cell invasion. Student’s T-test. ns: no significant, *P < 0.05, **P < 0.01, ***P < 0.001. O.D. Optical density

We then evaluated whether the functional consequences changes observed after IE8 are rescued by ectopic expression of MMP1. Thus, MDA-MB-231 and MDA-MB-436 IE8dis cells were transfected with either a functional MMP1 (hMMP1) or an inert catalytic mutant form of MMP1 (MMP1mut, Figure S10a). Ectopic expression of these vectors increased the MMP1 mRNA levels in both cell lines (Figure S10b). As expected, the MMP1 activity was only increased in the cells expressing hMMP1 (Figure S10c). Once the models were validated, we performed functional experiments including cell proliferation, colony formation, migration, and invasion. Ectopic expression of MMP of any of the two variants did not affect cell proliferation (Figure S10d). However, we determined that while hMMP1 increased in the clonogenic ability (Figure S10e) and wound healing rate (Figure S10f), the MMP1mut showed no differences with the IE8-disrupted cells. Beyond that, functional MMP1 overexpression on the IE8-disrupted clones triggered an enhanced invasion rate on collagen I-covered membranes (Figure S10g). Altogether, these results suggest that IE8 disruption diminishes invasiveness potential in the presence of collagen type I fibers by orchestrating a reprogramming of the MMP gene expression pattern.

MMP shift is associated with progression to invasion in early Breast cancer

Finally, we explored whether the ratios between pro-invasive and antitumorigenic MMPs are linked to the progression of ductal carcinomas in situ (DCIS) to invasive ductal carcinoma (IDC) regardless of the breast cancer subtype using the transcriptomic and clinicopathological data from the SCAN-B cohort [31]. We compared clinical specimens of normal breast tissue, pure DCIS, and IDC from different TNM stages. The MMP (3’MMP/5’MMP) ratio was significantly higher in DCIS and IDC when compared to normal breast tissue (PDCIS < 0.001, PIDC < 0.001) (Figure S11a). Importantly, DCIS that eventually progressed to invasive disease displayed a significantly higher MMP ratio than those that did not progress (Fig. 6a). Remarkably when considering the ratio between the pro-invasive MMP1 and the antimetastatic MMP8, the differences in DCIS that progressed to invasive disease were even more pronounced (P = 0.002, Fig. 6b). There were no changes between IDC stratified by tumor size high MMP ratios were also associated with shorter relapse-free survival in patients with DCIS (p-value < 0.001) considering both 3’MMPs/5’MMPs and MMP1/MMP8 ratios (Fig. 6c and d). Importantly, both MMP ratios significantly predicted which DCIS patients will progress to invasive disease (AUC 3’MMPs/5’MMPs ratio = 0.67, AUC of MMP1/MMP8 = 0.77; Fig. 6e and f), but not, but not 3’MMPs and 5’MMPs signatures tested individually (Supplementary figures S11b). Moreover, we found that MMP1 is significantly upregulated in DCIS from patients who had ipsilateral breast events in an independent cohort of DCIS patients studied at the TBCRC 038 clinical trial [33] (Figure S11c).

Fig. 6

Implications of MMPs in DCIS progression. (a) 3’MMPs/5’MMPs and (b) MMP1/MMP8 ratios of standardized gene expression DCIS and invasive ductal carcinoma. Mann Whitney test. ns: no significant, **P < 0.01. Kaplan-Meier curves of relapse-free survival of DCIS using (c) 3’MMPs/5’MMPs and (d) MMP1/MMP8 ratios. ROC curves displaying the performance of (e) 3’MMPs/5’MMPs and (f) MMP1/MMP8 ratios in in situ ductal carcinoma