Spatial Transcriptomic Landscape of Brain Metastases from Triple-Negative Breast Cancer: Comparison of Primary Tumor and Brain Metastases Using Spatial Analysis.

Introduction
Brain metastases (BM) are the most prevalent type of brain tumor in adults, outnumbering primary brain tumors [1]. Common primary sites for BMs include lung, breast, melanoma, and renal cell carcinoma [2]. Among metastatic breast cancer patients, 30%-50% will develop BM during their illness, with prognosis and incidence varying according to the molecular subtype [3]. In triple-negative breast cancer (TNBC), roughly 30% of patients develop BM, generally experiencing poor clinical outcomes [3].
TNBC is defined by its lack of expression of estrogen receptor (ER) and progesterone receptor, along with non-amplified human epidermal growth factor receptor 2 (HER2) status. It represents an aggressive subtype of breast cancer, accounting for 10%-15% of all breast cancer cases [4]. In TNBC, mortality is chiefly associated with distant metastasis: the 5-year survival rate for localized disease is 91%, which declines to 12% upon development of distant metastasis [5]. Despite its aggressive nature, TNBC typically responds well to chemotherapy and immunotherapy, significantly affecting the long-term prognosis.
TNBC is considered the most immunogenic subtype among breast cancers and has attracted significant attention as a target for immunotherapy [6]. The level of tumor-infiltrating lymphocytes (TILs) is higher in TNBC compared to other subtypes. In the large-scale Breast International Group (BIG) 02-98 trial, the median percentage of tissue area infiltrated by TILs was 20% in TNBC, compared to 15% in HER2-positive tumors and 10% in ER-positive, HER2-negative tumors [7]. Trials combining immune checkpoint inhibitor (ICI) targeting the programmed death protein 1 (PD-1) or programmed death-ligand 1 (PD-L1) in TNBC have shown benefits in some patients; however, the cancer and tumor microenvironment (TME) cell ecosystems have been minimally investigated [8-10]. Specifically, the distribution of immune cells in BM following ICI treatment holds significant implications for the future treatment of metastatic breast cancer [11]. As functional immune statuses and tumor characteristics in BM are expected to differ from those in the primary tumor, spatially resolved analysis would provide a more comprehensive understanding of BM.
Recent advances in spatial RNA sequencing (spRNA-seq) have enabled comprehensive analysis of tumors and their surrounding microenvironment by quantifying spatially resolved mRNA expressions. Unlike single-cell RNA sequencing (scRNA-seq), which lacks spatial context, spRNA-seq offers valuable insights into gene expression landscapes by analyzing tumors at the tissue level. In this study, we aimed to dissect the spatial architecture of TNBC and paired BM by compartmentalizing tumor regions and analyzing the TME to assess the genomic status of each region.

Materials and Methods

1. Specimen preparation
Between 2015 and 2022, three TNBC patients with available tissue from primary breast cancer (PC) and corresponding BM were selected at Gangnam Severance Hospital. All three patients underwent neoadjuvant chemotherapy for PC, developed BM post-mastectomy, and subsequently had surgical removal of the BM. The disease course for each patient is depicted in Fig. 1. Clinical data, including age at diagnosis, the interval between breast surgery and BM onset, treatment regimens, and responses, were extracted from the electronic medical records (Table 1).
Histopathologic assessments, such as nuclear grade, histologic grade, and triple-negative status, were conducted by a pathologist (YJC). For each patient, the most representative formalin-fixed paraffin-embedded (FFPE) blocks of one post-treatment residual TNBC and one BM were selected following a pathologic review. For each sample, a 10.5×10.5 mm section of FFPE tissue was cut to a thickness of 5 μm and mounted on a polarized slide for subsequent spatial transcriptomics analysis.

2. Visium data analysis
To acquire whole transcriptome spatially barcoded sequence data from serial sections, we adhered to the Visium Spatial Gene Expression for FFPE protocol (Demonstrated Protocol CG000520). Tissue sectioning adhered to the Visium Spatial Gene Expression for FFPE – Tissue Preparation Guide (Demonstrated Protocol CG000518). The formalin-fixed, paraffin-embedded tissue samples were sectioned at 5 μm thickness and placed on Superfrost Plus Microscope Slides (Fisherbrand). The sections were then deparaffinized and stained with hematoxylin and eosin (H&E). After imaging, the sections underwent decoverslipping, hematoxylin destaining, and decrosslinking.
The tissue sections were transferred to a Visium CytAssist Spatial Gene Expression slide using the Visium CytAssist instrument. This slide, featuring a capture area of 1.21 cm2, facilitated analyte transfer. Following this, probe extension and library construction were conducted according to the Visium for FFPE standard workflow. Sequencing libraries were prepared using paired-end dual-indexing with settings of 28 cycles for Read 1, 10 cycles for i7 and i5 indices, and 90 cycles for Read 2. Libraries were demultiplexed using bcl2fastq (Illumina), and the resulting FASTQ files were analyzed using the Space Ranger pipeline v2.0.1 (10x Genomics) with the GRCh38-2020-A reference genome.

3. Pathologic annotation
Using Loupe Browser v6.2.0, pathologist (YJC) annotated paired PC and BM samples. For PC samples, tissue areas were categorized into six clusters based on their cellular structures: invasive carcinoma, ductal carcinoma in situ (DCIS), necrosis, tumor stroma, non-tumor stroma, and non-tumor epithelium.
Invasive carcinoma was characterized by infiltrating tumor epithelial cells that lacked surrounding myoepithelial cells. DCIS was identified by tumor epithelial cells confined within ducts and surrounded by myoepithelial cells. Necrotic areas, consisting of karyorrhectic cell debris typical of tumor necrosis, were annotated to be excluded from further analysis.
Tumor stroma was defined as the non-epithelial regions of the tumor tissue, including immune cells, mesenchymal cells, and vasculature constituting the TME. Non-tumor stroma included areas outside the tumor regions, devoid of epithelial components. Non-tumor epithelium comprised non-cancerous epithelial cells such as ductal and acinar cells from normal mammary glands adjacent to the tumor area on the same slide.
For BM samples, tissue areas were categorized into four clusters: metastatic carcinoma, non-epithelial tumor area, hemorrhage, and necrosis. For all three patients, glial tissue was not included in the specimens. Metastatic carcinoma encompassed the tumor epithelial cells of metastatic breast cancer. The non-epithelial tumor area consisted of regions of the metastatic tumor excluding the tumor epithelial cells. Necrosis and hemorrhage were annotated to segregate these compartments for analysis. The necrosis cluster contained tumor cell debris, and the hemorrhage cluster consisted of red blood cells.

4. Zone discrimination of specimen area by Visium data
To systematically analyze regional transcriptomics, we divided the peritumoral and non-tumor spots into two distinct distance zones, termed proximal TME and distal TME, respectively. We calculated the distances between the central coordinates of the peritumoral/non-tumor spots and tumor spots, as indicated by “pxl_row_in_fullres” and “pxl_col_in_fullres” from the tissue positions file. The diameter of each spot was 65 μm. The distance from each non-tumor spot to the nearest tumor spot was measured. Spots located within a predefined range 120 μm of tumor spots were defined as the proximal TME, corresponding to the peritumoral area, while spots located more than a specified distance 120 μm away from tumor spots were classified as the distal TME, representing the non-tumor neighboring area.
We conducted differentially expressing genes (DEGs) analysis using the FindMarkers() function in the Seurat package ver. 4.1.0 to generate the DEG list based on the regional categorization. With three zones (tumor, proximal TME, distal TME), three DEG lists were derived from the regional comparisons.

5. Cell deconvolution via TIMER
To calculate immune infiltration for each spot across the total of six Visium datasets, spot-specific count data was analyzed using TIMER 2.0 [12]. Among the various tools provided by TIMER 2.0, xCELL, which exhibited no missing values across all datasets, was utilized. The xCELL deconvolution technique facilitated the estimation of immune cell infiltration, including CD8+ T cells, CD4+ T cells, B cells, M1 macrophages and M2 macrophages in each spot. The estimated immune cell counts for each spot were then applied to the pathology slides for visualization. Furthermore, the immune cell estimations were categorized by zone and compared, with the degree of increase evaluated for each patient.

6. Identification of relative modules across the samples
To identify the core transcriptional programs shared between PC and BM, we filtered out cancer cells from the integrated datasets and performed sub-clustering. Marker gene sets for each cluster were obtained using the FindAllMarkers function in the Seurat package. To ensure robust gene set clusters, shared genes between different clusters were excluded, and clusters with fewer than 50 genes were filtered out based on adjusted p-values < 0.05 and log2FC > 0.25. This process yielded 20 clusters, with 11 from BM samples and nine from PC samples. Expression scores for the relevant gene set were calculated using the AddModuleScore function in Seurat. Pearson correlation analysis was then conducted on the module scores for the 20 clusters, and the resulting matrix was utilized for hierarchical clustering with the ComplexHeatmap package. A correlation coefficient threshold of > 0.6 was employed to identify significant modules.

7. Functional enrichment analysis
Enrichr was used to confirm the biological pathways associated with shared genes across the clusters in each module [13]. The shared genes, identified using the FindAllMarkers function (adjusted p < 0.05 and log2FC > 0.25), were submitted to the Enrichr portal (https://maayanlab.cloud/Enrichr/). Hallmark gene sets from the molecular signatures database (MSigDB) were utilized for functional enrichment analysis [14].
Additionally, a protein-protein interaction (PPI) network was constructed using the STRING database [15]. The 109 genes from module 2 were uploaded, and disconnected nodes were eliminated from the network. Visualization of the PPI network was achieved using Cytoscape [16].

8. Survival prognosis analysis for modules
The survival map was characterized using the GEPIA2 portal to assess the overall survival (OS) across five modules (http://gepia2.cancer-pku.cn/) [17]. The shared genes between DEGs of clusters from each module served as inputs for survival analysis, with the group cutoff set at median values. A p-value significance threshold of < 0.05 was applied. Survival analysis was conducted using the Cox proportional hazards model.

9. Identification of candidate drugs
To identify potential drug candidates targeting module 2, we utilized the Drug-Gene Interaction Database (DGIDB) (https://www.dgidb.org/), a comprehensive tool for exploring drug-gene interactions [18]. We uploaded the 18 genes from module 2 into DGIDB and filtered the results to include only Food and Drug Administration (FDA)–approved drugs, ensuring clinical relevance.

Results

1. Pathologic annotation and RNA expression
Regions within the breast were categorized as DCIS, invasive carcinoma, necrosis, tumor stroma, non-tumor stroma, or epithelium (Fig. 2A). The spatial distribution and proportion of labeled spots were visualized for each patient, demonstrating that DCIS was present only in patient 1 (Fig. 2B). Utilizing these pathological annotations, spatial transcriptome datasets from the three patients were integrated and visualized using a uniform manifold approximation and projection (UMAP) plot (Fig. 2C). Invasive cancer cells were distinctly separated from adjacent stroma, except in spots identified as DCIS and necrosis. These cancer cells were found to express genes, including PRAME and ETV7, which are known to be highly enriched in TNBC (Fig. 2D). Additionally, the DEGs in the tumor stroma identified genes associated with epithelial-mesenchymal transition (EMT) and tumor progression, such as COL10A1 and MMP13.
Regional categories of BM included invasive carcinoma, necrosis, tumor stroma, and hemorrhage, similar to those of PC tissue except for the absence of DCIS and non-tumor stroma (Fig. 2E and F). The UMAP displayed closely located areas of necrosis and invasive carcinoma, reflecting their shared molecular features (Fig. 2G). Meanwhile, tumor stroma and hemorrhage formed distinct clusters from invasive carcinoma. A heatmap of the DEGs revealed that invasive carcinoma expressed genes associated with breast cancer BM, such as KRT5 and GBP1, compared to the tumor stroma (Fig. 2H). DEG analysis between tumor tissues from BM and PC revealed significant transcriptional differences. Several genes were markedly upregulated in PC compared to BM, with PIGR, CCL11, and CCL21 showing the highest expression levels. Other notable genes included COL6A6, LRRC15, TAC1, and CD1C, along with immunoglobulin-related genes such as IGHV6-1, IGHA1, and JCHAIN. Conversely, multiple genes were significantly upregulated in BM compared to PC. Among them, MAGEA10, DCD, and FTHL17 exhibited the highest differential expression. Additionally, GFAP and SLC18A3, genes associated with neural and glial function, were enriched in BM. Other significantly upregulated genes, including ARG1, SPPL2C, LY6L, and GRM3, suggest potential adaptations of metastatic tumor cells to the brain microenvironment (S1 Fig.). These findings confirm the accuracy and specificity of the pathological annotation of each spatially defined region.

2. Tumor-proximity TME analysis
We analyzed the transcriptomic profiles of peritumoral and non-tumor regions in relation to the tumor across paired primary and BM samples. Non-tumor spots were categorized into proximal and distal groups (referred to as proximal and distal TMEs) based on their distance from the tumor, with ≤ 120 μm defined as proximal and > 120 μm as distal (Fig. 3A). These criteria were superimposed on the H&E image, as shown in Fig. 3B. In PC samples, these three regions displayed distinct transcriptional profiles (Fig. 3C). The expression of TNBC-specific tumor genes, such as RAB25 and CLDN7, progressively decreased from the proximal to the distal TME. Interestingly, matrix metalloproteinase (MMP) family members exhibited differential expression between the proximal and distal TME. MMP1 increased while MMP2 decreased in the proximal TME, in line with MMP1’s involvement in extracellular matrix remodeling supportive of cancer invasiveness and metastasis. In contrast, this expression pattern was reversed in the distal TME, where genes commonly expressed in the fibroblasts of normal breast stroma such as TNXB, COL14A1, and FBLN1 were predominantly expressed. Region-specific expression was less evident in brain metastasis samples; however, the proximal TME of BM was still differentiated from the tumor, characterized by the expression of neuronal genes involved in myelination and axonal transport, including GFAP, PLP1, APOD, SFRP4 and CLDN11 (Fig. 3D-F, S2 Fig.). Compared to the distal TME, the proximal TME in these samples prominently expressed genes typically found in epithelial tissues, such as S100A14, S100A2, KRT16, KRT6B, KLK10, and POSTN, which are normally absent in the healthy brain (Fig. 3G and H, S3 Fig.).

3. Deconvolution
The TIMER algorithm was used to estimate the infiltration of immune cells by zone, specifically targeting CD8+ T cells, regulatory T cells, M1 macrophages, and M2 macrophages (S4 Fig.). The analysis scheme is depicted in Fig. 4A. Fig. 4B and E demonstrate the estimated zonal distribution of immune cells in the BM of patient 2 and patient 3, respectively. Spatial heatmaps visualize the relative infiltration of CD8+ T cells and M2 macrophages in patient 2 (Fig. 4C and D) and patient 3 (Fig. 4F and G). In the tumor region, the estimation of CD8+ T cells showed a significantly higher rate in the BM compared to the PC in patient 3, who underwent immunotherapy, than in patients 1 and 2, who were treated with conventional chemotherapy (p < 0.001) (Fig. 4H). Notably, in the spatial context for patient 2, there was no increase in CD8+ T cell estimation within the tumor region; rather, an increase was noted in the distal TME. Conversely, in patient 3, an increase in CD8+ T cells was observed in the tumor region (Fig. 4F and H). In the distal TME region, the estimation of M2 macrophages showed a decreased pattern in the BM compared to the PC in patients 1 and 2, while a significant increase was observed in patient 3 (p < 0.001) (Fig. 4I). A detailed regional distribution of different immune cells is presented in S4 Fig.

4. Tumor-specific modules and clinical relevance
To identify the core transcriptional signatures associated with BM in three patients, we first derived DEGs from cancer clusters found in primary and metastatic breast cancer datasets. Genes common to multiple clusters were excluded, and clusters with fewer than 50 DEGs were omitted from the analysis. We selected 20 clusters (9 from PC and 11 from BM) for integration. Hierarchical clustering and correlation analysis were applied to identify significant modules based on a correlation coefficient threshold of R > 0.6 (Fig. 5A). Five modules were confirmed, enriched with genes involved in hypoxia, glycolysis, immune response, cell cycle regulation, and EMT (Fig. 5B). Among them, the glycolysis-related program (M2) comprising 109 significant genes showed high expression across the six samples (p < 0.05 and average log2FC > 0.25) (Fig. 5C). To identify hub genes within M2, we constructed a PPI network using the 109 genes as nodes via the STRING database. Fourteen hub genes were identified in the network, overlapping with an 18-gene list based on adjusted p < 0.05 in M2 (Fig. 5D). Kaplan-Meier survival analysis revealed that these 18 genes in module M2 (adjusted p < 0.05 and average log2FC > 0.25) significantly impacted survival outcomes in the breast cancer cohort of The Cancer Genome Atlas (TCGA)/Genotype-Tissue Expression (GTEx) data, with a hazard ratio of 1.4 and a p-value of 0.045 (Fig. 5E). To identify potential drug candidates targeting the 18 genes in M2, we conducted drug repurposing analysis using the DGIDB database, validated as an effective drug screening tool. The analysis identified 25 potential drug candidates associated with four genes: VEGFA, RAC1, EGLN3, and CAMK1D (Fig. 5F).

Discussion
Despite decades of research aimed at understanding the molecular basis of BM in breast cancer, the underlying mechanisms remain unclear, and treatment outcomes are still poor. Single-cell technologies have revolutionized the field by providing high-resolution insights into various cellular phenotypes based on gene expression but have lacked spatial information. In this study, we employed spRNA-seq to analyze residual tumor in post-treated TNBC and paired BM tissues, aiming to investigate the biological nature of tumors across different sites and the molecular alterations in the adjacent non-tumor area at metastatic sites. Through zonal categorization based on proximity to tumor regions, we identified transcriptional differences in the TME, focusing on immune-related genomic features.
The proximal TME of BM exhibited heightened activity in central nervous system processes, such as myelination and trans-synaptic signaling compared to the tumor epicenter, thereby indirectly suggesting neuron-tumor interactions at the tumor margin. Venkataramani et al. [19] demonstrated that excitatory synapses form between neuronal cells and cells of metastatic melanoma and breast cancer. Additionally, Zeng et al. [20] reported that activation of the N-methyl-D-aspartate receptor by glutamate ligands is linked with poor prognosis in cancer. Similarly, we noted upregulation of GRIA4 (glutamate ionotropic receptor AMPA type subunit 4) in the proximal TME, providing additional clinical evidence that synaptic molecule upregulation may facilitate the colonization of metastatic TNBC cells in the brain.
The introduction of ICB has profoundly transformed systemic treatment approaches for various solid tumors. Within breast cancer, TNBC is identified as the subtype most likely to benefit from ICBs, owing to its relatively higher immunogenicity [21]. The IMpassion130 trial has established atezolizumab in combination with nab-paclitaxel as the preferred first-line treatment for patients with advanced TNBC who are PD-L1–positive (SP-142 ≥ 1% in immune cells) [22]. Elevated TILs level have been associated with improved responses to anti–PD-1/PD-L1 therapies in metastatic TNBC. In the IMpassion130 trial, patients with PD-L1 positivity and TILs ≥ 10% exhibited better progression-free survival and OS. The phase III KEYNOTE-119 trial further confirmed the superior therapeutic effect of pembrolizumab in metastatic TNBC patients with TILs ≥ 5% [21]. These findings are consistent with earlier results from the phase II KEYNOTE-086 trial and previous phase I studies [23,24]. Although high TILs level is a robust predictor of ICB response, evidence on whether ICBs directly enhance TILs levels in addition to reinvigorating their function remains limited [25]. In our study, longitudinal tissue analysis revealed increased CD8+ T cell infiltration in the BM following ICB therapy in patient 3, relative to breast cancer tissue collected before ICB treatment. However, conventional chemotherapy did not seem to recruit CD8+ T cells in patients 1 and 2. Further studies involving larger patient cohorts are needed to confirm the activation of systemic T cell responses following ICB treatment. Additionally, as immune cell distribution may be influenced by factors such as PD-L1 expression and other immune profiles, further investigation is necessary to elucidate these interactions.
We also identified core transcriptional programs common to both primary and metastatic breast cancer, including those associated with glycolysis, immune response, cell cycle, and EMT. Among these, key genes in the glycolysis-related module, such as VEGFA and RAC1, were linked to poorer OS in the breast cancer cohort of the TCGA/GTEx data. Previous studies have shown that breast cancer cells undergo metabolic reprogramming during EMT, increasing glycolysis and lipid metabolism to enhance their invasiveness [26]. While luminal subtypes display either the reverse Warburg effect or a null phenotype, TNBCs predominantly rely on glycolysis [27]. Furthermore, the upregulation of key glycolytic enzymes, such as HK2 and glucose transporter 3, has been observed in BM [28,29], indicating that glycolysis may serve as an alternative energy source to support the colonization and survival of metastatic cancer cells in the brain.
spRNA-seq has inherent limitations in deriving detailed cellular conclusions. Unlike scRNA-seq, it does not enable direct analysis of cell-to-cell interactions, identification of subpopulations within specific cell types, or measurement of dynamic cellular changes. However, its key strength lies in its ability to quantitatively assess spatial relationships between cells and define specific regions based on spatial proximity. Using this approach, we successfully distinguished non-tumor regions based on their distance from tumor spots. Despite the small sample size, we noted a significant genomic difference between the proximal and distal TMEs, highlighting the need for thorough consideration in pathological non-tumor areas. By analyzing the proximal and distal TMEs, we identified two distinct expression patterns relating to neuronal and epithelial genes. Furthermore, employing spot-specific count data for cell deconvolution techniques enables the quantitative analysis of immune cell infiltration. However, the spatial resolution of Visium, which utilizes 55-micron spot sizes, offers an average resolution of 1 to 10 cells per spot. This limitation indicates that with Visium data alone, we can only estimate immune cell proportions. To overcome this, platforms such as Xenium in situ spatial transcriptomics or multiplex immunohistochemistry could prove beneficial.
This study has several limitations. Firstly, the small number of patients and the heterogeneity in their treatment regimens may limit the generalizability of the findings. Additionally, the analysis was conducted without single-cell resolution, which could diminish the granularity of insights into the TME. Another limitation is the inability to analyze pretreatment biopsy samples, which could have provided valuable baseline information. Finally, the study focused on post-treatment residual disease, and all patients exhibited residual disease following neoadjuvant therapy, suggesting that these patients were already predisposed to poor outcomes, further complicating the interpretation of the results.
In conclusion, we found that the microenvironments of PC and BM can be transcriptionally distinguished based on the relative location of tissue components and their proximity to the tumor. In BM, gene sets related to synaptic connections with neurons may serve as a crucial niche for metastatic tumor cells. Additionally, alterations in immune-related genomic features, driven by ICB, could suggest potential therapeutic options for patients with specific immune profiles. While further research is needed to validate some of our findings, we have contributed to a better understanding of the spatial context within the BM TME.