Main

Cancers are complex and dynamic entities that are constantly reshaped by the interactions between neoplastic cells and their microenvironments4,5,6. Whole-genome sequencing (WGS) analysis of the average cancer detects thousands of somatic mutations and multiple genetically related yet distinct groups of cells termed ‘subclones’2,7,8. However, as genomic technologies typically assay DNA from dissociated tissues, the phenotypic consequences and the ecosystem pressures that are critical to fully understanding cancer evolution are lost9,10. Consequently, relatively little is currently known about the nature or causes of spatial patterns of cancer growth, phenotypic characteristics of distinct subclonal lineages or their interactions with tissue ecosystems11. Still, this information appears key because adverse cancer outcomes—growth, progression and recurrence—are properties of genetically distinct subclones3,12,13,14,15,16.

Lineage tracing using somatic mutations is a powerful tool for inferring the ancestral relationships between cancer subclones, but methods to perform this in preserved human tissue context are lacking3,14,17,18,19,20. Histology-driven sampling, such as laser capture microdissection (LCM)21, combined with low-input nucleic acid library sequencing or even single-cell sequencing goes some way towards resolving subclone spatial structure19. However, even the most exhaustive sampling strategy will struggle to provide an unbiased representation of the cancer clone territories, particularly across whole-tumour sections. Recently described spatial genomics approaches permit the de novo spatial detection of cancer clones with distinct copy number profiles, but this does not permit the detection of point mutations or quantitative read outs of intermixed clones22,23. It has previously been demonstrated that individual mutations can be detected in situ using in situ hybridization24 or mutation-specific padlock probes25,26,27,28. However, these approaches are limited by the number of available fluorophores. Given that every cancer and subclone therein is genetically unique, to reconstruct ancestral relationships in both space and time, we need to be able to trace multiple, cancer-specific somatic mutations simultaneously8.

To address this need, we developed a genetic clone mapping workflow that is centred around base-specific in situ sequencing (BaSISS) technology. We derived quantitative maps of multiple genetic clones in eight tissues from two multifocal breast cancers that span the main histological stages of early cancer progression: ductal carcinoma in situ (DCIS), invasive cancer and lymph node metastasis29. In a case of DCIS, clones exhibited co-existence and segregation patterns in different parts of the breast ductal anatomy. By integrating genetic clone maps with multimodal spatial data layers, we found that genetically similar regions can be scattered across wide areas yet maintain similar transcriptional and histological features and foster recurrent ecosystems. Finally, we found that genetic progression, which encapsulates the historical order of events, does not necessarily translate directly to transitions in histological state that are commonly assumed to reflect the stages of cancer progression, thus warranting a combined genetic and histological assessment of cancer evolution.

The BaSISS workflow

The BaSISS workflow is centred around fresh frozen tissue blocks that undergo serial cryosectioning to generate tissue for bulk WGS and z-stacked sections for in-tissue spatial clone mapping and spatial phenotyping (Fig. 1a). Following subclone detection from bulk WGS data, there are three core BaSISS steps. First, to facilitate detection of multiple clones of interest, BaSISS padlock probes with sequence-specific oligonucleotide target recognition arms are designed towards both mutant and wild-type alleles of clone-defining somatic variants. A unique 4–5 nucleotide reader barcode on each probe enables multiplexing27. BaSISS targets can take the form of any expressed somatic mutation, including point mutations and rearrangement breakpoints, and can be supplemented with copy number alterations (Supplementary Table 1). Second, BaSISS and transcript detection are performed as previously described for gene expression ISS using cyclical microscopy27,30 (Fig. 1a and Supplementary Methods).

Fig. 1: The BaSISS workflow to generate cancer clone maps.
figure 1

a, Following de novo mutation detection and subclone discovery in WGS data, the BaSISS workflow is performed as follows: (1) bespoke mutation-specific padlock probes are designed. (2) BaSISS transcripts are detected. To achieve this, BaSISS padlock probes hybridize to complementary DNA (cDNA) in situ. By virtue of a highly specific DNA ligase, only completely target-complementary padlock probes are ligated and form closed circles. Ligated probes are amplified through rolling circle amplification and their reader barcodes are detected in tissue space through sequencing by ligation with fluorophore-labelled interrogation probes and cyclical microscopy. (3) Mathematical modelling of BaSISS signals and the genotype of clones is then performed to derive clone maps. (4) Subsequent phenotype and microenvironment characterization of clones is then possible, by integrating clone fields with spatial datasets acquired from serial tissue sections. The BaSISS model and cell typing are described further in Extended Data Figs. 1 and 2. b, The two cases of multifocal primary breast cancer (PBC) used to develop the BaSISS approach. Coloured tiles report the histological features within each sample and the experiments performed. The number of clones identified by WGS and targeted by BaSISS are reported as white numerals. c, The traditional histological model of breast cancer progression. DCIS, ductal carcinoma in situ; H&E, haematoxylin and eosin; LN, lymph node; NST, invasive carcinoma of no special type; TME, tumour microenvironment.

Full size image

Third, continuous spatial subclone maps are generated using a statistical algorithm that exploits BaSISS signals as well as local cell counts (derived from the DAPI channel during the fluorescence microscopy of BaSISS) using two-dimensional Gaussian processes (Extended Data Fig. 1 and Supplementary Methods). The variational Bayesian model also accounts for unspecific or wrongly decoded signals and variable probe efficiency and is augmented by variant allele fractions in the bulk genomic sequencing data. In an optional, fourth characterization step, BaSISS clone maps can be aligned and integrated with additional layers of spatial phenotype data. In this study, we performed spatially resolved single-cell transcriptomics using targeted ISS (using a previously published 91 gene oncology, a novel 62 gene immune panel and drawing on published single-cell RNA sequencing data)30,31 and immunohistochemistry (IHC) staining (Extended Data Fig. 2a–c and Supplementary Methods). Additional sections were obtained to perform validation of our workflow using LCM and low-input WGS as previously described32 (Extended Data Fig. 3a).

Two cases of multifocal breast cancer

The cohort includes eight tissue blocks from two patients (P1 and P2) who underwent a surgical mastectomy for a multifocal breast cancer. These patients were selected to permit a comparison between genetic and histological progression models in early breast cancer development29 (Fig. 1b,c). P1 had two separate oestrogen receptor (ER)-positive, human epidermal growth factor receptor 2 (HER2)-negative primary invasive breast cancers (PBCs) within a 5-cm bed of DCIS; we used tissue blocks from both PBCs (samples P1-ER1 and P1-ER2) and three regions from DCIS (samples P1-D1, P1-D2 and P1-D3). P2 had two separate PBCs of the ‘triple-negative’ subtype (lacking the ER, progesterone receptor and HER2). We sampled both PBCs (samples P2-TN1 and P2-TN2) and an axillary lymph node that contained metastatic cancer deposits (sample P2-LN1) (Fig. 1b).

Accurate and reproducible maps of clones

To demonstrate that spatial BaSISS signal counts can provide a meaningful read out of the underlying somatic genotype, we first focused on three samples from P1 (P1-ER1, P1-ER2 and P1-D1) (Fig. 1b). Previous multiregional WGS experiments identified mutation clusters that equated to six phylogenetic tree branches, and these were present at different levels across the three samples3 (Fig. 2a,b, Extended Data Fig. 3b and Supplementary Table 3). To enable spatial detection of subclones, BaSISS padlock probes were designed towards 51 alleles that report on each branch of the phylogenetic tree: 25 single-base substitutions and the equivalent wild-type base, as well as an amplified oncogene (FGFR1) (black numbers; Fig. 2b and Supplementary Table 1). Subclones are referred to by a patient identifier and the colour of the corresponding node of the phylogenetic tree: P1-purple, P1-red, P1-grey, P1-orange, P1-green and P1-blue (Fig. 2b). A subclone genotype comprises the branch mutations accumulated as one moves from the tree root to the subclone node, therefore P1-green contains grey, blue and green branch mutations (Extended Data Fig. 3b). The bulk WGS-derived tree was corroborated by spatial co-occurrence of BaSISS signals and LCM–WGS validation data (Fig. 2c and Extended Data Figs. 3c,d and 4a).

Fig. 2: Converting BaSISS spatial signals into maps of clones.
figure 2

a, Bar plots of cancer cell fractions (CCFs) derived from bulk WGS of the P1 samples. b, Phylogenetic tree reconstructed from multiregional bulk WGS data from P1 (see Supplementary Methods for details). Each branch is labelled with the total number of WGS mutations defining the branch (grey text) and the number of BaSISS probes designed to target that branch (black text). c, Three heatmaps of variant allele fractions (VAFs) calculated using data derived from n = 11 regions of P1-ER1 and P1-ER2 (marked in d). Raw BaSISS VAFs (for each target mutation the number of mutant signals divided by total number of mutant plus wild-type signals) (top) and model-imputed BaSISS VAFs (middle) are derived from raw BaSISS signal data within these regions. In serial tissue cryosections, corresponding z-stack regions were identified and subjected to LCM–WGS. Resulting LCM–WGS VAFs are presented (bottom). Mean per-gene correlations are approximately 0.41 and 0.90 for BaSISS to LCM–WGS and model-imputed VAFs to LCM–WGS comparisons, respectively. Sample names are coloured according to the dominant BaSISS subclone in the sampled region. Each row represents a targeted mutation. The mutations plotted in d are labelled by their gene name; for PTEN there are two separate mutations. d, Spatial BaSISS detections of barcodes reporting on five selected mutations, coloured according to their targeted branch. White contours indicate LCM regions (relates to c). e, BaSISS clone maps in physical space projected on the DAPI image (nuclei are white), derived using BaSISS mathematical modelling of signals from 45 informative targets. Each clone has a different colour, and dominant clones are reported (shown if the CCF is more than 25% and the inferred local cell density is more than 300 cells per mm2). Scale bars, 2.5 mm (d,e).

Source data

Full size image

On average, 97% of detected BaSISS spot signals were converted into feasible barcodes33. The median target-specific coverage across 300 mm2 of breast tissue was 13,000-fold (Supplementary Table 4). BaSISS-derived variant allele fractions exhibited strong correlation across replicate experiments on serial tissue sections (R = 0.76–0.93, Pearson’s; Extended Data Fig. 4b), demonstrating quantitative reproducibility.

BaSISS signals coloured according to their subclonal mutation branch revealed a first, albeit noisy, visual glimpse of subclonal growth structure (Fig. 2d). Broad patterns were preserved in technical replicate experiments using adjacent tissue sections (Extended Data Fig. 4c). Although the number of signals detected per nucleus (n = 0.82) does not provide single-cell resolution of the somatic genotype, it is possible to aggregate information (1) spatially over areas of approximately 100 × 100 µm2, and (2) across alleles co-occurring in a particular subclone to infer the local clonal composition of different tumour clones and normal cells (Supplementary Table 4). This process generated detailed maps covering several squared centimetres of tumour tissues (Fig. 2e). Of note, the clone mapping algorithm also implicitly adjusts the observed allele frequencies for a range of systematic biases (Extended Data Fig. 4d), stemming from the use of RNA-derived signals, differential BaSISS probe sensitivity and allele confusion to produce highly consistent maps across replicates (Extended Data Fig. 4e). Although the raw BaSISS variant allele fractions of many probes were noisy owing to the aforementioned biases, the modelled allele frequencies were in highly accurate agreement with LCM–WGS validation data (Fig. 2c,g and Supplementary Table 5). This further corroborates the quantitative nature of BaSISS-derived clone maps that can be explored using an interactive web browser (https://www.cancerclonemaps.org/).

Charting histogenomic relationships

Histology-driven sampling of well-defined stages of cancer progression can uncover mechanisms and markers of disease progression10,19,29,34. Up to two-thirds of PBCs contain both invasive cancer and intermixed DCIS, a non-obligate precursor lesion. How these distinct ‘stages’ of cancer development might relate to genetic diversification within the same tissue is generally unknown35 (Fig. 1c). To demonstrate that BaSISS can chart these relationships across entire tissue sections, we examined three PBC samples with intermixed invasive and DCIS histology: P1-ER1, P1-ER2 and P2-TN1 (Fig. 3 and Extended Data Figs. 5 and 6a–c).

Fig. 3: Genetic clones mapped in histological context from three PBCs.
figure 3

a, BaSISS maps of two PBCs from P1 with intermixed DCIS and invasive cancer. The most prevalent genetic clone is projected as a coloured field (corresponds to b) on DAPI images (reported if the CCF is more than 25% and the inferred local cell density is more than 300 cells per mm2). Scale bar, 2.5 mm. Pie charts report the WGS-estimated clone composition of P1-ER1 and P1-ER2. Inset images (right) are regions of P1-ER2 (H&E-stained serial tissue sections) that represent three histological progression states. Scale bar, 250 µm. b, The phylogenetic tree was inferred from P1 multiregion WGS: branches are scaled according to and annotated with the number of WGS mutations and driver mutation-containing genes. Branches and nodes are coloured to reflect the clones mapped in a. Heatmaps report clone composition in 34 and 44 histologically annotated epithelial-containing microregions of P1-ER1 and P1-ER2, respectively. Microregions include individual ducts or randomly selected similarly sized regions of invasive cancer (see Extended Data Figs. 4b and 5b and the web browser https://www.cancerclonemaps.org/ for microregion details). HP, hyperplasia; N, normal ducts. c,d, IHC in P1-ER1 (c) and P1-ER2 (d) for the proliferative marker Ki-67 in six clone territories (indicated by contour colour); the percentage of nuclei staining positive (brown) is reported. Scale bars, 250 µm. e, As in a, but a clone map of P2-TN1. Scale bar, 2.5 mm. Mini-images report ISS-derived cell types (right) and H&E tissue section snapshots of the two cancer growth patterns (GP1 and GP2) reported in P2-TN1 (left). Scale bar, 250 μm. f, Phylogenetic tree for P2 and heatmap of 36 P2-TN1 microregions, as in b. Branches relating to clones not detected in this sample (that is, only found in P2-LN1) are shaded grey. The bottom heatmap is the estimate by the histopathologist and reports the contribution of different growth patterns to the microregion, defined by distinct nuclear and architectural features (Supplementary Methods).

Source data

Full size image

BaSISS detected 2–4 subclones per PBC in accordance with bulk WGS data. Clone maps (Fig. 3a,e) and the quantitative clonal composition of 73 individually annotated microregions (Fig. 3b,f and Extended Data Figs. 5a,b and 6a,b) revealed that individual subclones form spatial patterns that were, by varying degrees, related to the histological progression states. Normal tissue elements, including immune aggregates and histologically normal ducts, appear unstained consistent with a wild-type status for the targeted clones (green and yellow contours, respectively; Fig. 3a,e). In P1-ER2, an area of hyperplasia was predicted and confirmed by LCM–WGS to be genetically unrelated to the cancer (blue contour; Figs. 2c and 3a).

In each PBC, the genetic and histological progression models were broadly consistent, in which the invasive disease was mainly composed of cells from the most recently diverged subclone: P1-red, P1-purple and P2-purple in samples P1-ER1, P1-ER2 and P2-TN1, respectively (Fig. 3b,f). By contrast, earlier diverging clones colocalized entirely or in part to the histological pre-invasive lesion: DCIS. For example, in P1-ER2, BaSISS predicted that green branch mutations were completely absent from the invasive compartment, a conclusion that is supported by three separate microdissections (LCM–WGS) from distant regions of invasive cancer in P1-ER2 (Fig. 2c and Extended Data Fig. 5c).

However, in each PBC, there was a subclone that spanned both DCIS and invasive histology, revealing that disconnects between histological and genetic progression states can exist. This was the case for clone P1-red in P1-ER1 and clone P1-purple in P1-ER2. These DCIS-invasive spanning clones could be distinguished from each other by hundreds of private mutations, including different inactivating driver mutations in PTEN, indicating parallel evolution along these divergent lineages that resulted in two distinct instances of cancer invasion (total mutation numbers label the phylogenetic tree branches; Fig. 3b). The spatial predictions of the BaSISS model of intraductal acquisition of PTEN mutations and PTEN protein loss was confirmed by LCM–WGS and IHC, respectively (Fig. 2c and Extended Data Fig. 5d). In sample P2-TN1, the only predicted driver point mutation was a deleterious mutation in the tumour suppressor gene TP53, and this was detected in both DCIS and invasive compartments and was also present in all cancer regions of the second PBC, P2-TN2, consistent with an early onset in the development of this cancer (phylogenetic tree; Fig. 3e,f). These data therefore suggest that many, if not all, of the genetic events necessary to initiate the invasive transition in these three cancers were acquired within the ducts, and subsequently both intraductal expansion and stromal invasion ensued.

Phenotypic changes accompany progression

Next, by integrating additional layers of spatial data, we sought to establish how phenotypic changes relate to genetic-state and histological-state transitions. In P1-ER1 and P1-ER2, consistent with a more proliferative phenotype, PTEN-mutant clone regions exhibited denser Ki-67 IHC nuclear staining, than PTEN wild-type ancestral clone regions (false-discovery rate (FDR) = 0.004 P1-red versus P1-orange; and FDR = 0.03 P1-purple versus P1-green) (Fig. 3c,d and Extended Data Fig. 5e). However, for a given genetic clone, the Ki-67 score was similar irrespective of whether it occupied a DCIS or invasive state, indicating that upregulation of Ki-67 is temporally related to acquisition of a PTEN mutation and precedes invasion.

By contrast, cellular resolution spatial transcriptomics analysis of P1-ER2 revealed that epithelial cell expression of several genes—CLDN4 (encoding claudin 4), ACTB (encoding β-actin), KRT5 (encoding keratin 5) and CTSL2 (encoding lysosomal cysteine protease cathepsin V)—differed between DCIS and invasive compartments occupied by the same, P1-purple, clone (Extended Data Fig. 5f). These transcriptional changes might therefore be considered more closely linked to the histological transition rather than genetic changes traced by this approach. Expression of CLDN4 was consistently lower in the invasive compartment than to each DCIS clone. However, for some genes such as ACTB, expression patterns changed in opposing directions in the invasive cancer relative to the sampled DCIS clone (expression is higher than P1-green DCIS (FDR = 0.02) and lower than P1-purple DCIS (FDR = 0.013)) or were highly specific to a genetically more distant DCIS clone (Extended Data Fig. 5f).

Attempts to isolate the changes associated with invasive transition might also be confounded by heterogeneity within the invasive compartment. In P2-TN1, we therefore sought to examine whether the two genetically distinct invasive subclones (P2-blue and P2-purple) were phenotypically distinct. The two cancer clones exhibited distinct morphological (nuclear and architectural) features (P = 0.04, Fisher’s exact test) (H&E image insets; Fig. 3e,f) and occupied neighbourhoods with different stroma (FDR = 0.02) and immune cells such as myeloid cell densities (FDR = 0.08) (mini-image insets; Fig. 3e and Extended Data Fig. 6a–c). Transcriptional programs were also distinct, with statistically significant differences in gene expression for 12 of 91 genes between clones (Extended Data Fig. 6d). Together, these data indicate that the particular clones sampled can have a profound effect on attempts to identify the phenotypic changes implicated in driving or arising during histological progression.

Growth patterns of pre-invasive clones

To demonstrate that BaSISS can be used to chart growth patterns in relation to complex tissue structures, we turned our attention to three DCIS samples from P1 that spanned a tissue surface area of 224 mm2 (P1-D1, P1-D2 and P1-D3) (Fig. 4a and Extended Data Fig. 7a). The adult female breast comprises multiple, branching ductal systems, termed lobes, that extend from the nipple surface to the acini of the lobules, as illustrated in Fig. 4c36,37. DCIS arises from the duct epithelium and is considered a lobar disease as it typically involves the ducts and lobules of a single lobe38. Although DCIS is known to be genetically heterogeneous19, how DCIS clones are organized and grow through the wider duct system remains elusive39.

Fig. 4: Growth patterns and histological associations of DCIS clones.
figure 4

a, BaSISS maps of pure DCIS samples: P1-D1 and P1-D2. The most prevalent genetic clone is projected as a coloured field (which corresponds to b) on DAPI images (reported if the CCF is more than 25% and the inferred local cell density is more than 300 cells per mm2). Scale bar, 5 mm. The quantitative, continuous nature of these data can be examined via an interactive web browser (https://www.cancerclonemaps.org/). The pie chart reports the WGS-estimated clone composition of P1-D1. The white dashed contours delineate morphologically defined lobules. The beige contours mark 114 and 40 manually selected microregions in P1-D1 and P1-D2, respectively, the clonal composition of which is reported by the heatmaps in b. Microregions were manually selected and represent single or small groups of intimately related acini or ductules from the same lobule. b, The phylogenetic tree was inferred from P1 multiregion WGS: branches are scaled according to and annotated with the number of WGS mutations and driver mutation-containing genes. Branches and nodes are coloured to reflect the clones mapped in a. Only branches detected in P1-D1 and P1-D2 are coloured. WT, wild type. c, Cartoon of a lobe of the breast with normal anatomy (left) and DCIS (right), with lobules exhibiting monoclonal and polyclonal involvement. d, H&E images report representative subclone histological features in regions selected from a. Scale bars, 100 µm and 50 µm (vacuoles). e, Stacked bar plot summarizes histological features of microregions dominated by P1-green (n = 66) or P1-orange (n = 72). Nuclear pleomorphism is a measurement of the amount of variability in size and shape of the nuclei and is a major determinant of the histological grade.

Source data

Full size image

The clone maps generated for the three samples formed striking mosaics of mainly green and orange, and occasional blue and grey that localized to areas of histologically confirmed DCIS (Fig. 4a and Extended Data Fig. 7a). Immune clusters and occasional normal or hyperplastic ducts appeared white (unstained), consistent with a different genetic ancestry. In P1-D3, a 3-mm length of a large duct exhibited both a genetic and a histological transition from normal ductal epithelium to DCIS along its length, confirming that, although neoplastic involvement was extensive in this lobe, it was incomplete (Extended Data Fig. 7a). On dividing the glandular tissue into lobules (white dashed contours; Fig. 4a), it was apparent that a handful of lobules contained a single clone, but often multiple clones co-occurred. Indeed, we were surprised to observe that the same clones repeatedly co-existed within lobules that spanned centimetres of tissue. These appearances seem at odds with the traditional model of clonal competition in which a fitter clone generates localized monoclonal sweeps (Fig. 4c).

However, at finer, sublobular resolution, complete or near-complete clonal sweeps are the dominant pattern, as exemplified by assaying 146 representative microscopic regions that represent individual or small clusters of intimately related acini and ducts (beige contours; Fig. 4a). The existence of frequent clonal sweeps as inferred by BaSISS (Fig. 4b) was corroborated by LCM–WGS of additional microregions (Extended Data Fig. 7b). In some instances, including P1-D1-88 (Extended Data Fig. 7c) and P1-D2-0 (Fig. 4a,b,d and Extended Data Fig. 7d–f), clonal interfaces are directly observed within a continuous anatomical space. However, more commonly, rapid clone field transitions (see interactive maps (https://www.cancerclonemaps.org/)) coincided with the myoepithelial cell layer and/or basement membrane that define an acinus or ductule border. It thus transpires that the microanatomical structure of resident tissues can have, an as yet poorly understood, role in shaping observed subclonal architectures (Fig. 4a,c).

DCIS clone-specific phenotypes

Integration of histological and spatial gene expression data from serial sections revealed that the DCIS clones, P1-green and P1-orange, exhibit many phenotypic differences that are consistent across large tissue areas (Fig. 4d,e and Extended Data Figs. 7e,f and 8a,b). Histogenetic associations were very strong, with regions dominated by P1-green being more likely to have an intermediate rather than a low nuclear grade (P < 0.0001; Fisher’s exact test after Bonferroni correction), exhibit more nuclear pleomorphism (P < 0.0001), necrosis (P < 0.0001), vacuoles (P < 0.0001) and a non-solid architectural growth pattern (P < 0.0001) (Fig. 4d,e and Extended Data Fig. 7e,f).

Clone and cell type-resolved spatial gene expression analysis using targeted ISS further corroborated phenotype–genotype correlations. A total of 28 of 91 interrogated genes were differentially expressed by the two main clones (FDR < 0.1, fold change > 1.5 both ways; Extended Data Fig. 8a,b). Consistent with a higher nuclear grade, P1-orange epithelial cells exhibited higher expression of the cell-cycle regulatory oncogenes CCND1 and CCNB1 and the oncogene ZNF703, which have been linked to adverse clinical outcome40. Overall, architectural and nuclear appearances and gene expression profiles were remarkably lineage-specific, and it was particularly notable that these different patterns could also be appreciated spatially, in regions with sublobular, microscopic clone intermixing, adding weight to the clone composition predictions by the model (Extended Data Fig. 7d).

Metastatic clones in a lymph node

Lymph node metastasis is associated with higher rates of cancer mortality41. Whether it has an active role in facilitating cancer progression or simply reflects a more aggressive or distinct biology of certain clones is largely unknown. A substantial challenge is low cancer purity of diffusely infiltrated lymph nodes, which can make it difficult to separate cancer from immune cell-derived molecular signals. To demonstrate that BaSISS can facilitate the simultaneous study of cancer and immune compartments in such challenging cases, we analysed BaSISS, histological annotation and ISS targeted gene expression datasets from sample P2-LN1 (Fig. 5 and Extended Data Fig. 9).

Fig. 5: Intrinsic and extrinsic features of metastatic subclones in a lymph node.
figure 5

a, BaSISS map of P2-LN1, which relates to P2-TN1 (Fig. 3e) and P2-TN2 (Extended Data Fig. 6a,b). The most prevalent genetic clone colours are projected as coloured fields on the DAPI image (reported if the CCF is more than 25%; a threshold of 5% is used in regions of diffusely infiltrating blue to allow visualization in very high normal contamination regions). Scale bar, 2.5 mm. Coloured contours define microregions with distinct metastatic cancer growth patterns (M-GP1 and M-GP2); ‘+’ indicates the surrounding sinus epithelium. b, Plots of the genomic structures in P2-blue and P2-orange clones in the vicinity of the HER2 gene, derived from WGS data of P2-TN2 and P2-LN1. Vertical lines represent genomic rearrangement breakpoints coloured by the phylogenetic tree branch where the event occurred. Dots represent local (binned) copy number. HER2 amplification, CACNB1 fusion and HER2 mutation are BaSISS targets used to track this complex event. BFB, breakage fusion bridge. c, Representative areas of the two main growth patterns stained with H&E. Scale bar, 100 µm. d, Phylogenetic tree inferred from P2 multiregion WGS. Branch and node colours inform the clones mapped in a. The top heatmap reports the BaSISS clone contribution to 39 histologically annotated microregions from a (regions with 5% or more tumour cells are included); see https://www.cancerclonemaps.org/. The bottom heatmap reports microregion histological features. Pan-CK, pan-cytokeratin. e, Volcano plot of immune cell expression of the 62 genes in the ISS immune panel. f, Volcano plot of epithelial cell expression of the 91 genes in the ISS immune panel. Significantly (FDR > 0.1), differentially expressed (fold change of more than 1.5 both ways) genes are coloured. g, Violin plots depict clone-specific cell-type contribution posterior density of the generalized linear mixed model with region-specific random effect, and includes the 22 clone territories with a dominant clone fraction of more than 0.05 in P2-LN1. Significant comparisons were controlled for FDR using the Benjamini–Hochberg procedure.

Source data

Full size image

BaSISS in P2-LN1 targeted 13 trunk and branch alleles, including point mutations and an expressed novel internal fusion in the CACNB1 gene that was co-amplified with the clinically targetable breast cancer oncogene HER2 in a breakage fusion bridge event (Fig. 5b and Supplementary Data Table 1). The model detected two clones (P2-blue and P2-orange) that formed spatially segregated patterns in P2-LN1 (Fig. 5a,d). Only P2-blue was detected in primary breast tumours (P2-TN1 and P2-TN2) (Fig. 3e and Extended Data Fig. 6b).

Detailed histological annotation, blinded to the clone territories, was performed using a combination of H&E, CD45 and pan-cytokeratin IHC and identified multiple metastatic cancer growth patterns (coloured contours; Fig. 5a,c,d and Supplementary Table 2). Intersecting the clone maps and histological annotations revealed strong associations between the two detected clones and the two main histological growth patterns (P < 0.0001, Fisher’s exact test) (Fig. 5d). The P2-orange clone formed monotonous sheets of cancer cells, exhibited weak immunoreactivity for pan-cytokeratin and often occupied sinusoidal structures. By contrast, P2-blue cells stained more strongly for pan-cytokeratin and, when clustered, surround densely packed lymphocyte cores (Fig. 5c,d and Extended Data Fig. 9a–d).

We sought to determine whether transcriptional differences support the spatial inference of clones. Consistent with the known HER2 amplification, P2-orange expressed higher levels of HER2 (Fig. 5f and Extended Data Fig. 9c). A total of 17 of 91 genes were differentially expressed and many of these are implicated in critical biological cancer pathways and/or have recognized prognostic value, including CTSL2, VEGFA (encoding vascular endothelial growth factor receptor A) and CD24 (refs. 42,43) (Fig. 5f). Spatially plotting these genes confirmed that clone-specific expression patterns are recapitulated within multiple, spatially distinct expansions across more than 1 cm2 of tissue (Extended Data Fig. 9a–c).

Integration of spatial transcriptomics data also revealed that metastatic subclones occupied distinct immune microenvironments. Relative to P2-orange cells, P2-blue cells resided in neighbourhoods enriched for T cells and B cells (Fig. 5e,g). In fact, P2-blue cells frequently formed clusters around B cell-rich germinal-like centres, highlighting a potential clone-specific interaction with the adaptive immune system (Fig. 5c and Extended Data Fig. 9a,d). By contrast, P2-orange regions frequently resided inside the lymph node sinuses that were lined by endothelial cells expressing CD34 and PDGFRB (Fig. 5c and Extended Data Fig. 9f). Most of the immune cells in P2-orange regions were myeloid cells with expression profiles consistent with the existence of both M1 and M2 macrophages (CD163, CD68, HAVCR2 and FCGR3A), and the most highly enriched gene, CXCL8, is released by hypoxic macrophages44 (Fig. 5e). Indeed, relative to P2-blue, it emerges that P2-orange experienced more hypoxic conditions manifesting as higher cancer cell expression of VEGFA and necrotic regions (Extended Data Fig. 9e,f). Hypoxia signatures are associated with adverse clinical outcomes, probably because they reflect the emergence of environments that can select for hypoxia-tolerant clones and/or cancer proliferation rates outstrip neoangiogenesis45. Together, these data demonstrate how BaSISS clone maps allow one to spatially relate such variation in microenvironments to individual clones.

Discussion

Here we present BaSISS, a pipeline that combines a highly multiplexed fluorescence microscopy-based protocol and algorithms to map and phenotypically characterize the unique set of subclones of cancer. These maps served as the basis for further spatially and single-cell-resolved molecular and histological characterization of each clone. Applying BaSISS to a series of samples from the key stages of breast cancer progression—carcinoma in situ, invasive cancer and lymph node metastasis—it is notable that virtually every sample exhibited a spatial organization of clones, which warrants further investigation in larger cohorts. The fact that nearly all clones examined in this dataset displayed distinct clone-specific gene expression, stromal and immune microenvironments and microanatomical niches highlights the functional relevance of at least some subclonal diversification.

The ability to chart clonal growth patterns and clone-specific genetic underpinnings of the tumour microenvironment is likely to be instrumental in elucidating how different evolutionary processes operate and manifest across different cancer types—or even in histologically normal tissues46. Understanding the forces of malignant progression, especially invasion and metastasis, and how interactions with the tumour microenvironment shape clinical outcomes10 appear of particular importance. Detailing the functional and microenvironmental characteristics of different clones is also relevant as some part of subclonal diversity in tumours may be due to selectively neutral drift, but the exact extent remains debated.

Particular advantages of the technology are that it is capable of interrogating very large tissue sections on the scale of squared centimetres, which enables studying entire cross-sections of smaller tumours. It is also comparably cheap, unlike solely relying on sequencing-based methods47. The three main limitations of the approach are relatively low sensitivity, which currently precludes single-cell genotyping, a reliance on RNA with the resulting variation in gene expression levels of targeted transcripts, and the fact that clone-defining mutations need to be detected first by separate sequencing-based assays. Greater sensitivity and spatial resolution may be achieved by including more targets per clone and by favouring mutations with higher predicted expression levels, for example, in higher copy number states. A switch to hybridization-based sequencing and direct RNA-binding probes may also improve base-specific detection by several fold48,49. Further discussion of the implications of our observations and limitations of the method is provided in a Supplementary Note.

It is often stated that “nothing in biology makes sense except in the light of evolution”50, which is likely to be true for cancer biology. The ability to spatially locate and molecularly characterize different cancer subclones adds essential features to the spatial-omics toolkit. It provides a robust evolutionary framework that is necessary to interpret the biological relevance of many of the more plastic spatial characteristics of a cancer. Future widespread applications of spatial genomics approaches such as BaSISS will uncover how cancers grow in different tissues and allow us to track, trace and characterize the ill-fated clones that are responsible for adverse clinical outcomes.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.