Abstract

   Breast cancers present intricate microenvironments comprising
   heterotypic cellular interactions, yet a comprehensive spatial map
   remained to be established. Here, we employed the DNA nanoball-based
   genome-wide in situ sequencing (Stereo-seq) to visualize the geospatial
   architecture of 30 primary breast tumors and metastatic lymph nodes
   across different molecular subtypes. This unprecedented high-resolution
   atlas unveils the fine structure of the tumor vasculature, highlighting
   heterogeneity in phenotype, spatial distribution, and intercellular
   communication within both endothelial and perivascular cells. In
   particular, venular smooth muscle cells are identified as the primary
   source of CCL21/CCL19 within the microenvironment. In collaboration
   with ACKR1-positive endothelial cells, they create a chemokine-rich
   venular niche to synergistically promote lymphocyte extravasation into
   tumors. High venule density predicts increased immune infiltration and
   improved clinical outcomes. This study provides a detailed spatial
   landscape of human breast cancer, offering key insights into the
   venular regulation of tumor immune infiltration.

   Subject terms: Breast cancer, Cancer microenvironment, Immunoediting,
   Next-generation sequencing
     __________________________________________________________________

   This study utilizes Stereo-seq to comprehensively map the geospatial
   architecture of breast cancer, revealing venular niches that promote
   immune infiltration and better clinical outcomes, offering new insights
   into tumor microenvironment regulation.

Introduction

   Breast cancer is a highly heterogeneous collection of tumors that are
   characterized by distinct molecular subtypes^[86]1. Advances in the
   subtype-specific therapies have yielded significant survival benefits
   for patients, primarily through the administration of hormone and
   HER2-targeted therapy^[87]2. In recent years, there has been a growing
   interest in the combination of immune checkpoint inhibitors (ICIs) and
   chemotherapy as a treatment for breast cancer^[88]3. The success of
   several clinical trials, including IMpassion130^[89]4,
   KEYNOTE-355^[90]5, and KEYNOTE-522 trials^[91]6, has led to the
   approval of atezolizumab and pembrolizumab for triple-negative breast
   cancer patients with PD-L1 expression. Despite these promising
   achievements, the efficacy of ICIs in treating other molecular subtypes
   and PD-L1-negative patients remains limited, partly due to inadequate
   immune infiltration in the tumor microenvironment (TME)^[92]3,[93]7. A
   deeper understanding of TME heterogeneity and the mechanisms dictating
   immune infiltration is required to optimize immunotherapy efficacy.

   The post-capillary venules have been identified as the major channel
   for lymphocytes extravasation through the endothelial
   barrier^[94]8,[95]9. Mechanistically, circulating lymphocytes engage
   venular endothelial cells (ECs) under shear stress via the binding of
   endothelial selectins, P-selectin and E-selectin, leading to the
   tethering and rolling of lymphocytes along the vessel wall^[96]8.
   Rolling lymphocytes are then activated through the interaction between
   CCR7 and homeostatic chemokines, including CCL21 and CCL19, that are
   immobilized on the EC surface^[97]10. Subsequently, activation of
   lymphocytes mediates firm adhesion and enables lymphocytes to crawl on
   the luminal surface of venules until they transmigrate through the
   barriers of ECs, associated pericytes/vascular smooth muscle cells
   (SMCs), and the basement membrane^[98]11.

   Previous single-cell and experimental studies have provided valuable
   insights into the interaction between ECs and
   lymphocytes^[99]12–[100]15, but ECs present significant plasticity and
   spatial heterogeneity that are coordinated with different functional
   states^[101]12. The spatial organization of venular ECs and other EC
   functional subtypes within intact tumors from patients in a clinical
   setting remains unclear. More importantly, our understanding of their
   interaction and coordination with perivascular cells in mediating
   immune infiltration is still limited, primarily due to limitations in
   technical precision and a lack of suitable experimental
   models^[102]11,[103]16. A high-resolution, unbiased spatial map of
   the tumor ecosystem and its vasculature could help to elucidate the
   complex spatial and functional networks between the immune and vascular
   systems.

   The application of spatial transcriptomics (ST) technologies has
   revolutionized our ability to retrieve spatial information that is lost
   in single-cell RNA sequencing (scRNA-seq)^[104]17. Wu et al. performed
   ST on six breast cancers to reveal the coordination of tumor and host
   cell phenotypes^[105]18. However, conventional ST methods, such as
   Visium, exhibit insufficient resolution due to their reliance on 55 μm
   barcoded spots, with a center-to-center distance of 100 μm, that
   encompass dozens of cells, thus impeding precise localization of
   dispersed ECs and perivascular cells within the TME. An innovative
   technology, SpaTial Enhanced REsolution Omics-sequencing
   (Stereo-seq)^[106]19, has recently been developed that combines 220 nm
   DNA nanoball (DNB)-patterned arrays and in situ tissue RNA capture to
   create spatial maps at cellular resolution, making it an ideal tool to
   depict the intricate tumor vasculature^[107]20.

   In this study, we applied Stereo-seq to characterize the geospatial
   architecture and cell-cell communication in 30 surgically resected
   specimens of primary breast tumors and paired lymph node metastases. We
   unveiled that tumor vasculature displayed distinct heterogeneity of
   biological phenotype, spatial organization, and interactions with
   surrounding cells. CCL21+ venular SMCs and ACKR1+ venular ECs
   cooperated through multiple chemokines to promote lymphocyte
   infiltration in breast cancer.

Results

Depicting a high-resolution spatial atlas of human breast cancer

   To comprehensively characterize the spatial landscape of breast cancer,
   we employed Stereo-seq on 30 fresh frozen tissues derived from 23
   breast cancer patients, including 23 primary breast cancers and 7
   paired metastatic lymph nodes that collectively covered all four
   molecular subtypes (luminal A: n = 8; luminal B: n = 9; HER2: n = 8;
   triple-negative: n = 5). Adjacent tumor tissues were subjected to
   scRNA-seq to facilitate cell type identification in the ST slides. We
   further applied immunohistochemistry (IHC) and multiplex
   immunofluorescence (mIF) on the adjacent slides of all 30 ST slides and
   on the validation cohort to confirm our discovery (Fig. [108]1a).
   Detailed clinical and pathological information of patients is provided
   in the Supplementary Data [109]1.

Fig. 1. A spatially resolved atlas of human breast cancer slides across
different molecular subtypes and tumor sites.

   [110]Fig. 1
   [111]Open in a new tab

   a Schematic diagram illustrating the acquisition and analysis workflow
   of the ST maps. Stereo-seq was performed on 30 slides of primary breast
   tumors and paired lymph node metastases from 23 patients, yielding a
   total of 4,894,305,510 DNBs. ScRNA-seq was performed on the fresh tumor
   samples from 6 out of 23 patients and identified 78,848 qualified
   cells. Each block represents a ST slide and is colored by molecular
   subtype. Red, blue, orange, and green spots represent venular
   endothelial cells (vECs), venular smooth muscle cells (vSMCs), tumor
   cells, and immune cells, respectively, in the projection of scRNA-seq
   and ST data. Schematic plots were created using bioRender. b Uniform
   manifold approximation and projection (UMAP) plot of the scRNA-seq
   dataset. Each cluster is colored by its corresponding cell type. EC:
   endothelial cell; SMC: smooth muscle cell, including vascular SMC and
   pericyte; T/NK: T or natural killer cell. c Differential expression of
   canonical markers in each cell type cluster. Macro: macrophage; Fibro:
   fibroblast; Melano: melanocyte; Granulo: granulocyte; PB: plasmablast.
   d–f Hematoxylin and eosin (H&E) staining (d), gene count maps (e), and
   cell type annotation maps (f) of the ST slides of the primary tumor
   (PT) and metastatic lymph node (LN) from the patient NCC-P05. ST spots
   are mapped to the H&E staining of adjacent slides. Raw spatial
   expression matrix of each bin/spot is convoluted into pseudo 50 spot
   with 25 µm squares (50 × 50 bins/spot, bin50 spot for short). IDC:
   invasive ductal carcinoma; DCIS: ductal carcinoma in situ. g The cell
   type annotation, H&E staining, gene counts, and unique molecular
   identifiers (UMI) counts in the selected area of NCC-BR5 slide
   highlighted in Fig. 1d–f. h Relative proportions of cell types across
   30 ST slides. ST slides are grouped according to the molecular subtype
   and tumor site.

   After quality control (Supplementary Fig. [112]1a-c), the scRNA-seq
   analysis of 6 breast cancer samples identified 78,848 qualified cells,
   which were categorized into 13 dominant cell types based on the
   canonical cell markers, including breast tumor cell (EPCAM), T or
   natural killer (T/NK) cell (TRBC2), B cell (CD79A), plasma cell
   (JCHAIN), plasmablast (GZMB and JCHAIN), macrophage (CD163), fibroblast
   (DCN), EC (VWF), vascular SMC/pericyte (referred to as SMC; NOTCH3),
   basal cell (COL17A1), neuron (CNTN5), granulocyte (CPA3), and
   melanocyte (CDH19) (Fig. [113]1b, c). In order to validate the
   robustness of cell classification, we conducted cell clustering and
   cell type annotation on an external scRNA-seq dataset of 26 breast
   cancer samples consisting of 88,876 cells^[114]18, which produced
   similar clustering results and feature genes for the dominant cell
   types (Supplementary Fig. [115]1d, e).

   For the ST data, Stereo-seq of 30 slides yielded a total of
   4,894,305,510 DNBs (220 nm × 220 nm), with an average of 163,143,517
   DNBs per slide. Utilizing a resolution similar to that of Visium,
   Stereo-seq captured an average of 7565 pseudo-spots (bin200 spot,
   200 × 200 DNB bins) per slide. Each spot, on average, quantified
   approximately 6120 genes (Supplementary Fig. [116]1f). To gain a more
   detailed view of slides with sufficient genes for spatial annotation,
   we convoluted the raw spatial expression data of each slide into
   smaller pseudo-spots (bin50 spot, 50 × 50 DNB bins, 25 μm diameter)
   (Fig. [117]1d-g), with each bin50 spot covering approximately 3-5 tumor
   cells and having an average detection of 921 genes and 2058 UMIs per
   spot (Fig. [118]1g) (Supplementary Fig. [119]1g). Random sampling of
   10,000 bin50 spots taken from each slide showed few batch effects
   (Supplementary Fig. [120]1h).

   Using the internal and external scRNA-seq results as independent
   reference datasets, we utilized the Robust Cell Type Decomposition
   (RCTD) algorithm to infer the cell composition for each bin50
   spot^[121]21. Annotation of eight major cell types, including tumor
   cell, T/NK cell, B cell, macrophage, fibroblast, EC, SMC, and basal
   cell, produced highly correlated results in both references