Clonal evolutionary analysis reveals patterns of malignant transformation of Intraductal Papillary Mucinous Neoplasms of the pancreas

Domain: chemistry
Published: 2026-03-04
Editorial outcome: major_revision

Abstract
Intraductal papillary mucinous neoplasms (IPMNs) are critical precursors to pancreatic ductal adenocarcinoma, a highly lethal cancer characterized by late detection and rapid progression. Here we integrate multi-region whole-genome and transcriptome sequencing to trace the evolution of IPMN, constructing detailed phylogenetic trees to provide insights into subclonal architectures and progression pathways. Our analysis identifies two distinct evolutionary trajectories: one driven by a single ancestral clone, and another involving multiple independent ancestral clones. We further explore the roles of mutational signatures and structural variants in promoting clonal evolution and the emergence of new subclones. Complementing these genomic findings, our transcriptomic analysis reveals unique gene expression profiles and variations in the immune landscape that correlate with different progression stages of IPMN. These insights reveal the complex molecular dynamics of IPMN heterogeneity and progression, highlighting the need to refine early detection and treatment strategies.

Introduction
Pancreatic ductal adenocarcinoma (PDAC) is associated with a dismal prognosis and is projected to be the second leading cause of cancer related deaths soon1. This poor prognosis is largely attributed to the late-stage presentation, which is further compounded by rapid disease progression once diagnosed.

As one of the precursor lesions of pancreatic cancer, intraductal papillary mucinous neoplasm (IPMN) are cystic lesions that are radiologically detectable offering a unique opportunity for detection and intervention to prevent malignant transformation. Current clinical recommendations for surgical resection of IPMNs are primarily based on imaging criteria, leading to large numbers of patients undergoing annual imaging surveillance2. A better understanding of the critical molecular events that drive progression in premalignant lesions is urgently needed to better manage these patients to achieve early diagnosis and cancer prevention.

Genomic studies on PDAC have advanced our understanding of its molecular landscape3-6. However, these studies have often focused on single tumor regions without integrating detailed pathological or radiological data, resulting in cohort analyses that combine PDACs arising from various precursor lesions. Studies modelling progression from IPMN to PDAC, using single driver mutations, have revealed significant intratumoral heterogeneity, frequently identifying distinct KRAS and GNAS mutations within the same lesion, particularly in low-grade dysplasia regions7-9,10.

This heterogeneity suggests polyclonal origins in IPMNs, contrasting with the clonal homogeneity observed in advanced PDAC11, reflecting distinct evolutionary patterns underlying early and late stages of tumorigenesis. These findings highlight the need for comprehensive analyses of the genetic complexity within individual lesions, as the contributions of other genomic alterations, such as structural variants (SVs) and copy number alterations (CNAs), remain poorly understood.

This study employs multi-regional whole-genome and transcriptomic sequencing on surgically resected IPMN samples to comprehensively analyse intratumoral heterogeneity, clonal evolution, and associated molecular pathways. By integrating genomic and transcriptomic data, we provide insights into the evolutionary trajectories and signaling networks underlying IPMN progression, advancing our understanding of the pathways leading to pancreatic cancer.

Results
Patient characteristics
A cohort of 12 patients underwent pancreatic resection for high-risk IPMNs following clinical guidelines (Supplementary Data 1). Forty-seven tumor samples, plus matching normal tissue from each patient, were harvested from 12 surgical specimens (mean of 4 regions per tumor, range 2–6). Six cases (Cases 6, 9, 7, 4, 15, 16) provided samples representing different grades of dysplasia and histologically associated cancer within the same surgical specimen. Whole-genome sequencing (WGS) was successfully performed on 54 samples (42 tumor samples, mean coverage 80X, and 12 normal samples at 60X). Thirty-seven tumor samples were included in the clonal analysis, excluding 5 with low tumor purity or mutation counts. Transcriptome sequencing (RNAseq; average 100 million paired reads) was performed on 36 tumor samples, 32 of which had matching WGS data (Figure S1A2).

Mutational landscape of IPMN
A total of 66,724 single-nucleotide variants (SNVs), 4,683 insertions/deletions (Indels), and 2,447 structural variants (SVs) were identified across all tumor samples. These spanned 48 driver genes and 11,200 other genes across IPMN and associated PDAC samples. Tumor mutation burden (TMB) (4.11 vs 1.38, p = 0.0001; Figure S1B1a, Supplementary Data 2), as well as the mean number of SNVs (3,558 vs. 1,608, p = 0.0038), Indels (253 vs. 105.5, p = 0.009), and SVs (91.8 vs. 31.7, p = 0.0004), were significantly higher in PDAC compared to IPMN. Higher SV counts were also observed for large translocations, tandem duplications, deletions, and inversions in PDAC (Figure 1B and S1B1c, Supplementary Data 3).

Key driver mutations were observed in tumor suppressors (ARID1A, ATM, CDKN2A, LRP1B, NBEA, RNF43, SMAD4, TP53) and oncogenes (FAT3, GNAS, KMT2C, KRAS, ZNF521). KRAS was the most frequently mutated gene, with hotspot mutations (15/22 G12D, 5/22 G12V, 2/22 G12R) present in both IPMN and PDAC samples and no significant differences in term of cancer cell fractions (CCFs). Case 4 and 15 showed two distinct KRAS variants within the same lesion, suggesting intralesional heterogeneity (Figure S1B2-3, Supplementary Data 4). TP53 SNVs/Indels were predominantly seen in high-grade IPMN and PDAC. LRP1B mutations were exclusive to invasive cancer (0 in IPMN vs 7 in PDAC, p = 0.0056), while GNAS mutations were more frequent in low-grade IPMN (8 in IPMN vs 1 in PDAC, p = 0.0297) (Figure S1B4).

Driver SVs in exonic regions were detected in Case 7 and Case 16, including alterations in MUC4 in high-grade IPMN and PDAC, MDM2 in both low- and high-grade lesions, and MYB, THRAP3, and PDCD1LG2 in high-grade IPMN and PDAC (Figure 1C).

High-frequency CNAs were identified in 64 regions (33 gains, 18 losses, and 13 complex regions). The most frequent regions were 17p11.2 (gain 13/41, loss 12/41, complex 15/41), 20q11.21 gain (36/41), and 9p21.3 loss (25/41). Most CNAs (52/64) persisted from non-invasive IPMN to PDAC (Figure S1C.A), with four gains and four losses showing higher prevalence in PDAC. Interestingly, 17p11 gain was more common in IPMN than PDAC (Figure S1C.B, Supplementary Data 5). Losses in 11 driver genes were frequent, including CDKN2A (9p21.3), LAT1S1, ARID1B, SMAD4, U2AF1, and RNF43, with RNF43 and U2AF1 losses significantly associated with high-grade lesions (p = 0.0086 and p = 0.0049 respectively) (Figure S1B4).

Mutational signatures in IPMN
Mutational signature analysis identified eight single-base substitution (SBS) signatures, with SBS1 and SBS5 (age-related) predominant across all samples. AID/APOBEC-related SBS2 and SBS13 emerged only in clones of Case 7 progressing to cancer. Indel signatures ID1 and ID2, linked to mismatch repair deficiency, were present in cases 12 and 7 respectively, with five additional indel signatures detected. However, all cases retained expression of MMR proteins by immunohistochemistry, and no molecular evidence of microsatellite instability—including low overall tumor mutational burden—was observed, suggesting that ID1 and ID2 may develop through MMR-independent mechanisms12.

Structural signature SV4 was significantly enriched in PDAC samples (p = 0.011). Copy number signature CN9, indicating chromosomal instability, was observed in 31 samples, while a CN48B signature, marked by LOH segments, was identified in 21 samples from 10 cases (Figure S2A-B). Notably, CN48B activity showed a moderate, statistically significant positive correlation with ID1 and ID2 (R = 0.53 and 0.55; p = 0.003 and 0.0002, respectively), suggesting that CN48B may reflect broader defects in DNA repair or replication fidelity (Figure S2B–C).

Heterogeneity in genomic evolution of IPMNs
Mutations were clustered according to their CCFs to construct phylogenetic trees from multiregional samples of 11 cases, representing for each case the genetic relationships among subclones based on the hierarchy of oncogenic events. The analysis revealed two distinct evolutionary trajectories: (1) all IPMNs develop from one single most recent common ancestor (MRCA) (N = 8; Cases 3, 6, 7, 9, 10, 12, 13, 16), and (2) IPMNs develop from separate, independent MRCAs (N = 3; Cases 2, 4 and 15).

Further heterogeneity in tumor evolution was observed within both single and multiple independent MRCA trajectories. In the single MRCA cases, driver alterations (including SNVs, Indels, and CNAs) occurred clonally and subclonally in both oncogenes and tumor suppressor genes (Cases 3, 6, 7, 9, 10, and 12). In contrast, in Cases 13 and 16, driver alterations were confined to subclones.

Among the multiple independent MRCA cases, Case 15 exhibited IPMN_HGD and PDAC samples initiated from distinct MRCAs, both harboring driver mutations, while other clones without driver mutations remained inactive in subclonal progression. Similar independent evolutionary patterns were observed in Cases 2 and 4.

Across all samples, subclonal progression revealed uneven branch lengths in the phylogenetic trees of Cases 6, 7, 9, and 10. Notably, Case 6 displayed 1,932 mutations in the red branch compared to only 32 in the green branch, suggesting a late subclonal sweep in the red branch (Figure 2A). These findings highlight that critical genomic events occurred not only at the MRCA level but also in descendant clones, contributing to malignant transformation.

To increase statistical power and validate our findings using external PDAC data, we performed mutational timing analysis on the accessible PDAC samples from the PCAWG dataset using our in-house developed approach4,13,14. Although the PCAWG cohort is based on single-region sequencing and includes both conventional PDACs and those arising from IPMNs, the distribution of mutations across timing categories, particularly the distinction between clonal early and subclonal mutations, supports our key observation of pervasive clonal heterogeneity across all samples (Figure S3.A). When comparing pooled proportions of clonal early versus subclonal mutations, no statistically significant difference was observed (Figure 3S.B). However, four recurrent driver mutations, most notably KRAS and TP53, were predominantly clonal early, consistent with their role as early trunk events. In contrast, mutations in SMAD4 and LRP1B occurred both early and late, supporting a model of subclonal evolution (Figure 3S.C). These findings align with our current study (Figure 2.A), reinforcing the temporally ordered acquisition of driver mutations in PDAC.

Single vs multiple branch clonal evolution in IPMN
By investigating whether lineages in different patients are associated with specific genomic alterations, we identified two patterns of subclone expansion: (1) a linear progression where the ancestral clone gives rise to a single subclone, and (2) a branching progression where the ancestral clone generates multiple subclones diverging independently.

In the linear model (Cases 4, 12, 13, and 15), a dominant clone acquires successive alterations, generating a sequential series of subclones, each retaining the genetic features of its predecessor while gaining unique mutations. This stepwise progression highlights a hierarchical evolution of disease. In contrast, the branching model (Cases 2, 3, 6, 7, 9, 10, and 16) involves multiple subclones diverging independently from the ancestral clone, resulting in a complex architecture with distinct lineages, each harboring unique genetic alterations in addition to those inherited from the original clone.

While these patterns are distinct, the observed tree shapes may be influenced by random sampling of IPMN regions. Notably, cases with branching evolution exhibited significantly higher counts of intrachromosomal structural variations (median 85.2 vs 29.7, p = 0.0024), suggesting a potential role for SVs in driving the generation of diverging subclones (Figure 2B and 2C).

Temporal dynamics of mutational signatures
To study the mutational processes during clonal progression in IPMNs, we calculated the proportion of mutations assigned to each mutational signature across the branches of the tumor phylogenies (Figure 2D). Overall, SBS1, SBS5 (clock-like/age-related), and SBS40 were the most prevalent signatures, present in both the trunk and branches of the phylogenies. The AID/APOBEC-associated signatures SBS2 and SBS13 emerged in the MRCA of Case 7 but were restricted to the branch associated with invasive PDAC. Similarly, in Case 16, SBS28 was observed only in the PDAC subclone, and not in the pre-malignant IPMN branch.

ID1 and ID2 signatures, linked to DNA mismatch repair deficiency, were observed in the early stages of tumor initiation in most cases (2, 3, 4, 6, 7, 9, 15, 16) but were unexpectedly absent in subclones. A similar trend was noted for other indel signatures (ID5, ID6, and ID9), which appeared in earlier clones but not in subclones. This discrepancy may reflect the smaller number of mutations assigned to subclones compared to trunks and the overall lower prevalence of indels relative to SNVs, potentially limiting the detection of certain signatures due to insufficient power.

Modelling malignant transformation from IPMN to PDAC
In six cases with distinct grades of dysplasia and/or associated cancer, two (Cases 6 and 9) included IPMN with microinvasive PDAC (IPMN_HGD_PDAC), while four (Cases 7, 4, 15, and 16) provided samples of both IPMN and invasive cancer. In Case 2, WGS was performed on both IPMN and PDAC samples, but the analysis was successful only for the IPMN.

In the phylogenetic tree of these cases, the branch point where IPMN and PDAC diverge indicates the emergence of infiltrative clones.

In Case 9, the clonal architecture showed a branching pattern, with the ancestral clone (grey) giving rise to several independent subclones (Figure 2A). The MRCA clone, carrying KRAS, TP53, and LRP1B mutations, as well as the “blue” and “green” subclones were observed in both IPMN_HGD_PDAC and PDAC samples, suggesting their potential to progress to invasive cancer. In contrast, subclones found only in PDAC samples (“yellow, pink, and purple”) likely represent a secondary wave of clonal expansion specific to cancer. Notably, one PDAC sample (S3) had multiple intrachromosomal SVs and contained most subclones, suggesting that these SVs may have played a role in driving subclone development.

Similarly, in Cases 6 and 16, a single MRCA was shared by IPMN_HGD and PDAC samples, suggesting that the MRCA in HGD may have already acquired invasive potential. In Case 6, similar SV patterns were detected in both HGD and PDAC samples. In contrast, in Case 16, PDAC samples (S4 and S5) showed a higher number of intrachromosomal SVs than the HGD sample (S2), potentially driving invasive disease and PDAC-specific subclones (green).

In Cases 4 and 15, distinct ancestral clones were identified for IPMN and invasive cancer samples. Notably, in Case 15, HGD clones carried KRAS, GNAS, and ATM mutations, while invasive cancer clones displayed different mutations in KRAS and TP53 (Figure 2A), indicating separate evolutionary trajectories for high-grade lesions within the same patient.

Transcriptomic patterns in IPMNs
To explore the biological programs associated with IPMN progression, we examined gene expression dynamics across individual cases (Figure S4). Although some pathways showed differential enrichment between invasive and non-invasive samples, only adipogenesis, androgen response, and heme metabolism were significantly more enriched in non-invasive samples, consistent with metabolic reprogramming during progression (Figure S5 A – B).

Principal Component Analysis (PCA) of RNAseq data from samples matched to WGS revealed two main transcriptomic clusters: Cluster 1 (Cases 6, 9, 12, 13) and Cluster 2 (Cases 2, 4, 7, 10, 15, 16) (Figure 1D). Although histological features were comparable between the clusters, Cluster 1 was genomically enriched for LRP1B mutations, exhibited higher TMB (SNVs + indels, p = 0.02), and displayed lower expression of most cancer hallmarks, except for those related to development, bile acid metabolism, and KRAS signaling (downregulated target enrichment) (Figure 3A–C)15.

Principal Component Analysis (PCA) of RNAseq data from samples matched to WGS identified two distinct gene expression clusters: Cluster 1 (Cases 6, 9, 12, 13) and Cluster 2 (Cases 2, 4, 7, 10, 15, 16) (Figure 1D). Cluster 1 showed lower expression of most cancer hallmarks, except those related to development, bile acid metabolism, and downregulated KRAS signaling15, genomically, higher TMB (combining SNV and indels, p = 0.02) and was enriched for LRP1B mutations, with no discernible histological differences (Figure 3 A-C).

Following gene expression normalization, a transcriptomic gradient score was assigned to each sample and compared to reference PDAC cohorts (ICGC PACA-AU16 and TCGA-PAAD cohorts14). Higher scores indicated greater alignment with the squamous molecular subtype (Figure 2A and S6).

IPMNs showed overall lower gradient scores compared to PDACs, indicating a higher prevalence of the Classical subtype, reflective of earlier-stage, progenitor-like disease with better prognosis (Figure 1E). Classical subtype uniformity was observed in non-invasive IPMN cases (2, 10, 12), while significant intratumoral variation was evident in Case 4, where an LGD sample (S1) exhibited a squamous-like signature alongside increased SVs and CNAs compared to other samples from the same tumor (Figure 2A).

Although the two clusters exhibited similar transcriptomic gradient scores overall—with greater variability in Cluster 2— distinct expression patterns were observed. Cluster 1 displayed a transcriptional profile consistent with squamous subtype differentiation, marked by elevated expression of squamous-associated genes (TP63, KRT6A, S100A2) and reduced expression of classical transcription factors (GATA6, HNF4A) (Figure 3A–B). In contrast, Cluster 2 reflected a more classical/progenitor phenotype, characterized by higher expression of KRAS, GNAS, GATA6, and HNF4A, suggestive of a more differentiated transcriptional state and potentially less aggressive biological behavior16-18.

Tumor Microenvironment Analysis
Using EPIC cell proportion estimates19, we assessed the tumor microenvironment (TME) during IPMN progression (Figure 4A). Cluster 1 demonstrated higher CD8+ T cell populations compared to Cluster 2 (0.08 vs 0.01, p = 0.01), despite its higher squamous subtype scores, while Cluster 2 exhibited increased macrophage populations (0.006 vs 0.009, p = 0.04) (Figure 4B). PDAC samples showed significantly elevated CAF populations (median 0.24 vs 0.04, p = 0.04) and reduced epithelial cell populations (median 0.52 vs 0.66, p = 0.01) compared to premalignant IPMNs, indicating a stromal shift associated with invasive disease (Figure 4C, 4D). CD4+ and CD8+ T cell counts were generally lower in PDAC compared to IPMN across cases (7, 9, 15, 16), though these differences did not reach statistical significance.

In LGD_IPMN samples (Cases 10, 12, 15), cellular composition remained consistent, characterized by higher CD4+ T cells and reduced CAFs compared to high-grade lesions. ESTIMATE scores20, which summarize stromal and immune components, confirmed positive stromal scores in most PDAC samples, with the exception of Case 7, where histology revealed prominent inflammatory infiltration. Higher stromal scores were associated with increased CAF populations and squamous subtype enrichment (Figures S7, S8, S9). Correlation analysis revealed that higher Squamous Scores were moderately associated with increased ESTIMATE Scores, strongly associated with stromal scores, suggesting that squamous features in IPMNs are linked to a more stromal-rich tumor microenvironment. This pattern was validated in the ICGC cohort, where the squamous subtype showed elevated stromal and ESTIMATE scores, while immune infiltration remained weakly associated (Figure S8).

TME cell proportions were compared between two distinct evolutionary trajectories: single MRCA vs multiple independent clones pathway. No significant differences were observed in B cells, macrophages, or CD4⁺ T cells between the two trajectories. However, CD8⁺ T cell proportions and ESTIMATE Immune scores were significantly elevated in tumors that evolved through the single MRCA pathway (Figures S10).

Methods
Patient cohorts and multiregional samples collection
Primary tumor samples were collected from patients undergoing surgery at Glasgow Royal Infirmary, UK, and the Pancreas Institute, University of Verona, Italy, following clinical guidelines for IPMN management2,41,42. Fresh-frozen tissue was obtained for genomic analysis using a multiregional sampling protocol, dividing each lesion into 4 to 9 segments, with their positions recorded in a clockwise numeric order for spatial reconstruction. Expert pancreatic pathologists (FD, CL, AS) ensured samples were from the same lesion with no normal tissue in between, to avoid the inclusion of concomitant but separate tumors. All slides from fresh-frozen specimens were centrally reviewed to confirm and classify IPMN or PDAC according to morphology and assess molecular analysis suitability. Following DNA and RNA extraction and quantification from each tumor and matched normal sample, whole genome- (WGS) and RNA- sequencing were performed for comprehensive genomic and transcriptomic analysis. Specific immunohistochemical staining for P53 and SMAD4 was performed on FFPE diagnostic slides for validation of the genomic results. The materials used have been collected under Ethics Committee Approval (ECA) of both institutions (Glasgow: ref. 22/WS/0020; Verona: program 1885 protocol 52438 and amended with protocol 25982).

Read alignment and somatic variant calling
Paired sequence reads were aligned to the human reference genome (hg38). Somatic single nucleotide variants (SNVs) and insertions/deletions (InDels) were identified by comparing normal and tumor pairs, with filtering applied to generate high-confidence variant datasets. Structural variants were processed and categorized into different types for further analysis.

Detection of clonal and subclonal copy number alterations (CNAs)
Battenberg algorithm based-pipeline was employed to identify copy number and estimate tumor purity and ploidy as previously described43 (v2.2.9;https://github.com/Wedge-lab/battenberg). Among the CNAs found in the set of samples, regions altered at a significant high frequency were identified using GISTIC244 with the following modified parameters: -brlen 0.7 and -conf 0.99.

Mutation and copy number signature profiling
Mutation signatures were analysed by using SigProfilerExtractor based on a non-negative matrix factorization (NNMF) framework (v0.0.5.77). Signatures from De novo extraction and decomposition were profiled as single-base substitution (SBS96), double-base substitution (DBS78) and small insertion-deletion (ID83)36. Using ASCAT matrix generated from CNA datasets, copy number signatures were extracted by SigProfilerExtractor45. Based on the updated known signatures as in COSMIC database, SigProflerAssignment was utilised to retrieve decomposed signatures.

In order to minimise NNMF artefact potentials, initial signature extraction was performed with all available samples simultaneously46 and each single signature assignment per sample was finally determined by the probability matrix.

Mutation clustering and phylogenetic tree reconstruction
Subclonal structures were modeled by analyzing mutation allele fractions, copy number variations, and cellularity data. A previously described Bayesian Dirichlet process (DPClustering)47, (v2.2.8; https://github.com/Wedge-lab/dpclust/releases) was used to cluster mutations and reconstruct phylogenetic trees, applying principles to determine linear and branching relationships among subclones.

Annotation of the trees with mutations and signatures
To annotate each tree with oncogenic or putative oncogenic alterations, in-house programs were utilised to connect multiple genomic variants and identify specific significance, including SNV, Indels, SV, CNA, mutation signature and cluster assignment information from mutation caller, Battenberg, SigProfilerExtractor and multisample DPClust. All pipelines applied for WGS data analysis and their steps are illustrated in the flowchart (Figure. S1A2), with detailed descriptions provided in the Supplementary Methods.

RNAseq analysis
Scores were calculated for each sample for the following transcriptomic subtypes: Bailey subtypes16, Collisson subtypes48, and Moffitt subtypes49. Two further subtype scores, ‘Squamous score’ and ‘Classical score’, were calculated for each IPMN sample using a signature derived from RNAseq data from ICGC PACA-AU16. A transcriptomic gradient score was therefore obtained by subtracting the classical score from the squamous score. Samples were also scored for ten previously defined gene programmes associated with Bailey subtypes and for 50 Hallmark gene sets15.

Cell type proportions were estimated using the EPIC method as implemented in the immunedeconv19 package using TPM counts. The ESTIMATE20 method was used to calculate stromal, immune and overall ESTIMATE scores, implemented in the estimate (https://bioinformatics.mdanderson.org/ estimate/rpackage.html) package.

Statistics and reproducibility
All neoplastic and normal samples were sequenced once: whole-genome sequencing was performed for all neoplastic samples, and RNA sequencing for selected cases as described. No statistical method was used to predetermine sample size. No data were excluded from the analyses. The experiments were not randomized, and analyses of genomic and transcriptomic data were conducted independently.

Data availability
Whole-genome sequencing (WGS) and RNA-seq raw data have been deposited in the European Genome-phenome Archive (EGA) under accession number EGAS50000001182 (https://submission.ega-archive.org/studies/EGAS50000001182).
Processed datasets supporting the findings of this study are available at the following repositories:
- https://github.com/Wedge-lab/IPMNPDACApaperArchive/tree/main/IPMNPDAC_WGS/Data
- https://github.com/Wedge-lab/IPMNPDACpaperArchive/tree/main/IPMNPDAC_WGS/Data/sigDPC
- https://github.com/Wedge-lab/IPMNPDACpaperArchive/tree/main/IPMNPDAC_WGS/Data/svDriverCluster
- https://github.com/Wedge-lab/IPMNPDACpaperArchive/tree/main/IPMNPDAC_WGS/Data/cnvDriverCluster

Code availability
The Code and data used for the analysis are available at:
https://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/ [gitfront.io].

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Pea et al conduct multiregion whole genome and transcriptome sequencing analyses on a series of fresh frozen IPMNs (12 patients, 42 tumor samples, including 37 samples used for clonality analyses and 36 for RNA sequencing). The authors identify two main patterns of clonal evolution of IPMNs,one driven by a single ancestral clone (MRCA) with subsequent subclonal parallel or branching evolution, and the other of multiple independent ancestral clones (MRCAs) with parallel evolution. The authors identify that overall IPMNs progress from LGD to HGD and PDAC via stepwise accumulation of genomic alterations including both SNVs and copy number alterations. this study also found evidence of polyclonal evolution lasts into HGD and PDAC, especially in cases with independent MRCAs (although the numbers of cases on which this conclusion is made is very small). Transcriptomic features showed LGD IPMNs as being predominantly classical, while progression to PDAC was associated with more basal features, and changes in the TME (such as greater CAF density).\nOverall, there is a lot of work done here but the issue of novelty and addition of new knowledge to literature is modest at best. There is substantial prior literature on patterns of clonal evolution in IPMNs using multi-region sequencing (including papers published in this very journal). Examples include Fischer at el, Gastroenterology 2019, Noe et al, Nat Comm 2020, Fujikura et al GUT 2021, Omori et al Gastroenterology 2019). Inconsistencies such as higher proportion of RNF43 in HGD lesions (contradicting findings in Fischer at el, Gastroenterology 2019, Noe et al, Nat Comm 2020) and absence of KLF4 hotspot mutations (Fujikura et al GUT 2021) are barely discussed. This work seems more appropriate for a gastroenterology or pancreas focused journal.", "reviewer_2": "The study by Pea et al. explored evolutionary patterns of intraductal papillary mucinous neoplasm (IPMN), a precursor to pancreatic ductal adenocarcinoma (PDAC). By performing WGS and RNAseq of 54 regions from 12 patients, limiting to 37 regions for WGS and 36 for RNA-seq after QC filtering, the authors identify two evolutionary trajectories and transcriptional clusters. The first trajectory was defined by a single MRCA, while the second evolved from several MRCAs. This study also studied changes in tumour microenvironment (TME) during tumorigenic transformation, also highlighting both transcriptional and clonal diversity observed in these patient samples.\nThe manuscript is well-written, providing a characterization of the clonal evolution and transcriptional landscape of IPMN. The conclusions align with previous findings and expand on the evolutionary trajectories, which had not been explored using a multiregional sampling approach before.\nThe main limitation of the study is the lack of clear novelty, and the findings are mainly descriptive, but the study nevertheless provide a useful resource for studies on pre-malignant lesions. A more thorough integration with transcriptomic data and TME analysis will strengthen the appeal of the manuscript.\n- The data is well presented; however, several areas need improvement. The order of the figures does not align with the text, making it difficult for readers to follow. Additionally, some figures, such as Figures 3D, S3, and S5, are not referenced in the manuscript. The authors provide signature analysis, revealing a new copy number variant (CNV) signature, CNV48B, but they do not explore or discuss the aetiology or significance of this signature. When examining the transcriptional analysis, it would be valuable to compare tumours with different evolutionary paths. This comparison should also include components of the TME.\n- It would be interesting to determine whether the TME composition converges or if there are notable differences between the evolutionary trajectories.\n- The clonality analysis description in Methods is too limited - Section “Mutation clustering and phylogenetic tree reconstruction” - and should be expanded with filtering criteria, signature assignment to clusters etc. In general, code to reproduce the main results is missing apart from references to the basic software used in the manuscript\n- Known driver mutations are described and annotated but not integrated into TME, transcriptomic and clonality types. Obviously, this is a small cohort in the number of patients, making statistical inference challenging. Could external PDAC data on clonal vs subclonal mutations be utilized to gain additional power in identifying early trunk mutations vs more subclonal ones?\n- Indel signatures ID1 and ID2 are found in the clonal/trunk population but not in subclones. The authors suggest that this could be due to undersampling of subclonal mutations. It would be worthwhile to investigate if a large fraction of clonal indels are not identified in subclones of the same lineage. If so, this could suggest false positive calls.\n- More SVs are found in samples with multiple branches. Were samples with multiple branches also on average with higher tumour purity?\n- PCA analysis of the transcriptome analysis suggests two mixed clusters with case 2 clustering separately. Did the authors check for batch effects, which are usually the strongest effects, i.e. could some of the PCA separation be explained by confounders?\n- The transcriptomic and methylation analysis is kept separate from the clonal evolution analysis. Did the authors assess, whether evolution types or tree shapes (linear vs branched) matched any of the cellular deconvolution or clusters? Was RNA-seq and DNA for EPIC analysis taken from the sample sample or from adjacent tissue?", "reviewer_3": "This manuscript presents a high-resolution, multi-regional whole-genome and transcriptome sequencing study to elucidate the evolutionary trajectories of IPMN to PDAC. The study provides important insights into intra-tumoral heterogeneity and its role in malignant transformation, employing several lines of computational and genomic analyses. While the findings contribute meaningfully to the understanding of pancreatic neoplasia, several aspects require clarification, expanded methodological detail, and analytical refinements to strengthen the conclusions.\nMajor Concerns:\n1)For clonal phylogenies, the study presents two evolutionary trajectories (single vs. multiple MRCAs), which is conceptually compelling. However, there is no statistical validation or robustness assessment of the inferred phylogenies. The authors utilize DPClust, which is appropriate, but technical details on input parameters, mutation filtering, or confidence scoring are absent.\n2) There’s also no indication of whether alternative methods were tested (e.g. PyClone, or LICHeE) or if the tree topologies were robust to input variation or bootstrapping. Are the phylogenetic trees depicted in the paper robust across different clustering algorithms? A sensitivity analysis comparing alternative phylogenetic reconstruction methods would increase confidence in these findings.\n3) The study identifies SBS1, SBS5, SBS13, and various indel signatures. However, the biological implications of these signatures in the context of IPMN-to-PDAC progression require further discussion.\n4) The authors should expand and improve the relationship between mutational and transcriptomic findings as they currently are very disconnected. For example, SBS2 and SBS13 are associated with APOBEC mutagenesis. Is there any evidence of increased APOBEC expression in transcriptomic data to support this mechanistic link? Similarly, the ID1 and ID2 signatures suggest mismatch repair deficiencies. Were MSI-related genes (e.g., MLH1, MSH2) examined to determine whether any cases exhibited microsatellite instability?\n5) For transcriptomic clustering, the PCA plot in Figure 1D appears to support more than 2 (more likely three) distinct gene expression clusters, rather than two as claimed. In particular, the group of samples represented by yellow points occupies a spatially distinct region of the PC1–PC2 space, separate from both the more dispersed group at the bottom left and the tightly clustered samples on the right. This third potential cluster suggests transcriptional heterogeneity within what the authors define as Cluster 1. A re-evaluation of downstream transcriptomic results based on three clusters is necessary and the authors should discuss whether their claims hold with 3 clusters. Also reinvestigation of these clusters using unsupervised methods (e.g., hierarchical clustering, DBSCAN, or Gaussian mixture modeling) could better capture the underlying biological diversity and avoid oversimplification of the expression landscape.\n6) The biological relevance of the identified clusters is unclear. What pathways or signaling cascades are enriched in each cluster? How do these relate to clonal architecture or tumor grade? A deeper functional annotation is necessary to draw mechanistic insights.\n7) The study mentions a \"transcriptomic gradient score was assigned to each sample and compared to the ICGC PACA-AU20 and TCGA-PAAD cohorts\" but lacks a clear explanation of how it is computed. This should be more explicitly described.\n8) The study identifies novel CNVs and SVs, but their functional impact is not well explored. Have these structural variations been mapped onto known oncogenic pathways related to pancreatic cancer progression? The authors discuss LRP1B mutations and CNAs but does not provide sufficient functional interpretation. Could this gene’s role in WNT signaling or immune modulation be further elaborated?\n9) The authors reference previous work on IPMN evolution but should further contextualize their findings within existing models.How do the observed clonal architectures compare to those reported previously?\n10) The methods used to annotate phylogenetic trees with oncogenic alterations and mutational signatures are insufficiently described. The authors mention use of “in-house programs” to integrate SNVs, indels, SVs, CNVs, and mutational signatures from various tools, but do not provide adequate detail on the algorithms, decision criteria, or confidence thresholds used to assign mutations to specific branches or subclones ( no code, pipeline diagram, or methodological framework is described). This negatively impacts reproducibility and limits the reader's ability to assess the robustness of the phylogenetic annotations. For a high-impact publication in a journal like this, the authors should make the entire analysis code (for all analyses not just this one) available and provide a detailed methodological supplement explaining how different variant types were integrated and validated across the tree structure.\n11) The manuscript lacks details on how mutational signatures were assigned to tree branches. Were these based on clustered mutations, or aggregated per sample? Were exposures filtered for noise or assigned with confidence intervals?\n12) The description of subtype scoring is non-reproducible. It is not clear which gene sets were used or whether batch correction was applied before subtype classification.\n13) The methods describe numerous comparisons but there’s no mention of multiple hypothesis correction, effect size thresholds, or the statistical models used. This weakens the credibility of reported p-values, especially for associations with cluster identity or progression status.\n14) The figures need significant improvement in terms of making them publication ready. They could benefit from additional annotations, particularly in phylogenetic trees to highlight key driver mutations at branch points. Some figures appear disconnected from the rest (Fig 2C) or hard to interpret (Fig 2B) due to lack of more information. For instance, the pie charts lack numeric annotations or mutation counts, making it hard to assess signature burden or dominance. The authors should add mutation counts or a bar/scale representing the number of SNVs/indels at each node. Also for this same figure, it is unclear how signatures were assigned to branches. Are these inferred from the mutation cluster, or are they aggregated across constituent mutations?"}, "author_response": "We thank the reviewer for this thoughtful comment and fully agree that the cited studies—\nparticularly those by Fischer, Noe, Fujikura, and Omori—have substantially advanced our\nunderstanding of IPMN evolution, especially in demonstrating polyclonality and intratumoral\nheterogeneity. These foundational studies primarily utilized multi-regional WES or targeted\nsequencing.\nOur study builds upon this prior work in several key ways:\n1. We are the first to apply multi-regional whole-genome sequencing (WGS) combined with\nmatched transcriptomic profiling (RNA-seq) across the full histological spectrum of IPMN to\nPDAC. This enabled the detection of non-coding alterations, structural variants (SVs), and\nmutational signatures not captured in previous studies. SVs, in particular, were significantly\nenriched in high-grade lesions and branching evolutionary patterns—supporting a role in\nsubclonal divergence.\n2. We identified a novel copy number signature (CNV48B) and characterized SV patterns\nassociated with invasive progression, expanding the mutational landscape of IPMN beyond\nwhat has been previously reported.\n3. We introduced transcriptomic clustering, identifying two distinct transcriptional programs\naligned with classical and squamous-like PDAC subtypes. These clusters were associated with\ndifferences in LRP1B expression, WNT signaling, TMB, and immune contexture, providing\nnew insights into the phenotypic consequences of clonal evolution.\n4. We confirmed the persistence of polyclonal evolution even in high-grade dysplasia and PDAC,\ndemonstrating that genetically and transcriptionally distinct regions may coexist within the\nsame lesion.\nWe believe these findings offer important new insights into the molecular evolution of IPMNs. In the\nrevised Discussion, we explicitly acknowledge prior contributions while clarifying how our use of\nmulti-regional WGS and RNA-seq enhances resolution into clonal dynamics and tumor phenotype.\nInconsistencies such as higher proportion of RNF43 in HGD lesions (contradicting findings in\nFischer at el, Gastroenterology 2019, Noe et al, Nat Comm 2020) and absence of KLF4 hotspot\nmutations (Fujikura et al GUT 2021) are barely discussed.\nWe thank the reviewer for highlighting the importance of known driver genes such as RNF43 and\nKLF4, which have been previously implicated in IPMN.\nRegarding RNF43, we carefully reviewed the findings of Fischer et al., who reported RNF43 mutations\nin 12 of 20 cases: 6 in low-grade dysplasia (LGD), 2 in both LGD and high-grade dysplasia (HGD), and\n4 in HGD (with 2 of those also harboring LGD). Similarly, Noe et al. reported RNF43 mutations in 8 of\n18 cases, including both IPMN and MCN, with a higher proportion observed in HGD lesions—though\nthis trend may partly reflect a higher sampling of HGD and invasive regions in their cohort. Fujikura\net al. identified RNF43 mutations in 10 of 17 samples, with no significant difference in prevalence\nbetween LGD and HGD.\nIn our dataset, RNF43 mutations were present across multiple samples, and importantly, copy\nnumber losses of RNF43 were significantly associated with high-grade lesions. Thus, our study\nprovides complementary evidence that RNF43 inactivation in IPMN may occur not only through\nSNVs/indels, but also via chromosomal loss—a mechanism uniquely captured by WGS.\nAs for KLF4, Fujikura et al. reported hotspot mutations in 12 LGD samples and 7 HGD samples,\nsuggesting that KLF4 alterations may arise early and persist into higher-grade disease. Noe et al., by\ncontrast, did not report KLF4 mutations in their WES dataset. In our cohort, KLF4 mutations were\ndetected in two samples (15_S4, HGD; and 9_S3, PDAC). These findings have now been added to the\nrevised Figure 1. Given our study’s design—focused on sampling distinct regions of macroscopically\nvisible lesions rather than small LGD foci—it is possible that KLF4 mutations in early, microscopic\nprecursors may have been underrepresented. Nonetheless, our results suggest that KLF4 alterations\ncan be maintained during progression, consistent with prior observations.\nTogether, these observations reinforce the added value of WGS in identifying diverse mechanisms\nof driver gene alteration—including SNVs, indels, and CNAs—and provide a more comprehensive\nview of their role across spatial and histological contexts in IPMN evolution. We have added a\ndedicated discussion of RNF43 and KLF4, and clarified differences from prior studies in the revised\nDiscussion section.\nThank you for your valuable comments and insightful suggestions. We have made revisions below\nto address your concerns point by point:\nThe main limitation of the study is the lack of clear novelty, and the findings are mainly descriptive,\nbut the study nevertheless provide a useful resource for studies on pre-malignant lesions. A more\nthorough integration with transcriptomic data and TME analysis will strengthen the appeal of the\nmanuscript.\nWe thank the reviewer for this thoughtful and constructive feedback. While we agree that\ndescriptive studies of early lesions can serve as valuable resources, we respectfully emphasize that\nour findings go beyond description and offer several novel insights into the genomic and\ntranscriptomic dynamics of IPMN progression. To our knowledge, this is the first study to apply\nmulti-regional whole-genome sequencing (WGS) and matched RNA-seq to primary tissue samples\nspanning the full histological spectrum from non-invasive IPMN with varying grades of dysplasia to\ninvasive PDAC. This design allows high-resolution tracking of clonal evolution and\nmicroenvironmental changes, which have not been comprehensively addressed in previous WES- or\npanel-based studies.\nImportantly, our study also identifies early transcriptomic alterations that stratify IPMNs into two\nbiologically distinct clusters with gene expression programs aligned to classical and squamous PDAC\nsubtypes—entities known to differ in clinical behavior, immune infiltration, and therapeutic\nresponse.\nSpecifically, our work presents the following novel contributions:\n1. Comprehensive genomic profiling across spatially distinct regions, revealing the stepwise\naccumulation of SNVs, indels, CNAs, and structural variants (SVs) from non-invasive IPMN to\nPDAC.\n2. Discovery of novel mutational signatures, including a previously unreported copy number\nsignature (CNV48B) associated with clonal evolution and an SV signature (SV4) enriched in\ninvasive lesions.\n3. Frequent and spatially distinct mutations in key driver genes such as TP53, KRAS, LRP1B, and\nGNAS, supporting subclonal diversification during progression.\n4. Recurrent loss of tumor suppressors (CDKN2A, RNF43, U2AF1) via CNAs, emphasizing the\nadded value of WGS for detecting structural mechanisms of inactivation.\n5. Delineation of two major evolutionary trajectories—one involving a single most recent\ncommon ancestor (MRCA), and the other involving multiple independent ancestral clones—\noffering refinement of current models of IPMN progression.\n6. Integration of transcriptomic and immune deconvolution data, revealing associations\nbetween immune cell populations (e.g., CD8⁺ T cells) and evolutionary trajectories\n7. Identification of early transcriptional divergence that mirrors PDAC subtype specification,\nsuggesting biological commitment to distinct molecular programs even in non-invasive\nlesions.\nWe believe this integrative, multi-modal approach provides novel biological insights into early\ntumour evolution and may guide future efforts in early detection, patient stratification, and\ninterception strategies. As suggested, we have also clarified and strengthened the integration of\ntranscriptomic and tumor microenvironment data in the revised manuscript.\nThe data is well presented; however, several areas need improvement. The order of the figures\ndoes not align with the text, making it difficult for readers to follow. Additionally, some figures,\nsuch as Figures 3D, S3, and S5, are not referenced in the manuscript. The authors provide signature\nanalysis, revealing a new copy number (CNV) signature, CNV48B, but they do not explore or\ndiscuss the aetiology or significance of this signature. When examining the transcriptional\nanalysis, it would be valuable to compare tumours with different evolutionary paths. This\ncomparison should also include components of the TME.\nWe thank the reviewer for pointing this out. We have revised the figure numbering and ensured\nthat the order of all main and supplementary figures now aligns with their appearance in the text.\nIn addition, Figures 3D, S3, and S5 are now explicitly referenced in the revised manuscript and\nSupplementary Materials to improve clarity and readability.\nWe appreciate the insightful comment regarding the newly identified copy number signature,\nCNV48B. In the revised manuscript, we have expanded the discussion to explore the potential\nbiological and aetiological relevance of this signature.\nNotably, CNV48B is observed across both low- and high-grade lesions and is characterized by\nmarked loss of heterozygosity (LOH). We found a moderate and statistically significant positive\ncorrelation between CNV48B activity and the activities of ID1 and ID2 (R = 0.53 and 0.55; p = 0.003\nand 0.0002, respectively). As ID1 and ID2 are associated with DNA mismatch repair deficiency in\ncancer, this correlation suggests that CNV48B activity may be linked to defects in DNA repair\nprocesses.\nGiven its distinct pattern and association with LOH, we hypothesize that CNV48B may reflect\nunderlying chromosomal instability. However, we acknowledge that further functional validation in\nlarger cohorts will be essential to confirm its aetiology and clinical significance.\n- It would be interesting to determine whether the TME composition converges or if there are\nnotable differences between the evolutionary trajectories.\nThank you for your comment. We identified two distinct evolutionary trajectories, referred to as\npathway 1 and pathway 2. In pathway 1, tumour clones expanded from a single most recent\ncommon ancestor (MRCA), as observed in cases 3, 6, 7, 9, 10,12, 13 and 16. In contrast, pathway 2\ninvolved progression from multiple independent clones in case 2, 4 and 15.\nTo investigate whether the tumour microenvironment (TME) composition differs between these\ntrajectories, we used the EPIC deconvolution method to estimate immune cell proportions from\nbulk RNA-seq data. We have now integrated this analysis into our genomic evolution framework and\nbriefly discuss it in the revised manuscript.\nB cell, Macrophage and CD4⁺ T cell levels did not show significant differences between tumours\nfollowing pathway 1 and pathway 2. However, CD8⁺ T cell proportions were significantly higher in\ntumours that progressed via pathway 1.\nThis suggests that while helper T cell involvement may be consistent across both evolutionary\ntrajectories, cytotoxic T cell infiltration is enhanced in tumours following pathway 1. These\ndifferences may reflect variations in immunogenicity or tumour-immune interactions linked to the\nunderlying clonal architecture.\n- The clonality analysis description in Methods is too limited - Section “Mutation clustering and\nPhylogenetic tree reconstruction” - and should be expanded with filtering criteria, signature\nassignment to clusters etc. In general, code to reproduce the main results is missing apart from\nreferences to the basic software used in the manuscript\nWe have expanded the description of the clonality analysis in the \"Mutation Clustering and\nPhylogenetic Tree Reconstruction\" section to include detailed filtering criteria. Specifically, we have\nreferenced the exact thresholds and parameters used in the mutation clustering process\n( https://github.com/Wedge-lab/multidimDPClust).\nWe have also elaborated on the process of signature assignment to clusters. This includes the\nspecific mutational signature analysis tools used and how signatures were assigned to different\nmutation clusters based on their mutational profiles. We have referenced the algorithm used for\nsignature extraction and the statistical approaches to assign these signatures to the identified\nclusters.\nWe have summarized the relevant code for the analysis in our lab’s GitHub repository, including\ndetailed instructions on how to access and use it. We have provided links to the repository in the\nrevised manuscript and the supplementary section to ensure full transparency and facilitate\nreproducibility of our results.\nWe have also referred to our previously published methods where applicable to maintain\nconsistency and to give additional context to the methodology employed.\n( https://www.nature.com/articles/ncomms3997)\n- Known driver mutations are described and annotated but not integrated into TME,\ntranscriptomic and clonality types. Obviously, this is a small cohort in the number of patients,\nmaking statistical inference challenging. Could external PDAC data on clonal vs subclonal\nmutations be utilized to gain additional power in identifying early trunk mutations vs more\nsubclonal ones?\nTo address the concern regarding integration of known driver mutations with TME,\ntranscriptomic, and clonality subtypes, we sought to increase statistical power and validate our findings using\nexternal PDAC data. Specifically, we performed mutational timing analysis on the accessible PACA\nsamples from the PCAWG dataset using our in-house developed approach (described in\nSupplementary Methods). This method integrates cancer cell fraction (CCF) estimates with local\ncopy number information derived from Battenberg and DPClust ( https://github.com/Wedge-lab ).\nAlthough the PCAWG cohort is based on single-region sequencing and includes both conventional\nPDACs and those arising from IPMNs, the distribution of mutations across timing categories,\nparticularly the distinction between clonal early and subclonal mutations, supports our key\nobservation of pervasive clonal heterogeneity across all samples (Figure 3S.A). When comparing\npooled proportions of clonal early versus subclonal mutations, no statistically significant difference\nwas observed (Figure S3.B). However, four recurrent driver mutations, most notably KRAS and TP53,\nwere predominantly clonal early, consistent with their role as early trunk events. In contrast,\nmutations in SMAD4 and LRP1B occurred both early and later, supporting a model of subclonal\nevolution. These findings align with our current study (Figure 2A), reinforcing the temporally ordered\nacquisition of driver mutations in PDAC. These results have now been incorporated into the Results\nsection of the revised manuscript.\nUnfortunately, due to the lack of large-scale datasets that include both whole-genome sequencing\n(WGS) and matched transcriptomic or TME profiling in pancreatic cancer, we were unable to directly\nintegrate mutation timing with transcriptomic subtypes or TME features. This remains a key\nlimitation in the field, and future integrative studies with multi-omic data will be essential to fully\nresolve the interplay between genomic evolution, transcriptional programs, and the tumour\nmicroenvironment in PDAC.\n- Indel signatures ID1 and ID2 are found in the clonal/trunk population but not in subclones. The\nauthors suggest that this could be due to undersampling of subclonal mutations. It would be\nworthwhile to investigate if a large fraction of clonal indels are not identified in subclones of the\nsame lineage. If so, this could suggest false positive calls.\nThank you for your thoughtful suggestion. We appreciate your concern regarding the possibility of\nfalse-positive calls for clonal indels.\nBy definition, clonal mutations occur in every tumour cell in a sample, and they therefore occur in\nall subclonal populations of cells. However, we are not able to formally demonstrate that this is the\ncase in this study, in the absence of single cell sequencing. This would be an interesting technique\nfor use in future studies.\nWe have investigated whether indels that occur clonally in one or more samples occur subclonally\nin other samples. Our findings suggest that clonal indels predominantly appear in the trunk\npopulation (i.e. clonal in every sample). The consistency of calls across multiple samples is strongly\nsupportive of the great majority of these indel calls being true positives. On the other hand, some\nsubclones may acquire additional mutations over time. These mutations may be present at lower\nallele frequencies in later subclones, making them more difficult to detect — even with deep\nsequencing approaches. For this reason, we believe that false negative calls of subclonal mutations\nare much more likely than false positive calls of clonal mutations.\n- More SVs are found in samples with multiple branches. Were samples with multiple branches\nalso on average with higher tumour purity?\nNo, samples with multiple branches are not on average associated with higher tumour purity. We\ncomputed the average purities for both groups: the mean purity for single branch samples is 0.332,\nwhile for multiple branch samples, it's 0.316. A Fisher-Pitman permutation test showed no\nsignificant difference between the two groups (p = 0.77), suggesting that tumour purity does not\nsignificantly differ between samples with single versus multiple branches.\n- PCA analysis of the transcriptome analysis suggests two mixed clusters with case 2 clustering\nseparately. Did the authors check for batch effects, which are usually the strongest effects, i.e.\ncould some of the PCA separation be explained by confounders?\nWe appreciate the reviewer’s observation. The separation of Case 2 in the PCA plot is indeed\nnotable. However, all samples in this cluster originate from the same patient, and we interpret this\nas a patient-specific effect rather than a distinct transcriptional subtype. Batch correction was\nperformed during normalization, and we carefully evaluated additional confounders, including\nhistological grade (LGD, HGD, PDAC), epithelial subtype (pancreatobiliary vs. intestinal), and patient-\nlevel variables.\nDespite correction for these factors, some individual cases—including Case 2—continued to form\ndistinct clusters, which we believe reflects genuine biological heterogeneity rather than technical\nartifacts. We have now clarified the normalization and batch correction procedures, as well as our\nassessment of confounding variables, in the revised Supplementary Methods.\n- The transcriptomic and methylation analysis is kept separate from the clonal evolution analysis.\nDid the authors assess, whether evolution types or tree shapes (linear vs branched) matched any\nof the cellular deconvolution or clusters? Was RNA-seq and DNA for EPIC analysis taken from the\nsample sample or from adjacent tissue?\nWe thank the reviewer for this important observation. As specified in the Methods section, RNA and\nDNA were extracted from the same fresh-frozen samples, ensuring direct correspondence between\ntranscriptomic (including EPIC deconvolution) and genomic analyses.\nRegarding the relationship between evolutionary trajectories and the tumour microenvironment,\nthis was addressed in response to a previous comment. In brief, we observed a higher proportion of\nCD8⁺ T cells in cases with a single-MRCA evolutionary trajectory, while other immune populations\nshowed no significant differences across trajectories. We also examined tree shape (linear vs.\nbranched) and found no consistent differences in immune composition between these categories.\nThese findings are now included in the revised manuscript.\nWe acknowledge that robust validation of phylogenetic models is essential for ensuring the\nreliability of our results.\nIn this study, we employed DPClust to infer clonal phylogenies, following the same methodology\noutlined and validated in our previously published work (Noorani et al 2020, PMID 31907488). We\nchose DPClust due to its proven utility in accurately identifying clonal lineages and its ability to\nhandle the complexities of tumour evolution in high-dimensional datasets. We note that, in our\nprevious study, we performed extensive validation using simulated datasets, which demonstrated\nthe robustness of our approach.\nTo address the specific concerns raised, we have now included additional details in the revised\nmanuscript regarding the following:\nInput parameters: We used the same set of parameters as described in our prior work, with values\nchosen to balance sensitivity and specificity in the clonal inference.\nMutation filtering: We applied the same mutation filtering steps as previously published, focusing\non high-confidence somatic mutations.\nConfidence scoring: We have referenced further explanation regarding how the confidence in the\nclonal assignment was assessed.\nThank you for your insightful comment. Our analysis is based on a well-established pipeline that we\nhave developed and refined over the past decade, which has been rigorously validated and widely\nadopted in previous studies published in leading journals. This approach has been cited multiple\ntimes, demonstrating its reliability and robustness in phylogenetic reconstruction.\n(Nat Commum 2014 https://www.nature.com/articles/ncomms3997\nNature 2015 https://www.nature.com/articles/nature14347,\nNat Genet 2020 https://www.nature.com/articles/s41588-019-0551-3)\nGiven that our study focuses on biological insights rather than methodological development, we did\nnot perform an extensive comparison with alternative methods such as PyClone or LICHeE.\nHowever, our method has been benchmarked in prior work and has consistently produced\nbiologically meaningful and reproducible results. Moreover, the tree topologies presented in our\nstudy were carefully evaluated for robustness, ensuring that our conclusions remain valid.\nWhile a full sensitivity analysis across different clustering and phylogenetic reconstruction\nalgorithms is beyond the scope of this study, our confidence in the presented results is supported\nby the extensive validation of our pipeline in prior research.\nWe thank the reviewer for this insightful comment. In the revised manuscript, we have expanded\nthe discussion of the biological relevance of the observed mutational signatures.\nSBS1 and SBS5 were the most prevalent signatures in our dataset and are well-established \"clock-\nlike\" signatures, associated with spontaneous deamination of 5-methylcytosine and endogenous\nreplication-related mutagenesis, respectively. Their consistent presence across both trunk and\nbranch mutations supports the role of age-related mutational processes in the initiation and early\nmaintenance of IPMNs.\nSBS13, and to a lesser extent SBS2—both attributed to APOBEC enzyme activity—were observed\nexclusively in the PDAC subclone of Case 7. Notably, this subclone also exhibited increased structural\nvariation, while the matched HGD sample displayed a distinct SV profile, suggesting independent\nclonal evolution. This observation aligns with findings in other cancer types, where APOBEC activity\nhas been linked to chromosomal instability and subclonal diversification.\nIndel signatures ID1 and ID2, typically associated with mismatch repair deficiency (MMR-D), were\nidentified in early branches of several cases but were largely absent from subclones. This pattern\nmay reflect limited statistical power to detect indel signatures in low-mutation subclonal branches.\nSupporting a non-MSI phenotype, all cases retained expression of MMR proteins (MLH1, PMS2,\nMSH2, and MSH6) by immunohistochemistry, consistent with the absence of microsatellite\ninstability and the overall low tumor mutational burden observed across our cohort.\nAltogether, our findings support a model in which early IPMN genomes are shaped predominantly\nby age- and replication-related processes, while APOBEC-mediated mutagenesis and structural\ninstability emerge later, contributing to malignant transformation. These interpretations are now\ndiscussed in greater detail in the revised manuscript.\nWe thank the reviewer for this important point. Regarding mismatch repair (MMR) deficiency, we\nperformed immunohistochemistry for MLH1, PMS2, MSH2, and MSH6 on all available samples. All\ncases showed retained protein expression, indicating intact MMR function. This finding is consistent\nwith the lack of microsatellite instability on molecular analysis and the overall low tumor mutational\nburden observed across the cohort. These results have now been clarified in the revised Methods\nsection.\nAs for APOBEC-associated signatures (SBS2 and SBS13), these were identified in a limited number of\nsubclones (notably in Case 7). Given the small number of APOBEC-positive cases, and the limited\nstatistical power to detect differential expression, we did not perform formal transcriptomic\ncorrelation with APOBEC family gene expression.\nWe thank the reviewer for this thoughtful comment. We acknowledge that in the PCA plot, the\nyellow points appear spatially distinct. However, these samples all originate from a single patient\n(Case 2), and their separation is likely driven by patient-specific transcriptional variation.\nWe also performed unsupervised hierarchical clustering based on Euclidean distances using VST-\ntransformed gene expression data. As shown in the correlation heatmap and dendrogram, the\nsamples consistently segregated into two main clusters, in agreement with the PCA results. Three\nsamples with a high proportion of zero counts were excluded from the analysis.\nWe additionally examined the distribution of clinical and histological features—including grade,\nIPMN subtype, and sampling site—and found that these variables did not account for the observed\nclustering. These clarifications have been added to the revised Results and Supplementary Materials.\nWe thank the reviewer for raising this important point. Subtyping pancreatic neoplasms into\nsquamous and classical-like categories is clinically relevant due to well-established differences in\nprognosis, immune microenvironment, and therapeutic response. In our study, we identified two\ndistinct transcriptomic clusters. To clarify their biological significance, we expanded the\ntranscriptomic analysis and updated the Discussion accordingly.\nCluster 1 displayed higher expression of squamous-associated genes (e.g. TP63, KRT6A, S100A2),\nelevated LRP1B expression, increased TMB, and greater immune infiltration—features consistent\nwith a more immunogenic and genomically unstable phenotype. In contrast, Cluster 2 was marked\nby higher expression of classical transcription factors (GATA6, HNF4A), consistent with a more\nprogenitor-like and differentiated state. While most samples in our cohort exhibited lower\ntranscriptomic gradient scores compared to conventional PDACs in the ICGC and TCGA datasets—\nconsistent with an earlier-stage, classical phenotype typical of IPMN-derived pancreatic cancers—\nthe divergence between Clusters 1 and 2 highlights potential early commitment to distinct molecular\nprograms.\nThese changes are now reflected in the revised Results and Discussion sections. We believe this\nintegrated interpretation enhances the biological and translational relevance of the transcriptomic\nclustering analysis\nWe thank the reviewer for this helpful comment. In the revised Supplementary Methods, we have\nclarified the computation of the transcriptomic gradient score. Specifically, single-sample gene set\nenrichment analysis (ssGSEA) was applied to each RNA-seq sample using curated gene signatures\nfor the classical and squamous PDAC subtypes, as previously published by our group (Bailey et al.\n(Nature 2016) and Brunton et al. (Cell Reports 2020)). The final gradient score was calculated as the\ndifference between the enrichment scores for the squamous and classical signatures. This approach\nallows for a continuous estimate of subtype identity across the epithelial transcriptomic spectrum.\nThis information was added in the supplementary method section.\nWe thank the reviewer for this insightful comment. Our study was designed to address the\nunderexplored contribution of structural variations (SVs) and copy number alterations (CNAs) in\nIPMN progression using multi-regional WGS—a level of resolution not available in prior IPMN studies\nbased on WES or targeted sequencing.\nWe observed recurrent CDKN2A deletions, consistent with its known involvement in PDAC\nprogression, and frequent RNF43 losses associated with high-grade lesions. While RNF43 mutations\nhave been linked to dysregulated WNT signaling in IPMN, our identification of recurrent copy\nnumber loss extends these findings and supports the notion that biallelic inactivation may contribute\nto progression. U2AF1 deletions were also enriched in high-grade IPMNs, although their role in\npancreatic tumorigenesis remains to be clarified.\nRegarding LRP1B, this gene was both mutated and transcriptionally downregulated in Cluster 1—\nalongside higher tumor mutational burden, reduced WNT signaling activity, and increased CD8+ T\ncell infiltration. These findings are consistent with the immunogenic phenotype associated with\nLRP1B loss in other malignancies, and support a potential link between LRP1B inactivation, immune\nevasion, and transcriptional reprogramming in a subset of IPMNs. However, our study did not\nexperimentally validate the functional consequences of these alterations, and we agree this\nwarrants further investigation.\nNotably, although we did not identify specific SVs consistently enriched across high-grade IPMNs or\nPDACs, we observed a significantly higher SV burden in tumors with branching clonal evolution (as\nshown in Figure 2D). This suggests that SVs may play a key role in promoting subclonal diversification\nand genomic instability. This observation challenges prior models of IPMN progression largely based\non SNV/indel data from WES, and underscores the added insight provided by WGS in revealing the\nstructural complexity of these lesions.\nWe have therefore reinforced the potential role of SV in the discussion.\nWe thank the reviewer for this important comment. As outlined in prior responses and now\nexpanded in the revised Discussion, we have acknowledged the contributions of previous studies\nand clarified how our work builds upon and extends this body of literature.\nWhile earlier models based on targeted sequencing and WES demonstrated the polyclonal origin of\nIPMNs—primarily through the identification of distinct KRAS mutations within the same IPMN\nlesion—they largely inferred clonal heterogeneity from SNVs. These models proposed that early\nsubclonal diversification drives progression through clonal selection during the transition to high-\ngrade lesions, with substantial mutational heterogeneity observed only in LGD regions.\nOur study builds upon and evolves these models by applying multi-regional WGS, enabling unbiased\ndetection of structural variants SVs, copy number alterations, and non-coding mutations—classes of\nalterations not captured by WES. This approach allowed us not only to confirm polyclonality but also\nto demonstrate the complete absence of shared CNAs and SVs among genomically distinct clones\nwithin the same lesion, reinforcing a model of parallel evolution. Furthermore, we show that\npolyclonal architecture can persist into high-grade dysplasia and associated PDAC, suggesting that\nclonal divergence may extend further into tumorigenesis than previously appreciated.\nWe also note that prior studies employed Bayesian frameworks, such as the JAGS timing model\n(PMC7428044), to estimate the temporal dynamics of malignant progression. However, such models\nrely exclusively on SNVs and do not account for SVs or CNAs. Given the central role that SVs appear\nto play in IPMN evolution, as highlighted in our findings, we chose not to apply these timing models,\nas they may underestimate the complexity of clonal dynamics in the absence of whole-genome data.\nWe applied both published and in-house-developed pipelines for phylogenetic annotation. To\nenhance clarity and reproducibility, we have added a flowchart (Figure S1A2) detailing each step in\nthe Supplementary Methods section. Additionally, all analysis codes, including those for integrating\nSNVs, indels, SVs, CNAs, and mutational signatures, have been made available on\nhttps://github.com/Wedge-lab/multidimDPClust,\nhttps://gitfront.io/r/xtgitfhe2/eMoLG6vndJxF/IPMNpaperArchive/.\nWe assigned mutational signatures to each branch in the phylogenetic tree using the following steps:\n1. Mutation identification: we obtained the genomic coordinates of all somatic mutations from\nwhole-genome sequencing (WGS) data.\n2. Variant clustering: Using multidimDPClust, we performed multidimensional clustering of variants,\nincorporating confidence intervals for assignment. This process generated mutational clusters along\nwith their corresponding genomic coordinates. We included clusters containing at least 1% of the\ntotal mutations.\n3. Phylogenetic reconstruction: we reconstructed the phylogenetic tree based on the outputs from\nmultidimDPClust, assigning mutations (along with their genomic coordinates) to specific branches.\n4. Signature extraction: we identified single-nucleotide variant (SNV) signatures (SBS96) and\ninsertion-deletion (Indel) signatures (ID83) using COSMIC methodologies. The genomic coordinates\nof mutations contributing to each mutational signature were extracted.\n5. Assignment of signatures to tree branches: The probability of each mutational signature\ncorresponding to a clonal cluster (tree branch) was derived by merging mutation coordinates from\nsteps 3 and 4 for each multi-sample case. For each cluster, we used the probability of each mutation\ncontext being assigned to a signature and summed these probabilities across all mutations within\nthe cluster.\nAll methods are detailed in the Methods section and Supplementary section - Annotation of the\ntrees with mutations and signatures\n(21-Assignment of Signatures to Clustered Branches:\nhttps://gitfront.io/r/xtgitfhe2/eMoLG6vndJxF/IPMNpaperArchive/).\nWe thank the reviewer for this observation. In the revised Supplementary Methods section, we have\nclarified the subtype scoring methodology and batch correction procedures. Specifically,\ntranscriptomic gradient scores were calculated for each sample using single-sample gene set\nenrichment analysis (ssGSEA) based on previously defined classical and squamous gene signatures\n(Bailey et al., Nature 2016; Brunton et al., Cell Reports 2020). The final gradient score was computed\nas the difference between squamous and classical enrichment scores, with higher values indicating\ngreater alignment with a squamous-like transcriptional program.\nBatch correction was applied during normalization using the DESeq2 workflow, and potential\nconfounders such as histological grade, epithelial subtype, and sequencing batch were considered\nin downstream analyses. This information has been added to the revised Supplementary Methods.\nThank you for highlighting the need for greater clarity in our statistical methods. In the revised\nmanuscript, we now include details on the statistical models used for each analysis in the\nSupplementary method section: Data Processing and Statistical Analysis. Multiple hypothesis testing\nwas controlled using the Benjamini–Hochberg procedure, with FDR < 0.05 considered statistically\nsignificant. Specifically, we applied Fisher–Pitman permutation tests to compare the number of\nSNVs, indels, and structural variants (SVs) across tumor types and genomic regions; Fisher’s exact\ntests to compare the frequency of driver mutations (including SNVs, indels, SVs, and CNAs) and SV\nsignatures across tumor types; and permutation-based Wilcoxon tests for pairwise comparisons of\ntiming group proportions\nThank you for these helpful suggestions to improve the clarity and interpretability of the figures.\nWe have revised all figures to enhance their presentation quality and ensure they are publication\nready. In the phylogenetic trees, we have added detailed annotations at key branch points to\nhighlight the presence of known or putative driver mutations. This should help guide interpretation\nof evolutionary trajectories.\nFigure 2B and 2C clarifications: To address concerns about Figure 2B, we have added numeric\nannotations and mutation counts within the pie charts to clearly reflect the signature burden and\ndominance.\nFor Figure 2C, we have clarified its relevance and integration with the rest of the analysis in the\nrevised figure legend and main text, ensuring it is well-connected to the study’s narrative.\nSignature assignment clarification: Regarding the assignment of mutational signatures to branches,\nwe have clarified this in both the figure legend and Methods section. Signatures were inferred based\non aggregated mutational profiles from each cluster of mutations assigned to the corresponding\nbranch.\nWe believe these updates address the reviewer’s concerns and significantly improve the clarity and\nimpact of the figures."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed my queries including the concerns about novelty and how this study differs from prior publications in this area.", "reviewer_2": "The authors have addressed most of my concerns.\nConsistent points:\n- Some of the figures are still low resolution, part of the text is not visible (with different fonts and sizes), and overall they look \"rushed\" - which is odd and unusual to observe in a high-impact journal.\n- The pathway analysis is rather superficial - Fig S4/S5 - is there anything significantly different or conclusions to be drawn from this - instead of \"notable differences\"?\n- I could not access the code in the notebooks on the github (linked through gitfront), e.g.\nhttps://github.com/xtgithubhe/IPMNpaperArchive/blob/main/ipmnPaperCode.ipynb", "reviewer_3": "I appreciate the authors’ revisions and thoughtful responses to the prior round of review. The addition of new analyses, expanded methodological descriptions, and improved code sharing through repositories represent meaningful progress.\nHowever, after detailed examination of the revised manuscript, code repositories, supplementary materials, several important issues remain. These concerns relate to both scientific interpretation and reproducibility, and should be addressed to meet the expectations of a high-impact journal.\n-Reproducibility and code sharing:\nWhile the authors have shared parts of their pipeline through different repos, the full analytical workflow is not yet reproducible. Critical components remain unavailable:\na) Tree Annotation (SNVs, SVs, CNVs, Signatures): No public scripts are provided for mapping mutations and mutational signatures onto tree branches. The process is described narratively but is not executable by readers.\nb) Transcriptomic and Immune Integration: The analyses of transcriptomic clustering, subtype assignment, and immune deconvolution are described but not accompanied by code.\nc) Shared code through different repos reference local paths with no standardized data input structure or templates provided for external reproducibility.\nOverall, while parts of the pipeline are transparent, the complete analytical workflow is not yet reproducible. For a study of this scope, particularly one centered on complex evolutionary modeling, end-to-end reproducibility is essential and everything should be unified in a single repo for users instead of pointing readers to different repos for different parts.\n-Integration of Genomic, Transcriptomic, and Immune Data:\nAlthough new immune analyses were added, the study still lacks mechanistic integration between genomic evolution, transcriptomic subtypes, and immune microenvironment features. The relationships remain correlative. For example, there is no test of whether specific mutational patterns predict transcriptional or immune shifts.\n-Figures are partially improved but some, especially Fig 2, lacks critical quantitative features and full integration of genomic events:\na) Fig 2a: No branch length scaling is provided; trees are topological, not quantitative.\nb) Fig2b: The SV counts across branching architectures are presented, but no mapping of specific SVs onto the phylogenetic trees is performed.\nc) Fig2c: the circos plots are of extremely low resolution, and labels are unreadable. Furthermore, SVs are displayed per sample but are not linked back to the clonal architecture in Fig 2D.\nd)Fig2d: Branch lengths are absent. For a study claiming to map clonal trajectories and mutational processes over time, branch lengths are essential. Also, SVs and CNVs are not mapped onto the trees, despite their importance to the clonal evolution narrative."}, "author_response": "We thank the reviewer for their positive feedback and are grateful that our revisions\nsatisfactorily addressed their concerns.\nWe thank the reviewer for this observation. All figure files have been carefully checked for\nconsistency in format and layout to ensure clarity and reproducibility, and have now been re-\nexported and uploaded at high resolution (300 dpi). We trust that the revised figures now meet\nthe standards of the journal.\nWe thank the reviewer for this valuable observation. While some pathways appeared more\nenriched in invasive compared to non-invasive samples, only adipogenesis, androgen response,\nand heme metabolism were significantly enriched in non-invasive samples. We have revised the\ntext accordingly and replaced Figure S5B with box plots highlighting these significantly enriched\npathways.\nThank you for your comment regarding GitHub accessibility, and we apologize for any\ninconvenience caused by the link being provided without explanation. We have updated the\nGitFront link. All files in the repository, including ipmnPaperCode.ipynb, can be directly viewed\nonline or cloned using the GitFront link provided in the manuscript:\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/\nor\ngit clone https://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/\nThey contain exact duplicates of the previous private GitHub repository\n(https://github.com/xtgithubhe/IPMNpaperArchive/blob/main/ipmnPaperCode.ipynb).\nOnce the manuscript is accepted, the GitHub xtgithubhe repository will be made public.\nThank you for this useful comment. We have now updated the methods and provided the\nmethod and executable code for tree annotation (SNVs, SVs, CNVs, and mutational signatures).\nTo ensure reproducibility, we have shared both the running code and relevant example datasets\nin the same repositories:\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/\nAs outlined step by step:\n1. Subclonal clustering: We applied multidimDPClust to the IPMN–PDAC WGS data, which is\nspecifically designed for multi-region WGS samples (https://github.com/Wedge-\nlab/multidimDPClust).\n2. Phylogenetic tree inference: Based on the identified SNV/indel clusters, we determined the\nmost likely tree branches using the “pigeon-hole” principle (PHP), the “sum rule,” and the\n“crossing rule,” as detailed in references 10–12. We have highlighted this procedure in the\nSupplementary Methods (section: Mutation clustering and phylogenetic tree reconstruction).\n3. Tree annotation of SV: SV clusters were identified by matching their genomic coordinates with\nthose of SNV/indel clusters from the same sample. SV drivers were then assigned to the\ncorresponding branches (see revised section Tree annotation of SV and CNV of the\nSupplementary Methods and analysis code at:\n24-Assignment of SV (gene-driver) to Clustered Branches,\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/ ).\n4. Tree annotation of CNV: Similarly, CNV clusters were assigned to tree branches by genomic\ncoordinate matching to SNV/indel clusters (see revised section Tree annotation of SV and CNV\nof the Supplementary Methods and analysis code at:\n25-Assignment of CNV (gene-driver) to Clustered Branches,\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/ ).\n5. Mutational signature annotation: In addition to the description in the Annotation of the trees\nwith mutations and signatures section of the Supplementary Methods, we have now updated\nand shared the analysis code for assigning signatures to tree branches:\n(21-Assignment of Signatures to Clustered Branches:\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/).\nTogether, these updates provide readers with a full method to reproduce the tree annotation\nsteps.\nWe thank the reviewer for this important comment. The full script for the transcriptomic and\nimmune analyses has now been uploaded at\nhttps://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/ and includes data input,\nclustering, ComBat-Seq correction, subtype assignment, cancer hallmark analysis, and immune\ndeconvolution (EPIC).\nWe are grateful to the reviewer for highlighting the importance of reproducibility. We have\nunified all analysis code into a single GitHub repository:\n(https://gitfront.io/r/xtgitfhe2/8f2rTdyZs1S9/IPMNPDACpaperArchive/)\nThis repository contains a standardized input directory along with relevant datasets required to\nreproduce the results. It also includes wrapper scripts covering the entire method from data\nprocessing and statistical analysis to evolutionary modeling, mutation assignment to\nphylogenetic branches, and figure generation.\nFurthermore, we have updated the integration analysis of genomic, transcriptomic, and\nimmune data within the same repository. This includes dataset preparation, clustering, ComBat-\nSeq–corrected expression data, subtype assignment, cancer hallmark annotation, and immune\ndeconvolution (EPIC).\nWe believe these updates address the reviewer’s concerns and ensure end-to-end\nreproducibility of our study, while also providing a more integrated view of genomic evolution,\ntranscriptomic subtypes, and the immune microenvironment.\nWe thank the reviewer for this constructive comment. In the revised version, we strengthened\nthe integration between genomic evolution and immune contexture by explicitly comparing\nimmune cell composition across clonal trajectories. We found that CD8⁺ T-cell infiltration and\nESTIMATE Immune Scores were significantly higher in tumors evolving from a single MRCA\ncompared to those with multiple independent clones, suggesting that clonal architecture\nmodulates immune surveillance.\nIn the Discussion, we now interpret this finding by proposing that single-MRCA tumors retain\nuniform clonal alterations that may favor immune recognition, whereas multiple independent\nclones generate distinct, lineage-specific alterations that fragment and dilute antigenic visibility,\nthereby contributing to immune evasion. Although we did not directly assess neoantigen\npresentation, this interpretation is consistent with prior studies showing that tumors enriched\nfor clonal rather than heterogeneous neoantigens are associated with stronger immune\nresponses. We acknowledge that this remains a hypothesis, which we plan to address directly in\nfuture studies through integrated neoantigen prediction and functional validation.\nThank you for your valuable comments. To improve the figures, we have made both scientific\nand technical revisions by adding quantitative features and optimizing visibility.\na) Figure 2a: We have now provided branch length scaling. Specifically, the legend at the bottom\nindicates the length scale based on the total number of SNVs and indels (with units of ~150,\n200, and 1000 mutations). Additionally, for clarity, we have labeled the number of mutations on\nthe longest branch in each tree, which will help readers estimate the relative lengths of all other\nbranches.\nb) Thank you for this suggestion. We have now mapped cancer driver SVs onto the\ncorresponding branches of the phylogenetic trees in Figures 2A and 2D, thereby directly linking\nthe SVs to the clonal architecture.\nc) We have improved the quality of the circos plots in Figure 2C by substantially increasing their\nresolution and adjusting the labeling for better readability. Furthermore, to integrate these data\nwith the clonal evolution framework, we have mapped the cancer driver SVs, including those\ndisplayed in the circos plots, onto the phylogenetic trees shown in Figure 2D.\nd) We appreciate your comment and fully agree on the importance of branch lengths, CNVs, and\nSVs in illustrating clonal evolution. Fig. 2d is intended as an unscaled schematic focusing on\nmutational signature profiles within the clonal trajectories and is largely derived from Fig. 2b.\nIn Fig. 2b, we have already included scaled branch lengths and mapped CNVs and SVs. For\nimproved readability in Fig. 2d, we have now labelled the key driver events related to CNVs and\nSVs on each branch and indicated in the legend that branch lengths can be referenced from Fig.\n2b. We hope this modification addresses your concern."}, {"round": 3, "reviewer_comments": {"reviewer_2": "I have no further comments", "reviewer_3": "The authors have made substantial progress in addressing previously raised concerns. The paper is now closer to meeting the standards of reproducibility and clarity expected for publication in a high-impact journal.\nHowever, I strongly recommend the following minor revisions prior to acceptance:\n1. Provide a README.md or tutorial describing end-to-end use of the shared repository, including required data formats, example commands, and environment setup.\n2. Consider a supplementary figure or model linking mutational features to immune phenotypes in a mechanistic way.\nAssuming these minor improvements are made, I am supportive of publication and commend the authors for their thoughtful revision."}, "author_response": "We thank the reviewer for this helpful suggestion. We have added a comprehensive\nREADME.md file to the shared repository. This document includes detailed instructions for\nenvironment setup, input data formats, and example commands to enable full end-to-end\nreproducibility of the analyses presented in the manuscript.\nWe appreciate this suggestion. Figure 2A now displays all genomic alterations—including\nSNVs/indels and structural variants—mapped onto the phylogenetic branches for each case, and\nthe case order has been aligned with Figure 3A, which defines the transcriptomic clusters. This\nalignment ensures consistency between figures and text, and allows readers to visually integrate\nclonal architecture, genomic alterations, and the corresponding transcriptomic and immune\nphenotypes within a coherent narrative framework."}]

Decision letter