Abstract

   Breast cancer categorized into hormone receptor-positive (HR+),
   HER2-positive (HER2+), and triple-negative (TNBC) subtypes, exhibits
   varied outcomes based on the number of tumor-infiltrating lymphocytes
   (TILs). To explore the divergent roles of TIL levels across different
   subtypes, we employed single-cell RNA sequencing on 31 patients with
   breast cancer. HR+ breast cancer with high TIL levels (TIL-high)
   revealed increased SPP1+ macrophages, increased SPP1 expression in
   other monocytes/macrophages (mono/macro) subgroups, and enriched
   pathways associated with extracellular matrix (ECM) remodeling in
   mono/macro. Moreover, cell–cell interaction analyses revealed enhanced
   SPP1, MIF, and FN1 signaling in the interaction between SPP1+
   macrophages and T-cells in TIL-high HR+ breast cancer. Spatial
   transcriptomics data highlighted the close proximity of SPP1+
   macrophages, CD8+ T-cells, and CD4+ T-cells in TIL-high HR+ breast
   cancer. Our findings unveil the novel influence of SPP1+ macrophages on
   T-cells in TIL-high HR+ breast cancer, potentially explaining the poor
   prognosis and offering insights for targeted interventions.

   Subject terms: Breast cancer, Cancer microenvironment, Tumour
   heterogeneity

Introduction

   Breast cancer is the second most prevalent type of cancer globally,
   with an age-standardized incidence rate of 46.8 per 100,000 people in
   2022^[42]1. It can be classified into three subtypes based on the
   expression of hormone and HER2 receptors: hormone receptor-positive
   (HR+), HER2-positive (HER2+), and triple-negative breast cancer
   (TNBC)^[43]2. A critical prognostic factor in breast cancer is the
   level of tumor-infiltrating lymphocytes (TILs), which can be assessed
   through histologic examination^[44]3. TILs are defined as mononuclear
   cells, (including plasma cells) that infiltrate the stromal tissue
   within the invasive tumor area^[45]4. TILs exhibit subtype-specific
   characteristics, wherein HER2+ breast cancer and TNBC typically
   manifest higher TIL levels than those of HR+ breast cancer^[46]5. In
   HER2+ breast cancer and TNBC, increased TIL levels are associated with
   improved pathological complete response (pCR) rates after neoadjuvant
   chemotherapy (NAC), longer disease-free survival, and overall
   survival^[47]6. Conversely, in HR+ breast cancer, higher TIL levels are
   linked to shorter overall survival^[48]7,[49]8.

   Numerous immune cell types can modulate TIL function within the tumor
   microenvironment (TME)^[50]9. Among these, tumor-associated macrophages
   (TAMs) represent a significant proportion of breast cancer immune cells
   and are closely related to unfavorable prognoses^[51]10,[52]11.
   Single-cell RNA sequencing (scRNA-seq) studies have revealed specific
   TAM subsets linked to T-cell infiltration and immunosuppression in
   patients with breast cancer^[53]12,[54]13.

   One of the specific TAM subsets, SPP1+ macrophages, plays a crucial
   role in the TME of various cancers. For example, in Colorectal cancer
   (CRC), the presence of SPP1+ macrophages is closely correlated with
   tumor-specific FAP+ fibroblasts, resulting in a desmoplastic fibrosis
   that limits T-cell infiltration and reduces the efficacy of
   immunotherapy^[55]14. Moreover, in head and neck squamous cell
   carcinoma, non-small lung cancer, and CRC, macrophage polarity defined
   by the CXCL9 and SPP1 expression ratio is strongly associated with the
   disease outcome, serving as a relevant prognostic indicator^[56]15.
   These findings underscore the significance of understanding the role
   played by various TAMs in shaping the TME and may be indicative of
   potential therapeutic strategies targeting the interaction between TAMs
   and other cellular components to enhance treatment efficacy in various
   cancers. However, our knowledge of differences in immune cell
   composition, including TAMs, according to TIL levels is limited, and
   the influence of such differences on immune cell infiltration remains
   poorly understood.

   Different subtypes of breast cancers have different TMEs^[57]16. HR+
   breast cancer often exhibits not only a low infiltration of TILs but
   also higher number of TAMs, which are typically immune-suppressive and
   associated with poor survival^[58]17. The emergence of novel
   technologies, including single cell and spatial transcriptomics, has
   enabled more detailed characterization of the differences in TME among
   breast cancer subtypes.

   Here, we stratified patients with breast cancer based on their TIL
   levels (‘TIL-high’ and ‘TIL-low’), considering each subtype (HR+,
   HER2+, and TNBC). We then analyzed the differences in the composition
   of CD45+ immune cell subtypes using scRNA-seq. Furthermore, we
   elucidated the differences in immune cell interactions within the TMEs
   between TIL-high and TIL-low breast cancers by employing cell–cell
   interaction analysis and spatial transcriptomics.

Results

Heterogeneity of tumor microenvironment in breast cancer

   To identify the characteristics of TME in breast cancer and its
   differences based on TIL levels, we conducted scRNA-seq using isolated
   CD45+ cells from 21 primary tumor tissues with accompanying TIL data.
   These TIL levels were evaluated by a certified pathologist using H&E
   stained slides. We also included 10 samples from the Genome Sequence
   Archive (GSA, HRA000477) dataset^[59]18. In total, we analyzed 31
   samples (21 samples from this study and 10 from the GSA database) using
   scRNA-seq, comprising 16 TNBC (10 and 6), 11 HR+ (8 and 3), and 4 HER2+
   samples (3 and 1). Among the HR+ breast cancer cases, seven samples
   were classified as luminal B type and four as luminal A type. Samples
   with TIL values ≥10% were classified as the TIL-high group (n = 18,
   (TNBC n = 10, HR+ n = 5, and HER2+ n = 3)) and those with <10% as the
   TIL-low group (n = 13, (TNBC n = 6, HR+ n = 6, and HER2+ n = 1))
   (Supplementary Table [60]1). Quality control was conducted to remove
   low-quality cells, after which 110,059 cells were left for analysis
   (Supplementary Fig. [61]1a). After removing the batch effect, cells
   were distributed evenly unaffected by sample type, breast cancer
   subtype, or TIL levels (Fig. [62]1a, b). To annotate CD45+ cells, we
   used the PBMC reference in Azimuth, a reference-based mapping
   algorithm^[63]19,[64]20. Based on the celltype.l2 annotated by Azimuth
   (Supplementary Fig. [65]1b, c), we classified immune cells into four
   categories: T-cells (T-cell, 62,489 cells), B cells (B-cell, 14,097
   cells), dendritic cells (DC, 6637 cells), and monocytes/macrophages
   (mono/macro, 26,836 cells) (Figs. [66]1c, [67]S1d). The annotated major
   cell types were validated using canonical markers (Fig. [68]1d). In
   addition to celltype.l2, we conducted a more detailed annotation of
   immune cell subgroups, notably T-cells and mono/macro, through
   re-clustering these cells extracted from the entire cell population
   (Figs. [69]1e and [70]S2a–e). The expression pattern of marker
   genes^[71]20 in detailed cell types is presented in Fig. [72]1f. To
   investigate a difference in the proportion of cell types between
   samples, we generated plots depicting the proportion of cell types
   relative to TIL levels and subtypes (Fig. [73]1g, h). Regardless of TIL
   levels and subtypes, breast cancer showed a high degree of
   heterogeneity in the proportion of immune cell subtypes (Fig. [74]1h).
   The proportions of major immune cell subtypes did not significantly
   differ among HR+, HER2+, and TNBC (Fig. [75]1i).

Fig. 1. Heterogeneity of immune cell composition in breast cancer.

   [76]Fig. 1
   [77]Open in a new tab

   a–c UMAP plot shows breast cancer subtypes, TIL levels, cell types, and
   detailed cell types in 31 samples. d Dot plot shows each major cell
   type and their expression of selected marker genes. e UMAP plot shows
   detailed cell types in 31 samples. f Heatmap depicts the expression of
   marker genes for detailed cell types, arranged according to major cell
   types. g Bar plot shows TIL value. TIL levels of each sample are
   represented in colors. h Bar plot shows the proportion of cell types in
   31 samples. i Bar plot shows the proportion of cell types in breast
   cancer subtypes. The proportion of cell types in each subtype was
   determined using mean values. TNBC Triple-negative breast cancer, HR+
   Hormone receptor-positive breast cancer, HER2+ HER2-positive breast
   cancer, TIL tumor-infiltrating lymphocytes, DC dendritic cells,
   Mono/macro monocytes, and macrophages.

Increased proportion of SPP1+ macrophages in the TIL-high group of HR+ breast
cancer

   To identify variation in immune cell composition according to TIL
   levels, we examined immune cell composition in more detail. Considering
   the diverse biological characteristics exhibited by different breast
   cancer subtypes, compositions within the same subtype were compared. In
   the immune cell distribution based on TIL levels, each subtype
   exhibited different patterns (Figs. [78]2a–c, [79]S3a–c). We initially
   compared major cell type compositions across TIL levels. Major cell
   type proportions showed no significant differences in each subtype
   (Figs. [80]2d, [81]S3d–f). Then, we divided cells by major cell
   branches (T, B, and dendritic cells and mono/macro) and compared the
   proportions of cells within each branch, further dissected into
   detailed subgroups, based on TIL levels (Supplementary Fig. [82]3g–i).
   Two cell subgroups exhibited statistically significant differences in
   HR+ breast cancer according to TIL levels. In the TIL-low group, there
   was a notable enrichment of mucosal-associated invariant T-cells (MAIT)
   in T-cell (p = 0.028) (Supplementary Fig. [83]3k). However, due to
   insufficient cell counts for meaningful comparative analysis (87 in the
   TIL-low group and 1 in the TIL-high group), our investigative focus
   shifted towards the TIL-high group. Within this group, a remarkable
   increase in the proportion of SPP1+ macrophages in mono/macro was
   observed (p = 0.0043) (595 in the TIL-high group and 54 in the TIL-low
   group) (Figs. [84]2e, f, [85]S3j). Sub-clustering of HR+ breast cancer
   mono/macro (a total of 3407 cells in the TIL-high group and 1374 cells
   in the TIL-low group) highlighted differences in the distribution of
   SPP1+ macrophages based on the TIL levels (Fig. [86]2g).

Fig. 2. Comparison of cell type proportion according to TIL levels revealed
increased SPP1+ macrophages in TIL-high HR+ breast cancer.

   [87]Fig. 2
   [88]Open in a new tab

   a–c UMAP plot shows detailed immune cell subgroups according to TIL
   levels and breast cancer subtypes. d Bar plot shows the cell type
   proportion of the breast cancer subtype in accordance with TIL levels
   (left), and the corresponding sample size (right). The proportion of
   cell types in each group was determined using mean values. e Bar plot
   shows the proportion of subgroups of mono/macro in HR+ breast cancer in
   accordance with TIL levels. The proportion of cell types in each group
   was determined using mean values. f Boxplot shows the proportion of
   SPP1+macro in accordance with TIL levels (Wilcoxon rank sum test). The
   box boundaries represent the top and bottom quartiles, encompassing the
   interquartile span. The whiskers stretch to 1.5 times the interquartile
   span, with any data points beyond them being labeled as outliers. g
   UMAP plot shows the proportion of mono/macro in HR+ breast cancer in
   accordance with TIL levels. Red circles indicate the differences in the
   number of SPP1+ macro. TNBC Triple-negative breast cancer, HR+ Hormone
   receptor-positive breast cancer, HER2+ HER2-positive breast cancer, TIL
   tumor-infiltrating lymphocytes, DC dendritic cells, Mono/macro
   monocytes, and macrophages.

SPP1+ macrophages and mono/macro in TIL-high group enriched in ECM
remodeling-associated pathways in HR+ breast cancer

   To explore the characteristics of SPP1+ macrophages in HR+ breast
   cancer, we performed pathway enrichment analysis with DEGs of SPP1+
   macrophages with other mono/macro. From the analysis of DEGs, we
   identified 24 genes with up-regulated expression with an average log[2]
   fold change greater than 1 and a p-value less than 0.05, including
   MMP12, SPP1, and NURP1. These genes were subsequently selected for
   pathway enrichment analysis (Fig. [89]3a). In this analysis, pathways
   related to ECM remodeling, such as ECM degradation, collagen
   degradation, and ECM organization, were significantly enriched in SPP1+
   macrophages from HR+ breast cancer (Fig. [90]3b). In other subgroups of
   mono/macro, the CTSK+ macrophages subgroup tended to increase in the
   TIL-high group without statistical significance, enriched in similar
   pathways related to ECM remodeling (Supplementary Fig. [91]4a).

Fig. 3. Pathways related to ECM remodeling process and SPP1 are increased in
TIL-high HR+ breast cancer mono/macro.

   [92]Fig. 3
   [93]Open in a new tab

   a Volcano plot of DEGs between SPP1+ macro and other mono/macro
   subgroups of HR+ breast cancer. Red dots indicate DEGs. b Bar plot
   shows enriched pathways of genes upregulated in SPP1+ macro. c Volcano
   plot of DEGs between HR+ mono/macro of TIL-high and TIL-low group. Red
   dots indicate DEGs. d Bar plot shows enriched pathways of genes
   upregulated in TIL-high HR+ mono/macro. e Feature plot shows level of
   SPP1 expression in HR+ Mono/macro in accordance with the TIL levels. f
   Violin plot shows level of SPP1 expression in HR+ Mono/macro in
   accordance with the TIL levels. SPP1+ macro SPP1+ macrophages, DEGs
   differentially expressed genes, macro macrophages, mono/macro
   monocytes, and macrophages, TIL tumor-infiltrating lymphocyte, NS not
   significant, avg_Log2FC average log[2] fold change, p_val p-value,
   p.adjust adjusted p-value.

   Due to our interest in differences based on TIL levels, we also
   conducted an analysis comparing the TIL-high and -low groups within
   SPP1+ macrophages and other mono/macro in HR+ breast cancer.
   Especially, in 8 out 10 mono/macro cell subgroups, including SPP1+
   macrophages, the ECM degradation pathway was enriched in the TIL-high
   group (Supplementary Fig. [94]4b). To establish whether pathway
   enrichment related to ECM remodeling was a characteristic of mono/macro
   from TIL-high HR+ breast cancer, we used 32 genes upregulated in
   mono/macro from TIL-high HR+ breast cancer compared to TIL-low for
   pathway enrichment analysis (Fig. [95]3c). Indeed, pathways related to
   ECM remodeling including ECM degradation and organization, were
   enriched in the TIL-high group (Fig. [96]3d).

   In the analysis of DEGs within HR+ mono/macro, SPP1 gene showed
   up-regulated expression in the TIL-high group (adjusted
   p-value = 6.128E-273, log[2] fold change = 2.963) (Fig. [97]3c). Both
   feature and violin plots confirmed higher SPP1 expression in the
   TIL-high group compared to the TIL-low group, and this was not specific
   to the SPP1+ macrophages subgroup (Fig. [98]3e, f). These trends were
   not observed in other breast cancer subtypes (Supplementary Fig.
   [99]4c, d).

   To examine whether changes in TAM and TIL subsets are associated with
   ECM rearrangement, Masson’s trichrome staining was performed on the HR+
   breast cancer samples. An increase in collagen fiber bundles in the
   intratumoral stroma area of TIL-high HR+ breast cancer was observed
   compared to that in TIL-low HR+ breast cancer (Supplementary Fig.
   [100]5a–h).

   In addition, we compared the expression of genes involved in
   angiogenesis in SPP1+ macrophages and other mono/macro to investigate
   whether SPP1+ macrophages played an angiogenic role, as previously
   reported^[101]21. However, we did not observe any significant increase
   in the expression of angiogenic genes within SPP1+ macrophages
   (Supplementary Fig. [102]6a–c). Similarly, no positive correlation of
   SPP1 expression and angiogenic gene expression was observed in
   mono/macro (Supplementary Fig. [103]6d–f). Taken together, TIL-high HR+
   breast cancer mono/macro exhibited higher expression of SPP1,
   particularly in SPP1+ macrophages subgroup, and displayed gene
   enrichment associated with ECM remodeling, leading to an increase in
   collagen deposition in TME.

SPP1+ macrophages in TIL-high HR+ breast cancer interact more with T-cells
via SPP1, MIF, and FN1 signaling

   Cell–cell interaction analyses were performed to study the interactions
   of SPP1+ macrophages with other immune cells using CellChat^[104]22. We
   initially compared TIL-high and low HR+ breast cancer with differential
   number of interactions and interaction strength (Fig. [105]4a, b). The
   TIL-high group displayed a higher number of interactions within
   mono/macro, CD4T, and CD8T, as well as enriched interaction numbers
   among mono/macro, CD4T, and CD8T (Fig. [106]4a). Notably, the
   interaction strength, especially from mono/macro to CD8T and CD4T,
   exhibited greater intensity in the TIL-high group (Fig. [107]4b).
   Consequently, we proceeded to analyze which signaling interactions were
   involved in the direction from mono/macro to CD8T and CD4T. When
   analyzed for each cell subgroup, the interaction number and strength
   were most pronounced in TIL-high when CD8 naive and CD8 proliferating
   were the target cells (Supplementary Fig. [108]7a). Comprehensive
   analysis of mono/macro towards CD8T and CD4T within the HR+ subtype
   identified significant enrichment of 19 signaling pathways in the
   TIL-high group (Fig. [109]4c). When analyzing the interaction in the
   directions of CD8T and CD4T separately, the signaling pathways enriched
   in TIL-high did not differ substantially from the previous results (19
   per 19 pathways matched in CD8T, and 12 per 19 pathways matched in
   CD4T) (Supplementary Fig. [110]7b). Therefore, we conducted a combined
   analysis of interactions from mono/macro to both cell groups. Among the
   19 signaling pathways examined, 7 pathways (RESISTIN, SN, VCAM, PARs,
   IL10, IL16, and PECAM1) showed enrichment exclusively in the TIL-high
   group, while 12 pathways (CD99, SPP1, CD86, CLEC, MIF, FN1, NECTIN,
   CCL, APP, PD-L1, ITGB2, and GALECTIN) were enriched in both groups but
   demonstrated greater significance in the TIL-high group (Fig. [111]4c).
   Subsequently, to investigate which cell subgroup exerts an influence,
   we analyzed each signaling pathway (Figs. [112]4d–f, [113]S7 c, d). The
   analysis of each signaling pathway in the TIL-high and -low groups
   revealed the involvement of multiple cell subgroups. Particularly,
   signaling pathways such as SPP1, MIF, FN1, CD99, CD86, NECTIN, CLEC,
   GELECTIN, and ITGB2 indicated the participation of most mono/macro cell
   subgroups in interactions with CD8T and CD4T. Especially, in the case
   of the SPP1, MIF, and FN1 signaling pathways, SPP1+ macrophages in the
   TIL-high group acted as the major sender cell with the highest
   communication probability compared to other mono/macro subgroups (Fig.
   [114]4d–f). In the analysis of all significant interactions from SPP1+
   macrophages to CD4T and CD8T within the TIL-high group, the following
   three signaling pathways emerged as major contributors (Supplementary
   Fig. [115]7e): SPP1, MIF, and FN1. Moreover, when comparing signaling
   of SPP1+ macrophages towards CD8T and CD4T between TIL-high and -low
   groups of the HR+ subtype, these pathways were more enriched in the
   TIL-high group (Supplementary Fig. [116]7f). These findings highlight
   the substantial involvement of SPP1+ macrophages, especially within the
   three identified signaling pathways, in the interactions between
   mono/macro and CD8T/CD4T within the TIL-high context. This implies a
   crucial role for SPP1+ macrophages in mono/macro and T-cell
   interactions in the TIL-high group.

Fig. 4. Cell-to-cell interaction comparison analysis between TIL-high and
TIL-low groups in HR+ breast cancer.

   [117]Fig. 4
   [118]Open in a new tab

   a Differential number of interactions between TIL-high and TIL-low.
   Arrow represents direction, colored edges indicate increased signaling
   in TIL-high (red) or TIL-low (blue) and line width reflects the
   differential number of interactions. b Differential interaction
   strength between TIL-high and TIL-low, Arrow represents direction,
   colored edges indicate increased signaling in TIL-high (red) or TIL-low
   (blue) and line width reflects the differential interaction strength. c
   Bar plot shows all significant signaling pathways from mono/macro to
   CD8T and CD4T that are ranked based on their differences of overall
   information flow between TIL-high and TIL-low groups in HR+ breast
   cancer. The signaling pathways colored by red are more enriched in the
   TIL-high group, the black ones are equally enriched in the TIL-high and
   TIL-low groups, and the blue ones are more enriched in the TIL-low
   group. d Communication probability of SPP1 signaling network of source
   (sender, mono/macro) (rows) across the target (CD8T/CD4T) (columns) in
   TIL-low (left) and TIL-high (right) groups. e Communication probability
   of MIF signaling network of source (sender, mono/macro) (rows) across
   target (CD8T/CD4T) (columns) in TIL-low (left) and TIL-high (right)
   groups. f Communication probability of FN1 signaling network of source
   (sender, mono/macro) (rows) across the target (CD8T/CD4T) (columns) in
   TIL-low (left) and TIL-high (right) groups. TIL tumor-infiltrating
   lymphocyte, mono/macro monocytes, and macrophages, DC dendritic cells
   CD4T CD4+ T-cells, CD8T CD8+ T-cells, NK NK T-cells, OtherT other
   T-cells.

   Additionally, we examined the expression of immunosuppression and
   proliferation-related markers in T-cells based on the expression of
   SPP1 in mono/macro to evaluate the impact of SPP1 signaling on T-cells.
   The analysis revealed no significant differences (Supplementary Fig.
   [119]8a–c). However, CD8T/CD4T cells in TIL-high HR+ breast cancer were
   enriched in interleukin −4, −10, and −13 signaling pathways, which are
   related to the regulation of immune response in T-cells (Supplementary
   Fig. [120]8d, e). Collectively, within the HR+ subtype, mono/macro and
   CD8T/CD4T engage in interactions through various signaling pathways in
   the TIL-high context, particularly through signaling pathways such as
   SPP1, MIF, and FN1 involving SPP1+ macrophages.

Spatial colocalization of SPP1+ macrophages and CD8T/CD4T in TIL-high HR+
breast cancer

   Through cell–cell interaction analysis, we demonstrated differential
   interactions between mono/macro and CD8T/CD4T in TIL-high and -low HR+
   breast cancer. In addition, among the signaling pathways enriched in
   TIL-high, 9 pathways (SN, PECAM, CD99, CD86, CLEC, NECTIN, APP, ITGB2,
   and PD-L1) were related to cell–cell adhesion^[121]22 (Figs. [122]4d–f,
   [123]S7 c, d). Furthermore, we identified that SPP1+ macrophages
   interacted more with CD8T/CD4T through three signaling pathways in
   TIL-high group. Therefore, direct interactions between the two cell
   groups could be inferred.

   To explore the spatial proximity of these cells for potential
   interactions, we used two spatial transcriptomics datasets of HR+
   breast cancer^[124]23. They were classified as TIL-high (CID4535) and
   TIL-low (CID4290) based on H&E images (Figs. [125]5a, c, [126]S9 a, b).
   We confirmed appropriate pathologic annotations were performed to
   determine the distribution of invasive cancer, stroma, lymphocyte, and
   other components in deconvoluted classification. (Fig. [127]5b, d).

Fig. 5. SPP1+ macrophages and CD8T/CD4T are spatially colocalized in TIL-high
HR+ breast cancer.

   [128]Fig. 5
   [129]Open in a new tab

   a H&E image of the TIL-low (CID4290) tissue section analyzed using
   VISIUM. b Classification of the VISIUM public dataset^[130]23. c H&E
   image of the TIL-high (CID4535) tissue section analyzed using VISIUM. d
   Classification of the VISIUM public dataset. e Signature score of
   mono/macro, SPP1+macro, CD8T, and CD4T in TIL-low (CID4290). f
   Signature score of mono/macro, SPP1+macro, CD8T, and CD4T in TIL-high
   (CID4535). g Correlation of signature scores between SPP1+macro and
   CD8T (left), and between SPP1+macro and CD4T (right) presented as
   Spearman’s correlation coefficient. SPP1+macro SPP1+ macrophages, CD8T
   CD8+ T-cells, CD4T CD4+ T-cells, TIL tumor-infiltrating lymphocytes.

   Subsequently, we employed signature scoring of cell subtypes using the
   top 25 genes to annotate the spots obtained through spatial
   transcriptomics. These scores were then used for conducting correlation
   analyses across the scores in each spot. Signature scores of SPP1+
   macrophages, CD8T, and CD4T were calculated using our scRNA-seq data of
   CD45+ breast cancer (Fig. [131]5e, f). To validate the legitimacy of
   the signature score of our CD45+ sorted dataset, we employed a public
   non-sorted scRNA-seq dataset^[132]23 ([133]GSE176078) to establish
   additional signature score criterion. Comparing the distribution of
   signature scores calculated from the public dataset with our dataset,
   we found no significant differences (Supplementary Fig. [134]10c–f).
   Then, we analyzed the correlation of signature scores in each spot of
   interest. To mitigate the impact of variations in TIL quantity on the
   scores, we normalized each score. After confirming that there were no
   statistically significant differences in the distribution of scores
   between the two samples (Supplementary Fig. [135]10g), we proceeded
   with the analysis. The correlation between SPP1+ macrophage scores and
   CD8T scores (R^2 = 0.44) and between SPP1+ macrophage scores and CD4T
   scores (R^2 = 0.38) was also higher in the TIL-high samples than that
   in the TIL-low samples (R^2 = 0.04 and 0.016, respectively) (Fig.
   [136]5g). Additionally, the correlation between mono/macro scores and
   CD8T scores (R^2 = 0.52) and between mono/macro scores and CD4T scores
   (R^2 = 0.43) was higher in the TIL-high sample compared to the low
   sample (R^2 = 0.045 and 0.017, respectively) (Supplementary Fig.
   [137]10h). These results indicate that in the TIL-high sample, there
   may be a higher probability of colocalization not only between SPP1+
   macro and CD8T/CD4T but also between mono/macro and CD8T/CD4T, than in
   the TIL-low sample.

   During the analysis of TIL-high samples, we aimed to elucidate
   differences in T-cell behavior influenced by the proximity of SPP1+
   macrophages. For this purpose, using the signature scores of T-cells
   (CD4T, CD8T) and SPP1+ macrophages, we classified the spots into high-
   and low-proximity categories and conducted comparative analyses between
   them (Supplementary Fig. [138]11a). Through differentially expressed
   gene (DEG) and pathway enrichment analyses, specifically in high
   proximity spots, we found enrichment not only in ECM
   remodeling-associated pathways in high proximity spots but also in
   pathways such as immunoregulatory interactions between a lymphoid and a
   non-lymphoid cell and signaling through the interleukins (Supplementary
   Fig. [139]11c). Additionally, in the DEG analysis, we observed
   up-regulation in the expression of the FN1 gene associated with FN1
   signaling, which is involved in the interaction between T-cells and
   SPP1+ macrophages, as identified through single-cell analysis
   (Supplementary Fig. [140]11b). These findings suggest a possible role
   of SPP1+ macrophages in inhibiting T-cell response to tumor cells
   through ECM remodeling and immunoregulatory signaling.

   Lastly, to determine if SPP1+ macrophages are associated with
   prognosis, we conducted a survival analysis using bulk RNA sequencing
   data from dataset [141]GSE96058^[142]24, focusing on SPP1+ macrophage
   markers. First, SPP1+ macrophage markers were identified based on their
   absolute log[2] fold change values from our dataset (Supplementary
   Table [143]2). Then, we explored the expression of these markers in the
   breast cancer scRNA-seq dataset^[144]23 ([145]GSE176078) using CZ
   CELLxGENE Discover^[146]25 to select markers that are either
   exclusively expressed or highly expressed in macrophages (Supplementary
   Fig. [147]12a). Among the markers, we used SPP1 and MS4A6A, which are
   highly expressed in mononuclear phagocytes and have the highest
   absolute log[2] fold change (avg_log2FC 2.02 and -0.98). Since bulk RNA
   sequencing includes all cell types, the expression of
   macrophage-specific genes may be diluted. To address this, we utilized
   the ratio of SPP1 (high expression in SPP1+ macrophages) to MS4A6A (low
   expression in SPP1+ macrophages) to find the prognostic effect of SPP1+
   macrophages. The results demonstrated that a high SPP1/MS4A6A ratio
   (indicating higher SPP1 expression relative to MS4A6A) was associated
   with poor overall survival in ER+ breast cancer (HR = 1.44, logrank
   P = 0.012) (Supplementary Fig. [148]12b).

Discussion

   Here, we highlighted the abundance of SPP1+ macrophages, and their
   heightened interaction with colocalized T-cells through SPP1, MIF, and
   FN1 signaling in TIL-high HR+ breast cancer.

   The SPP1 gene encodes osteopontin, a matricellular protein previously
   associated with promoting angiogenesis in breast cancer^[149]26.
   Single-cell studies have indicated that SPP1+ macrophages interact with
   FAP+ fibroblasts in colon cancer, exhibiting a pro-angiogenic
   signature^[150]14,[151]21. However, our study revealed that SPP1+
   macrophages in breast cancer did not exhibit higher expression of
   angiogenic marker genes, and SPP1 expression in macrophages was not
   correlated with the expression of angiogenic marker genes. This pattern
   remained consistent across all breast cancer subtypes. Thus, in breast
   cancer, SPP1+ macrophages may serve roles other than angiogenesis.

   SPP1+ macrophages have previously been associated with fibrogenesis in
   chronic inflammation^[152]27,[153]28 and a higher metastatic rate in
   cancer^[154]29,[155]30. To examine how macrophage function affects
   TILs, we conducted analyses of DEGs. Both SPP1+ macrophages and
   monocytes/macrophages in TIL-high HR+ breast cancer showed up-regulated
   expression of genes related to ECM degradation compared to other
   monocytes/macrophages and those in the TIL-low group. Furthermore, the
   amount of collagen fibers was increased in TIL-high HR+ breast cancer
   compared to that in TIL-low HR+ breast cancer; furthermore, through
   spatial transcriptomics, the elevated expression of genes related to
   ECM remodeling in areas where SPP1+ macrophages and T-cells were
   present in high proximity was determined. These findings suggest an
   enhanced ECM remodeling process in TIL-high HR+ breast cancer,
   potentially facilitating tumor invasion and metastasis^[156]31.
   Although the precise mechanism underlying this up-regulated ECM
   remodeling in macrophages in TIL-high HR+ breast cancer remains unclear
   and the spatial transcriptomics analysis utilized only one sample per
   group, it could potentially contribute to the unfavorable prognosis
   associated with this subtype. Further validation in larger cohorts is
   required for the generalizability of our findings.

   In this study, we proposed a novel role of SPP1 within tumor-associated
   macrophages in breast cancer, shedding light on their primary
   interaction mediated by SPP1 signaling with T-cells. The SPP1 signaling
   pathway was associated with cancer-promoting pathways in gastric
   cancer^[157]32 and was linked to TAM infiltrations and poor prognosis
   in glioma^[158]33. However, these studies mainly focused on the role of
   SPP1+ macrophages in their interaction with tumor epithelial cells. In
   a recent study involving ovarian cancer, SPP1 expression was associated
   with higher TIL levels, but also with increased immune checkpoint
   marker gene expression and poor prognosis, potentially indicating
   immune tolerance^[159]34. Unfortunately, this study did not specify the
   cell types responsible for these expression changes. We could not
   identify any significant correlation between heightened SPP1 expression
   in macrophages and immune checkpoint marker gene expression or cell
   cycle progression scores in T-cells. Although T-cells in TIL-high HR+
   breast cancer were enriched in several interleukin signaling pathways
   associated with immunosuppressive TMEs like Interleukin −4, −10, and
   −13, they might help in attacking the tumors by promoting CD8T cell
   activation^[160]35,[161]36. Therefore, it is plausible that the impact
   of macrophage-derived SPP1 on T-cells may not be associated with immune
   suppression, and further research is needed to ascertain their precise
   role.

   The dynamic crosstalk between TAMs and TILs plays a pivotal role in
   orchestrating immune responses in various cancers. In this study, we
   highlighted three key signaling pathways (SPP1, FN1, and MIF) that
   mediate crucial interactions between TAMs and T-cells in TIL-high HR+
   breast cancer. MIF, identified as a contributor to immune-tolerant
   phenotypes in malignant tumors^[162]37, inhibits T-cell activation in
   the tumor context^[163]38,[164]39. Fibronectin, encoded by the FN1
   gene, is related to TAM infiltration and suppression of TIL
   infiltration^[165]40,[166]41, yet it can induce T-cell migration in
   dermal interstitium^[167]42,[168]43. In the context of
   macrophage-T-cell interactions, FN1 signaling might be implicated in
   shaping the immune-tolerant tumor microenvironment. Macrophages
   producing fibronectin contribute to the ECM structure, can influence
   T-cell adhesion and migration to TAMs, hindering their attack on tumor
   cells. Collectively, crosstalk between TAMs and T-cells in TIL-high HR+
   breast cancer may be associated with immune-tolerant conditions,
   despite the lack of correlation between macrophage-derived SPP1 and
   immune suppression in our results.

   In contrast to other breast cancer subtypes, patients having HR+ breast
   cancer with higher TIL levels are associated with worse
   prognosis^[169]7,[170]8. However, after receiving chemotherapy,
   patients having higher TIL levels show with higher pCR rates along with
   longer disease-free survival^[171]8,[172]44. While our study suggests
   increased SPP1 expression in macrophages as a unique characteristic of
   TIL-high HR+ breast cancer that may impact TILs, ECM remodeling, and
   prognosis, all patients with HR+ breast cancer included in our study
   were untreated, which limits our ability to assess the impact of
   chemotherapy on the TME and its association with prognosis. In
   addition, due to the small number of patients included in the study,
   there were no cases of breast cancer recurrence in HR+ breast cancer
   patients during an average follow-up period of 73.6 months (n = 7).
   Notably, recent research has indicated an elevated risk of recurrence
   in tamoxifen-treated HR+ breast cancer patients with high SPP1 mRNA
   expression in bulk tumor tissue^[173]45. This association may be
   attributed to the previously mentioned immunosuppressive and ECM
   remodeling effects of SPP1+ macrophages^[174]31,[175]34. Therefore, a
   potential strategy to enhance the prognosis of patients with TIL-high
   HR+ breast cancer could involve combining drugs that target SPP1 with
   the current chemotherapy regimen.

   Currently, there are no commercially available drugs for SPP1
   inhibition. However, based on our investigation using the Drug Gene
   Budger (DGB)^[176]46, daunorubicin was identified as a potential agent
   capable of reducing SPP1 expression in acute myeloid leukemia cell
   lines such as SKM1. Additionally, a review of the literature revealed
   that the TRK inhibitor Entrectinib can bind and inhibit SPP1, as
   demonstrated in a previous study^[177]47. Nonetheless, further research
   is imperative for drug development, as daunorubicin did not exhibit
   similar effects in other cell lines within the DGB database.

   In conclusion, we demonstrated that increased SPP1+ macrophages in
   TIL-high HR+ breast cancer may directly influence TILs, potentially
   enhancing ECM remodeling and T-cell interactions. Therefore, targeting
   SPP1+ macrophages could offer a promising approach to overcome the poor
   prognosis associated with TIL-high HR+ breast cancer.

Methods

Ethics

   This study received approval from the Institutional Review Board of
   Asan Medical Center (2016-0935) and adhered to the principles of the
   Declaration of Helsinki of 1975, as revised in 1983. All patients
   provided written informed consent.

Histologic analysis

   For histological examination, formalin-fixed and paraffin-embedded
   (FFPE) blocks were obtained from surgically resected breast specimens.
   The TIL value was defined as the proportion of the area occupied by
   mononuclear inflammatory cells within the tumor stroma^[178]4. In
   brief, a full assessment of average TILs in the stromal area within the
   borders of the invasive tumor (at least one FFPE block produced for
   every 1 cm increase in tumor size) was calculated by pathologists in
   one representative section per one FFPE block of 5-µm thickness
   (magnification 200–400×). Subtyping of breast cancer was based on
   immunostaining or silver in situ hybridization results of ER, PR, HER2,
   and Ki-67.

   To examine collagen deposition, FFPE sections were stained with
   Trichrome III Blue stain (Ventana Medical Systems, Tucson, AZ.) using a
   standard protocol for the NexES® automated slide staining system
   (Ventana Medical Systems). Then, the stained slides were scanned at x20
   with Philips Pathology Scanner SG300 (Philips Healthcare, Amsterdam,
   Netherlands.). The scanned image was converted to a tiled TIFF file and
   downloaded using Philips IntelliSite Pathology Solution 5.1. (Philips
   Healthcare). Then, the area with invasive carcinoma was selected by
   pathologists and classified as tumor, stroma, and glass/empty area
   using the DenseNet algorithm by HALO (v. 3.6.4134.137, Indica Labs,
   Albuquerque, NM) and HALO AI (v. 3.6.4134, Indica Labs). Positive
   stained areas in the stroma were measured using HALO Area
   Quantification v.2.4.3 through categorizing the intensity into three
   grades. H-score was obtained by summing up the positive stained area
   for each category (H-score = 3 × P[3] + 2 × P[2] + 1 × P[1]; where
   P[3], P[2], and P[1] denote the percentages of strongly, moderately,
   and weakly stained areas, respectively).

Single-cell dissociation

   Breast cancer tissues were placed in RPMI 1640 medium (11875093, Thermo
   Fisher Scientific, Waltham, MA, USA) and promptly transported to the
   laboratory within two hours post-resection. After washing with
   Dulbecco’s Phosphate Buffered Saline (DPBS, L0615, Biowest, Nuaillé,
   France) containing 1× ZellShield anticontaminant agent (Minerva
   Biolabs, Berlin, Germany), the tumor tissues were minced into
   1-mm-diameter pieces. The minced breast cancer tissue was incubated for
   1 h in digestion buffer (DMEM-F12) (11330032, Thermo Fisher
   Scientific), 2% FBS (Corning, VA, USA), 1% penicillin/streptomycin
   (15140-122, Thermo Fisher Scientific), 10 μg/mL insulin (51500056,
   Thermo Fisher Scientific), and 10 ng/mL epidermal growth factor
   (PHG0311, Thermo Fisher Scientific) supplemented with 1× collagenase/
   hyaluronidase (CA094-002, Gendepot, Barker, TX, USA). The digested
   samples were centrifuged at 80 × g for 30 s, and the resultant
   supernatant was filtered through a 70-μm nylon mesh strainer (352350,
   Corning Life Sciences, Tewksbury, MA, USA). The pellets from the
   digested samples were resuspended in 0.25% trypsin/EDTA (25200072,
   Thermo Fisher Scientific) and incubated for 3–5 min to release single
   cells. Once a single-cell suspension was obtained, cold Hank’s balanced
   salt solution (L0612, Biowest) containing 2% FBS was added, and cells
   were recovered through centrifugation at 300 × g for 5 min.
   Subsequently, the cells from the supernatants and the single cells from
   the pellet were combined, counted, and cryopreserved.

CD45+ cell isolation

   CD45+ cells were isolated either through magnetic bead sorting or
   fluorescence-activated cell sorting (FACS). In the magnetic separation
   method, dead cells were initially removed using a dead cell removal kit
   (130-090-101, Miltenyi Biotec, Auburn, CA, USA), followed by cell
   labeling with CD45-microbeads (130-045-801, Miltenyi Biotec) and
   enriched using a magnetic-activated cell sorting (MACS) column
   (130-042-901, Miltenyi Biotec) and MACS separator (130-091-051,
   Miltenyi Biotec). For FACS sorting, single cells were stained with an
   anti-human CD45 antibody (304016, Biolegend, San Diego, CA, USA) and
   resuspended in 4′,6-diamidino-2-phenylindole solution (D1306, Thermo
   Fisher Scientific). CD45+ live cells were sorted using FACSAria (BD
   bioscience, San Jose, CA, USA). The isolated cells were resuspended in
   CS10 freezing media (100-1061, Stemcell Technologies, Bothell, WA, USA)
   and stored in an NL2 tank until the next use or were immediately
   sequenced.

scRNA-seq library preparation and sequencing

   scRNA-seq library preparation was conducted using the Chromium
   Single-Cell 3′ Kit v2 and the Chromium Next GEM Single-Cell 3′ Kit v3
   platform (10x Genomics, Pleasanton, CA, USA). Cell suspension volumes
   were determined to attain a target capture of either 5000 or 10,000
   cells per sample. This process involved the generation of gel beads in
   emulsion (GEMs), barcoding, GEM-reverse transcription (RT) clean-up,
   complementary DNA amplification, and library construction, as per the
   manufacturer’s specified protocol. Library quality assessment was then
   performed using Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA).
   Subsequently, sequencing was executed using the NovaSeq 6000 System
   (Illumina, Inc., San Diego, CA, USA), generating 150 base pairs (bp)
   paired-end reads.

scRNA-seq data processing

   To increase the sample size and analyze independent datasets, we
   performed joint analysis of our dataset with a public dataset (GSA,
   HRA000477^[179]18), which was a single-cell RNA sequencing dataset
   where TIL values were measured and human CD45+ cells were sorted from
   breast cancer tissue. Analysis was conducted starting from the Cell
   Ranger processing stage with raw fastq files. The 10x Genomics Cell
   Ranger 6.1.2^[180]48 was used for mapping sequencing reads against the
   GRCh38 human reference (2020-A, GENCODE v32/Ensembl 98), unique
   molecular identifier (UMI) counting, and cell identification.
   Pre-processing using the Seurat package (version 4.1.1)^[181]19 was
   performed in R (version 4.1.3)^[182]49. In total, 31 samples were
   merged. Cells that had >20% mitochondrial genes or unique feature
   counts over 6000 or less than 200 were excluded. To normalize the
   dataset, we employed the “NormalizeData” procedure. Scaling was applied
   to the dataset with 2000 features selected through
   “FindVariableFeatures”. We then performed principal components analysis
   (PCA) on the scaled data using “RunPCA” (dims = 1:30). Harmony (version
   0.1.0)^[183]50 was employed for batch effect correction, using each
   sample as a batch assignment. Harmony embedding was used in downstream
   analysis. Uniform manifold approximation and projection (UMAP) was
   performed using “RunUMAP”.

Cell type annotation

   Azimuth package (version 0.4.5)^[184]19,[185]20 (reference = PBMC) was
   used for cell type annotation for CD45+ immune cells. The metadata
   column named ‘predicted.celltype.l2’ determined by Azimuth was utilized
   to isolate immune cells into four groups: T-cell
   (predicted.celltype.l2 = “CD4 CTL”, “CD4 Naive”, “CD4 Proliferating”,
   “CD4 TCM”, “CD4 TEM”, “CD8 Naive”, “CD8 Proliferating”, “CD8 TCM”, “CD8
   TEM”, “dnT”, “gdT”, “MAIT”, “NK”, “NK Proliferating”, “NK_CD56bright”,
   and “Treg”), B-cell (predicted.celltype.l2 = “B memory”, “B
   intermediate”, “B naive”, and “Plasmablast”), DC
   (predicted.celltype.l2 = “ASDC”, “cDC1”, “cDC2”, and “pDC”), and
   mono/macro (predicted.celltype.l2 = “CD14 Mono” and “CD16 Mono”).
   Further annotation was conducted for T-cells and mono/macro. Memory
   T-cells (predicted.celltype.l2 = “CD8 TCM”, “CD8 TEM”, “CD4 TCM”, and
   “CD4 TEM”) were defined based on marker expression. CCR7 and SELL
   double-positive cells (UMI count >0) were annotated as TCM, while
   others were annotated as TEM. Mono/macro cells
   (predicted.celltype.l2 = “CD14 Mono” and “CD16 Mono”) were separated
   and clustered (resolution 0.3). Each cluster was named based on
   expressed markers except mitochondrial genes. Expressed markers were
   identified using the “FindAllMarkers” function in Seurat with the
   parameter “pct = 0.5”. Following cell type annotation, each breast
   cancer subtype was individually analyzed downstream.

Pathway enrichment analysis

   Differentially expressed genes (DEGs) were computed for each group
   using the “FindMarkers” function in Seurat with default parameters.
   Specifically, we considered genes as having up-regulated expression if
   they exhibited an average log[2] fold change greater than 1 and a
   p-value less than 0.05, for further analysis. Pathway enrichment
   analysis was subsequently conducted using the “enrichPathway” method
   from the ReactomePA package (version 1.38.0)^[186]51 in R. Bar plots,
   visualizing significant pathways based on an adjusted p-value threshold
   of less than 0.05 (utilizing the Benjamini–Hochberg method), were
   generated using ggplot2^[187]52.

Cell–cell interaction analysis

   The CellChat package (version 1.6.1)^[188]22 was used to analyze
   cell–cell interactions. The CellChat object was created using the
   “createCellChat” function with the CellChatDB human dataset. T-cells
   were categorized into four groups: CD4T (“CD4 CTL”, “CD4 Naive”, “CD4
   Proliferating”, “CD4 TCM”, and “CD4 TEM”), CD8T (“CD8 Proliferating”,
   “CD8 TCM”, “CD8 TEM”, and “CD8 Naive”), NK (“NK”, “NK Proliferating”
   and “NK_CD56bright”), and OtherT (“dnT”, “gdT”, “MAIT”, and “Treg”). To
   identify differential interactions in cell–cell communication and to
   explore conserved and context-specific signaling pathways across
   different TIL levels, we compared TIL-high and TIL-low using
   “netVisual_diffInteractions” and “rankNet”. Then, the
   “netVisual_heatmap” function was employed to visualize each signaling
   pathway in different cell groups of varying TIL levels. Identification
   of ligand-receptor pairs of all the significant interactions was
   performed using “netVisual_chord_gene”.

Spatial transcriptomics data processing and signature scoring

   Seurat objects for each sample were created using the filtered count
   matrix and spatial data from the Zenodo data repository^[189]23 via the
   Seurat package. To normalize each sample, we employed the “SCTransform”
   function. TIL levels in each sample in H&E images were assessed by a
   certified pathologist.

   To identify cell type signatures, markers for each cell group were
   analyzed in scRNA-seq data using the “FindMarkers” function with
   default parameters. We selected the top 25 markers as signature genes.
   Signature scoring was conducted through the “AddModuleScore” method
   with a list of signature genes. To minimize the influence of variations
   in TIL quantity on the scores, we normalized each score. We scaled and
   centered the signature scores in R and subsequently employed them for
   correlation coefficient analysis.

   To calculate the relative proximity between T-cells and SPP1+
   macrophages, spots with both CD8T and CD4T scores above the median were
   selected. Among these spots, those with an SPP1+ macrophage score above
   the median were defined as proximity high spots, while those with
   scores below the median were designated as proximity low spots.

Cell cycle scoring

   To assign cell cycle scores to each cell, we used G2/M and S phase
   markers. This calculation was performed using the “CellCycleScoring”
   function, available in the Seurat package.

Survival analysis

   Survival analysis was carried out with Kaplan–Meier Plotter^[190]53
   ([191]https://kmplot.com/analysis). The overall survival was analyzed
   based on the ratio of two genes. The cut-off value was determined using
   the auto-select best cut-off mode.

Statistical analysis

   The Wilcoxon and Kruskal–Wallis rank sum tests were performed using the
   R stats package (ver.3.6.2)^[192]49. Additionally, Spearman’s rank
   correlation coefficient was calculated using the “stat_cor” function in
   the R package ggpubr (ver.0.6.0)^[193]54.

Supplementary information

   [194]Supplementary Information^ (12.6MB, pdf)

Acknowledgements