Dynamics and variability in the pleiotropic effects of adaptation in laboratory budding yeast populations

Evolutionary adaptation to a constant environment is driven by the accumulation of mutations which can have a range of unrealized pleiotropic effects in other environments. These pleiotropic consequences of adaptation can influence the emergence of specialists or generalists, and are critical for evolution in temporally or spatially fluctuating environments. While many experiments have examined the pleiotropic effects of adaptation at a snapshot in time, very few have observed the dynamics by which these effects emerge and evolve. Here, we propagated hundreds of diploid and haploid laboratory budding yeast populations in each of three environments, and then assayed their fitness in multiple environments over 1000 generations of evolution. We find that replicate populations evolved in the same condition share common patterns of pleiotropic effects across other environments, which emerge within the first several hundred generations of evolution. However, we also find dynamic and environment-specific variability within these trends: variability in pleiotropic effects tends to increase over time, with the extent of variability depending on the evolution environment. These results suggest shifting and overlapping contributions of chance and contingency to the pleiotropic effects of adaptation, which could influence evolutionary trajectories in complex environments that fluctuate across space and time.

hundreds of diploid and haploid laboratory budding yeast populations in each of three 23 environments, and then assayed their fitness in multiple environments over 1000 generations of 24 evolution. We find that replicate populations evolved in the same condition share common 25 patterns of pleiotropic effects across other environments, which emerge within the first several 26 hundred generations of evolution. However, we also find dynamic and environment-specific  Our results allow us to investigate differential roles for chance and contingency over

100
To study the dynamics of the pleiotropic consequences of adaptation, we experimentally evolved After completing the evolution, we revived populations from generation 0, 200, 400, 600, 800, 119 and 1000. We then conducted parallel bulk fitness assays (2 technical replicates) to measure the 120 fitness of each population at each timepoint across five environments (the three evolution 121 environments plus YPD + 0.4M NaCl at 30°C (transfers every 24 hours) and YPD at 21°C 122 (transfers every 48 hours), environments which exposed the populations to unique osmotic and 123 temperature stresses). In each bulk fitness assay, we pooled all populations from a given 124 generation along with a small number of common reference clones and propagated them for 50 125 generations ( Fig 1B). We then sequenced the barcode locus at generation 10, 30, and 50, and we 126 inferred the fitness of each population from the change in log frequency of each corresponding 127 barcode. By exploiting the fact that each population is uniquely barcoded, these bulk fitness  Based on the measured fitness of the generation 0 ancestral populations, we found that some 132 diploid populations had substantially higher ancestral fitness in certain assay environments, 133 likely because they acquired mutations prior to the start of the evolution. To clarify our 134 downstream analyses, we excluded 19 outlier diploid populations whose ancestors differed from 135 the mean ancestral fitness by at least 4% in at least one environment, leaving us with 133 diploid 136 populations (43 YPD at 30°C, 48 YPD + acetic acid, and 42 YPD at 37°C) and 20 haploid 137 populations (153 populations total). However, we note that the results of all our analyses are very 138 similar when we consider the entire dataset with outliers included (see Figure Supplements).

140
Adaptation to the home environment leads to consistent fitness gains and pleiotropic effects 141 While there is modest variability between replicate populations, adaptation in each environment 142 leads to a consistent increase in fitness in that "home" environment ( Fig. 2, subplots with bold 143 black borders). As observed in earlier experiments, this fitness increase is largely predictable, 144 and follows a characteristic pattern of declining adaptability: early rapid fitness gains that slow 145 down over time (Couce and Tenaillon 2015). This declining adaptability trend is less obvious 146 among populations evolved at 37°C, possibly because the fitness gains in this environment were 147 generally minimal, but we do observe declining adaptability in the handful of diploid populations 148 at 37°C that experienced larger-than-average fitness gains.

150
Adaptation in each evolution environment also led to fitness changes in most other environments 151 (Fig. 2). In general, these fitness changes tend to have a consistent direction over time for each 152 environment pair. For example, populations adapted to YPD + acetic acid and YPD at 37°C 153 steadily gained fitness in the YPD at 30°C and YPD + 0.4M NaCl environments over time, with 154 the average fitness across populations largely following the same trend seen at home: initial rapid 155 fitness gains followed by slower increases over time. In other instances, fitness gains at home 156 correspond to fitness declines in away environments. For example, populations evolved in YPD 163 To visualize how these pleiotropic effects change over time, we plot these fitness trajectories 164 across pairs of environments (Fig. 3). This representation of the data shows clear but sometimes 165 subtle differences in patterns of pleiotropy depending on evolution environment and ploidy. For 166 instance, while almost all populations gained fitness in both YPD at 30°C and YPD + NaCl, the 167 dynamics of fitness change differed based on evolution environment: populations evolved at 168 37°C (orange lines in Fig. 3) initially made substantial fitness gains in YPD + NaCl sometimes 169 followed by more significant gains in YPD at 30°C, whereas the populations evolved in YPD at 170 30°C (cyan lines) and YPD + acetic acid (green lines) only gained substantial fitness in YPD + 171 NaCl after initial fitness increases in YPD at 30°C (Supplementary File 1). Separately, plotting 172 fitness in YPD + acetic acid against fitness in YPD at 21°C reveals trajectories that segregate not 173 only by evolution environment, but also by ploidy (Supplementary File 1).

175
Characteristic environment-and ploidy-specific pleiotropic profiles emerge over time 176 To understand the diversity of fitness trajectories across environments, we treated the fitness of 177 each population across all five assay environments as a single "pleiotropic profile." We then 178 conducted principal component analysis across all these pleiotropic profiles to characterize 179 variation between replicate populations, across different evolution environments, and over time.   To more formally quantify the emergence of characteristic pleiotropic profiles over time in 202 Figures 4A and B, we developed a simple clustering metric, which counts how many of a given 203 population's five nearest neighbors belong to the same evolution condition on average. We see 204 that the degree of clustering in this two-dimensional space rises appreciably until the 600-205 generation mark, at which point it plateaus (Fig. 4D). The observed clustering from generation 7 200 onward is much greater than expected by chance, as is clustering for the total-data PCA 207 shown in Figure 4B (compared to a null expectation constructed by randomly permuting the 208 evolution condition assigned to each population; p < 0.001). Note that this trend is consistent  We find that these patterns of variability are structured, with specific evolution conditions

252
The role of stochasticity and temporal shifts in pleiotropic dynamics also can be seen in the 253 relative non-monotonicity of fitness trajectories in away environments compared to home 254 environments. To assess non-monotonicity, we interpolated fitness at 500 generations for each 255 population in each assay environment and compared the 0-to-500-generation and 500-to-1000- Despite these characteristic patterns, we also observe ample variability within these trends.   , it is therefore likely that fluctuations on longer timescales (e.g., longer 308 than 600 generations in this system) will lead to qualitatively different outcomes than 309 fluctuations on shorter timescales. Our data show that both the average and variance in these 310 outcomes will also depend critically on the specific sequence of environments experienced by a 311 population.

313
These results underscore the need for further empirical and theoretical work to understand   primer sequences). This 20 µL 10-cycle PCR reaction was performed using Q5 polymerase (NEB 453 M0491L) following the manufacturer's guidelines, using 10 µL (~250 ng) of gDNA as template, 454 annealing at 54°C, and extending for 45 seconds. The first-round PCR products were then purified 455 using one equivalent volume of DNA-binding beads (Aline Biosciences PCRCleanDX C-1003-5) 456 and eluting in 33 µL 10 mM Tris pH 8.5. In the second-round PCR, the remainder of the Illumina 457 adapters and sample-specific Illumina indices were appended to the first-round PCR products (see 458 Supplementary Table 5 for primer sequences). The second round PCR was performed using Kapa µL reaction, using 17.25 µL of first round PCR product, annealing at 63°C and extending for 30 461 seconds for 26 cycles. The second-round PCR products were then purified using one equivalent 462 volume of DNA-binding beads and eluting in 33 µL 10 mM Tris pH 8.5. Following bead cleanup, 463 the concentration of the PCR products was quantified using the Accugreen High Sensitivity  Once fastq files were concatenated, barcode information was extracted as described below. 10, which were longer than the indices associated with epochs 0, 2, and 4, we allowed up to 3 478 mismatches. 479 Then, as with the barcode association mapping, we used a previously described "deletion-error-  We removed all apparently cross-contaminating reads from our analysis.

489
Then, we summed reads associated with all barcodes in a given population, since some 490 populations contained more than one unique barcode (or, in the case of diploids, more than two 491 unique barcodes). In addition, some barcodes were present in the BFAs that could not be 492 confidently assigned to a single well, representing 0.3% of all reads. These were summed 493 together and retained in the dataset. If not, these intervals were excluded from subsequent analysis.

504
To determine s values for each population in each environment at each generation, interval-505 specific s estimates were averaged. Then, s estimates from each of the two technical replicates 506 were averaged, producing a final s estimate. The standard error of this final s estimate was 507 calculated from the two technical replicate s estimates.

508
To clarify our downstream analyses, we excluded 19 outlier diploid populations whose ancestors 509 differed from the mean ancestral fitness by at least 4% in at least one environment. We believe   To triangulate the position of each barcode across the eight plates, for each error-corrected 557 barcode that appeared in the sequencing data, we tabulated which barcodes were present in 558 which libraries, and how many reads were associated with each barcode in each library. These 559 data allowed us to determine the wells in which barcodes belonged.

561
Fitness variability was examined by plotting box-and-whisker plots of population mean fitness 562 values, where the line, box, and whiskers represent the median, quartiles, and data within 1.5xIQR 563 of each quartile, respectively, and outlier populations beyond whiskers are shown as points (Fig.   564   5A). To compare the resulting IQR for various evolution conditions and fitness assay 565 environments, 95% confidence intervals of the IQR were calculated from bootstrapped interval-566 specific replicate s measurements (Fig. 5B). 567 To evaluate whether home environment fitness variance was less than away environment fitness 568 variances at each evolution timepoint, we applied a Brown-Forsythe test (Brown and Forsythe 569 1974). Since this test is typically a two-tailed test, and we wanted instead to employ a one-tailed 570 test, we used the z scores from the Brown-Forsythe test to arrive at a two-tailed t-statistic. We   (Figs. 4A,B). The clustering metric plotted in Fig. 4D is the number of five 584 nearest neighbors that belong to the same evolution condition as the focal population, averaged for 585 each evolution condition. Error bars represent 95% confidence intervals of the mean clustering 586 metric, which were calculated by performing the PCA and clustering analysis on bootstrapped 587 interval-specific replicate s measurements. The null expectation for populations to cluster by 588 evolution condition was computed by permuting the evolution condition 1000 times and 589 performing the clustering analysis as described above. The permuted clustering metrics were then 590 compared to the true mean clustering metric by a two-sided Student's t-test (using the Scipy.stats 591 ttest_ind_from_stats function).

592
Non-monotonicity analysis 593 To assess non-monotonicity, we linearly interpolated fitness at 500 generations for each 594 population in each assay environment. We achieved the interpolated standard errors in fitness by top-left or lower-right quadrant --corresponding to an increase followed by a decrease, or a 607 decrease followed by an increase, over the 1000-generation experiment --these were considered 608 to be "clearly non-monotonic." We applied a χ2 test to evaluate the significance in the difference 609 in the frequency of non-monotonicity in home versus away trajectories.        using all fitness data from across the 1000 generations. The first two PCs are plotted and explain 30% and 22% of the variance, respectively. (C) Plots of fitness trajectories in all 5 assay environments for 8 example populations (a-h, identified as points in (B)). (D) Population clustering in PCA by evolution condition over time. Clustering of each population was quantified as the number of five nearest neighbors that share the same evolution condition, for each 200-generation interval, and across all intervals. Clustering metrics were averaged for each evolution condition to calculate point estimates; error bars represent 95% confidence intervals of the mean clustering metric, estimated by performing PCA on bootstrapped replicate fitness measurements.      Figure 4A.  Figure 4B.     Brown-Forsythe significance test results for differences between variance at home and away. Each panel shows, for each of the 5 assay environments, the change in fitness over the first 500 (x-axis) and second 500 (y-axis) generations of evolution of each population in a given evolution environment. Populations that fall in shaded quadrants have trajectories that are non-monotonic. Points corresponding to fitness in the home environment are colored more opaquely than points corresponding to fitness in away environments, and panel borders have been colored to match the home environment. Fitness at generation 500 has been interpolated. (B) Each panel corresponds to a given evolution environment and shows the proportion of populations evolved in that environment that exhibit clearly nonmonotonic fitness trajectories in (A). "Clearly non-monotonic" trajectories are those populations (points) in (A) that fall in the grey quadrants and whose error bars (1 standard error in either direction) do not span either the x-or y-axis. As in (A), bars corresponding to the home environment are colored more opaquely than bars corresponding to away environments.         Figure 4D quantified for (A) ten and (B) three nearest neighbors. Clustering metrics were averaged for each evolution condition to calculate point estimates; error bars represent 95% confidence intervals of the mean clustering metric, estimated by performing PCA on bootstrapped replicate fitness measurements.

. Non-monotonicity in evolutionary trajectories for unfiltered data (outliers included). (A)
Each panel shows-for each of the 5 assay environments-the change in fitness over the first 500 (x-axis) and second 500 (y-axis) generations of evolution of each population in a given evolution environment. Populations that fall in shaded quadrants have trajectories that are nonmonotonic. Points corresponding to fitness in the home environment are colored more opaquely than points corresponding to fitness in away environments, and panel borders have been colored to match the home environment. Fitness at generation 500 has been interpolated. (B) Each panel corresponds to a given evolution environment and shows the proportion of populations evolved in that environment that exhibit clearly non-monotonic fitness trajectories in (A). "Clearly non-monotonic" trajectories are those populations (points) in (A) that fall in the grey quadrants and whose error bars (1 standard error in either direction) do not span either the x-or y-axis. As in (A), bars corresponding to the home environment are colored more opaquely than bars corresponding to away environments. As with the outliers-excluded data, populations exhibit clearly non-monotonic trajectories in away environments much more commonly than in home environments (p < 0.0001), with most of these reflecting initially positive pleiotropic effects.