Introduction

Defining the genetic basis and evolutionary dynamics of adaptation is a central goal in evolutionary biology. Mutations underlying adaptation or biological innovation can depend on multiple factors including genetic backgrounds, phenotypic states, and genome architecture (Blount et al., 2008, 2012). One important class of mutation mediating adaptive evolution are copy number variants (CNVs) which comprise duplications or deletions of genomic sequences that range in size from gene fragments to whole chromosomes. Quantifying the rates at which CNVs occur, the factors that influence their formation, and the fitness and functional effects of CNVs is essential for understanding their role in evolutionary processes.

CNVs play roles in rapid adaptation in multiple contexts and are an initiating event in biological innovation. In laboratory evolution experiments, a spontaneous tandem duplication captured a promoter for expression of a citrate transporter and resulted in Escherichia coli cells, typically unable to use citrate, to start metabolizing citrate as a carbon source (Blount et al., 2012). CNVs can be beneficial in cancer cells, promote tumorigenesis (Ben-David & Amon, 2020), enhance cancer cell adaptability (Rutledge et al., 2016), and accelerate resistance to anti-cancer therapies (Lukow et al., 2021). Over longer time scales, CNVs serve as substrate from which new genes evolve (Ohno, 1970; Taylor & Raes, 2004) as duplicated genes redundant in function can accumulate mutations and evolve to acquire new functions. For example, the globin gene family in mammals arose from rounds of gene duplication and subsequent diversification (Storz, 2016). CNVs also contribute to macro-evolutionary processes and thereby contribute to species differences, such as between humans and chimpanzees (Cheng et al., 2005) and speciation (Zuellig & Sweigart, 2018).

Mutations, including CNVs, occur in part because of errors made during DNA replication or DNA repair. Two general processes underlie CNV formation: (1) DNA recombination-based mechanisms and (2) DNA replication-based mechanisms (Harel et al., 2015; Hastings, Lupski, et al., 2009; Malhotra & Sebat, 2012; Pös et al., 2021; Zhang, Gu, et al., 2009). Recombination-mediated mechanisms of CNV formation include non-allelic homologous recombination (NAHR) and nonhomologous end joining. NAHR occurs via recombination between homologous sequences that are not allelic. As such, NAHR occurs more frequently with repetitive sequences due to improper alignment of DNA segments and can occur either between (interchromosomal) or within (intrachromosomal) a chromosome (Harel et al., 2015). One prevalent class of repetitive sequence are retrotransposons and both full length and partial sequences, such as long terminal repeats (LTR), are substrates for homologous recombination generating gene amplifications (Avecilla et al., 2023; Dunham et al., 2002; Gresham et al., 2008; Lauer et al., 2018; Spealman et al., 2022). DNA replication-based mechanisms include fork stalling template switching (FoSTes) and microhomology mediated break-induced repair (MMBIR) (Carvalho et al., 2013; Gu et al., 2008a; Hastings, Ira, et al., 2009; Lee et al., 2007). During FoSTes and MMBIR, the DNA replication fork stalls due to a single strand nick and a replication error occurs in which the lagging strand switches to an incorrect template strand mediated by microhomology. Reinitiation of DNA synthesis at the incorrect site can form CNVs. A particular type of DNA replication-based error is Origin-Dependent Inverted Repeat Amplification (ODIRA), in which short inverted repeats near an origin of DNA replication enable template switching of the leading strand to the lagging strand. Subsequent replication generates an intermediate linear DNA molecule that can recombine into the original genome to form a triplication with an inverted middle copy (Brewer et al., 2015; Martin et al., 2024).

In microbes, CNVs can mediate rapid adaptation to selective conditions imposed through nutrient limitation in a chemostat. Selected CNVs often include genes encoding nutrient transporters that facilitate import of the limiting nutrient (Dunham et al., 2002; Gresham et al., 2008; Horiuchi et al., 1963; Payen et al., 2016; Sonti & Roth, 1989), likely as a result of improved nutrient transport capacity due to increased protein production. Previous studies have found amplification of the general amino acid permease gene, GAP1, when Saccharomyces cerevisiae populations are continuously cultured in glutamine-limited chemostats (Gresham et al., 2010; Lauer et al., 2018). Amplification of GAP1 confers increased fitness in the selective environment (Avecilla et al., 2023). Sequence characterization of these CNVs revealed that a diversity of de novo CNV alleles arose and were selected for including tandem duplications of GAP1, complex large CNVs, aneuploidies, and translocations. However, little is known about the molecular basis of this genetic diversity.

Local genome elements are likely to be an important determinant of CNV formation rates and mechanisms. Genomic context can influence multiple properties including mutation rate, epigenetic regulation, chromatin state, transcription levels, DNA replication, and recombination rate (Arndt et al., 2005; Chuang & Li, 2004; Lang & Murray, 2011; Lercher & Hurst, 2002; Matassi et al., 1999; Nishant et al., 2009; Wolfe et al., 1989). Prior work has shown that CNVs occur more frequently in repetitive regions in the genome (Harel et al., 2015; Pentao et al., 1992; Stankiewicz et al., 2003; Turner et al., 2008). However, little is known about the role of local genomic architecture and organization on CNV formation rates, the types of CNVs that are generated, their associated fitness effects, and ultimately the paths taken during adaptive evolution.

Here, we aimed to directly investigate the effect of local genome architecture elements on de novo GAP1 CNV formation and selection dynamics during adaptive evolution of Saccharomyces cerevisiae. We hypothesized that sequence elements proximate to GAP1 potentiate CNV formation. The GAP1 locus, which is located on the short arm of chromosome XI, consists of two flanking Ty1 long terminal repeats (LTRs) that share 82% sequence identity and an origin of DNA replication or autonomously replicating sequence (ARS) (Figure 1A). Both LTRs and ARS may facilitate GAP1 CNV formation due to their proximity. First, the flanking LTRs can undergo inter-chromatid NAHR to form tandem duplications of GAP1 on a linear chromosome (Lauer et al., 2018; Spealman et al., 2022). Second, intra-chromatid NAHR between the flanking LTRs can form an extrachromosomal circle containing GAP1 and an ARS able to self-propagate and integrate into the genome (Gresham et al., 2010). Finally, GAP1 triplications can form through ODIRA using short inverted repeats and the proximate ARS (Brewer et al., 2015; Lauer et al., 2018; Martin et al., 2024). These elements are thought to facilitate a high rate of GAP1 amplification, estimated to be on the order of 10^-4 per haploid genome per generation (Avecilla et al., 2022). To test our hypothesis we used a CNV reporter, wherein a constitutively expressed fluorescent GFP gene is inserted adjacent to GAP1 (Lauer et al., 2018). We engineered strains that lacked either the ARS (ARSΔ), both flanking LTRs (LTRΔ), or all three elements (ALLΔ) (Figure 1A). We performed experimental evolution using wildtype (WT) and genomic architecture mutant populations in glutamine-limited chemostats for 137 generations and quantified GAP1 CNVs using flow cytometry (Figure 1). Surprisingly, we find that the proximate DNA elements are not required for GAP1 CNV formation as GAP1 CNVs were identified in all evolving populations. We used neural network simulation-based inference (nnSBI) to infer the CNV formation rate and selection coefficient (Avecilla et al., 2022). We find that genomic architecture mutants have significantly reduced CNV formation rates relative to WT and significantly lower selection coefficients even though GAP1 CNVs repeatedly form and sweep to high frequency in all strains. We performed genome sequence analysis to define the molecular mechanisms of CNV formation for 177 CNV lineages and found that 49% of GAP1 CNVs are mediated by ODIRA regardless of background strain. In the absence of the local ARS, a distal ARS facilitates CNV formation through ODIRA. We also find that homologous recombination mechanisms still mediate gene amplification in the absence of LTRs in part initiated by de novo insertion of retrotransposon elements at the locus. Our study reveals the remarkable plasticity of the genome and that DNA replication errors are a predominant source of adaptive CNVs.

Figure 1.

Results

Accurate estimation of CNV allele frequencies remains challenging using molecular methods such as DNA sequencing and qPCR. To address this challenge we developed a CNV reporter comprising a constitutively expressed fluorescent gene inserted upstream of GAP1 and observed recurrent amplification and selection of GAP1 in glutamine-limiting conditions (Lauer et al., 2018). Subsequently, we showed that a high rate of GAP1 CNV formation and strong fitness effects explain the highly reproducible evolutionary dynamics (Avecilla et al., 2022). Noncoding sequence elements proximate to GAP1, including flanking LTRs in tandem orientation and an ARS, contribute to GAP1 CNV formation (Gresham et al., 2010; Lauer et al., 2018). Many studies have shown that repetitive sequence regions and origins of replications are hotspots of CNVs (Arlt et al., 2012; Cardoso et al., 2016; Di Rienzi et al., 2009; Gresham et al., 2010; Lauer et al., 2018; Martin et al., 2024). Thus, we hypothesized that the local genomic architecture of GAP1 facilitates its high rate of CNV formation.

To test the role of proximate genomic features we engineered strains deleted for each element and thus differing from our wildtype strain (WT) containing a GAP1 CNV reporter by a single modification. Specifically, we constructed ARSΔ, a strain lacking the single ARS, LTRΔ, a strain lacking the flanking LTRs, and ALLΔ, a strain lacking all three elements (Figure 1A). All strains contain the CNV reporter at the identical location as the WT strain. We confirmed scarless deletions of genetic elements using Sanger and whole-genome sequencing.

Local genomic architecture contributes to GAP1 CNV evolutionary dynamics

We founded independent populations with each of the three engineered strains lacking proximate genomic features and a WT strain. We studied GAP1 CNV dynamics in populations maintained in glutamine-limited chemostats over 137 generations (Figure 1). For each of the four strains, we propagated 5-8 clonal replicate populations, each originating from the same inoculum (founder population) derived from a single colony. Approximately every 10 generations, we measured GFP fluorescence of sampled populations using a flow cytometer and quantified the proportion of cells containing GAP1 CNVs (Methods). We observed similar CNV dynamics across independent populations within each strain (Figure S3B). Therefore, we summarized CNV dynamics for each strain using the median proportion of the population with a GAP1 CNV. In every strain, GAP1 CNVs are generated and selected resulting in qualitatively similar dynamics in WT and mutant strains (Figure 1B).

Deletion of the ARS, but not the flanking LTRs alters CNV dynamics

We quantified three phases of CNV dynamics 1) the time to CNV appearance, defined by the inflection point before rise in frequency (Fig 1C), 2) the selection phase, corresponding to the increase in proportion of CNVs per generation during the initial expansion of CNVs (Figure 1D), and 3) the equilibrium phase, corresponding to the plateau (Figure 1E). The time to CNV appearance (Figure 1C) and the CNV selection phase (Figure 1B) does not differ between WT and LTRΔ populations (pairwise wilcoxon test, adjusted p = 1, pairwise t-test p = 1, respectively). In the WT and LTRΔ populations, GAP1 CNVs appear at generation 50 and increase in frequency at similar rates, ∼15% per generation in WT and ∼18% per generation in LTRΔ. The two strains both reach their equilibrium phase at the same point around generation 75 (pairwise t-test, adjusted p =1). The absence of a significant difference in CNV adaptation dynamics between the two strains suggests that the LTRs are not a major determinant of GAP1 CNV evolutionary dynamics.

By contrast, in ARSΔ and ALLΔ populations, we observe a delay in the time to CNV appearance. In both of these strains, CNVs are first detected at generations 65-80, whereas in WT and LTRΔ populations CNVs are first detected at generation 50 (ARSΔ vs. LTRΔ, wilcoxon pairwise test, adjusted p = 0.0059, ALLΔ vs. LTRΔ, wilcoxon pairwise test, adjusted, p = 0.0124) (Figure 1). Thus, the local ARS contributes to the initial GAP1 CNV dynamics. Similarly, CNV selection is significantly different between the LTRΔ (18%) and ALLΔ (13%) (pairwise t-test, adjusted p-value = 0.0026). Finally, we also observe a significant delay (ANOVA, p = 0.00833) in the generation at which the CNV frequency reaches equilibrium in ARSΔ (∼generation 112) compared to WT (pairwise t-test, adjusted p = 0.05) (Figure 1E). These observations suggest that absence of the ARS in the ARSΔ and ALLΔ strains delays the appearance of GAP1 CNVs compared with the presence of the ARS in WT or LTRΔ strains.

GAP1 amplifications can occur without CNV reporter amplification

In both WT and LTRΔ populations we observed that GAP1 CNV abundance stabilized around 75% during the equilibrium phase (Figure 1B). In these populations, inspection of raw flow cytometry data revealed a persistent single-copy GFP subpopulation (Figure 2A). These data could be explained by two possible scenarios: (1) the existence of a non-GAP1 CNV subpopulation comprising beneficial variation at other loci with fitness effects equivalent to GAP1 CNVs or (2) lineages with GAP1 CNVs without co-amplification of the CNV reporter. To resolve these two possibilities we sequenced clones from the single-copy GFP subpopulation in each of the five WT populations (Supplementary table 1) and identified the presence of GAP1 amplifications without co-amplification of the CNV reporter (Supplementary figure S2A). The presence of GAP1 CNVs in four out of five populations lacking the CNV reporter amplification suggests it may have occurred in the founder population. We found eight distinct CNV breakpoint pairs (Supplementary figure S2A) indicating at least eight independent events occurred, either in the founder population or just after being propagated to independent populations. All clones from one of the five populations (i.e. population 3) contained one copy of GFP and one copy of GAP1, suggesting these clones have a beneficial mutation elsewhere in the genome allowing them to stably coexist with the GAP1 CNV subpopulation.

Figure 2.

Incorporating unreported pre-existing CNVs in an evolutionary model

To quantify the evolutionary parameters underlying empirically measured CNV dynamics (Figure 1) we built an evolutionary model with the goal of performing network simulation-based inference (nnSBI) (Avecilla et al., 2022; Cranmer et al., 2020; Gonçalves et al., 2020). Previously, our evolutionary model assumed the GAP1 CNV reporter allowed us to detect all GAP1 CNVs (Avecilla et al., 2022). However, our new experimental results indicate the existence of a small subpopulation of unreported GAP1 CNVs present at the beginning of the experiments (Figure 2A). Therefore, we expanded the evolutionary model to include φ, the proportion of cells with GAP1 CNVs without co-amplification of the reporter, at the commencement of the experiment (i.e. generation 0) (Figure 2B). The remaining model parameters are δ_C, the rate at which GAP1 duplications form, and δ_B, the rate other beneficial mutations occur. We find that this updated evolutionary model can accurately describe the observed dynamics (Supplementary Figure S3B), which are clearly affected by the value of φ. When the total CNV proportion is very different from the reported proportion, e.g., when φ≫δ_C > δ_B, a reduced CNV formation rate results in a greater discrepancy between reported and total CNV proportions (Figure 2C)

Decreased CNV formation rates in mutants suggest adjacent elements can drive GAP1 CNV formation

We used nnSBI to estimate CNV formation rates and selection coefficients from evolutionary dynamics observed in glutamine-limited chemostats. We trained a neural density estimator using evolutionary simulations (Methods). This neural density estimator then allows us to infer posterior distributions over model parameters (formation rate, δ; selection coefficient, s) from the CNV dynamics in single replicates and a collective posterior distribution from a set of replicates of the same strain. δ_C is the GAP1 CNV formation rate. s_C is the GAP1 CNV selection coefficient, wherein the fitness effect is 1 + s. We estimated the confidence of our inference approach on synthetic simulations by computing its coverage, i.e., the probability that the true parameter falls within the 95% highest density interval (HDI) of the posterior distribution. Our neural density estimator is slightly over-confident for φ with a coverage of 0.934, and under-confident for GAP1 CNV selection coefficient and formation rate with a coverage of 0.992 for s_C and 0.995 for δ_C, respectively. (Supplementary Table 2). Despite this under-confidence, the posterior distributions are narrow in biological terms: the 95% HDI represents less than an order of magnitude for both s_C and δ_C. Thus, we did not apply post-training adjustments to the neural density estimator, such as calibration (Cook et al., 2006) or ensembles (Caspi et al., 2023; Hermans et al., 2022).

We find that the individual maximum a posteriori (MAP) estimates vary across strains and replicates (Supplementary Figure S3A). Overall, the CNV selection coefficient, s_C, ranges from 0.1 to 0.22 (with one exception of 0.3) whereas the CNV formation rate, δ_C, ranges from 10⁻⁶ to 10⁻⁴ (with one exception of 10⁻³ and two of 10⁻⁷); and the proportion of pre-existing GAP1 CNVs that do not amplify the reporter φ ranges from 10⁻⁶ to 10⁻² (with two exceptions of 10⁻⁸). We found that the MAP estimates for replicate populations within the same strain cluster together with some outliers (Supplementary Figure S3A). We performed posterior predictive checks, drawing parameter values from the posterior distributions and simulating the frequency dynamics (Supplementary Figure S3H), which agree with the observed data (Supplementary Figure S3B). For each strain, we estimate the collective posterior distribution based on all individual posteriors, which generates a posterior distribution conditioned on all observations, P(θ|X₁, …, X_n) (Methods). The collective posterior allows us to estimate whether there is a difference in CNV formation rate and fitness effects across the four strains.

Collective posterior HDIs are very narrow (Figure 3A) and fit empirical observations (Figure 3C). The collective MAP estimates of the CNV selection coefficient are similar for the WT and LTRΔ (0. 182). For ARSΔ and ALLΔ, the selection coefficient is estimated to be lower, with values of 0.146 and 0.126, respectively. However, these selection coefficients are still large, consistent with these populations containing GAP1 CNVs that are highly beneficial under glutamine-limitation. The collective MAP estimate for the CNV formation rate in WT is 4. 5⋅10⁻⁵. In contrast with selection coefficients, the CNV formation rate is markedly lower in all strains ranging from 1⋅10⁻⁵ for LTRΔ and ALLΔ to 2. 4⋅10⁻⁶ in ARSΔ. These results are consistent with our hypothesis that proximate sequence features facilitate GAP1 CNV formation.

Figure 3.

In all strains, collective parameter estimations are highly correlated, as expected for joint estimation of selection coefficients and beneficial mutation rates (Gitschlag et al., 2023). Thus, the collective MAP predictions can be interpreted as a range of parameter values for each strain. Indeed, other than the very final time point for the WT population, all collective MAP predictions lay within the interquartile ranges (Supplementary Figure S3B). The observed GAP1 CNV frequency stabilizes at different levels in the different experiments (Figure 1B). This can be explained by pre-existing unreported CNVs with frequencies estimated to be between φ=4⋅10⁻⁶ to 1. 6⋅10⁻⁴ by the collective MAPs (Supplementary Figure S3C). Indeed, our model predicts that the total (reported and unreported) final CNV frequency is nearly one in all cases (Supplementary Figure S3G).

Next, we estimated the de novo CNV diversity. Previous work showed a diversity of CNV allele types formed under glutamine-limited selection including tandem duplications, segmental amplification, translocations, and whole chromosome amplification (Lauer et al., 2018), and that lineage richness decreases rapidly over the course of evolution due in part to competition and clonal interference (Lauer et al., 2018; Levy et al., 2015; Nguyen Ba et al., 2019). Our model does not include competition, clonal interference, or recurrent CNV formation. Therefore, diversity calculations are likely overestimations. Nonetheless, a comparison of diversity between strains is informative of whether proximate genome elements affect CNV allele diversity. Therefore, for each strain, we used its collective MAP to simulate a posterior prediction for the genotype frequencies (Figure 3C), which we then used to predict the posterior Shannon diversity (Jost, 2006). We predict the set of CNV alleles to be highly diverse: the final predicted Shannon diversity ranges from 1. 6⋅10⁴ in ARSΔ to 3. 2⋅10⁵ in WT (Figure 3B). Our model predicts that the diversity increases rapidly during the selection phase and stabilizes in the equilibrium phase. This is because CNV alleles that form towards the end of the experiment would have a low frequency with a minor effect on diversity. We observe the greatest diversity in WT populations with lower diversity in the three genomic architecture mutants. Moreover, diversity saturates faster in WT populations. This suggests that the WT strain is able to form more unique CNVs allele types earlier compared to the other three strains (Figure 3B). Shannon diversity is lower in LTRΔ and further lower in ALLΔ and ARSΔ (Figure 3B). This rank order is the same as the CNV formation rates (Figure 3A).

We used a modified version of the evolutionary model with the inferred parameters to simulate an evolutionary competition between WT and the three architecture mutant strains over 116 generations, a point in which CNVs have reached high frequency in the original experiment. To “win” these competitions, the competitor strains need to adapt to glutamine-limitation by producing CNVs (Supplementary Figure S3J and supplementary information). The results of the simulated evolutionary competitions predict that the WT would dominate over all three strains, as its predicted final proportion almost always exceeds its initial frequency of 0.5. The average predicted proportion of WT cells when competing with LTRΔ is 0.717. In contrast, ARSΔ and ALLΔ are predicted to be almost eliminated by generation 116, as the average predicted WT proportion is 0.998 and 0.999, respectively. These simulated competitions further suggest that the ARS is a more important contributor to adaptive evolution mediated by GAP1 CNVs.

Inference of CNV mechanisms in genome architecture mutants

Contrary to our expectations, removal of proximate genomic elements from the GAP1 locus does not inhibit the formation of GAP1 CNVs. We sought to determine the molecular basis by which GAP1 CNVs form in the absence of these local elements. Therefore, we isolated ∼40 GAP1 CNV-containing clones from each population containing the four different strains at generations 79 and 125 and performed Illumina whole-genome sequencing. Using a combination of read depth, split read, and discordant read analysis, we defined the extent of the amplified region, the precise CNV breakpoints, and GAP1 copy number. On the basis of these features, we inferred the CNV-forming mechanisms for each GAP1 CNV (Methods). Over the 177 analyzed GAP1 CNVs, we observed tandem amplifications, tandem triplications with an inverted middle copy, intra-chromosomal translocations, aneuploidy, and complex CNVs. GAP1 copy numbers range from two to six in any given clone. Each of the four strains is able to produce a diversity of CNV alleles ranging from small (tens of kilobases) to large (∼hundreds of kilobases) segmental amplifications (Figure 4). We identified six major CNV-forming mechanisms across the four strains: ODIRA, LTR NAHR, NAHR, transposon-mediated, complex CNVs, and whole chromosome duplication (aneuploidy) (Figure 4A and Methods).

Figure 4.

ODIRA is a predominant mechanism of CNV formation

We inferred GAP1 CNVs formed through ODIRA in all four genotypes at high frequencies: 22 out of 37 WT clones (59%), 42 out of 52 LTRΔ clones (81%), 11 out of 42 ARSΔ clones (26%), and 12 out of 46 ALLΔ clones (26%). Considering the set of all CNVs in all strains, ODIRA is the most common CNV mechanism comprising almost half of all CNVs (87/177, 49%). The second most common mechanism is NAHR between flanking LTRs (38/177, 21%), which generate tandem amplifications. In the WT background, ODIRA (22/37) and NAHR between LTRs (11/37) account for 89% of GAP1 CNVs.

In LTRΔ, GAP1 CNVs form via ODIRA, chromosome missegregation, and NAHR using other sites. As expected, in LTRΔ clones we did not detect NAHR between LTRs in 52 clones and as a result focal amplifications were not detected (Figure 4C). In LTRΔ, CNVs are formed predominantly by ODIRA (42/52, 81%) (Supplementary Table 3), a significant increase relative to WT clones (chi-sq, p = 0.02469) (Figure 4C). By contrast, aneuploidy (5/52), complex CNV (3/52), and NAHR (2/52) account for less than 10% of GAP1 CNVs in LTRΔ. Consequently, we observe an increase in GAP1 CNV size in LTRΔ relative to WT (Figure 4B) as there is an increased prevalence of segmental amplifications and aneuploidy (Figure 4E).

Aneuploidy was rarely observed. Whole amplification of chromosome XI was detected in six out of 177 clones (3.4%) (Figure 4C). It was detected only in 2 strains: WT and LTRΔ (Figure 4D).

ODIRA generates CNVs using distal ARS

Whereas removal of proximate LTRs abrogates the formation of CNVs through NAHR, removal of the local ARS does not prevent the formation of GAP1 CNVs through ODIRA (Figure 4C). In the absence of the proximate ARS, distal ones are used to form ODIRA as all amplified regions of ODIRA clones contain a distal ARS (Figure 4D), with one exception (see Methods, Supplementary Figure S4B). We observe an increase of LTR NAHR in the ARSΔ clones (27/52, 52%) relative to WT clones (11/37, 39%) (Figure 4C, chi-sq, p = 0.03083). In ARSΔ, we find two CNV size groups (Figure 4B). Small amplifications are formed via NAHR of the flanking LTRs (Figure 4E) and large segmental amplifications are formed via ODIRA and by NAHR between one local and one distal LTR (Figure 4E).

Novel retrotransposition events potentiate GAP1 CNVs

CNVs in the ALLΔ clones form by two major mechanisms: 1) ODIRA using distal ARS sites to form large amplifications and 2) LTR NAHR following novel Ty LTR retrotransposon insertions to form focal amplifications (transposon-mediated, Figure 4). ALLΔ clones have larger amplifications formed by ODIRA than ODIRA-generated amplification in WT and LTRΔ (Figure 4E) because they encompass distal ARS and inverted repeats (Figure 4D). Surprisingly, we detected novel LTR retrotransposon events that generated new LTRs that subsequently formed GAP1 CNVs through NAHR with a pre-existing LTR in the genome or an LTR from a second novel retrotransposition (Figure 4E). Regions upstream of tRNA genes are known to be hotspots for Ty retrotransposons (Ji et al., 1993; Mularoni et al., 2012). We find the novel retrotransposons insert near one or both of the previously deleted LTR sites (Supplementary File 1), which flank GAP1 and are downstream of tRNA genes (Figure 1A). We only detected novel retrotranspositions in ALLΔ populations. In total we detected 15 unique Ty retrotransposon insertion sites of which eight were upstream of the deleted LTR, YKRCδ11, and four were downstream of deleted LTR, YKRCδ12 (Supplementary File 1). The remaining two insertions were distal to the GAP1 gene: one on the short arm and the second on the long arm of chromosome XI. Every insertion was upstream of an tRNA gene, consistent with the biased preference of Ty LTR insertions (Ji et al., 1993; Mularoni et al., 2012). Recombination between a new and preexisting LTR produces large amplifications whereas recombination between two newly inserted Ty1 flanking the GAP1 gene forms focal amplifications of the GAP1 gene (Figure 4E).

Discussion

In this study we sought to understand the molecular basis of repeated de novo amplifications and selection of the general amino acid permease gene, GAP1, in S. cerevisiae continuously cultured in glutamine-limited selection. We hypothesized that a high formation rate of GAP1 CNVs is due to the unique genomic architecture at the locus, which comprises two flanking long terminal repeats and a DNA replication origin. We used genetic engineering, experimental evolution, and neural network simulation-based inference to quantify de novo CNV dynamics and estimate the CNV formation rate and selection coefficient in engineered mutants lacking the proximate genome elements. We find that removal of these elements has a significant impact on de novo CNV dynamics, CNV formation rate, and selection coefficients. However, CNVs are formed and selected in the absence of these elements highlighting the plasticity of the genome and diversity of mechanisms that generate CNVs during adaptive evolution.

Despite their proximity to GAP1 we found that flanking LTRs are not an essential driver of CNV formation. The de novo CNV dynamics of WT and LTRΔ populations are similar and we find that although the CNV formation rate is reduced, the effect is small. In contrast, a significantly decreased CNV formation rate and delayed CNV appearance time was observed in the absence of the ARS in ARSΔ and ALLΔ populations, which suggests that the local ARS is a major determinant of GAP1 CNV-mediated adaptive dynamics. Indeed, ODIRA was identified as the predominant CNV mechanism in sequence-characterized clones revealing that DNA replication errors are a major source of CNV formation during adaptive evolution.

The prevalence of ODIRA is a result of many replication origins and the pervasiveness of inverted repeat sequences throughout the chromosome (Figure 4). In particular, breakpoint analysis of LTRΔ CNV clones show that ODIRA produces a continuum of CNV sizes along the short arm of chromosome XI. Downstream breakpoints range from near the GAP1 gene (∼3 kilobases) to all the way to the right telomere of chromosome XI (153 kilobases) (Figure 4E). The S. cerevisiae genome contains a high number, 1-1000 per kilobase, of inverted repeats ranging from 3bp to 14bp (Martin et al., 2024). The longer repeats are more likely to be used in ODIRA (Martin et al., 2024). The ubiquity of inverted repeats is in stark contrast to the relative paucity of LTR sequences, which are dispersed throughout the genome. Thus, ODIRA supplies a diverse and high number of gene amplifications for selection to act on, setting the stage for genome evolution and adaptation.

Consistent with previous reports of increased Ty insertions in S. cerevisiae under stress conditions (Morillon et al., 2000, 2002), we observed novel retrotransposon insertions in populations evolved in glutamine-limited chemostats. Transposon insertions could be harmful and lead to loss-of-function mutations but are also substrate for generating beneficial alleles including CNVs (Blanc & Adams, 2003; Dunham et al., 2002; Gresham et al., 2008; Wilke & Adams, 1992). We only detected novel Ty insertions in the ALLΔ strain. This is likely because regions upstream of tRNA genes unoccupied by LTRs are predisposed to transposition. Our detection of novel retrotransposon insertions is consistent with a previous experimental evolution study that suggested that Ty insertions were rare under constant nitrogen-limitation and substantially more common under fluctuating nitrogen limitation, in which cells experience total nitrogen starvation periodically (Hays et al., 2023). We observe a substantially lower frequency of novel Ty retrotransposition events, 15/177 clones, compared to 898/345 clones, with multiple insertions per genome, isolated from serial transfer ammonium-limited experimental evolution (Hays et al., 2023). Importantly, the role of Ty differs in the two studies, as in our case beneficial CNV formed after novel retrotransposition through recombination of newly introduced LTR sequences, whereas Hays et al. found Ty associated null alleles that are beneficial in nitrogen-limited conditions. Our results reveal the prevalence of Ty insertions even under constant nitrogen-limiting conditions and demonstrate how new retrotransposition events can contribute to CNV formation.

Aneuploidy was not a major source of adaptation in our experiments as it was infrequently detected (n = 6/177). This contrasts with studies suggesting aneuploidy is a rapid and transient route to adaptation over short evolutionary time scales (Chen, Bradford, et al., 2012; Chen et al., 2015; Chen, Rubinstein, et al., 2012; Pavelka et al., 2010; Selmecki et al., 2015; Yona et al., 2012).

However, aneuploidy incurs a fitness cost (Robinson et al., 2023; Tsai et al., 2019; Yang et al., 2021) and therefore can be outcompeted by slow-forming but less costly beneficial mutations in large populations (Kohanovski et al., 2024). Our observations of higher frequencies of focal and segmental amplifications may be because they are less costly than whole-chromosome amplifications.

A variety of DNA replication errors generate CNVs. Replication slippage at palindromic DNA and DNA repeats can cause fork stalling and downstream CNV formation (Lee et al., 2007; Zhang, Khajavi, et al., 2009). DNA repeats can form secondary structures like R loops, cruciforms, non-B DNA structures, and hairpins which stimulate CNV formation (Gu et al., 2008b). Untimely replication, faulty fork progression, S-phase checkpoint dysfunction, defective nucleosome assembly, and DNA repeat sites including LTRs are sources of replication-associated genome instability (Aguilera & García-Muse, 2013). However, DNA replication alone might not be the only contributor. The GAP1 gene is highly transcribed under glutamine-limitation (Airoldi et al., 2016). Transcription-replication collisions may fuel ODIRA CNV formation at this locus (Lauer & Gresham, 2019; Wilson et al., 2015). CNV formation can also be stimulated by transcription-associated replication stress and histone acetylation (Hull et al., 2017; Salim et al., 2021; Whale et al., 2022) and replication fork stalling at tRNA genes (Osmundson et al., 2017; Yeung & Smith, 2020). Testing the role of transcription in promoting the formation of adaptive CNVs warrants further investigation.

Recent work has proposed that ODIRA CNVs are a major mechanism of CNVs in human genomes (Brewer et al., 2015; Martin et al., 2024). Studies of human and yeast genomes have typically considered homologous recombination as the predominant mechanism of CNV formation (Lupski & Stankiewicz, 2005). CNV hotspots identified in the human (Chance et al., 1994; Lupski, 1998; Lupski & Stankiewicz, 2005; Pentao et al., 1992) and yeast genomes are indeed mediated by NAHR of long repeat sequences (Green et al., 2010; Gresham et al., 2010). However, a focus on recombination-based mechanisms as a means of generating copy number variation may be the result of ascertainment bias or the comparative ease of studying the effect of long repeat sequences over short palindromic ones. Our study demonstrates that experimental evolution in yeast is a useful approach to elucidating the molecular mechanisms by which DNA replication errors generate CNVs.

Methods

Strains and Media

Each of the three architecture mutants were constructed independently starting with the GAP1 CNV reporter strain (DGY1657). To construct each deletion strain, we performed two rounds of transformations both using PCR amplified donor templates designed for homology-directed repair. The first transformation used a repair template containing a nourseothricin resistance cassette to replace the pre-existing kanamycin resistance cassette and GAP1 gene. The repair template was designed to also deleted the elements of interests (ie. ARS, both flanking LTRs, or both LTRs and ARS). The second transformation replaced the nourseothricin cassette with a kanamycin resistance cassette and GAP1 gene thus yielding a genomic architecture Δ strain. We confirmed scarless deletions with sanger sequencing and whole-genome-sequencing. Final identifiers are DGY1657 for the WT architecture strain, DGY2076 for the LTRΔ strain, DGY2150 for the ARSΔ strain, and DGY2071 for the ALLΔ strain.

DGY1, DGY500, and DGY1315 are zero-, one-, two-, copy GFP controls, respectively, described in Lauer et al, 2018 and Spealman et al. 2023 (Lauer et al., 2018; Spealman et al., 2023). Briefly, GFP under the ACT1 promoter was inserted at neutral loci that do not undergo amplification in glutamine-limited continuous culture. 400μM glutamine-limited media is described in Lauer et al. 2018.

Long-Term Experimental Evolution

We performed experimental evolution of 30 S. cerevisiae populations in miniature chemostats (ministats) for ∼137 generations under nitrogen limitation with 400μM glutamine as in Lauer et al. (2018). Of the 30 populations, there were three controls: one control population with no fluorescent reporter (DGY1), one with one GFP fluorescent reporter (DGY500), one with two GFP fluorescent reporters (DGY1315). The remaining 27 populations have the GAP1 CNV GFP reporter. Of these, five populations are WT (DGY1657), seven are LTRΔ (DGY2076), seven are ARSΔ (DGY2150), and eight are ALLΔ (DGY2071). We inoculated each ministat containing 20ml of glutamine-limited media with 0.5 ml culture from its corresponding genotype founder population. The founder population was founded by a single colony grown overnight in glutamine-limited media at 30°C. Replicate populations of the same strain were inoculated from the same founder population derived from a single colony. Strains were randomized among the 30-plex ministat setup to account for the possibility of systematic position effects. After inoculation, populations were incubated in a growth chamber at 30°C for 24 hours with the media inflow pump off. After 24 hours, the populations had reached early stationary phase and we turned on the media inflow pump and waited 4 hours for the populations to reach steady-state equilibrium, at which the population size was ∼10⁸ cells. This was generation zero. Ministats were incubated in a growth chamber at 30°C with a dilution rate of 0.12 culture volumes/hr. Since the ministats had a 20ml culture volume, the population doubling time was 5.8 hours. Approximately every 10 generations, we froze 2-ml samples of each population in 15% glycerol stored at -80°C. Approximately every 30 generations, we pelleted cells from 1-ml samples of each population and froze them at -80°C for genomic DNA extraction.

Flow Cytometry analysis to study GAP1 CNV dynamics

To track GAP1 CNV dynamics, we sampled 1-ml from each population approximately every 10 generations. We sonicated cell populations for 1 minute to remove any cell clumping and immediately analyzed samples on the Cytek Aurora flow cytometer. We sampled 100,000 cells per population and recorded forward scatter, side scatter, and GFP fluorescent signals for every cell. We performed hierarchical gating to define cells, single cells, unstained (zero-copy-GFP control) cells, cells with one copy of GFP (GAP1), and two or more copies of GFP (GAP1) (Spealman et al., 2023). First we gated for cells (filtered out any debris, bacteria) by graphing forward scatter area (FSC-A) against side scatter area (SSC-A). Second, we gated for single cells by graphing forward scatter area against forward scatter height and drawing along the resulting diagonal. Finally, we drew non-overlapping gates to define three subpopulations: zero copy, one copy, and two or more copies of GFP by graphing B2 channel area (B2-A), which detects GFP (excitation = 516 λ, emission = 529 λ), against forward scatter area (FSC-A). We note that the one copy and two copy events overlap some, which is a limitation in this experiment (Spealman et al., 2023).

We found that two architecture mutants, DGY2150 and DGY2071, had strain-specific GFP fluorescence even though they only harbored one copy of GFP. DGY2150 and DGY2071 had slightly higher fluorescence than the one copy GFP control strain, DGY500, but less than that of the two copy GFP control strain, DGY1315. The third architecture mutant, DGY2076, had the same GFP fluorescence as the one-copy GFP control strain (DGY500). We ruled out that they were spontaneous diploids by looking at forward scatter signals. The forward scatter signal was not different from that of the one copy control (a haploid) and was not as high as a diploid. Therefore due to strain-specific fluorescence, we decided to perform strain-based gating, ie. one set of gates for the WT strain, a second set of gates for the LTRΔ strain, and so on. Since the controls are also a strain of their own, they were not used to set universal gates for one-copy or two-copy. Thus, for each strain, we chose the basis of our one-copy gate as the timepoint per strain in aggregate with the lowest median cell-sized normalized fluorescence. The two-or-more-copy (CNV) gate was drawn directly above and non-overlapping with the one-copy gate.

Quantification of dynamics

To obtain the proportion of CNVs for each population at each timepoint, we applied gates that correspond to zero-,one-, and two-or-more copy subpopulations. Using such proportion per population per timepoint, we summarize population CNV dynamics as follows (Lauer et al., 2018; Spealman et al., 2023). We calculate the generation of CNV appearance for each of the evolved populations. We defined CNV appearance as the generation where the proportion of CNV-containing cells first surpasses a threshold of 10% for three consecutive generations. Next, modified from Lang et al. (2011) (Lang et al., 2011) and Lauer et al. (2018) (Lauer et al., 2018), we calculate the percent increase in CNVs per generation for each evolved population. We compute the natural log of the proportion of the population with CNVs divided by the proportion of the population without CNVs for each timepoint. These proportions were obtained previously by gating. We plot these values across time and perform linear regression during the initial increase of CNVs. The slope of the linear regression is the percent increase in CNVs per generation. Finally, we calculate the time to CNV equilibrium, as defined by the generation at which a linear regression results in a slope < 0.005 after the selection phase.

Neural network simulation-based inference of evolutionary parameters

Evolutionary model

We developed a Wright-Fisher model that describes the evolutionary dynamics, similar to our previous study (Avecilla et al., 2022). In this study we have shown that a Wright-Fisher model is suitable for describing evolutionary dynamics in a chemostat. This is a discrete-time, non-overlapping generations model with a constant population size. Every generation has three stages: selection, in which the frequency of genotypes with beneficial alleles increases; mutation, in which genotypes can gain a single beneficial mutation or CNV; and drift, in which the population of the next generation is generated by sampling from a multinomial distribution. Our model follows the change in frequency of four genotypes (Figure 2B): A, the ancestor genotype; B, a cell with a non-CNV beneficial mutation; C⁺, a genotype with two copies of GAP1 and two copies of the CNV reporter; and C⁻, a genotype with two copies of GAP1 but only a single copy of the CNV reporter. CNV and non-CNV alleles are formed at a rate of δ_C and δ_B and have a selection coefficient of s_C and s_B, respectively. The frequency of genotype i is X_i. Unlike X_B and X_C⁺, which may increase due to both mutation and selection, we assume that C⁻ is not generated after generation 0 (as experimental results suggest that the reporter is working properly). Hence, the C⁻ genotype only increases in frequency due to selection, with s_C as its selection coefficient. We assume C⁻ has an initial frequency φ. Model equations and further details are in the supplementary information.

Simulation-based inference

We use a neural network simulation-based inference method, Neural Posterior Estimation or NPE (Papamakarios, 2019) to estimate the joint posterior distribution of three model parameters, s_C, δ_C and φ, while the other parameters, s_B and δ_B are fixed to a specific value(Table 1). Inferring all five model parameters resulted in similar prediction accuracy and s_C and δ_C estimates.

Table 1.

We applied NPE implemented in the Python package sbi (Tejero-Cantero et al., 2020) using a masked autoregressive flow (Papamakarios et al., 2018) as the neural density estimator: an artificial neural network that “learns” an amortized posterior of model parameters from a set of synthetic simulations. Posterior amortization allows us to estimate the posterior distribution P(θ|Χ) for a new observation Χ without the need to re-run the entire inference pipeline, i.e., generating new simulations and re-training the network (as is the case in sampling-based methods such as Markov chain Monte Carlo or MCMC).

We generated 100,000 synthetic observations simulated from our evolutionary model using parameters drawn from the prior distribution (Table 1). The neural density estimator was trained using early stopping with a convergence threshold of 100 epochs without decreases in minimal validation loss (the default in sbi is 20). Using 100 epochs as a threshold resulted in improved predictions. We validated that this improvement in prediction accuracy is not a result of over-fitting (Supplementary Figure S3F).

We validated the trained neural density estimator by measuring the coverage property: the probability of parameters to fall within the estimated posterior marginal 95% HDI. Then, we used the distribution of (Supplementary Figure S3D) and posterior predictive checks (Supplementary Figure S3B) as quantitative and qualitative measures of prediction accuracy, respectively.

Collective posterior distribution

NPE estimates a posterior distribution per observation, i.e., P(θ|X). Given n observations X₁, …X_n generated from the same model distribution P(θ), where each observation is a time-series of GAP1 CNV frequency, NPE infers n individual posterior distributions, each conditioned on a single observation, P(X_i). We estimate the collective posterior distribution based on n individual posteriors, that is, a posterior distribution conditioned on all observations,

This can be computed using the individual posteriors P(X_i) and the prior P(θ) (see supplementary information for derivation).

However, as P(X_i) could be infinitesimally small, a single observation could potentially reject a parameter value that is likely according to other observations. We want the collective posterior to be robust to such non-representative observations. Therefore, we define P_∈(X_i) = max(∈, P(X_i)) and use this quantity instead of P(X_i) in eq. 1. For a correct choice of ∈, the collective posterior mode should reflect a value with high posterior density for multiple observations, rather than a value that no individual posterior completely rejects. Thus, we implemented an amortized collective posterior: given a set of observations, X₁, …X_n, an amortized posterior for individual observations, P(θ|X), a prior, P(θ), and ∈, a collective posterior distribution P(θ|X₁, …, X_n), can be computed. This distribution can then be used to evaluate the posterior distribution of specific parameters given the observations, generate random samples from the posterior distribution (e.g., using rejection sampling), or find the collective MAP.

We set ∈ = e⁻¹⁵⁰ based on a visual grid-search. To find the normalizing factor (denominator in eq. 1), the integral is approximated by a dense Riemann sum (300³ points). Maximizing the distribution, i.e., finding the collective MAP, is implemented using scipy’s minimize method with the Nelder-Mead algorithm.

Genetic diversity

Using our evolutionary model with the inferred parameters, we can estimate the diversity of CNV alleles in the experiments. For each strain, we used samples from its collective posterior to simulate a posterior prediction for the CNV allele frequencies (Figure 3C), which we then used to compute the posterior Shannon diversity (Jost, 2006), as detailed in the supplementary information.

Whole genome sequencing of isolated clones

Clones were isolated from archived populations and verified to harbor a GAP1 CNV by measuring GFP fluorescence signal consistent with two or more copies. Populations of each strain from generation 79 were streaked out from the -80°C archive on YPD and incubated at 30°C for 2 days. Plates containing single colonies were viewed under a blue light to view GFP fluorescent colonies by eye. Relative to the fluorescence of the 2 copy control strain, we picked single colonies that fluoresced as bright or brighter, reasoning that these colonies would likely contain GAP1 CNVs. Single colonies were used to inoculate cultures in glutamine-limited media and incubated at 30°C for 18 hours. The cultures were analyzed on the Cytek Aurora to verify they indeed harbored two or more copies of GAP1 based on GFP fluorescence signal. For Illumina whole genome sequencing, genomic DNA was isolated using Hoffman-Winston method. Libraries were prepared using a Nextera kit and Illumina adapters. Libraries were sequenced on Illumina NextSeq 500 platform PE150 (2 x 150 300 Cycle v2.5) or Illumina NovaSeq 6000 SP PE150 (2x 150 300 Cycle v1.5). We also used custom Nextera Index Primers reported in S1 Table Baym et al. 2015 (Baym et al., 2015).

Breakpoint analysis and CNV mechanism inference in sequenced clones

Reference genomes

We created a custom reference genome for each of the genomic architecture mutants. The custom reference genome containing the GAP1 CNV reporter in Lauer et al. 2018 (NCBI assembly R64) was modified to delete the flanking LTRs, single ARS, or all three elements.

Copy number estimation by read depth

The estimation of GAP1 copy number from read depth used is described in Lauer et al. 2018, except we searched for ≥ 1000 base pairs of contiguous sequence. CNV boundaries were refined by visual inspection.

Structural variation calling and breakpoint analysis

Whole genome sequences of clones were run through CVish, a structural variant caller (Spealman, 2019). Structural variant calling was also done on each of the ancestor genomes: WT, ARSΔ, LTRΔ, and ALLΔ. Output .bam files containing split reads and discordant reads of evolved clones and its corresponding ancestor were visualized on Jbrowse 2 to confirm locations of de novo CNV breakpoints and orientation of sequences at the novel breakpoint junctions. Novel contigs relative to the reference genome were outputted in addition to the supporting split reads that generated the contig. Blastn was used to verify orientation of contigs, namely inverted sequences used to define ODIRA (see Definitions of Inferred CNV Mechanisms).

Definitions of Inferred CNV mechanisms

We used the following liberal classifications for each CNV category. We called a clone ODIRA if we found inverted sequences in at least one breakpoint (Figure 4A). We define LTR NAHR as having both breakpoints at LTR sites (Figure 4A). This is evidence of recombination between the homologous LTR sequences. This mechanism typically forms tandem amplifications. In some cases, we find the hybrid sequence between two LTRs, but this is hard to retrieve in short-read-sequencing.

We define NAHR as having breakpoints at homologous sequences, with at least one breakpoint not at an LTR sequence (Figure 4A). We define transposon-mediated as a clone having a breakpoint at a novel LTR retrotransposon insertion and the other breakpoint at a different LTR site (Figure 4A). This supports that the newly deposited LTR sequence recombined with other LTR sequences (either pre-existing or introduced by a second de novo retrotransposition) to form CNVs. Sometimes we are able to recover the hybrid sequence between LTR sequences. We define complex CNV as having more than two breakpoints on chromosome XI and a read depth profile that looks like more than one amplification event occurred (multi-step profile). For the complex CNV clones, we were not able to resolve the CNV mechanisms due to the limitations of short-read sequencing.

Supplement

Supplementary Table 1.

Supplementary figure S2A.

Figure S3A.

Figure S3B.

Figure S3C.

Figure S3D.

Supplementary Table 2.

Figure S3F.

Fig S3G.

Fig S3H.

Figure S3J.

Supplementary Figure S4A.

Supplementary Figure S4B.

Supplementary Table 3.

Supplementary File 1. Ty-associated clones and locations of novel Ty insertions.

Supplementary File 2. CNV Clone Sequencing Analysis

Data availability

Sequencing data is available at SRA PRJNA1098800.

Other associated data are available here: https://osf.io/js7z8/.

Source code repository simulation-based inference: https://github.com/yoavram-lab/chuong_et_al.

Scripts for flow cytometry-based evolutionary dynamics and analysis of CNV clones: https://github.com/GreshamLab/local_arch_variants.

Acknowledgements

We thank all the members of the Gresham lab for helpful discussions, NYU Gencore for sequencing samples, and NYU High Performance Cluster for computing and storage. We thank Joshua Caleb Macdonald, Saharon Rosset, Uri Obolski, and Adi Stern for discussions and advice. This work was supported by NSF GRFP DGE1839302 (JNC), NIGMS T32GM132037 (JNC), NIGMS R01GM134066 (DG), R01GM107466 (DG), and R35GM153419 (DG), NIAID R01AI140766 (DG), NSF 1818234 (DG), Israel Science Foundation (ISF, Y.R. 552/19), US–Israel Binational Science Foundation (Y.R. & D.G. 2021276), Minerva Center for Live Emulation of Evolution in the Lab (Y.R.) fellowship from the Edmond J. Safra Center for Bioinformatics at Tel-Aviv University (N.B.N), and fellowship from the AI and Data Science Center at Tel-Aviv University (N.B.N).