An in vitro system for telomerase-mediated escape from natural telomere crisis

To gain a clearer insight into the nature of SVs that arise during telomere crisis, we developed an in vitro system in which we could reproduce telomere crisis and generate a large number of post-crisis clones. MRC5 human lung fibroblasts were chosen to model telomere crisis since they lack telomerase activity and as a consequence have a well-defined in vitro replicative potential determined by telomere attrition. To bypass senescence, the Rb and p21 pathways were inactivated by infecting the population of MRC5 cells with retrovirus-bearing shRNAs targeting the respective transcripts (Supplementary Fig. ^2A). This population of MRC5/Rbsh/p21sh was then endowed with an inducible CRISPR activation system (iCRISPRa) to activate the TERT promoter and induce telomerase expression (Supplementary Fig. ^2B). The iCRISPRa system employed a doxycycline-inducible nuclease-dead Cas9 fused to a tripartite transcriptional activator (VP64-p65-Rta)³² and four gRNAs targeting the TERT promoter (Fig. 2a and Supplementary Fig. ^2B). The addition of doxycycline (dox) to MRC5/Rbsh/p21sh/iCRISPRa-TERT cells resulted in induction of TERT mRNA within 96h, whereas without dox, TERT transcripts are undetectable in this cell line (P<0.001, Fig. 2b). A similar dox-induced increase in mRNA expression was noted upon the introduction of sgRNAs to a control gene (Supplementary Fig. ^2C). Induction of telomerase activity was readily detectable in a TRAP (telomerase-repeated amplification protocol) assay (Fig. 2c). However, the induced TERT mRNA levels and the TRAP activity were significantly lower than in telomerase-positive control cell lines. The relatively weak telomerase activity in this system harmonizes with recent work showing that cancer-associated TERT promoter mutations initially result in low levels of telomerase activity that is not sufficient to maintain bulk telomere length³³.

Open in a separate window

Fig. 2

An in vitro system for a telomerase-mediated escape from natural telomere crisis.

a Immunoblot for dCas9-VPR (using a Cas9 Ab) in MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with or without doxycycline treatment for 96h (see also Supplementary Fig. ^{2A, B}). The blot shown is representative of at least two experiments. b qPCR of TERT mRNA expression in RPE-1, HCT116, U2OS (n=2), and MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with and without doxycycline treatment (n=6). Values are normalized to β-actin mRNA. Error bars represent means±SDs; P value from two-tailed Student’s t test; ****P<0.0001. c TRAP assay on MRC5 and MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with and without doxycycline treatment for indicated time periods. HCT116 and 293T (Phoenix) cells are included as positive controls. IC=internal control PCR product at 36bp. The gel shown is representative of at least two experiments. d Growth curve of parental MRC5 cells, MRC5/Rbsh/p21sh cells, and MRC5/Rbsh/p21sh/iCRISPRa-TERT cells grown with or without doxycycline. Arrows indicate when each construct was introduced. Days in culture represent total time in culture from parental MRC5 cells to late passage MRC5/Rbsh/p21sh/iCRISPRa-TERT cells. Time points for telomere analysis (presented in Fig. 3) and the approximate onset of senescence in the parental MRC5 cells are indicated. e STELA of XpYp telomeres in MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with or without doxycycline treatment at 70 and 150 days of culture. f Quantification of band intensity in e, with background signal subtracted. Data from two independent experiments (see also Supplementary Fig. ^2D) were analyzed with two-way ANOVA with multiple comparisons; all points at day 70 are not significant, day 150; 5–6kb P=0.0005; 4–5kb P=0.0199; 2–3kb P=0.005; 1.5–2kb P=0.0025. Biological replicates represent cells at approximately the same days in culture (±5 days). g Genomic blot of telomeric MboI/AluI fragments in MRC5/Rbsh/p21sh/iCRISPRa-TERT cells grown with or without doxycycline at the indicated time points. The blot shown is representative of at least two experiments.

At approximately 120 days after the start of the experiment (55 days with dox), the MRC5/Rbsh/p21sh/iCRISPRa-TERT population was proliferating faster than their untreated counterparts (Fig. 2d). Inspection of individual telomere lengths using single-telomere-length analysis (STELA³⁴) revealed that although telomerase expression was sufficient to allow the cells to proliferate, it was not sufficient to maintain bulk telomere length (Fig. 2e). After 150 days of continuous culture, the majority (86%) of XpYp telomeres in induced MRC5/Rbsh/p21sh/iCRISPRa-TERT cells were between 1 and 4kb compared to 40% in uninduced cells (Fig. 2f and Supplementary Fig. ^2D). Consistent with this, genomic blotting showed bulk telomere shortening in both induced and uninduced cells (Fig. 2g). These telomere dynamics are consistent with the expectation that in the culture without telomerase, cells with critically short telomeres will preferentially be lost, leading to a surviving population with relatively longer telomeres. In contrast, cells in the induced culture with (low) telomerase activity have the ability to elongate the shortest telomeres. As a result, the induced cells are expected to tolerate telomere attrition better and present with overall shorter telomeres at later time points.

Dissipating telomere crisis in MRC5/Rbsh/p21sh/iCRISPRa-TERT cells

To confirm that the MRC5/Rbsh/p21sh/iCRISPRa-TERT cells experienced telomere crisis before the induction of telomerase increased their proliferation rate, we investigated cells at various time points from the start of the experiment (Fig. 2d). Metaphase spreads showed both induced and uninduced MRC5/Rbsh/p21sh/iCRISPRa-TERT cells contained dicentric and multicentric chromosomes (Fig. 3a) and genomic blots showed high-molecular weight telomere bands consistent with fused telomeres (Fig. 2g). As expected from the ability of telomerase to counteract the formation of critically short telomeres, at 125 days after the start of the experiment, induced cells had significantly fewer fusions than untreated cells (21% vs. 40%, P<0.05; Fig. 3b and Supplementary Fig. ³). PCR-mediated detection of fusions between the Tel Bam 11 family of telomeres^35,36 confirmed these dynamics (Fig. 3c). Quantification of the fusion frequency showed a significant reduction in the number of fusions per haploid genome in the induced population (day 110, P<0.01; Fig. 3d). Consistent with telomerase-mediated genome stabilization, there was a trend toward a lower level of 53BP1-marked DNA damage foci at later time points (Fig. 3e, f) and the percentage of cells with micronuclei (an indicator of genome instability) was significantly reduced at day 110 (P<0.05; Fig. 3g). Taken together, these data indicate that after a period of genomic instability induced by critically short telomeres, iCRISPRa-mediated telomerase activation is sufficient to partially stabilize the genome and allow the MRC5/Rbsh/p21sh/iCRISPRa-TERT cells to navigate the deleterious effects of telomere crisis.

Open in a separate window

Fig. 3

Dissipating telomere crisis in MRC5/Rbsh/p21sh/iCRISPRa-TERT cells.

a Metaphase spreads from MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with and without doxycycline at day 95. Telomeres are detected with a telomeric repeat PNA probe (TelG, red), and centromeres are detected with a probe for CENPB (green). DNA was stained with DAPI (gray). Chromosome fusions are indicated by white arrowheads. b Quantification of the percentage of metaphase spreads with at least one fusion after the indicated days of continuous culture for MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with and without doxycycline (see also Supplementary Fig. ^3A), two-tailed Student’s t test; ns not significant, *P=0.0422. c Gel showing products of telomere fusion PCR on MRC5/Rbsh/p21sh/iCRISPRa-TERT cells cultured with and without doxycycline for the indicated time. Each lane represents an independent replicate PCR reaction. Telomere fusion products are detected by hybridization with a probe for the 21q telomere (see “Methods”), and the control XpYp PCR product is detected with ethidium bromide staining. d Quantification of the number of telomere fusion products per haploid genome using the assay shown in panel c. Each dot represents a single PCR reaction. Reactions from two independent biological replicates are shown, two-tailed Student’s t test; ns not significant; **P=0.0061. e Detection of micronuclei (arrowheads) and DNA damage foci using indirect immunofluorescence for 53BP1 (red) in the indicated cells. DNA is stained with DAPI (blue); scale bar (white)=10µm. f Quantification of the percentage of cells with >10 53BP1 foci at the indicated time points; two-tailed Student’s t test; ns not significant. g Quantification of the percentage of cells with micronuclei after the indicated days in culture, two-tailed Student’s t test; ns not significant; *P=0.0157. In panels b, f, and g, error bars indicate means and standard deviations from three independent biological replicates.

Genomic screening of post-crisis clones

To assess the genome structure of proliferating post-crisis cells, single-cell clones were isolated from induced MRC5/Rbsh/p21sh/iCRISPRa-TERT cells at day 120 (“Y clones”) and day 150 (“Z clones”) (Fig. 4a). The clonal yield at day 150 was greater than at day 120 in induced cells, but no clones could be isolated from the uninduced population at either time point. The lower clonal yield at day 120 may be due to incomplete stabilization of the telomeres since clones from this time point showed a higher burden of fused telomeres than those derived from day 150 (Supplementary Fig. ^4A). Post-crisis clones from both time points showed evidence of ultrashort telomeres and reduced telomere length (Supplementary Fig. ^{4B, C}). Telomerase activity in post-crisis clones was comparable to the parental induced population, indicating that clone viability was not due to selection for increased telomerase activity (Supplementary Fig. ^4D). To generate control clones that had not passed through a period of telomere crisis, early-passage MRC5 cells were infected with a retrovirus expressing hTERT and single-cell clones were isolated (Supplementary Fig. ^4E). Genome profiling with a low-pass (~5X) WGS was performed on eight hTERT-expressing control clones (CT clones), 36 Y clones from day 120, and 82 Z clones from day 150 (Supplementary Table ²).

Open in a separate window

Fig. 4

Genomic screening of post-crisis clones.

a Growth curve of MRC5/Rbsh/p21sh/iCRISPRa-TERT cells with and without doxycycline, indicating the time points at which single-cell clones were derived (day 120 and day 150). b Circular heatmap showing genome-wide binned purity- and ploidy-transformed read depth (in units of CN across 118 low-pass WGS-profiled clones. Heatmap rows correspond to concentric rings in the heatmap. Clones are clustered with respect to genome-wide copy number profile similarity (see “Methods”). c Zoomed-in portion of chromosomes 12 and 21 that underwent copy number alterations in a majority of the clones, clustered based on their coverage across these regions. Clusters are named with respect to their consensus copy number pattern, and on the basis of high-depth WGS analyses presented in Fig. 5. Chromosome 21 gain n=6 clones; unrearranged n=38; chromothripsis-like n=1; arm loss n=6; early BFB-like n=20; BFB-like n=47.

Analysis of genome-wide read depth across 118 clones from both day 120 (Y clones) and day 150 (Z clones) demonstrated predominantly diploid genomes with a striking enrichment of clones with DNA loss on most of chromosome 12p (63%, 74/118, Fig. 4b). Within the other 44 samples, we observed a subset of clones (5%, 6/118) with gains of chromosome 21q. As expected, control CT clones showed no evidence of SVs or copy number variants (Supplementary Fig. ^4F). Hierarchical clustering of all clones by their coverage on chromosomes 12p and 21q revealed six distinct clusters (Fig. 4c). A minority of clones were diploid on chromosomes 12 and 21 and elsewhere in the genome and are therefore designated as “unrearranged” (32% of clones, 38/118). Of note, the unrearranged group was enriched in day 120 (Y) samples compared to day 150 (Z) samples (P=1.79×10⁻⁹, odds ratio 14.7, Fisher’s exact test; Fig. 4c), suggesting that these clones may have largely avoided crisis prior to telomerase induction. The cluster of clones with 21q gain was diploid on 12p.

The remaining 74 clones (63%) all showed a heterogeneous pattern of copy number alterations targeting 12p (Fig. 4c). One out of the 118 clones (0.8%) displayed the singular pattern of distinct interspersed losses that resembled chromothripsis. Complete loss of one copy of 12p (“arm loss”) was found in a cluster of six clones (6/118, 5%). The second cluster of 67 clones all shared a breakpoint near the distal end of 12p and a large deletion starting ~9Mbp from the centromere. These clones were differentiated into two clusters by the presence or absence of an amplification ~8–9Mbp from the 12p telomere. In the 47 clones that contained this amplification, aggregated consensus read-depth profiles revealed stepwise gains at the distal end of 12p, a pattern reminiscent of BFB cycles (Supplementary Fig. ^5A). This cluster was therefore labeled “BFB-like”, a designation that is further supported by the data presented below. The 20 clones (17%) that lack the amplicon ~8–9Mbp harbored varying boundaries of the shared larger deletion; based on the analysis described below, we designate these as “early BFB-like”. In summary, these low-pass WGS copy number profiles indicated a limited set of distinct lineages surviving telomere crisis, with at least two lineages independently converging on 12p.

High-resolution reconstruction and lineage of post-crisis genomes

To gain further insight into structural variant evolution along these lineages, we chose 13 representative clones spanning the five clusters with rearrangements involving 12p for high-depth WGS to a median read depth of 50× (range: 30–88). Phylogenies derived from genome-wide SNV patterns demonstrated a median branch length of 551 SNVs (range: 9–2409), a low mutation density (<1 SNV/Mbp) that is consistent with previous WGS studies of clones in cell culture³⁷. This analysis revealed four major clades (Fig. 5a). These clades had good concordance with copy number alteration and rearrangement junction patterns in the same 12p region, suggesting these clones represent distinct post-crisis evolutionary lineages (Fig. 5b).

Open in a separate window

Fig. 5

High-resolution reconstruction and lineage of post-crisis genomes.

a SNV-based phylogeny inferred across 13 high-depth WGS clones and heatmap of variant allele fractions (VAF) for SNVs detected among two or more clones. For simplicity, private SNVs (those found only in a single clone) are not shown. b Heatmap of chromosome 12p copy numbers and variant junction patterns in chromosome 12 (see the text and “Methods”). c Proposed tree showing distinct trajectories of structural variant evolution following 12p attrition and subsequent telomere crisis. Each terminal node in the tree is associated with a unique 12p profile comprising a representative binned read-depth pattern (bottom track) from one or more clones mapping to an identical junction-balanced genome graph (second track from bottom). The top track in each profile represents a reconstruction of the rearranged allele. Each allele is a walk of genomic intervals and reference/variant junctions that, in combination with an unrearranged 12p allele (not shown), sum to the observed genome graph (see “Methods”). Two distinct arrows linking Y11 and Y15 demonstrate that these clones are distinct lineages (based on divergent SNV patterns, see panel a), that converge to identical WGS 12p CN profiles (although with likely distinct breakpoints inside the 12p centromere unmappable by WGS, see the text).

In order to further reconcile the shared and distinct rearrangement junctions present in the evolution of these clones, we carried out a local assembly of rearrangement junctions and junction balance analysis (see “Methods”³), which revealed seven distinct junction-balanced genome graphs spanning 12p (Fig. 5c). With the exception of the chromothriptic lineage (see below), each of these distinct lineages was represented by more than one post-crisis clone.

To reconstruct a set of linear alleles that parsimoniously explain these different genome graph patterns³ (Fig. 5c), we applied gGnome to the data (see the section “Joint Reconstruction of allelic evolution in MRC5”). We constrained our model to contain one intact allele of chromosome 12 for the following reasons: (1) karyotypes and chromosome painting showed a single copy of chromosome 12 was altered in the post-crisis clones (see Fig. 6); and (2) rearrangement of one allele is more likely than rearrangement across two alleles. Application of this constraint to the full set of MRC5 clones in a joint inference revealed a parsimonious set of rearranged alleles that explained the observed collection of clonally related junction-balanced graphs (Fig. 5c).

Open in a separate window

Fig. 6

Resolution of BFB cycles in telomere crisis.

a DAPI banded karyotypes of post-crisis clones Z43 and Y8 showing a rearranged chromosome 12 (red star) and loss of one copy of chromosome 21 (dashed box) (see also Supplementary Fig. ^6A). b Representative metaphase spreads of clone Z43 and Y8 hybridized with whole-chromosome pairs for chromosomes 12 (green) and 21 (red). DNA was stained with DAPI (gray). Insets show enlarged images of the 12–21 derivative marker chromosome and intact copies of the sister alleles (see also Supplementary Fig. ^6B). c Images of derivative chromosome 12:21 from representative clones from each branch of the evolution of chromosome 12 post-crisis (according to the analysis in Fig. 5c). Metaphases were hybridized with whole-chromosome pairs for 12 (green) and 21 (red). DNA was stained with DAPI (gray).

Analyzing the clonal evolution of these rearranged 12p alleles, we identified eight clones demonstrating progressive stages of a BFB cycle. This complex variant evolved after a long-range inversion junction (j1) joined a distal end of 12p to its peri-centromere. This junction was followed by subsequent fold-back inversion junctions (j2, j3, j4), clustered at the 8–9 Mbp focus on 12p, which are present in two different sets of post-crisis clones (Early BFB, BFB, Fig. 5c). The earliest of the fold-back inversion junctions (j2) in the BFB lineage was associated with a cluster of 3G or C mutations within 2kbp of each other, consistent with APOBEC-mediated mutagenesis²⁵ (Supplementary Fig. ^5B). The most complex locus in the BFB lineage (Z43, late BFB, Fig. 5c), contained six variant junctions in cis, including two late tandem duplications (j5, j6). Although j6, which connects the distal portion of 12p to the 12p centromere, was not directly observed in the short-read WGS data, it was imputed (dashed line, j6, Fig. 5c) to resolve the duplication of j1 in clone Z43, as well as two allelic ends in the genome graph. Remarkably, the vast majority (97%) of SNVs detected in this BFB lineage (Fig. 5a) were either shared by all clones or private to a single clone, indicating that these stages of BFB evolution occurred rapidly in the history of the experiment.

We confirmed a chromothripsis event in an independent lineage (Y8), which lacked j1 and all subsequent junctions of the BFB lineage, further supporting the idea that this is an independent lineage (Fig. 5c). Integration of copy number data with the SNV phylogeny showed clones from the unrearranged lineage (Y1 and Y4) and one of the 12p arm loss clones (Y11) to be mutationally distant (>2000 SNVs) from the chromothripsis (Y8) and BFB lineages, which shared over 1583 SNVs (Fig. 5a). Supporting this, a small (~21.5 kbp) simple deletion junction was shared across Y8, Y15, and all the BFB lineage samples, yet was absent in Y11 (Supplementary Fig. ^5C).

This comparison established that the 12p loss in Y11 could not have occurred after j1 and indicates that a second independent arm loss must have given rise to Y15. Interestingly, the Y15 arm loss clone was clustered in the BFB/Y8 clade in the SNV phylogeny, sharing 30 SNVs with the BFB lineage which it did not share with Y8 (Fig. 5a). This indicates that the 12p arm loss in Y15 may have arisen either before or after j1. Although the breakpoints of the Y11 and Y15 arm losses could not be mapped due to their location in the 12 centromeric region, based on the SNV phylogeny, they likely represent distinct events. Taken together, these results support a model whereby at least three lineages independently rearranged a previously wild type 12p during telomere crisis (Fig. 5c). Our data appear to have captured sequential steps in the formation of an increasingly complex BFB-like event. Each of these stages must represent a stabilized allele since the post-crisis lines are clonal, and multiple clones share the same rearrangement junctions (Fig. 5b). This necessarily raises the question as to what caused the on-going instability, and how and where these complex alleles are terminated.

Resolution of BFB cycles in telomere crisis

Analysis of junction-balanced genome graphs allows for the nomination of “loose ends” (or allelic ends), representing copy number changes that cannot be resolved through assembly or mapping of short reads. We identified three distinct loose ends across the four variant graphs spanning the eight clones in the BFB lineage (Fig. 5c and Supplementary Fig. ^5D). Each of these loose ends was placed at the terminus of their respective reconstructed allele, and we posit they represent the new “ends” of the derivative alleles of the BFB lineage. Distinct ends for each of these rearranged lineages suggest the derivative 12p allele could have been stabilized independently. We did not observe telomere repeat-containing reads mated to these loose ends, arguing against neo-telomere formation at these loci. Instead, loose reads represented highly repetitive unmappable sequences which may be a result of the junctions being in close proximity to centromeric regions (see below).

To resolve the genomic architecture at these loci, we generated karyotypes from metaphase spreads for representative rearranged clones (Fig. 6a and Supplementary Fig. ^6A), which revealed that in the BFB and chromothripsis (Y8) lineages, the chromosome 12 derivative was likely linked to a copy of chromosome 21 with an intact long arm (Supplementary Fig. ^6A). These observations were confirmed with chromosome painting, demonstrating a derivative chromosome transitioning between 12 and 21 (Fig. 6b, c and Supplementary Fig. ^6B). Two possible events can explain these findings: the 12–21 fusion could have occurred as an early event during telomere crisis, preceding the divergence of Y8 (chromothriptic) and the BFB lineage; alternatively, independent 21 fusion events stabilized the derivative chromosome 12 following the formation of the distinct junction lineages in Fig. 5c. We consider the first possibility unlikely since the creation of the long-range inversion (j1) and subsequent fold-back junctions in the BFB lineage would require the formation of interstitial 12p breaks on a 12p-21 derivative chromosome. Such breaks are predicted to result in the loss of 21, which would be distal to these junctions on the fusion allele. Furthermore, the acrocentric nature of chromosome 21 would make it more likely to stabilize the overall chromosome architecture, suggesting that an early 12–21 derivative chromosome would be unlikely to engage in the additional SV events observed in the BFB lineage. We, therefore, consider it likely that each of the BFB cycles and chromothripsis clones was independently resolved through subsequent fusion to 21 (Fig. 6c).

Unlike the BFB and chromothripsis clusters, one of the two 12p arm loss lineages (Y11) did not appear to be fused to chromosome 21. In this clone, the derivative chromosome 12 appears to contain a distinct fusion (with a longer p-arm) (Supplementary Fig. ^6C). This is consistent with our analysis of the SNV phylogeny, showing Y11 to be mutationally distant from the BFB lineage (Fig. 5a and Supplementary Fig. ^5C). We were unable to further resolve the nature of the stabilization event in this clone. It would be necessary to perform long molecule DNA sequencing across different lineages in order to confirm the distinct nature of the fusion junction in each of the post-crisis clones.

Immortalized cell line panel

Details of the post-crisis immortalized cell line panel are provided in Supplementary Table ¹. HA-1M cells were a kind gift of Silvia Bacchetti^29,49, SW13/26/39 cells were a kind gift of Jerry Shay²⁸, and Bet-3B/3K and BFT3B/G/K cells were a kind gift of Roger Reddel²⁷.

Cloning and plasmids

A dual-shRNA vector LM2PshRB.698-p21.890-PURO was used to knockdown Rb and p21⁵⁰. The inducible dCas9-VPR (pCW57-dCas9-VPR) construct was created by Gibson assembly of the dCas9-VPR insert from SP-dCas9-VPR (Addgene #63798) into pCW57-MCS1-P2A-MCS2-Neo (Addgene #89180). Retroviral pLVX-hTERT was a kind gift of Teresa Davoli. Activating TERT gRNAs were targeted up to 1000bp upstream of the TERT promoter transcriptional start site and designed using online software from the Broad Institute (portals.broadinstitute.org/gpp/public/analysis-tools/sgrna-design). gRNA sequences were cloned into a modified version of lentiGuide-Puro (Addgene #52963) in which the selection cassette had been swapped for Zeocin resistance. Activating TERT gRNA sequences are shown in Supplementary Table ³. TTN gRNA sequences were used as described³².

Viral gene delivery

Retroviral constructs were transfected into Phoenix amphitropic cells using calcium phosphate precipitation. Lentiviral constructs were transfected with appropriate packaging vectors using calcium phosphate precipitation into 293-FT cells. Viral supernatants were collected and filtered before addition to target cells, supplemented with 4μg/ml polybrene. For activating gRNA constructs, multiple viral supernatants were collected and concentrated using PEG-it Virus Precipitation Solution (System Biosciences LV810A-1). Cells were infected two to three times at 12-h intervals before selection in the appropriate antibiotic.

Immunoblotting

For immunoblotting, cell pellets were directly lysed in 1× Laemmli buffer (2% SDS, 5% β-mercaptoethanol, 10% glycerol, 0.002% bromophenol blue, and 62.5mM Tris-HCl pH 6.8) at a concentration of 10⁷ cells/ml. Lysates were denatured at 100°C, and DNA was sheared with a 28½ gauge insulin needle. Lysates were resolved on SDS/PAGE gels (Life Technologies), transferred to nitrocellulose membranes, and blocked with 5% milk in TBS with 0.1% Tween-20. Primary antibodies (anti-Cas9 7A9-3A3, Cell Signaling Technology #14697S 1:1000, anti-γ-tubulin Sigma #T5326 1:1000, anti- Human Retinoblastoma protein BD Pharmigen #554136 1:500, anti-p21 F-5 Santa Cruz sc-6246 1:200) were incubated overnight, before membrane washing and incubation with appropriate HRP-conjugated secondary antibodies (Amersham, NA934 and NXA931, 1:20,000) and detection with SuperSignal ECL West Pico PLUS chemiluminescence (ThermoFisher).

Immunofluorescence

Cells were grown on glass coverslips and fixed in 3% paraformaldehyde and 2% sucrose. Coverslips were permeabilized in 0.5% Triton-X-100/PBS, and blocked in goat block (0.1% BSA, 3% goat serum, 0.1% Triton-X-100, 2mM EDTA) in PBS. Primary and secondary antibodies (Rabbit anti-53BP1 Abcam #ab-175933 1:1000, F(ab’)2-goat anti-rabbit IgG (H+L) Cross-adsorbed Alexa Fluor 488 ThermoFisher A-11070, 1:500) were diluted in goat block. Slides were counter-stained with DAPI and mounted using prolong gold antifade medium. Images were acquired on a DeltaVision microscope (Applied Precision) equipped with a cooled charge-coupled device camera (DV Elite CMOS Camera), with a PlanApo ×601.42 NA objective (Olympus America), and SoftWoRx software. Images were analyzed for foci numbers using a custom-made algorithm written for FIJI, courtesy of Leonid Timashev⁵¹.

Metaphase spread preparation and staining

Metaphase spreads were prepared by treatment of cells with 0.1µg/ml colcemid (Roche) for 3h before trypsinization and swelling at 37°C for 5–10min in 0.075M KCl. Cells were fixed in a freshly prepared 3:1 mixture of methanol to acetic acid and stored at 4°C overnight or longer. Spreads were prepared by dropping cell solution onto cold glass slides exposed to steam from a 75°C water bath, flooding slide with acetic acid, before exposure of the dropped cells for 3–5s in steam. Slides were dried overnight before storage in 100% ethanol at −20°C. For visualization of fusions, slides were rinsed in PBS, fixed in 4% formaldehyde/PBS for 5min, and dehydrated in an ethanol series before co-denaturation of the slide and PNA probes (TelG-Cy3 PNA Bio F1006, CENPB-AF488 PNA Bio F3004) for 3min at 80°C in hybridization solution (10mM Tris-HCl pH 7.2, 70% formamide, 0.5% Roche 11096176001 blocking reagent). Hybridization was carried out for 2h at RT in the dark, before washing twice in 10mM Tris-HCl pH 7.2, 70% formamide and 0.1% BSA, then washing three times in 0.1M Tris-HCl pH 7.2, 0.15M NaCl, 0.08% Tween-20. DAPI was included in the second wash. Slides were dehydrated through an ethanol series before mounting with Prolong Gold antifade medium (Invitrogen).

For chromosome painting, slides were prepared as above for chromosome fusions, and co-denaturation of chromosome-specific paints (XCP-12 Metasystems D-0312-050-FI, XCP-21 Metasystems D-0321-050-OR) was carried out at 75°C for 2min, before hybridization overnight at 37°C. Post hybridization washes were 0.4× SSC for 2min at 72°C, 2× SSC, 0.05% Tween-20 for 30s, followed by counterstaining in DAPI for 15min, and a rinse in ddH₂O before mounting in Prolong Gold antifade medium (Invitrogen). For karyotyping, slides were prepared as above, and analysis was carried out on DAPI stained chromosomes.

BAC probes

To identify individual chromosomes on metaphase spreads, BAC probes were ordered from BACPAC Genomics (Chr.12p11.2 RP11-90H7, Chr.18q12.3–21.1 RP11-91K12, Chr.6p21.2–21.3 RP11-79J17). Probe DNA was nick-translated with either Digoxigenin-11-UTP or Biotin-16-UTP (Roche) using DNase I (Roche) and DNA polymerase I (NEB) overnight at 15°C. Probes were precipitated with Cot1 human DNA (Invitrogen) and salmon sperm DNA (Invitrogen) and resuspended in 50% formamide, 2× SSC, and 10% dextran sulfate before denaturation for 8min at 80°C. Metaphase spreads were prepared as above, and slides were denatured with 70% formamide, 2× SSC for 2min at 80°C before dehydration through an ethanol series. Slides were co-denatured for 2min at 80°C with TelG-647 (PNA Bio F1014) in hybridization solution (10mM Tris-HCl pH 7.2, 70% formamide, 0.5% Roche 11096176001 blocking reagent) followed by a 2-h hybridization at RT. Denatured BAC probes were then applied and hybridized overnight at 37°C. Slides were washed for 3×5min in 1× SSC at 60°C, followed by a blocking step in 30μg/ml BSA, 4× SSC and 0.1% Tween-20 for 30min at 37°C. BAC probes were detected with anti-digoxigenin–rhodamine (Roche 11207750910, 1:400) and Avadin-FITC antibodies (VWR CAP21221, 1:400) in 10μg/ml BSA, 1× SSC, and 0.1% Tween-20 by incubating for 30min at 37°C, before washing twice for 5min in 4× SSC and 0.1% Tween-20 at 42°C. Counterstaining with DAPI was carried out for 15min at RT, before a further wash at 42°C in 4× SSC, 0.1% Tween-20, and mounting in Prolong Gold antifade (Invitrogen). Images were acquired on the Deltavision microscope equipment detailed above. Images were analyzed using FIJI⁵². Briefly, individual chromosomes were detected with the BAC probes. Intensity measurements of the TelG signal were quantified for each p-arm and q-arm of identified chromosomes. Measurements were also taken for all other telomeres in the same spread. Background subtracted measurements for all telomeres were compared to the shortest (lowest intensity) of the Chr.12p (or Chr.18p or Chr.6p) telomeres on each spread, and the results expressed as a ratio. For identification of fusions containing chromosome 12, the same protocol was carried out using a centromere probe (CENPB-AF488, PNA Bio F3004).

qPCR

RNA was isolated from cell pellets using a Qiagen RNeasy kit, according to the manufacturer’s instructions. cDNA was synthesized using Superscript IV first-strand synthesis (ThermoFisher). qPCR was carried with SYBER Green reagents (ThermoFisher) and run on a Life Technologies QuantStudio 12K machine. qPCR primer sequences are shown in Supplementary Table ³. Expression was quantified using the standard ΔΔCT method relative to β-actin.

TRAP assay

Telomerase activity was assessed using the TRAPeze kit (EMD Millipore S7700) according to the manufacturer’s instructions. Amplification products were resolved on 12% PAGE gels and visualized with EtBr staining.

STELA, fusion PCR, and telomeric blots

High-molecular-weight DNA was extracted from cell pellets using a MagAttract HMW DNA kit (Qiagen) and solubilized by overnight digestion with EcoRI (for STELA and Fusion PCR) or a combination of AluI and MboI (for telomeric blots). STELA was carried out essentially as described³⁴. Briefly, 10ng of DNA was ligated to a mixture of six telorette linkers (Supplementary Table ³) overnight at 35°C, before dilution with water to a concentration of 200pg/μl. Multiple PCR reactions for each sample were carried out with 200pg of annealed DNA using the XpYpE2 and teltail primers (Supplementary Table ³) and FailSafe PCR reagents (Epicenter). PCR conditions were as follows: 94°C for 15s, 27 cycles of 95°C for 15s, 58°C for 20s, 68°C for 10min, and a final extension at 68°C for 9min. PCR products were resolved on 0.8% TAE gels, denatured, and transferred to the Hybond membrane via southern blotting. Products were detected with a randomly primed α-³²P DNA probe created by amplification of the telomere-adjacent region of the XpYp telomere (using XpYpE2 and XpYpB2 primers, Supplementary Table ³). For quantification, FIJI was used to measure the relative signal between indicated molecular weight markers relative to the background signal for each sample.

Fusion PCR was carried out essentially as described^11,35. Subtelomeric primers (Supplementary Table ³) used for amplification of telomeric fusions were XpYpM, 17p6, and 21q4. The control primer XpYpc2tr was included for control amplification of XpYp subtelomeric DNA and detected using EtBr. Fusion products were detected with a random primed α-³²P-labeled (Klenow) DNA probe (21q probe) specific for the TelBam11 telomere subfamily^36,53, which was created with the 21q4 primer and 21q-seq-rev2. The number of fusions per haploid genome (6pg) is calculated based on the amount of input DNA in each PCR reaction.

Telomere length was assessed using telomeric restriction fragment analysis. Briefly, AluI/MboI-digested genomic DNA was run on 0.8% TAE gels, before denaturation, neutralization, and transfer onto a Hybond membrane according to standard Southern blotting procedures. Telomeric DNA was detected using a TTAGGG repeat primed α-³²P (Klenow)-labeled Sty11 telomeric repeat probe⁴⁵.

WGS library preparation

Genomic DNA was extracted from cell pellets using a QIAGEN QIAamp DNA mini kit and sheared using a Covaris Ultrasonicator (E220) to ~300bp fragments. DNA concentration was measured using Qbit 4.0 reagents (ThermoFisher), and 200ng of fragmented DNA was used for library preparation. End repair and A-tailing was carried out with NEBNext End repair reaction enzyme mix and buffer (E7442), and KAPA dual-indexed adapters (Roche) were ligated using the T4 DNA ligase kit from NEB (M0202). Post-ligation size selection was performed with AMPure XP beads (Beckman Coulter) before washing two times in 80% ethanol. Libraries were amplified using KAPA HiFi HotStart ready mix (Roche) and P5 and P7 primers (IDT). PCR program was as follows: 98°C for 45s, five cycles of 98°C for 15s, 60°C for 30s, 72°C for 30s, and a final extension at 72°C for 5min. A further size selection and washing step was carried out after library amplification, and library quality was confirmed on Bioanalyzer chips (Agilent) and using a KAPA Library Quantification kit (Roche). Libraries were pooled and submitted for sequencing on NovaSeq 6000 at the New York Genome Center.

WGS basic data processing

Reads were aligned to GRCh37/hg19 using the Burroughs-Wheeler aligner (bwa mem v0.7.8, http://bio-bwa.sourceforge.net/)⁵⁴. Best practices for post-alignment data processing were followed through use of Picard (https://broadinstitute.github.io/picard/) tools to mark duplicates, the GATK (v.2.7.4) (https://software.broadinstitute.org/gatk/) IndelRealigner module, and GATK base quality recalibration.

Variant rearrangement junctions were identified using SvABA³⁰ (https://github.com/walaj/SvABA) and GRIDSS³¹ (https://github.com/PapenfussLab/gridss) with standard settings. For MRC5 samples, the somatic variant setting of each tool was used, with the ancestral MRC5 line as the matched normal. For post-crisis SV40 transformed cell lines the respective pre-crisis clone was used as the matched normal. 1-kbp binned read depth was computed and corrected for GC and mappability using fragCounter (https://github.com/mskilab/fragcounter). Systematic read-depth bias was subsequently removed using dryclean (https://github.com/mskilab/dryclean)⁵⁵.

Low-pass WGS clustering

Genome-wide binned read depth was aggregated across 118 low-pass WGS clones across 10-kbp bins by taking the median of 1-kbp binned normalized read depth from dryclean (see above). To minimize read-depth noise in unmappable regions, recurrent (>10% of the cohort) low-quality coverage regions (defined by Hadi et al.³) were combined with regions bearing consistently high variance in our high-pass sequencing dataset (standard deviation >0.3 for bin value over the mean in 100-kbp windows). Hierarchical clustering was then applied on the genome-wide Euclidean distance of bins, with “method=ward.D2” option. Six clusters were identified following dendrogram inspection.

Junction balance analysis

Preliminary junction-balanced genome graphs were generated for MRC5 and SV40T cell lines from binned read depth and junction calls (see above) using JaBbA (https://github.com/mskilab/JaBbA)³. Briefly, 1-kbp binned read-depth output from dryclean was collapsed to 5kbp and JaBbA was run with slack penalty 500 for MRC5 clones and 100 for SV40T cell lines. gGnome (https://github.com/mskilab/gGnome) was used to identify complex structural variant patterns. Genome graphs and corresponding genomic data (e.g., binned coverage, allelic bin counts) were visualized using gTrack (https://github.com/mskilab/gTrack).

Joint inference of junction balance in MRC5

To chart structural variant evolution across sub-clades of MRC5 clones, a procedure was developed to jointly infer junction-balanced genome graphs in a lineage (e.g., BFB lineage in Fig. 5c). This co-calling algorithm augmented the existing JaBbA model, described in detail by Hadi et al.³, enabling the application to a compendium of genome graphs by minimizing the total number of unique loose ends assigned a nonzero copy number across the graph compendium.

To describe this algorithm, we extend the notation introduced in Hadi et al.³. Formally, we define a collection {Gⁱ}_i  1,…,_n of identical genome graphs across n clones, each a replica of a “prototype” genome graph G⁰. The mapping p maps each vertex v  V(Gⁱ) and edge e  E(Gⁱ), i  1,…,n to its corresponding vertex p(v)  V(G⁰) and edge p(e)  E(G⁰) in the prototype graph. We then jointly infer unique copy number assignments κⁱ to the vertices and edges of each genome graph Gⁱ by solving the mixed-integer program:

where xⁱ and $J^i$ represent the binned read-depth data and bin-node mappings for clone i and $\mathcal{V}\left(G^i,{\kappa }^i,x^i,J^i\right)$ is the read-depth residual for genome graph i, analogous to Hadi et al.³. An additional term in this new joint formulation is uⁱ:E(Gⁱ)→{0,∞}, which is a data-derived mapping that constrains the upper bound of each edge e  E(Gⁱ), e.g., on the basis of whether that junction has read support in clone i. In addition, a joint complexity penalty $\mathcal{R}$ couples the collection {κⁱ}_i  1,…,_n of copy number mappings across the collection of graphs {Gⁱ}_i  1,…,_n to each other by jointly penalizing loose ends at all vertices that map to the same prototype graph vertex v  V(G⁰). Formally,

λ in Eq. ⁽¹⁾ controls the relative contribution of the read-depth residual and complexity penalty to the objective function. It is important to note that while each of the graphs Gⁱ have an identical structure, the constraints imposed by the upper bounds uⁱ and bin profiles xⁱ couple each graph to its junction and read-depth data, and hence lead to a unique fit κⁱ on the basis of this data. The ${\ell }_0$ penalty (defined using the Iversion bracket〚〛operator) in Eq. ⁽²) couples the solutions κⁱ by adding an exponential prior on the number of unique loose ends across the entire graph compendium, where uniqueness is defined by the mapping p to the prototype graph G⁰.

This joint mixed-integer programming model in Eq. ⁽¹) is implemented in the “balance” function of gGnome. The model was applied to a collection of genome graphs representing the structure of chromosome 12 across 13 clones. The prototype graph for this genome graph collection was built from the disjoint union of intervals of the 13 preliminary graphs (via the GenomicRanges “disjoin” function) and the union of junction calls fit across those graphs (via gGnome “merge.Junction” function). Each graph was associated using the read-depth data and bin-to-node mappings as per Hadi et al.³. The mapping uⁱ for each reference edge was set to ∞ while variant edges were assigned ∞ on the basis of bwa mem realignments of read pairs in each clone.bam file to the corresponding junction contig via rSeqLib (https://github.com/mskilab/rSeqLib)⁵⁶, otherwise they were assigned 0.

Equation ⁽¹) was then solved using the IBM CPLEX (v12.6.2) MIQP optimizer within the gGnome package after setting the hyperparameter λ to 100. This value was chosen after a parameter sweep observing for the visual concordance of genome graphs, loose ends, and read-depth profiles in the region.

Joint reconstruction of allelic evolution in MRC5

Evolving 12p alleles were jointly reconstructed across 13 MRC5 clones through the analysis of junction-balanced genome graphs (Gⁱ,κⁱ) (see “Joint inference of junction balance in MRC5” section above). The procedure for joint allelic phasing described in Hadi et al.³ was extended to identify the most parsimonious collection of linear and/or cyclic walks and associated walk copy numbers that summed to the vertex and edge copy numbers in the compendium (Gⁱ,κⁱ).

Formally, the subgraph of vertices and edges with a nonzero copy number in each (Gⁱ,κⁱ) were exhaustively traversed to derive all minimal paths and cycles Hⁱ, where for each walk h  Hⁱ maps to subsets V(h)  V(Gⁱ) and E(h)  V(Gⁱ) of vertices and edges in the graph Gⁱ. The nodes and vertices of these walks were then projected via the mapping p to define a unique set of walks H⁰ in the prototype graph G⁰. We extend our notation p (see the previous section) so that for a walk h  Hⁱ the mapping p(h)  H⁰ denotes the walk formed by projecting the vertices and edges of h via p to H⁰. With these definitions, the single graph haplotype inference defined in Hadi et al.³ was extended to a joint inference by solving the following mixed-integer linear program to assign a copy number ϕⁱ(h)  $\mathbb{N}$ to each walk h  Hⁱ.

where the function δ(v,h) and δ(e,h) is 1 if vertex v and edge e belong to walk h and 0 otherwise. The Iverson bracket (〚〛) operator in the objective function Eq. ⁽³) minimizes the total number of unique walks used across the compendium, hence identifying a jointly parsimonious assignment of copy number to walks across the compendium of graphs. Equation ⁽³) was solved using the IBM CPLEX (v12.6.2) MIQP optimizer within the gGnome package. Variant cycles and paths from the resulting solution were manually combined to yield a set of consistent linear paths, i.e., somatic haplotypes, to yield allelic reconstructions in Fig. 5c.

Loose-end classification

Each loose end in each MRC5 genome graph was analyzed to identify a clone-specific (i.e., absent in the ancestral MRC5 line) origin for the mates of high mapping quality (MAPQ=60) reads mapping to the location and strand of the loose end. These mates were assessed for neo-telomeric sequences by counting instances of 11 permutations of a 12-bp telomere repeat motif (TTAGGGTTAGGG) (using the R/Bioconductor Biostrings package) in the mates. The mates were also assembled into contigs using fermi⁵⁷ aligned using bwa mem⁵⁸ via the RSeqLib R package⁵⁶ to hg38 (which contains a more highly resolved centromere build) to characterize novel repeat (e.g., centromere) fusions (GATK Human reference genome, hg38, data bundle including Homo_sapiens_assembly38.fasta, gs://gcp-public-data--broad-reference). The loose-end loci were also assessed through overlap with the hg19 repeatMasker database (human_g1k_v37_decoy.repeatmasker) for the presence of reference annotated repeats that might explain the absence of a mappable junction explaining the copy number change.

SNV phylogeny

To compute an SNV phylogeny across MRC clones, we first identified SNV that were acquired in MRC5 clones relative to the ancestral MRC5 line using Strelka2⁵⁹ (https://github.com/Illumina/strelka) under paired (i.e., tumor/normal) mode with the clone as the “tumor” and the MRC ancestral line as the “normal” sample and default parameters and GATK hg19 resource bundle (Genome Analysis Toolkit GATK Resource Bundle for hg19; gs://gatk-legacy-bundles). Acquired SNVs were first filtered according to the Strelka2 PASS filter as well as additional filters (MQ=60, SomaticEVS>12, total ALT count >4) yielding 27,220 total unique variants across the 13 MRC5 clones. Reference and variant allelic read counts were assessed at each SNV site (via the R/Bioconductor Rsamtools package, version 3.6.1, http://www.r-project.org/) across all 13 clones. We then further required a >0.5 posterior probability of a variant being present in a sample, by assuming the Binomial likelihood of variant read count and using the aggregated allele frequency in all samples as the prior, resulting in the final 14,970 unique mutations. The binary matrix of clones by SNV loci was then used to derive a neighbor-joining phylogenetic tree using the R/Bioconductor package ape. Following tree construction, we associated each SNV with its most likely phylogenetic tree branch by comparing the binary incidence vector associated with each SNV with the binary incidence vector associated with each tree branch, and finding the closest branch using Jaccard distance, only linking SNV to branches when the SNV was within <0.1 Jaccard distance of the closest branch, thus producing the groupings of SNVs in Fig. 5a.

Parental SNP allelic phasing and imbalance

Germline heterozygous sites in the parental MRC5 line were identified by computing allelic counts at HapMap sites (GATK human reference genome, hg19 data bundle, hapmap_3.3.b37.vcf) and identifying loci with variant allele fraction >0.3 and <0.7. Y11, a clone with loss of a single allele at 12p, was chosen to phase parental SNPs on 12p. At each locus, the allele (reference or alternate) with a 0 read count was assigned to the “L” (lost) haplotype and the other allele was assigned to the “R” (retained) haplotype. (All heterozygous SNP loci in the region contained exactly one allele with a 0 read count). L and R allelic counts were then computed at these sites across all 13 high-pass WGS and 131 low-pass WGS samples. These counts were divided by the genome-wide mean of heterozygous SNP allele counts (in these 100% pure and nearly diploid samples) to derive the absolute allelic copy number⁶⁰.

SNV clustering

Inter-SNV distances were computed for all pairs of reference adjacent acquired SNVs associated with each MRC5 clone and visualized as rainfall plots. Runs of two or more SNVs with inter-SNV distances <2kbp were nominated as clusters. Two distinct SNV clusters were identified on chromosome 12p across the 13 clones.

Statistical analysis

Statistical analysis for in vitro experiments was carried out using Prism software (GraphPad Software). All relevant statistical experimental details (n numbers, SD) are provided in the figure legends. Statistically significant associations between binary variables were determined using a two-tailed Fisher’s exact test. Significance was assessed on the basis of Bonferroni-corrected P values <0.05. Effect sizes (odds ratios) are reported alongside 95% confidence intervals for each test. Details of all other quantitative analyses (e.g., read-depth processing, clone clustering, genome graph inference, allelic reconstructions, phylogenetic reconstruction, SNV clustering, parental SNP phasing) are described above.

We thank Jerry Shay (UTSW), Silvia Bacchetti, AdVec, and Roger Reddel (CMRI) for the generous provision of cell lines used in this study. We also thank Keshav Sharma for assistance with cell culture maintenance. S.M.D is funded jointly by a NIH/NCI grant (5R35CA210036) and Melanoma Research Alliance grant (577521) both awarded to T.d.L. T.d.L. has additional funding support from NIH/NIA (5R01AG016642-21A1), Starr Cancer Consortium (I9-A9-047), Breast Cancer Research Foundation and a Glenn Foundation Award from the Glenn Foundation. M.I., X.Y., J.B., and H.T. are supported by Burroughs Wellcome Fund Career Award for Medical Scientists, Doris Duke Clinical Foundation Clinical Scientist Development Award, Starr Cancer Consortium Award, and National Institutes of Health (NIH) U24-CA15020 to M.I., as well as Weill Cornell Medicine Department of Pathology Laboratory Medicine startup funds.