Prostate cancer cells converge to an inflammatory-like state upon metastatic dissemination

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

Abstract
Identifying drivers of cancer progression to guide treatment selection is hindered by our limited understanding of tumor heterogeneity and its impact on tumor evolution. Here, we delineate the phenotypic variability across ~300,000 cells collected from multiple tumor loci in primary prostate and matched locoregional metastases using single-cell chromatin accessibility and gene expression sequencing. We find inter-patient heterogeneity to be confined to malignant populations. Within individual tumor loci, we see phenotypic heterogeneity among malignant cell populations despite a shared clonal genotypic architecture. We also observe that malignant cell populations disseminating to locoregional lymph nodes mirror the clonal architecture and phenotypic heterogeneity across primary tumor loci, while shifting from canonical prostate-cancer states to non-canonical inflammatory-like states. Our findings suggest a bottleneck imposed during the dissemination process, funneling prostate cancer cells toward an inflammatory-like cell state. These insights into the interplay between phenotypic identity and clonal architecture refine our understanding of prostate cancer progression and suggest that convergence of cancer cells towards an inflammatory-like state underlies dissemination to lymph nodes, offering a critical framework for future studies into prostate cancer metastatic potential.

Introduction
Understanding how tumor heterogeneity impacts cancer evolution has been a significant obstacle in our efforts to identify the drivers of cancer development and progression. This is particularly true for high-risk prostate cancer, with a 48–80% likelihood of progression within 5 years compared to overall survival rates of 90–99% in lower-risk groups1,2. This disparity underscores a need to uncover the dimensions of tumor heterogeneity that fuel tumor progression, inclusive of inter-patient, intra-tumor, intra-tissue and inter-tissue.

Most studies to date have relied on investigating single primary or metastatic biopsies using bulk sequencing, guiding our understanding on the existence of intra- and inter-tissue heterogeneity. For example, bulk genetic sequencing has showcased distinct clonal populations residing within different regions of the same tissue to harbor unique genetic aberrations, which could be indicative of outcome and potential treatment sensitivities3–9. Additionally, non genetic sources of variability, such as epigenetic reprogramming and transcriptional variability, play critical roles in reconfiguring cellular states to adapt to adverse conditions of metastasis and treatment resistance10–17. Phenotypic changes are observed throughout the onset and progression of prostate cancer to mCRPC, where non genetic variability is showcased with prostate-specific transcription factors (TFs) recruitment to chromatin variants to drive de-differentiation and promote progression9,18–26. Collectively, these findings underscore the necessity to characterize phenotypic diversity within and across tumors to fully understand the functional impact of heterogeneity towards cancer progression and treatment response.

In this work, we tackle the complexity of the multiple dimensions of tumor heterogeneity in prostate cancer progression by integrating multi-regional sampling of primary and metastatic sites with high-resolution single-nucleus chromatin accessibility and transcriptomic profiling to characterize the multi-faceted nature of heterogeneity and its role in progression. We find that inter-patient variability is restricted to malignant populations, that malignant cells within individual loci exhibit phenotypic divergence despite shared clonal origins, and that metastatic dissemination involves an inflammatory-like malignant state, revealing a bottleneck that funnels prostate cancer cells toward this noncanonical phenotype.

Results
Inter-patient heterogeneity is restricted to malignant epithelial as opposed to normal tumor microenvironment populations
We characterized 38 prostate cancer specimens derived from five high-risk human patients with ISUP ≥ 3 pathology who underwent radical prostatectomy and pelvic lymph node dissection at surgery as part of the International Cancer Genome Consortium (ICGC) late-onset cohort27 (Supplementary Dataset 1 and Fig. 1a). All patients were treatment-naïve before prostatectomy with pelvic lymphadenectomy, except one (ICGC_PCAL25), who received hormonal therapy 6 weeks before surgery. Using single-nuclei multiome sequencing (scMultiome- gene expression and chromatin accessibility from RNA- and ATAC-seq, respectively), we characterized 282,956 nuclei, with a median of 6 primary tumor loci from the prostate gland and 3 locoregionally lymph node metastases across all patients. Here, we use the term “loci” to refer to tumors sampled from spatially distinct sites of the prostate or across lymph nodes.

Leveraging both RNA and ATAC-seq data from all samples with a weighted nearest neighbor model28, 41, populations were identified on the phenotypic characteristics across all cells (Fig. 1b and Supplementary Fig. 1a). Cell type annotations identified 30 malignant epithelial populations that accounted for 86.4% of the cells. The remaining populations included major tumor microenvironment (TME) cell types, such as fibroblasts, macrophages, B and T cells (Fig. 1c). The cell types we identified displayed expected gene expression markers, such as AR, KRT18, and KRT8 in malignant epithelial cells and PTPRC in lymphocytes (Supplementary Fig. 1b). Additionally, we computed transcription factor DNA motif enrichments using ChromVar, which is the deviation score relative to an average chromatin accessibility signal computed across DNA motifs29. We observed significant TME-specific transcription factor DNA recognition motifs exemplified by the enrichment of the ETS1 motif in lymphocytes, the CEBPB motif in macrophages and the depletion of these motifs from accessible chromatin regions in malignant epithelial cells (Fig. 1d). Similarly, DNA recognition motifs of prostate-specific TFs, namely FOXA1 and HOXB13, were depleted in TME populations (Fig. 1d).

To quantify “inter-mixing” of populations, we utilized the Shannon entropy metric (H). We refer to population inter-mixing where populations include cells that come from multiple patients (Hp) or tissue site (Ht). In all cases, an entropy score of 0 indicates no inter-mixing with respect to the variable of interest. Conversely, maximum Hp or Ht reflects a cell population where cells from all patients or tissue sites contribute to the population’s composition with equal probability. We found that individual malignant epithelial populations were mainly composed of cells originating from a single patient, with a median Hp=0 across malignant epithelial populations, indicating no patient inter-mixing (Fig. 1e). Fibroblasts and macrophages were the most inter-mixed TME populations (Hp=2.15 and 1.5), with fibroblasts being closest to the maximum entropy value possible (Hpmax=2.32). Similarly, we identified patient inter-mixing in the basal epithelial population (Hp=1.2) with no inter-mixing according to tissue site (primary tumor vs lymph node; Ht=0) (Fig. 1e and Supplementary Fig. 1b–d). This agrees with prostate basal cells being a unique cell type to the prostate gland, and showcases the validity of biological signals revealed by single-cell multiome. Lastly, B cell populations were the only normal cell type showing limited patient inter-mixing, with median Hp=0.24 across B cell populations (Supplementary Fig. 1d), representative of patient-level phenotypic diversity previously reported in B cells30, while having patient-specific driver alterations (Supplementary Fig. 1e). Overall, these observations suggest biological constraints on most TME cell populations towards a shared identity across different patients. In contrast, malignant populations are patient-specific and suggest patient-dependent tumor evolutionary paths, warranting analyses on a patient-by-patient basis.

Malignant cell populations exhibit intra-locus heterogeneity, yet are phenotypically similar within a prostate gland
Single primary tumor biopsies do not capture the diverse histology and genetics reported within the prostate gland and are unlikely to reflect the phenotypic diversity of prostate cancer cells in a patient31. Thus, to better understand phenotypic heterogeneity within a prostate but across tumor loci from spatially distinct locations, we analyzed the prostate gland on a patient-level basis (3–8 tumor sections per patient) (Fig. 2a and Supplementary Fig. 2). A median of 88% of cells were assigned to malignant epithelial populations across patients (Fig. 2a). Malignant epithelial populations exhibited almost no inter-mixing of cells from different loci according to sample entropy measurements in contrast to TMEs (median normalized Hl=0.06 versus Hl=0.8 across patients, Mann–Whitney U test P=3.05 × 10−6) (Fig. 2b). A median of 3 malignant epithelial populations contained more than 90% of cells from the same tumor locus (Supplementary Fig. 3). Inter-mixing across loci in malignant epithelial populations was skewed to those engaged in cell division (Fig. 2c) (median normalized Hl=0.80) suggesting that cell cycle signal can in part override locus-specificity of most malignant epithelial populations and drive inter-mixing of cells from different loci. Overall, the low degree of inter-mixing across loci replicated in each patient suggests a detectable degree of phenotypic heterogeneity in cancer cells across loci within the prostate gland.

To further assess the heterogeneity across malignant epithelial populations, we looked at correlations across populations based on their ATAC or RNA read-outs. Pan-cancer analysis of ATAC- and RNA-seq data reported chromatin accessibility at non-promoter elements to reflect cancer types with increased specificity compared to accessibility at gene promoters or gene expression levels32. Hence, we generated a patient-specific repertoire of chromatin-accessible DNA elements based on the normalized ATAC signal at promoters and non-promoter elements in each population. Pearson correlations had a wider range of values across non-promoter elements (median: range=0.96, N elements=120,942) compared to promoter elements (median: range=0.44, N elements=38,545) and transcripts (median range=0.13; N transcripts=36,601) (Supplementary Fig. 3b–d). These results highlight the specificity of chromatin accessibility at non-promoter elements as a key source of variation across populations, accounting for intratumoral heterogeneity.

Accessible non-promoter elements can serve to regulate gene expression through the recruitment of TFs recognizing the DNA recognition motif they harbor. We reasoned that the variance in TF-DNA recognition motif enrichments across the patient-specific repertoire of accessible DNA elements can reveal sources of heterogeneity in malignant populations. Principal component analysis in each patient revealed motifs for the FOS and JUN-related TFs, recently linked to stem cell-like properties in castration resistant prostate cancer26, consistently emerging as top contributors to the major axis of variation across malignant epithelial populations in all patients (Fig. 2d and Supplementary Fig. 4). Nuclear receptor, Forkhead, and GATA transcription factor families were also found as strong sources of variation, in line with their role as drivers of prostate cancer progression33. Variance across malignant epithelial populations was associated with neuroendocrine-like TFs, such as ASCL1, SNAI1, and TCF4, supporting a neuroendocrine-based source of heterogeneity within one patient (ICGC_PCAL34)26,34,35. Additionally, we observed PAX-family TF motifs, such as PAX6, which is implied in neuro-ectoderm development but not well characterized in prostate cancer36. To further reproduce and contextualize this, we linked TFs, accessibility, and gene expression using single-cell enhancer-gene regulatory networks analysis37 on tumors from patients ICGC_PCAL34 and ICGC_PCAL03. The highly specific TFs to individual malignant epithelial clusters included FOX, HOX, FOS/JUN, GRHL, and NR-based TFs, aligned with our orthogonal analysis (Fig. 2d and Supplementary Dataset 2). Performing pathway enrichment analysis on the lists of genes specific to a given cluster, regulated by the reported TFs, revealed distinct and overlapping pathways (FDR<0.05). The commonly enriched pathways included androgen and estrogen response, whereas more specific pathways involved EMT, hypoxia, and TGF-β signaling in ICGC_PCAL34, along with p53, TGF-β, PI3K/AKT/mTOR Signaling, mTORC1 and DNA repair in ICGC_PCAL03 (Supplementary Fig. 5 and Supplementary Dataset 2: Hallmark Pathways). Overall, these results suggest that prostate lineage-specific TFs, along with JUN/FOS, differentially engage with accessible chromatin, driving heterogeneity across malignant epithelial populations.

To contextualize malignant epithelial population heterogeneity, we computed Pearson correlations across malignant epithelial and TME populations. TME-to-TME comparisons offer a reference point for expected variability between well-defined cell types. Across patients, malignant epithelial populations were more correlated with each other than TME populations were among themselves across patients (Fig. 2e) (Mann–Whitney U test P: ICGC_PCAL10=8.77 × 10−17, ICGC_PCAL03=0.02, ICGC_PCAL25=1.83 × 10−09, ICGC_PCAL34=3.76 × 10−09, ICGC_PCAL37=1.21 × 10−06). Thus, while malignant populations exhibit distinct regulatory programs driving their heterogeneity between and within tumor loci, their core molecular identity is more preserved than what is observed across populations from the TME.

Locus-specific malignant epithelial populations are composed of similar clonal architectures
High-risk prostate cancer is known for its genotypic diversity with a high frequency of structural variants and copy number variations (CNVs)38. For this reason, we inferred CNV profiles from the scMultiome data using NUMBAT39 to link malignant epithelial populations to genetic clones within a patient. We detected typical prostate cancer CNV events such as chromosome 6 and 8p loss, and 8q gains (Supplementary Dataset 3), with notable variation across patients in line with prior reports40,41. The genetic architecture of each patient included a median of 6 clones per patient (Fig. 3a) with varying degrees of similarity between malignant clones. We measured similarity between clones using the Hamming metric, where 0 represents the same CNV label across all genomic segments between groups, and 1 indicates no shared CNV states over segments called per patient.

Patients ICGC_PCAL03 and ICGC_PCAL37 had the smallest difference between clones (mean distance between clones: ICGC_PCAL03=0.13, ICGC_PCAL37=0.06), whereas ICGC_PCAL10, ICGC_PCAL34, and ICGC_PCAL25 were more diverse (mean distance between clones: ICGC_PCAL10=0.39, ICGC_PCAL34=0.22, ICGC_PCAL25=0.3) (Supplementary Fig. 6a). In the latter patients, the most distinct clones contained cells from a single tumor locus, suggesting genotypic diversity to be locus-restricted. For instance, ICGC_PCAL10’s clone 3 was the most distinct, with a mean distance of 0.9 from other clones, where 98.3% of clone 3 cells were from the Tu5 locus (Fig. 3b, c). This clone was the exception in ICGC_PCAL10, as the rest of the clones in ICGC_PCAL10 were similar to each other (mean distance=0.08). Similarly, in ICGC_PCAL34, clones 4 and 6 accounted for 90% of all Tu6 locus malignant epithelial cells used in CNV inference (Fig. 3b, c). In ICGC_PCAL25, clone 4 and its subclones (5–8) were largely composed of Tu6 locus malignant epithelial cells, averaging 87% (Fig. 3b, c). Thus, the most distinct clones relative to the most common patient CNVs were restricted to a single sampled tumor section and not found in multiple loci.

Next, we wanted to understand the mapping between phenotype-based populations and genotype-based clones. We used the Jensen–Shannon divergence score to compare how malignant populations relate to each other with respect to the overlap in the probability distribution of their genotype-based clonal membership. Specifically, a score of 1 indicates that two malignant populations are composed of the same underlying sets of clones with identical proportions. Conversely, a JSD score of 0 represents no overlap in the probability distributions of clones between a pair of malignant populations. We found that despite the intra-tissue and intra-locus specificity of malignant epithelial populations based on chromatin and expression data, these populations had similar clonal architectures (median clonal distribution distance between populations across patients=0.07) (Fig. 3 and Supplementary Fig. 6b). The similar clonal composition of phenotypically heterogeneous malignant epithelial populations from multiple tumor loci further supported a population diversity relying on nongenetic variation embedded with chromatin accessibility and transcriptional programs.

Malignant cells disseminating to locoregional lymph nodes undergo an inflammatory-like cell state transition despite genotypic similarity to primary cancer
Locoregional lymph nodes are the earliest and most common site of malignant epithelial cell disseminations in prostate cancer and are linked to worse prostate cancer outcomes2,42. Therefore, their involvement offers crucial insights into tumor heterogeneity and its role in the evolution towards a metastatic state. For this reason, we utilized multiple locoregional lymph nodes collected in four of four prostate cancer patients (median 3 metastatic lymph nodes per patient) to characterize their cell state composition from scMultiome assays (Fig. 4a and Supplementary Fig. 7). Except for patient ICGC_PCAL34, lymph node samples had a similar clonal architecture to each other and to their matched primary tumor loci (Fig. 4b and Supplementary Fig. 7). Additionally, locus-restricted intraglandular clones, such as clone 3 in ICGC_PCAL10, were not detected in multiple matched lymph nodes (Fig. 4b and Supplementary Fig. 7). From a phenotypic angle, malignant epithelial populations continued to show low inter-mixing of tissue sites and samples (Fig. 4c), with median Ht=0 across all patients’ malignant epithelial populations. ICGC_PCAL37, which had the lowest Gleason score (7) was an exception, with malignant populations being most inter-mixed across tissues in our cohort (median Ht=0.4). To quantify the degree of phenotypic differences, we compared Pearson correlations of non-promoter DNA accessibility profiles between malignant epithelial populations containing >90% of cells from a single locus in either the lymph node or the prostate. In most patients, correlation ranges for inter-tissue comparisons (lymph node-specific vs prostate-specific) were not significantly different from intra-tissue comparisons (prostate-only), indicating that dissemination did not introduce a greater range of variability between populations than already present within the primary tumors (Mann–Whitney U test, P: ICGC_PCAL34=0.47, ICGC_PCAL10=0.28, ICGC_PCAL03=0.06, ICGC_PCAL37=6.78 × 10−13) (Fig. 4d). Notable exceptions included ICGC_PCAL34, where lymph node-specific populations were more similar to each other than to prostate-specific populations (P=0.046), and ICGC_PCAL37, which had only a single lymph node-specific population, precluding intra-tissue correlation estimates. Clinically, ICGC_PCAL37 had no evidence of bone metastases and remained stable for nearly 5 years before initiating systemic treatment (Supplementary Dataset 1), consistent with reduced molecular heterogeneity in lower Gleason score disease. Together, these data indicate that the stability of patient-level clonal architecture, paired with comparable intra- and inter-tissue phenotypic correlations, points to a role for phenotypic variability that may account for the tissue-specificity of populations, changes in cell state not explained by genetics alone and are likely induced by the microenvironment, as malignant prostate cancer cells adapt to colonize lymph nodes.

Given the genotypic similarity but phenotypic-specificity of lymph node malignant epithelial populations, we wanted to understand if malignant cells in the lymph nodes are undergoing a coordinated phenotypic shift across patients. Thus, we extended the comparison of malignant epithelial cells at the tissue level (lymph node versus primary tumor) across all patients to determine the biological basis of disease progression. Differential gene expression analysis on malignant cells identified statistically significant upregulated pathways in lymph node tumors to be associated with the immune system, such as allograft rejection (normalized enrichment score (NES)=2.86), inflammatory response (NES=2.05), IL6 JAK-STAT3 signaling (NES=2.02), as well as the epithelial-mesenchymal transition (NES=2.27) (Supplementary Dataset 4). In contrast, pathways upregulated in primary tumor malignant epithelial cells associated with the cell cycle (E2F targets (NES=−1.71), G2M checkpoints (NES=−1.56)) and hormone-driven cancers (cholesterol synthesis (NES=−1.54) and estrogen response (NES=−1.46)) (Fig. 4e). To ensure the differential signal observed across tissues is unlikely to be due to ambient RNA, we applied ambient RNA correction on our samples43. Ambient RNA correction was marginally helpful in increasing TME-related signal strength, which belonged to cells that were excluded in the differential expression analyses (Supplementary Fig. 8a, b). Repeating differential gene-expression analysis on all patients using ambient RNA corrected results revealed the same trends, and the majority (8/11) of the GO terms were in agreement (Supplementary Fig. 8c).

Analysis of DNA recognition motif enrichment on differentially accessible chromatin regions detected the motifs for the ETS, IRF and RUNX TFs as the top statistically significant motifs in lymph node malignant epithelial cells (Supplementary Dataset 4). These immune-related TFs44,45 were implicated in immune evasion of melanoma cancer cells in the lymph nodes46. In contrast, nuclear receptors and Forkhead DNA recognition motifs, such as those recognized by the androgen receptor (AR) and FOXA1, respectively, were enriched in differentially accessible regions of primary tumor-specific malignant epithelial cells (Supplementary Dataset 4 and Fig. 4f). We further validated the tissue-specificity of phenotypic changes in malignant epithelial populations on an intra-patient level. We computed the expression signature of the statistically significant genes per tissue at the single cell level and ranked all populations based on their median enrichment37. When using lymph node-upregulated genes, lymph node-specific malignant epithelial populations ranked higher than primary tumor-specific malignant epithelial populations across all patients, except for our lowest Gleason Score patient (ICGC_PCAL37) (Supplementary Fig. 9a). We observed the opposite trend for differentially upregulated genes at the primary tumor site (Supplementary Fig. 9b). In addition, fibroblasts, then immune populations ranked the highest in lymph node differential gene enrichment and surpassed all malignant epithelial populations (Supplementary Fig. 9a). These findings demonstrated inter-patient commonalities in metastatic versus primary malignant epithelial cell states, despite significant inter-patient genotypic variabilities. Taken together, malignant epithelial cells found within locoregional lymph nodes have a similar clonal architecture to the primary gland, suggesting limited constraint in the dissemination of clones. Furthermore, the malignant epithelial cells from all patients undergo shared phenotypic change toward an immune/inflammatory-like state across tissues.

Our own cohort precluded contextualizing these gene sets with respect to drug treatment, as the only androgen deprivation therapy (ADT)-treated patient (ICGC_PCAL25), the LN samples had insufficient quality for downstream single-cell profiling. Thus, we turned to external single-cell RNA sequencing data from metastatic prostate cancer47. This dataset covers major treatment stages of pre-treatment, on-ADT and/or anti-PD-1 combination therapy and recurrence in metastatic settings. Focusing on the relative enrichment across treatment stages of our signatures above, showcased LN-up genes to increase over treatment relative to pre-treatment (Fig. 4g). In our model, this suggests that transcriptional programs associated with lymph node metastases can emerge as the disease progresses, suggesting that an immune/inflammatory-like state may confer advantages beyond the lymph node microenvironment/metastatic niche. Concomitantly, primary prostate-upregulated genes, representing canonical prostate pathways, were also relatively enriched in post-treatment recurrence, which aligns with a renewed reliance on oncogenic AR signaling during disease relapse (Fig. 4g). Together, these results support a model of parallel pathway exploration within the primary site along the axis of progression in prostate cancer, reflecting coexisting adaptive strategies.

Methods
Biospecimens and quality control
Patients. Biospecimens were collected from five patients diagnosed with prostate cancer as part of the German ICGC cohort. Written informed patient consent and an ethical vote (Institutional Reviewing Board) were obtained according to the ICGC guidelines (see http://www.icgc.org). The local ethics board commission in Hamburg approved the study (PV3552 and PV4679), and the study was conducted following ethical guidelines (Declaration of Helsinki). No statistical methods were used to determine sample sizes since this study was exploratory.

Inclusion criteria were: primary prostate cancers, high-risk, histologically adenocarcinomas, multiple biopsies from a locally advanced specimen with a unifocal disease, and matching pelvic lymph node metastases.

Patients did not receive any neo-adjuvant radiotherapy, ADT, or chemotherapy before the prostatectomy with pelvic lymphadenectomy, except for one patient (ICGC_PCAL25) who received neoadjuvant hormone therapy with LH-RH agonist 6 weeks before surgery. One patient (ICGC_PCAL03) had a known oligometastatic disease stage at the time of operation (M1b stadium in conventional imaging with bone metastases).

All prostate specimens were analyzed according to a standard procedure, including a complete embedding of the entire prostate for histological analysis72. Clinical metadata were retrieved from the patient’s records, and follow-up data were collected with questionnaires.

Pathology review
Specialized pathologists dissected each prostate immediately after surgery as described before27. In brief, dissection followed a pre-defined scheme to represent the position of each block relative to the entire prostate. This procedure resulted in 60–150 pieces of tissues depending on the size of the prostate. After dissection, each tissue block was placed on a separately labeled cork plate covered with a special compound for cryopreservation (OCT) before the tissue was frozen to −20 °C and subsequently stored at −80 °C. Accordingly, the dissection of the lymph node packages proceeded. For each patient, 10 distinct tumor biopsies were taken, and matching lymph node metastases were included in this study (Supplementary Dataset 1). In total, we had a sample collection of 66 samples from 5 different high-risk locally advanced prostate cancer patients (50 primary gland biopsies and 16 matching lymph node biopsies). Gleason grading per individual biopsy was reported (Supplementary Dataset 1). Tumor cell content exceeded 80% for each sample according to pathology reviews.

Sample and library preparation for 10x single-cell multiome ATAC and gene expression
Single-cell sequencing can be prone to batch effects. Given the complexity of our cohort (multiple patients, samples, and tissues), we implemented the following strategy to ensure we are minimizing potential confounding variables of wet-lab library preparation batches and sequencing batches outlined in (Supplementary Dataset 1). Most library-preparation batches (4 samples at a time) incorporated samples from a distinct patient while including randomly chosen samples from another patient. We aimed to randomize samples derived from the same tissue sites across batches as well, such that not all metastatic samples were processed together. This approach was chosen to reduce library preparation batch effects and to better reveal biological variability from technical noise in downstream analyses.

Nuclei isolation
The nuclei were isolated following the adapted protocol for clinically relevant frozen tissues described previously73. We kept the samples on dry ice until use and weighed them individually before starting. The range in weight was between 7 and 24 mg. To perform single-cell multiome ATAC+gene expression analysis, we added 0.2 U/µl of Protector RNAse-Inhibitor (Sigma Aldrich PN-3335399001) to the homogenization buffer (according to the 3’ GEX assay). In brief, we thawed the frozen tissue in a pre-chilled douncer with 1 ml homogenization buffer for 5 min. We then dounced with pestle A for 10 strokes on ice and continued with pestle B for another 10 strokes. We next filtered the solution through a 40 µm tip into a prechilled 1.5 ml Eppendorf tube and spun the nuclei (5 min, 500 × g, 4 °C) in a fixed-angle centrifuge. Then, we removed all but the pellet with 50 µl of supernatant and resuspended it in 350 µl HB to have a total volume of 400 µl. We transferred it to a 2 ml Lo-Bind Eppendorf tube and added 1 volume (400 µl) of 50% iodixanol solution to give a final concentration of 25% iodixanol and mixed by pipetting. Following, we slowly layered 600 µl of 30% iodixanol on top of 600 µl prepared 40% iodixanol and put 800 µl of 25% iodixanol/sample mixture carefully on top. After using a swinging bucket centrifuge (spinning for 20 min at 3000 × g with the brake off), we collected 200 µl of nuclei from the 30%–40% interface into a chilled 1.5 ml tube. We then aliquoted 32,500 nuclei (to aim eventually for 10,000 nuclei) into a 1.5 ml tube. We again spun (10 min, 600 × g, 4 °C) and discarded the supernatant. We next resuspended the pellet in ~6 µl diluted nuclei buffer+RNase inhibitor (see user guide, Chromium Next GEM Single Cell Multiome ATAC+gene expression from 10x).

Library preparation
Single-nucleus libraries were generated using the 10x Chromium single-cell Instrument and Next GEM Single cell multiome ATAC+gene expression kit (10x Genomics) according to the manufacturer’s protocol (Chromium Next GEM single cell multiome ATAC+gene expression Reagent Kits User Guide; CG000338 Rev E). To minimize PCR biases during ATAC library construction, we performed an extra short quantitative PCR (qPCR) to assess how many cycles are needed to amplify our libraries without reaching saturation (1: 72 °C for 5 min, 2: 98 °C for 30 s, 3: 98 °C 10 s, 4: 63 °C for 30 s, 5: 72 °C for 1 min, 6: go to 3–20 times). At the end of the qPCR run, we determined the number of cycles that should be used by looking at the plot of fluorescence vs. the cycle numbers to find out what the maximum fluorescence was, and at what cycle number ¼ of that maximum fluorescence was reached and continued with the ATAC library construction by PCR with 7–9 cycles.

Multiome (snATAC-seq and snRNA-seq) sequencing
Final libraries were quantified using the Qubit dsDNA HS Assay Kit (Life Technologies). The fragment size of every library was analyzed using the Bioanalyzer high-sensitivity chip. All 10x scMultiome ATAC and gene expression libraries were sequenced on the NovaSeq 6000 instruments (Illumina) according to the manufacturer’s protocol (Chromium Next GEM Single Cell Multiome ATAC+gene expression reagent kits user guide; CG000338 Rev E). BCL files were converted into fastq using the command cellranger_arc mkfastq with the default parameter. The generated fastq files were processed using cellranger-arc v2.0.0 count function for read trimming, alignment, duplicate marking (ATAC), and UMI counting (GEX) (Cell Ranger ARC, 10x Genomics). The reads were aligned to the pre-built human reference genome GRCh38—v2020-A-2.0.0 (3 May 2021) provided by 10x Genomics.

Library preparation exclusion criteria and quality control
Nuclei were evaluated after staining with 4′,6-diamidino-2-phenylindole (DAPI) and visualized under a high-powered fluorescent microscope (minimum magnification: 40x). Only samples containing high-quality nuclei, defined as clump-free, debris-free, and exhibiting intact membranes, were selected for downstream experimental procedures to ensure data integrity.

After ATAC-library amplification, a quality control (QC) step was performed using qPCR to assess chromatin accessibility. Fold enrichment was calculated by comparing two accessible chromatin regions to two inaccessible chromatin regions. A minimum fold enrichment of 10-fold in accessible regions was required for sample inclusion. The following primer pairs were used:

• Accessible regions (POS).
 ∘ GAPDH: (F) 5′-GCCAATCTCAGTCCCTTCCC-3′, (R) 5’-TAGTAGCCGGGCCCTACTTT-3′
 ∘ KAT6B: (F) 5′-GAAGAGGCGGACCCAGCGGT-3′, (R) 5′-TTCCTGCCGGTCATCTCGCTT-3′

• Inaccessible regions (NEG).
 ∘ SLC22A3: (F) 5′-GGAGAGGGTGGACAGATTGA-3′, (R) 5′-TCAGCCTTGCTGCTACAGTG-3′
 ∘ QML93: (F) 5′-CACTGGTTGTCTTTGCAGGA-3′, (R) 5′-CCTGGGTCATATTGGGACAC-3′

Samples failing to meet these criteria—such as those with poor nuclei integrity or insufficient fold enrichment in accessible over inaccessible chromatin regions—were excluded from further analysis. Therefore following samples were excluded: ICGC_PCAL03_Tu5, ICGC_PCAL03_Tu6, ICGC_PCAL03_Tu8, ICGC_PCAL03_Tu9, ICGC_PCAL10_LN3, ICGC_PCAL10_Tu8, ICGC_PCAL10_Tu9, ICGC_PCAL25_LN1, ICGC_PCAL25_LN2, ICGC_PCAL25_Tu2, ICGC_PCAL25_Tu3, ICGC_PCAL25_Tu4, ICGC_PCAL25_Tu5, ICGC_PCAL25_Tu7, ICGC_PCAL25_Tu9, ICGC_PCAL25_Tu10, ICGC_PCAL34_Tu2, ICGC_PCAL34_Tu3, ICGC_PCAL34_Tu4, ICGC_PCAL34_Tu7, ICGC_PCAL34_Tu10, ICGC_PCAL37_Tu6. Samples with fewer than 1500 cells remaining after cellranger-arc v2.0.0 and QC (described below) were also removed: ICGC_PCAL03_LN2 (838 cells post QC), ICGC_PCAL10_Tu7 (52 cells post QC), ICGC_PCAL10_Tu10 (1339 cells post QC), ICGC_PCAL34_Tu1 (319 cells post QC), ICGC_PCAL37_LN1 (670 cells post QC), ICGC_PCAL37_Tu1 (317 cells post QC). In total, 38 samples (27 individual primary gland tumor biopsies and 11 matching lymph node biopsies) passed all QC criteria for inclusion in the analyses presented in this manuscript (Supplementary Dataset 1 batches and QC).

scMultiome data pre-processing
For each sample, output files from Cell Ranger ARC were processed with Seurat and Signac, with the following versions: R=4.1.0, Signac=1.4.0, Seurat=4.1.0. Cells were filtered if one of the following criteria was met: nFeature_ATAC>500, nFeature_RNA>500, percent.mt<20, TSS.enrichment>2. Cells were removed if both of the following conditions were met: nFeature_RNA>8000 and nFeature_ATAC>80000. Notably, we computed TSS enrichment using the Signac function (TSSEnrichment) and found that median TSS scores reported by ArchR on the same cells are on average, 2.5× higher than those reported by Signac, showcased by recomputing TSS scores using ArchR and Signac across 11 samples from ICGC_PCAL34 (Supplementary Dataset 1). 99% with TSS over 2 in Signac had TSS over 4 using ArchR, and the remaining 1% had a median of Signac TSS 3.6. Lastly, we did not opt for doublet correction and chose to filter on the top percentile of counts for RNA and ATAC, as doublet correction methods on RNA and ATAC individually did not produce results that agreed across samples74,75. Lastly, ATAC data from each sample were then reprocessed using a binned genome approach with 5 kb window sizes across hg38 using the Signac FeatureMatrix function. Downstream analyses were performed on the resulting Seurat object, containing an RNA and tiled genome matrix and were merged across patients.

Inter-patient analysis
Inter-patient data analysis was done using Python (version 3.9.13). Due to the large size of the patient-wide dataset containing all samples, we transitioned from using R to Python-based packages. Subsequent to QC steps and cell filtering described above, Seurat RNA and ATAC assays from each patient were converted to a Python-readable format using R packages MuDataSeurat and SeuratDisk. We then used the Muon package to load and merge each assay across all samples76 and proceeded with downstream processing.

For RNA, we removed mitochondrial and ribosomal transcripts prior to cluster/population calling. We had already used mitochondrial-based metrics to remove low-quality cells. However, mitochondrial and ribosomal reads can represent technical variability unrelated to the biological variability of interest77. Thus, for a more conservative approach in understanding inter-patient heterogeneity, we decided against the inclusion of mitochondrial and ribosomal reads across patients to minimize potential technical variation across patients influencing population-calling (clustering) results. We used scanpy default functions to normalize and scale data (normalize_total), retrieve and highly variable genes (highly_variable_genes), perform dimensionality reduction with PCA, and build the nearest-neighbor graph with defaults.

For ATAC, we normalized with TF-IDF and retrieved the top 5000 highly variable features, performed dimensionality reduction with LSI, and built the nearest neighbor graph on the top 50 LSI dimensions, excluding dimensions with higher than 0.7 Pearson correlation coefficient with depth (nCount_ATAC). We then used the Muon implementation of the weighted-nearest neighbor approach to build an integrated graph of the RNA and ATAC data, which captures the complementary information from both gene expression and chromatin accessibility in each cell28,76. We used Leiden clustering resolution with 1.3 on the obtained multi-modal graph. After clustering, we called peaks for each population using the MACS3 algorithm and built an inter-patient-wide feature set as described below (pseudobulk population calling). We then repeated dimensionality reduction steps for ATAC, recomputed WNN, clustering, and population assignments. Lastly, we did not proceed with batch correction, as nonmalignant cell populations showed inter-mixing across patients. This is indicative of the clustering algorithms preserving genuine biological differences. Batch correction could be potentially detrimental to removing meaningful differences across malignant epithelial populations78.

Intra-patient analyses
All intra-patient analyses were performed in R version=4.1.0, Signac=1.4.0, and Seurat=4.1.0. For prostate-specific analyses (Fig. 2), we merged samples coming from the primary tumor sections within a prostate gland of a given patient. For intra-patient analyses across tissues (Fig. 4), we merged all samples from metastatic lymph nodes and primary prostate tumor sections. For each merged object, we processed RNA, ATAC, and the subsequent Multiome object as described above, but used the equivalent R-based Seurat/Signac implementations. For RNA processing, we used SCTransform to normalize the filtered count matrix and PCA to run louvain graph-based clustering on the first 50 principal components, and for ATAC, we used RunTFIDF to normalize the data, found highly variable regions using FindTopFeatures, and RunSVD to project ATAC data to lower dimensions. We built the nearest neighbor graph on the top 50 PCA and top 50 LSI dimensions, excluding dimensions with higher than 0.7 correlation with ATAC depth (ATAC_nCount). The weighted nearest neighbor graph was built using FindMultiModalNeighbors. Populations were called on this graph using FindClusters for each patient with multiple resolutions, we tested and annotated clusters across resolutions 0.4, 0.6, and 0.8; and found 0.4 to cover the major cell types across all patients.

Pseudobulk peak calling
We wrote a Snakemake pipeline for calling and merging peaks for any group of cells in our dataset. First, we generated a fragment TSV file for every cell in a given group (in this case, our called populations) per sample. Second, MACS3 was then used to call peaks on the fragment files with the following parameters: -f BED -q 0.01 -g 2.7e9 --keep-dup all --nomodel --nolambda --shift −75 --extsize 150 --bdg --call-summits79 and peaks were iteratively merged80. We then re-computed the FeatureMatrix in Signac and either transferred to Python for inter-patient analysis or added them to the Seurat objects for intra-glandular/intra-patient analyses. Third, the new FeatureMatrices were used to obtain signal over promoter and non-promoter elements and compute ChromVar-based transcription factor enrichment scores29. For each peak catalog, we assigned whether a peak was a promoter element if it was within +−1 kb TSS in gene code transcript annotations, otherwise, the peak was annotated as non-promoter, which was used in subsequent Pearson correlation analyses described below. We used a merged peakset across all samples for inter-patient and patient-specific peaksets for intra-patient analyses.

Population annotations
We used a multi-faceted manual and automated approach to go about population annotations, independently utilizing RNA and ATAC data. For every population, we computed the marker genes and chromVar transcription factor motif enrichments by looking at the top upregulated genes relative to the entire sample of interest (Seurat’s FindAllMarkers). For every cell in our dataset, we computed the singleR assignments using the publicly available default datasets and a prostate-specific dataset81,82. We also computed each cell’s ModuleScore across gene expression markers reported in relevant papers24,82–84 using Seurat’s AddModuleScore, with additional cell-cycle markers obtained from https://github.com/hbc/tinyatlas. We then visualized all singleR’s cell type assignment proportions per population and visualized ModuleScores to facilitate manual cell type assignments. For every population, we checked the alignment of the highly expressed genes with expected genes per cell type with literature search and available atlases (protein atlas.org). We also utilized ATAC-based transcription factor-deviation scores per population, as ATAC signal tends to be less prone to biological variability such as cell cycle, and TF-motif enrichment computations rely on more global genome-wide statistics. We cross-referenced TFs based on their enrichment relative to the expected average (which would contain more malignant epithelial cell TFs) and annotated populations according to TFs. Furthermore, more than one expert was exposed to the population markers and visualizations and independently annotated populations according to the steps above to ensure agreement. This rigorous, multi-layered approach provided confidence that our findings represent biological cell types within and across patient samples.

Assessing the degree of phenotypic heterogeneity in locus-specific populations
We assessed heterogeneity between populations in two ways. First, we reported how inter-mixed populations were based on variables of interest, such as patient, tumor locus, and tissue, using Shannon entropy. Additionally, to compare scores across patients when reporting the median across patients, we normalized the scores to be between 0 and 1 when 0 indicating no entropy and 1 indicating maximum entropy relative that can be achieved given the number of loci and cells per loci. Second, we derived Pearson correlation values on DESeq2 normalized pseudobulked signal across the three feature sets of non-promoter DNA elements, as well as promoter elements and gene expression. We obtained normalized pseudobulked values per population by summing raw counts of each population according to the three feature sets and normalizing them using DESeq2 vst85. To assess how transcription factor families drive variance, we conducted PCA on median chromVar motif deviations per population. We then binned the feature loadings into 30 bins across 633 motifs, covering the top and bottom ~5th percentiles of features contributing to a given axis of variation.

To contextualize biological programmes driving intra-prostate tumor heterogeneity, we used the SCENIC+ pipeline on tumors from ICGC_PCAL34 and ICGC_PCAL0337. Pre-processing of ATAC data for the SCENIC+ pipeline was done using the polars_1xx branch of the pycisTopic package for scalability purposes. ScRNA preprocessing and the SCENIC+ pipeline were run according to recommendations highlighted at https://scenicplus.readthedocs.io/en/latest/. eRegulon specificity scores for each malignant epithelial cluster were calculated using sceniecplus.RSS.regulon_specificity_scores() to select the top 10 most specific and confident eRegulons per malignant epithelial cluster. The genes of the highly specific eRegulons were used to create gene lists for pathway enrichment using the gseapy Python package.

CNV inference
We used the NUMBAT algorithm39,86 to perform CNV inference. We performed NUMBAT per sample to ensure tumor vs. normal calls were properly cross-referenced with our RNA/ATAC-based cell type annotations. We then used NUMBAT on a patient-wide basis to get a sense of all CNV segments with higher accuracy. However, due to the large number of cells from each patient, we kept only the cells annotated as malignant epithelial and had to further subsample the number of cells used in the per-patient analysis as follows, for the pipeline to run. To allow computational feasibility, we kept a maximum of 3000 cells per sample in ICGC_PCAL03, ICGC_PCAL10, ICGC_PCAL37, and 6000 malignant epithelial cells per sample in ICGC_PCAL34 and ICGC_PCAL25. Lymphocyte populations from ICGC_PCAL10 and ICGC_PCAL37 were used as background, since they had the largest number of cells from the TME. For visualization purposes (Fig. 2b), we used the most common CNV call in all cells assigned to a given clone. To compute the difference between clones, we used the Hamming metric, which measures the distance between two vectors of equal length. The CNV calls per cell were then used to compute this metric between all cells and clones, to broadly assess the difference between clones in a given patient.

Assessing the degree of clonal architecture heterogeneity in populations
We assessed how similar the clonal makeup of each population was according to the Jensen–Shannon (JS) divergence metric. We used JS divergence as it provides a symmetrical and finite measurement of the difference between two probability distributions. For every malignant epithelial population with more than 200 cells, we computed the proportion of cells that were assigned to each clone, and the JS distance was calculated per population pairs.

Differentially accessible regions and transcripts identification and analyses
To identify differentially accessible peaks and transcripts, we pseudobulked raw counts from malignant cells from a given patient per sample across all features. Only cells annotated as malignant epithelial were kept. We then used DESeq2’s ~patient+tissue to retrieve differential peaks and transcripts of primary prostate and metastatic lymph nodes, FDR thresholded at <0.1. We then computed transcription factor DNA motif enrichment scores using HOMER on differentially accessible peaks with absolute LFC > 1 and q < 0.187.

Ambient RNA correction
To gain confidence in our differential analyses, we repeated our pipeline but used ambient RNA corrected count matrices using CellBender43. We optimized parameters per sample, except for the FPR, which we had to keep consistent across all, according to the authors (FPR=0.05). Overall, CellBender retained a median of 40% more cells than Cell Ranger and did not have an apparent benefit to the malignant epithelial cell expression patterns (Supplementary Fig. 8a, b). We kept only the cells that passed QC according to our scMultiome data pre-processing as outlined, and repeated our differential analysis as described above.

Validation of gene signatures in an external single-cell RNA seq dataset
To validate the AUCell enrichment scores of gene signatures, we utilized an external dataset47. Specifically, we downloaded the Seurat eRegulons object (primecut.integrated.geneexp.seurat.rds) from the publicly available repository hosted at https://data.mendeley.com/datasets/5nnw8xrh5m/1. This dataset contains gene expression profiles of castration-sensitive prostate cancer (mCSP) cells from various metastatic sites. We subset the dataset to exclude TME cell types, isolating a total of 3976 malignant epithelial mCSPC cells. These cells were distributed across four metastatic sites: lymph nodes (n=1096), bone (n=1263), lung (n=641), and liver (n=976). We calculated AUCell enrichment scores for our lymph node up- and down-regulated genes across the treatment stages in all metastatic sites to assess the dynamics of our gene signatures in response to treatment. We calculated AUCell enrichment scores for lymph node up- and down-regulated genes across treatment stages in all metastatic sites to assess gene signature dynamics. Scores were normalized between 0 and 1, and a Kruskal–Wallis test determined whether at least one treatment group differed significantly. Pairwise comparisons were conducted using the Wilcoxon rank-sum test. This approach allowed us to investigate both site-specific and systemic changes in gene signature dynamics during mCSPC progression and treatment.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have addressed my comments effectively.", "reviewer_2": "Authors have partially addressed my concerns related to the manuscript. Still major issues remain.\nFor CNV analysis authors could consider more optimized (e.g. python) implementation of the algorithms. None of the less, authors do demonstrate data to alleviate concern regarding subsampling. However, how the choice of reference cell type impacts the CNV calls has not been properly addressed. As different cell types expressed different set of genes this can have fundamental impact on inferred CNVs and thus on conclusions.\nWhile authors do present innovative comparisons through TF motif enrichment, I do not agree that this would remove the need to address biology behind sites through standard approaches. While motif analysis has novelty, results need to be presented and put into context of differential expression (and possibly accessibility) analysis and pathway activation, that I would consider as baseline benchmark for any more complex analysis.\nAuthors have not properly addressed their claims of plasticity. In my view, raising plasticity as a key concept, will require authors to properly address the transcriptomic (and/or epigenetic) states (and provide definition of cell state) in relation to CNV based clones. Current revision, where non-formal definition of plasticity is given at the end of discussion is far from sufficient.", "reviewer_3": "The authors clarified a sufficient amounts of critical points and I am convinced of the robust and rigorous analysis of the data.\nMaybe, given the lack of further validation that prostate cancer cells converge to a phenotype similar to an inflammatory-like status in the LNs, the title of the manuscript could be slightly modified.", "reviewer_4": "The authors have answered my questions."}, "author_response": "The authors have approached the question of tumour heterogeneity in a novel way, resulting in interesting findings regarding phenotypic heterogeneity in malignant epithelial cells and its absence in the microenvironment cells. The methodology provides an interesting perspective to answer biological questions about prostate cancer progression. It appears technically robust and well interpreted. The finding that malignant cells acquire an inflammatory cell-like state in lymph node metastases is interesting (TME-mimicry) – if confirmed in larger cohorts, this could be a therapeutic target. The study is limited however by the small number of patients (which the authors acknowledge), making the findings less generalisable.\nI have the following comments for discussion:\nC1.1. Shannon Entropy analysis to measure population intermixing is interesting; the finding that fibroblasts are most inter-mixed raises the question of whether cancer associated fibroblasts are all the same regardless of the genetic alterations in the malignant cells..?\nAlthough our primary analyses focused on malignant epithelial cell heterogeneity, we were able to detect fibroblast populations in the dataset. However, these fibroblasts represented fewer than 1% of total cells in four out of five patients (ICGC_PCAL03=2.6%, ICGC_PCAL10=0.01%, ICGC_PCAL25=0.02%, ICGC_PCAL34=0.01%, ICGC_PCAL37=0.01%) (Rebuttal Figure 1.1A). The overall low abundance of fibroblasts (median = 61 cells per sampled region) limited our ability to perform robust spatial comparisons.\nHowever, in one patient (ICGC_PCAL10), two regions Tu3 and Tu6 contained a higher number of fibroblasts (141 and 278 cells, respectively). Comparing ambient-RNA corrected fibroblast cells of these two sites suggested that Tu3 and Tu6 fibroblasts likely correspond to iCAFs and myoCAFs, respectively. This classification was based on the enriched hallmark pathways (see citation 1 supplementary figure S3) 1, such as IL6-JAK STAT3 , IL2 STAT5, complement, TNFalpha, interferon gamma and inflammatory responses in iCAF, as well as epithelial-mesenchymal transition (EMT) in myoCAFs (Rebuttal Figure 1.1B).\nAlthough we are underpowered to comprehensively investigate how the TME maintains a pro-tumorigenic environment, these preliminary observations warrant future studies specifically designed to comprehensively profile the TME.\nRebuttal Figure 1.1 A) Proportions of cell type annotations across all cells isolated from each patient B) Hallmark gene set enrichment analysis between fibroblast cells in ICGC_PCAL34 Tu3 and Tu6 regions corresponding to iCAFs and myoCAFs respectively (Figure S3 from 1)\nC1.2. Fig S4 suggests that cells with the same genotype can have multiple different phenotypes (gene expression clusters) – but if the microenvironment is not causing this difference (because the microenvironment is largely the same across patients and genomic clones, except B cells), what explains the phenotypic heterogeneity within genomic clones? Could it be a limitation in the resolution of CNAs in identifying genomic clones? Might it be that they are different genomic clones at the SNV level? Could it be epigenetic ‘subclones’?\nNon-genetic yet inheritable drivers of phenotypic heterogeneity are well established in developmental biology and gaining recognition in cancer 2,3. Specifically, differentiating cells arise from the acquisition of epigenomic variation, such as changes in chromatin accessibility, over narrow (~500bp) genomic loci 4. Importantly, epigenomic variation can replicate and segregate in a Mendelian manner with DNA during cell division 5. Furthermore, transient cell state switches can be observed from both epigenomic and/or gene expression variation in response to temporal environmental stimulation 6.\nIn this manuscript, phenotypic heterogeneity is evident in both chromatin accessibility and transcriptional profiles. Specifically, the cluster-calling algorithm integrates both mRNA expression and chromatin accessibility data to construct a graph of cellular relationships 7, leveraging the complementary information from these modalities to define cell populations.\nThe association between microenvironmental changes, such as colonization of a new organ, with convergent phenotypic plasticity observed across malignant cells sampled from multiple tumor regions and patients may reflect pressures exerted by microenvironmental influences. However, this remains difficult to directly measure, as our experimental design and current technologies are underpowered to fully map the diversity and complexity of interactions between malignant cells and their niche.\nAs we mention in our response to R1.1, although we are underpowered to comprehensively investigate how the TME maintains a pro-tumorigenic environment, future studies specifically designed to comprehensively profile the TME are warranted.\nThe advent of spatial epigenomic and transcriptomic sequencing technologies is of particular interest to deciphering how the microenvironment may direct malignant cell plasticity, which we now mention in the discussion section of our revised manuscript on page 11 line 6-20:\n“Phenotypic heterogeneity in malignant cells, both between and within patients, suggests that primary tumor formation is shaped by patient-specific tumor evolutionary processes. In contrast, TME cell types, such as immune and stromal cells, are relatively conserved across patients. The similarity in TME composition aligns with prior single-cell studies for other malignancies 49–51 and supports targeting the TME as an alternative therapeutic strategy in high-risk prostate cancer 52. While our Shannon entropy analysis revealed relative conservation of TME composition across patients and spatial sites, it does not preclude the possibility that local, spatially restricted interactions contribute to phenotypic heterogeneity between malignant populations 53–55. Whether TME and malignant epithelial populations with similar phenotypes congregate in distinct spatial pockets 54,55 and if local cellular interactions result in genetically distinct cells to adopt shared phenotypes remain to be explored. Addressing these will require a different setup, such as high-resolution spatial sequencing technologies, to more precisely map the interplay between the TME and malignant populations.”\nC1.3. A limitation with the methodology is that the spatial relationships between malignant cells and their associated TME cells are lost using the single cell approach – hence, it is not possible to know, for example, whether specific genomic clones within a locus are associated with specific TME subtypes. Perhaps the authors could acknowledge this in the discussion?\nSee response to comment C1.2 – We have acknowledged this in our manuscript discussion section lines page 11 line 6-20,\nMinor\nC1.4. In Fig 3A, I think it’s easier for the reader if the CNV segments are ordered by chromosome (1 to 22) – the current sorting order seems unusual\nWe have updated Figure 3 by ordering CNV segments by chromosomes.\nC1.5. “malignant populations exhibited more heterogeneity across patients, whereas the TME cell types are more similar across patients and tissues.” – The authors clearly mean phenotypic heterogeneity, but could help to state this more explicitly\nWe have revised the sentence in the discussion section (See C1.2).\nC1.6. Repetitive text – for example, lines 367 and 378 basically describe the same thing, don’t they?\nWe have reformatted our discussion section to minimize redundancies, see C1.2\nReviewer #1 (Remarks on code availability):\nFIX: The github repo link is not accessible (perhaps still private?)\nThe Codeocean link is accessible, and code appears to be available to reproduce all the figures in the manuscript. Software and environment requirements seem to be well documented.\nWe have updated our GitHub link.\nReviewer #2, expertise in prostate cancer tumour heterogeneity, lineage plasticity and single cell sequencing (Remarks to the Author):\nThe authors present a single cell multiome dataset from 38 unique samples collected through multiple primary tumour and pelvic lymph node metastasis biopsies of 5 high-risk prostate cancer patients. Altogether, the study comprises the transcriptomic and chromatin accessibility profiles of 282,956 unique cells, boasting an impressive amount of high-resolution data accompanied by a commendable amount of details regarding each sample and patient, including biopsy site, copy number variation profiles, driver mutations, clinical baseline, follow-up data, and histological features. Given all this information, the authors are well-positioned to assess genomic, epigenomic, and transcriptomic leading to metastatic dissemination. However the manuscript has some major shortcomings in this regard.\nMajor:\nC2.1. The main issue in the manuscript regards the copy number variation analysis. Accurate assessment of CNV profiles is a prerequisite for finding matching cancer clones in the primary tumor and the metastatic sites. For this task the authors have used Numbat, a tool that combines count data with and haplotype information to infer the copy number profiles. However the described workflow contains some ambiguity and plain ad hoc solutions. On rows 727-735 the authors state the following: \"...due to the large number of cells from each patient, we kept only the cells annotated as malignant epithelial and had to further subsample the number of cells used in the per-patient analysis... To allow computational feasibility, we kept a maximum of 3000 cells per sample... We found ICGC_PCAL10’s lymphocyte to work best as background in all patients except for ICGC_PCAL37, where its own lymphocyte populations worked best.\" It is advised that the authors completely revise this analysis with sufficient computational resources and a more meticulous computational pipeline. Furthermore, it would be useful that the authors present a complete \"raw\" inferred copy number changes for all epithelial cell populations and autosomes at some point in the manuscript, e.g. the supplementary (as seen in e.g. Fig 4a of the Numbat publication https://www.nature.com/articles/s41587-022-01468-y).\nWe acknowledge the importance of accurate CNV inference and thank them for pointing out ambiguities in the methods section. Numbat automatically chooses the background based on the correlation of supplied pseudobulked background cells with individual malignant epithelial cells. We have revised our sentence in the manuscript for clarification (Methods: CNV Inference page 22, line 2). Our rationale for subsampling was motivated both by computational feasibility and statistical considerations.\nSubsampling was necessary due to R’s memory ceiling, as Numbat computes a pairwise distance matrix (ape::nj) which grows quadratically with the number of supplied cells, and we found the limit to be around ~35,000 cells. We estimated the maximum number of cells which enables computational feasibility of per-patient CNV inference to be 3,000, with each run taking ~1TB of RAM and ~7-10 days to complete depending on the patient. Thus, higher input sizes were not feasible for per-patient runs. Furthermore, we reasoned that for each patient’s CNV inference, would mitigate the risk of overfitting clonal inference to any one spatial location by subsampling the number of cells obtained from each region to roughly an equal number. Finally, in response to the reviewer’s request and technical considerations, we showcase results for ICGC_PCAL34 when subsampling to a maximum of 3,000 versus 6,000 cells per spatial region. We found little difference in segment detection between the two runs, with 93% of basepairs detected as segments in the 6,000 cells per region run are detected as segments in the 3,000 cells per region run (Rebuttal Figure 2.1A). We then computed the weighted hamming distance of CNV labels between cells, which represents the difference between CNV segment calls between 6,000 versus 3,000 across the genome weighted by the length of the CNV ranging from 0 to 1. 0 implies no difference in CNV labels and 1 implies all CNV labels differing in the 3,000 cell from 6,000 cell run. Each CNV segment is weighted by its length in basepairs divided by total number of basepairs covered by all segments called in the larger 6,000 cells/region run. We saw a median hamming distance of 0.1 across shared (N=15,000) cells, indicating little difference across the called CNV segments (Rebuttal Figure 2.1B). We have now included segment calls on the single cell level in Supplementary Figure 5B.\nRebuttal Figure 2.1: CNV calls across two computationally feasible experimental set ups A) Segments detected by Numbat across the two runs. Top and bottom bars represent 6,000 and 3,000 cells per region (total N = 30,000 and 15,000 respectively) B) Weighted hamming distance between shared cells (N=15,000)\nAssociated manuscript revisions:\nSupplementary Figure 5.\nClonal architecture across sampled loci\nA) Pie charts depicting distribution of separately inferred genotypic clones for cells coming from each primary tumor loci and each matched lymph node metastases per patient. B) Segment CNV calls across the human autosome per single cells according to their maximum likelihood estimation (MLE, left) and maximum a posteriori (MAP, right). Annotation bars represent sampled regions and inferred clones.\nC2.2. The analyses regarding tumor heterogeneity between different sampling locations is limited to descriptive metrics such as Shannon's entropy and Pearson correlation. A more fruitful analysis would be to identify the relevant biology varying between these sampling locations through a standard differential expression and gene set enrichment analysis. This is also related to the clonality analysis as its results would guide the relevant comparisons.\nWe would like to clarify that, beyond descriptive metrics, we performed a comprehensive analysis to identify biologically relevant variation across different sampling locations. We reasoned that heterogeneity across sites can be captured by evaluating the variance in transcription factor (TF) motif enrichments across accessible non-promoter elements, which have been shown to reflect increased cancer-type specificity compared to accessibility at gene promoters or gene expression levels 15.\nFurthermore, given that genome-wide TF motifs were not utilized to compute the cell graphs used for intra-patient population calling (as opposed to ATAC and RNA information), we reasoned that TF motifs would be a more robust and orthogonal metric for assessing population variability and avoids the double-dipping problem (when the same features used for clustering are used for differential analysis). To this end, we performed principal component analysis on TF motif enrichments across multiple samples, revealing motifs for FOS and JUN-related TFs as top contributors to the major axes of variation across malignant populations (Figure 2D, Supplementary Figure 3). We now provide differential motifs across all clusters (Rebuttal Figure 2.2), although we found summarizing this information across patients and groups to be more biologically meaningful in deciphering sources of heterogeneity.\nAdditionally, we performed differential chromatin accessibility analysis between malignant epithelial cells from primary prostate (Tu) and lymph node (LN) sites, followed by Hallmark pathway enrichment and TF motif enrichment analysis of differentially accessible regions to identify biologically relevant differences across sampling sites (Figure 4E-F). Together with Pearson correlation analyses of normalized pseudo-bulk chromatin accessibility across sites (Figure 4D), these analyses align with the reviewer’s suggestion, providing both a standard and chromatin accessibility-informed approach to identifying relevant biology underlying heterogeneity across different tumor locations.\nRebuttal Figure 2.2: Marker transcription factors visualized per clusters called in primary tumor foci of each patient.\nC2.3. The authors reference phenotypic plasticity throughout the manuscript. Cell plasticity represents the ability of cells to be reprogrammed and to change their fate and identity. Any formal analysis regarding plasticity is lacking from the manuscript, as there is no temporal data that would indicate a change in phenotype. To substantiate the claims regarding plasticity with these data, the authors would need to 1) robustly define the subclonality of cancer cells (e.g. pseudotime) and 2) show divergent gene expression and chromatin accessibility profiles in these subclones. In it's current form, the differential gene expression and chromatin accessibility analyses are limited to simply comparing the lymph node metastasis and primary tumor populations. The prospect of being able to place subclone populations on a timeline ranging from the primary tumor to the metastasis would add a lot of value to the manuscript and the field overall.\nWe thank the reviewer for their comment and share their interest in understanding potential ancestral relationships of malignant epithelial populations. We refer to plasticity as it points toward phenotypic changes between cells that cannot be explained by genetics alone.\nOther definitions of plasticity have been suggested to distinguish microenvironmental and cell-intrinsic sources of heterogeneity according to phenotypic read-outs 3, and we have now contextualized our findings accordingly in the discussion section (page 12, lines 26-38).\nOur intention in generating this dataset was to clarify how single-cell-based populations could be utilized for identifying common ancestors that may exhibit migration capabilities. We considered methods for population tracking and found the CNV-based approach applied in our manuscript to be the most technically-sound.\nHowever, inferred CNVs did not support a genetic (i.e.: clonal) history between regions, which led us to focus on expanding our study to the phenotypic comparisons between lymph nodes and primary tumors.\nWe agree with the reviewer on the value of placing cells on a timeline, and explain our rationale for excluding pseudotime or velocity-based approaches from the main figures below:\n1. Algorithmic ancestral relationships identification (pseudotime and velocity-based approaches): Population dynamics of single cells can be modelled using velocity-based algorithms, which rely on the observed shifts in intronic and exonic read ratios. While we implemented this approach in our dataset, we found that it was not an ideal method for solid tumor samples. These methods assume relatively uniform transcriptional dynamics and consistent splicing behavior across cells and would be correctly utilized for single-sample analysis. We must note that velocity-based approaches, or other graph-based algorithms, can be sound and useful for more well-defined systems of development combined with more reductive/controlled models with well-known directional patterns in cell fate, where the entire space of intermediate cell states are captured within a single experiment 16,17. However, the assumptions made by these algorithms are often not met in tumor tissues and are exacerbated by pre-processing of multiple samples, which leads to distortions in algorithmic inference of trajectories 18. This is especially prominent in velocity algorithms, which are unsuitable for our complex multi-sampled cohort.\nLastly, in cancer settings, single-cell sequencing provides a snapshot of phenotypic behavior during the time of surgery. Algorithmic inference of trajectories using graph-based models requires a knowledge of directionality and sufficient similarity between populations to avoid false positives. This is complex in the case of cancer. For example, in normal development, well-researched stemness markers can be used for guidance in assessing directionality and avoiding false-positives. Whereas in cancer settings they may represent progenitor populations 20,21 or growing daughter cells as a result of cancer progression 22. Thus, in the case of our cohort, we found identifying ancestral relationships and population dynamics using the available phenotype-based algorithms to be unsuitable for defining the relationship between malignant populations.\nWhile acknowledging these limitations, we performed pseudotime analysis on one patient ICGC_PCAL34, with four tumour regions and a single lymph node, presents the simplest case in our cohort for defining trajectories using Monocle3 19. We used Monocle3 to preprocess the data subsetted on all cells belonging to malignant epithelial clusters (Rebuttal Figure 2.3 A). We observed partitioning to resemble inter-loci heterogeneity, where the best possible path between cell states are localized to single loci, with the exception of Tu5 and Tu8. While relaxing partitioning by Monocle3 can reveal continuous paths between sampled regions (Rebuttal Figure 2.3 B-C), it can lead to spurious paths that may not be representative of the biological differentiation path and is not recommended practice. Thus, pseudotime analysis is constrained by patient-specific inter-loci phenotypic states rather than inter-loci transitions. Additionally, the Monocle-based graph and called populations were similar to those already captured by the weighted nearest neighbor graph-derived populations, underscoring the robustness of our integrated clustering approach (Rebuttal Figure 2.3 D-E).\n2. Copy number variation(CNV)-based approaches: Our intention for incorporating CNVs was to (i) more confidently annotate populations and (ii) utilize CNVs as surrogates for lineage tracing in our single-cell sequencing analysis as CNVs have been used to track evolution in cancer progression 23 and address the complexities of tumor heterogeneity, which render the methods above unsuitable. While we initially utilized the CNV-based approach to identify and isolate malignant populations in the prostate with migration capabilities to the lymph nodes, our findings suggested a more homogeneous clonal architecture among the malignant epithelial populations, an observation which we elaborated in our paper. Thus, we could not confidently place subclones on a genetic timeline representing the seeding to the lymph nodes. This underscored that there may be phenotypic changes that occur without detectable clonal divergence in CNV-based markers, so we turned to primary and metastatic comparisons.\nRebuttal Figure 2.3: Pseudotime reveals constrained cell state trajectories to specific loci. Pseudotime trajectories overlaid on Monocle3- derived UMAPs on RNA expression data from ICGC_PCAL34 malignant epithelial cells labelled by A) Sampled region, B) Pseudotime trajectory on partitions, root cells identified algorithmically C) Pseudotime trajectory irrespective of partitioning, root cells identified algorithmically D) Monocle3 identified populations E) Weighted-nearest neighbor (RNA+ATAC) populations.\nMinor\nC2.4. Fig. 1a states there are \"5x treatment-naïve patients\" while methods describe that \"one patient (ICGC_PCAL25) who received neoadjuvant hormone therapy with LH-RH agonist 6 weeks before surgery\".\nWe have fixed figure 1 to accurately reflect the clinical parameters of this patient cohort:\n5x patients with locally advanced high-risk prostate cancer\nC2.5. Fig. 1d dotplot is confusing as the markers are not very cell type specific (e.g. HOXB13 and FOXA1 do not appear to score highly in any cell type).\nThe dot plot displays the deviation score relative to an average chromatin accessibility signal computed across DNA motifs. Since this reference profile is primarily shaped by AR-dependent malignant populations, we generally expect malignant epithelial cells to exhibit minimal deviation. To clarify this, we have revised the figure legend accordingly. Furthermore, we showcase this phenomena by subsampling up to 1,000 cells per cell type in one patient, where malignant epithelial cells show more pronounced increased deviations of FOXA1/HOXB13, which would otherwise be obscured by dominant AR-driven signals in the non-subsampled average (Rebuttal Figure 2.5)\nRebuttal Figure 2.5: Motif deviation scores re-computed in subsampled cell types in ICGC_PCAL37\nC2.6. Fig 2d. TF family annotation does not add anything to the plot. Column annotations are also missing although they presumably describe patient of origin.\nWe have updated Figure 2D to include the column annotations for patient identifiers. We retained the TF family annotation on the y-axis to make TF groupings apparent to the reader, which may imply their shared structural and functional characteristics. For example, clusters of FOS/JUN, FOX or HOX family members showed consistent activity in specific patients which could suggest lineage-related programs or motif-similarity of the reported factors.\nC2.7. Remarks about the tumor microenvironment components being similar across different tumors is to be expected and little value to the manuscript (lines 125-140, 367-379).\nWe report the inter-mixing of tumor microenvironment components across patients and organs, as it serves as a reference point for the expectations for the inter-mixing of malignant epithelial cells. The inter-mixing of tumor microenvironment components is additionally used to assess the need for batch effect correction. For example, we observed basal-epithelial cells, which are biologically constrained to the prostate and not the lymph nodes, to contain higher inter-mixing across patients, but not across organ sites in our data.\nReviewer #2 (Remarks on code availability):\nCodeocean link asks to log in with email and github link gives 404 error\nReviewer #3, expertise in single cell multi-omics and prostate cancer (Remarks to the Author):\nPerforming elegant analyses of single nuclei multiome data from matched clinical speci in mens from the prostate and from lymph node dissections, this work by Keshavarzian et all., provides some evidence of cell plasticity and interindividual ethereogenicty of such plasticity during prostate cancer evolution and metastatic seeding from the prostate to the lymphnodes. According to the performed analysis some of the processes identified, specifically linked to inflammation, could be exquisitely associated with prostate cancer dissemination.\nThe analysis is conducted rigorously and appropriate statistical methods are used. The manuscript is clearly written and the content of general relevance. Some questions remains, that can be addressed with additional analyses of the dataset produced, which is also a great resource for the community and should be easily accessible.\nMajor comments\nC3.1. The phenotypic plasticity is only inferred by the uniformity of the clonality analysis across the specimens from the same individuals.\nLines 277-9 “Thus, the patient-level clonal architecture stability suggests a predominant role for phenotypic plasticity as malignant prostate cancer cells evade the primary site to colonize the lymph node.” To confirm phenotypic plasticity, additional anaysis should be performed using independently the ATAC and RNASeq data and then possibly the data combined. Phenotypic plasticity can be referred as cellular trandifferentiative states dictated by both chromatin and transcriptional plasticity. This can be assessed by presenting the analysis for the 4 patients separately and computing the specimens data per patient.\nWe appreciate the reviewer’s insightful comment. Our initial inference of phenotypic plasticity was based on the observation that malignant cells in lymph node metastases maintain clonal architectures resembling primary tumors, yet exhibit distinct gene expression and chromatin accessibility states when we analyze differential phenotypes on the cohort level. Thus, allowing us to capture the non-genetic phenotypic changes that cells undergo as a result of colonizing a new organ. Specifically, this is showcased with (i) differential gene expression and chromatin accessibility between malignant cells of the primary and lymph node sites, with patients as covariates (Table 3), (ii) lymph-node upregulated genes are associated with chromatin variants with immune/inflammatory transcription factor DNA binding motifs (Figure 4E-F). This analysis is done on gene expression and chromatin accessibility independently.\nAdditionally, given that gene expression read-outs can capture transient phenotypic heterogeneity such as transcriptional bursting and cell cycle, and chromatin accessibility read-outs can include more heritable and stable sources of heterogeneity, we used the integrated weighted-nearest neighbor approach which considers data from both RNA and ATAC to build the underlying cell-graph used for clustering to effectively capture the variability between cell states. When focusing on a single modality at a time and then combined (RNA, ATAC, and WNN), the overall structure of the graphs are conserved (Rebuttal Figure 3.1). We have revised the manuscript to clarify how we conducted the analysis and to further justify the contribution of different sources of non-genetic variation, namely epigenomic and transcriptional variation as described in our response to comment C1.2.\nRebuttal Figure 3.1: UMAP plots generated based on weighted-nearest neighbor considering information on ICGC_PCAL34 from both ATAC & RNA (top), ATAC-only (middle), and RNA only (bottom), coloured based on sampled region (left column) and cell type annotation (right column)\nC3.2. Can the authors show directionality of differentiation from the primary tumor to the different sites of colonization? Can a pseudotime RNA velocity analysis be performed?\nPlease refer to C2.2\nC3.3. The analysis of the clonality is still very interesting: is there a relationship between the lymphnods anatomical regions and the clonality? Can the authors comment on whether the high clonality is derived from high er proximity of the sites investigated? How far apart are the sites in primary tumors and LNs? And how many cells per loci were on average retrieved? Are the co-annotated cells types providing some evidence that those would also influence the colonization of the Lymphnode?\nDue to the nature of pelvic lymphadenectomy, lymph nodes are removed en bloc, which prevents precise reconstruction of their anatomical locations or spatial proximity.\nPathology images are available for primary tumor sites and schematic representation of the inferred spatial distribution of tumor loci within the prostate gland of each patient are shown in Figure 2A. Detailed pathology reports have been added to the rebuttal for reference (Rebuttal Figure 3.3). Regarding spatial relationships within the primary tumor, we sampled anatomically distinct foci (see Methods), but distances between primary and lymph node sites were not quantifiable in our study. Thus, we are not able to derive a relationship between phenotypic and genotypic diversity and physical proximity. However, we observed that, with the exception of patient ICGC_PCAL34, lymph node samples exhibited similar clonal architecture both to one another and to their matched primary tumors. In contrast, some intraglandular clones, such as clone 3 in ICGC_PCAL10, were restricted to specific regions of the primary tumor and not detected in any lymph nodes, suggesting their limited dissemination (Figure 4B, Supplementary Figure 5). These findings indicate that major clones can disseminate across lymph nodes, while some subclones remain confined to specific primary tumor regions. We cannot exclude the possibility of spatially proximal regions of the prostate tumors disseminating to the lymph node which may share a similar clonal architecture to our sampled regions.\nThe number of cells retrieved per locus is summarized in Table 1 (\"batches_qc\"), which also provides an overview of sampling depth across primary and metastatic sites (Methods: Sample and library preparation for 10x single-cell multiome ATAC and gene expression).\nRegarding whether the co-annotated TME cell types may influence colonization of the Lymphnode: while the majority of cells in our dataset are malignant epithelial cells (see response to C1.1), we do detect a range of TME cell types. However, as noted in our response to comment C1.1, the limited sample size and proportion of non-malignant cells precludes a robust investigation of specific TME-malignant crosstalk axes (e.g., fibroblast-epithelial interactions). We have further clarified this point in the discussion (Page 11, lines 6-20): “Phenotypic heterogeneity in malignant cells, both between and within patients, suggests that primary tumor formation is shaped by patient-specific tumor evolutionary processes. In contrast, TME cell types, such as immune and stromal cells, are relatively conserved across patients. The similarity in TME composition aligns with prior single-cell studies for other malignancies 49–51 and supports targeting the TME as an alternative therapeutic strategy in high-risk prostate cancer 52. While our Shannon entropy analysis revealed relative conservation of TME composition across patients and spatial sites, it does not preclude the possibility that local, spatially restricted interactions contribute to phenotypic heterogeneity between malignant populations 53–55. Whether TME and malignant epithelial populations with similar phenotypes congregate in distinct spatial pockets 54,55 and if local cellular interactions result in genetically distinct cells to adopt shared phenotypes remain to be explored. Addressing these will require a different setup, such as high-resolution spatial sequencing technologies, to more precisely map the interplay between the TME and malignant populations.”\nRebuttal Figure 3.3 (New Supplementary Figure 2): Detailed pathology map of biopsied regions within the primary tumor. Annotated pathology map displaying the precise localization of all biopsied regions within the primary prostate tumor. A comprehensive description of the pathology review and region annotation criteria is provided in the Methods section (Pathology Review), and detailed histomorphological features, including Gleason scores for each region, are provided in Table 1\nC3.4. Most of the patients analyzed had tumor specimens from multiple cores of the prostate and from the Lymphnodes but the analysis does not provide any hint on how prostate cancer cells seem to metastasize. The results could be masked by the fact that several sites in the primary tumors are taken into account in the analysis. Can the author perform the analysis taking the prostate cores and each LN site separately? Can the author infer that high “phenotypic” plasticity can be an account for, or evidence of monoclonal spreading and seeding ?\nOur analysis was performed at the individual sample level, that is, each regionally distinct sample was used to prepare the single-cell sequencing library independently. This design allows us to preserve and investigate site-specific variation in both clonal and phenotypic composition. In Supplementary Figure 5, we show per-sample pie charts depicting the distribution of genotypically inferred clones across each primary tumor locus and matched LN metastasis. To improve clarity, we have now revised the figure legend to read: “Pie charts depicting distribution of separately inferred genotypic clones for cells coming from each primary tumor loci and each matched lymph node metastases per patient.”\nSince the clonal architecture in LN metastases generally resembled that of the primary tumors, it suggested that multiple clones can seed and survive dissemination without undergoing a strong evolutionary bottleneck from a genetic perspective. This is consistent with previous findings from bulk sequencing 24–26. Thus, we sought out to identify common phenotypic shifts that may be shared across patients by considering data from all patients. This reinforces a role for non-genetic mechanisms, such as epigenetic variation, to support phenotypic plasticity for cancer cells to survive in changing environments and tissues. These observations are consistent with DNA barcode-based lineage tracing studies 27–2927,2829, which have shown that prostate cancer cells can undergo common phenotypic transitions under metastatic and therapeutic pressure without clonal selection. (Manuscript: page 11 lines 21-37: “Despite phenotypic and genotypic inter-patient variability, we found that a common phenotypic shift occurs in malignant epithelial cells disseminating to the lymph nodes across patients. The clonal architecture in lymph node metastatic cancer cells resembled that of the malignant epithelial cells in primary tumor loci in each patient. This suggests that clones can seed and survive the dissemination process in lymph nodes without undergoing an evolutionary bottleneck based on their CNV composition, consistent with bulk sample genetic sequencing studies 4,6,56. Additionally, barcode-based studies support a model of tumor evolution relying on phenotypic plasticity57,58. For instance, individual malignant prostate cancer cells, tagged with unique DNA barcodes, were shown to withstand selective pressures, such as metastatic dissemination and treatment resistance, by undergoing phenotypic transitions 59. While paired deep whole-genome sequencing data may reveal rare mutations that could impact phenotypic transitions, our observations support a contribution for non-genetic variation, including changes in chromatin accessibility in non-coding DNA elements, in defining the phenotype of disseminating malignant epithelial cells adapted to survive in sites beyond the primary tumor.”\nC3.5. Why AR and other TFs e.g. FOXA1 TFs were not enriched in the tumor cells in the Lymphnodes compared to primary regions?\nFigure 4 represents transcription factor motif enrichment within differentially accessible chromatin between primary and lymph node sites. AR (NR) family and FOXA1 (Forkhead family) are still detected as enriched factors (Figure 4F). Our findings imply increased accessibility around motifs of immune-based transcription factors such as the ETS family, rather than the lack of AR signalling in malignant cells in the LN. This also supports the differential transcriptional programs observed with RNA (Figure 4E).\nGiven that AR and FOXA1 are still observed in LN cells, our findings reflect an increase in immune-signalling co-occurring with continued reliance on AR as the core driver of prostate cancer. (Manuscript: page 12 lines 3-7 “Our investigation revealed phenotypic plasticity between primary and metastatic lymph node malignant cells were marked by a transition from prostate-specific pathways to immune/inflammatory-like states within the cancer cells themselves, suggesting engagement of epithelial-intrinsic immune programs during dissemination. ”)\nC3.6. The reviewer understands that the 4/5 patients were treatment näive but one of them was ADT pretreated. Is there any noticable difference between the tumor cells retrieved in this ADT treated patient lymphnodes compared to the tumor cells retrieved from the lymph nodes of the other treatment näive patients?\nWe thank the reviewer for this important question. Indeed, among the five patients in our cohort, one had received short-course neoadjuvant ADT prior to surgery. We agree that comparing the lymph node metastases from this patient to those from treatment-naïve individuals would be highly informative in understanding treatment-related effects on tumor phenotypes and clonal dynamics. While lymph node tissue was available for this ADT-treated patient, the nuclei isolation procedure unfortunately yielded insufficient quality for downstream single-cell profiling. As a result, only the primary tumor samples from this case were included in the current single-cell analysis. Therefore, our present dataset does not support a direct comparison of lymph node metastases between treatment-naïve and ADT-pretreated patients. (Manuscript Page 10 lines 22-26: “Our own cohort precluded contextualizing these gene sets with respect to drug treatment, as the only androgen deprivation therapy (ADT)-treated patient (ICGC_PCAL25), the LN samples had insufficient quality for downstream single-cell profiling. Thus, we turned to external single-cell RNA sequencing data from metastatic prostate cancer 48. [...]”\nC3.7. The authors focus their analysis on the epithelial-like cells but other cell types present in the specimens seem to be ignored. In light of the tumor (immune) microenvironment shaping the metastatic potential, can the author distinguish more immune active vs colder tumors and foci and their propensity and features leading to tumor cell dissemination? Are these just simply the result of higher inflammation on the primary site?\nRecent studies have classified tumors into immune “hot” versus “cold” states based on T cell infiltration, PD-L1 (CD274) expression, and tumor mutational burden (TMB) 30. Prostate cancer is considered an immune-cold tumor type 31, and this is supported by our own single-cell data, where T cells constituted only ~1.95% of all profiled cells, limiting robust T cell-based stratification across patients. With this caveat in mind, we decided to assess CD274 (PD-L1) expression within malignant epithelial cells across tumor and lymph node (LN) foci. Gene expression was log-normalized using the LogNormalize method in Seurat before visualization. Overall, CD274 expression was low across all sites from all patients (Rebuttal Figure 3.7A), consistent with the immune-cold phenotype of prostate cancer.\nRebuttal Figure 3.7: PD-L1 (CD274) expression across tumor loci. Violin plots displaying the mean log-normalized expression of PD-L1 (CD274) across all tumor loci for each individual patient.\nReviewer #3 (Remarks on code availability):\nAt this stage I did not review the code.\nReviewer #4, expertise in prostate cancer epigenetics and mechanisms of lineage plasticity (Remarks to the Author):\nIn this manuscript, Keshavarzian et al employs scMultiome to dissect intra- and inter-tumoral heterogeneity in high-risk prostate cancer across spatially distinct primary tumor loci and matched locoregional lymph node (LN) metastases. This research show that Inter-patient heterogeneity is largely restricted to malignant epithelial cells and LN metastases mirror the clonal diversity of primary tumors but undergo a transcriptional shift toward inflammatory-like states. This work advances the understanding of prostate cancer progression by linking clonal architecture to phenotypic plasticity. The inflammatory transition paradigm also provides a framework for targeting non-genetic adaptation mechanisms in metastasis. The following points need to be considered.\nMajor comments:\nC4.1. The cohort size and clinical annotation are limited in this study which includes only 5 patients, with one receiving neoadjuvant therapy. The author can explicitly discuss how cohort size limits statistical power for inter-patient comparisons and include clinical outcomes (e.g., recurrence, survival) if available to contextualize the observed heterogeneity.\nWhile this study was designed as an in-depth exploratory analysis of spatial heterogeneity using single-cell multiomic profiling in a rare and clinically relevant setting (paired primary and locoregional metastatic tissues from high-risk prostate cancer patients), we agree that a larger cohort size would have been beneficial. Our current cohort includes five patients, the extensive spatial sampling per patient (up to 10 biopsies per individual) and the high resolution of ~300,000 single-nucleus profiles which enabled robust intra-patient comparisons and allowed us to identify conserved phenotypic trends across patients despite inter-individual variation.\nNonetheless, we acknowledge that the small number of patients limits the statistical generalizability of inter-patient comparisons. We now explicitly discuss this limitation in the revised Discussion section and caution against over-interpreting patient-specific findings as broadly representative. Clinical annotations, such as histopathologic features (e.g., Gleason score), were leveraged to contextualize some findings, such as ICGC_PCAL37 in the Results Section “Malignant cells disseminating to locoregional lymph nodes undergo an inflammatory-like cell state transition despite genotypic similarity to primary cancer”; page 9 lines 3-8 and page 10 lines 7-10. Our cohort, which is limited in size and restricted to lymph node metastasis, constraints our ability to expand our findings to distal metastases and contextualize clinical outcomes.\nWe agree that future studies with larger cohorts and longitudinal follow-up data will be essential to link the observed cellular heterogeneity to patient outcomes and treatment response (raised on page 12, lines 21-25)\nC4.2. ICGC_PCAL37, which has the lowest Gleason score, has some specificity in some aspect, such as lower tumor heterogeneity and the unique inter-site difference. These may indicate that the stage of tumor progression can have an impact on the inter- and intra-tumor heterogeneity. The author can discuss how tumor progression relates to the heterogeneity of prostate tumor inter- and intra-site through the analysis of wider range of statistic data from their own or in the public database\nExtensive prior work has shown that prostate cancer displays marked inter-patient heterogeneity across histo-morphological and molecular dimensions in both primary and metastatic disease settings 32–37. Large-scale genomic studies have further delineated molecular subtypes of prostate cancer that differ at the genetic level despite sharing similar histopathological grades 38–44.\nIn our cohort, patient ICGC_PCAL37, the only case with a Gleason score of 7 (3+4), stood out in several respects. This individual exhibited comparatively lower tumor heterogeneity and unique inter-site differences. Clinically, this patient had no evidence of bone metastases (similar to ICGC_PCAL10), experienced biochemical recurrence relatively early, and remained stable for an extended period following salvage radiation therapy before initiating systemic treatment nearly five years later (Table 1: baseline_metadata, therapy_follow_up). These characteristics may reflect a more indolent disease course and support the idea that earlier-stage tumors may harbor lower levels of molecular heterogeneity.\nWhile our dataset includes only five patients, limiting the statistical power to formally compare Gleason score with heterogeneity metrics, previous studies support this link. For example, high-risk prostate cancers with elevated Gleason scores tend to exhibit more pronounced somatic copy number alterations (SCNAs) and greater genomic instability 45, Figure 2). Similarly, The Cancer Genome Atlas44 reported that tumors with high SCNA burden are associated with higher Gleason scores and elevated PSA levels (Figures S3B–S3D).\nFurthermore, the degree of structural rearrangement has been shown to increase across Gleason grades: tumors with any Gleason pattern 4 component (e.g., 3+4 or 4+3) demonstrate significantly higher rearrangement rates than Gleason 3+3 tumors 44, Figure 1). These studies collectively support the hypothesis that tumor progression correlates with increasing inter- and intra-tumor heterogeneity.\nC4.3. Bone metastases are very common in prostate cancer. Does prostate bone metastases also show an immune/inflammatory-like states transition similar with LN metastases? The author can discuss the potential universality of this inflammatory-like states.\nAs our cohort includes only lymph node metastases, we are unable to directly assess whether bone metastases exhibit a similar immune/inflammatory-like state transition. However, prior work by others 46 has examined progression biopsies from mCRPC patients, including bone metastases, and similarly reported enrichment of immune/inflammatory-related pathways, such as interferon signaling and IL6/JAK/STAT—during treatment resistance. This suggests that such inflammatory-like states may represent a shared, site-independent mechanism associated with lineage plasticity. While distant metastatic organs such as bone and liver are considered aggressive metastatic disease, locoregional lymph nodes may similarly capture this evolving aggressive biology due to ongoing clonal seeding from the primary site prior to therapy 24. Thus, although our analysis is limited to lymph node metastases, the inflammatory-like transcriptional states we observe may reflect a generalizable mechanism of metastatic adaptation and resistance. sWe have added a brief statement in the Discussion to reflect this point (page 12, lines 3-14): “These immune/inflammatory-like pathways, including JAK-STAT signaling, have also been observed in mCRPC progression samples across diverse metastatic sites, including bone 63. JAK-STAT signalling is also implicated as a driver of plasticity and dedifferentiation from luminal-epithelial cell states in prostate cancer models 64,65. This suggests that activation of these pathways may represent a mechanism of progression and potentially contribute to treatment resistance, though this hypothesis warrants further validation. ”\nLastly, we used our chromatin variants to get their global enrichment over mCRPC bone metastases profiled using bulk ATAC-sequencing. We observed the primary Tu chromatin variants to score lower than the LN chromatin variants (Wilcoxon P=0.008), with respect to the average chromatin accessibility profile across mCPRC samples. However, these samples were profiled using bulk ATAC-sequencing which implies that the computed signal is inclusive of cells in the tumor microenvironment. For this reason, we turned to performing the same analyses on over prostate cancer cell lines which are composed of purely malignant cells. We saw that LN chromatin variants were more enriched in aggressive cell lines (PC3s and Du145, Wilcoxon P=0.029) relative to Tu chromatin variants which were more pronounced in less aggressive hormone dependent cell lines (LNCaP, 22RV1, VCaPs, Wilcoxon P=0.008).\nRebuttal 4.3 LN chromatin variant signatures are more enriched relative to primary Tu in aggressive prostate cancers. A) Chromatin variant enrichments in 28 bulk ATAC-seq mCRPC bone metastasis samples. Enrichment scores computed using ChromVar, which assesses the relative enrichment of all differential LN (N=5544, Log Fold Change > 0.5, adjusted P-val < 0.1) and Tu (N=5312, Log Fold Change > 0.5, adjusted P-val < 0.1) peaks on a consensus mCRPRC peakset. B) Chromatin variant enrichments computed in a similar manner on a consensus cell line peakset, with two replicates per cell line with the exception of VCaPs. C) Depiction of individual cell line data and replicates in (B) as a heatmap, color-scale is limited to a minimum and maximum of -10 and +10 Z-scores.\nC4.4. While the transcriptional shifts in LN metastases are compelling, the functional consequences of inflammatory reprogramming remain speculative. In vitro assays (e.g., cytokine exposure) or spatial proteomics could help in validating inflammatory signatures as well as the cellular function of this transition.\nWhile the functional consequences of inflammatory reprogramming remain to be fully elucidated, we partially addressed this by validating our lymph node (LN)-upregulated transcriptional signature using an independent single-cell RNA-seq dataset from metastatic prostate cancer across treatment stages 47. We observed an enrichment of LN-up genes during disease progression and treatment resistance (Figure 4g), supporting their biological relevance. Similar pathways were also shown to be upregulated in mCRPC progression post-Enzalutamide treatment 46, and aggressive prostate cancer cells expressing immune-like markers were previously observed in mouse models 48.\nWe agree that spatial transcriptomic and proteomic technologies could provide further insight into the local context and function of these states. Similarly recent technological advances, such as spatial scATAC-seq and in situ tagmentation, highlight promising future directions for dissecting spatial epigenomic heterogeneity 49,50.We believe our findings help lay the groundwork for such future studies.\nWe also acknowledge the potential value of in vitro cytokine exposure assays to interrogate the functional consequences of inflammatory reprogramming. However, given the broad scope of our transcriptional signatures, encompassing multiple pathways and hundreds of genes, designing targeted perturbation experiments would require substantial additional work to deconvolute pathway-specific effects. As with IHC-based validation, these experiments, while valuable, fall beyond the scope of the current study.\nC4.5. In this research, sampling focuses on primary tumors and LN metastases but lacks the adjacent normal tissue comparisons to define tumor-specific changes. Spatial transcriptomics may help in mapping phenotypic transitions at the invasion front.\nWhile our manuscript was focused on investigating the primary and metastatic heterogeneity, we did not have access to matched adjacent normals to define tumor-specific changes. However, we utilized the identified tumor-microenvironment cell types as a reference point as discussed in our response to comment C2.7. We have additionally covered how our data highlights spatial transcriptomics to be imperative to our understanding of prostate cancer heterogeneity and metastasis (See response to comment C1.2 and revised Discussion page 11 line 6-20)\nMinor comments:\nC4.6. In line 335, the author claims that “lymph node-specific malignant epithelial populations consistently ranked higher than primary tumor-specific malignant epithelial populations across all patients”. However, in ICGC_PCAL37 the rank of the two population looks the same.\nPatient ICGC_PCAL37 had the most inter-mixing across sites, with a single lymph-node specific malignant epithelial population (as reported in Results page 9 lines 3-8), hence, as the reviewer points out, we were not able to draw strong conclusions from this patient. We have revised our paper to more accurately reflect our finding (page 10 lines 7-10) “When using lymph node-upregulated genes, lymph node-specific malignant epithelial populations ranked higher than primary tumor-specific malignant epithelial populations across all patients, except for our lowest Gleason Score patient (ICGC_PCAL37) (Supplementary Fig. 8a). ”\nC4.7. Supplementary figure 7 is incorrectly labeled as supplementary figure 5.\nThank you for bringing this to our attention, we have now revised our figure numbers.\nC4.8. The label of the lateral axis of figure 1D at the top of the chart is a little confusing. The author should explain more clearly.\nWe have modified our figure caption to better convey what Z-score deviations represent.\nRevised Figure 1 Caption:\n“D) Selected transcription factors’ motif deviations according to Chromvar (right panel), values represent the average z-score of motif accessibility deviations across all cells of a given cell type. Higher z-scores indicate greater enrichment of TF motifs in accessible chromatin, compared to background accessibility.”\nC4.9. In the third panel of supplementary figure 5 showing ICGC_PCAL34, the population marked pink does have a corresponding caption on the left.\nWe have corrected the figure in the revised version.\nReferences:\n1. Elyada, E. et al. Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts. Cancer Discov. 9, 1102–1123 (2019).\n2. Laisné, M., Lupien, M. & Vallot, C. Epigenomic heterogeneity as a source of tumour evolution. Nat. Rev. Cancer 25, 7–26 (2025).\n3. Whiting, F. J. H., Househam, J., Baker, A.-M., Sottoriva, A. & Graham, T. A. Phenotypic noise and plasticity in cancer evolution. Trends Cell Biol. 34, 451–464 (2024).\n4. Stergachis, A. B. et al. Developmental Fate and Cellular Maturity Encoded in Human Regulatory DNA Landscapes. Cell 154, 888–903 (2013).\n5. Stewart-Morgan, K. R., Petryk, N. & Groth, A. Chromatin replication and epigenetic cell memory. Nat. Cell Biol. 22, 361–371 (2020).\n6. Lupien, M. et al. Growth factor stimulation induces a distinct ER(alpha) cistrome underlying breast cancer endocrine resistance. Genes and Development 24, 2219–2227 (2010).\n7. Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e29 (2021).\n8. Liu, Z. et al. Single-Cell Chromatin Accessibility Analysis Reveals the Epigenetic Basis and Signature Transcription Factors for the Molecular Subtypes of Colorectal Cancers. Cancer Discov. 14, 1082–1105 (2024).\n9. Regner, M. J. et al. A multi-omic single-cell landscape of human gynecologic malignancies. Mol. Cell 81, 4924–4941.e10 (2021).\n10. Terekhanova, N. V. et al. Epigenetic regulation during cancer transitions across 11 tumour types. Nature 623, 432–441 (2023).\n11. Bejarano, L., Jordāo, M. J. C. & Joyce, J. A. Therapeutic Targeting of the Tumor Microenvironment. Cancer Discov. 11, 933–959 (2021).\n12. Arora, R. et al. Spatial transcriptomics reveals distinct and conserved tumor core and edge architectures that predict survival and targeted therapy response. Nat. Commun. 14, 5029 (2023).\n13. Wang, X. et al. Spatial transcriptomics reveals substantial heterogeneity in triple-negative breast cancer with potential clinical implications. Nat. Commun. 15, 10232 (2024).\n14. Oliveira, M. F. de et al. High-definition spatial transcriptomic profiling of immune cell populations in colorectal cancer. Nat. Genet. 57, 1512–1523 (2025).\n15. Thurman, R. E. et al. The accessible chromatin landscape of the human genome. Nature 489, 75–82 (2012).\n16. Li, C. et al. Single-cell brain organoid screening identifies developmental defects in autism. Nature 621, 373–380 (2023).\n17. Weinreb, C., Wolock, S., Tusi, B. K., Socolovsky, M. & Klein, A. M. Fundamental limits on dynamic inference from single-cell snapshots. Proc. Natl. Acad. Sci. U. S. A. 115, E2467–E2476 (2018).\n18. Gorin, G., Fang, M., Chari, T. & Pachter, L. RNA velocity unraveled. PLoS Comput. Biol. 18, e1010492 (2022).\n19. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386 (2014).\n20. Eppert, K. et al. Stem cell gene expression programs influence clinical outcome in human leukemia. Nat. Med. 17, 1086–1093 (2011).\n21. Bonnet, D. & Dick, J. E. Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell. Nat. Med. 3, 730–737 (1997).\n22. Pomerantz, M. M. et al. Prostate cancer reactivates developmental epigenomic programs during metastatic progression. Nat. Genet. 52, 790–799 (2020).\n23. Derrien, J. et al. Acquired resistance to a GPRC5D-directed T-cell engager in multiple myeloma is mediated by genetic or epigenetic target inactivation. Nat. Cancer 4, 1536–1543 (2023).\n24. Warner, E. W. et al. Multiregion sampling of de novo metastatic prostate cancer reveals complex polyclonality and augments clinical genotyping. Nat Cancer 5, 114–130 (2024).\n25. Singhal, U. et al. Integrative multi-region molecular profiling of primary prostate cancer in men with synchronous lymph node metastasis. Nat. Commun. 15, 4341 (2024).\n26. Reiter, J. G. et al. Lymph node metastases develop through a wider evolutionary bottleneck than distant metastases. Nat. Genet. 52, 692–700 (2020).\n27. Rehman, S. K. et al. Colorectal Cancer Cells Enter a Diapause-like DTP State to Survive Chemotherapy. Cell 184, 226–242.e21 (2021).\n28. Berthelet, J. et al. The site of breast cancer metastases dictates their clonal composition and reversible transcriptomic profile. Sci Adv 7, (2021).\n29. Serio, R. N. et al. Clonal Lineage Tracing with Somatic Delivery of Recordable Barcodes Reveals Migration Histories of Metastatic Prostate Cancer. Cancer Discov. (2024) doi:10.1158/2159-8290.CD-23-1332.\n30. Wu, B., Zhang, B., Li, B., Wu, H. & Jiang, M. Cold and hot tumors: from molecular mechanisms to targeted therapy. Signal Transduct. Target. Ther. 9, 274 (2024).\n31. Kwon, W.-A. & Joung, J. Y. Immunotherapy in prostate cancer: From a ‘cold’ tumor to a ‘hot’ prospect. Cancers (Basel) 17, 1064 (2025).\n32. Network, The Cancer Genome Atlas Research et al. The Molecular Taxonomy of Primary Prostate Cancer. Cell 163, 1011–1025 (2015).\n33. Stopsack, K. H. et al. Oncogenic Genomic Alterations, Clinical Phenotypes, and Outcomes in Metastatic Castration-Sensitive Prostate Cancer. Clin. Cancer Res. 26, 3230–3238 (2020).\n34. Mateo, J. et al. Genomics of lethal prostate cancer at diagnosis and castration resistance. J. Clin. Invest. 130, 1743–1751 (2020).\n35. Løvf, M. et al. Multifocal Primary Prostate Cancer Exhibits High Degree of Genomic Heterogeneity. Eur. Urol. 75, 498–505 (2019).\n36. Tang, D. G. Understanding and targeting prostate cancer cell heterogeneity and plasticity. Semin. Cancer Biol. 82, 68–93 (2022).\n37. Haffner, M. C. et al. Genomic and phenotypic heterogeneity in prostate cancer. Nat. Rev. Urol. 18, 79–92 (2021).\n38. Cancer Genome Atlas Research Network. The Molecular Taxonomy of Primary Prostate Cancer. Cell 163, 1011–1025 (2015).\n39. Stopsack, K. H. et al. Oncogenic Genomic Alterations, Clinical Phenotypes, and Outcomes in Metastatic Castration-Sensitive Prostate Cancer. Clin. Cancer Res. 26, 3230–3238 (2020).\n40. Mateo, J. et al. Genomics of lethal prostate cancer at diagnosis and castration resistance. J. Clin. Invest. 130, 1743–1751 (2020).\n41. Small, E. J. et al. Characterization of neuroendocrine prostate cancer (NEPC) in patients with metastatic castration resistant prostate cancer (mCRPC) resistant to abiraterone (Abi) or enzalutamide (Enz): Preliminary results from the SU2C/PCF/AACR West Coast Prostate Cancer Dream Team (WCDT). J. Clin. Orthod. 33, 5003–5003 (2015).\n42. Gerhauser, C. et al. Molecular Evolution of Early-Onset Prostate Cancer Identifies Molecular Risk Markers and Clinical Trajectories. Cancer Cell 34, 996–1011.e8 (2018).\n43. Woodcock, D. J. et al. Genomic evolution shapes prostate cancer disease type. Cell Genom 4, 100511 (2024).\n44. Fraser, M. et al. Genomic hallmarks of localized, non-indolent prostate cancer. Nature 541, 359–364 (2017).\n45. Espiritu, S. M. G. et al. The evolutionary landscape of localized prostate cancers drives clinical aggression. Cell 173, 1003–1013.e15 (2018).\n46. Westbrook, T. C. et al. Transcriptional profiling of matched patient biopsies clarifies molecular determinants of enzalutamide-induced lineage plasticity. Nat. Commun. 13, 5345 (2022).\n47. Hawley, J. E. et al. Anti-PD-1 immunotherapy with androgen deprivation therapy induces robust immune infiltration in metastatic castration-sensitive prostate cancer. Cancer Cell 41, 1972–1988.e5 (2023).\n48. Morel, K. L. et al. Low tristetraprolin expression activates phenotypic plasticity and primes transition to lethal prostate cancer in mice. J. Clin. Invest. 135, (2024).\n49. Llorens-Bobadilla, E. et al. Solid-phase capture and profiling of open chromatin by spatial ATAC. Nat. Biotechnol. 41, 1085–1088 (2023).\n50. Noberini, R. et al. Spatial epi-proteomics enabled by histone post-translational modification analysis from low-abundance clinical samples. Clin. Epigenetics 13, 145 (2021).\nReviewer #1 (Remarks to the Author):\nThe authors have addressed my comments effectively.\nResponse: Thank you for acknowledging our efforts to address your previous comments effectively, which enabled us to enhance the scientific content of our manuscript.\nReviewer #3 (Remarks to the Author):\nThe authors clarified a sufficient amounts of critical points and I am convinced of the robust and rigorous analysis of the data. Maybe, given the lack of further validation that prostate cancer cells converge to a phenotype similar to an inflammatory-like status in the LNs, the title of the manuscript could be slightly modified.\nResponse: We thank the reviewer for recognizing the rigor and robustness of our work. We have revised the manuscript (title and other sections) as suggested.\nReviewer #4 (Remarks to the Author):\nThe authors have answered my questions.\nResponse: Thank you for recognizing the quality and impact of our study.\nReviewer #2 (Remarks to the Author):\nAuthors have partially addressed my concerns related to the manuscript. Still major issues remain.\n1. For CNV analysis authors could consider more optimized (e.g. python) implementation of the algorithms. None of the less, authors do demonstrate data to alleviate concern regarding subsampling. However, how the choice of reference cell type impacts the CNV calls has not been properly addressed. As different cell types expressed different set of genes this can have fundamental impact on inferred CNVs and thus on conclusions.\nResponse: We agree with the reviewer that the majority of CNV calling algorithms are based solely on RNA expression and are therefore sensitive to the choice of reference. We did not experience this to be a problem for Numbat because it relies on both haplotype frequency (reference-independent) and molecular abundance (reference-dependent).\nTo this point, a very recent CNV-calling benchmarking paper showcases that with the exception of Numbat, most other CNV callers are sensitive to the choice of reference (https://doi.org/10.1038/s41467-025-62359-9). This paper demonstrated Numbat’s high performance in\n1. Outperforming all other methods at detecting the correct CNVs regardless of the reference cell type used\n2. Separating tumor and normal cells\n3. Identifying the right clonality\nThere are four ways that references can be chosen that deliver equivalent performance when using Numbat, including the Within-dataset, different normal cell types. Our results are based on this mode of reference choice, by using normal cell types found in our own dataset. This mode is recommended over using external-dataset modes, even if the external datasets have matched normal cell types. We had also consulted with the Numbat first author, and their reference choice for their datasets are also from different cell types coming from the same dataset, as their internal experimentations also proved this to be a more effective method. The benchmarking paper recommends this choice above matched normal cell types from external datasets as well.\n- External-dataset, matched normal cell types: According to the benchmarking paper, utilizing a within-dataset, different normal cell type is preferred. Additionally, we were unable to find publicly released single-nucleus based RNA ‘normal’ epithelial prostate cells. The released datasets were on single-cell (sc) RNA sequencing, which has a different read- distribution than single-nucleus (sn) RNA sequencing. sn versus sc conducted by the same groups on healthy cells show a moderate difference between the two experimental setups (https://doi.org/10.1371/journal.pone.0209648, https://doi.org/10.1371/journal.pgen.1011301). Thus, introducing a non-trivial difference in signal due to external datasets. Hence why we chose to use internal-dataset, different cell types.\n- External-dataset, different normal cell types: We did not test this as we believe it to not align with neither ours or reviewer’s interest, although this data is reported in the benchmarking paper.\nWe appreciate the reviewer’s concern in the choice of the reference, as other non-haplotype-dependent CNV inference methods are highly sensitive to the reference choice. However, given the strong evidence from benchmarking papers, our internal experimentation, and consultations with the main Numbat group, and the benchmarking paper, we have adopted a scientifically sound setup for our analyses.\nWith respect to a more optimized (e.g. Python) implementation of Numbat, we have thoroughly considered it and gone through the entire code-base. To clarify, we previously mentioned that the pairwise distance matrix computation for phylogenetics (ape::nj) grows quadratically with the number of supplied cells. The function’s core is implemented in C++ rather than R, which is significantly faster than both Python and R. Given that the highest performing functions in Python are also C++ based, it is simply a compute limitation despite allocating 2TB of memory. Other parts of Numbat utilize R’s data.table which has a generally faster performance than much of Python’s data processing libraries. Furthermore, according to the benchmarking paper, Numbat’s resource uptake (in memory/GB) is better than most CNV-calling algorithms.\n2. While authors do present innovative comparisons through TF motif enrichment, I do not agree that this would remove the need to address biology behind sites through standard approaches. While motif analysis has novelty, results need to be presented and put into context of differential expression (and possibly accessibility) analysis and pathway activation, that I would consider as baseline benchmark for any more complex analysis.\nResponse: To address the reviewer’s request for a deeper integration of differential gene expression, chromatin accessibility, pathway activity, and TF motif enrichment, we expanded our analyses using SCENIC+ (single-cell enhancer–gene regulatory network) inference (doi:10.1038/s41592-023-01938-4). As now described in the revised Methods section (“Assessing the degree of phenotypic heterogeneity in locus-specific populations,” p. 21, line 31–p. 22, line 4), we applied SCENIC+ to seven tumor samples (ICGC_PCAL34: tumors 1, 2, 4, 10; ICGC_PCAL03: tumors 6, 8, 9).\nThis algorithm links accessibility, genes, and transcription factors as triplet regulons, and derives the most specific, high-confidence regulons for each cluster. Running this analysis on the seven tumors individually revealed transcription factors and genes specific to low entropy (Fig. 2, Supplementary Fig. 3) malignant epithelial clusters with a sufficient number of cells (N>1000). Since the malignant clusters are specific to each tumor, the high specificity scores of the listed transcription factors and genes for each cluster is representative of intra-tumor heterogeneity (manuscript, newly added Table 2, Cluster_eRegulons). The reported transcription factors by this algorithm include FOX, HOX, FOS/JUN, GRHL, and NR-based TFs.\nThis reproducibly confirmed our reported sets of TFs driving variance using an orthogonal approach (Fig. 2d, Supplementary Fig. 4).\nGiven that SCENIC+ also links TFs/accessibility to genes, we also performed pathway enrichment analyses on the lists of genes specific to a given cluster from a tumor (newly added Table 2, Cluster_eRegulons). Performing pathway enrichment analysis on the lists of genes specific to a given cluster, regulated by the reported TFs, revealed distinct and overlapping pathways (FDR < 0.05) (newly added Supplementary Fig. 5, Table 2: Hallmark Pathways). The commonly enriched pathways included androgen and estrogen response, whereas more specific pathways involved EMT, hypoxia, and TGF-β signaling in ICGC_PCAL34, along with p53, TGF-β, PI3K/AKT/mTOR Signaling, mTORC1 and DNA repair pathways in ICGC_PCAL03 (Results, p. 6, line 31-p. 7, line 3).\n3. Authors have not properly addressed their claims of plasticity. In my view, raising plasticity as a key concept, will require authors to properly address the transcriptomic (and/or epigenetic) states (and provide definition of cell state) in relation to CNV based clones. Current revision, where non-formal definition of plasticity is given at the end of discussion is far from sufficient.\nResponse: We thank the reviewer for highlighting the need to better define and contextualize plasticity in relation to transcriptomic and epigenomic states across CNV-defined clones. We have revised the manuscript to more clearly articulate this relationship and to restrain our use of the term plasticity throughout.\nWhile our data, relying on lineage-tracing using CNV-defined clones combined with single-cell RNA and ATAC data from the primary site of tumour development and a secondary site of colonization, supports a role for phenotypic adaptation, we acknowledge the lack of a genetically evolving barcode-based in vitro model that would allow definitive single-cell resolution lineage tracing. Such an approach lies beyond the scope of the present study. Accordingly, we have removed the term “phenotypic plasticity” from the title and all other sections of the manuscript, retaining it only in the discussion when referring to relevant manuscripts by others demonstrating that metastatic potential in prostate cancer can emerge through transcriptional and epigenetic reprogramming rather than de novo clonal selection (manuscript refs 4,6,56–58)."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Authors have addressed all my concerns in their revision."}, "author_response": ""}]

Decision letter