Spatial transcriptomics in breast cancer reveals tumour microenvironment-driven drug responses and clonal therapeutic heterogeneity

Abstract

Breast cancer patients are categorized into three subtypes with distinct treatment approaches. Precision oncology has increased patient outcomes by targeting the specific molecular alterations of tumours, yet challenges remain. Treatment failure persists due to the coexistence of several malignant subpopulations with different drug sensitivities within the same tumour, a phenomenon known as intratumour heterogeneity (ITH). This heterogeneity has been extensively studied from a tumour-centric view, but recent insights underscore the role of the tumour microenvironment in treatment response. Our research utilizes spatial transcriptomics data from breast cancer patients to predict drug sensitivity. We observe diverse response patterns across tumour, interphase and microenvironment regions, unveiling a sensitivity and functional gradient from the tumour core to the periphery. Moreover, we find tumour therapeutic clusters with different drug responses associated with distinct biological functions driven by unique ligand-receptor interactions. Importantly, we identify genetically identical subclones with different responses depending on their location within the tumour ducts. This research underscores the significance of considering the distance from the tumour core and microenvironment composition when identifying suitable treatments to target ITH. Our findings provide critical insights into optimizing therapeutic strategies, highlighting the necessity of a comprehensive understanding of tumour biology for effective cancer treatment.

Graphical Abstract

Open in new tabDownload slide

Introduction

Cancer is a group of genetic diseases that has been classically approached with a one-size-fits-all strategy, which often resulted in treatment inefficacy and adverse drug reactions (1), underlining the existence of intertumour heterogeneity (2). Personalized precision oncology is a novel therapeutic strategy that aims to provide tailored treatments to cancer patients according to the molecular characteristics of individual tumours (3,4). Breast cancer research has contributed to the foundations of personalized precision oncology by identifying biomarkers to guide patient stratification and selection of targeted therapies. In the clinical setting, an immunohistochemical analysis determines the expression of the oestrogen receptor (ER) and the progesterone receptor (PR) and the overexpression or amplification of the human epidermal growth factor receptor 2 (HER2). The status of these three hormone receptors guides the stratification of breast cancer patients into three groups with different prognosis and treatment strategies: luminal (ER+, PR±, HER2−); HER2+ (ER±, PR±, HER2+) and triple-negative breast cancer (TNBC; ER-, PR-, HER2−) (5). ER+ and HER2+ patients are treated with targeted therapies against these biomarkers, including hormonal therapies like tamoxifen or HER2 inhibitors. On the contrary, TNBC patients are usually prescribed chemotherapies since targeted therapies are still unavailable (6). Since the adoption of this strategy, breast cancer patient outcomes have improved. However, there is still a high number of treatment failure cases and relapses, especially in TNBC patients (7).

Several cancer cell subpopulations coexist within a single lesion and display different epigenomics, genomics and transcriptomics profiles. This variability, usually referred to as intratumour heterogeneity (ITH), reflects the existence of distinct tumour subpopulations that may harbour different sensitivities to anticancer therapies (8) and have been linked to treatment failure and relapse (9). Until recently, ITH research has mainly focused on cancer cells and their genomic and gene expression variability, thus studying ITH at the subclone or transcriptional cell state levels (10,11). In the last few years, the role of the tumour microenvironment (TME) on ITH and response to treatment has been acknowledged. The TME, which aggregates cellular and non-cellular components, such as immune cells or the extracellular matrix (ECM) (12), influences transcriptomics changes in the tumour compartment. At the same time, the tumour hijacks the TME to favour its own growth (13). Thus, this active cross-talk within the tumour ecosystem is a source of selective pressure for cancer evolution and, ultimately, a key factor in treatment response (14). Despite recent efforts, there is a lack of studies relating drug response with the spatial organization of the tumour and the cross-talk with the TME.

Single-cell technologies, particularly single-cell RNA-seq (scRNA-seq), have become a prominent method for studying ITH and TME. scRNA-seq methods require sample disgregation, causing cell stress (15) and death and not preserving the spatial organization of the tissue. However, carcinomas grow and evolve within a spatial context (16), so studying the microanatomical niches within a tissue can illuminate unexplored ITH sources. In this scenario, an incipient sequencing-based technology known as spatial transcriptomics (ST) is becoming more relevant since it maps high-resolution transcriptomics data on top of tissue slides, allowing the study of the distribution of cell types and their cell-to-cell communications, the border between the tumour and TME compartments and the localization of different tumour subpopulations (17).

In this work, we acquired ST data from nine patients with invasive adenocarcinomas stratified into luminal, HER2+ or TNBC subtypes. We predicted the sensitivity to >1200 drugs to explore the therapeutic heterogeneity of breast cancer while accounting for the spatial context and the interaction between the tumour and TME compartments. For this purpose, we have exploited a new version of Beyondcell (18), a tool for identifying tumour cell subpopulations with different drug response patterns in scRNA-seq data that now allows users to work with 10x Visium ST data.

Materials and methods

Visium ST data preprocessing

We created Seurat v4.3.0.1 (19) objects from raw counts, the corresponding metadata and tissue images. Spots marked as ‘Artefact’ were filtered out, as well as low-quality spots that did not meet any of the thresholds specified in Supplementary Table S1. Next, we normalized the filtered raw counts with the SCTransform function (vst.flavor=‘v1’). We assessed the cell cycle effect using the CellCycleScoring function and plotting a Principal Component Analysis (npcs = 50). Based on this plot, we regressed the cell cycle when needed using the SCTransform function with vars.to.regress = c(‘S.Score’, ‘G2M.Score’). In multi-region patients (CID44971, CID4465 and CID4535), the effect of the region was also regressed.

Spot deconvolution

We used Robust Cell Type Decomposition v2.2.1 (20) to predict the proportion of different cell types contained within each spot. First, we constructed three subtype-specific deconvolution references using annotated scRNA-seq data coming from 26 primary tumours (11 luminal, 5 HER2+ and 10 TNBCs) (5). Labelled cell types comprised B cells, cancer-associated fibroblasts (CAFs), cancer cells, endothelial cells, myeloid cells, normal epithelial cells, plasmablasts, perivascular-like (PVL) cells and T cells. As part of quality control preprocessing, we removed single cells that met any of the following criteria: (i) percentage of transcripts mapping to mitochondrial genes ≥15%; (ii) percentage of transcripts mapping to ribosomal genes ≥40%; (iii) number of transcripts ≤250 or ≥50 000; (iv) number of expressed genes ≤200 or ≥7500. Then, we merged the scRNA-seq datasets by breast cancer subtype and input the filtered raw counts to the Reference function to compute subtype-specific deconvolution references using a maximum of 500 cells per cell type. Finally, we deconvoluted the spots in each sample using the corresponding reference, the filtered raw counts, the coordinates and the default parameters of create.RCTD and run.RCTD functions.

Estimation of tumour purity and spot categorization

To estimate the tumour purity of each spot, we computed ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumours using Expression data) v1.0.13 scores (21) using the SCTransform-normalized and regressed counts obtained with Seurat v4.3.0.1. We scaled this score between 0 and 1 across all patient spots.

Then, we categorized each spot as either tumour or TME based on the agreement of three annotation sources: (i) the histopathological annotations when available (spots labelled as ‘Ductal carcinoma in situ’ or containing the ‘cancer’ or ‘carcinoma’ keywords), (ii) the scaled ESTIMATE score (≤0.4), which is inversely proportional to tumour purity, and (iii) the proportion of deconvoluted cancer cells (≥0.6) in each spot. The spots that did not meet at least two of these criteria were labelled as TME and assigned the non-malignant cell type with maximum deconvoluted proportion.

Subclone inference

We determined sample-wise clonal structures from inferred copy-number alteration (CNA) profiles using Single CEll Variational Aneuploidy aNalysis v1.0.1 (22). Briefly, we input filtered raw counts to pipelineCNA function (beta_vega = 0.5, SUBCLONES = TRUE, ClonalCN = TRUE, plotTree = TRUE) using the spots labelled as TME as reference. TME spots were relabelled as tumour if their inferred CNA profile was similar to other tumour spots. Moreover, tumour spots with no inferred CNAs were marked as diploid.

Breast sensitivity signature collection generation

The Breast sensitivity signature collection (SSc breast) contains 1372 gene signatures that reflect the transcriptional differences between sensitive and resistant breast cancer cell lines before drug treatment. In order to obtain these signatures, we performed a differential expression analysis against the area under the curve with limma v3.54.0 (23) for all compounds tested in ≥10 different breast cancer cell lines. We selected the top 250 upregulated and downregulated genes in sensitive versus resistant cancer cell lines, ranked by the t-statistic, to create bidirectional gene signatures of 500 genes each. The area under the curve was used to measure drug response because, contrary to IC50, it can always be estimated without extrapolation from the dose-response curve and has shown more accuracy in predicting drug response (24).

Expression and drug response data were retrieved from three independent pharmacogenomics assays: the Cancer Therapeutics Response Portal v2 (25–27), the Genomics of Drug Sensitivity in Cancer v2 (28–30) and the PRISM (31,32) repurposing compendium through the DepMap portal v22Q4 (33). As these sources are independent, several signatures refer to the same compound. Consequently, the 1372 transcriptomic signatures that form the SSc breast reflect the predicted response to >1200 drugs. To verify that the cancer cell lines included in the signature generation analyses were representative of breast cancer patients, we used the corrected lineage reported in the Celligner project (34).

Beyondcell score calculation and correlation

Beyondcell Scores (BCS) constitute a spot-wise enrichment metric that is normalized to penalize spots with many zeros or outliers. A positive BCS indicates enrichment in the upregulated genes that comprise the signature. Conversely, a negative BCS designates spots enriched in the downregulated genes. BCS close to 0 indicate uncertainty on the enrichment directionality. As the SSc breast was computed comparing treatment-naive sensitive versus insensitive cells, the BCS reflects the predicted sensitivity (normalized BCS > 0) or insensitivity (normalized BCS < 0) to the corresponding drug. Similarly, a normalized BCS > 0 indicates enrichment in the phenotype represented by a functional signature.

BCS were computed using the SCTransform-normalized counts and the bcScore function (expr.thres = 0.1) from the beyondcell v2.2.0 package (18). The non-available values in the normalized BCS matrix were transformed to 0. Then, we combined all patient-wise results into a single object and regressed the effect of patient identity with bcRegressOut function.

Therapeutic cluster computation

We computed therapeutic clusters (TCs) with the bcUMAP function from the beyondcell v2.2.0 package. This function performs a dimensionality reduction via Principal Component Analysis (npcs = 50) and constructs a K-Nearest Neighbors graph (k.neighbors = 20) based on the Euclidean distance in the Principal Component Analysis space. The top principal components (pc) for clustering are determined by drawing an elbow plot, and the spots are grouped using the Louvain algorithm, specifying different resolutions (res). The TCs and the Uniform Manifold Approximation and Projection computed by bcUMAP were visualized with the bcClusters function. We selected a res = 0.5 and pc = 20 for all spots and a res = 0.15 and pc = 10 for tumour spots.

To identify groups with similar drug sensitivity, we computed a vector with the mean BCS per drug signature for each patient and (i) major TC (for all spots) or (ii) tumour TC (for tumour spots), excluding groups with less than 50 spots or mixed TCs. Then, we calculated the Pearson correlation coefficients between vectors and clustered the results according to Ward’s method.

Neighbourhood analysis

For each patient, we selected the spots located at the inner and outer edge of each patient’s major TC using the RegionNeighbors function (mode=‘inner_outer’) from the semla v1.1.6 package (35). Then, we performed a neighbourhood enrichment analysis to assess whether the spots belonging to two major TCs co-localized more than expected by chance. RunNeighborhoodEnrichmentTest transforms the Visium data in a network, with each spot constituting a node with a maximum of six edges with other spots. For each pair of major TCs, this function randomly permutes the class labels a fixed number of times (n.perm = 1000), calculates a null distribution with the number of edges between spots from different classes and returns a z-score computed with the number of observed edges and the mean and standard deviation of the null distribution. A z-score around 0 can be interpreted as random spot localization. A positive z-score indicates an over-representation of the label pair co-localization. In contrast, a negative z-score can be viewed as a spatial repellent effect of the label pair. Finally, we aggregated the results across samples by computing the mean z-scores of each major TC pair.

Radial distances

We computed the radial distances to the centre of the tumour-rich TCs, which we designated as the tumour core, with the RadialDistance function from the semla v1.1.6 package. This metric is calculated from the border of the region of interest (ROI). Thus, radial distances are equal to 0 at the tumoural margin, and they become more positive or negative further from or closer to the tumour core, respectively.

We tested the correlation between functional or drug sensitivity BCS and the radial distance to the tumour core in each sequenced region. Negative Pearson correlation coefficients indicate increased functional enrichment or drug sensitivity at the proximity of the tumour core. In contrast, positive correlation coefficients reflect increased functional enrichment or drug sensitivity at the tumour periphery. The P-values were adjusted using false discovery rate (FDR) correction for multiple testing.

Functional enrichment analysis

To compare spots from different patients, we performed a differential expression analysis with the log-normalized counts of the groups of interest (major TCs, tumour TCs or expression clusters from patient V19L29). We compared each group against the rest using the FindMarkers function (min.pct = 0, logfc.threshold = 0) from the Seurat v4.3.0.1 package, which implements a non-parametric Wilcoxon rank sum test. This function returns the log fold-change of the average expression of each gene between the two groups, which was used to rank all genes in the expression matrix. For each comparison, we performed a pre-ranked gene set enrichment analysis (GSEA) with the fgsea function (minSize = 15, maxSize = 500, nPermSimple = 10 000) from the fgsea v1.26.0 package (36) using the Hallmark and Reactome v2023.1 collections (37–39), the cancerSEA gene sets (40) and other functional signatures (41–47).

Differential sensitivity analysis

To determine the drugs that specifically target the groups of interest (non-mixed major TCs or tumour TCs), we compared each group against the rest using the bcRanks function from the beyondcell v2.2.0 package. For each comparison, this function ranks the drugs in the SSc breast according to the residual’s mean and the switch point. The switch point is a signature-specific Beyondcell metric that reflects the directionality and homogeneity of the drug response throughout the experiment. An SP = 0 indicates a homogeneously sensitive response to the drug, whereas an SP = 1 indicates a homogeneously resistant phenotype. Intermediate SPs represent a heterogeneous response.

We kept the drugs with intermediate switch points (between 0.4 and 0.6) and lowest or highest residual’s mean (5th percentile and 95th percentile for major TCs; 1st percentile and 99th percentile for tumour TCs). These drugs display a heterogeneous sensitivity pattern across the experiment and are the most specifically effective or ineffective against the group of interest.

Cell-cell communication analysis

We inferred spatially proximal cell–cell ligand-receptor interactions using CellChat v2.1.0 (48). We created a spatial CellChat object with the SCTransform-normalized expression counts, metadata and coordinates from all patients, grouping the tumour spots by TCs and the TME spots by cell type. Then, we performed a cell–cell communication analysis using the Secreted Signaling, ECM-Receptor and Cell–Cell Contact databases with the computeCommunProb function (type=‘truncatedMean’, trim = 0.1, distance.use = FALSE, interaction.range = 250, contact.knn = TRUE, contact.knn.k = 6, population.size = TRUE). We filtered out results for groups with <10 spots and retrieved the significant (thres = 0.05) interaction strengths for ligand-receptor pairs of interesting signalling pathways [signaling = c(‘PD-L1’, ‘IGF’, ‘DESMOSOME’)] and cell groups [sources.use = targets.use = c(‘TC1.1″, ‘TC1.2″, ‘TC2’, ‘TC3’, ‘T cells’, ‘CAFs’, ‘B cells’)] using the subsetCommunication function. Finally, we scaled the interaction strengths of each ligand-receptor pair between 0 and 1.

Tumoural ROI definition

Using SCTransform-normalized expression counts and FindNeighbors(dims = 1:30), we applied the FindClusters function from the Seurat v4.3.0.1 package five times, progressively increasing the resolution settings (0.1, 0.25, 0.5, 0.75 and 1). To visualize the stability and relationships between clusters across these different resolutions, we employed the clustree v0.5.1 package (49), as illustrated in Supplementary Figure S5C. Since we aimed to identify clusters corresponding to tumoural ROIs, we coloured the nodes based on the proportion of tumour spots within each cluster. We selected a resolution of 0.25 as the most suitable, as it identified seven clusters with high tumour purity that remained stable across increasing resolutions.

All expression clusters but number 7 were confined within unique spatial regions. Thus, we subdivided cluster 7 into two subclusters using the FindSubCluster(resolution = 0.075) function. We defined eight initial ROIs from the expression clusters overlapping globular tumour ducts (3, 4, 5, 6, 7_1, 7_2, 8 and 9). Next, we refined each ROI by splitting spatially disconnected spots using the DisconnectRegions function from the semla v1.1.6 package and keeping the main group of connected spots as ROI. Finally, we selected ROIs with at least 100 spots from any subclone.

For each tumoural ROI, we defined the edge and inner regions with the RegionNeighbors function (mode=‘inner’) from the semla v1.1.6 package. The spots detected at the inner border of each ROI were labelled as ‘edge’, whereas the rest were categorized as ‘inner’. Then, we tested the differential TC proportion between the edge and inner regions for each ROI and subclone using Fisher’s exact test. The P-values were adjusted using FDR correction for multiple testing.

Statistical analysis and visualization

Statistical analysis, data manipulation and visualization were carried out using the R v4 programming language (R Core Team, 2023; https://www.r-project.org/) and the packages specified in Supplementary Table S2.

Results

Tumour and TME dissection in ST breast cancer samples

In order to spatially dissect the therapeutic and functional heterogeneity in breast cancer samples, we analysed Visium ST data retrieved from public repositories (5, https://www.10xgenomics.com/resources/datasets [last accessed 16 Jun 2023]) coming from nine breast cancer patients stratified into luminal (n = 2), HER2+ (n = 3) and TNBC (n = 4) subtypes. Of these nine patients, seven had histopathological annotations, and sample V19L29 (HER2+) contained two consecutive slides. In addition, two TNBC patients (CID4465 and CID44971) and one ER+ patient (CID4535) had two or three regions sequenced into the same slide.

To classify each spot as either tumour or TME, we used four independent annotation sources with high overlap (Figure 1A and B): the pathologist annotations (Figure 1C), the ESTIMATE score, which is inversely proportional to the tumour purity (Figure 1D), the proportion of cancer cells per spot (Figure 1E) and the clonal composition of each sample based on CNA profiles (Figure 1F). In total, we successfully annotated 18 207 malignant and 11 858 non-malignant spots (Figure 1G and Supplementary Figure S1A) coming from all nine patients. Notably, deconvolution analysis revealed that epithelial cancer cells, CAFs and myeloid cells were present in all samples (Supplementary Figure S1B–D).

Spatial compartmentalization of drug response in the tumour ecosystem

To predict the drug sensitivity of each spot, we created the SSc breast (see the ‘Materials and methods’ section), which contains 1372 gene signatures that reflect the transcriptional change between treatment-naive sensitive and resistant breast cancer cell lines to >1200 compounds. With Beyondcell (18), we computed a matrix of enrichment scores for each spot and signature that was subsequently used to obtain 16 TCs with different predicted drug sensitivities (Figure 2A).

We defined three major groups of TCs according to the observed proportion of tumour spots: tumour-rich (tumour proportion ≥0.95), mixed (tumour proportion >0.65 and <0.80) and TME-rich (tumour proportion ≤0.15) (Figure 2A and B). Interestingly, in mixed TCs, tumour and TME spots were grouped together, suggesting that therapeutic heterogeneity is not entirely driven by cell type identity (Supplementary Table S3). We confirmed that these three major TCs were present in all patients, implying similar drug responses independent of the breast cancer subtype (Figure 2C and Supplementary Table S4).

Next, we studied the spatial organization of these major TCs. We observed that mixed TCs tended to co-localize with both tumour-rich and TME-rich TCs. Moreover, we noticed a clear repellent effect between spots in tumour-rich and TME-rich TCs (Figure 2D–E). These results indicate that the drug responsiveness defined by these major TCs is spatially organized and might overlap with the main compartments of the tumour ecosystem: the tumour and the TME regions and the interphase that physically separates them.

To assess whether drug response patterns were conserved across patients, we generated a correlation matrix from the BCS within the tumour and TME regions of each patient (Figure 2F). Interestingly, these two spatial regions displayed opposite drug response patterns within the same patient, as indicated by their high anticorrelation, revealing the existence of intratumour therapeutic heterogeneity. At the same time, we identified two correlation clusters that perfectly matched our tumour and TME-rich annotations, suggesting recurrent response patterns within the same region across the whole cohort. Nevertheless, within the same cluster, the correlations that involved TNBC samples were lower than the rest, implying a higher interpatient therapeutic heterogeneity in TNBC samples.

Altogether, these results point to the existence of therapeutic heterogeneity between the spatial compartments of the tumour ecosystem. These different drug response patterns are conserved in all breast cancer patients independently of the subtype. Nevertheless, interpatient therapeutic heterogeneity exists and is particularly evident in TNBC tumours, in concordance with previous knowledge (50).

A functional and drug sensitivity gradient exists between the tumour and TME compartments

In order to functionally characterize the major TCs, we first deconvoluted their cellular composition and confirmed that the tumour-rich compartment was almost exclusively formed by cancer cells, together with a minority of CAFs and myeloid cells (Figure 3A). The interphase also contained endothelial cells, while the TME region displayed the highest cellular diversity. Overall, T cells appeared in the interphase and TME compartments, but interestingly, patient CID4465 (luminal) showed lymphocyte infiltration in the tumour region.

We applied Beyondcell to compute enrichment scores for expression biomarkers (Figure 3B), cancer-related drugs like doxorubicin (Figure 3C) and well-known tumoural biological processes such as cell proliferation (Figure 3D). Taking advantage of the spatial information associated with each spot, we also calculated the radial distance to the core of the tumour-rich compartment (Figure 3E). Positive radial distances correspond to spots far away from the tumour core, whereas negative radial distances indicate proximity to the centre of the tumoural mass, and a distance of 0 denotes the tumoural margin. We then computed the region-wise correlation between radial distances and functional or sensitivity BCS (Supplementary Tables S5–S7), discovering a functional and therapeutic gradient from the tumour core to the periphery of the tumoural regions (Figure 3F).

Notably, the spots within the tumour region were enriched in a ductal-like phenotype and predicted to be sensitive to standard hormonal therapies and HER2 inhibitors such as tamoxifen, afatinib and allitinib. We also observed an enrichment in proliferation and stemness functions closer to the tumour core. Accordingly, this region was more sensitive to cell cycle arrest agents (BI-2536 and KW-2449), telomerase inhibitors (MST-312), PI3K/Akt/mTOR inhibitors (piperlongumine, PI-103, PF-4 708 671) and WNT inhibitors (LGK974). Moreover, the tumour core was more hypoxic and paradoxically more sensitive to VEGFR inhibitors (BMS-690 514, nintedanib, foretinib and sunitinib) since hypoxia triggers the expression of VEGF, a key regulator of angiogenesis. Interestingly, we observed a decreased sensitivity to DNA-related agents and olaparib within the tumour region, despite displaying a higher enrichment in DNA damage and BRCAness. This lower sensitivity can be due to an enrichment in DNA repair functions, which may help cancer cells repair DNA damage efficiently, contributing to treatment resistance. Indeed, upregulation of DNA repair mechanisms and hypoxia have been linked to chemoresistance in breast cancer (51,52). On the other hand, the tumour periphery was more enriched in an invasive phenotype, inflammation, ECM stiffness and epithelial-to-mesenchymal transition (EMT). Consistently, these spots displayed increased sensitivity to JAK-STAT signalling inhibitors, statins and immunosuppressants or anti-inflammatory agents, drugs that mainly target the immune component and the inflammatory response of the TME.

Compartment-wise functional analysis and drug ranking confirmed these functional and therapeutic differences (Supplementary Figure S2A and B, and Supplementary Tables S7–S9). We assessed that the interphase displayed medium sensitivities to all these drugs [(mean BCS and standard deviation for Tumour-specific therapies: Tumour-rich = 7.74 (7.67), TME-rich = −7.27 (7.64), Interphase = −0.52 (4.74); mean BCS and standard deviation for TME-specific therapies: Tumour-rich = −7.92 (7.25), TME-rich = 7.44 (6.70), Interphase = 0.53 (5.18)], suggesting the existence of a therapeutic continuum from the tumour core to the TME and vice versa.

Using Beyondcell, we checked these predicted trends in three examples of breast cancer clinical subtypes. As expected, only a normal epithelial tissue patch in region A of the TNBC sample expressed ESR1 (ER gene), PGR (PR gene) and ERBB2 (HER2 gene) biomarkers (Figure 3B). The triple-negative tumour-rich region displayed specific sensitivity to first-line treatment doxorubicin (Figure 3C), which can be explained by an enrichment in cell cycle activity in the same region (Figure 3D). Conversely, in luminal and HER2+ samples, the high expression of ESR1 and ERBB2 biomarkers in the tumour-rich region was correlated with the predicted response to standard-of-care therapies like tamoxifen or HER2 inhibitors (Supplementary Figure S3A–C). Furthermore, for each TME spot, we computed a single sample bidirectional enrichment score (53) depicting its protumour (>0) or antitumour (<0) status, according to the BCS ranking of 22 pan-cancer microenvironment signatures (54). We observed TME heterogeneity among samples (Supplementary Figure S3D). The TME of the luminal example was predominantly antitumoural and thus did not need chemical inhibition. In contrast, the HER2+ patient displayed a dual TME whose protumoural region could also be targeted with the HER2 inhibitor afatinib.

The interaction with the TME influences cancer drug response

To determine whether there was therapeutic heterogeneity within the tumour region, we reclustered the spots labelled as tumour with Beyondcell. From this unsupervised clustering, we obtained five TCs (Figure 4A) and, to determine whether these drug response groups were conserved across patients, we once again generated a correlation matrix from the BCS within the tumour TCs of each patient (Figure 4B, see the ‘Materials and methods’ section). TC4 and TC5 were patient-specific and displayed high correlations with TC3. Since these two clusters likely resulted from over-clustering of TC3, we regrouped them into TC3. On the other hand, we noticed that TC1 could be subdivided into two correlation subclusters. After this refinement process (Figure 4C), we obtained three TCs (TC1, TC2 and TC3), one divided into two subclusters (TC1.1 and TC1.2), that were conserved across patients (Supplementary Figure S4A and Supplementary Table S10).

We verified that these TCs were different in terms of proximity to the tumour core and compartment affiliation (Supplementary Figure S4B and C). TC2 was the closest to the centre of the tumoural mass and displayed the highest percentage of tumour-rich spots (TC1.1 = 68%, TC1.2 = 33%, TC2 = 90% and TC3 = 4%). In contrast, TC3 was the more distal cluster, mainly composed of interphase spots (TC1.1 = 25%, TC1.2 = 66%, TC2 = 10% and TC3 = 84%). Thus, we hypothesized that TC2 constituted the tumour core, whereas TC3 comprised the tumoural margin.

Functional enrichment analysis (Figure 4D, and Supplementary Tables S7 and S11) revealed that TC2 was the least invasive TC (e.g. downregulation of invasion gene set or collagen formation). Furthermore, ligand-receptor analysis (Supplementary Figure S4D) confirmed that TC2 spots were in close contact via DSC2/DSG2 interactions. These cadherins are primary constituents of the desmosome, a cell–cell junction that maintains the integrity of normal epithelia. Thus, these results convey that the tumour core represented by TC2 is in a pre-invasive state. The remaining TCs were enriched in functions that suggested communication with the TME. TC3 was characterized by B cell activation, phagocytosis and triggering of MAPK signalling. TC1.1 was defined by the triggering of the complement cascade, cytokine signalling (TNFα and IFNα, β and γ) and functions that indicated a T-cell suppressive state via PD-1 signalling. Conversely, TC1.2 displayed an increased translational activity triggered in response to starvation, as well as upregulation of Nonsense-mediated mRNA decay (NMD), Slit/Robo signalling and collagen formation, which might be related to ECM reorganization, contributing to its invasive phenotype. The upregulation of these functions points to nutrient deprivation due to inadequate blood supply, which may favour a hypoxic niche. Some authors have proposed intron retention, an aberrant splicing event which can introduce premature stop codons and lead to NMD, as a possible mechanism by which hypoxia can promote angiogenesis, cancer cell migration and invasiveness in breast cancer (55). Ligand-receptor analysis (Supplementary Figure S4D) confirmed that TC1.1 interacted with suppressed T cells via PD-1/PD-L1 binding and TC1.2 with CAFs via IGF signalling, which has been related to invasiveness, EMT and resistance to EGFR inhibitors such as gefitinib and chemotherapies like gemcitabine and paclitaxel in breast cancer (56). Collectively, these results imply that each TC exhibits distinct functional states conditioned by interactions with different components of the TME.

Using the SCSubtype gene signatures from Wu et al. (5), we classified the tumour spots into intrinsic subtypes (Luminal A, Luminal B, HER2+ or Basal) based on the maximum BCS, with ties resulting in an ‘Undetermined’ annotation. Next, we identified drugs specifically targeting each tumour TC, but no TC was enriched in a particular intrinsic subtype (Figure 4E and Supplementary Table S12). The tumour core depicted by TC2 displayed sensitivity to non-apoptotic cell death agents and kinase inhibitors, whereas TC3 appeared to be differentially sensitive to MAPK inhibitors. Interestingly, both TC1 subclusters displayed lower sensitivity to their specific drugs compared with the other clusters [mean BCS and standard deviation for specific therapies: TC1.1 = 2.35 (6.81), TC1.2 = 1.08 (4.26), TC2 = 6.15 (5.77), TC3 = 7.95 (4.66)]. Still, cancer cells in an immunosuppressive microenvironment defined by PD-1/PD-L1 binding or interacting with CAFs via IGF signalling exhibited specific sensitivity to cell cycle arrest agents and PI3K/AKT/mTOR inhibitors. Our analysis also indicated that TC1.1 and TC1.2 were sensitive to SKI-II, a sphingosine kinase (SK) inhibitor and simvastatin. Upregulation of sphingosine kinase 1 (SK1) has been associated with invasiveness and chemoresistance in breast cancers, and targeting SK1 could have therapeutic potential (57). Moreover, simvastatin has been proposed to attenuate breast cancer metastases and recurrence (58,59). Of note, none of these subclusters was predicted to be sensitive to the standard-of-care tamoxifen, whereas the rest of the TCs displayed high sensitivity to this drug. This result highlights the relevance of studying ITH to identify pre-resistant tumour subpopulations that, if overlooked, can lead to treatment failure and relapse. We hypothesize that the interacting TME influenced the low drug sensitivity observed in TC1. We confirmed that TC1.1 was significantly enriched in the tumour inflammation signature gene set (44), which reflects a suppressed adaptive immune response and has been associated with sensitivity to anti-PD-1 therapies (44,60) (Supplementary Figure S4E). On the other hand, TC1.2 was significantly enriched in ECM stiffness compared to the rest of the TCs (Supplementary Figure S4F). ECM stiffness is linked to CAF activity, constitutes a mechanical barrier for treatment diffusion and promotes treatment resistance by triggering different signalling pathways in tumour cells (61).

Altogether, these results support the idea that intratumoural therapeutic heterogeneity does exist in breast cancer, and some niches may display decreased sensitivity, or even insensitivity, to standard-of-care treatments. These differences may be related to distinct transcriptomics states triggered by the interaction with specific components of the TME, which can sensitize to immunotherapies or protect cancer cells from drugs.

Subclones display heterogeneous responses depending on their location in tumour ROIs

We asked ourselves whether different subclones drove therapeutic heterogeneity within the tumour region. In order to explore this hypothesis, we selected patient V19L29 (HER2+) because it was the sample with more sequenced tumour spots (n = 4713; Supplementary Figure S1A), granting higher statistical power. We defined six tumoural ROIs based on well-delimited globular ducts (Figure 5A and Supplementary Figure S5A). We confirmed that ROIs had a biological identity, as they were identified as different expression clusters (Supplementary Figure S5B and C) enriched in different functions (FDR < 0.25), such as invasion in ROI5 or hypoxia in ROI4 among others (Supplementary Figure S5D, and Supplementary Tables S7 and S13). These ROIs contained six cancer subclones defined by CNA profiles (Figure 5B), grouped into TC1.1, TC2 and TC3 (Figure 5C). Subclones showed ROI-specificity (i.e. were localized in one or few ROIs), but the same subclone belonged to different TCs, suggesting intra-subclonal therapeutic heterogeneity (Figure 5D). We asserted that, for some subclones and ROIs, the proportion of TCs was significantly different (FDR < 0.05) between the edge and inner regions (Figure 5E), TC2 being more abundant in the inner region and TC1.1 and TC3 overrepresented at the edge. This result implies that the spatial location of spots, even with the same genetic background, conditions their transcriptomic profile and might ultimately result in different response patterns, highlighting the importance of studying ITH in its spatial context. When analysing all patients, we confirmed that tumour TCs consisted of mixtures of subclones sharing similar interactions with the TME (Supplementary Figure S5E), indicating that TME interactions may shape varying transcriptional profiles affecting drug sensitivity. Therefore, we propose that therapeutic heterogeneity is more likely driven by convergent transcriptional phenotypes shaped by specific TME interactions rather than distinct mutations defining subclones.

In the clinic, HER2+ patients are typically classified based on average ERBB2 overexpression or amplification. However, some authors have described molecular subtype heterogeneity between cancer cells within the same tumour. Notably, subclones 1 and 6 within ROI5 were diploid for ERBB2 and exhibited lower ERBB2 expression than the rest of the tissue (Figure 5F and G). HER2+ patients typically receive trastuzumab as a first-line treatment. However, based on a publicly available gene signature (62), Beyondcell predicted that none of the tumour spots of patient V19L29 would respond to this therapy (Supplementary Figure S5F). For patients unresponsive to trastuzumab, current guidelines recommend tyrosine kinase inhibitors (63). Beyondcell reported insensitivity to the HER2 inhibitor varlitinib in ROI5, whereas the rest of the tumour was predicted to be responsive (Figure 5H). If patient V19L29 were given HER2 inhibitors as monotherapy, ROI5 would survive to regenerate a HER2 inhibitor-resistant tumour. Beyondcell identified doxorubicin as a complementary therapy specifically targeting ROI5 (Figure 5I). Thus, we propose a combination of varlitinib and doxorubicin to target the major subpopulation of HER2+ cancer cells and the minor HER2− subclones to eliminate the tumour and avoid relapse.

We present another relevant case study of a TNBC patient in Supplementary Figure S6. For TNBC patients, anthracycline and taxane-based therapies are the only proven adjuvant systemic treatments preventing recurrence and improving survival (64). Patient CID44971 exhibited differential predicted responses to the taxane docetaxel, likely due to subclonal populations with distinct regional distribution within the primary tumour. Specifically, subclones 2 and 4 in region C showed decreased sensitivity to docetaxel but increased sensitivity to oxaliplatin. This platinum-based chemotherapy has proven beneficial in TNBC patients resistant to taxanes or anthracyclines (65). These two case studies illustrate pre-existing treatment resistance and highlight the value of using ST data to dissect ITH for therapy design.

Discussion

Personalized precision oncology faces a significant challenge in addressing ITH to boost treatment efficacy and avoid relapse (8,66). Extensive research has focused on characterizing cancer subpopulations and identifying drug vulnerabilities across samples (18,67–69). However, the impact of the TME on the tumour ecosystem has been usually overlooked. Current research indicates that the interactions with the TME are critical for tumour progression, favouring vascularization, ECM remodelling, immune exclusion and suppression (70). The selective pressure exerted by the local TME leads to the diversification of malignant and non-malignant cell subpopulations within the same tumour, thus enhancing ITH and increasing the chances of therapeutic resistance (71). Moreover, the TME components exhibit some conservation across breast cancer subtypes, providing a comprehensive therapeutic target (72). Sequencing-based ST provides new opportunities to study tissue organization and cancer–TME interactions (73). Leveraging ST and Beyondcell (18), this study identifies therapeutic niches based on drug sensitivity predictions across over 1200 drugs. Results reveal compartmentalization of drug response in breast cancer, unveiling sensitivity and therapeutic gradients relative to the distance from the tumour core. Moreover, we dissected therapeutic heterogeneity within the tumour region, which appears to be related to different TME interactions conserved across patients. Notably, we observed different response patterns within the same cancer subclone depending on its spatial location, underscoring the importance of considering spatial context in therapeutic decision-making.

We reported enrichment in cell proliferation, stemness and hypoxia near the tumour core, whereas inflammation, invasiveness, and EMT were enriched in the periphery. This functional gradient could result from the co-existence of diverse molecular gradients, including oxygen, nutrients, pH, chemokines and growth factors (74–77). They arise from the leaky vasculature of solid tumours and the rapid cancer cell proliferation. Anaerobic glycolysis is a cellular adaptation to hypoxia and low nutrient concentrations in regions far from blood vessels. In this metabolic process, pyruvate is converted to lactate and released to the extracellular space, acidifying the medium and generating a pH gradient. The secretion of growth factors and chemokines by TME cells generates additional biochemical gradients that, together with hypoxia and an acidic medium, favour cell migration and EMT (78). Moreover, these gradients also influence drug susceptibility. Hypoxia has been associated with increased DNA repair activity and subsequently reduced efficacy of multiple chemotherapeutic agents. Moreover, hypoxia induces the expression of ABC transporters that decrease the intracellular concentration of chemotherapeutic agents, favouring resistance (79). In addition, an acidic microenvironment can ionize weakly basic chemotherapies such as mitoxantrone, reducing their cellular uptake and cytotoxicity (80,81). Accordingly, we observed an enrichment in DNA repair activity in the hypoxic tumour core and a decreased sensitivity to different chemotherapies, including mitoxantrone. This hypoxic state could also explain the increased sensitivity to VEGFR inhibitors near the centre of the tumour. On the other hand, the enrichment in inflammation at the periphery of the tumour is coherent with an increased sensitivity to immunosuppressants and anti-inflammatory agents. These sensitivity and functional gradients defined three therapeutic regions that overlapped with the three main compartments of the tumour ecosystem: the tumour, the TME, and the interphase that physically separated the previous two.

We also dissected the therapeutic heterogeneity of spots labelled as tumour and characterized the resultant TCs, which were conserved across patients, linking function and ligand-receptor interactions to drug response within the tumour compartment. The tumour core, characterized by a non-invasive phenotype and active metabolism, was predicted to be sensitive to kinase inhibitors and non-apoptotic cell death agents. In contrast, cancer cells at the region’s periphery and close to B lymphocytes were sensitive to MAPK inhibitors. The remaining TC was subdivided into subclusters with different biological functions and TME interactions. TC1.1 cancer cells inhabited an immunosuppressive microenvironment promoted by PD-1/PD-L1 interactions, suggesting sensitivity to immune checkpoint inhibitors (70) that we confirmed. TC1.2 displayed increased collagen production and an invasive phenotype, potentially due to interactions with CAFs via IGF signalling. This interaction may contribute to ECM remodelling, which can reduce the sensitivity to targeted therapies by (i) creating a protective niche that physically impedes the delivery of drugs, (ii) mechano-transduction of proliferation or anti-apoptotic signals mediated by integrins, or (iii) metabolic changes induced in cancer cells due to the internalization of ECM components. These mechanisms have been observed in breast cancer in vivo and in vitro models, which displayed increased resistance to HER2 and multi-kinase inhibitors (82). Accordingly, TC1.2 was significantly enriched in ECM stiffness, which may explain the overall low sensitivity of this subcluster to the entire collection of tested drugs. Our findings, in line with previous studies (83,84), suggest that therapeutic heterogeneity within a single tumour is driven by distinct TME interactions that drive conserved biological functions across patients. Arora et al. identified the leading edge of oral squamous cell carcinoma as invasive, neighbouring CD8 + T cells and CAFs, while the tumour core retained an epithelial phenotype. In our study, TC2 mirrors the core, while TC1’s subclusters resemble the invasive leading edge, which was conserved across different tumour types.

Finally, based on the results of a single patient, we propose that the TME also modulates the expression of cancer cells with identical genetic backgrounds. Although we would need a higher sample size to extract robust conclusions, therapeutic heterogeneity in single-cell-derived clones has been previously linked to different cell states before drug administration (85). Other authors have noted that spatially distinct cancer cell states are not primarily driven by distinct subclones but rather by different TME interactions (83). Thus, studying the subclonal composition alone without considering the interacting TME may not be enough to dissect and understand therapeutic heterogeneity. Patient outcomes are likely driven by the least sensitive subclone (86). Thus, to address therapeutic heterogeneity and avoid relapse, some researchers have proposed adaptive therapy, which allows sensitive cells to survive, compete and eliminate resistant subpopulations (87), or clonetherapy, which targets both major cancer subpopulations and minor pre-resistant cells using combination therapy (88). While inferring clinical trajectories remains challenging, our findings suggest that therapeutic heterogeneity is lower than ITH. Focusing on heterogeneity from a therapeutic perspective could improve our understanding of cancer evolution, help predict resistance, and guide more effective treatments.

This analysis highlights several areas for improvement and future considerations. Spots in 10x Visium technology correspond to groups of one to ten cells, necessitating scRNA-seq references for cell type deconvolution and spot labelling according to arbitrary proportion thresholds. Moreover, these spots are discrete, thus leaving a large portion of the tissue unmeasured. Here, we categorized spots as tumour or TME based on the agreement of four highly overlapping annotation sources. Although tumour spots contained minor proportions of non-cancer cells, we detected ITH and therapeutic heterogeneity using tools designed for single-cell technologies. Nevertheless, leveraging new spatial sequencing technologies that allow single-cell resolution (89,90) would enhance the precision of these methodologies. Moreover, the spatial organization of drug response that we observed encourages the use of spatial-informed clustering methods (91–93) in future analyses. Additionally, integrating data from adjacent (horizontal) and serial (vertical) tissue sections may provide a comprehensive understanding of 3D tumour organization and therapeutic differences between the tumour’s core and periphery. We expect these 3D models to aid in studying how complex microanatomical structures, like blood vessels or lymph nodes, may influence drug sensitivity. Still, observing spatial organization solely through expression clustering implies that gene expression alone can explain most spatial variability. Thus, a potential application of these ST analyses is to extract geographic expression signatures that can be remapped onto scRNA-seq data to clarify expression differences between clusters of the same cell type. Furthermore, ITH is a multi-level phenomenon that includes genomics, epigenomics, transcriptomics, proteomics and the interaction with the TME. Multi-omics analysis facilitated by tools like Seurat v5 (94) would further our understanding of ITH and its relationship with therapeutic response. Finally, our collection of drug signatures derives from screenings in cancer cell lines with no tissue context and we lack specific drugs to target the TME. Incorporating drug signatures derived from screenings in co-cultures could elucidate the TME’s role in drug action. On top of that, we must recognize the importance of non-cellular constraints like ECM stiffness in understanding treatment diffusion and efficacy. Thus, future efforts to characterize the effect of non-cellular TME components on cellular function and drug response would benefit therapeutic heterogeneity dissection in breast and other cancer types.

Overall, our results provide valuable insights into the complex dynamics of therapy response, which depends on the genetic background and transcriptional state of cancer cells, their position in the tissue and their interactions with the TME. This study showcases the potential of integrating ST data into clinical decision-making to guide personalized therapy selection based on the tumour molecular characteristics and microenvironmental interactions.

Data availability

We downloaded luminal and TNBC Visium ST data (CID4290, CID4535, 1142243F, 1160920F, CID4465 and CID44971) (5) from Zenodo (DOI: 10.5281/zenodo.4739739). HER2+ Visium ST data (738811QB, 1168993F and V19L29) were downloaded from the 10x Genomics Datasets portal [https://www.10xgenomics.com/resources/datasets (last accessed on 16 Jun 2023)]. All these data included raw counts and matching tissue images. Wu et al. dataset also included metadata information with pathologist annotations for each spot. 10x Genomics provided an image with pathologist annotations for patient 738811QB, and we manually labelled each spot using Loupe Browser v5.0.1 (10x Genomics, 2021; https://www.10xgenomics.com/products/visium-analysis#loupe-browser). Additionally, we downloaded scRNA-seq data (5) from the Single Cell Portal (ID: SCP1039) to construct the spot deconvolution reference.

The SSc breast and Beyondcell objects generated in this work and the code used to create and analyse them can be found on Zenodo (DOI: 10.5281/zenodo.14247036). The beyondcell (https://github.com/cnio-bu/beyondcell) and ggseabubble (https://github.com/mj-jimenez/ggseabubble) packages are available on GitHub.

Supplementary data

Supplementary Data are available at NAR Cancer Online.

Acknowledgements

We thank all CNIO Bioinformatics Unit members for their support. We also thank Coral Fustero-Torre for her literature suggestions and fruitful discussions about the state of the art.

Author contributions: M.J.J.-S. designed the study, performed the analyses and wrote the manuscript under the supervision of G.G.-L. and F.A. S.G.-M. computed the SSc breast and reviewed the manuscript. M.R.-F. analysed patient V19L29 independently, providing insights that influenced the principal analysis. Beyondcell was updated and tested by M.J.J.-S. and maintained by M.J.J.-S. and S.G.-M. All authors read and approved the final paper.

Funding

Agencia Estatal de Investigación, Instituto de Salud Carlos III (ISCIII), European Regional Development Fund (ERDF), “A way of making Europe' [IMP/00019]; Ministerio de Ciencia e Innovación / Agencia Estatal de Investigación (10.13039/501100011033), European Union, European Regional Development Fund (ERDF-EU), “A way of making Europe' [PID2021-124188NB-I00]; Comunidad de Madrid, European Structural and Investment Funds [P2022/BMD-7379]; Horizon Europe Framework Programme [101058427]; Fundación científica de la Asociación Española Contra el Cáncer [PRYCO234528VALI]; ‘la Caixa’ Foundation [HR23-00051]; ISCIII, European Union, the Recovery, Transformation and Resilience Plan (PRTR) through Next GenerationEU recovery funds (MRR) [PMP22/00064].

Conflict of interest statement. None declared.

References