Tumor-Promoting Crosstalk of MMP11⁺ Fibroblasts and SPP1⁺ Macrophages Drive Poor Prognosis in Breast Cancer: Integrated Multi-Omics Analysis

Background

Breast cancer (BRCA) remains one of the most prevalent malignancies and the leading cause of cancer-related mortality worldwide.¹ Despite unprecedented breakthroughs have been achieved over the past decade, persistent challenges, including immune evasion and chemoresistance, continue to pose formidable clinical obstacles.^2,3 Tumor microenvironment (TME) heterogeneity has been recognized as a critical contributor to therapeutic failure.⁴ Therefore, a deeper understanding of the BRCA TME may unveil novel opportunities for precision oncology.

TME is a dynamically evolving ecosystem that orchestrates complex crosstalk between the tumor, immune, and stromal cells.⁵ Extensive research has established that the TME plays a pivotal role in driving tumor progression, modulating therapeutic response, and fostering immunosuppression landscapes.⁴ As dominant stromal constituents, cancer-associated fibroblasts (CAFs) can promote tumor progression by including angiogenesis, facilitating epithelial-mesenchymal transition (EMT),⁶ establishing immunosuppression niches,⁷ and mediating metabolic reprogramming.⁸ Tumor-associated macrophages (TAMs), key immune components with the TME, can orchestrate immune escape and malignant progression through the secretion of immunosuppressive cytokines and dynamic stromal-immune crosstalk.^9,10 Recent advances in single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) now enable precise dissection of the TME and its intricate interaction networks, thereby providing new insights for mechanistic exploration and therapeutic targeting.

CAFs and TAMs, as two pivotal cellular populations in the TME, engage in extensive bidirectional crosstalk. CAFs secrete various chemokines, such as CCL2 and TGF-β, which mediate monocyte recruitment and drive their differentiation toward an M2-polarized phenotype.¹¹ Polarized TAMs subsequently accumulate at the tumor periphery and actively exclude cytotoxic T lymphocytes (CTLs) through extracellular matrix remodeling, thereby fostering an immune-privileged niche and conferring chemoresistance.¹² Reciprocally, TAMs can activate CAFs via the paracrine secretion of IL-6 and SDF-1, establishing a feedforward loop that promotes EMT-mediated tumor progression and metastatic dissemination.¹³ However, the mechanisms of CAFs-TAMs crosstalk remain elusive, largely due to TME complexity and previous technological limitations. Recent advances in scRNA-seq and ST now provide powerful tools to decipher TME heterogeneity and elucidate such cellular interactions.

In this study, we integrated scRNA-seq with ST to investigate heterogeneity and cellular interactions within the tumor microenvironment. Our analysis identified five distinct CAF subpopulations and six macrophage subpopulations. Among these, MMP11⁺ CAFs and SPP1⁺ TAMs were significantly enriched in BRCA tissues and correlated with poor clinical outcomes. Further in-depth analysis revealed that MMP11⁺ CAFs and SPP1⁺ TAMs were localized in adjacent regions. Cell-cell communication analysis indicated the most extensive communication between the CAFs and macrophages. Collectively, our study delineated the complex cellular interplay between MMP11⁺ CAFs and SPP1⁺ TAMs, highlighting their potential as promising therapeutic targets for future BRCA treatments.

Method

Data Sources

ScRNA-seq data for BRCA were obtained from the Gene Expression Omnibus (GEO) database, including GSE176078 (26 samples) and GSE167036 (8 samples). For bulk transcriptome analysis, gene expression profiles and corresponding clinical data were sourced from from TCGA, METABRIC, and GEO databases (GSE20685 and GSE7390). Only primary breast cancer samples with complete follow-up information and follow-up time greater than zero days were included in our analysis. To increase statistical power, the GSE20685 and GSE7390 datasets were merged into a combined cohort, with batch effects were removed using the ComBat algorithm from the sva R package. A detailed summary of all datasets used in this study, including sample sizes, molecular subtype distributions, platforms, and preprocessing tools, was provided in Table S1.

scRNA-Seq Data Processing

ScRNA-seq data were processed using the Seurat (version 4.3.0) R package. GSE176078 served as the training cohort, and GSE167036 served as the validation cohort. Each dataset was processed separately using an identical analytical pipeline.

To ensure data quality, low-quality cells were excluded based on the following criteria: nFeature RNA < 300 or the percentage of mitochondrial reads > 10%. The gene expression matrix was log-normalized and scaled using the default parameters. Subsequently, high variable gene were identified using the FindVariableFeature function. Principal component analysis (PCA) was performed for dimensionality reduction. To correct for technical variations across individual samples, batch effects were removed using the RunHarmony function from the harmony package (version 1.2.3), with group.by.vars set to “orig.ident”. Batch correction effectiveness was evaluated via UMAP visualization (Figure S1). Cell clusters were identified using FindClusters function with a resolution parameter of 0.2. The resulting clusters were visualized in two-dimensional space using a Uniform Manifold Approximation and Projection (UMAP). Finally, major cell lineages were annotated based on the expression of canonical marker genes.

Identification of Subtypes

To further resolve cellular heterogeneity of CAFs and myeloid cells, we extracted and re-clustered these populations. The cluster algorithm was employed to determine the optimal clustering resolution for each lineage (0.2 for CAFs and 0.35 for myeloid cells). Differentially expressed genes in each subpopulation were identified using the FindAllMarkers function. Subsequently, these subclusters were annotated based on their top differentially expressed genes.

Functional Enrichment Analysis

To assess pathway activities, we performed enrichment analysis was performed using irGSEA R package (version 3.3.2). Enrichment scores for distinct cell subpopulations were calculated using four methods, including AUCell, UCell, single-score, and ssgsea functions, followed by rank aggregation for statistical evaluation. Additionally, we employed the AddModuleScore function from the Seurat R package to quantify the activity of three critical pro-tumorigenic signatures across subpopulations: (1) M1/M2 macrophage polarization, (2) angiogenesis, and (3) epithelial-mesenchymal transition (EMT) activity. The M1 and M2 macrophage signatures were derived from published studies¹⁴ (Supplementary File 1 and Supplementary File 2). The angiogenesis and EMT signatures were downloaded from GSEA database (https://www.gsea-msigdb.org/gsea/msigdb, HALLMARK_ANGIOGENESIS and HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION).

Pseudotime Trajectory Analysis

To investigate the developmental trajectories and dynamic transitions among cell subpopulations, pseudotime trajectory analysis was performed using the Monocle2 R package (version 2.34.0). First, a pseudo-timing analysis objects was created using the “newCellDataSet” function. Subsequently, the “differentialGeneTest” function was employed to select 2000 genes that were most relevant to the pseudo-timing model. Finally, Dimensionality reduction was subsequently performed using the DDRTree algorithm, and cells were ordered along pseudotime based on their transcriptional similarity. The pseudotime trajectory was visualized using the plot_cell_trajectory function.

Intercellular Communication Analysis

To decode intercellular communication within the tumor microenvironment, we applied the CellChat R package (version 2.1.2). The scRNA-seq data were processed using the following CellChat functions: “identifyOverExpressedGenes”, “identifyOverExpressedInteractions”, “computeCommunProb”, and “aggregateNet” functions to process scRNA-seq. The resulting cell-cell communication networks were visualized using the FacenetVisual heatmap function.

Cell Type Deconvolution of ST Data by cell2location

ST data derived from two BRCA samples (CID4290 and CID4465) were obtained from the Zenodo Data Repository (DOI: 10.5281/zenodo.4739739). To infer cell-type abundances at each spatial spot, we applied cell2location (version 0.1.3), a computational framework that integrates spatial transcriptomics with scRNA-seq data to enable spatially resolved cell-type mapping. The model was trained for 50,00 epochs using the following parameters: N_cells_per_location = 20 and detection_alpha = 200.

Spatial Cells Distance Analysis

To investigate proximity-dependent cell–cell interactions within the spatial transcriptomics data, we used the mistyR R package (v1.16.0), which implements the Multi-view Interpretable Spatio-Temporal (MISTy) framework. This approach analyzes cell–cell interactions based on spatial distances between cells, enabling the identification of local and regional communication patterns within the tissue microenvironment.

Survival Analysis and Immune Landscape Characterization

To evaluate the prognostic significance of cell subpopulations, we performed Kaplan–Meier (KM) survival analysis. Infiltration levels of cell subpopulation were quantified in bulk RNA-seq data using the single-sample gene set enrichment analysis (ssGSEA) method, based on the top 15 marker genes from scRNA-seq data. Patients were stratified into high- and low-infiltration groups using the optimal cutoff value determined by the surv_cutoff function from the survminer R package. Overall survival (OS) differences between differential groups were assessed using KM curves with Log rank tests.

To evaluate the combined prognostic effect of MMP11⁺ CAFs and SPP1⁺ macrophages, patients were further stratified into three groups based on co-infiltration levels: low (MMP11⁺ CAFˡ^oʷ/SPP1⁺ macrophageˡ^oʷ), intermediate (MMP11⁺ CAFʰ^igʰ/SPP1⁺ macrophageˡ^oʷ or MMP11⁺ CAFˡ^oʷ/SPP1⁺ macrophageʰ^igʰ), and high (MMP11⁺ CAFʰ^igʰ/SPP1⁺ macrophageʰ^igʰ), with survival differences assessed by Log rank test for trend, and HR (95% CI) derived from univariate Cox regression models.

Tumor purity, immune and stromal scores were calculated using ESTIMATE. Immunotherapy response was predicted by TIDE, and CD8⁺ T cell proportions were estimated by CIBERSORT. Comparisons across groups were performed using Kruskal–Wallis test, and response rate differences by chi-square test for trend.

Immunofluorescence Staining

FFPE tissue section (4 μm) from BRCA patients were deparaffinized, rehydrated, and subjected to antigen retrieval in EDTA buffer using a pressure cooker. Two staining panels were performed.

For three-color staining, sections were blocked with 5% BSA for 1h at room temperature. Subsequently, sections were incubated overnight at 4°C with mouse anti-human FAP (1:200, Cat#222271, Zenbio) and rabbit anti-human MMP11 (1:100, Cat#ET1611-33, HUABIO). After washing, sections were incubated with Coralite594-conjugated goat anti-rabbit IgG (1:300, Cat#SA00013-4, Proteintech) and Coralite488-conjugated goat anti-mouse IgG (1:300, Cat#SA000013-1, Proteintech) for 1h at room temperature in the dark. Nuclei were counterstained with DAPI (1μg/mL) for 5 min.

For five-color multiplex staining, tyramide signal amplification (TSA, Cat#RK05904P, ABclonal) was employed. Endogenous peroxidase was blocked with 3% H₂O₂. Sections were then blocked with 5% BSA for 1 h at room temperature. Sections underwent sequential rounds of staining with primary antibodies against SPP1 (1:200, Cat#83341-1-RR, Proteintech), CD68 (1:300, Cat#ab955, abcam), MMP11 (1:100, Cat#ET1611-33, HUABIO), and FAP (1:200, Cat#222271, Zenbio), each detected by HRP-conjugated secondary antibody (1:500, Cat#SA00001-2/ Cat#SA00001-1, Proteintech) and TSA with distinct fluorophores. Microwave stripping was applied between rounds. Nuclei were counterstained with DAPI.

All slides were mounted with antifade mounting medium and imaged using a fluorescence microscope.

Statistical Analysis

All analyses were performed using R (version 4.5.1) and Python (version 3.10). The flow diagram of this study is shown in Figure 1. Functional analysis of subtypes was performed using DEseq2 and clusterProfiler R packages. Comparisons between two groups were performed using the Wilcoxon rank-sum test for continuous variables and Fisher’s exact test for categorical variables. For analyses involving multiple comparisons, P-values were adjusted using the Benjamini–Hochberg false discovery rate (FDR) method, with raw P-values provided in Table S2. Statistical significance was set at P < 0.05.

Figure 1 Study design and workflow schematic.

Result

Single-Cell Transcriptomic Atlas and Cell-Cell Interactions in the Breast Cancer Microenvironment

A total of 26 breast cancer single-cell transcriptome profiles were included in this analysis. After applying quality control thresholds, 83,033 cells were retained for downstream analysis. Unsupervised clustering analysis identified 11 distinct cell clusters. Based on canonical markers, these clusters were annotated as eight major cell types (Figure 2A-2C): T/NK cells (n = 30,936; identified by CD3D, CD3E, GNLY, NKG7), epithelial cells (n = 119,044; identified by EPCAM, KRT19, KRT8, KRT18), cancer-associated fibroblasts (CAFs, n = 6024; identified by DCN, COL3A1), perivascular cells (PVL, n = 4713; identified by RGS5, PDGFRB), endothelial cells (n = 6565; identified by VWF, PLVAP), B cells (n =4031; identified by CD79A, MS4A1), plasma cells (n = 3139; identified by IGHG1, MZB1), and myeloid cells (n = 8397; identified by LYZ, CD68, CD14).

Figure 2 Single-cell atlas of human breast cancer. (A) UMAP visualization of all cells clustered into 8 cell clusters. (B) Violin plot showing the expression of maker genes across the 8 cell types. (C) Expression levels of selected known marker genes of 8 main cell types. (D) Dot plots showing incoming and outgoing interaction strengths among cell types in BRCA TME. (E) Heatmaps revealed cell-cell interaction strengths within BRCA TME.

To investigate cellular interactions within the tumor microenvironment, we performed a cell-cell communication analysis. Our results showed that CAFs were the predominant signal-sending population, whereas myeloid cells exhibited the highest signal-receiving activity (Figure 2D). Notably, extensive crosstalk has been observed between CAFs and myeloid cells in the BRCA microenvironments (Figure 2E), suggesting that CAFs-myeloid cells crosstalk might constitute a key hallmark of BRCA microenvironment.

MMP11⁺ CAFs Correlate with BRCA Progression and Poor Prognosis

Given that our cell-cell communication analysis identified CAFs as the dominant signaling hub in the BRCA microenvironment, we next sought to deconvolute the heterogeneity of this critical population. Re-clustering analysis revealed that CAFs could be classified into five distinct clusters (Figure 3A). These subpopulations were named according to their most significantly and uniquely overexpressed marker genes: MMP11⁺ CAFs, PLA2G2A⁺ CAFs, SRGN⁺ CAFs, CCL19⁺ CAFs, and APOD⁺ CAFs (Figure 3B). The differentially expressed genes (DEGs) for each cluster are presented in Supplementary File 3. Notably, MMP11⁺ CAFs and APOD⁺ CAFs constituted the predominant subpopulations, representing 30.96% and 39.03% of total CAFs, respectively (Figure 3C).

Figure 3 Characteristics of CAF subtypes in BRCA. (A) UMAP visualization of CAF subpopulation clustered into 5 distinct subtypes. (B) Volcano plots delineating differential gene expression patterns across distinct CAF subpopulations. Statistical significance was determined using the Wilcoxon rank-sum test with Benjamini–Hochberg FDR adjustment; Subpopulation names are derived from the most significantly upregulated gene within each cluster. Each subpopulation is characterized by a unique set of marker genes. (C) Pie chart depicting proportional distribution of CAF subpopulations. (D) Heatmap plot showing the terms of the hallmark gene set enrichment pathway for each CAF cluster (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001). (E) Quantification of EMT pathway across each CAF clusters. (F) Quantification of angiogenesis pathway across each CAF clusters. (G) Left panel: the pseudotime of the developmental trajectory inferred by Monocle2. Right panel: CAF transition along the pseudotime trajectory. (H) EMT score transition along the pseudotime trajectory. (I) Angiogenesis transition along the pseudotime trajectory. (J) Box plots showing infiltration scores of distinct CAF subpopulations in normal and tumor tissues from the TCGA-BRCA dataset. Statistical comparisons between normal and tumor tissues for each CAF subset were performed using the Wilcoxon rank-sum test (***P < 0.001). (K and M) Kaplan–Meier survival curves showing overall survival in the TCGA-BRCA cohort (K), METABRIC cohort (L), and GEO cohort (M) stratified by MMP11⁺ CAF infiltration levels. High infiltration was consistently associated with poor prognosis across all cohorts. Survival comparisons were performed using the Log rank test; hazard ratios (HR) with 95% confidence intervals (CIs) are indicated on each plot.

Functional analysis revealed extensive heterogeneity across these fibroblast populations. MMP11⁺ CAFs were specifically enriched for tumor-promoting processes, including tumor-related epithelial-mesenchymal transition (EMT) and angiogenesis (Figure 3D). Quantitative pathway analysis demonstrated that MMP11⁺ CAFs exhibited the highest enrichment scores for both EMT and angiogenesis (Figure 3E-3F). These results suggested that MMP11+ CAFs might be associated with tumor-promoting EMT and angiogenesis. Pseudotime trajectory analysis positioned MMP11⁺ CAFs at the terminal state of the differentiation trajectory, suggesting their critical role in advanced disease progression (Figure 3G). Notably, both EMT and angiogenesis scores exhibited a consistent distribution pattern with MMP11⁺ CAFs, further indicating their importance in the later stages (Figure 3H-3I).

We next applied single-sample gene set enrichment analysis (ssGSEA) to evaluate the association between MMP11⁺ CAFs and clinical outcomes in BRCA. MMP11⁺ CAFs were significantly enriched in tumor tissues compared with adjacent paracancerous tissues (Figure 3J). Moreover, high infiltration of MMP11⁺ CAFs was associated with poor overall survival in the TCGA-BRCA cohort (Figure 3K, p =0.034). These findings were consistently validated across multiple independent cohorts (Figure 3L-3M, METABRIC cohort: p = 0.042, GEO cohort: p = 0.0029), confirming the robust prognostic significance of MMP11⁺ CAFs in BRCA.

Validation of CAF Subpopulation Across Independent scRNA Datasets

To validate the characterization of the CAF subpopulation, we performed previous analyses using an independent scRNA-seq dataset from 8 BRCA tissues. By clustering 41,751 single cells, we identified eight main cell types based on canonical marker gene expression (Figure 4A-4B): T/NK cells (n = 18,008), epithelial cells (n = 11,866), CAFs (n = 1878), PVL (n = 1096), endothelial cells (n = 2322), B cells (n = 2108), plasma cells (n = 2603) and myeloid cells (n = 1870).

Figure 4 Validation of CAF subtypes in an independent BRCA scRNA-seq dataset. (A) UMAP visualization of all cells clustered into 8 major cell clusters in the independent validation dataset. (B) Violin plot showing the expression of maker genes across the 8 cell types. (C) UMAP visualization of CAF subpopulations clustered into 5 distinct subtypes. (D) Volcano plots delineating differentially expressed genes across distinct CAF subpopulations. Statistical significance was determined using the Wilcoxon rank-sum test with Benjamini–Hochberg FDR adjustment. (E) Quantification of EMT pathway scores across each CAF clusters. (F) Quantification of angiogenesis pathway scores across each CAF clusters. (G) Quantification of MMP11⁺ CAF infiltration. Statistical comparison between tumor and adjacent normal tissues was performed using the Wilcoxon rank-sum test (***P < 0.001).

Subsequent re-clustering CAFs identified five distinct subtypes in the validation dataset (Figure 4C-4D): MMP11⁺ CAFs, PLA2G2A⁺ CAFs, SRGN⁺ CAFs, CCL19⁺ CAFs, and APOD⁺ CAFs. Quantitative pathway analysis suggested that MMP11⁺ CAFs consistently exhibited the highest enrichment scores for EMT and angiogenesis scores in this independent cohort (Figure 4E-4F). Moreover, MMP11⁺ CAFs were also significantly enriched in the BRCA tissues (Figure 4G). To further validate the presence of MMP11⁺ CAFs at the protein level, we performed immunofluorescence staining on breast cancer tissue sections. Co-staining for MMP11 and the CAF marker FAP confirmed that MMP11⁺ cells were localized within the FAP⁺ stromal compartment, with significant enrichment in tumor tissues compared to adjacent normal tissues (Figure S2). Collectively, these results demonstrated that MMP11⁺ CAFs represented a distinct fibroblast subpopulation. The enrichment in BRCA tissues and association with poor prognosis suggested a potential role in BRCA pathogenesis.

SPP1⁺ Macrophage Correlate with BRCA Progression and Poor Prognosis

Our cell-cell communication analysis revealed that myeloid cells exhibited the highest signal-receiving activity in the BRCA microenvironment, with extensive crosstalk observed between myeloid cells and CAFs. Given their prominent role as signal receivers, we next performed unsupervised clustering of myeloid cells to delineate their heterogeneity and identify functionally distinct subsets. Re-clustering myeloid cells identified eight distinct subsets (Figure 5A-5B), including one dendritic cell population (DCs, identified by CD1C and CD1E), one monocyte cell population (mono, identified by S100A8 and S100A9), and six macrophage cells (SPP1⁺ macrophages, CCL4⁺ macrophages, CXCL10⁺ macrophages, SEPP1⁺ macrophages, HSPA6⁺ macrophages, and STMN1⁺ macrophages). The differentially expressed genes (DEGs) for each cluster are presented in Supplementary File 4.

Figure 5 Characteristics of myeloid cells subsets in BRCA. (A) UMAP visualization of eight myeloid subsets identified by re-clustering: DCs, monocytes, and six macrophage subsets (SPP1⁺, CCL4⁺, CXCL10⁺, SEPP1⁺, HSPA6⁺, STMN1⁺). (B) Volcano plots showing average log2 fold-change of significantly upregulated (Orange) and downregulated (purple) genes across distinct myeloid cell subpopulations. Statistical significance was determined using the Wilcoxon rank-sum test with Benjamini–Hochberg FDR adjustment. (C) M1 and M2 module scores across macrophage subsets. (D) Module scores for EMT signatures. (E) Module scores for angiogenesis signatures. (F) ssGSEA analysis showing significant enrichment of SPP1⁺ macrophages in tumor tissues compared to adjacent normal tissues (***P < 0.001, Wilcoxon rank-sum test). (G–I) Kaplan–Meier survival curves in TCGA-BRCA (G), METABRIC (H) and GEO (I) cohorts. High SPP1⁺ macrophage infiltration was consistently associated with poor prognosis across all cohorts. Survival comparisons were performed using the Log rank test; hazard ratios (HR) with 95% confidence intervals (CIs) are indicated on each plot.

To investigate the polarization status of macrophage subsets, we examined the expression of canonical M1 and M2 markers. As shown in Figure S3, M1-associated genes (CXCL9, CXCL10) were predominantly expressed in the CXCL10⁺ subset, while M2-associated genes (CD163, MRC1, TGFB1) showed moderate expression across multiple subsets. Notably, no macrophage subset displayed exclusive enrichment of either M1 or M2 markers, reflecting the complexity of TAM polarization in vivo.

Given this mixed expression pattern, we next used module scores to quantitatively assess polarization bias. As shown in Figure 5C, SPP1⁺ macrophages exhibited the highest M2 module score and low M1 module score among all subsets, indicating a relative M2-like polarization bias, despite the mixed expression of individual markers.

Beyond polarization status, we further characterized the functional features of SPP1⁺ macrophages by examining their enrichment for pro-tumorigenic pathways. SPP1⁺ macrophages demonstrated the highest enrichment scores for EMT and angiogenesis among all subsets, supporting their potential pro-tumoral role (Figure 5D- 5E).

We next evaluated the clinical significance of SPP1⁺ macrophages using ssGSEA. SPP1⁺ macrophages were significantly enriched in BRCA tissues compared to adjacent normal tissues (Figure 5F). High infiltration of SPP1⁺ macrophages was associated with poor overall survival in the TCGA-BRCA cohort (Figure 5G, p = 0.025), and this finding was consistently validated across multiple independent cohorts (Figure 5H-5I, METABRIC: p < 0.001; GEO: p = 0.024), suggesting a potential prognostic role for SPP1⁺ macrophages in BRCA.

Collectively, these results indicate that SPP1⁺ macrophages, characterized by functional enrichment for EMT and angiogenesis and relative M2-like bias, may be associated with breast cancer progression.

scRNA-Seq and ST Revealed Interaction Between MMP11⁺ CAFs and SPP1⁺ Macrophages

Recent studies have demonstrated that CAFs and macrophages engaged in reciprocal interactions that synergistically facilitated tumor progression.^15,16 Given our previous identification of MMP11⁺ CAFs and SPP1⁺ macrophages as key functionally distinct subsets within the BRCA microenvironment, we next sought to investigate the potential crosstalk between these two populations using cell–cell communication analysis.

Applying CellChat to scRNA-seq data, we found that MMP11⁺ CAFs were the most prominent signal-sending population, whereas SPP1⁺ macrophages exhibited the highest signal-receiving activity (Figure 6A). Moreover, extensive interactions were observed between MMP11⁺ CAFs and SPP1⁺ macrophages (Figure 6B-6C). Specifically, MMP11⁺ CAFs directly interacted with SPP1+ macrophages, potentially via the collagen signaling pathway involving COL1A1/COL1A2-CD44 signal axis (Figure 6D). Reciprocally, SPP1⁺ macrophages might signaled to MMP11⁺ CAFs through the SPP1 pathway involving the SPP1-ITGAV/ITGB1 axis (Figure 6E).

Figure 6 Cell-cell communication in BRCA TME. (A) Scatter plot showing outgoing versus incoming interaction strength for different cell types. Each point represents a cell type, with color indicating cell identity and size reflecting the number of interactions. (B and C) Heatmaps depicting the number of interactions (B) and interaction strength (C) between CAF subtypes and macrophage subtypes. Side bars indicate total interactions per cell type. (D and E) Signaling pathway-specific communication networks: (D) Collagen signaling pathway from CAF_MMP11 to Macro_SPP1, and (E) SPP1 signaling pathway from Macro_SPP1 to CAF_MMP11. Dot color represents communication probability, and dot size indicates statistical significance. (F and G) Spatial transcriptomic maps showing infiltration levels of Macro_SPP1 and CAF_MMP11 in two representative breast cancer samples (#1 in F, #2 in G). Color intensity represents infiltration levels. (H–J) Spatial niches and cell-cell interactions identified by MISTy analysis: (H) Heatmap of mean importance scores from MISTy analysis, (I and J) Cell–cell interaction networks derived from MISTy-inferred multiview models in two representative samples (#1 in I, #2 in J). Nodes represent cell types, colored by major functional clusters, and edges denote predicted interactions.

To validate the crosstalk between MMP11⁺ CAFs and SPP1 macrophages, we analyzed spatial transcriptomic data using cell2location framework to infer the distribution of cell populations within the BRCA tumor microenvironment (Figure 6F-6G). Subsequently, we applied MISTy analysis to estimate proximity-dependent cell-cell interaction. Spatial distribution mapping revealed significant intra-view colocalization. Notably, the ST results were consistent with our CellChat predictions from scRNA-seq, suggesting that MMP11⁺ CAFs spatially colocalized with SPP1⁺ macrophages (Figure 6H-6J). To further validate this spatial interaction, we performed multiplex immunofluorescence (mIF) staining on three breast cancer tissue specimens. Representative images confirmed close spatial proximity between MMP11⁺ CAFs and SPP1⁺ TAMs (Figure S4), providing protein-level validation of the cellular crosstalk identified by our transcriptomic analyses.

High MMP11⁺ CAFs and SPP1⁺ Macrophage Infiltration Correlated with Poor Prognosis and Reduced Therapeutic Responsiveness

Given the significant crosstalk observed between MMP11⁺ CAFs and SPP1⁺ macrophages in the BRCA microenvironment, we next sought to investigate its clinical relevance. Patients were stratified into three molecular subtypes based on the infiltration levels of these two cell populations: low (MMP11⁺ CAFs^low/SPP1⁺ macrophages^low), intermediate (MMP11⁺ CAFs^high/SPP1⁺ macrophages^low or MMP11⁺ CAFs^low/SPP1⁺ macrophages^high), and high (MMP11⁺ CAFs^high/SPP1⁺ macrophages^high). Kaplan-Meier analysis revealed that patients with dual-high infiltration of MMP11⁺ CAFs and SPP1⁺ macrophages had significantly worse survival outcomes across three independent cohorts (Figure 7A-7C; TCGA-BRCA: p = 0.011; METABRIC: p = 0.0012; GEO: p = 0.0029), suggesting that these two cell types may synergistically promote BRCA progression.

Figure 7 Prognostic significance of MMP11⁺ CAF and SPP1⁺ TAM co-infiltration. (A–C) Kaplan–Meier survival curves of overall survival in breast cancer across three independent cohorts: (A) TCGA-BRCA, (B) METABRIC, and (C) GEO cohort. Patients were categorized into high, intermediate, and low groups combined MMP11⁺ CAF and SPP1⁺ macrophage infiltration levels. High co-infiltration was consistently associated with poor prognosis across all cohorts. Survival comparisons were performed using the Log rank test for trend, with hazard ratios (HR) and 95% confidence intervals (CIs) derived from univariate Cox regression models. (D) Top 20 KEGG pathways enriched in MMP11+ CAFshigh/SPP1+ macrophageshigh groups. P-values were adjusted using the Benjamini–Hochberg FDR method. The dot size represents -log10 (adjusted p-value), and the x-axis shows the number of genes involved in each pathway. (E–G) Comparison of tumor microenvironment features across high, intermediate, and low CAF signature groups: (E) immune score, (F) stromal score, and (G) tumor purity. Statistical comparisons were performed using the Kruskal–Wallis test. (H–K) Immunotherapy response-associated analyses: (H) TIDE score distribution (Kruskal–Wallis test for statistical comparisons), (I) waterfall plot of TIDE scores with predicted responders (red) and non-responders (blue), (J) percentage of predicted responders in each CAF signature group (chi-square test for trend). (K) CIBERSORT-inferred proportion of CD8+ T cells. Statistical significance was determined using Kruskal–Wallis test.

We next explored the functional and immunological characteristics associated with this high-risk subtype. Functional analysis showed enrichment of multiple immune-related pathways in the high-infiltration subtype (Figure 7D), including Th17 cell differentiation, NF-κB signaling pathway, Th1/Th2 cell differentiation, IL17 signaling pathway, PD-L1 expression and PD-1 checkpoint pathway in cancer, as well as antigen processing presentation. These findings point to distinct immune microenvironment characteristics across the different molecular subtypes.

Consistent with this, ESTIMATE algorithm analysis indicated higher immune and stromal scores, accompanied by lower tumor purity, in the high-infiltration subtype (Figure 7E). Furthermore, TIDE analysis revealed substantial heterogeneity in immunotherapy response across these subtypes (Figure 7F). Notably, the high-infiltration subtype showed the lowest predicted response rate (25.7%), whereas the low-infiltration subtype showed the most favorable response profile (74.2%).

Discussion

BRCA represents the most prevalent malignancy among women worldwide, posing significant threats to female health.^1,17 Systematic characterization of TME in BRCA can uncover patient-specific molecular signatures and facilitates precision therapy.^18,19 Therefore, there is an urgent need to employ scRNA-seq and ST technologies to for a comprehensive analysis of the BRCA TME. Here, we employed an integrated multi-omics approach to systematically delineate the BRCA TME and investigate the interplay between CAFs and TAMs.

CAFs, one of the most abundant stromal components, have been shown to promote tumor progression and immune evasion.²⁰ Recent pan-cancer scRNA analysis revealed that CAFs can be classified into eight functionally distinct CAF subtypes.²¹ These CAF subpopulations demonstrate tumor type-dependent heterogeneity in both composition and functional states.²² Our analysis delineated five functionally heterogeneous CAF subpopulations in BRCA, with APOD⁺ and MMP11⁺ CAFs emerging as the two predominant subpopulations. Among these CAF subpopulations, MMP11⁺ CAFs were significantly enriched in BRCA and were associated with worse overall survival. This subpopulation demonstrated a pronounced myofibroblastic CAF (myCAF) signature, including overexpression of fibrillar collagens (COL1A1 and COL11A1) and enhanced fibronectin production (FN1). MyCAFs are associated with multiple protumorigenic processes.²³ Furthermore, MMP11⁺ CAFs exhibited the highest scores for EMT and angiogenesis compared to the other CAF subpopulations. This finding aligns with that of previous studies. CAF-derived MMP11 promotes EMT in pancreatic ductal adenocarcinoma via the PI3K/AKT pathway.²⁴ Wu et al demonstrated that MMP11+ myCAFs could enhance the proangiogenic activity of ESM1+ tumor endothelial cells (tECs) through CCL11/CCL2 signaling.²⁵ These findings suggested that MMP11+ CAFs are critical mediators of BRCA progression.

TAMs demonstrated remarkable phenotypic flexibility with the TME and was capable of dynamic interconversion between tumoricidal M1-like and pro-tumoral M2-like activation states.^26–28 Studies have revealed that TAMs orchestrate immunosuppression through multiple mechanisms, including PD-L1 expression,²⁹ Treg recruitment,^30,31 and metabolic reprogramming.^32,33 However, accumulating evidence suggests that TAMs in human tumors do not strictly adhere to the conventional M1/M2 binary classification but rather display mixed polarization states.^34,35 Consistently, our analysis revealed moderate expression of M2-associated markers across multiple macrophage subsets, with no subset showing exclusive enrichment of M1 or M2 markers, highlighting the complexity of TAM polarization in vivo. Among these subsets, SPP1⁺ macrophages exhibited a relative M2-like bias, were highly enriched for pro-tumorigenic pathways including epithelial-mesenchymal transition (EMT) and angiogenesis, and were significantly associated with poor prognosis in BRCA. The pro-tumoral role of SPP1⁺ macrophages have been increasingly recognized across multiple cancer types, such as prostate,³⁶ lung,³⁷ and pancreatic cancers,³⁸ where their presence correlates with aggressive tumor features and adverse clinical outcomes. Mechanistically, SPP1⁺ macrophages promote tumor progression through multiple pathways, including angiogenesis, extracellular matrix remodeling, and the establishment of an immunosuppressive niche.^39,40 Our findings extend these observations to breast cancer, revealing a potential link between SPP1⁺ macrophages and tumor immunity. Specifically, high infiltration of SPP1⁺ macrophages was associated with reduced CD8⁺ T cell infiltration and impaired response to immunotherapy, consistent with prior reports demonstrating their negative correlation with cytotoxic T cell presence and contribution to immune checkpoint blockade resistance.^34,41 Collectively, our study suggests that SPP1⁺ macrophages are characterized by a relative M2-like bias and pro-tumorigenic features, potentially contributing to tumor progression and immunotherapy resistance. Further protein-level validation of their polarization state and functional role in T cell exclusion is warranted.

Extensive crosstalk between CAFs and TAM within TME has been established.^42–44 For instance, FAP⁺ CAFs can interact with SPP1+ TAMs to regulate fibrotic structure formation and promote immune exclusion.⁴⁵ Xu et al identified GREM1+ fibroblasts and SPP1+ macrophages in gastric cancer through signaling pathways, including the COLLAGEN and SPP1 signaling pathways.⁴⁶ In our study, we characterized a unique interaction axis between MMP11⁺ CAFs and SPP1⁺ macrophages in the breast cancer TME. Ligand-receptor analysis revealed bidirectional signaling between these two populations: MMP11⁺ CAFs signal to SPP1⁺ macrophages via the COL1A1/COL1A2-CD44 axis, while SPP1⁺ macrophages respond through the SPP1-ITGAV/ITGB1 signaling pathway. A recent study in bladder cancer reported that MMP11⁺ CAFs recruit SPP1⁺ TAMs via the CCL11/CCL2 chemokine axis,²⁵ suggesting a conserved pro-tumoral role for this CAF subset across cancer types. Distinct from this chemokine-mediated recruitment, the collagen- and SPP1-mediated ligand-receptor interactions identified in our study represent a unique feature in breast cancer. As summarized in Table S3, our findings reveal a previously underappreciated interaction module.

Drug repurposing has emerged as a cost-effective strategy for developing novel anticancer therapies, offering advantages such as reduced development timelines, lower costs, and established safety profiles of approved drugs.⁴⁷ Given the potential involvement of the MMP11⁺ CAF–SPP1⁺ macrophage axis in pro-tumorigenic inflammation and immune exclusion, targeting this crosstalk through drug repurposing may represent a promising therapeutic strategy in breast cancer. Metabolic agents such as metformin and statins can reprogram TAMs,^47,48 thereby potentially reducing the pro-tumorigenic activity of SPP1⁺ macrophages, while thalidomide analogs may further suppress SPP1⁺ TAM function by modulating the stromal cytokine milieu.⁴⁹ Concurrently, targeting the structural barriers orchestrated by MMP11⁺ CAFs is critical for enhancing therapeutic efficacy; anti-fibrotic agents such as losartan can decompress the tumor microenvironment (TME) by reducing collagen deposition and CAF activation,^50,51 potentially downregulating MMP11 expression and facilitating T-cell infiltration, whereas tranilast has been shown to inhibit TGF-β-induced fibroblast activation,⁵² thereby potentially disrupting the MMP11⁺/SPP1⁺ interaction circuit. Collectively, integrating these repurposed immunomodulatory agents with standard-of-care therapies may offer a multifaceted strategy to “heat up” the TME, converting immune-excluded tumors into immune-inflamed phenotypes. Future experimental validation is needed.

Precision medical approaches enable more effective therapeutic interventions for breast cancer.⁵³ In our analyses, we identified three molecularly distinct breast cancer subtypes based on MMP11⁺ CAFs and SPP1⁺ macrophages infiltration patterns. Patients with MMP11⁺ CAFs^high and SPP1⁺ macrophages^high displayed the most unfavorable clinical outcomes, a distinctive tumor microenvironment characterized by prominent immune/stromal infiltration and poor immunotherapy response. Although the MMP11⁺ CAFs^high/SPP1⁺ macrophage^high subtype displayed an “immuno-hot” signature, the co-enrichment of these stromal and immune components may create a dysfunctional tumor ecosystem that favors immune evasion. Recent studies have demonstrated that MMP11⁺ CAFs contribute to immunotherapy resistance through extracellular matrix remodeling, which might physically exclude cytotoxic T cells, whereas SPP1+ macrophages likely promote an immunosuppressive niche via multiple pathways.⁵⁴ Our findings provide a novel framework for personalized therapy.

Although our findings provide novel insights, several limitations should be considered. First, the absence of matched multi-omic samples (bulk, scRNA, spatial) limits our ability to fully resolve spatial architecture. Second, the MMP11⁺ CAFs-SPP1⁺ macrophages crosstalk lacks functional validation, despite preliminary mIF support. Third, reliance on public datasets introduces heterogeneity in sample annotation. Fourth, defining cell phenotypes from transcriptomic data remains ambiguous, as macrophage polarization and CAF subtyping represent continuous spectra. Fifth, we did not assess GZMB⁺ and IFN-γ⁺ T cell subsets, which would be valuable for understanding the impact of this niche on anti-tumor immunity. Moreover, the predominance of luminal tumors may skew pathway analyses toward hormone receptor-associated signatures. Thus, validation in larger, balanced cohorts is warranted. Future directions should include: (1) functional validation using co-culture systems or organoid models; (2) protein-level characterization of T cell subsets (GZMB⁺ and IFN-γ⁺) via mIHC; (3) cytokine profiling pre- and post-drug treatment.

Conclusion

In this study, we conducted an integrative multi-omics analysis combining bulk RNA sequencing, single-cell RNA sequencing and spatial transcriptomics to characterize the complex landscape of the breast cancer tumor microenvironment. Notably, both MMP11⁺ CAFs and SPP1⁺ macrophages were significantly associated with poor prognosis. Moreover, we uncovered a sophisticated crosstalk between MMP11⁺ CAFs and SPP1⁺ macrophages, which may serve as a promising therapeutic target for future breast cancer treatment strategies.