Reconstructing single-cell karyotype alterations in colorectal cancer identifies punctuated and gradual diversification patterns

Central to tumor evolution is the generation of genetic diversity. However, the extent and patterns by which de novo karyotype alterations emerge and propagate within human tumors are not well understood, especially at single-cell resolution. Here, we present 3D Live-Seq—a protocol that integrates live-cell imaging of tumor organoid outgrowth and whole-genome sequencing of each imaged cell to reconstruct evolving tumor cell karyotypes across consecutive cell generations. Using patient-derived colorectal cancer organoids and fresh tumor biopsies, we demonstrate that karyotype alterations of varying complexity are prevalent and can arise within a few cell generations. Sub-chromosomal acentric fragments were prone to replication and collective missegregation across consecutive cell divisions. In contrast, gross genome-wide karyotype alterations were generated in a single erroneous cell division, providing support that aneuploid tumor genomes can evolve via punctuated evolution. Mapping the temporal dynamics and patterns of karyotype diversification in cancer enables reconstructions of evolutionary paths to malignant fitness.

A neuploidy-defined as a cell having an abnormal number of chromosomes or sub-chromosomal fragments-is among the most common features of human cancers 1 . The prevalence of aneuploidy and the detection of recurrent copy-number alterations (CNAs) across cancer types indicate that genome-wide aneuploidy plays an active role in promoting malignant phenotypes 2,3 .
How an aneuploid tumor genome evolves over time remains incompletely resolved. Although single-cell sequencing from patient samples captures genetic diversity in human tumors 4,5 , the evolutionary history of an aneuploid genome is difficult to determine, as the two critical parameters that fuel its progression-ongoing chromosomal instability (CIN) and selection pressures-are dynamically intertwined 3,6-8 . Consequently, general models of tumor evolution are still under debate, in part because the karyotype alteration rate and the patterns by which de novo karyotypes emerge and propagate within human tumors are not well understood 6,8,9 . These parameters in particular are difficult to extract from patient material since it provides limited to no information on genetic intermediates and the timescales that separate individual tumor subclones (the number of cell generations). For instance, punctuated bursts of karyotype alterations in human cancers have been described [10][11][12][13] , but the exact timescale and prevalence of such events remain elusive 14 .
Patient-derived tumor organoids (PDTOs) offer a unique opportunity to address these fundamental questions. PDTOs are three-dimensional (3D) mini-organs derived from primary tumor tissue. PDTOs maintain the histopathological features of the native tumor, including high concordance in somatic mutations and transcriptome and drug response profiles [15][16][17] . Crucially, PDTO models are the closest representatives of human tumors that are compatible with high-spatiotemporal-resolution imaging of cells, enabling live microscopy of clonal tumor outgrowth. stable) and PDTO-19b (microsatellite unstable) 15 . Sequencing data from both PDTO lines displayed de novo whole-chromosome as well as sub-chromosomal CNAs among cells (Fig. 1b,c and Extended Data Figs. 1 and 2), in agreement with our previous study showing that both PDTOs are in CIN 7 . Capturing the entire population of cells from individual clonal PDTOs has the critical advantage of detecting reciprocal copy-number gains and losses between cells, indicating that the two lineages are progeny of the same ancestral CNA event (Fig. 1b,c, red boxes; >1 cell with the same CNA indicates propagation). Across datasets, CNA events occurred at a rate of approximately one in ten cell divisions for both PDTO-9 and PDTO-19b (Fig. 1d). Notably, we identified instances where identical CNAs emerged from independent events within a single PDTO-9 organoid, which were discriminated by the co-occurrence of additional CNAs that prohibited the assembly of a coherent single phylogeny ( Fig. 1b (lineages I and II) and Extended Data Fig. 1b (lineages I and II)). Indeed, we frequently detected cell divisions where multiple CNAs arose simultaneously. Extreme cases were represented by a subset of PDTO-19b cells with highly aberrant karyotypes that could be classified as hopeful monsters 18,19 , defined as a genome-wide set of changes in ploidy that are likely to substantially alter fitness. Two of these cells displayed extensive reciprocal CNAs across their genomes (Fig. 1c, lineage III), indicating that they were the progeny of a single catastrophic cell division. The sum of both karyotypes pointed to a genome-duplicated ancestor, in line with models suggesting that extreme aneuploid states often evolve from unstable tetraploid intermediates 20 (Extended Data Fig. 3a).
Collectively, these data indicate that de novo CNAs readily arise and propagate during PDTO outgrowth. Furthermore, we detected several instances of punctuated karyotype alterations that involved multiple chromosomes in a single-cell division, with hopeful monsters representing the most extreme cases.
Integrating live-cell imaging and single-cell sequencing data using 3D Live-Seq. To further investigate the dynamics of de novo CNAs as they arise and propagate across cell generations, we developed 3D Live-Seq-a protocol inspired by the LookSeq strategy 21 . 3D Live-Seq incorporates confocal live-cell imaging data of a growing PDTO structure and single-cell sequencing data of the entire imaged cell population, to allow a precise reconstruction of evolving tumor karyotypes across consecutive cell generations. The expression of transgenic H2B-Dendra2 (ref. 22 ) is sufficiently bright to support long-term confocal imaging with remarkable spatiotemporal resolution. Importantly, cells of interest can be photoconverted  1  2  3  4  5  6  7  8  9  10  11  12  14  13  15  17  16 18 19  . The dataset size is indicated within each circle. Events include reciprocal CNAs, non-reciprocal whole-chromosome gains or losses and non-reciprocal sub-chromosomal CNAs represented in more than one cell. Hopeful monster karyotypes were excluded from this analysis. CNAs that occurred in the same cell division were considered as one event.
from green to red Dendra2 fluorescence to serve as a reference landmark for the direct integration of imaging and sequencing datasets.
Chromatin errors during mitosis are phenotypes of de novo CNAs. We first implemented 3D Live-Seq to investigate whether de novo CNAs correlate with specific chromatin errors during mitosis. Collectively, we generated high-resolution imaging data of 64 cell divisions, each from individual tumor organoids (PDTO-9), to score their chromatin phenotype. We then sequenced both daughter cells along with a bulk sample consisting of the remaining cells of the dissociated PDTO structure (Fig. 2a). As expected, de novo CNAs were almost exclusively generated among daughter cells of cell divisions that involved either a chromatin bridge or lagging chromatin ( Fig. 2b and Supplementary Note 1). Furthermore, chromatin bridges generally resulted in sub-chromosomal CNAs, whereas whole-chromosome missegregations were mainly associated with lagging chromatin. This bias was consistent with reports showing that chromatin bridges primarily represent the stretching of fused or dicentric chromosomes 23 , likely to result in sub-chromosomal fragments when broken. To exclude that the expression of transgenic H2B substantially exacerbates CIN phenotypes, we generated a CRISPR knock-in 24 of Dendra2 at the carboxy terminus of the HIST1H2BC gene (Extended Data Fig. 4a,b). In a direct comparison, both lines showed clear CIN phenotypes with similar error-type frequency distributions (Extended Data Fig. 4c).
The lower overall CIN rate in the knock-in line could be attributed to the substantially reduced brightness of the knock-in, which hampers effective scoring of subtle chromatin errors, as well as to clonal differences between lines 7 (Extended Data Fig. 4d). Furthermore, single cells isolated from genetically unmodified PDTOs (using the DNA dye Syto 11) that had not been exposed to live-cell imaging procedures, displayed a similar number of de novo CNAs compared with other PDTO datasets (Extended Data Figs. 1d and 2b,c). While these data confirm that chromatin errors are distinct phenotypes of karyotype diversification, chromatin bridges in particular often failed to generate detectable CNAs (Fig. 2b). In contrast with artificially induced chromatin bridges 25 , naturally occurring chromatin bridges in PDTOs are likely to reflect a more complete spectrum of underlying causalities, of which a subset may be resolved at high fidelity, while other cases may give rise to mutations that are below our detection threshold (Supplementary Note 2). 3D Live-Seq captures karyotype alterations across cell generations. Next, we applied 3D Live-Seq to study evolving tumor genomes across multiple consecutive cell divisions. To showcase the feasibility of our protocol, we recorded the unperturbed outgrowth of an individual PDTO-9 organoid across three cell generations (from two to 13 cells) and reconstructed the true mitotic tree with a detailed characterization of mitotic fidelity for each cell division (Fig. 2c). Before single-cell isolation, we selected one cell of interest for photoconversion based on a lagging chromatin phenotype during its preceding cell division (Fig. 2c, cell number 1.1.1.1). Single-cell sequencing data from 12 out of 13 cells displayed several de novo CNAs across three lineages (Fig. 2c, lineages I-III). Using the sequencing result of the photoconverted cell as a landmark and the mitotic tree as a structural constraint on potential phylogenetic solutions of lineage I, we mapped a consecutive missegregation of chromosome 7 to the highlighted branch within the mitotic tree (Fig. 2c, branch 1.1). In agreement with our previous data, divisions at which chromosome 7 missegregated displayed a lagging chromatin phenotype migrating toward the daughter cell that gained a copy of chromosome 7 (Fig. 2c, red arrows). The remaining de novo CNAs could not be mapped with high accuracy as their potential phylogenetic solutions matched multiple branches of the mitotic tree (Fig. 2c, lineages II and III). Lineage III displayed an unusual CNA of chromosome 6q (chromosome 6q24.2-ter) with one cell carrying five copies of chromosome 6q while its sister cell, identified by the reciprocal CNA of chromosome 1q, displayed only one copy of chromosome 6q. Given the structural constraints set forth by the mitotic tree, any conventional phylogenetic interpretation that could explain the origin of additional chromosome 6q copies appeared unlikely.
Replication and collective missegregation of acentric chromosomal fragments. We encountered a similar amplification of chromosome 9q (chromosome 9q21.33-ter) within the photoconverted lineage of another PDTO-9 3D Live-Seq dataset (25 cells; 100% recovery), which allowed accurate mapping of its phylogenetic origin across four cell generations (Fig. 3a, lineage I). The amplification of chromosome 9q in one of the photoconverted cells and the reciprocal loss across all other cells within the lineage could only be reconciled by two cycles of replication and collective missegregation of an acentric chromosome 9q fragment (Supplementary Note 3). Chromatin phenotypes along the photoconverted branch suggested that replicated chromosome 9q fragments were missegregated during consecutive divisions 1.1 and 1.1.1, but were probably shielded from a third round of replication before division 1.1.1.1 by micronuclear containment (asterisk), which are known to be associated with replication defects [26][27][28] . Next, we performed additional deep sequencing of cell 1 (five copies), cells 2 and 3 (one copy) and a cell with the core karyotype (diploid), to extract single-nucleotide variants (SNVs) unique for the chromosome 9q21.33-ter region of each parental chromatid. As expected, the cell containing an amplification of chromosome 9q displayed a biased variant allele frequency of SNVs in favor of one parental allele and at a degree consistent with two rounds of replication (Fig. 3b). Notably, no indications of chromothripsis 21,29 were found within the amplified chromosome 9q fragment (Extended Data Fig. 5), suggesting that micronuclear envelope integrity was maintained 30 . Similar sub-chromosomal amplifications were observed in additional datasets (Extended Data Fig. 1 (lineages III and IV) and Extended Data Fig. 6, (lineage I). In each case, the amplified sub-chromosomal fragments were acentric, suggesting a model where a lack of spindle microtubule attachment results in non-disjunction of the replicated acentric fragments. During anaphase, the unattached replicated fragments are collectively displaced to one daughter cell, probably resulting in a stochastic nuclear or micronuclear containment. While the impact of amplified acentric fragments on tumorigenesis and tumor evolution is speculative, it is conceivable that amplified acentric chromosomal fragments act as substrates for chromothriptic events, contribute to the formation of extrachromosomal DNA or eventually fuse to chromosomes, resulting in a stable inheritance of the amplified segment.
Gross karyotype alterations result from multipolar spindle defects. While CNA events generally involve one or a few chromosomes, we identified two cells that displayed extensive reciprocal CNAs across their genomes, and their karyotypes can be classified as hopeful monsters (Fig. 3a, lineage II). In contrast with the population of hopeful monsters that we previously detected in PDTO-19b, these cells displayed an unbalanced distribution of the core karyotype of lineage I without signs of previous genome duplication, demonstrating that near-haploid and near-triploid genomes can be generated in a single division without tetraploid intermediates 20 . The uneven distribution of chromosomes matched the substantial difference in chromatin mass seen between daughter cells of division 1.1.2.2. Furthermore, the irregular shape of the metaphase plate indicated a multipolar spindle defect as a cause of the unbalanced chromosome distribution. Multipolar spindle defects are associated with a loss of spindle pole integrity 31 or supernumerary centrosomes. Indeed, the distribution of a replicated genome among excess spindle poles is likely to instigate a genome-wide misallocation of chromosomes.
Since centrosome amplification is frequently associated with whole-genome doubling 32 , we hypothesized that genome-duplicated cells, resulting from mitotic entry and exit defects, could serve as efficient substrates for the generation of hopeful monster karyotypes by frequently instigating multipolar spindle defects. Indeed, we captured a tripolar division following a re-replication event, which generated three daughter cells with reciprocal hopeful monster karyotypes ( Fig. 4a    To investigate cell fate after mitotic slippage, we recorded confocal imaging data of PDTO-9 organoids for 72 h to identify and trace mitotic slippage events. Although apoptosis was a more frequent outcome, we scored a high number of cells entering mitosis following a previous mitotic slippage event. Multipolar spindle defects were common among these downstream cell divisions, with 3 out of 12 cells displaying three or more spindle poles (Extended Data Fig. 7a,b). To confirm that multipolar spindle defects resulting from upstream mitotic slippage generate progeny with hopeful monster karyotypes, we photoconverted cells after mitotic slippage and isolated the progeny after their next   The VAF of SNVs located on the missegregated chromosome 9q21.33-ter region is 0% for cell numbers 2 and 3 (only one copy of chromosome 9q) and 50% for a cell with the core karyotype. Cell number 1 shows nearly 80% VAF for these SNVs, indicative that four of the five chromosome 9q21.33-ter copies originated from the same parental allele. The dashed red line indicates the median.
cell division. As expected, when two daughter cells were generated with equal chromatin mass, we detected polyploidy (including a few reciprocal CNAs) but no genome-wide misallocation of chromosomes (Fig. 4b, lineage I). In contrast, when downstream division resulted in three daughter cells, indicating a multipolar spindle defect, the isolated cells displayed hopeful monster karyotypes (Fig. 4b, lineage II).
Collectively, these data demonstrate that genome duplication events, irrespective of their origin, can act as unstable intermediates in the generation of progeny with gross genome-wide karyotype alterations. In addition, since the majority of cell divisions that follow mitotic slippage appeared fairly normal (bipolar) (Fig. 4b (lineage I) and Extended Data Fig. 7a,b), mitotic slippage may represent an important pathway for whole-genome doubling during carcinogenesis and tumor progression, which is a feature of many human cancers 2,3,33 .
Cells with novel karyotypes frequently remain proliferative. As presented here, 3D Live-Seq of PDTOs enables the reconstruction of karyotype alterations across multiple consecutive cell generations. It reveals the immediate genomic consequences of CIN (that is, the temporal dynamics and patterns by which new karyotypes  are constantly being generated, upon which selection pressures act to shape the genomic evolution of advanced human cancers). Our 3D Live-Seq datasets suggest that intrinsic negative selection against de novo karyotypes, including grossly altered hopeful monster karyotypes, is often not absolute or instantaneous. To derive quantitative evidence for the strength of selection experienced by novel karyotypes, we used a stochastic branching process model to simulate CNA evolution and clonal selection within organoids using cell cycle parameters and measures of karyotype diversity extracted from our datasets. The results were best explained by a neutral drift model, although the power to detect subtle selection pressures was restricted due to the limited dataset size (Extended Data Fig. 8a-d).
To further explore selective forces in 3D Live-Seq data, we developed a likelihood ratio test to assess whether there was evidence for a change in birth rate following mitotic events. Using this approach, we were able to aggregate all four mitotic trees of PDTO-9 and found evidence for a decrease in birth rate after a mitotic error (P = 0.029; Supplementary Note 4). Taken together, this implied that many de novo karyotypes remain proliferative over a timescale of a few cell generations and experience neutral or subtle negative fitness effects.

Gross deviations of the core karyotype are generally less fit.
To investigate cell intrinsic selection pressures against de novo karyotypes in closer detail, we captured a high number of individual cell lineages in PDTOs during early stages of outgrowth (days 4-6; up to ~30 cells in size) and mature stages (days 9-11; hundreds of cells). Overall, the level of CIN remained constant from young to mature organoids (Fig. 5a). Apoptosis was readily detected during both stages of outgrowth and, in agreement with previous reports 7 , we detected elevated levels of apoptosis among the progeny of erroneous cell divisions, in particular after multipolar spindle defects (Fig. 5b,c and Extended Data Fig. 7c). For cells that remained proliferative, the apparent slight shift towards longer cell cycle durations immediately after an error was not significant (Fig. 5d; Wilcoxon test, P = 0.64). However, a comprehensive analysis of mitotic tree structures using our derived likelihood ratio test was, again, able to indicate evidence for a decrease in birth rate following mitotic errors (P = 5 × 10 −6 ; Supplementary Note 4).
Thus, a subset of de novo karyotypes are subject to stringent intrinsic negative selection, with large deviations from the core karyotype originating from multipolar spindle defects being least   The difference in apoptotic rates was statistically significant, as determined by chi-squared test (P = 0.014 (young) and P = 0.0026 (mature)). c, As in b, but the fate analysis is shown for normal (N), chromatin bridge (CB), lagging chromatin (LC) and multipolar cell divisions (MP). d, Ratio between the cell cycle duration before a division and the progeny's cell cycle duration (averaged) after that division. Data from young and mature PDTOs were pooled and subdivided by mitotic fidelity. The difference in cell cycle duration ratios between normal (94 branches) and erroneous divisions (76 branches) was not significant (data not normally distributed; P = 6.104 × 10 −13 (Saphiro); P = 0.64 (two-sided Wilcoxon test)). In the box and whisker plots, the boxes represent quartiles 2 and 3, the horizontal lines represent median values and the whiskers represent minimum and maximum values within 1.5× the interquartile range. The data points indicate outlier values that deviate by more than 1.5× the interquartile range. e, Top: schematic of karyotype analysis of successfully formed organoids (bulk) grown from single (fit) cells derived from a parental (clonal) PDTO. Bottom: karyotype heatmaps of sequenced PDTOs (bulk; similar sized) as a proxy for the karyotype of the seeding cell. In total, 67 organoids shared the reference PDTO-9 core karyotype, seven were polyploid and five had de novo localized CNAs.
fit. In contrast, the majority of cells with localized karyotype alterations that result from chromatin bridges or lagging chromatin phenotypes clearly remain proliferative, although these cells seem to experience a subtle fitness cost in the aggregate.

Cells with novel karyotypes can seed new PDTOs.
To functionally probe the fitness of novel karyotypes, we assessed the organoid-forming capacity of single tumor cells derived from one parental clonal PDTO. Subsequent karyotype analysis of the newly formed individual organoids (bulk) acts as a proxy for the karyotype of the seeding single cell and can confirm sustained fitness of de novo CNAs. Indeed, cells with (sub-)chromosomal karyotype alterations and even whole-genome duplications were able to develop into new PDTOs (Fig. 5e and Extended Data Fig. 8e), in line with our general observations, demonstrating that from a cell intrinsic and short-term perspective, some cells with novel karyotypes maintain a high level of fitness.
Punctuated and localized karyotype alterations are ongoing in mature tumors. To investigate whether advanced CRCs generate a similar degree of karyotype diversity in vivo before the reshaping of the karyotypic landscape as a consequence of selection, we adapted our sequencing methodology to measure karyotype diversity within clonal lineages of freshly frozen human CRC samples. To maximize recent ancestry within the sampled cell population, we excised small fragments from individual CRC glands, which are analogous to clonal intestinal crypts of the normal epithelium 34 . Subsequent single-cell isolation and prospective genome analyses displayed substantial karyotype diversity between cells, with almost 30% of cells carrying CNAs that deviated from the core gland karyotype-a number that was in agreement with our PDTO data ( Fig. 6a and Extended Data Fig. 9). In contrast, only a few single cells from a wild-type intestinal crypt displayed de novo CNAs, supporting the sensitivity of our assay, in agreement with an earlier report on basal CIN levels in healthy colon tissue 35 (Fig. 6b). Notably, by capturing a large fraction of clonal tumor cells (44%), we detected the emergence and propagation of a de novo CNA in vivo (Fig. 6a, lineage I). Moreover, we readily identified cells with hopeful monster karyotypes (Fig. 6a, lineage II), further demonstrating that the generation of cells with genome-wide karyotype alterations is ongoing at advanced cancer stages.

Discussion
Our work reveals the genomic consequences of CIN across consecutive cell generations. By reconstructing intermediate genomic states, we demonstrate that both genome-wide (punctuated) and localized (singular, Extended Data Fig. 10) karyotype alterations are ongoing and prevalent in advanced CRC. The substantially reduced fitness of many hopeful monster karyotypes is not surprising given the extent of their karyotype alterations. Nevertheless, low probabilities of increased fitness can still be of potential impact considering that tumors can consist of hundreds of billions of cells 9 . Thus, while most diversity provides incremental variation around a proven-fit aneuploid genome, millions of hopeful monster karyotypes represent extreme, potentially adaptive phenotypes within a fitness landscape. The ongoing generation of karyotype diversity in advanced cancers is in agreement with recent analyses of tumor genomes from large patient cohorts 2,3 . However, as most of these studies measure the net outcome of evolution, it remains difficult to disentangle the rate at which karyotype diversity is generated and the selection pressures that act on them to shape an aneuploid landscape over time. Therefore, our insights into the alteration rate and diversification patterns of tumor karyotypes provide valuable parameters to understand the historic trajectories and temporal dynamics of an evolving aneuploid genome. In addition, our work suggests that aneuploid tumor karyotypes can, in principle, be generated in a single erroneous cell  division, providing support for punctuated evolution of malignant karyotypes during the earliest stages of carcinogenesis 13,14,36 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41588-021-00891-2.
Karyotyping. PDTOs were treated with 0.1 µg ml −1 colcemid (Thermo Fisher Scientific) in culture medium for 12 h. Then, PDTOs were dissociated into single cells using TrypLE Express, incubated in a hypotonic 27 mM trisodium citrate solution at 37 °C for 10 min and fixed in methanol:acetic acid solution (3:1). Fixed cells were dropped on a microscope slide, mounted with DAPI-containing Vectashield and imaged on a Zeiss Axio Imager Z1 microscope with a 63× objective. In total, 36 metaphase spreads were imaged for each PDTO line and quantified by manual chromosome counting (Extended Data Fig. 10a).

Confocal live-cell microscopy and image analysis.
To support live-cell microscopy of cell divisions and subsequent photoconversion of daughter cells, PDTO cultures expressing transgenic H2B-Dendra2-IRES-Puromycin were passaged 5-7 d before imaging. PDTOs were harvested 24 h before imaging and resuspended in an ice-cold mix of culture media containing 50% vol/vol BME. The organoid suspension was then seeded in an ice-cold glass-bottom WillCo dish (WillCo Wells) coated with a thin film of BME. PDTOs were allowed to settle on ice before BME polymerization at 37 °C and the addition of culture media. Cells about to go into mitosis were identified by chromosome condensation and imaged until completion of mitosis with a confocal laser-scanning microscope (Leica SP8X, z step 1 µm, 30-s frame rate) equipped with atmospheric and temperature control. Daughter cells of imaged cell divisions were marked via photoconversion of Dendra2 (405 nm laser; 1 min; 10% laser power; per region of interest). Imaging data were analyzed and 3D rendered with Imaris version 9.3 image analysis software (Oxford Instruments).
To support long-term imaging of organoid outgrowth, a dissociated PDTO suspension was filtered twice through a CellTrics 10-μm sieve (Sysmex) to obtain a pure single-cell suspension. Single cells were seeded in BME-coated glass-bottom WillCo dishes as described previously. At 30 h after seeding, the outgrowth of cell doublets was imaged (3-min frame rate; z step = 1 µm). The culture media was refreshed each day to prevent excessive evaporation. Raw data were converted to videos using an ImageJ macro as described 39,40 . Lineage traces were generated manually.
Spinning disk live-cell imaging. To support long-term imaging of organoid outgrowth using a spinning disk confocal system (Nikon), 3-d-old organoids carrying a HIST1H2BC-Dendra2 knock-in or expressing transgenic H2B-Dendra were cultured and seeded in BME-coated glass-bottom WillCo dishes as previously described. The outgrowth was captured for 72 h on days 4-6 and days 9-11 (4.5-min frame rate; z step = 1.4 µm). The imaging data were analyzed with Fiji (ImageJ).

Single-cell isolation of PDTOs.
To generate clonal PDTOs, a dissociated PDTO suspension was filtered twice through a CellTrics 10-μm sieve to obtain a pure single-cell suspension and seeded at medium density in BME. Clonal PDTOs were grown into structures of 60-100 cells before harvesting. Before single-cell isolation, a single PDTO was transferred to a WillCo dish containing phosphate-buffered saline (PBS) (and 1:1,000 Syto 11 (S7573; Thermo Fisher Scientific) for PDTOs lacking H2B-Dendra2-IRES-Puromycin expression) for image acquisition (LSM510; Zeiss). Raw Z-stack data were analyzed with Fiji (ImageJ) to obtain the exact cell number, and a 3D render of the imaged organoid was generated using Imaris version 9.3 image analysis software. To isolate single cells, the imaged organoid was mouth-pipetted to a drop of PBS on a siliconized object slide (Sigmacote; Sigma-Aldrich) and mounted on a fluorescence stereo microscope (SMZ18; Nikon) equipped with a Sola LED (Lumencor) for fluorescence-assisted single-cell picking. The organoid was transferred to a 20-μl drop of trypsin-EDTA (Lonza) for dissociation (30 s). After removing trypsin-EDTA, PBS containing 0.5% bovine serum albumin (Sigma-Aldrich) was added immediately. The organoid structure was further dissociated by repeated mouth-pipetting using a 30-μm glass capillary (custom made). Single cells were transferred to PCR tubes containing either 30 μl mineral oil (Affymetrix) or 4 μl H 2 O and each tube was immediately spun down and stored ice-cold. Photoconverted daughter cells were isolated as previously described and visualized with a custom filter set on SMZ18 (FF01-433/24, FF458-Di02 and FF01-482/35). The remaining cells were pooled and transferred to a single PCR tube as a bulk sample. An overview of the single-cell datasets is presented in Supplementary Table 2.
Organoid-forming capacity of PDTO-9 cells. PDTO-9 organoids were cultured and dissociated as previously described. Single cells obtained from a single parental clonal organoid were reseeded and cultured for 9 d as previously described.
Organoids of similar size were isolated and individually bulk sequenced as a proxy for the karyotype of the seeding cell.
Fresh-frozen CRC tissue specimens. Fresh-frozen tissue samples were collected from University College London Hospitals under ethical approval 11/LO/1613 and via the University College London Hospital Biobank (15/YH/0311) and processed as previously described 41 . All surgically resected samples were collected from patients who had given informed consent.
Single-cell isolation of CRC glands. Individual tumor glands were mechanically separated from thawed CRC tissue pieces using a fluorescence stereo microscope (SMZ18; Nikon) and transferred to a WillCo dish containing PBS and 1:1,000 Syto 11 (S7573; Thermo Fisher Scientific) for image acquisition (LSM510; Zeiss). A small fragment of individual CRC glands was resected for single-cell isolation. Raw Z-stack data were analyzed with Fiji (ImageJ) to obtain the exact cell number of the resected fragment, and a 3D render of the imaged gland was generated using Imaris version 9.3 image analysis software. Single cells of the resected fragment were subsequently isolated as described for PDTOs. The remaining part of the gland was transferred to a single PCR tube as a bulk sample.
Single-cell sequencing. Single-cell sequencing was performed on clonal PDTOs with and without photoconverted daughter cells and on CRC glands.
Single-cell sequencing of clonal PDTOs and photoconverted daughter cells. For cells expressing H2B-Dendra2-IRES-Puromycin, single-cell lysis and whole-genome amplification (WGA) were performed using the REPLI-G single-cell kit (Qiagen) according to the manufacturer's instructions. A positive control (multiple organoids or HEK293T cells) and negative control (3 µl H 2 O) were included during each WGA reaction. Amplified DNA was purified by phenol:chloroform:isoamyl alcohol (25:24:1) treatment, followed by precipitation with ethanol. To assess the quality of amplified DNA, a multiplex PCR was used to simultaneously amplify nine loci across the human genome in each single cell: SCYL1, PPP5C, JMJD6, ACTR10, ROCK2, SCAP, SCAMP3, SCARB2 and XPOT 42 . DNA libraries of each single cell, with two or more amplified products in the multiplex PCR, were constructed with the TruSeq Nano DNA library preparation kit (Illumina). The concentration of all libraries was quantified using a Qubit dsDNA HS Assay kit. Library fragment size profiling was performed using the Agilent TapeStation system using Agilent High Sensitivity D1000 ScreenTapes (Agilent Technologies), and subsequently all libraries with unique indices were pooled equimolarly.
Single-cell sequencing of clonal PDTOs and CRC glands. For cells stained with Syto 11, cell lysis was performed with 0.63 mUA protease (Qiagen) and 1× NEBuffer 4 at 50 °C for 2 h, followed by a protease inactivation incubation for 20 min at 80 °C. DNA libraries were prepared from individual cells with a modified KAPA HyperPlus kit protocol (one-seventh of all volumes; Roche). Cleaned-up pooled single-cell libraries were quantified using a Qubit dsDNA HS Assay kit and analyzed using the Agilent 2100 Bioanalyzer HS kit.
Low-coverage whole-genome single-end 75-base pair sequencing was performed on an Illumina NextSeq 500 aimed at generating ~5 million reads per single cell for the WGA TruSeq libraries and more than 100,000 reads for the KAPA libraries. For a subset of WGA TruSeq libraries (cell numbers 1-3 and 15 of the clonal PDTO dataset presented in Fig. 3), whole-genome paired-end sequencing (2 × 150 base pairs) was performed on an Illumina NovaSeq 6000 to a coverage of 15× (cell numbers 2 and 3) or 30× (cell numbers 1 and 15). Genomic analysis. Sequence reads were aligned to the human genome reference (hg19/GRCh37) using the Burrows-Wheeler Aligner mapping tool (BWA-MEM; version 0.7.15). Duplicated sequence reads were marked by Sambamba (version 0.6.5) and realigned using the Genome Analysis Toolkit IndelRealigner (GATK; version 3.8), and sequence read quality scores were recalibrated with GATK BaseRecalibrator.
The copy-number status of each single cell was analyzed using Ginkgo 43 . Briefly, the created BAM files were converted to BED format files using BEDtools (version 2.25.0). Then, the aligned reads were binned into 1-Mb variable-length intervals across the genome, normalized and corrected for GC biases, and segmented in regions with an equal copy-number status. Thereafter, a final integer copy-number profile was assigned to each single cell. Single cells that exhibited low-quality sequencing data, defined as high variation between copy-number ratios of segmented regions, were excluded. One exception was made for cell number 9 of dataset 04032020, as this cell clearly belonged to the lineage with segmental loss of chromosome 1p.
To identify high-confidence copy-number variations, regions of at least 25 Mb (excluding bins greater than 5 Mb) with copy-number values deviating >0.6 from the average ploidy were considered to indicate losses or gains. Copy-number variations smaller than 25 Mb were only included if both the loss and gain of that particular region were identified in the clonally expanded organoid, photoconverted daughter cells or CRC gland. To avoid missing copy-number variations, less confident variations (>0.3 and ≤0.6 below or above the average ploidy) with >70% reciprocal overlap between high-confidence variations were reviewed. Heatmaps were generated using the R (version 3.4.3) package ggplot2. All de novo CNA events observed in PDTO datasets are shown in Extended Data Fig. 10b using the R package karyoploteR.
For cell numbers 1-3 and 15 of the clonal PDTO dataset presented in Fig. 3 with 15× or 30× coverage, SNV and indel calling was performed using GATK HaplotypeCaller. The allele frequency of variants located on chromosome 9 (genomic location = 89502430-141213431) in sequences 1 (SC11-20190409_1) and 15 (SC15-20190409_1) (core) were calculated by subtracting all of the variants present on the allele of sequences 2 (SC10-20190409_1) and 3 (SC17-20190409_1) (Fig. 3b). Density plots of the allele frequencies of variants were generated using the R package ggplot2. Manta (version 0.29.5) was used with standard settings to detect structural variants. To exclude structural variants introduced by the WGA (short-range chimeras of the inverted orientation), we established a threshold of structural variant length for which a similar count in inverted and non-inverted structural variants was observed. Next, we estimated the frequency of structural variants (subdivided based on threshold structural variant length) on all control chromosomes and on the missegregated chromosome to determine enrichment of structural variants in the missegregated chromosome. Parameter estimation and model selection. We used a stochastic birth-death branching process to model the outgrowth of PDTOs (Supplementary Note 4). Since apoptotic events were rare in the 3D Live-Seq datasets, we modeled only cell birth with rate b (per day). During each cell division, n 1 reciprocal chromosome-level CNAs and n 2 reciprocal arm-level CNAs were introduced to daughter cells, where n 1 = ~Pois(μ 1 ) and n 2 = ~Pois(μ 2 ). We assumed two models of evolution: a neutral model where all cells have the same fitness; and a selection model in which cells with novel mutation(s) have a selective (dis) advantage, s, given by bd bp = 1 + s, where b p and b d are the birth rate of the parent and daughter cells respectively. The model parameters, θ = (μ 1 ,μ 2 ,b,s), were inferred using approximate Bayesian computation with sequential Monte Carlo sampling 44,45 . Using simulated data, we validated that the true parameters were inferred accurately. For datasets without recorded mitotic trees, we used single-cell copy-number profiles to estimate mutation rates under the neutral model (assuming b = 0.5). For datasets with mitotic trees, we estimated all parameters under both models and selected the model that best fit the data using deviance information criteria 46 . We also examined the power of our approach to detect negative selection at varying selection strengths and cell population sizes. This was done by simulating 100 datasets for each parameter setting, with the parameters estimated from real data (μ 1 = 0.1, μ 2 = 0.1 and b = 0.4) using the approximate Bayesian computation rejection algorithm, and computing deviance information criteria under the neutral and selection models for each dataset.
Statistics and reproducibility. Imaging stills in the manuscript represent (time-compressed) representations of the entire captured imaging data. Representative snapshots of organoid outgrowth were extracted from imaging data at evenly distributed time points across the imaging span. The numbers of samples analyzed are indicated when applicable.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
BAM files of the single-cell sequencing data are available through controlled access at the European Genome-phenome Archive (EGA), which is hosted by the European Bioinformatics Institute and Centre for Genomic Regulation (https:// ega-archive.org) under accession number EGAS00001003812. Data access requests will be evaluated by the University Medical Center Utrecht Department of Genetics Data Access Board (EGAC00001000432) and transferred on completion of an agreement and authorization by the medical ethical committee at the University Medical Center Utrecht at the request of Hubrecht Organoid Technology to ensure compliance with the Dutch Medical Research Involving Human Subjects Act. The copy-number segment calls of all sequenced single cells are publicly available on Zenodo (https://doi.org/10.5281/zenodo.4732372). Fig. 3 | Copy-number profiles of hopeful monster karyotypes. a, Reciprocal copy-number profiles of two PDTO-19b cells with hopeful monster karyotypes are shown (Fig. 1c, III). Since the ancestor karyotype of these two cells is unclear, reciprocity between hopeful monster karyotypes refers to the observation that all chromosomal segments add up to an even number. Grey dots indicate bin ratios (Bin size: 1 Mbp). Integer copy-number states are indicated with black (no deviation of a tetraploid state), blue (loss) and orange (gain) lines. Below: The copy-number profile of the merged sequencing data of the two hopeful monsters mimics a duplicated core karyotype of the sequenced PDTO-19b organoid. Orange boxes indicate CNAs in the merged profile that deviate from a duplication of the core karyotype, indicating that these CNAs were present in the ancestor cell. b, Reciprocal copy-number profiles of three PDTO-9 cells with hopeful monster karyotypes are shown (Fig. 4a). Grey dots indicate bin ratios (Bin size: 1 Mbp). Integer copy-number states are indicated with black (no deviation of a tetraploid state), blue (loss) and orange (gain) lines. Adding up the copy-number state of chromosome segments across all three hopeful monster karyotypes results in a twice replicated genome ('ploidy 8', not shown). Below: The copy-number profile of the merged sequencing data mimics the core karyotype carrying a de novo loss of 1pter-p34.2 (orange box).