Patterns of selection reveal shared molecular targets over short and long evolutionary timescales

Standing and de novo genetic variants can both drive adaptation to environmental changes, but their relative contributions and interplay remain poorly understood. Here we investigated the dynamics of drug adaptation in yeast populations with different levels of standing variation by experimental evolution coupled with time-resolved sequencing and phenotyping. We found a doubling of standing variation alone boost the adaptation by 64.1% and 51.5% in hydroxyuea and rapamycin respectively. The causative standing and de novo variants were selected on shared targets of RNR4 in hydroxyurea and TOR1, TOR2 in rapamycin. The standing and de novo TOR variants map to different functional domains and act via distinct mechanisms. Interestingly, standing TOR variants from two domesticated strains exhibited opposite resistance effects, reflecting lineage-specific functional divergence. This study provides a dynamic view on how standing and de novo variants interactively drive adaptation and deepens our understanding of clonally evolving diseases.


24
Connecting allele frequency, phenotype and fitness change in a causally cohesive manner is 25 challenging in both natural and clinical populations, but feasible in experimental populations.

26
Experimental evolution can reveal the molecular determinants of adaptation across a wide range 27 of biological systems with unprecedented resolution (Long et al., 2015). It can be initiated from 28 populations with known levels of standing variation, evolved under fixed selection regimes and 29 preserved ad infinitum as frozen fossil records that can be revived and studied in detail. Clonal 30 evolutions of initially isogenic populations have confirmed key theoretical predictions, notably 31 how competing clones carrying different beneficial mutations interfere with each other (clonal 32 interference), and how neutral or slightly deleterious mutations can hitchhike to higher 3 frequencies on the same clone as beneficial mutations (Barrick et al., 2009;Gerrish & Lenski, 1 1998;Herron & Doebeli, 2013;Kvitek & Sherlock, 2013;Lang et al., 2013;Levy et al., 2015; 2 Payen et al., 2016;Venkataram et al., 2016). Causal relationships in heterogeneous populations, 3 usually derived from sexual crosses of diverged parents, are much more challenging to pinpoint, 4 because of the number of variants that segregate in these populations and the linkage between 5 them. Nevertheless, experimental evolution using heterogeneous budding yeast, fly and Virginia 6 chicken populations have shown that standing variation alone can drive adaptation (Burke et al., 7 2010;Burke, Liti, & Long, 2014;Parts et al., 2011;Sheng, Pettersson, Honaker, Siegel, & 8 Carlborg, 2015) with no need for de novo mutations to emerge and spread. We recently 9 performed experimental evolution using heterogeneous populations derived from diverged West 10 African (WA) and North American (NA) natural yeast strains (hereafter referred to as "two-11 parent population"), in two different drugs rapamycin (RM) and hydroxyurea (HU) (Vázquez-12 García et al., 2017). RM is an inhibitor of the eukaryotic serine/threonine kinase TOR and HU is 13 an inhibitor of DNA replication. Specifically, budding yeast contains two TOR genes -TOR1 and 14 TOR2. They form two different complexes termed TOR complex 1 (TORC1) and TORC2. The 15 former contains either TOR1 or TOR2 and is uniquely sensitive to RM while the latter specifically 16 contains TOR2 and is insensitive to RM (Loewith et al., 2002). HU impair DNA synthesis by 17 inhibiting ribonucleotide reductase and preventing the reduction of ribonucleotides to 18 deoxyribonucleotides (Koç, Wheeler, Mathews, & Merrill, 2004). In the two-parent population,

19
we found drug-specific adaptive contributions of standing and de novo variants. Selection on 20 standing variation explained more of growth rate increases in HU (51%) than selection on de 21 novo mutations (23%) but less in RM (22% vs 70%).

22
Overall, the relative contribution of standing and de novo variants to adaptation depends on 23 multiple factors, including the degree of standing variation, the typical fitness effects of standing 24 and de novo variation, the selective constraints imposed by the environment and the relevant time 25 scales (Long et al., 2015). Theory predicts early adaptation in heterogeneous populations to be 26 faster because beneficial standing variants are immediately available and less likely to be lost by 27 drift (Barrett & Schluter, 2008). Standing variants are predicted to disproportionately drive 28 adaptation when de novo beneficial mutations are rare, have small selection coefficients, or when 29 the duration of selection is short (Hermisson & Pennings, 2005). However, strict experimental 30 comparisons of adaptation on standing and de novo variants are scarce. In particular, it remains to 31 be explored: (1) how the degree of standing variation affects the adaptation rate and yield, (2) 32 whether standing and de novo variants are selected in a shared target. These questions have a 33 direct bearing on our understanding of the evolution of resistance to chemotherapy and 4 antimicrobials (Palmer & Kishony, 2013;Turner & Reis-Filho, 2012). To this end, we evolved 1 highly-heterogeneous yeast populations derived from intercrossing four diverged parents over 12 2 consecutive meiotic generations (Cubillos et al., 2013) (hereafter referred to as "four-parent 3 population", Figure 1A) to fixed concentrations of RM and HU. In comparison to the two-parent 4 population, the four-parent population has approximately twice the genetic diversity segregating

5
(1 SNP/120bp vs. 1 SNP/230bp), indicating higher level of standing variation (Cubillos et al., 6 2013). We tracked the adaptation of these four-parent populations to the two drugs at high 7 resolution, comparing the molecular and phenotypic changes to that of the isogenic populations of 8 the four parental strains as well as the published two-parent populations (Vázquez-García et al., 9 2017). We found that the four-parent populations adapted earlier and faster than the two-parent 10 populations. Resistant standing and de novo variants were selected on shared mutational targets 11 (RNR4, TOR1 and TOR2). However, the standing and de novo variants of the TOR paralog genes 12 occur in different domains and conferred RM resistance via distinct mechanisms.  Figure 1A, Tables S1-S2). Two four-parent populations (F12_1 and F12_2) were 21 independently derived from the four parents by 12 rounds of intercrossing (Cubillos et al., 2013) 22 and were therefore highly heterogeneous at the onset of selection. We evolved two replicates of 23 each isogenic parental population and eight replicates of four-parent populations in batch-to-batch 24 selection regimes, storing a subsample of each batch (T0 to T14 in HU and T0 to T15 in RM) to 25 create a dense fossil record.

26
To track the adaptation dynamics comprehensively, we revived the frozen subsamples of all the 27 isogenic, four-parent populations and the previously published two-parent populations (Vázquez-28 García et al., 2017) across all the time points (Tables S1-S2). We estimated their fitness related 29 properties by both precise measurement of their doubling time and spotting assay (Figures 1B, 30 S1-S3). Over the whole RM experiment, the adaptive gain between four-parent and two-parent 31 populations was similar (45.3% vs. 42.6% of doubling time reduction, Mann-Whitney U-test, p = 32 0.96). However, the early adaptive gain (T0 to T2) was larger in the four-parent populations 5 (19.2% vs. 11.2% of doubling time reduction, Mann-Whitney U-test, p = 0.038), highlighting the 1 advantage of higher level of standing variation in driving expeditious adaptation. There was no 2 substantial late stage (last three time points) adaptation in either four-parent or two-parent 3 populations (5.1% of doubling time increase and 1.1% reduction respectively), reflecting 4 exhaustion of adaptive potentials within the experimental timescale. In HU, the adaptation was 5 slow, gradual and persisted to the end in both the four-parent and two-parent populations but with 6 seemingly greater adaptive gains in the four-parent populations (20.4% vs. 12.3% of doubling 7 time reduction, Mann-Whitney U-test, p = 0.06). Therefore, a doubling of segregating diversity 8 in the four-parent populations translated into more rapid and more remarkable adaptive gains in 9 both RM and HU. No observable adaptation to control condition (no drug) was observed (Figure 10 S1).

11
To measure the adaptive gains in individuals independently of their background population, we 12 isolated > 2,600 random clones from ancestral and a subset of the endpoint populations (Table S2) 13 and measured their doubling time separately. Before selection (T0), the variability in doubling 14 time between individuals of the four-parent population was much greater than that in the two-

27
Mann-Whitney U-test, p < 2.2 × 10 -16 ). Individuals drawn from the NA, SA and WE endpoint 28 populations reached the same level of adaptation as those from the four-parent populations, 29 whereas those from the WA populations adapted more slowly, which is consistent with their 30 weaker initial growth. Remarkably, only the NA managed to adapt to HU (28.2% of doubling 31 time reduction, Mann-Whitney U-test, p < 2.2 × 10 -16 ). Even though NA individuals failed to 32 reach the same adaptation level of four-parent individuals ( Figure 1C; mean endpoint doubling 33 time 3.16 vs 2.62 hours, Mann-Whitney U-test, p < 2.2 × 10 -16 ). The SA and WE individuals 6 grew worse at the end of HU selection than their respective ancestral states (6.5% and 4.2% of 1 doubling time increase). The WA individuals went complete extinction after two cultivation 2 rounds (T2), suggesting the lowest evolvability. In summary, we found that higher level of 3 standing variation positively impacts the rate of adaptation, the absolute adaptive gains and the 4 endpoint performance, with exact effects depending on the selection regime.

6
De novo mutations in TOR1 and FPR1 drive rapamycin adaptation in isogenic populations 7 To lay a solid foundation for understanding the adaptation of the highly heterogeneous four-8 parent populations, we sequenced the initially isogenic populations at multiple time points 9 (Tables S1-S3). As expected, de novo mutations drove adaptive evolution in isogenic populations.

10
In RM, we detected recurrent mutations in TOR1 and FPR1 ( Figure 2A). TOR1 mutations (six 11 mutations in three sites) emerged in all the eight populations, indicating TOR1 as a background 12 independent source of RM resistance. In contrast, FPR1 mutations (frame shift and start codon 13 disruption in two sites) emerged only in the two NA populations. Surprisingly, all the NA clones 14 carrying FPR1 mutations became haploids during selection. This may be a consequence of NA 15 diploids being highly prone to sporulate even in relatively rich medium (Cubillos, Louis, & Liti, 16 2009) and a strong selection for haploids carrying loss-of-function FPR1 mutations given that 17 they are fully recessive (Vázquez-García et al., 2017). The frequency increase of TOR1 and FPR1 18 mutations agreed well with the doubling time reduction of the populations in which they emerged 19 ( Figure 2A). This supports that they are true drivers of adaptation, rather than hitchhikers or 20 drifters, and that adaptation is genetic, rather than initially epigenetic and later genetically 21 assimilated (Gjuvsland et al., 2016).

22
To quantify the individual contributions of TOR1 and FPR1 mutations to RM adaptation, we 23 isolated and estimated the doubling time of individual clones carrying these mutations ( Figure 2B, 24 Table S4). Except for FPR1 Met1Ile, the doubling time reduction conferred by each individual 25 mutation in the relevant state (heterozygous or homozygous) equaled (e.g. TOR1 S1972I in WE), 26 or approached (>90%, e.g. TOR1 W2038L and S1972I in NA) the total doubling time reduction 27 of the population in which it emerged ( Figure 2B). Nearly all the clones from the evolved 28 populations carried one of these mutations; thus, they were capable of explaining almost the 29 complete adaptive gains (Heitman, Movva, & Hall, 1991). All RM-adapted populations 30 performed equally well in presence and absence of RM; thus RM adaptation had plateaued 31 ( Figure S4). TOR1 mutations recurrently emerging in different genetic backgrounds (Ser1972Ile 32 in NA, SA and WE and Trp2038Leu in WA and WE) consistently gave complete tolerance to 7 RM ( Figure 2B). The larger adaptive gain conferred by FPR1 Ile11X frame shift than by FPR1 1 Met1Ile stop in the NA background (69.8% vs. 32.9% of doubling time reduction, Mann-2 Whitney U-test, p = 2.7 × 10 -6 , Figure 2B) agreed with the near fixation of FPR1 Ile11X in 3 NA_RM_1 and the low frequency of FPR1 Met1Ile in NA_RM_2 (Figure 2A). Given that both 4 should be complete loss-of-function mutations, this distinction is intriguing. In the WE 5 background, the TOR1 Ser1972Ile homozygous clones grew faster than those with the 6 heterozygous mutation ( Figure 2B, 68.1% vs. 59.8% of doubling time reduction, Mann-Whitney 7 U-test, p = 1.9 × 10 -4 ), giving them a competitive edge and suggesting that continued selection 8 should drive the homozygote state to fixation. Such homozygous mutations should have occurred 9 via loss of heterozygosity, as demonstrated in our previous study (Vázquez-García et al., 2017).

10
Whole-genome population sequencing uncovered no copy number changes, except in both SA 11 populations where the sequencing depth of chromosome IX (chrIX) increased under RM selection 12 ( Figure 2C, Table S3). Copy number qPCR confirmed that the RM-evolved diploid SA clones 13 carried three or four rather than the normal two copies of chrIX. Given that extra chrIX copies 14 conferred dramatically increased heat sensitivity ( Figure S5F), we estimated that ~12.5% and 15 ~8.3% of the evolved population (SA_RM_2_T15) carried three and four copies chrIX copies 16 respectively based on the frequencies of heat-sensitive clones. This is roughly in agreement with 17 the estimates based on the sequencing depth analysis ( Figure 2C). All the SA clones with extra 18 chrIX copies carried the TOR1 Ser1972Ile heterozygous driver mutation. Mutated TOR1 clones 19 carrying three copies of chrIX grew faster in RM than those with two or four copies (Mann-20 Whitney U-test, p = 6.90 × 10 -5 and p = 3.43 × 10 -3 respectively) ( Figure 2B). To better 21 understand the interplay between the TOR1 Ser1972Ile mutation and the chrIX amplification, we 22 constructed a cross grid of diploid strains with all possible combinations of TOR1 (wild type or 23 mutated) and chrIX copy number (2-4 copies) and estimated their doubling time ( Figure 2D).

24
Overwhelmingly, the TOR1 Ser1972Ile mutation is the major contributor to RM resistance 25 (53.2% and 56.1% of doubling time reduction for heterozygous and homozygous mutation 26 respectively), with extra chrIX copies being marginally beneficial in the TOR1 Ser1972Ile clones.

27
In sharp contrast to RM selection, isogenic populations propagated under HU almost uniformly 28 failed to generate and maintain detectable de novo variants. No de novo driver mutations were 29 detected in WE, WA or SA populations, which probably explains the evolution failure of these 30 populations throughout the 32-day experiment (Figures S1, S3-S4). The RNR4 mutations 31 (Arg34Ile and Lys114Met) were detected in the NA background. The clones carrying these 32 mutations in heterozygous state showed a mean population doubling time reduction of 31.8%. 8 This adaptive gain was directly comparable to those of NA endpoint populations, and thus 1 capable of fully explaining their adaptive gains. ( Figure S4B).

3
De novo mutations in TOR1, TOR2 and FPR1 drive rapamycin adaptation in heterogeneous 4 populations 5 Both de novo and standing variants could contribute to adaptation in the four-parent populations.

6
Given that they were derived from four parents, the frequency spectrum of each parental allele is 7 centered around a median frequency of 0.21 (WA), 0.26 (NA), 0.26 (WE) and 0.26 (SA) at T0 8 (Cubillos et al., 2013). In comparison, the initial frequency of de novo mutations is extremely low, 9 arising during the crossing or selection phases (Vázquez-García et al., 2017). We called de novo 10 driver mutations in genes that were recurrent mutation targets in the eight four-parent populations 11 and found that FPR1, TOR1 and TOR2 harbor such mutations ( Figure 3A). Identical FPR1 12 mutations occurred in all the replicated populations derived from one intercrossed population 13 F12_1 and the same TOR1 mutations were found in all the replicated population derived from the 14 other intercrossed population F12_2; thus these drivers emerged during the shared crossing phase 15 and then expanded independently during the selection phase. Validating this assumption, the 16 same haplotype blocks increased in frequency in the replicated populations derived from the same 17 intercrossed population, reflecting expansion of the same clones present at T0 ( Figure S6).

18
Similar to the isogenic lines, we found recurrent mutations at the TOR1 1972 and 2038 amino 19 acid sites, further confirming that these were the primary RM selection targets and that they arise 20 independently of the genetic context. In one population (F12_2_RM_4) we also found a TOR2 21 Ser1975Ile mutation rising to high frequency ( Figure 3A). TOR2 Ser1975 is located in the RM-22 binding domain and is paralogous to TOR1 Ser1972, implying that RM driver targets are 23 conserved between TOR1 and TOR2 (Helliwell et al., 1994). Isolation and genotyping of single 24 clones from the population containing TOR2 Ser1975Ile showed that heterozygous and 25 homozygous clones co-exist. This explains its frequency higher than 0.50 in the population. The  In stark contrast to RM, we did not detect any de novo driver mutations in the four-parent 1 populations under HU selection. Nevertheless, adaptation to HU is obvious ( Figures 1B-C, S4D), 2 and genome-wide frequency changes of parental alleles showed broad jumps at later time points 3 ( Figure S7), indicating resistant clones rising to high frequencies. We therefore conjectured that 4 standing variation largely drove the adaptation to HU in the four-parent populations.

6
Standing variation provides multiple selection targets to drive adaptation in heterogeneous 7 populations 8 We next investigated how the standing variation in the four-parent populations contributed to RM 9 and HU adaptation. We searched for genomic regions (quantitative trait loci, QTLs) with a 10 consistent change in the frequency of one or more alleles across both time points and replicated 11 populations. At later time points (T4 to T15), we observed strong shifts of allele frequencies over 12 large genomic regions, reflecting drug-resistant clones rising to high frequency in both selection 13 regimes ( Figures S7-S8). Therefore, we analyzed allele frequency changes before the clones arose 14 (T0-T4 for HU and T0-T2 for RM) to map QTLs using 99% and 95% quantiles cut-offs (see 15 Materials and Methods).

16
In HU, two QTLs passed the 99% quantile cut-off and seven more QTLs passed the 95% cut-off 17 ( Figure 4A and Table 1) with a median size of 22 kb and containing on average 10 genes. The  (Table S5, 29 Figure S9) while the allele frequency changes of the other three QTLs wore off before the end of 30 selection. Given that there were no detectable de novo driver mutations, the latter was probably 31 due to overwhelming competition from clones carrying the beneficial versions of the stronger 32 QTLs (e.g. RNR4 WE , ENA SA ).
10 Similarly, we mapped QTLs for RM resistance by analyzing allele frequency changes from T0 to 1 T2. We identified four QTLs at the 99% quantile cut-off (Table 1, Figure S10A). The two 2 strongest QTLs (52 and 26 kb respectively) covered the TOR1 and TOR2 genes respectively.

3
Interestingly, the WE and SA alleles of TOR1 and TOR2 showed opposite allele frequency 4 changes: TOR1 SA and TOR2 WE were selected for while TOR1 WE and TOR2 SA were selected against 5 ( Figure 5A). We validated such parental-specific allele preference by reciprocal hemizygosity 6 (Figures 5B, S5C-D). Clones carrying the strong TOR1 SA showed significantly shorter doubling 7 time and higher yield than the ones with the weak TOR1 WE allele (Mann-Whitney U-test, p = 3.1 8 × 10 -4 and p = 1.5 × 10 -4 respectively). Clones carrying strong TOR2 WE allele showed significantly 9 higher yield than clones carrying TOR2 SA (Mann-Whitney U-test, p = 1.5 × 10 -4 ). Nine additional 10 QTLs passed the 95% quantile cut-off. While we have not experimentally validated their effects, 11 we considered SNQ2, NPR3, KOG1 and CFT8 to be strong candidates for driving these QTLs 12 based on previous studies. Among them, CFT8 also contributed to RM resistance in the two-  Figure S9), despite the frequency increase of clones carrying de novo driver mutations.

23
Shared selection targets between standing and de novo variants in RNR4, TOR1 and TOR2

24
The multi-hit de novo mutations and QTLs identified in the same genes (RNR4 in HU and TOR1, 25 TOR2 in RM) showed a pattern of selection on shared molecular targets over short and long 26 evolutionary timescales. To understand why this pattern arose, we compared the standing variants 27 with de novo mutations identified in isogenic, two-parent and four-parent populations (Table 2).

28
The HU-resistant RNR4 WE allele had a single derived amino acid change, Ala161Thr, located 29 within the ribonucleotide reductase domain; this substitution was predicted to be functional 30 critical by sequence conservation analysis (see Materials and Methods). The RNR4 de novo driver 31 mutations emerging in both NA and two-parent populations were in the same domain but the 32 exact sites differed (Arg34Ile, Arg34Gly, and Lys114Met). All the de novo mutations in the 33 11 TOR1, TOR2 paralogs occurred in the highly conserved RM-binding domain, where they 1 prevented the binding of the FKBP12-RM complex and thereby conferred RM resistance. In 2 sharp contrast, none of the TOR1 and TOR2 standing variants mapped in the RM-binding domain 3 and occurred outside any characterized functional domains with the exception of TOR1 WE 4 Phe1640 ( Table 2). Three derived amino acid changes were unique to the weak TOR2 SA allele 5 (Glu122Gly, Ile1369Met, Ile1872Leu) and were all predicted to be deleterious (Table 2). To 6 expand our understanding of natural genetic variation of these shared selection targets, we 7 compared the sequences of RNR4, TOR1 and TOR2 across >1,000 S. cerevisiae natural isolates 8 (http://1002genomes.u-strasbg.fr/). All the three genes were well conserved ( Figure S11). A total 9 of 9, 79 and 73 amino acid sites of RNR4, TOR1 and TOR2 respectively were predicted to be 10 functionally critical, based on sequence conservation (Table S6). All nine RNR4 sites were in the 11 ribonucleotide reductase domain in which the standing and de novo variants driving HU 12 adaptation were located. About 38.0% (30/79) of the amino acid sites in TOR1 and 46.6% (34/73) 13 in TOR2 were located in known domains, including four TOR1 and two TOR2 sites in the RM-14 binding domain. We experimentally confirmed that natural alleles TOR1 His2000 and TOR2 15 Leu2047 in the RM-binding domain conferred RM resistance ( Figure S5E). This adds additional 16 support to that drug resistance can emerge through selection on existing natural variants that 17 prevent drug binding, with no need for de novo mutations to emerge.

18
Given that the TOR1, TOR2 standing and de novo variants often co-existed in the same 19 population, we further investigated their potential interactions. We genotyped the local genetic 20 background of TOR1 de novo mutant clones drawn from different endpoint populations and found 21 their background genotypes can be different, i.e. the TOR1 de novo mutations had arisen on 22 different clones (Table S4). Thus, TOR1 mutations conferred strong adaptive gains across 23 different clone backgrounds. This lack of a specific interplay between TOR1 de novo driver 24 mutations and their genetic background is consistent with the observation that TOR1 de novo 25 mutations also emerged and reached high frequency in all the isogenic populations. The interplay 26 between TOR2 mutations and their genetic background is particularly evident in a TOR2 clone 27 from population F12_2_RM_4. This clone carried the weak TOR2 SA allele, whose frequency 28 dropped from 0.29 (T0) to 0.04 (T8). However, the frequency was nevertheless enough for one 29 TOR2 SA clone to acquire a compensatory TOR2 Ser1975Ile mutation in a heterozygote state.

30
Consequently, the growth performance drastically increased and the TOR2 SA allele frequency was 31 driven to 0.46 at the end of selection ( Figure S12B). Thus, the emergence of this de novo TOR2 32 mutation in the TOR2 SA allele compensated for the sensitivity of the TOR2 SA allele by preventing 12 RM binding. This indicated that standing and de novo variants of TOR probably act towards RM 1 resistance via distinct mechanisms.

3
Functional consequences of TOR natural variants 4 TOR1 and TOR2 are master regulators of growth, controlling yeast performance in many 5 environments of relevance to industry, particularly in alcoholic beverage production. In this 6 industrial context, we were particularly intrigued by the opposite RM resistance phenotype of the 7 SA and WE TOR1, TOR2 alleles, as they occur in two lineages independently domesticated for 8 alcoholic beverage production (Fay & Benavides, 2005). We therefore further characterized the 9 TOR1 and TOR2 alleles in these two genetic backgrounds to determine their respective impact on 10 yeast performance in environments of industrial and medical interest. Potentially, such benefits 11 could also explain their distinct evolutionary trajectories.

21
Next, to understand the effects of TOR1, TOR2 standing variants on TORC1 activity, we used a 22 highly specific commercial antibody to measure phosphorylation of the ribosomal protein S6 23 (Rps6) under RM exposure. Rps6 phosphorylation is regulated by TORC1 and used as a specific 24 in vivo assay for TORC1 activity (González et al., 2015). Rps6 phosphorylation increased in 25 strains with the strong TOR1 SA and TOR2 WE alleles ( Figure 5D). Thus, the SNPs distinguishing 26 these alleles enhance TORC1 activity. This was quite surprising, first because a majority of SNPs 27 in these alleles occur outside functional domains and second because not even mutations in the 28 RM-binding domain affect TORC1 activity (González et al., 2015). This further underscored that 29 the standing and de novo variants of TOR1 and TOR2 cause RM resistance by distinct 30 mechanisms.

31
RM is an unlikely selection pressure on natural yeast alleles; however, real ecological constraints 32 such as nitrogen limitation do affect cell growth in a TOR-dependent functions (Loewith & Hall, 33 1 TOR2 allele preferences between SA and WE backgrounds and illuminate the underlying 2 mechanism, we measured their effect on doubling time in 18 relevant environments, including 3 nitrogen-limitations and synthetic wine must ( Figure 5E). As expected, the WE strain grew the 4 fastest in synthetic wine must, consistent with its niche-specific domestication history. Overall, 5 the removal of one TOR allele tended to result in growth defects in nitrogen-limited 6 environments, and the removal of a WE allele was generally worse than the removal of a SA 7 allele. For example, hybrids carrying TOR1 WE grow faster than those carrying TOR1 SA on 8 methionine and threonine; and hybrids with TOR2 WE grow faster than those with TOR2 SA in 9 tryptophan, threonine, serine, methionine, isoleucine, asparagine and adenine.

10
Finally, we investigated TOR gene essentiality in SA, WE, WA and NA genetic backgrounds by 11 knocking out TOR1 or TOR2. Previous studies in the laboratory strain S288C showed that TOR1 12 was non-essential, whereas TOR2 was essential (Liu et al., 2015;Winzeler et al., 1999). As 13 expected, TOR2 could not be deleted in WE, NA or WA. Surprisingly, however, we successfully 14 deleted the TOR2 gene in the SA haploid. The tor2Δ SA strain was able to grow on synthetic 15 complete medium (SC), although with marked growth defects, but not on YPD ( Figures 6A-B).

16
Because the TORC1 activity of tor2Δ SA remained unaltered upon RM treatment ( Figure 6C), 17 the SA background is either able to make up for the TOR2 loss by compensatory induction or by 18 complex incorporation of TOR1, or do not use TOR2 in TORC1 at all. We further dissected ~900 19 spores from WE/SA TOR2 reciprocal hemizygous deletions, as well as WE/SA wild type on both 20 YPD and SC medium. On SC, the spore viability was 83.5% for the homozygote TOR2 cross and 21 55.3% for the hemizygote cross. Thus, TOR2 was essential in a fraction of the recombined 22 WE/SA offspring. By tracking the deletion marker, we estimated that 18.5% of the tor2Δ spores 23 carrying recombinants survived on SC (Table S7) (Table S8). Taken together, the divergent 29 functions of the WE and SA alleles on TOR1 and TOR2 alleles impacted not only on RM 30 resistance but also on chronological aging, TORC1 activity, nitrogen control of growth and

32
Towards the mid and later phase of selection, highly resistant clones emerged and arose to high 33 frequency in both HU and RM. Nevertheless, the genetic make-up and origin of these clones 15 differed dramatically between the two selection regimes. Standing variants appeared to drive HU 1 adaptation all the way to the end, implying that beneficial de novo mutations are either too rare or 2 too weak to compete against the bulk dynamics driven by the standing beneficial alleles (i.e. the 3 nine QTLs). It could also be partially explained by negative or sign epistasis weakening the 4 effects of beneficial de novo alleles (Khan, Dinh, Schneider, Lenski, & Cooper, 2011). An 5 additional explanation is that the RNR driver mutations appear to be strongly background 6 dependent. This is evident from the isogenic populations, with only one background (NA) that 7 acquired RNR4 mutations and evolved. In contrast, mid to late adaptation to RM was consistently 8 driven by clones with de novo mutations in TOR1, TOR2 and FPR1 emerging and overtaking 9 other competing bulk subpopulations. This was consistently true in all genetic contexts and at all 10 levels of standing variation. These highly penetrant de novo mutations in members of the TOR 11 pathway have long been known to prevent their interaction with RM (Heitman et al., 1991; 12 Helliwell et al., 1994), which is manifested again by our experiments.

13
In the widest sense, we found strong examples of convergent selection on standing and de novo 14 variants to confer RM resistance -TOR1 and TOR2. This was not given a priori. First, strong for this drastic difference is that TOR2 is under stronger selection constraints likely reflecting its 30 unique, essential role in the TOR complex 2 (TORC2).

31
The standing WE and SA variants of TOR1 and TOR2 have opposite effects on RM resistance, 32 reflecting lineage-specific functional divergence after the gene duplication in their shared 33 ancestor. Although we cannot stringently reject a purely neutral explanation, the directly 16 opposing effects of these alleles on growth and survival corresponds to an evolutionary trade-off 1 between the two key determinants of fitness and makes it tempting to invoke selection to explain 2 this divergence. Domestication to distinct human made niches (Sake and grape-wine), including 3 different substrates of fermentation (Giudici & Zambonelli, 1992;Sasaki et al., 2014), may be the 4 ultimate explanation for this divergence with drug resistance as a side-effect caused by TOR

18
We previously performed two independent intercrosses to generate the F12 populations (four-19 parent populations -F12_1 and F12_2) that were derived from four diverged parents: 20 DBVPG6044 (West Africa, "WA"), DBVPG6765 (Wine European, "WE"), Y12 (Sake, "SA") 21 and YPS128 (North America, "NA"). The strain information is listed in Table S9. Here 22 experimental evolution was initiated from random subsamples of F12_1 and F12_2, with each 23 subsample comprised of 10 7 -10 8 cells. In parallel, experimental evolution was also initiated from 24 clonally expanded, near isogenic parental populations of similar size. Cells were evenly spread on 25 YPD agar plates (2% peptone, 1% yeast extract, 2% glucose, 2% agar) with hydroxyurea (10 26 mg/ml) or rapamycin (0.025 µg/ml), and incubated at 23°C. Every 2-3 days, all the cells were 27 collected from each plate into 1 ml distilled water. Ten percent of the cell suspension was 28 transferred to a freshly made plate while the rest were kept in 25% glycerol at -80°C. The 29 selection experiment lasted for 32 days. The detailed timeline and population specifics are listed 30 in Tables S1-S2. For each drug, there are four independently evolving replicates derived from 31 F12_1 and F12_2 respectively, as well as two replicates for each of the four parents. Besides, 32 17 there were two replicates derived from F12_1 and F12_2 respectively using drug-free YPD as 1 control. Procedures were identical to those used for generating and evolving the previously 2 published two-parent population (Vázquez-García et al., 2017). DNA was extracted from 3 populations of T0, T1, T2, T4, T8 and the last transfer using "Yeast MasterPure" kit (Epicentre, 4 USA). The samples were sequenced with Illumina TruSeq SBS v4 chemistry, using paired-end 5 sequencing on Illumina HiSeq 2000/2500 at the Wellcome Trust Sanger Institute. Sequence data 6 is deposited to NCBI SRA database with accession number for BioProject PRJEB4645.

11
Given our allele frequency estimates, we used a 10-kb sliding window with a 2-kb step size to 12 localize quantitative trait loci (QTLs). For each heterogeneous population, we compared the allele 13 frequency change in a window ! between time point ! and T0 (e.g.
If there is selection on standing variation, the absolute frequency change 15 of a parental allele in regions under selection is expected to be higher than in neutral regions and 16 to increase gradually as selection proceeds. On this basis, for each earlier transfer, we calculated 17 z-score of allele frequency changes compared with T0 in each population:

27
We searched for regions with z-score square higher than 99% or 95% quantile for each earlier 28 time point. If the regions were able to pass the cut-off at T1, T2 for RM and at T1, T2, T4 for HU,

29
and not pass the same cut-off in control (drug-free condition), they are assumed to be QTLs

30
(Figures 4A and S10, Table 1). We excluded regions located near chromosome ends, which could 31 be false positives due to repetitive sequences. The discrepancy of QTL numbers between the two-32 parent and the four-parent populations cannot be attributed to the different approaches to perform 33 1 only few weak QTLs were mapped ( Figure S13) including CTF8.
2 QTLs could be either maintained until later time points or be hijacked by the spread of clones 3 with beneficial mutations. We define whether a QTL is maintained by counting the replicates in 4 which the strong allele keeps increasing or the weak allele keeps decreasing until T4, T8 and the 5 end. If the number of such replicates is more than six (of a total of eight), we defined the QTL as 6 maintained until the later time points ( Figure S9, Tables 1 and S5). We randomly selected thousands of isolates from the initial and final populations (Table S2), bulk 11 population from the isogenic, two-parent and four-parent populations at each serial transfer of the 12 experimental evolution (Table S1) and strains with gene deletion (

24
We also used the Tecan Infinite 200 PRO plate reader to measure growth curves in small scale.

25
We pre-cultured the cells overnight and diluted the saturated culture 100 times into fresh medium.

26
We measured OD 600 every 15 minutes for at least 3 days in drugs and control. The raw OD 600 27 values were corrected and then used to generate growth curves. Doubling time and yield were 28 extracted using the online tool "PRECOG" (Fernandez-Ricaud, Kourtchenko, Zackrisson,

31
We did serial dilution and spotting of the cells to visualize the adaptation on population level 32 visually ( Figure S1) as well as the growth phenotypes of gene deletions ( Figure S5). Cells were 20 pre-cultured in YPD overnight to saturation. Then 5 µl of the culture was taken to do spotting 1 assay in the condition of interest. There were a total of six 1:10 dilutions from left to right on the 2 plate.

3
We also did spotting assay of 48 isolates drawn from the SA population evolved in RM 4 (SA_RM_2_T15) in heat condition (40ºC). We pre-cultured cells in YPD overnight. Then 5 µl 5 cells of 1,000-fold dilution from saturation were taken to put on YPD and incubated at 40ºC. The 6 plates were scanned after two days.  We isolated diploid SA clones from a RM-evolved population at T15. We validated the copy 3 number of chrIX (three or four copies), induced sporulation (in 2% KAc) and dissect spores. We 4 genotyped the spores of the mating type, TOR1 mutation or wild type and chrIX copy number 5 (one or two copies). With these genotypes, we crossed spores to create an array of diploids where 6 all possible genotypes were combined ( Figure 2D).     genomes. The first haploid pseudo-genome of each isolates was used for our downstream analysis.

12
We extracted the CDS regions of the RNR4, FPR1, TOR1 and TOR2 genes from these haploid

22
Genome evolution and adaptation in a long-term experiment with Escherichia coli.

27
Revealing the genetic structure of a trait by sequencing a population under selection.      showed Rps6 phosphorylation in SA wild type and tor2Δ strains shows that TORC1 activity is 11 not altered by the TOR2 deletion. All conditions are similar to the one reported in Figure 5D. (D)

12
Representative plates acquired 4 days after tetrad dissection on SC and YPD for WE/SA wild 13 type and its TOR2 reciprocal hemizygotes. The red circles indicate viable tor2Δ strains.
14 15 Figure S1. Spotting assay of all the four-parent and isogenic populations in this study.

16
Populations were sampled at T0, T1, T2, T4, T8 and T14 (HU) or T15 (RM and control). Each  Table S1. WA in HU died out after T2. The second replicate of WE in RM was 16 17 Figure S10. Identification of QTLs in RM and control conditions. We used allele frequency 18 from T0 to T2 to identify QTLs for RM resistance (A) and T0 to T4 for drug-free condition (B).

19
Dashed and solid lines indicate 99% and 95% quantile cut-off respectively. The red and black 20 labels respectively indicate the QTLs passing 99% and 95% cut-off.

30
Circle indicates the positions of standing variants and star indicates de novo mutations (Table 2).    We aligned the sequences of RNR4, TOR1 and TOR2 from the four parents, the S. cerevisiae 3 reference strain S288C and S. paradoxus reference strain CBS432 and extracted all the amino 4 acid changes. The unique amino acid change among the four parents are shown in red. We also 5 show de novo mutations identified in isogenic, two-parent and four-parent populations. We 38 predicted their functional impact by the online tool (mutfunc.com). The predicted scores are listed 1 in #Score and #IC columns.

2
The column description is as follows: 3 # Domain: predicted by Pfam 4 # Score: SIFT score, any mutation with a score below 0.05 occur in a highly conserved site and 5 are predicted to be deleterious.

6
# IC: information content at each position of the alignment. Here, a high value indicates strong 7 conservation, where the maximum value is 4.32.
8    1 * The isolates were from the ancestral or the final populations. The population of "WE_RM_2" 2 was contaminated at T15 and was replaced by T8.
3 4 ** Driver mutations identified by whole-genome population sequencing in the final population.