Clonal evolution and transcriptional plasticity shape metastatic dissemination routes in prostate cancer

Domain: biology
Published: 2025-12-23
Editorial outcome: minor_revision

Abstract
Prostate cancer is a highly heterogeneous disease, driven by genomic and transcriptional changes that impact disease progression and metastatic potential. The interplay between clonal evolution, transcriptional plasticity, and tumour microenvironment is, however, poorly understood. Here, we leverage and integrate single-nuclei RNA sequencing and whole-genome sequencing from 43 spatially distinct tumour samples from five patients with locally advanced prostate cancer to reconstruct clonal evolution trajectories and transcriptional changes driving metastasis at single-cell resolution. We find extensive clonal heterogeneity, including both monophyletic and polyphyletic metastatic dissemination, and ongoing clonal evolution in the primary tumour after metastatic spread. Metastatic seeding converges on disease trajectories involving both genomic and transcriptional changes, including androgen receptor independence and activation of estrogen-, WNT- and JAK-STAT-pathway activity, in spatially distinct areas. Our findings suggest an intricate interplay between clonal evolution and cellular plasticity driving metastatic seeding and point toward more integrative prognostic markers for improved patient management.

Introduction
Prostate cancer is characterised by substantial molecular and clinical heterogeneity1,2. While most patients do not experience fatal outcomes, approximately one in eight patients progresses to aggressive metastatic disease. High intra- and inter-tumoral genetic diversity complicates efforts to unravel the evolutionary trajectories driving metastasis. Effective management of high-risk primary prostate cancer requires a deep understanding of the intricate genetic, phenotypic, and tumour microenvironment (TME) factors that contribute to heterogeneity and metastatic potential, particularly in the context of multimodal therapy concepts, including neoadjuvant and adjuvant treatments as well as targeted systemic therapies.

Our understanding of prostate cancer evolution has primarily come from bulk genome sequencing studies. While most of these studies have focused on single-biopsy tumour material, multi-region bulk sequencing studies have revealed significant intratumoural heterogeneity linked to distinct tumour microenvironments3,4. Different metastatic seeding patterns have been identified from locally advanced prostate cancer, mainly involving monoclonal dissemination (a single primary subclone seeds all metastatic lesions), but multi-region sampling approaches have shown that polyclonal dissemination is not uncommon1,5–7.

While it is believed that high clonal diversity can serve as a strong indicator of future tumour progression2, the metastatic dissemination in prostate cancer is, nevertheless, often driven by a single dominant lineage3,8. The characteristic properties of seeding clones in the primary tumour are poorly understood, but both intrinsic and extrinsic factors play an important role, including tumour microenvironment (TME)9 and signalling pathway activities10. This highlights the need for more comprehensive approaches to unravel the intricate clonal architecture and population evolutionary dynamics of both primary and metastatic lesions.

The cancer cells in prostate cancer are thought to arise from a luminal epithelial cell of origin, situated with different cell types, including luminal, basal, hillock, and club cells, as well as stromal and immune cells11–13. Recent advancements in single-cell and spatial transcriptomic sequencing technologies have advanced our ability to study this cellular and molecular heterogeneity, providing insights into the spatial relationship between diverse genetic and phenotypic profiles of both tumour cells and TME in primary prostate cancer14–17. Despite these advancements, there has been a relative paucity of studies focusing on both the genomic and transcriptional evolutionary aspects of prostate cancer at single-cell resolution. Understanding the evolutionary dynamics of prostate cancer is crucial, as it can reveal the mechanisms driving tumour progression, metastasis, and resistance to therapy.

Here, we integrate single-nuclei transcriptomics with whole-genome DNA sequencing (WGS) from multiregional sampling to explore the evolutionary pathways of prostate cancer and elucidate convergent clonal and phenotypic routes and processes involved in metastatic dissemination.

Results
Pathology-guided whole genome sequencing and single-cell transcriptomics of advanced prostate adenocarcinoma

We examined 43 spatially distinct specimens from five patients with locally advanced high-grade prostate adenocarcinoma with a pathology-guided single-cell transcriptomics workflow to study tumour heterogeneity and metastatic potential (Fig. 1a). Four patients were treatment-naive, while patient PCAL25 received neoadjuvant androgen deprivation hormone therapy six weeks before surgery. Additionally, patient PCAL03 had synchronous bone metastasis (Supplementary Fig. 1A). Following radical prostatectomy and pelvic lymph node dissection, the whole prostate specimens were processed following a whole prostate cryopreservation approach18. Ten regions with high cancer cell content (>80%) were selected and marked in a topographic pathology report (Fig. 1b). Tissue punches were obtained from these areas and, along with corresponding regional lymph node metastasis samples (LNM), subjected to nuclei isolation, followed by single-nuclei RNA sequencing (snRNA) and low-pass bulk WGS (Fig. 1c). The spatial distances from the centre of all prostate samples per patient (centroid) to the analysed area were calculated based on the topographic pathology report.

After quality-based filtering, including the removal of doublets, nuclei with high mitochondrial content, and ambient RNA (“Methods”), we retained more than 356,000 single nuclei from 43 samples (Supplementary Data 1). Each patient contributed between 5 and 9 areas and at least one LNM sample (Fig. 1d, e and Supplementary Fig. 1B). Malignant cells were predicted using Numbat19, and cell type assignment was determined by consensus among several reference-based cell-type prediction approaches, including scPred20 and SELINA21. As expected, the majority of cells captured (64.9%) were cancer cells (Fig. 1d). Two distinct cancer cell populations (Ciliated and Proliferating cells) comprising 0.46% and 1.0% of cancer cells, were defined by motile cilia markers (FOXJ1, DNAI1, DNAH11) and proliferation markers (MKi67), respectively (Fig. 1f). In addition, the proliferating cells exclusively expressed EZH2, a component of the polycomb repressive complex 2 (PRC2), which is implicated in lineage plasticity and cell cycle progression in prostate cancer22. We also detected normal luminal- (LE, 8%), hillock- (HE, 1.29%), club- (CE, 0.6%), and basal epithelial cells (0.1%). The immune compartment consisted of mononuclear phagocytic cells (MNPs, 3.04%), T cells (3.95%), B cells (2.2%) and mast cells (0.08%), while the stromal compartment consisted of fibroblasts (11%) and endothelial cells (3%) (Supplementary Data 2). The PCAL25 tumour exhibited a high degree of fibrosis, as indicated by a substantial stromal compartment contribution of 33.9%, compared to the other patients (median 9.27%, range 4.71–20.4%, Fig. 1e and Supplementary Fig. 1D). In contrast, PCAL37 contained the highest proportion of tumour cells (81%), except for the lymph node sample (55 tumour cells). This patient also exhibited minimal immune compartment contribution compared with the other patients (averaging between 5.9 and 13.5%). We found a high degree of heterogeneity at the tumour microenvironment (TME) level, both between and within patient samples, with TME entropy ranging from 0.3 to 2.7 (indicating low and high TME heterogeneity, respectively) (Fig. 1e). To estimate the spatial distribution of cell types, we computed their Gini coefficient (Fig. 1g), a metric for inequality and cellular heterogeneity. The Gini coefficients for most cell types were normally distributed across patients (Supplementary Fig. 1C), indicating similar variability in their spatial distribution among patients. Malignant and ciliated cells were evenly distributed across areas, as indicated by their low Gini coefficient (Gini coefficient 0=perfect equality, 1=perfect inequality), whereas LE and B cells were among the cell types exhibiting higher Gini coefficients, reflecting greater cellular heterogeneity.

Reconstructing tumour phylogenies at spatially distinct single-cell resolution reveals extensive multi-lineage branching

We examined intra-tumour heterogeneity at the genetic level by integrating our multi-region, high-resolution snRNA-seq and WGS data. To reconstruct the phylogenetic trees for each tumour, we first inferred copy number alteration (CNA) profiles for all areas. We integrated genomic segments from both area-specific snRNA-seq and bulk low-pass WGS data (Supplementary Figs. 2–4), to build phylogenetic trees representing local tumour evolution (Fig. 2a–c, and Supplementary Fig. 5). For comparison, we built phylogenetic trees from single-biopsy bulk-WGS data (Supplementary Fig. 5, 6, and Supplementary Data 3)23. All patient tumours exhibited extensive genetic heterogeneity at the single-cell level, with numerous intra-prostatic lineages (Supplementary Fig. 5). Our single-biopsy bulk-WGS tree structures were overall in agreement with the single-cell multi-region phylogenetic trees, but only captured a minor fraction of the clonal diversity. We note that the branch length of our high-resolution phylogenetic trees is based on the number of CNAs, and is, therefore, only a simplistic estimate of the actual evolutionary timing. In the case of PCAL37, clone PCAL37-C18 represented the most recent common ancestor (MRCA), characterised by copy number losses of chromosomes 8p (NKX3-1), 17p (TP53), 18, and smaller deletions on chromosomes 1q, 2p, and 4q, along with an 8q arm gain including MYC (Supplementary Fig. 5A, B). PCAL37 diverged into two main branches, one marked by a 19q13 amplification, which further evolved into multiple lineages. Additionally, we found support for a subclonal shattering of chromosome 19 in this clade (Supplementary Fig. 4). We identified several lines of evidence for convergent evolution, for example, with clones PCAL37-C11 and PCAL37-C6 both acquiring independent subclonal 12p (CDKN1B) losses. The PCAL37-C1 clone, responsible for metastatic seeding to the lymph node, constituted the majority of the PCAL37 tumour (cancer cell fraction (CCF) of 0.498). Areas T05 and T08, which were spatially close in the prostate (1 unit pairwise distance, with a median 2.6 units, and a standard deviation of 1.2 units for all areas from patient PCAL37), were also similar at the expression and clonal composition level (Fig. 1e, Supplementary Fig. 5A, B and 7A, B). Clones detected in more than one area were from areas in close spatial proximity, as expected. Interestingly, the MRCA clone (PCAL37-C18) was located at the basis of the prostate, far from the largest tumour lesion, suggesting extensive intra-tumour dissemination.

We found numerous subclonal and recurrent copy number gains in PCAL34 (Supplementary Fig. 5C, D). Truncal events included losses of chromosomes 8p and 13, 8q gain, and small deletions on chromosomes 2, 5, and 6. The earliest branching event was characterised by whole-genome doubling (WGD). The other clade remained diploid but further split into multiple lineages, all characterised by large genomic gains of chromosomes 1q21.2-44, 3, and 7. Tumour cells in area T09 contained both tetraploid and diploid clades while a distinct lineage, comprising clones PCAL34-C6 and PCAL34-C5, was restricted to area T06 (Supplementary Fig. 7C). This clone acquired additional copy number gains across the genome, including on chromosomes 3, 6, 12, 14, 17q, 19 and 21 (Supplementary Fig. 4). The observed convergent evolution of independent copy number gains of chromosomes 3 and 7 suggests a selective advantage in this prostate tumour. This is supported by the retention of four copies of these chromosomes after WGD, while most of the genome was reduced to a 3N state (Supplementary Fig. 5D). The PCAL34-C14 clone from area T10 seeded the lymph node and acquired additional gains of chromosomes 5p and 14.

PCAL03 was the largest tumour (Supplementary Fig. 1A). Despite an early WGD event, PCAL03 surprisingly showed the lowest clonal diversity (Supplementary Fig. 5E, F). We found 8p deletion and gains of chromosomes 8q and 20q to be early events, with WGD followed by several whole chromosome losses, present in the majority of the clones, which was supported by our DNA-based WGS copy number profiles across the examined areas (Supplementary Fig. 3A). The largest clone PCAL03-C2 (CCF 0.579) spanning four areas (Supplementary Fig. 7E, F), acquired additional rearrangements on chromosome 11p, and gave rise to smaller subclones PCAL03-C3, PCAL03-C5 and PCAL03-C7, CCF of 0.043, 0.276 and 0.127, respectively (Supplementary Fig. 5E, F). A branching event gave rise to one leaf clone spanning areas T09 and T10 and another leaf giving rise to metastatic seeding, with an additional 4q focal gain and 7p loss following lymph node dissemination. We note that two different LNMs were subjected to single-biopsy deep WGS and to snRNA-seq from PCAL03. Interestingly, while the phylogenetic trees were comparable in the primary tumour, the single-biopsy WGS-based phylogenetic tree inferred a different seeding clone (Supplementary Fig. 7C), suggesting polyphyletic metastatic dissemination, where several distinct clones in the primary tumour independently seed metastases.

Acquired mismatch repair deficiency and JAK-STAT signalling in metastasising lineage

Genomic instability is associated with aggressive disease and metastatic competency24,25. We found a particularly striking example of copy number instability in PCAL10, with early truncal deletions of chromosomes 6q13-23.3, 8p, and 18, followed by WGD in one branch, which was not detected in the single-biopsy bulk WGS tree, while the other branch remained diploid (Fig. 2a, c). In the diploid lineage, clone PCAL10-C2, with the least percentage of genome altered (PGA), emerged as the largest clone (CCF = 0.459), spanning three areas with the highest proportion of proliferating tumour cells (Fig. 1, Supplementary Fig. 8). We also identified monophyletic seeding from PCAL10-C2 to the lymph node. Expression-based clustering identified T04 as the most likely metastatic seeding area (Fig. 2b–d). Interestingly, we found the lowest Gleason grade in this area (Fig. 1e), suggesting that pathology-based tumour grading should be accompanied by genomic-based analysis to identify clones with metastatic potential. A second branching event in the diploid lineage gave rise to clones with major genomic losses/gains, all confined to areas T09 and T10. Moreover, both areas contained clones from diploid and tetraploid lineages (Fig. 2d), pointing to local migration and intermixing of distinct clonal lineages within an individual prostate tumour. To further investigate somatic alterations in these regions, we performed deeper WGS on areas T02 and T09. Area T02, which contained the PCAL10-C2 clone, exhibited a biallelic deletion at the MSH2/6 locus (2p21-p16.3, expression log2 fold-change=-3.7, T02 versus T09) and concomitant increased levels of the DNA mismatch repair (MMR) mutational signatures SBS15 and SBS21, and a microsatellite instability (MSI) phenotype compared with area T09 (Supplementary Fig. 9).

WGD is a prevalent event observed across various cancer types and is strongly associated with genomic instability, metastatic disease and poor clinical outcomes26,27. Relying solely on bulk WGS makes it challenging to accurately determine the subclonal composition of a sample. For example, our multi-region WGS alone indicated area T05 to be 4N, with subclonal single-copy losses of chromosomes 1, 3 and 20 (Supplementary Fig. 10), but our snRNA-seq-based clonality revealed two distinct cell clusters in almost equal proportions, with divergent loss of chromosomes 1 and 3 in clone PCAL10-C5 and chromosome 20 in clone PCAL10-C4.

Metastatic seeding is a crucial step in the progression of localised prostate cancer to lethal metastatic disease. Therefore, we aimed to identify pathway activity unique or shared across subclones within the entire tumour, as well as specific to the seeding clone, using decoupleR28. As expected, pathway activity was similar between clones derived from the same lineage within an area (Fig. 2e). For instance, clones PCAL10-C4 and PCAL10-C5 clustered together based on pathway activity, in contrast to PCAL10-C2 from areas T02, T04, and T06. We found enriched JAK-STAT activity, associated with lineage plasticity and resistance to Androgen Deprivation Therapy (ADT)29, exclusive to the PCAL10-C2 clone, regardless of the area (Fig. 2e, f). Notably, the seeding clone and metastatic clones (in the LNM) exhibited low androgen activity. Based on mean androgen activity, we grouped clones into androgen-positive and androgen-negative, which was further confirmed by inferred AR transcription factor (TF) activity (Fig. 2g and Supplementary Fig. 11A). AR-positive clones had higher NCOA3 activity, a steroid receptor coactivator30, while AR-negative clones were enriched in SHOX2 and MECP2, both of which promote cell growth and proliferation31,32. These findings suggest that AR-negative copy number unstable lineages can drive the clonal evolution towards metastatic seeding.

Early chromothripsis followed by polyphyletic metastatic seeding to regional lymph nodes

Chromothripsis is a one-off genomic shattering process, often encompassing simultaneous disruption of several cancer genes and associated with aggressive disease33. Although the clonality analysis included a limited number of samples (see Methods), our clonal reconstruction of PCAL25 primary tumour identified clone PCAL25-C3 in area T10 as the MRCA, with TP53 loss and chromothripsis of chromosome 6 (Fig. 3a and Supplementary Fig. 3C). We found the MRCA to branch off, leading to parallel multi-lineage dissemination to two regional lymph nodes, LNM1 and LNM2. Single-cell analysis revealed that metastatic tumour clades originated from polyphyletic clonal lineages within the primary tumour (Fig. 3b). In contrast, only a single monophyletic metastatic spread was detectable in the single biopsy bulk WGS tree (Fig. 3a). We observed two clades seeding LNM1, a direct descendant of PCAL25-C3 in T10, and another with acquired WGD in areas T02 and T05. This was followed by multiple subclonal amplifications in chromosomes 7, 10, 11, 13, and 20, as well as large deletions in chromosomes 2, 4, 9, 15, 16q, and 22q (Supplementary Fig. 5H). We traced a third metastatic seeding to LNM2 back to area T10 (Fig. 3a, b).

Next, we examined common pathway activities for clones of different lineages within each lymph node and seeding clones in the primary tumour. Again, we found low androgen activity in seeding clones (Fig. 3c, d). The PCAL25-C3 metastatic lineage displayed transcriptional plasticity, with exclusive WNT activity in the metastatic LNM1, while AR and JAK-STAT pathway activity was increased in the other metastatic area, LNM2. In support, we found differential enrichment of TFs in each LNM (Fig. 3e, Supplementary Fig. 11B), with LNM1 and the seeding area having higher activity of SATB2 and TCF7, a TF negatively regulated by AR34. Both TFs are implicated in activating the WNT pathway and epithelial-mesenchymal transition (EMT)34,35. Samples with low tumour purity are often excluded from WGS analysis, as low CCF can hinder the detection of somatic mutations and lead to inaccurate reconstruction of clonal structures36. Using our snRNA-seq approach, we estimated LNM1 CCF to be 0.172. Despite low tumour purity, we were able to detect multi-lineage seeding events. In total, we identified 14% of the tumour cells to be diploid (114 cells) and 86% of tetraploid lineages (689 cells) (Figs. 3f, g, 2N cells (PCAL25-C2; C3; C4) in orange. 4N cells (PCAL25-C1) in blue), providing further support for WGD following the initiating chromothripsis event on chromosome 6. The single-cell resolution also allowed us to explore further clonal heterogeneity (Fig. 3h and Supplementary Fig. 12). We detected heterogeneous expression of long non-coding RNAs, implicated in prostate cancer tumorigenesis, such as ARLNC1 and SChLAP137. All 4N clones exhibited higher expression of ONECUT2, a TF implicated in suppression of AR and progression to lethal disease38. Notably, PCAL25-C2 from the 2N lineage was the most distinct from all LNM1 resident clones, while PCAL25-C3/C4 clones clustered together with a subset of C1 cells.

These findings suggest convergent transcriptional programmes across different lineages, including AR suppression, JAK-STAT and WNT activity, even with distinct genotypes.

The impact of spatial organisation on clonal evolution and transcriptomic profiles in primary prostate cancer

The clonal composition of adjacent areas was more similar than distant ones, as expected (Figs. 2, 3). We also found a positive correlation between the phylogenetic distance between tumour cells and the TME composition (P=0.01, R=0.25, Pearson correlation, Supplementary Fig. 13A). To complement our centroid-based distance analysis (Fig. 1), we also computed pairwise Euclidean distance between each tumour area, to investigate the spatial relationships between individual clones. Neither the phylogenetic distance nor TME were proportional to the pairwise distance between areas (Supplementary Fig. 13B, C), suggesting a connection between TME and tumour cell composition at shorter distances, but also that tumour and normal cell types follow non-linear distributions at larger distances within a prostate tumour. We observed tumour cells preferentially in the peripheral and central zones, with higher clonal diversity (entropy) in the central zone (Supplementary Fig. 13E, D). In contrast, the transition zone was occupied by fewer tumour cells with lower clonal diversity (Supplementary Fig. 13F), suggesting that clonal expansions may be more constrained in this region.

Our spatially resolved analyses implied a strong interplay between transcriptional programmes and the clonal evolution in prostate cancer. AR signalling is the most common driver for prostate cancer growth, but accumulating evidence suggests AR-independent transcriptional reprogramming can drive lineage plasticity in tumour cells29,39. Overall, pathway activities exhibited varying degrees of intra- and inter-tumour heterogeneity, with, for example, Oestrogen, androgen and JAK-STAT pathway activity having high and low intra-tumour heterogeneity, respectively (Gini index, Supplementary Fig. 14F), and with the majority of pathways exhibiting patient-specific activities (pairwise variance test, Supplementary Data 4 and Supplementary Fig. 15). For example, PCAL25 (ADT-treated) was among the tumours with the lowest AR activity, as expected. To investigate the contribution of transcriptional programmes to the clonal evolution and metastatic seeding areas across all patients samples, we integrated and compared TME components, tumour cell states, which represent gene modules derived from a pan-cancer single-cell RNA-sequencing study40 including prostate cancer, and key pathways within and between patient samples using principal component analysis (PCA) (Fig. 4a). We found PC1 (20.79%) to be explained primarily by (positively correlated with) tumour content and negatively correlated with MNPs, TME diversity, Grade group and T cell proportion. PC2 (12.66%) was explained primarily by lymphocyte content and normal luminal prostate cells. The presence of areas with low tumour content, high immune cell content and high Gleason grade (Fig. 4a and Supplementary Fig. 16) suggests immunologically ‘warm’ areas in the prostate with tumour-immune interactions driving more aggressive cell phenotypes. High immune cell content also correlated with MSH2 homozygous deletion, whereas TP53 loss (17p) and MYC gain (8q) were both linked with immune cell content depletion and low JAK-STAT signalling (Supplementary Fig. 16B), supporting the previously established immune suppressive roles of TP53 and MYC mutations in cancer41,42. We also found a strong contribution of proliferating tumour cell content and JAK-STAT activity in PC2 (Fig. 4a). JAK-STAT activity, which we also found elevated in the seeding area of PCAL10 (Supplementary Fig. 14C–F), can drive immune evasion and increased lineage plasticity in prostate cancer in an AR-independent manner29. Changes in EMT signalling and mast cells correlated with PC3, whereas PC4 was mainly explained by ciliated tumour cells, endothelial and fibroblasts, as well as metal response, which is linked with homoeostasis, cell differentiation, and proliferation40. PC5 was explained by stress signalling and tumour areas toward the periphery of the tumour (high centroid distance).

Interestingly, we found metastatic seeding areas correlating with PC6 and PC7 together with ciliated and proliferating tumour cells. To further characterise these cells, we investigated their spatial properties and key signalling pathways. We found ciliated cells to be located preferentially near the centre of the prostate gland (anti-correlation with centroid distance), whereas the metastatic seeding areas were located preferentially towards the periphery (Fig. 4b, Supplementary Fig. 16). Ciliated tumour cells were also enriched for TRAIL (Tumour necrosis factor-related apoptosis-inducing ligand) pathway and depleted for WNT signalling. In agreement, WNT signalling is known to repress TRAIL-induced apoptosis and immune response modulation43. Hormone signalling plays a pivotal role in prostate cancer dissemination. Whereas androgen response genes were reduced at the transcriptional level in the metastatic seeding area, we found significant enrichment of oestrogen pathway response genes (P=0.017, R=0.39, Pearson correlation, Fig. 4b), which are important for progression in primary untreated prostate cancer44,45. We also found seeding areas to exhibit a trend towards increased WNT and TRAIL signalling, both markers of cell plasticity.

These findings imply spatially distinct areas with different clonal and TME compositions, and suggest that seeding clones, phylogenetically closely related to the metastasis, are localised in spatially confined regions in the primary tumour. These regions tend to localise towards the periphery of the prostate, populated by tumour cells with elevated levels of oestrogen and WNT signalling and reduced androgen pathway signalling (Fig. 4c), accompanied by extensive and recurrent rewiring of transcriptional programmes in primary tumour cells.

Methods
Patient material

Sample collection, consent and clinical data. Five prostate cancer patients with high-risk, locally advanced primary prostate adenocarcinomas were selected from the German ICGC cohort23. Written informed consent and ethical approval from the patients were obtained in accordance with ICGC guidelines. The study was approved by the Ethics Committee of the Ärztekammer Hamburg (Hamburg Medical Association) under protocols PV3552 and PV4679. The study followed the Declaration of Helsinki guidelines and did not use statistical methods to determine sample size due to its exploratory nature. All patients underwent radical prostatectomy with lymph node dissection as the primary treatment for prostate cancer, except for patient PCAL25, who received neoadjuvant androgen deprivation therapy (Gonadotropin-Releasing-Hormone (GnRH) analogue) 6 six weeks before surgery. Patient PCAL03 presented with oligometastatic disease with bone metastasis at the time of surgery. Two different synchronous LNMs were obtained from this patient, one for the initial single biopsy WGS and another LNM for multi-region snRNA-seq.

After surgery, following a predefined scheme as described in Gerhauser et al.23, specialised pathologists dissected each prostate. A single biopsy was isolated and subjected to deep bulk WGS23. The precise location of this sample in the 3D coordinate system of the prostate is unknown. For multi-region sample-isolation, tissue blocks were embedded in a tissue-freezing medium, cryopreserved and stored at −80°C. 10 punches (1.5 mm in diameter) were taken from different areas of the primary tumour and matching lymph node metastases from each patient. In total, we obtained 66 samples: 50 primary tumour biopsies, 16 lymph node biopsies. The selected tumour areas had more than 80% of tumour cell content, according to the pathologist. Gleason grading was reported for each biopsy.

Single-nuclei transcriptomics

Tissue dissociation. Fresh-frozen prostate cancer tissue punches were processed using an adapted version of Pfisterer et al. protocol53. Shortly, frozen tissue punctures were mechanically dissociated using dounce-homogenizer: 5–10 strokes with pestle A, following 20–25 strokes with pestle B, in an ice-cold homogenization buffer (HB): NIM with 1 mM DTT (Invitrogen), 1X Protease inhibitors (Roche), 0.4 U/µl RNAse inhibitors (Takara), 0.2 U/µl Superasin (Invitrogen) and 0.1% (v/v) Triton X-100. To make HB we first made a nuclei isolation medium (NIM): 250 mM sucrose, 25 mM KCl, 5 mM MgCl, and 10 mM tris buffer, pH 8. Homogenised tissue suspensions were filtered through a 40 µm mesh filter and centrifuged at 1000 × g for 8 min, 4°C. To remove debris, samples were subjected to equilibrium density centrifugation: sample pellets were first resuspended in 250 µl HB and then mixed 1:1 with 50% iodixanol solution (50% v/v OptiPrepTM (Sigma), 25 mM KCl, 5 mM MgCl, 10 mM tris buffer, pH 8.0, 1 mM DTT, 1X Protease inhibitors, 0.4 U/µl RNAse inhibitors, 0.2 U/µl Superasin). 500 µl of the resulting suspension was layered over 500 µl 29% iodixanol solution (29% v/v of 50% iodixanol, diluted with NIM, supplemented with DTT (1 mM), Protease inhibitors (1X), RNAse inhibitors (0.4 U/µl), Superasin (0.2 U/µl)) in a 2.2 ml Ultra-Clear Centrifugation Tube (Beckman Coulter). Samples were centrifuged at 14,000 × g for 22 min, 4°C. Then, the remaining nuclei pellet was resuspended in 1X PBS, supplemented with 2.5 mM MgCl, 1 mM DTT, 0.50% BSA and 0.2 U/µl RNAse inhibitors. Sample quality was assessed under a microscope and nuclei were counted with a hemocytometer.

snRNA-seq library preparation and sequencing. Single-nuclei suspensions containing 16,000 nuclei were loaded to the Chromium Next GEM chip (10X genomics), aiming to capture 10,000 nuclei per area. snRNA-seq libraries were constructed either with Chromium Next GEM Single Cell 3′ v3 or v. 3.1 Reagent Kit, following manufacturer’s protocol. 10X single-nuclei libraries were sequenced on Illumina NovaSeq-6000 platform, using NovaSeq S4 flow cell with paired-end 150 bp mode (28×10x10×281 bp configuration), at the Department of Genomic Medicine at Rigshospitalet, Denmark. We sequenced on average 46,000 reads (range 22.5K-92K) per cell (Supplementary Data 1).

snRNA-seq processing. We used 10X Cell Ranger software (version 5.0.0) with default parameters for alignment to human genome reference hg38 (including intronic sequences), demultiplexing, barcode processing and gene x cell count quantification. To remove technical artefacts from ambient RNA contamination, CellBender54 (version 0.2.0) was used with the following parameters: 10,000 expected cells, 25,000 total droplets included, a low count threshold of 15, and all remaining parameters set to default. In samples with higher background RNA contamination, the low count threshold parameter was increased accordingly.

The resulting corrected and raw gene x cell count matrices were analysed with scCustomize package (R 4.2.0) to create Seurat objects (version 4.3.0). Cells with fewer than 500 genes and genes present in fewer than five cells were excluded from further analysis. The SCP package (version 0.5.1, available from: https://github.com/zhanghao-njmu/SCP) was used to compute quality control metrics and perform initial filtering. Doublets were predicted using scDblFinder (available from: https://github.com/plger/scDblFinder) with default parameters, and cells with mitochondrial transcript content exceeding 1% were excluded.

An additional filtering step was implemented to remove low-quality or dying cells, based on the assumption that mitochondrial reads should be minimal in single-nuclei data. Apoptosis scores for each cell were calculated as the ratio of mitochondrial percentage to total RNA counts. The mean apoptosis score was then calculated for each cluster, and clusters with a mean apoptosis score exceeding one standard deviation above the mean were excluded from further analysis.

The filtered and corrected snRNA-seq data were normalised using the “SCTransfrom” function from Seurat (version 4.3.0) with vst.flavour= “v2”, while regressing out mitochondrial percentages, keeping all other parameters at their default settings. We used the “FindVariableFeatures” function to identify 3000 variable genes, the “RunPCA” function for dimensionality reduction with npcs=20, and the “FindNeighbors” function to identify nearest neighbours, exploring a range of resolutions between 1 and 2. Dimensionality reduction was further performed using the “RunUMAP” function based on the first 20 principal components. All the above functions are from the Seurat package (version 4.3.0), with default settings for the remaining parameters.

Cell annotation and Integration. scPred, version 1.9.220 was used for single-nuclei preliminary data annotation with publicly available prostate cancer datasets49,55. SELINA.py tool, version 0.121, with “--disease” mode and default parameters, was used for cell-type prediction. We labelled cell types using each dataset as a reference, which included a literature-based reference set13,16,56. In addition, “normal prostate” and “NSCLC” references, provided by SELINA, were used for prediction. For cell cluster identity assignment, we took the consensus of all predictions, except for the assignment of proliferating and ciliated cells. We calculated the G2M score and a custom literature-based ciliated gene-set using JASMINE (v1, available from: https://github.com/NNoureen/JASMINE). Cells that had a high ciliated or proliferating score (>threshold=median+one standard deviation) were assigned as ciliated or proliferating, respectively. To integrate all patient data, we utilised STACAS, version 2.2.1, a semi-supervised method developed by Andreatta et al.57. We employed cell type annotation as input, using the default parameters. Cancer cell fraction was subset, normalised with “SCTransform” normalisation (Seurat v5) on raw counts, and clustered for visualisation.

Pathway, transcription factor activity and cell states scores analysis. We computed pathway activity inference on single-cell count data with decoupleR-py28, version 1.8.0, together with a curated collection of pathways from PROGENy58 and transcription factor (TF) database from CollecTRI59. For each pathway set, we used the top 500 genes ranked by p-value and Multivariate Linear Model (mlm) method for pathway activity analysis. TF enrichment scores were inferred by running the Univariate Linear Model (ulm) method. The resulting activity scores represent the predicted regulatory activity of each pathway or TF in individual cells, based on the expression of their downstream targets. These scores are relative within each dataset, making them comparable between patient samples. Scores reflect the direction and magnitude of pathway activity compared to other cells in the same patient’s sample, rather than providing an absolute measure across patients. Mean pathway/TF activity per cell group was calculated and visualised with pheatmap using R version 4.3.3., python version 3.9.16.

Cell state scores were calculated on normalised single-cell expression data using the AddModuleScore function from Seurat (version 4.3.0.1). Gene sets defining cell states were obtained from Barkley et al.40. We calculated the proportions of cells exhibiting cell state activity per sample, using two different cutoff strategies based on the overall distribution of activity values. For each area and each cell state, if the mean percentage of active cells for a given signature was <70%, a binary cutoff of 0 was applied. Otherwise, the signature’s mean activity was used as the threshold.

Gini index, Shannon entropy, weighted phylogenetic dissimilarity, cell population dissimilarity and physical distance. We assigned each biopsy sample a 3D coordinate on a Cartesian map and calculated pairwise Euclidean distances. The pathology-provided tumour diameter was used to estimate the unit conversion, with 1 unit corresponding to 1 cm. The resulting pairwise distances were as follows: PCAL37 (median: 2.6 units

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting study that effectively integrates single cell transcriptomics with whole genome sequencing to produce insights into prostate cancer spread. The study appears to be methodologically sound, and the results are very nicely presented. The finding that two distinct cancer cell populations exist, ciliated and proliferating, potentially in distinct zones of the prostate is interesting. The limitation, as the authors acknowledge, is the small number of patients.\nMajor comments:\nMetastatic seeding pockets – This is a big (and interesting) claim, but with insufficient evidence. It must be noted that in the absence of data from the entire prostate (which is impractical), one cannot comment on whether a metastatic clone is arising from a specific spatial location based on phylogeny. One can only claim that the specific clones in the primary tumour site(s) are phylogenetically linked to the lymph node mets. In actuality, the metastasis may have occurred from a different site in the prostate with a similar clonal composition.\nSpatial information – Encoding location as distance from the centroid is novel, but ultimately too simplistic and not tracking biology in my opinion. It helps describe the spread of the tumour in the different prostate zones, but the authors themselves state that “…tumour and normal cell types follow non-linear distributions at larger distances within a prostate tumour.” Tumours can follow convoluted routes within the prostate in 3D and Euclidean distance from the centroid may not be the best way to encode this information. For example, regions with a similar distance to the centroid could also be far away from each other if they are situated on opposite ends (anterior vs. posterior) of the prostate. Perhaps the authors could also incorporate direction, treating the distance as a vector instead of a scalar?\nThe study is descriptive, owing to the small number of patients studied. This in itself shouldn’t be a problem, but broader conclusions can’t be drawn from this small set of samples. For example, while the association with JAK-STAT pathway in lymph node metastasis is interesting, it can’t be generalised.\nThe authors state that “To date, our understanding of prostate cancer evolution has primarily been derived from bulk genome sequencing studies. These studies typically rely on single biopsies from the primary tumour, which may not fully capture the extent of genetic heterogeneity present within the entire tumour.” – but I believe this is not the case, there have been several studies (dating back to Cooper et al., Nature Genetics 2015) that used multi-region sequencing to study prostate cancer evolution.\n- Hong et al https://www.nature.com/articles/ncomms7605\n- Warner et al https://www.nature.com/articles/s43018-023-00692-y\n- Singhal et al https://www.nature.com/articles/s41467-024-48629-y\n- Rao et al https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-024-01302-x\n- Nurminen et al https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-023-01242-y\n- Boutros et al https://www.nature.com/articles/ng.3315", "reviewer_2": "The manuscript “Clonal Evolution and Transcriptional Plasticity Shape Metastatic Dissemination Routes in Prostate Cancer” by Mikutenaite et al. describes the clonal evolution of 5 prostate cancers in relation to the transcriptional changes. The study includes multiple well-portrayed examples about why single-sample WGS biopsies are not enough to depict the total breadth of clonal evolution in prostate cancer.\nMajor issues:\nPresented work needs to be tied more clearly with current literature and put into context of state of the art. The manuscript refers to a limited number of studies in the introduction and discussion. For example, several studies report metastatic seeding patterns from locally advanced prostate cancers based on spatially distinct genome sequencing. These works should be acknowledged more clearly and putative novelty of the presented work articulated based on complete literature.\nTranscriptional plasticity or the results of other prostate cancer single cell transcriptomics (or spatial transcriptomics) papers are not addressed in introduction. Mentioning generally “signalling pathway activities” does not tell much. On the other hand, e.g. in discussion they mention “a recent study identified four major groups of androgen-independent transcriptional programmes (28)”, but do not elaborate in relation to their results.\nAs there is WGS data on every patient (and sample?), could you report the driver mutations in these cases, and how they are related to the inferred clones? Why are they not involved in tumor evolution reconstruction?\nWhile the manuscript is well-written, the text is occasionally a bit difficult to follow, especially keeping in mind which patient the results address, or if a clone in question has a special role, e.g. in seeding. Would it be possible to call the clones with patient name like PCAL03-C1 throughout the text?\nMinor issues:\n- It should be stated in the abstract that these are locally advanced prostate cancers.\n- The Figure 1B is not clear. What are the smaller areas on the left and right side?\n- Is the PCAL37 LNM1 sample a metastasis? it seems to contain only normal cells based on Figure 1E.\n- Authors should check the references to supplementary figures, they might be off in some cases, e.g. page 6.\n- Not all data is shown, e.g. “subclonal shattering of chromosome 19 in this clade”\n- The WGDs are not marked in all trees in Figure S4.\n- Why are not all tumor areas included in the evolutionary reconstruction of PCAL25?\n- In case PCAL03, the authors state that “The metastatic seeding occurred from a branching event derived from areas T09 and T10”. However, this does not seem to be the case based on the evolutionary tree in Figure S4E.\n- Is MSH2 expressed in the T02 area with biallelic deletion (case PCAL10?)?\n- On page 8, last sentence, how does the copy number instability drive the evolution towards AR-negative lineage?\n- Page 15, “tumour cell states” are mentioned for the first time, what are these states?\n- Is it really shown that the metastatic seeding occurred before the WGD in PCAL25? They are in different branches and there’s no timing analysis presented.\n- The project clearly has generated a lot of results and it is understandable that not all are brought up in the manuscript text or showed in the supplementary figures. However, some observations that could be elaborated. For example,\n• area T09 in PCAL34 has both diploid and tetraploid cells?\n• single biopsy tree has wrongly inferred a different seeding clone than single cell data, in case PCAL03.\n• low pass WGS results are not shown (Figure S2 and S3 legends say they are from deep sequencing)", "reviewer_3": "The manuscript by Mikutenaite et al. investigated clonal relationships in five locally advanced, high-risk prostate cancer (PCa) patients using single-nucleus RNA sequencing (snRNA-seq) integrated with low-pass whole genome sequencing (WGS). All patients underwent radical prostatectomy with lymph-node (LN) dissection; one patient (PCAL25) had received neoadjuvant androgen deprivation therapy (ADT), while the others were treatment-naïve. The manuscript effectively comprises five case studies in which tumor samples from multiple spatial regions including LN metastases subjected to phylogenetic reconstruction. The use of snRNA-seq layered on genomic phylogenies provides higher-resolution inference than standard bulk WGS and allows for mapping of transcriptional states onto clonal architectures.\nWhile the genomic and phylogenetic findings themselves (e.g., examples of both monophyletic and polyphyletic metastatic seeding) are not novel in the context of prostate cancer, the integration of transcriptomic data with clonal evolution is an interesting and valuable aspect. For example, the observation that metastases tend to arise from clones with reduced androgen receptor (AR) activity and enhanced expression of pathways associated with lineage plasticity is potentially clinically relevant. However, the manuscript currently lacks sufficient depth in both the analysis and interpretation of these biological and clinical implications.\nWith limited number of patients including only one treated patient, the study is limited in the scope. Further development of the analysis and more insightful contextualization of the transcriptomic findings and addition of treated patients will strengthen the finding.\nMajor comments\n1. The authors should discuss how the findings for PCAL25 relate to the context of prior exposure to ADT. For example, the increased evidence of fibrosis in this patient is probably a consequence of inflammation and regression of the prostate in response to ADT. It is also surprising to see several clones with apparently high AR activity in a newly diagnosed PCa patient after 6 weeks of ADT – is this because pathway/activity scores are only relative to other clones in the same patient, meaning that some clones must always be considered “high” and others “low” in every patient? If samples from different patients were run in multiplex for snRNA-seq (please clarify this) to minimize batch effects between patients, pathway/activity scores should be made comparable across patients if not already. Are the high Androgen/AR scores in PCAL25 similar to the high scores in other patients? Also, do the authors have any pre-ADT biopsies from PCAL25 that they could add to their analysis?\n2. It is difficult to picture how the samples in the various transcriptomic heatmaps line up with the phylogenies because the sample orders are different, and closely related clones do not always have similar colours (e.g. in Fig 2c, there are blue clones in three different branches). Although heatmaps clustered by transcriptomic profiles/scores are informative and should be kept (possibly as Supplementary Figures), it would be more helpful to have the samples and dendrograms in transcriptome-based heatmaps to be arranged and oriented as in the genomic phylogenies. This would better illustrate whether closely related clones have similar transcriptomes and whether more distantly related clones display convergent evolution.\n3. Most of the biologically interpretable transcriptomic data is presented at the resolution of clones. However, there is clearly evidence for transcriptomic heterogeneity within individual clones (e.g. Fig 3F,H). The authors should explore intra-clone transcriptomic heterogeneity in more depth. For instance, they could investigate how much transcriptomic heterogeneity exists within individual clones compared to variability between clones.\n4. The authors should discuss whether transcriptional profiles reflect key genomic aberrations. For example, do clones with gain of MYC or FOXA1 have higher activity scores for those transcription factors? Do losses of specific tumour suppressors, e.g. TP53, associate with activation of plasticity pathways? The authors say in the Discussion that they “observed tumour areas containing clones with homozygous MSH2 loss to have higher lymphocyte infiltration”, but this isn’t presented anywhere in the Results despite being an important, albeit expected, observation. And so on.\n5. The authors need to provide more details about the single-biopsy bulk WGS. In each patient, which tumour site was that biopsy taken from? What were the methods for the analysis? Etc.\n6. The phylogenetic trees from single-biopsy WGS and multiregional snRNA-seq for PCAL03 (Fig S4E) show completely different clones being inferred as the clone that seeded the LN metastasis, while truncal whole-genome doubling was inferred in the single-biopsy WGS-based tree but nowhere in the multiregional tree. How do the authors explain these discrepancies, and do they affect biological interpretation of results?\n7. The authors should discuss how their work fits within the broader field of prostate cancer genomics."}, "author_response": "Nature Communications NCOMMS-25-31194.\nWe would like to thank the reviewers for their appreciation of our work and the highly constructive\ncomments and suggestions, which we have responded to point-by-point below. Briefly, we have\ni) improved clarity of the interpretation of metastatic seeding, rephrasing our claims to emphasize\nphylogenetic links rather than spatial origin, ii) enhanced our spatial analysis, adding further\nsupport to our use of centroid-based and pairwise Euclidean metrics, iii) further contextualised\nour findings, citing the important literature that led up to our work and the context with which our\nresults should be interpreted, iv) performed new analysis supporting genomic-transcriptomic\nintegration analysis, v) improved clarify and labelling of the clonal trees, vi) performed a new intra-\nclonal and transcriptomic heterogeneity analysis, adding further support to the intra-tumour\nheterogeneity at both genomic and transcriptional level, vii) clarified sampling and methodology\napproaches, in particular our usage of whole-genome sequencing and sampling strategy.\nWe believe that our changes have substantially improved the manuscript, both in terms of\ninterpretability and scientific significance.\nREVIEWER COMMENTS\nReviewer #1\nexpertise in prostate cancer genomics and clonal evolution (Remarks to the Author):\nThis is an interesting study that effectively integrates single cell transcriptomics with whole\ngenome sequencing to produce insights into prostate cancer spread. The study appears to be\nmethodologically sound, and the results are very nicely presented. The finding that two distinct\ncancer cell populations exist, ciliated and proliferating, potentially in distinct zones of the prostate\nis interesting. The limitation, as the authors acknowledge, is the small number of patients.\nMajor comments:\n1. Metastatic seeding pockets – This is a big (and interesting) claim, but with insufficient evidence.\nIt must be noted that in the absence of data from the entire prostate (which is impractical), one\ncannot comment on whether a metastatic clone is arising from a specific spatial location based\non phylogeny. One can only claim that the specific clones in the primary tumour site(s) are\nphylogenetically linked to the lymph node mets. In actuality, the metastasis may have occurred\nfrom a different site in the prostate with a similar clonal composition.\nResponse: We fully agree that without exhaustive spatial sampling of the entire prostate, it is not\npossible to pinpoint the exact spatial origin of metastatic clones. Our phrasing in the manuscript\naimed to convey that the metastatic clones we observed are indeed phylogenetically closely\nlinked to clones present in the sampled primary tumour sites.\nTo clarify this point and avoid overinterpretation, we have revised the relevant sections of the\nmanuscript to emphasize that our findings demonstrate phylogenetic relationships between\nmetastases and specific primary tumour regions sampled, and that our proposed framework for\nunderstanding metastatic seeding is of a hypothetical nature. Additionally, we have added a brief\ndiscussion of this limitation in the manuscript to ensure transparency and guide future studies that\nmight explore more comprehensive spatial sampling.\n2. Spatial information – Encoding location as distance from the centroid is novel, but ultimately\ntoo simplistic and not tracking biology in my opinion. It helps describe the spread of the tumour in\nthe different prostate zones, but the authors themselves state that “…tumour and normal cell\ntypes follow non-linear distributions at larger distances within a prostate tumour.” Tumours can\nfollow convoluted routes within the prostate in 3D and Euclidean distance from the centroid may\nnot be the best way to encode this information. For example, regions with a similar distance to\nthe centroid could also be far away from each other if they are situated on opposite ends (anterior\nvs. posterior) of the prostate. Perhaps the authors could also incorporate direction, treating the\ndistance as a vector instead of a scalar?\nResponse: We appreciate the reviewer’s insightful comment and suggestions regarding our use\nof distance from the centroid to encode spatial location within the prostate tumour. We agree that\nthis metric, while useful as a first approximation, does not fully capture the, sometimes, convoluted\nspatial relationships present within prostate tumours.\nTo clarify, we used two distance metrics with distinct properties: a centroid-based metric (Figure\n1), which enabled us to compare the relative spatial distance between each 3D patient tumour,\nand a pairwise Euclidian distance analysis (Figure S12), comparing the spatial distance between\nall pairs of tumour foci. We realise that these two metrics were not sufficiently well explained, and\nwe have now expanded the description and clarified the analysis conducted in the revised\nmanuscript.\nWe have also expanded the discussion to emphasize these points and the limitations of the\ncurrent approach, acknowledging that the complexity of tumour spread within the prostate may\nrequire more detailed 3D spatial data in future studies. Our analysis of the directional relationship\nbetween clones and anatomical zones allowed us to investigate the spatial heterogeneity beyond\nsimple radial distances. This warrants further investigation, which we now emphasise in the\nrevised discussion.\n3. The study is descriptive, owing to the small number of patients studied. This in itself shouldn’t\nbe a problem, but broader conclusions can’t be drawn from this small set of samples. For example,\nwhile the association with JAK-STAT pathway in lymph node metastasis is interesting, it can’t be\ngeneralised.\nResponse: We appreciate the reviewer’s point and agree that our study, based on a limited\nnumber of patients, is inherently descriptive. While we believe that the observations, such as the\nassociation of JAK-STAT signalling with lymph node metastasis, are interesting and hypothesis-\ngenerating, we fully acknowledge that they cannot be generalized at this stage.\nTo address this, we have revised the relevant sections of the manuscript to clarify that our findings\nare exploratory and may not reflect broader trends without further validation in larger patient\ncohorts. We have also adjusted the language used in our conclusions and explicitly stated that\nfuture studies are needed to assess the generalisability of our findings.\n4. The authors state that “To date, our understanding of prostate cancer evolution has primarily\nbeen derived from bulk genome sequencing studies. These studies typically rely on single\nbiopsies from the primary tumour, which may not fully capture the extent of genetic heterogeneity\npresent within the entire tumour.” – but I believe this is not the case, there have been several\nstudies (dating back to Cooper et al., Nature Genetics 2015) that used multi-region sequencing\nto study prostate cancer evolution.\n- Hong et al https://www.nature.com/articles/ncomms7605\n- Warner et al https://www.nature.com/articles/s43018-023-00692-y\n- Singhal et al https://www.nature.com/articles/s41467-024-48629-y\n- Rao et al https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-024-01302-x\n- Nurminen et al https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-023-\n01242-y\n- Boutros et al https://www.nature.com/articles/ng.3315\nResponse: Thank you for pointing out these relevant publications. We acknowledge that the\nexisting literature was not properly cited, and we have now included a new paragraph in the\nintroduction, where we properly cite this relevant literature.\nReviewer #1 (Remarks on code availability):\nCode is available to reproduce the figures and appears well documented.\nReviewer #2\nexpertise in prostate cancer tumour heterogeneity, lineage plasticity and single cell sequencing\n(Remarks to the Author):\nThe manuscript “Clonal Evolution and Transcriptional Plasticity Shape Metastatic Dissemination\nRoutes in Prostate Cancer” by Mikutenaite et al. describes the clonal evolution of 5 prostate\ncancers in relation to the transcriptional changes. The study includes multiple well-portrayed\nexamples about why single-sample WGS biopsies are not enough to depict the total breadth of\nclonal evolution in prostate cancer.\nMajor issues:\n1. Presented work needs to be tied more clearly with current literature and put into context of state\nof the art. The manuscript refers to a limited number of studies in the introduction and discussion.\nFor example, several studies report metastatic seeding patterns from locally advanced prostate\ncancers based on spatially distinct genome sequencing. These works should be acknowledged\nmore clearly and putative novelty of the presented work articulated based on complete literature.\nResponse: We agree with the reviewer and apologize for the brevity of our prior introduction (this\npoint was also raised by reviewer #1. Please see response 1.4). In response, we have revised\nand expanded the introduction, properly citing the relevant literature and providing a proper\nbackground to the study.\n2. Transcriptional plasticity or the results of other prostate cancer single cell transcriptomics (or\nspatial transcriptomics) papers are not addressed in introduction. Mentioning generally “signalling\npathway activities” does not tell much. On the other hand, e.g. in discussion they mention “a\nrecent study identified four major groups of androgen-independent transcriptional programmes\n(28)”, but do not elaborate in relation to their results.\nResponse: Thank you for pointing out the missing context. We have now expanded the\nintroduction and appropriate discussion section, which is particularly interesting and relevant to\nour study due the orthogonal finding of a WNT-subtype.\n3. As there is WGS data on every patient (and sample?), could you report the driver mutations in\nthese cases, and how they are related to the inferred clones? Why are they not involved in tumor\nevolution reconstruction?\nResponse: We thank the reviewer for the constructive suggestion. Although prostate cancer is\noverall considered a “C-type” cancer, driven mainly by structural variations and somatic copy\nnumber alterations, driver point mutations such TP53 and SPOP (PMID: 34716314) are also\nimportant. We would like to clarify that we used shallow WGS to identify large CNAs concordant\nwith our scRNA-seq based CNA inference. The shallow WGS does, unfortunately, not allow us to\nidentify point mutation drivers. Instead, we performed deep WGS from a single area of each\ntumour, which is the basis of our single-biopsy WGS phylogenetic trees (see e.g. Figure 2). From\nthe deep WGS, we have identified both CNA, SV and SNV drivers, which we have now included\nin a new supplementary Table 3. We now also provide the CNA profiles from the shallow WGS of\nall the areas in a new Supplementary Figure 4. We have also further clarified our usage of shallow\nWGS in the revised manuscript.\n4. While the manuscript is well-written, the text is occasionally a bit difficult to follow, especially\nkeeping in mind which patient the results address, or if a clone in question has a special role, e.g.\nin seeding. Would it be possible to call the clones with patient name like PCAL03-C1 throughout\nthe text?\nResponse: We thank the reviewer for the good suggestion to improve the clarity of the complex\nclonal and transcriptional relationships. In response, we have revised the manuscript to label\nclones with patient-specific identifiers (e.g., PCAL03-C1) throughout the text, to improve clarity\nand traceability.\nMinor issues:\n5. It should be stated in the abstract that these are locally advanced prostate cancers.\nResponse: Thank you for pointing this out. We have revised the abstract to clearly state that the\nstudy focuses on locally advanced prostate cancers.\n6. The Figure 1B is not clear. What are the smaller areas on the left and right side?\nResponse: Thank you for your valuable feedback. We have revised Figure 1B to improve clarity.\nSpecifically, we have updated the schematic to include annotations that more clearly show the\nsmaller areas on both sides of the pathologist’s report correspond to specific sectioning locations\nwithin the prostate. We believe this enhancement will help readers better understand the spatial\ncontext of the pathologist’s map.\n7. Is the PCAL37 LNM1 sample a metastasis? it seems to contain only normal cells based on\nFigure 1E.\nResponse: We thank the reviewer for the careful observation. PCAL37 LNM1 is indeed a lymph\nnode metastasis; however, as noted in the text, it contains predominantly normal (non-cancer)\ncells. To clarify this, we have now added the cell number for PCAL37 LNM1 in the main text and\na supplementary table listing all samples along with their cell type proportions and total cell\nnumbers (Supplementary Table 2).\n8. Authors should check the references to supplementary figures, they might be off in some cases,\ne.g. page 6.\nResponse: We thank the reviewer for the attention. We have carefully verified all references to\nmain and supplementary figures and tables.\n9. Not all data is shown, e.g. “subclonal shattering of chromosome 19 in this clade”\nResponse: We thank the reviewer for pointing out the lack of supportive data for this statement.\nWe have now provided the appropriate supportive data including copy number profiles from the\nshallow WGS of all the areas, including patient PCAL37, in a new Supplementary Figure 4.\n10. The WGDs are not marked in all trees in Figure S4.\nResponse: We thank the reviewer for bringing our attention to this mistake. We have now\nupdated Supplementary Figure 4 (new Supplementary Figure 5), and added the presence or\nabsence of WGD events consistently across all trees. This has also been clarified in the figure\nlegend.\n11. Why are not all tumor areas included in the evolutionary reconstruction of PCAL25?\nResponse: We thank the reviewer for addressing this point. We performed a careful quality\nassessment of all snRNA-seq data prior to our integrative analysis. We found that not all tumour\nareas from PCAL25 met our objective quality control criteria and were therefore excluded from\nthe clonality analysis. We have now clarified this point in both the main text and the Methods\nsection.\n12. In case PCAL03, the authors state that “The metastatic seeding occurred from a branching\nevent derived from areas T09 and T10”. However, this does not seem to be the case based on\nthe evolutionary tree in Figure S4E.\nResponse: We realise that the sentence was not clearly formulated, and have now revised the\ntext to “A branching event gave rise to one leaf clone seeding in areas T09 and T10 and another\nleaf giving rise to metastatic seeding.”.\n13. Is MSH2 expressed in the T02 area with biallelic deletion (case PCAL10?)?\nResponse: This is a very good point. In response, we performed differential analysis between\nT02 (MSH2 biallelic inactivation) and T09 (MSH2 diploid), which showed a substantial\ndownregulation of MSH2 (log2 fold change = -3.7), supporting the loss of MSH2 in this area.\n14. On page 8, last sentence, how does the copy number instability drive the evolution towards\nAR-negative lineage?\nResponse: We thank the reviewer for this interesting observation. While there is evidence for\ncopy number instability driving AR-negative cancers, such as neuroendocrine prostate cancer,\nwe acknowledge that the order of events remains to be determined. We have, therefore,\nrephrased the sentence to reflect this uncertainty, and now write “...suggest that AR-negative\ncopy number unstable lineages can drive the clonal evolution towards metastatic seeding”.\n15. Page 15, “tumour cell states” are mentioned for the first time, what are these states?\nResponse: We apologize for the lack of clarity and have now added a description and reference\nto the original definition in the revised manuscript.\n16. Is it really shown that the metastatic seeding occurred before the WGD in PCAL25? They are\nin different branches and there’s no timing analysis presented.\nResponse: We thank the reviewer for the insightful observation. Our clonal reconstruction shows\nthat the diploid clone PCAL25-C3 seeded the diploid PCAL25-C2 and PCAL25-C4 clones in\nLMN1. A WGD event occurred in another C3-derived branch (C9), which seeded C1 in LMN1 (Fig\n3A and F). However, as the reviewer correctly pointed out, we cannot deduce the relative timing\nof events, and therefore cannot exclude the possibility that the WGD-acquired clade (“C1”) was\nthe first metastatic seeding event. We have now revised the text to clarify that we can not deduce\nthe relative timing of events.\n17. The project clearly has generated a lot of results and it is understandable that not all are\nbrought up in the manuscript text or showed in the supplementary figures. However, some\nobservations that could be elaborated. For example, area T09 in PCAL34 has both diploid and\ntetraploid cells?\nResponse: Thank you for highlighting this interesting observation. We have now elaborated our\ndescription of different ploidy clades, including this interesting observation in the revised\nmanuscript and referenced it in Supplementary Figure 7C, where the presence of both diploid and\ntetraploid cells in area T09 of PCAL34 is shown.\n18. single biopsy tree has wrongly inferred a different seeding clone than single cell data, in case\nPCAL03.\nResponse: We appreciate the reviewer’s point, and we should have clarified the origin of the\nsamples better. In the case of PCAL03, the single-biopsy WGS-based phylogenetic tree indeed\ninferred a different seeding clone compared to the single-cell data. In the revised manuscript, we\nemphasize that while the primary tumour was the same between the single-sample WGS biopsy\nand the multi-region scRNA-seq analysis, the single-sample metastatic biopsy for deep WGS was\nobtained from a different lymph node than the one used for scRNA-seq. Interestingly, the single-\nbiopsy WGS-based phylogenetic tree inferred a different seeding clone, suggesting polyphyletic\nseeding. We have now clarified this sampling difference and revised the interpretation to make\nthis point clearer.\n19. low pass WGS results are not shown (Figure S2 and S3 legends say they are from deep\nsequencing)\nResponse: Thank you for pointing out the missing data. We have revised Supplementary Figure\n4, so that it now contains the copy number profiles from the shallow WGS of all the areas.\nReviewer #2 (Remarks on code availability):\nProvided code acts as a documentation of the performed analysis, but it is not directly executable\nin different environments due to e.g. specific library dependencies. As these are analysis scripts,\nI find this level of documentation sufficient.\nReviewer #3\nexpertise in prostate cancer omics and mechanisms of cell plasticity (Remarks to the Author):\nThe manuscript by Mikutenaite et al. investigated clonal relationships in five locally advanced,\nhigh-risk prostate cancer (PCa) patients using single-nucleus RNA sequencing (snRNA-seq)\nintegrated with low-pass whole genome sequencing (WGS). All patients underwent radical\nprostatectomy with lymph-node (LN) dissection; one patient (PCAL25) had received neoadjuvant\nandrogen deprivation therapy (ADT), while the others were treatment-naïve. The manuscript\neffectively comprises five case studies in which tumor samples from multiple spatial regions\nincluding LN metastases subjected to phylogenetic reconstruction. The use of snRNA-seq layered\non genomic phylogenies provides higher-resolution inference than standard bulk WGS and allows\nfor mapping of transcriptional states onto clonal architectures.\nWhile the genomic and phylogenetic findings themselves (e.g., examples of both monophyletic\nand polyphyletic metastatic seeding) are not novel in the context of prostate cancer, the\nintegration of transcriptomic data with clonal evolution is an interesting and valuable aspect. For\nexample, the observation that metastases tend to arise from clones with reduced androgen\nreceptor (AR) activity and enhanced expression of pathways associated with lineage plasticity is\npotentially clinically relevant. However, the manuscript currently lacks sufficient depth in both the\nanalysis and interpretation of these biological and clinical implications.\nWith limited number of patients including only one treated patient, the study is limited in the scope.\nFurther development of the analysis and more insightful contextualization of the transcriptomic\nfindings and addition of treated patients will strengthen the finding.\nMajor comments\n1. The authors should discuss how the findings for PCAL25 relate to the context of prior exposure\nto ADT. For example, the increased evidence of fibrosis in this patient is probably a consequence\nof inflammation and regression of the prostate in response to ADT. It is also surprising to see\nseveral clones with apparently high AR activity in a newly diagnosed PCa patient after 6 weeks\nof ADT – is this because pathway/activity scores are only relative to other clones in the same\npatient, meaning that some clones must always be considered “high” and others “low” in every\npatient? If samples from different patients were run in multiplex for snRNA-seq (please clarify this)\nto minimize batch effects between patients, pathway/activity scores should be made comparable\nacross patients if not already. Are the high Androgen/AR scores in PCAL25 similar to the high\nscores in other patients? Also, do the authors have any pre-ADT biopsies from PCAL25 that they\ncould add to their analysis?\nResponse: We thank the reviewer for the insightful comments and suggestions. The reviewer is\ncorrect that ADT is expected to lead to lower levels of AR activities and that the increased\nproportion of fibroblasts in PCAL25 could suggest a fibrotic response.\nTo clarify, the pathway activities are normalised across patients, making them directly\ncomparable, and we do observe a lower level of AR activity in this patient. We have now included\na new analysis, comparing the pathways between patient samples, presented in a new\nSupplementary Figure 14. We have also included a sentence describing the new analysis in the\nrevised manuscript. As PCAL25 is the only patient sample with prior ADT, we have deliberately\nrefrained from conducting any comparative analyses. Instead, we emphasise that differences in\ncell composition and pathway activities in this particular patient may be influenced by prior\ntherapy.\n2. It is difficult to picture how the samples in the various transcriptomic heatmaps line up with the\nphylogenies because the sample orders are different, and closely related clones do not always\nhave similar colours (e.g. in Fig 2c, there are blue clones in three different branches). Although\nheatmaps clustered by transcriptomic profiles/scores are informative and should be kept (possibly\nas Supplementary Figures), it would be more helpful to have the samples and dendrograms in\ntranscriptome-based heatmaps to be arranged and oriented as in the genomic phylogenies. This\nwould better illustrate whether closely related clones have similar transcriptomes and whether\nmore distantly related clones display convergent evolution.\nResponse: We thank the reviewer for the excellent suggestion to improve the clarity and\ninterpretability of our findings. In response, we have re-clustered the transcriptomic heatmaps,\nusing the phylogenetic structure instead of the hierarchical clustering. We believe that this\nalignment has improved the illustration of the relationship between genomic relatedness and\ntranscriptomic profiles.\n3. Most of the biologically interpretable transcriptomic data is presented at the resolution of clones.\nHowever, there is clearly evidence for transcriptomic heterogeneity within individual clones (e.g.\nFig 3F,H). The authors should explore intra-clone transcriptomic heterogeneity in more depth.\nFor instance, they could investigate how much transcriptomic heterogeneity exists within individual\nclones compared to variability between clones.\nResponse: We would like to thank the reviewer for the excellent and insightful suggestion. In\nresponse, we have performed a more detailed analysis of the transcriptional relationship between\nand within clones. For Figure 3H, for example, we have added additional support for clonal\ntranscriptional heterogeneity in our new Supplementary Figure 12, which adds further support to\nthe transcriptional heterogeneity between clones. We have also performed a new analysis,\ninvestigating the transcriptomic clonal heterogeneity at the pathway level, which is now presented\nin a new Supplementary Figure 15.\n4. The authors should discuss whether transcriptional profiles reflect key genomic aberrations.\nFor example, do clones with gain of MYC or FOXA1 have higher activity scores for those\ntranscription factors? Do losses of specific tumour suppressors, e.g. TP53, associate with\nactivation of plasticity pathways? The authors say in the Discussion that they “observed tumour\nareas containing clones with homozygous MSH2 loss to have higher lymphocyte infiltration”,\nbut this isn’t presented anywhere in the Results despite being an important, albeit expected,\nobservation. And so on.\nResponse: We thank the reviewer for the interesting and constructive suggestion. In response,\nwe have expanded our analysis comparing somatic genomic and transcriptional correlates (see\nnew Supplementary Figure 16). This includes major genomic events such as chromosome 8q\namplification (MYC gain), 21q loss (ERG), 13q loss (RB1), 17p loss (TP53). Our expanded\nanalysis supports our initial analysis (presented in our PCA analysis), with a consistent\nrelationship between areas with MSH2 homozygous loss and high lymphocyte infiltration. In the\nnew analysis, we also find support for immune cell depletion in areas with TP53 loss or MYC gain,\nin agreement with prior literature.\nWe note that our analysis is limited to analysing large copy number alterations, which precludes\nanalysis of driver point mutations in, for example, FOXA1, which is primarily hotspot base\nsubstitutions and focal duplications.\n5. The authors need to provide more details about the single-biopsy bulk WGS. In each patient,\nwhich tumour site was that biopsy taken from? What were the methods for the analysis? Etc.\nResponse: We thank the reviewer for pointing this out, and we apologize for the lack of clarity in\nour original manuscript. In response, we have expanded the description of the single-biopsy bulk\nWGS sample and the method used for copy number segmentation and clonal reconstruction in\nthe revised manuscript.\n6. The phylogenetic trees from single-biopsy WGS and multiregional snRNA-seq for PCAL03 (Fig\nS4E) show completely different clones being inferred as the clone that seeded the LN metastasis,\nwhile truncal whole-genome doubling was inferred in the single-biopsy WGS-based tree but\nnowhere in the multiregional tree. How do the authors explain these discrepancies, and do they\naffect biological interpretation of results?\nResponse: We thank the reviewer for bringing our attention to this special case, which was\ninsufficiently explained in the initial manuscript. PCAL03 was a special case, where the lymph\nnode metastasis from the single biopsy sample was, in fact, different from the lymph node\nmetastasis sampled for our multi-region single-cell analysis. We have expanded the description\nand explanation in the Results section and Methods, lending support to our findings of polyphyletic\ndissemination routes. We have also added a supportive whole-genome doubling mark on the\nmultiregional tree.\n7. The authors should discuss how their work fits within the broader field of prostate cancer\ngenomics.\nResponse: We thank the reviewer for the constructive suggestion. We have now included a\nperspective on how our findings may fit and impact the field in the Discussion of the revised\nmanuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of the points raised by me. However, I am still not convinced that the data shown here is evidence for metastatic 'seeding pockets' - I believe that such a claim, even in its softer version after the changes to the text, requires a more thorough investigation. For example, let us take Fig 3B. T10 and T05 harbour 'metastatic' clones (i.e., clones that are also found in the lymph node mets. However, there is also plenty of tumour extending into the central region at those Z levels of the prostate. If you sequence such regions, and you find that the 'metastatic' clones can be found in those regions too, where is the metastasis being seeded from? I see the fundamental issue as _definite_ phylogenetic relationships being extrapolated to a _potential but by no means definite_ spatial relationship. It is valid to make a claim about a clone being the seeding clone (let us say clone A) because even if deeper sequencing or sequencing from more regions shows that there is a subclone A' that is in fact the 'real' seeding clone, the original claim (that clone A is the seeding clone) broadly holds true because A' is a subset of A. However, the same cannot be done for spatial relationships without evidence from more regions - claiming that region X is the seeding region would mean all other regions (which haven't been profiled) are not seeding regions. Apologies if I am missing something here, I remain open to be convinced.", "reviewer_2": "Authors have address all the issues I raised and I'm happy with their revision.", "reviewer_3": "Although the authors responded to all reviewer comments."}, "author_response": "Nature Communications NCOMMS-25-31194.\nRebuttal version 2\nReviewer #1 (Remarks to the Author):\nThe authors have addressed most of the points raised by me. However, I am still not convinced\nthat the data shown here is evidence for metastatic 'seeding pockets' - I believe that such a claim,\neven in its softer version after the changes to the text, requires a more thorough investigation.\nFor example, let us take Fig 3B. T10 and T05 harbour 'metastatic' clones (i.e., clones that are also\nfound in the lymph node mets. However, there is also plenty of tumour extending into the central\nregion at those Z levels of the prostate. If you sequence such regions, and you find that the\n'metastatic' clones can be found in those regions too, where is the metastasis being seeded from?\nI see the fundamental issue as _definite_ phylogenetic relationships being extrapolated to a\n_potential but by no means definite_ spatial relationship. It is valid to make a claim about a clone\nbeing the seeding clone (let us say clone A) because even if deeper sequencing or sequencing\nfrom more regions shows that there is a subclone A' that is in fact the 'real' seeding clone, the\noriginal claim (that clone A is the seeding clone) broadly holds true because A' is a subset of A.\nHowever, the same cannot be done for spatial relationships without evidence from more regions\n- claiming that region X is the seeding region would mean all other regions (which haven't been\nprofiled) are not seeding regions. Apologies if I am missing something here, I remain open to be\nconvinced.\nResponse: We appreciate the reviewer’s comment to refine our interpretation. We agree that\nother areas of the primary prostate may indeed also represent similar seeding zones. Our findings\nare based on our comprehensive multi-region sampling, where we present the spatial\nrelationships that support certain areas to be putative seeding zones, relative to other regions. To\nmake this point clearer and to address the reviewer’s comment, we have changed our wording in\nthe abstract and discussion to reflect the point raised by the reviewer.\nWe have revised the statement in the abstract:\nThese findings imply spatially distinct areas with different clonal and TME compositions, and\nsuggest that seeding clones, phylogenetically closely related to the metastasis, are localised in\nspatially confined regions in the primary tumour.\nThe relevant Discussion section:\nWe find tumour areas that are phylogenetically closely related to the metastasis to exhibit\ndistinct transcriptional and genomic characteristics, including androgen independence.\nThis suggests that seeding clones are located in spatially confined regions, often located\ntoward the periphery of the prostate. Despite our comprehensive multi-region sampling\napproach, we acknowledge that the spatial resolution of our sampling is, nevertheless,\nlimited, and that the presence of metastatic clones in a given region does not exclude\ntheir existence elsewhere in the tumour.\nAuthors have address all the issues I raised and I'm happy with their revision.\nReviewer #2 (Remarks on code availability):\nCore provides documentation for reproducing the key results.\nReviewer #3 (Remarks to the Author):\nAlthough the authors responded to all reviewer comments.\nReviewer #3 (Remarks on code availability):\nAlthough the number of the patient is limited, the authors responded adequately to the reviewer\ncomment.\n2"}]

Decision letter