Rates and spectra of de novo structural mutations in Chlamydomonas reinhardtii

SM detection

We performed PacBio HiFi or continuous long-read (CLR) sequencing of MA lines of the C. reinhardtii strains CC-1952 (N = 4 MA lines, all HiFi) and CC-2931 (N = 8, six HiFi and two CLR). These geographically distinct strains were selected on the basis of their relatively low (CC-1952, overall μ = 4.05 × 10⁻¹⁰ per site per generation) and high (CC-2931, μ = 15.6 × 10⁻¹⁰) SNM and indel rates in the MA lines (Ness et al. 2015), which had been maintained for approximately 1050 generations by single-cell descent (Morgan et al. 2014). To perform strain-specific detection of SMs, we first produced near-complete reference assemblies for the two ancestral strains. The 17 C. reinhardtii chromosomes were assembled into totals of 50 and 39 contigs with N50s of 4.25 and 3.81 Mb, for CC-1952 and CC-2931, respectively (Supplemental Fig. S1; Supplemental Table S1). We subsequently defined ∼98% of each ∼111 Mb ancestor genome as “callable” (i.e., sites where SMs could be called with high confidence; see Methods). This represents a substantial increase over the ∼71% obtained in our previous study using short-read technology (Supplemental Fig. S2; Ness et al. 2015). We sequenced MA lines at sufficient depth of coverage to produce highly contiguous assemblies (Supplemental Fig. S1; Supplemental Table S1), enabling us to call SMs using three approaches: directly from MA line PacBio read alignments against the appropriate ancestral reference using Sniffles (Sedlazeck et al. 2018), from MA line PacBio assembly alignments against the reference using MUM&Co (O'Donnell and Fischer 2020), and from Cactus pangenomes using vg (Garrison et al. 2018; Armstrong et al. 2020). Sniffles and MUM&Co were run individually for each MA line, whereas all MA lines and the ancestor for each strain were analyzed collectively with vg from a single Cactus alignment. We subsequently collated and manually curated all variant calls using a combination of read visualization and mapping approaches (Fig. 1A).

We classified our curated data set of SMs (a total of 120 in CC-1952 and 443 in CC-2931; Supplemental Dataset S1) into eight categories: expansions and contractions of tandemly repeated sequence (e.g., in microsatellites or satellite DNA, collectively termed tandem repeat mutations [TRMs]), duplications, deletions, insertions and excisions of mobile elements, and inversions and translocations (Fig. 1B). The different callers varied substantially in their ability to identify different SM types (Fig. 1C). Only 19.2% of SMs were called by all three tools, highlighting the importance of combining approaches. vg was most successful overall, calling 79.4% of SMs and 22.4% uniquely, although it called only 20.0% of inversions and translocations. MUM&Co called only 16.7% of TRMs, although it identified 64.5% of SMs in other categories. Sniffles called 47.2% of SMs and was generally outperformed by the assembly-based approaches, although it was superior in identifying duplications (calling 62.7% of duplications and 39.0% uniquely). The assembly-based methods failed to call duplications >30 kb (i.e., longer than the reads) because they were generally collapsed and absent in the assemblies. We also called 7.3% of SMs manually from alignment files and read visualization, most of which were TE insertions located at the breakpoints of other complex SMs (e.g., inversions and translocations).

The three approaches also differed markedly in the proportions of rejected calls (i.e., variants that could not be classified as genuine SMs), and all three returned more rejected calls than confirmed SMs (reaching 3× as many for vg) (Supplemental Fig. S3). Rejected calls generally fell in three categories: calls made within regions that we had defined as uncallable, which were not considered further; MUM&Co or vg calls that received assembly support but not read support (i.e., “valid” calls introduced by assembly errors); or calls that received support from neither assemblies nor reads (unsupported variants). Considering the ratio of rejected calls to confirmed SMs in CC-2931, Sniffles performed best (1.14), followed by MUM&Co (1.81) and vg (3.42). However, these ratios varied considerably depending on the type and sequence context of the rejected call. All categories of rejected calls were associated with repetitive sequence (Supplemental Fig. S4). Uncallable regions correspond to the most repetitive regions of the genome, generally long satellite arrays, including some centromeric and subtelomeric regions (Supplemental Dataset S2). Although they contain only ∼2% of sites, a substantial proportion of all variant calls were made in these regions (26.4%). Although we expect that the majority of calls in the uncallable regions were false positives, the enrichment of these regions for tandem repeats may have led to an underestimation in the rates of expansions and contractions in our analysis. Considering callable regions, tandem repeats of all lengths were a major source of both assembly errors and unsupported variants. After excluding all tandem repeats, the ratio of unsupported variants to confirmed SMs improved for all callers, falling from 0.21 to 0.17 for MUM&Co and from 0.91 to 0.20 for vg, although substantial proportions of calls attributed to assembly errors remained (Supplemental Fig. S3). Consequently, although our results show that aligners and callers are capable of detecting SMs in tandemly repeated regions (expansions and contractions) (Fig. 1C), attempting to do so risks introducing many false positives in fully automated pipelines. Rejected calls are further discussed in the Supplemental Material.

Rates and spectra of SMs

The rates and spectra of SMs were markedly different between the two strains (Fig. 2A,B). CC-2931 MA lines experienced ∼85% more SMs than CC-1952 MA lines, and overall SM rates were significantly different between the strains (μ_{SM (CC-1952)} = 2.58 × 10⁻¹⁰ and μ_{SM (CC-2931)} = 4.30 × 10⁻¹⁰ per site per generation; W-test, P = 4 × 10⁻³). However, the within-strain variance in μ_SM among MA lines was ∼12% higher than between them (ANOVA test, F = 10.8, P = 8 × 10⁻³). In terms of the number of bases affected, CC-2931 MA lines experienced larger SMs than CC-1952 MA lines and also experienced a greater variety of SM types (Fig. 2A). We observed only three SMs >20 kb in CC-1952 MA lines, whereas almost all large chromosomal rearrangements were found in CC-2931; that is, there were 1.75 inversions (median ∼243 kb) and 2.50 translocations per CC-2931 MA line (collectively 4.05 × 10⁻³ rearrangements per genome per generation) (Supplemental Fig. S5), compared with a single 5.2-kb inversion in the CC-1952 MA lines. Deletions were also rare in CC-1952 (less than one per MA line on average), as were mobile excisions, explained by a relatively low frequency of active cut-and-paste DNA transposons (see below).

Although there were clear differences between the strains, some SM properties were shared by CC-1952 and CC-2931. TRMs were the most common category of SMs <3 kb in length in both strains (representing 60.8% and 35.0% of SMs in CC-1952 and CC-2931, respectively) and occurred at similar rates in the two strains (median μ_{TRM (CC-1952)} = 1.68 × 10⁻¹⁰, μ_{TRM (CC-2931)} = 1.60 × 10⁻¹⁰; W-test, P = 0.93). Mobile insertions dominated the spectra of SMs >3 kb in length (representing 21.7% and 37.9% of all SMs in CC-1952 and CC-2931, respectively), although the median rate of mobile insertions was ∼2.5× higher in CC-2931 (14.08 × 10⁻¹¹ per site and generation) than in CC-1952 (5.15 × 10⁻¹¹). Duplications were also relatively common in both strains (12.5% and 9.9% of all SMs in CC-1952 and CC-2931, respectively), although duplications in CC-1952 (median length 0.9 kb) were significantly shorter than in CC-2931 (median length 12.9 kb; W-test, P = 3 × 10⁻⁶). Duplications were more frequent than deletions in both strains (W-test, P < 1.5 × 10⁻²) and also significantly longer in CC-2931 (median length, 12.9 kb vs. 1.7 kb, respectively; W-test, P = 1.81 × 10⁻²). Across all SM types, we found a bias toward increasing genome size: CC-2931 MA lines gained 171.7 bp, and CC-1952 MA lines gained 50.4 bp, per generation, on average. After excluding TEs, the rate of genome size expansion fell by approximately half in CC-2931 and approached zero in CC-1952, highlighting the role of TEs in driving genome size expansion in C. reinhardtii MA lines.

In addition to identifying SMs, we used the PacBio HiFi reads to call SNMs and indels <50 bp in length. This analysis suggested that SMs represent ∼6% of the total mutations in CC-1952 and ∼12% in CC-2931. However, most mutations affecting short lengths of sequence were called in tandem repeats (∼60%). Given the uncertainty of variant calling in these sequences, we recalculated mutation rates, excluding these regions. This had a minor effect on the relative contributions of SMs to the total rate (Fig. 2C). We did not find a significant correlation between μ_SM and either SNM or indel mutation rates among MA lines in either strain. Disregarding covariance terms between SM, SNM, and indel mutation rates, variance in μ_SM explained 4.3% of total variance attributed to the total mutation rate in CC-1952 and 9.7% in CC-2931. Hence, our results suggest that SMs, including those that could have functional consequences, occur at a rate approximately 10-fold lower than the rate of SNMs and indels combined, yet their rates of occurrence are highly variable between strains. The broader improvements to short mutation detection brought by PacBio HiFi sequencing are briefly discussed in the Supplemental Material.

Genomic distribution of SMs

To explore the distribution of SMs across the genome and, in particular, their distribution relative to functional sequences, we generated ∼20 Gb of RNA-seq data for the CC-2931 ancestor and annotated coding sequences, introns, and untranslated regions (UTRs). Because >70% of genes are separated by <1 kb in the compact C. reinhardtii genome, we divided intergenic regions into “proximal” (within 500 bp of a gene) and “distal” sequences (>500 bp from a gene), the latter category largely capturing long, highly repetitive intergenic regions, including the centromeres and subtelomeres. We then used two approaches to compare the observed distribution of SM coordinates to the expected distribution based on random sampling. First, we considered the entire span of each SM. Second, we considered the specific coordinates of the SM breakpoints. Although the former approach takes into account the wider possible effects of mutation (e.g., the duplication of entire genes), the latter approach explores whether the breakpoints of SMs were enriched relative to any particular genomics features.

When considering the whole length of duplications, deletions, and inversions, these mutations intersected with genomic annotations, as expected under a random distribution (Fig. 3A). As a result, their proportions of overlap with each functional annotation resembled the actual proportion of the corresponding annotation in the genome. This included coding sequences, suggesting that many of these SMs could have large fitness effects. Furthermore, the breakpoints of duplications and deletions were also distributed randomly (Fig. 3B). These results are consistent with a near absence of selection in our MA experiment, as expected for populations subject to regular bottlenecks. As detailed below, mobile insertions were underrepresented in coding sequences and introns and were enriched in 5′ UTRs and intergenic sequences (Fig. 3B). The coordinates of inversion and translocation breakpoints had qualitatively similar distributions, a result that we attribute to the association of these rearrangements with TEs (see below).

Active TEs

With the exception of one putative mobile satellite in CC-2931, all mobile elements were TEs. We found that the two strains differed markedly in the number and diversity of active TE families. There were 12 active TE families from seven subclasses in CC-2931 and only three active families from two subclasses in CC-1952. We also observed considerable heterogeneity in insertion rates among MA lines and among TE families (Figs. 2B, 4A).

The most active retrotransposons in CC-2931 were an autonomous (Chlamys-9_cRei) and nonautonomous (Chlamys-N4_cRei) pair of Penelope-like elements (PLEs), which generally caused very short insertions of median length 128 bp as a consequence of 5′ truncation. All Chlamys insertions were intronic, a pattern that is broadly consistent with the underlying distribution of their (C)_n microsatellite target (Fig. 4B,C; Craig et al. 2021b). The remaining CC-2931 TEs comprised various types of DNA transposons. The cut-and-paste DD(E/D) transposons EnSpm-3_cRei (autonomous) and EnSpm-N3_cRei (nonautonomous) were active in all lines. EnSpm insertions were significantly enriched in intergenic regions distant from genes, although unlike other TE families, they were not underrepresented in coding sequences (Fig. 4C). Consistent with a higher efficiency of nonautonomous transposons relative to their longer autonomous counterparts (Han et al. 2013), we observed a net increase in copies of EnSpm-N3_cRei (i.e., three ancestral copies compared with a mean of 3.6 copies in the MA lines) and a net decrease in copies of EnSpm-3_cRei (i.e., three ancestral copies compared to a mean of 2.4 copies in the MA lines). The increase in EnSpm-N3_cRei copies is presumably a consequence of transposition occurring during DNA replication (Feschotte and Pritham 2007). It should be noted that the estimated number of EnSpm insertions is a minimum estimate (Fig. 4A) because several cut-and-paste transpositions may have occurred during MA that were not captured.

We observed 12 insertions of copy-and-paste Helitrons, including the particularly long 20.4-kb autonomous Helitron2-7_cRei (Fig. 4B). We also observed an unusual pair of copy-and-paste transposons, the autonomous Replitron-1 and nonautonomous Replitron-N1, which are the founding elements of a new group of eukaryotic transposons named Replitrons that will be described elsewhere (Craig 2022). Like Helitrons, Replitron-1 encodes an HUH endonuclease of the Rep class, although it features only one catalytic tyrosine (i.e., Y1 rather than Y2) in the Rep domain and does not feature a C-terminal helicase domain (Fig. 4B). The Replitrons inserted upstream of “RG” target sequences (R indicating purine), causing variable length target site duplications. Replitron-N1 was the second most active TE, causing insertions in all lines and a maximum of 22 insertions in L13 (Fig. 3B). Replitron-1 and Replitron-N1 insertions were significantly underrepresented within coding sequence and introns (Fig. 4C).

We observed two families of autonomous Cryptons, including CryptonF-1_cRei, which was the most active TE in CC-2931. To our knowledge, these are the first observations of active Cryptons. Goodwin et al. (2003) proposed a model of Crypton insertion via site-specific recombination between a short donor sequence at one terminus and a near-identical target sequence at the integration site, catalyzed by the Crypton-encoded tyrosine recombinase. Our data were consistent with this model, because CryptonF-1_cRei terminates in the motif “CACCG” and targeted “CAYCG” (Y indicating pyrimidine) (Fig. 4B). However, it was also proposed that Cryptons would undergo excision. Although we observed clean Crypton excisions that left behind only the “CAYCG” target, Cryptons increased in copy number, suggesting a complex mode of transposition. CryptonF-1_cRei insertions were enriched at gene-proximal intergenic sequences and 5′ UTRs.

In contrast to CC-2931, only the LINE retrotransposon Dualen-4b_cRei was active in multiple CC-1952 MA lines, and the number of insertions ranged from four insertions in MA line 1 (L1) to nine in L15 (Fig. 4A). This family was also active in CC-2931, but it only caused one insertion in L11 (not shown in Fig. 4A). We also observed two cut-and-paste DD(E/D) transposons in CC-1952 MA lines: a single excision and a single insertion of a P element (P-2_cRei) in L1, and an excision of a giant single-copy 31.7-kb Zisupton DNA transposon (Zisupton-3_cRei) in L3, resulting in its extinction from the genome.

Incidentally, we did not observe any gene retrocopy events in either CC-1952 or CC-2931, presumably because of an absence of active polyadenylated LINEs (e.g., L1 elements in mammals) (Kaessmann et al. 2009) and autonomous LTRs (Tan et al. 2016). Although LTR elements are abundant in the genomes of both strains, we also observed no “deletions” of LTR elements (i.e., formation of solo LTRs by homologous recombination).

TE-mediated SMs

We found that TE insertion and excision events were associated with many other SM types. The number of TE insertions was positively correlated with the combined number of deletions, duplications, inversions, and translocations across all MA lines (Pearson's product-moment correlation, t₁₀ = 2.64, r = 0.64, P = 0.02). More specifically, 70.6% of inversions and translocations in the CC-2931 MA lines involved either CryptonF-1_cRei or Replitron elements at one or both of the breakpoints (Table 1). One exception to this pattern involved truncated copies of both Replitron-N1 and CryptonF-1_cRei at a single breakpoint. We also attributed a far smaller number of SMs to homology-mediated double-strand break (DSB) repair mechanisms (Table 2). These are described in more detail below.

Translocations mediated by Replitron-1 and Replitron-N1 generally featured breakpoints coinciding precisely with the right end of the transposons and were frequently repaired so that one of the derived chromosomes had two elements in a tail-to-tail arrangement. For example, this outcome was observed in a reciprocal translocation between Chromosomes 14 and 17 in CC-2931 L9 (Fig. 5A). Although we do not yet understand the transposition mechanism of these newly discovered TEs, the presence of target site duplications suggests that they may cause DSBs, and their simultaneous insertion at different genomic regions could lead to aberrant repair and rearrangements. CryptonF-1_cRei-mediated rearrangements were associated with both insertions and excisions. For example, CC-2931 L11 experienced a reciprocal translocation between Chromosomes 9 and 11, involving insertions at each breakpoint (Fig. 5B). Other events involved only one CryptonF-1_cRei insertion, although in all cases we observed a “CAYCG” target site at the breakpoint without insertion. Notably, the high insertion rates of both CryptonF-1_cRei and Replitrons in CC-2931 L13 (Fig. 4A) resulted in a highly derived karyotype in this line (Supplemental Fig. S6A). This included a cluster of nonreciprocal translocations involving four chromosomes, mediated by three Replitron-N1 and one Replitron-1 insertions that may have occurred simultaneously.

CryptonF-1_cRei was associated with six inversions on Chromosome 16 in different MA lines (Fig. 5C). In the CC-2931 ancestor, Chromosome 16 features two CryptonF-1_cRei copies in opposite orientation separated by ∼233 kb. However, we did not observe this ancestral state in any of the eight MA lines. The ∼233-kb region was inverted in four MA lines, with either none, one, or both of the copies excised. In L2 and L14, longer inversions were mediated between one of the CryptonF-1_cRei copies and independent “CAYCG” target sites elsewhere on the chromosome. Overall, there were seven different mutated states relative to the ancestor (Fig. 5C), evidencing the hypermutability of Crypton activity. The presence of target sites at each breakpoint in both inversions and translocations may suggest a role for the site-specific recombination activity of the Crypton tyrosine recombinase enzyme in mediating rearrangements.

Finally, we observed duplications associated with the LINE Dualen-4b_cRei, which was most active in CC-1952 (Fig. 4A) but also caused a single insertion and associated duplication in one CC-2931 MA line. More than half of the Dualen-4b_cRei insertions were associated with duplications >50 bp, which represented 87% of duplications in CC-1952 (Table 1). These duplications (of median length 900 bp) (Supplemental Fig. S7) resembled the variable length target site duplications that flank insertions of non-LTR elements (LINEs and PLEs); that is, the duplicated sequence flanked either side of the Dualen insertions. Such target site duplications are caused by resolution of the DNA nicks introduced during insertion, the distance between cleavage sites corresponding to the target site duplication length. These putative Dualen-4b_cRei target site duplications are considerably longer than other large target site duplications reported previously, that is, the 126-bp target site duplications observed in R9 LINEs of rotifers (Gladyshev and Arkhipova 2009). Active Dualens have not previously been observed, and exceptionally long target site duplications are possibly mediated by the dual action of RLE and APE endonucleases (Fig. 4B), which are uniquely present together in the Dualen clade (Kojima and Fujiwara 2005).

Homology-mediated SMs

We next attempted to identify homology-based DSB repair mechanisms that may have mediated the remaining SMs. These included most deletions and duplications in CC-2931 (Table 2). Many mechanisms associated with SMs involve DSB repair, which typically proceeds via two distinct pathways: homologous recombination and canonical nonhomologous end joining. Homologous recombination can induce SMs via recombination between interspersed paralogous sequences, that is, nonallelic homologous recombination. Nonhomologous end joining can also mediate SMs. For example, two DSBs that are present on the same chromosome could be repaired aberrantly, yielding a deletion or an inversion, or two DSBs that are present on different chromosomes could yield a translocation. Several other DSB repair pathways that involve varying lengths of homology tracts also exist (So et al. 2017). For example, microhomology-mediated end joining requires ∼1–16 bp of sequence homology, and single-strand annealing requires macrohomology >30 bp. Both mechanisms generate deletions at single DSBs because they involve DNA end resection, whereby homology is revealed at DNA ends when degraded to single-strand sequences, and this is followed by deletion of any sequence overhanging the homology tract (Sfeir and Symington 2015).

We found evidence of macrohomology (>30 bp) in only ∼5% of duplication events, one event in CC-1952 and two in CC-2931 (Table 2; Supplemental Fig. S8). These duplications could potentially have been caused by nonallelic homologous recombination. In all these cases, the paralogous sequence did not feature TEs. We found no evidence of macrohomology-mediated deletions, suggesting little role for the single-strand annealing repair pathway causing SMs. In contrast, we identified microhomologies in 28% of deletions, four in CC-2931 and one in CC-1952 (Table 2; Supplemental Fig. S9). One of these deletions (in CC-2931 L6), was flanked by several clustered SNMs and indels (Supplemental Fig. S10), a phenomenon that has been observed flanking DSBs repaired by microhomology-mediated end joining (Sinha et al. 2017). Although similar hypermutability was not observed flanking the other microhomology-associated deletions, these deletions were generally shorter than deletions showing no evidence of homology (median length, 225 bp vs. 2,945 bp; W-test, P = 3.77 × 10⁻²). This is consistent with the involvement of DNA end resection in microhomology-mediated end joining.

We also found three duplications in CC-2931 that had five, three, and two nucleotides of sequence homology between their respective breakpoints. Although these duplications are consistent with microhomology-mediated mechanisms, the presence of such short stretches of homology might also be explained by chance. Finally, we found no evidence of either micro- or macrohomology at the breakpoints of the inversions and translocations not associated with active TEs. However, we did observe short deletions at several inversion and translocation breakpoints, which is consistent with repair via nonhomologous end joining. Notably, 40% of non-TE-associated inversions and translocations occurred in a single CC-2931 line, L15 (one inversion and three translocations) (Supplemental Fig. S6B).

Tandem repeat mutations

Tandemly repeated sequences are known to be hypermutable and evolve via mechanisms that include replication slippage and unequal exchange (Lower et al. 2018). As mentioned above, unlike most other SMs, TRMs (grouping expansions and contractions >50 bp) occurred at similar rates in the CC-2931 and CC-1952 MA lines. This could be owing to the independent underlying mechanisms of mutation compared with other SMs. Approximately 22% of TRMs occurred in microsatellites (defined as tandem repeats with monomers <10 bp), and most of the remainder occurred in satellite DNA. The C. reinhardtii genome contains many introns featuring tandem repeats (Craig et al. 2022), and 30.3% and 64.7% of satellite and microsatellite TRMs, respectively, were located in genic sequences. It is important to note that the rate of TRMs per site per generation is several times higher than the genome-wide rate displayed in Figure 2B, because TRMs can only occur within the ∼12% of the genome that features tandem repeats, effectively resulting in a smaller “callable” genome.

We also observed numerous expansions and contractions in the centromeric and subtelomeric regions. The centromeres mostly consist of an L1 LINE retrotransposon, ZeppL-1_cRei (Craig et al. 2021a), and we observed expansions and contractions of these sequences (generating 7.0% of TRMs), consistent with centromere evolution via satellite-like mechanisms such as unequal exchange rather than by active transposition. Subtelomeres feature a unique satellite called Sultan, and we observed expansions and contractions of the Sultan monomer within the same subtelomere, as hypothesized by Chaux-Jukic et al. (2021). A small number of TRMs involved the expansion or contraction of tandemly repeated gene families. Examples included 5S ribosomal RNA arrays and clusters of the large Chlamydomonas-specific NCL gene family on Chromosome 15, which encode RNA-binding proteins and appear to be experiencing rapid and ongoing evolution in C. reinhardtii (Boulouis et al. 2015).

Finally, we observed one example of a “mobile” satellite, which caused four insertions ranging from ∼0.5 to >21 kb in length in three different CC-2931 MA lines. This satellite, MSAT-11_cRei, consists of a ∼1.9-kb monomer and does not feature any characteristics typical of a TE. We have recently observed mobile insertions of MSAT-11_cRei in other C. reinhardtii strains (CC-1690 and CC-4532) (Craig et al. 2022). Although the mechanisms mediating satellite dissemination are not well understood, the phenomenon has been observed in other species and may be an important mechanism of satellite evolution (Ruiz-Ruano et al. 2016).

Whole-chromosome duplication and terminal deletions

We explored the possibility of whole-chromosome duplications using a genome coverage analysis and found evidence for one whole-chromosome duplication in CC-2931 L15 (Supplemental Fig. S11), which gives an estimated rate of 1.2 × 10⁻⁴ per genome per generation in the CC-2931 genotype. The existence of the duplication is supported by nearly 2× coverage of both the PacBio and Illumina reads relative to their median whole-genome coverage. Furthermore, we identified “heterozygous” mutations of all types along this chromosome, which presumably accumulated after the chromosome duplication event. One example of a heterozygous SM is a reciprocal translocation between one copy of the duplicated Chromosome 1 and the single copy of Chromosome 11 (Supplemental Fig. S11). Aneuploidy events, such as chromosome gain, can likely be attributed to chromosome segregation errors and are generally expected to be deleterious (Krasovec et al. 2022), although they can also facilitate adaptation in some instances (Gilchrist and Stelkens 2019). We did not observe chromosome loss, which would be lethal. However, we did observe the deletion of a chromosome terminus in two MA lines (CC-1952 L15 and CC-2931 L1), followed by de novo telomere addition (i.e., “telomere healing”) (Supplemental Fig. S12). These events resulted in terminal deletions of 94 kb and 163 kb, in CC-1952 L15 and CC-2931 L1, respectively.

Calling SMs from long reads and assemblies

De novo mutations are inherently rare events, and highly accurate methods are therefore required for their detection. Our ability to detect SMs was aided by major advances in long-read sequencing technology, long-read and whole-genome aligners, and structural variant callers. Nonetheless, our results show that the detection of structural variants remains much more challenging than the detection of shorter variants, making it necessary to use a combination of approaches. Because there was only a partial overlap between variants identified by different callers, it is possible that some SMs were undetected in our study. Although a read-based caller (Sniffles) was required for the identification of many variants, especially duplications, this approach failed to detect the full range of SMs. We found that the pangenome approach implemented by Cactus and vg was particularly successful, calling ∼80% of all curated SMs.

However, in all cases, our results show that structural variant callers are likely to yield high rates of false positives, even when the samples analyzed are nearly isogenic, as in our experiment. Although we were able to detect many genuine SMs in tandem repeats, we found that these regions were responsible for the majority of false-positive calls. Accurately calling SMs in tandem repeats may therefore require specific alignment and variant calling tools to be developed, and given that manual curation of variants is unlikely to be manageable in larger and more complex genomes than that of C. reinhardtii, masking of tandem repeats may be appropriate in automated analyses. The substantial contribution of assembly errors to rejected calls may also warrant the use of multiple assemblers in variant calling, or the development of methods that combine assembly-based structural variant detection with read-based verification. Overall, we recommend sequencing samples with sufficient coverage to enable de novo assembly, followed by both assembly and read-based SM or variant calling. If possible, manual verification with visualization tools such as IGV should also be performed.

Mechanisms underlying SM and between-strain heterogeneity of rates and spectra

Excluding expansions and contractions of tandem repeats, ∼79% of SMs were associated with TEs. Of these, 84% were TE insertions and excisions, which therefore formed major components of the SM spectra in both strains. The remaining 16% of TE-associated SMs included most inversions and translocations in CC-2931 and almost all of the duplications in CC-1952. The C. reinhardtii genome contains more than 200 diverse TE families, almost all of which show evidence of recent activity (Craig 2021). However, none of the three TEs implicated in causing rearrangements or duplications have previously been observed as active elements. Two of these active elements, Cryptons (Goodwin et al. 2003) and Dualen LINEs (Kojima and Fujiwara 2005), were first described from multicopy repeats in genetic and genomic data, whereas the active Replitrons found here have not been described previously.

Beyond TEs, we found little role for homology-based mechanisms of DSB repair, with the possible exception of the microhomology-mediated end joining pathway, which putatively caused deletions of moderate lengths. This is in contrast to SM studies in yeast, in which nonallelic homologous recombination has been shown to be the predominant mechanism mediating deletions, duplications, and rearrangements (Sui et al. 2020). Similarly, 10%–20% of de novo SMs in humans are thought to be mediated by nonallelic homologous recombination (Parks et al. 2015). Instead of homology-based mechanisms, the nonhomologous end joining repair pathway may have been involved in many of the SMs, particularly those not mediated by TEs. This is consistent with the very low rates of homologous recombination observed in C. reinhardtii under vegetative growth, where the species is haploid and nonhomologous end joining is the dominant DSB repair pathway (Ferenczi et al. 2021). The relative rates at which species repair DSBs via either nonhomologous end joining or homologous recombination have been implicated in many aspects of genome evolution, such as the evolution of base composition (Weissman et al. 2019) and intron density (Farlow et al. 2011), and it is likely that variation in the activity of DSB repair pathways also leads to substantial variation in SM spectra among species.

We found substantial variation in the rates of SNMs and indels between the CC-1952 and CC-2931 strains, both in our previous analysis with Illumina sequencing (Ness et al. 2015) and in the reanalysis reported herein using PacBio HiFi. There is substantial nucleotide diversity among C. reinhardtii strains (Flowers et al. 2015; Craig et al. 2019), and some of the within-species variation in mutation rates that we have observed may be caused by the presence of mutator alleles in certain strains. Consistent with our findings for SNM and short indel mutations, the SM rate in CC-1952 MA lines was significantly lower than in CC-2931. Furthermore, the SM spectra differed substantially between the strains. In particular, CC-2931 has a higher overall rate of transposition than CC-1952, involving a more diverse array of TEs. TE suppression is not well understood in C. reinhardtii but likely occurs at transcriptional and post-transcriptional levels (van Dijk et al. 2006) via mechanisms including repressive histone modifications (Jeong et al. 2002; Zhang et al. 2002) and RNA interference (Casas-Mollano et al. 2008). Because the CC-2931 and CC-1952 genomes both harbor many more potentially active TE families than we observed to be active in our study, we infer that all but a few families are silenced effectively. TE suppression may differ between CC-2931 and CC-1952 owing to genetic variation in genes involved in silencing pathways. Alternatively, environmental factors could impact TE activity. Almost nothing is known about local adaptation in C. reinhardtii, although natural isolates differ in their growth rates under laboratory conditions (Morgan et al. 2014; Kraemer et al. 2017). CC-1952 and CC-2931 were sampled ∼1600 km apart, from Minnesota and North Carolina, respectively, and they may differ in the extent of their adaptation to the highly artificial laboratory environment. Such differences could potentially cause stress-related interactions with transposition, which can result in TE activation or repression (Horváth et al. 2017).

In addition to a higher rate of transposition, CC-2931 MA lines also appear to have experienced a higher rate of DSBs than CC-1952 MA lines, which may explain the higher rates of duplications, deletions, and rearrangements in MA lines of this strain, even after TE-mediated SMs have been accounted for. DSBs are generally considered to be the most mutagenic DNA lesions (So et al. 2017) and are induced by intrinsic cellular factors (e.g., replicative and oxidative stresses) and by exogenous sources (e.g., mutagens). The rate of DSBs in C. reinhardtii is, however, not well understood, and differences between the strains could arise from genetic or environmental factors. Overall, we infer that the CC-2931 genome appears to be less stable than that of CC-1952, although the reasons for this are currently obscure.

Evolutionary implications of high and variable rates of SM

Population-level long-read sequencing projects have generally found that structural variants are common but segregate at low frequencies, implying that most are strongly deleterious (Chakraborty et al. 2019; Weissensteiner et al. 2020). Although we have not explored the relationship between SMs and fitness, the genomic distribution of SMs, in many cases overlapping coding sequences, suggests that many have large fitness effects. In contrast, a growing body of evidence suggests that coding sequences are less mutable than intergenic sequences with respect to SNMs and indels (Lee et al. 2012; Krasovec et al. 2017; Belfield et al. 2018; López-Cortegano et al. 2021; Monroe et al. 2022). We also found a substantial number of tandem repeat expansions and contractions in genic regions, a class of mutations that has been implicated as a major source of human genetic disease (Hannan 2018). Furthermore, the genomic distribution of TE insertions was not random. Considering the most active elements in CC-2931, Replitron insertions were enriched in intergenic sequences distant from genes, whereas CryptonF-1_cRei insertions were overrepresented in 5′ UTRs and gene proximal intergenic regions. Similar insertion biases in other TE families are well documented, and for example, P elements in Drosophila and Mu elements in maize show biases toward transcription start sites of highly expressed genes (Zhang et al. 2020). Gene proximal TE insertions can have important effects on gene expression, via the disruption of regulatory sequences, regional effects of transcriptional silencing, or the deposition of new regulatory elements (Cridland et al. 2015; Uzunović et al. 2019; Rech et al. 2022).

Consistent with the association between TEs and chromosomal rearrangements, the breakpoints of inversions and translocations also shared a similar genomic distribution to that of TE insertions. We expect that many of these events may have large fitness effects. Although not a factor in clonally propagated MA lines, many of the observed rearrangements are likely to cause meiotic incompatibilities. Furthermore, it is very likely that we have underestimated the true rate of translocations; we only observed translocations that retained a single centromere per chromosome, likely owing to highly reduced fitness following the formation of acentric chromosomes. Overall, the high fraction of the total mutation rate explained by SMs, together with their genomic distribution, suggests that SMs contribute substantially to the mutation load but also to the potential for adaptive evolution.

Notably, we also observed a strong bias toward genome size expansion in MA lines, which was driven by TE insertions and duplications. In some cases, the ancestral copy number of a TE increased severalfold in an MA line genome over the course of the experiment. The C. reinhardtii genome is highly compact, and TE copy number is generally low within TE families (Craig 2021), suggesting that the observed tendency toward genome expansion in MA lines is counteracted by purifying selection in nature. Based on short-read data, similar expansion biases have been reported for duplications in MA lines of Caenorhabditis elegans (Konrad et al. 2018) and for TE insertions in MA lines of Daphnia magna (Ho et al. 2021). Briefly, it is important to note that copy-and-paste TEs may effectively act as mutators in an MA experiment because each new insertion has the potential to create a new active copy and subsequently increase transposition rate. Thus, our observed rates of TE insertion, which were averaged over almost 1050 generations, may be overestimated relative to those occurring in nature.

Our results highlight the prevalence and importance of TEs among SMs and support a prominent evolutionary role for TEs. The heterogeneity in the rate and identity of TE insertions between CC-2931 and CC-1952 contributed substantially to the overall differences in SM rates and spectra between the strains. Our results suggest that species, populations, and even individuals may differ considerably in their SM spectra as a result of their active TE repertoire. More work would be required to investigate the rate of new SMs in other species and to elucidate the generality of results observed here. In particular, it will be important to test whether similar within-species variation in SM rates and spectra exists in other taxa. Our results add to the weight of evidence supporting the importance and prevalence of SMs and should further encourage structural variant discovery via assembly-based methods.

Biological samples and nucleic acids extraction

The MA line ancestors were C. reinhardtii wild strains CC-1952 (from Minnesota, 1986) and CC-2931 (North Carolina, 1991), which were originally obtained from the Chlamydomonas Resource Center (https://www.chlamycollection.org/). The MA experiment was conducted by Morgan et al. (2014). Briefly, MA lines were initiated from the ancestor strains and cultured on Bold's medium agar plates under white light at 25°C. MA lines were bottlenecked at regular intervals of 3–5 d by randomly picking single colonies and transferring them from one plate to another. MA lines were maintained for estimated averages of 1066 and 1050 generations for CC-1952 and CC-2931, respectively, after which Illumina sequencing was performed. The original ancestors and MA lines from the last transfer of the experiment were cryopreserved in liquid nitrogen.

For this study, we reconditioned the CC-1952 and CC-2931 ancestors along with several MA lines after cryopreservation and grew all samples in liquid Bold's medium before transferring to agar slants to produce stock cultures. Four CC-1952 MA lines (L1, L3, L6, and L15) and eight CC-2931 MA lines (L1, L2, L6, L9, L11, L13, L14, and L15) were sequenced together with the two ancestors. MA lines were randomly selected after excluding lines with combined SNM and indel mutation rates greater than or less than 1.5 times the interquartile range for that strain (Ness et al. 2015), because these lines may have accumulated mutations that modify the ancestral mutation rate. Cells were inoculated in six-well plate liquid cultures and grown for 4 d under constant light to produce sufficient biomass for DNA extraction. High-molecular-weight genomic DNA was extracted using a cetyltrimethylammonium bromide (CTAB) and phenol:chloroform protocol, following the method of Craig et al. (2021a). RNA was extracted in triplicate from independent cultures of the CC-2931 ancestor grown in liquid Bold's medium under constant light via a Maxwell RSC 48 instrument.

Nucleic acid sequencing

All of the CC-1952 samples and six of eight CC-2931 MA lines were sequenced on the PacBio sequel II platform with a 30-h movie, using the circular consensus sequencing (CCS) mode to generate HiFi reads. Samples were multiplexed, with between four and six samples per SMRT cell. Library and barcoding preparation, sequencing, CCS analysis, and demultiplexing were performed at the Earlham Institute (Norwich, United Kingdom). The mean read length N50 was 20.2 kb per sample, and mean coverage was about 27× per sample.

The CC-2931 ancestor and CC-2931 MA lines L2 and L9 were sequenced on individual SMRT cells using the PacBio sequel I platform with a 10-h movie, using the CLR mode. Library preparation and sequencing was performed at Edinburgh Genomics (Edinburgh, United Kingdom). Mean read length N50 was 19.3 kb, and mean coverage was about 54× per sample.

CC-2931 ancestor RNA-seq library preparation was conducted with the NEB mRNA stranded library preparation kit and Illumina NovaSeq preparation. The three replicate samples were sequenced using Illumina 100-bp paired-end sequencing. Library preparation and sequencing were performed by Génome Québec (Montreal, Canada).

Genome assembly

A different assembly strategy was followed for MA lines and ancestors. MA lines were assembled de novo to the contig-level. Flye v2.8.2 (Kolmogorov et al. 2019) was selected for assembly because it produced assemblies that were most representative of the haploid state (i.e., other assemblers yielded redundant contigs in repetitive regions). A genome length of 111.1 Mb was assumed (“-g 111.1m” in Flye; command line arguments are shown here and elsewhere in double quotation marks). Postprocessing was performed with purge_dups (Guan et al. 2020). Error correction was only performed for the two CC-2931 MA lines (L2 and L9) sequenced using CLR by applying two iterative rounds of the Arrow algorithm (Hepler et al. 2016). To facilitate their use as reference genomes, the ancestor genomes were assembled de novo to the chromosome-level using a combination of assembly methods and extensive manual assembly (see Supplemental Material).

Assembly metrics were quantified using QUAST v5.0.2 (Gurevich et al. 2013). Assembly completeness was estimated using BUSCO v4.0.6 (Manni et al. 2021), which was run in genome mode using the chlorophyta_odb10 data set (“augustus_species chlamy2011”). A list of all samples with a summary of their assembly approaches and quality metrics can be found in Supplemental Table S1.

Genome annotation

Structural annotation of the CC-2931 ancestor assembly was performed using BRAKER2 v2.15 (Brůna et al. 2021). RNA-seq data were mapped to the assembly using STAR v2.7.9 (“‐‐alignIntronMax 5000 ‐‐twopassMode Basic”) (Dobin et al. 2013). BRAKER2 was first run on an assembly softmasked for all repeats (see below) using existing AUGUSTUS parameters for C. reinhardtii (“‐‐species = chlamy2011 ‐‐skipAllTraining”). A second BRAKER2 run was performed with UTR prediction (“‐‐stranded = +, ‐‐UTR = on”), where the input BAM alignments were split into forward- and reverse-strand read sets using SAMtools v1.9 (Danecek et al. 2021). Because the run without UTR prediction returned a superior BUSCO score (“protein” mode), we designated this the primary annotation but also attempted to add UTRs where possible. For gene models in which coding sequence coordinates had a one-to-one correspondence between the two BRAKER2 runs (with variation permitted at only one exon for models with more than two exons), the model with UTRs was introduced as a replacement. Gene models were filtered if their coding sequence intersected in >30% with TE sequence or >70% with simple repeats (i.e., microsatellites) identified with RepeatMasker v4.0.9 (https://www.repeatmasker.org).

Repeat annotation was performed for the CC-2931 and CC-1952 ancestor assemblies. A custom library of C. reinhardtii TEs and satellites (Craig 2021) was first passed to RepeatMasker. A small number of TE families newly identified in this study were manually curated and added to the library (see below). The RepeatMasker annotations were supplemented with additional microsatellites and satellites identified by Tandem Repeats Finder (“2 7 7 80 10 50 1000 -f -d -m -ngs”) (Benson 1999). A final set of tandem repeats for each assembly was produced by combining simple repeats and satellites from RepeatMasker, satellites and microsatellites from Tandem Repeats Finder, and manually curated centromeric, subtelomeric, and ribosomal DNA array coordinates. Centromeres were identified based on the span of the constituent LINE ZeppL-1_cRei (Craig et al. 2021a), and subtelomeres were identified as sequence between the telomere and the characteristic spacer sequence, which coincides with the transition from the subtelomeric repeats (Chaux-Jukic et al. 2021). Tandem repeat sequences were considered to be microsatellites if they had a monomer length of <10 bp, and as satellite DNA if their monomer was >10 bp.

Mapping and alignment

To produce the input files for Sniffles, MA line PacBio read mapping was performed against the appropriate ancestor reference assembly using pbmm2 align (https://github.com/PacificBiosciences/pbmm2), adding the “‐‐preset CCS” flag for HiFi samples. For input to DeepVariant (see below), additional processing with SAMtools view (Danecek et al. 2021) run with the option “-F 256” was performed to remove reads with a flag indicating secondary alignment.

For visualization and curation of SMs, MA line assemblies were aligned to the appropriate reference assembly using minimap2 (“-x asm5”) (Li 2018). We also used unimap v0.1-r41 (https://github.com/lh3/unimap), a tool derived from minimap2 and optimized for assembly alignments, which provided superior alignments across many tandem repeats.

A pangenome alignment was produced for each ancestor and its MA lines using Cactus v1.3.0 (Armstrong et al. 2020). All assemblies were first softmasked for repeats using RepeatMasker and Tandem Repeats Finder, as described above. Minigraph v0.15-r426 (“-xggs”) (Li et al. 2020) was used to produce the graphical fragment assembly (GFA) file that was provided as input for cactus-graphmap together with softmasked assemblies. The resulting PAF file was then passed to cactus-align (“‐‐pangenome ‐‐pafInput ‐‐outVG”) in order to produce the final pangenome in variation graph (VG) format.

Callable sites

Callable sites were defined as genomic sites where SMs could be called with high confidence, based on visualization of the sequencing and assembly data. We previously defined callable sites based on short-read mapping parameters (Ness et al. 2015; López-Cortegano et al. 2021). Here, we used two criteria based on the de novo genome assemblies of the ancestors and MA lines. First, the ancestor assembly was aligned against itself with minimap2 (“-x asm5”), and genomic regions that were absent in the resulting PAF file (i.e., that were unmapped) were deemed uncallable. This first criterion was used because these regions are essentially unmappable even as isogenic sequences (at least with minimap2) and are hence inaccessible to variant calling. Second, a similar procedure was followed for each MA line by aligning their assemblies to the ancestor genome and extracting unmapped genomic coordinates. The unmapped coordinates extracted from all MA lines per ancestor were then intersected using BEDTools “intersect” (Quinlan and Hall 2010), and any regions that were present in at least two MA lines were defined as uncallable, because an unmapped region in a single MA line could be an SM such as a large deletion. This second criterion was adopted because these regions are prone to assembly breaks across multiple lines, even though they are assembled in the ancestor references. Taking the output from both criteria, uncallable regions separated by <30 kb were merged. Finally, coordinates corresponding to active cut-and-paste DNA transposons were manually reincluded as callable, because the excisions at these sites could otherwise be classified as uncallable if the same TE was excised in multiple lines. Because the uncallable regions are enriched in tandem repeats (see Supplemental Material; Supplemental Dataset S2; for the example dotplot, see Supplemental Fig. S4), we expect that we may have underestimated the mutation rate of TRMs. However, there was no overrepresentation of other SM types in callable tandem repeats, suggesting that the callable regions of the ancestor assemblies provide near-complete references for detecting most SMs genome-wide.

To compare callable sites between our previous Illumina sequencing and the PacBio sequencing described here, callable site coordinates from Ness et al. (2015) were converted to correspond to our new ancestor assemblies. A whole-genome alignment of the v5 reference assembly and the ancestor assemblies was generated using Cactus (Armstrong et al. 2020). Coordinates were then lifted over to the relevant ancestor assembly using the HAL tools command halLiftover (Hickey et al. 2013).

SM identification

SMs were called using three different variant callers, each of which relied on a different underlying alignment tool. Sniffles v1.0.12b (Sedlazeck et al. 2018) was used to call SMs based on the pbmm2 read alignments described above. BAM files were preprocessed using SAMtools-calmd to generate the MD tag, which provides information on mismatching positions (i.e., variable coordinates in the reads). Sniffles was first run on each MA line individually, and the resulting VCF files were merged using SURVIVOR v1.0.7 (Jeffares et al. 2017). Following the pipeline recommended for population calling (https://github.com/fritzsedlazeck/Sniffles/wiki/), Sniffles was then run again with the merged VCF as input and the option “‐‐Ivcf.” This population calling enables consistent presence or absence calls for SMs across all MA lines within a strain. SURVIVOR was used again to generate a multisample VCF.

MUM&Co v3 (O'Donnell and Fischer 2020) was used to call SMs from individual alignments of MA line assemblies to their ancestral reference, setting a genome size of 110 Mb (“-g 110000000”). MUM&Co calls variants based on alignments produced by MUMmer v4 (Marçais et al. 2018), which is performed as part of a single script. Variants were obtained as TSV and VCF files.

The variation graph tool (vg) (Garrison et al. 2018) was used to call variants directly from the pangenome alignments using the deconstruct command (“‐‐path-traversals”). The resulting VCF file for each strain was reduced to variants >50 bp.

All called variants in callable regions were manually curated via visualization of read and assembly alignments using the Integrative Genomics Viewer (IGV) (Robinson et al. 2011). SMs were rejected if they were not supported unambiguously by the read alignments. Read support for very large SMs was visualized via Ribbon v1.1 (Nattestad et al. 2021), which enables the visualization of reads mapping to discordant genomic regions. Supplemental Figures S12–S26 provide examples of SM visualization and curation. Most variants were entirely spanned by the reads, leading to simple visual confirmation in IGV, but variants >30 kb in length (approximately the upper limit of read lengths), including large inversions and translocations, required additional curation. In addition to read support from Ribbon, these rearrangements were traced in the MA line assemblies by manually assessing the discordant mapping of MA line contigs in the PAF alignment files (see Supplemental Fig. S23). Complex SMs, including large rearrangements and duplications, were further visualized using Ribbon v1.1 (Nattestad et al. 2021).

Duplications and deletions were curated as tandem repeat expansions or contractions if they involved the duplication or deletion of one or more monomers of a tandem repeat. Most fell within existing tandem repeat annotations, that is, satellites and microsatellites, whereas a small number required manual inspection of indel flanks by self-vs-self dotplots generated using the MAFFT v7 online server (Katoh et al. 2019). Deletions that perfectly intersected with TEs annotated by RepeatMasker in the ancestor genome were called as mobile excisions. Mobile insertions for described TE families were identified as cases in which the inserted sequence had a near-perfect BLASTN match (Camacho et al. 2009) to the Chlamydomonas repeat library (Craig 2021). These hits all had expected length distributions; LINE and PLE insertions frequently only contained the 3′ end owing to 5′ truncation, whereas insertions of other TEs corresponded to the entire length of the TE. In cases in which an inserted sequence had no match to an existing TE model, we queried the insert sequence against the ancestor genome, extracted and aligned hits, and manually curated new consensus sequences following established protocols for mobile element annotation (Goubert et al. 2022). All insertions unambiguously matched either the existing or newly produced consensus sequences and could be neatly defined to specific mobile element families. The one exception to this pattern was the duplications mediated by Dualen LINEs, where the sequence called as an insertion partly matched Dualen-4b_cRei and partly matched the sequence immediately flanking the insertion. These Dualen-mediated duplications were manually split to two called SMs: one mobile insertion and one duplication of the appropriate lengths.

When curating inversions and translocations, we noticed that many events featured additional insertions at the rearrangement breakpoints that were not specifically detected by the variant callers. As above, these insertions were compared to the annotated TEs and defined as mobile insertions of specific TE families. Five rearrangements could not be fully characterized because one of the breakpoints was clearly supported, but the other was in an uncallable region. These were arbitrarily classified as translocations.

SNMs and indels

DeepVariant 1.1.0 (Poplin et al. 2018) was used for calling SNMs and indels based on read alignments, using the option “‐‐model_type = PACBIO.” DeepVariant was run on individual MA line pbmm2 alignment files. The resulting VCF files were merged using GLNexus 1.3.1 (“ ‐‐config Deepvariant_unfiltered”) (Yun et al. 2020). The merged VCF file was further processed to retain only high-quality calls (QUAL ≥ 20) of sites called as homozygous genotype, because DeepVariant assumes diploid genomes and biallelic variants, with a minimum read depth of eight. In addition, only variant calls that were unique to a single MA line were retained as mutations. All SNMs and indels were confirmed visually in IGV.

Genomic distributions of SMs and TE insertions

The coordinates of SM breakpoints in CC-2931 were intersected with genomic annotations (coding sequence, introns, etc.) using BEDTools. Additionally, we also calculated the intersection of deletions, duplications, and inversions based on the entire span of these SMs. The observed distributions of SMs were compared against null expectations based on random sampling of the callable genome. For the analysis of individual TE families, the random expectation was adjusted to account for insertion target sequences (see Fig. 4B). The target motifs were identified using SeqKit v2.1.0 (Shen et al. 2016), and these motifs were then sampled from the callable genome. TRMs and mobile excisions were not analyzed.

Sequence homology

To explore the role of different DSB repair pathways in SM generation, we searched for patterns of micro- and macrohomology at SM breakpoints, following the method of Belyeu et al. (2021). We first searched for evidence of macrohomology between the breakpoints of deletions, duplications, inversions, and translocations that were not previously associated with TEs (≥95% sequence identity detected by the megablast BLASTN algorithm) (Camacho et al. 2009). Query and target sequences were the 100 bp upstream of and downstream from each breakpoint. In addition to macrohomology, we looked for patterns of microhomology. Given that microhomology mechanisms such as microhomology-mediated end joining can require as little as 2 bp of homology, we manually compared the 20 bp of sequence surrounding breakpoints.