The evolutionary maintenance of sexual reproduction has long challenged biologists as the majority of species reproduce sexually despite inherent costs. Providing a general explanation for the evolutionary success of sex has thus proven difficult and resulted in numerous hypotheses. A leading hypothesis suggests that antagonistic species interaction can generate conditions selecting for increased sex due to the production of rare or novel genotypes that are beneficial for rapid adaptation to recurrent environmental change brought on by antagonism. To test this ecology-based hypothesis, we conducted experimental evolution in a predator (rotifer) - prey (algal) system by using continuous cultures to track predator-prey dynamics and in-situ rates of sex in the prey over time and within replicated experimental populations. Overall, we found that predator-mediated fluctuating selection for competitive versus defended prey resulted in higher rates of genetic mixing in the prey. More specifically, our results showed that fluctuating population sizes of predator and prey, coupled with a trade-off in the prey, drove the sort of recurrent environmental change that could provide a benefit to sex in the prey, despite inherent costs. We end with a discussion of potential population genetic mechanisms underlying increased selection for sex in this system, based on our application of a general theoretical framework for measuring the effects of sex over time, and interpreting how these effects can lead to inferences about the conditions selecting for or against sexual reproduction in a system with antagonistic species interaction.

Experimental evolution: We conducted an experimental evolution study that used continuous cultures (chemostats) for evolving replicate (n = 3) algal populations with and without (control) predators in order to follow predator-prey dynamics, predation intensity (predator/prey ratios) and the rate of sex in algal prey populations over time (Fig. 1A-B). Each chemostat had a volume of 850 mL, was continuously stirred, ventilated, and fresh media [80 μM nitrogen] supplied at a constant flow-through rate of 0.11 d^-1. Cultures were maintained at 26°C. To allow zygospore maturation, we set the cultures to an 8h L:16h D cycle. We inoculated all 6 chemostats with 1x10⁶ cells of each mating type, followed by 300 rotifers in 3 of the chemostats 3 days later. Chemostat sampling began 4 days after rotifer inoculation.

Predator-prey population dynamics, predation intensity, in-situ rate of sex, mating type allele frequencies:  The experimental evolution study lasted for 62 days, and each chemostat was subsampled (1) every other day (with a 9-day gap between the last two sampled timepoints) by extracting 50 mL samples to enumerate predators, prey and zygospores using a Neubauer hemocytometer and compound microscope (Fig. 1B; Fig. 2A-F), (2) every week to track the frequency of mating types over time in each algal population using ddPCR (Fig. 1C; Fig. 4A-F), and (3) every 4 days by streaking out subsamples onto agar plates to isolate algal clones for subsequent fitness assays to test for changes in selection and the effects of sex (Fig. 1D; Fig. 3; Fig. 6). Zygospore counts are used to estimate the in-situ rate of sex in the prey (Fig. 2A-F). Predation intensity was calculated as the ratio of predator/prey densities and scaled to the maximum value for each corresponding population of the predation treatment (Fig. 2G-I). 

Assays measuring changes in selection and the effects of sex: From the predation treatment, we randomly isolated 10 asexual clones from each of the 15 stored algal populations per replicate chemostat and kept them on agar plates in dim light. To obtain sexual clones from each of the corresponding populations, we transferred the populations from the agar plate to nitrogen-free media to induce sex and random mating. Mating reactions form sticky pellicles of aggregated zygotes that can be easily separated from the liquid and unmated cells. We plated the pellicles and followed standard protocol in carrying out the sexual cycle. After germination was complete, we were able to streak out progeny to obtain 10 random, sexually-produced clones. The 20 asexually and sexually-produced offspring served as a representative population subsample from each timepoint. 

Before each assay, algal clones were washed to remove excess nitrogen and pre-conditioned on media with 80 μM nitrogen for approximately 12 hours. Fitness assays lasted 24 hours and were conducted in 24-well plates containing 1.5 mL media with 80 μM nitrogen. Well-plates were continuously illuminated and each well contained a 1 mm stir bar that was set to stir for 1 minute every hour. We used a Tecan Infinite 200 PRO plate reader to obtain optical density measurements (absorption wavelength of 680 nm) and calculate clone growth rates (day^-1). For all clones and assays, we obtained 3 technical replicates per measurement. 

To confirm changes in selection for competitive and defended prey over time, we used the absolute growth rate (no predation) and defense data. For the assays measuring changes in selection, each well was inoculated with an algal suspension of ~500 µL to obtain a starting OD of 0.1 and either 0 (no predation = growth rate) or 10 rotifers (defense). We estimated algal defense indirectly by counting the number of rotifers present in each well at the end of the assay and calculated the defense as 1/total number of rotifers. For plotting purposes, we used reciprocal values (1/x) for the defense data in Fig. 3A-C.

Estimating the short and long-term effects of sex requires the measurement of fitness of sexually and asexually-produced offspring in the environment that they were produced in. Therefore, to account for fluctuations in selection for defense and growth, we measured algal growth rates under low predation (‘L’) or high predation pressure (‘H’), which simulates the selection conditions that the clones evolved under. We used prey growth rates, with growth rates in the (L)ow predation environment representing selection for competitive ability and growth rates in the (H)igh predation environment representing defense, because prey defense was measured indirectly through predator fitness. For this reason, we chose 1 trait which was directly obtained from the prey in order to evaluate the effects of sex in the prey. At each timepoint for each replicate population, whether fitness was measured in the low or high predation environment was selected based on the predator-prey ratios (Fig. 2G-I; see Table S2). Specifically, we assigned timepoints with high/increasing predator densities and high predation intensity as times selecting for defense in the prey, and timepoints with low/decreasing predator densities and low predation intensity as times selecting for prey growth (i.e. competitive ability). We used these designations to determine from which assay environment (L, H) data would be plotted for the short and long-term effects of sex (Fig. 6; Table S2: Selection environment). For the assays measuring the short and long-term effects of sex, each well was inoculated with an algal suspension of ~500 µL to obtain a starting OD of 0.1 and either 5 (low predation) or 10 (high predation) rotifers.

Estimating rates of sex: To quantify sex, we measured (1) in-situ rates of sex, (2) allele frequencies of mating types, and (3) conducted a test for evolutionary changes in the rate of sex. The first two measurements were done for all control and experimental populations (n = 6), while the last test was only conducted for the prey populations within the predation treatment. Zygospore counts and changes in the frequency of mating types over time were used as estimates of the in-situ rate of sex, with the latter using the weekly population subsamples and digital droplet polymerase chain reaction (ddPCR) (see Supplementary Methods) (Koch et al. 2015). We conducted an additional assay under standardized conditions testing for the time that is required for the production of a mating reaction (pellicle formation) comparing algal populations from day 62 to the ancestors (see Supplementary Methods). We used the time that is required for the production of a mating reaction (pellicle formation) based on the assumption that the ancestor switches only in a nitrogen-free environment or when all nitrogen is consumed. In contrast, the genotypes isolated from the predation treatment, which likely evolved to switch in the absence of the trigger, were predicted to switch earlier as compared to the ancestors. Pellicle assays can be used as an end-point detection method for quantifying a rate of sex in Chlamydomonas because gametes are too difficult to distinguish, visually, from vegetative cells (Harris 2009). Pellicles are the result of aggregated zygotes forming a film on the air-water surface and are easily seen by the naked eye, which eliminates the need to continuously sample and disturb the mating reaction to look for gametes or zygotes. Rotifers were not included in this assay.

Pellicle assay: From timepoint 15 (end of experiment, day 62), subsamples of each replicate population of the predation treatment (n = 3) were grown in separate culture flasks containing 80 mL media with 800 μM nitrogen for 5 days on a linear shaker in bright light. The same was done for the ancestors. Cultures were concentrated by centrifuging at 3000 rpm for 10 minutes at room temperature. Resulting algal pellets were washed 3 times with media containing 0 μM nitrogen. Half of each pellet was re-suspended in 4 mL of media containing 0 μM nitrogen and the other half in 4 mL media with 80 μM nitrogen. Since nitrogen deprivation is the sexual cycle cue, mating reactions with 0 μM nitrogen served as the controls and those with 80 μM nitrogen represent the evolved environment. Each algal suspension was inoculated into a 24-well plate with 1 mL per mating reaction and 4 technical replicates per treatment and nitrogen concentration. Rotifers were not included. Mating reactions were left stationary in bright light until every mating reaction was complete, as evidenced by the formation of a pellicle, which is biofilm-like layer of aggregated zygotes on the liquid surface. Every 30 minutes, mating reactions were poked with a pipette tip to test for the presence of the sticky pellicle, at which point the time was recorded. We followed the protocol of Werner and Mergenhagen (1998) for using PCR to confirm the presence of both mating types in each mating reaction of the third chemostat replicate that did not produce a pellicle. Since in the control populations one mating type took over, we could not conduct a pellicle assay for these populations nor include them in this analysis.

Selection differentials: To estimate selection differentials, we first plotted a linear regression of the phenotypic data (fitness: growth rate versus defense) for each of the 20 clones to obtain a slope for each timepoint (n = 15) and replicate population (n = 3). These selection differentials (absolute values) were then plotted over time to show the changes in selection for each replicate chemostat, along with the frequency of zygospores (in-situ rate of sex). 6 outliers (out of 900 values) were removed. 

Statistical analyses: All statistical analyses were conducted using R (R Core Team 2015), version 3.1.3. To estimate the period of the cycles in the algal prey populations (Figs. 2D-F; S2), predator-prey ratios (Figs. 2G-I; S2), changes in selection (Fig. 3A-C), and the short and long-term effects of sex (Figs. 6A-C; S3), we used wavelet analyses (WaveletComp package in R (Roesch 2014)).  For all wavelet analyses, we followed standard time-series analysis practices and used detrended (pracma package (Borchers 2015)) and smoothed time-series data (spline function in R, excluding day 62). We excluded day 62 because there was a 9-day gap between the last 2 sampled timepoints. Generally, wavelet analysis decomposes a time-series into time/frequency space simultaneously and provides information on both the amplitude of any periodic signal within the series and how this amplitude varies with time. The wavelet analysis used here calculates the wavelet power by applying the Morlet wavelet. Significances of periodicity were assessed by testing the null hypothesis that a period is not relevant at a certain time of the time-series by using a simulation algorithm representing white noise (default methods in the WaveletComp package (Roesch 2014); we also calculated significance levels for simulations where the time-series was shuffled but there were no differences compared to the default method of white noise). From the wavelet analysis, we estimated the dominant period by using the period when the average wavelet power is maximized and significant at the 0.05 level (wavelet power averages are shown in Figs. S2; S3). We further tested for a significant relationship between oscillations in the predator-prey ratios (Fig. 2G-I) and population mean defense, population mean growth, and the frequency of defended phenotypes (Fig. 3D-F). We did this by identifying the dominant and significant phase shifts between these time series. Next we used wavelet coherence analyses to measure the local correlation between two series over a specific period (WaveletComp package in R (Torrence & Compo 1998; Roesch 2014)). The value of wavelet coherence is ‘0’ when there is no relation between the two oscillations (no phase coupling) and ‘1’ when there is a full correlation (perfect phase coupling) between the two oscillators. We extracted from these analyses the significant phase shifts (at the 0.05 level and outside the cone of influence). Phase shifts are expressed in days and show the number of days the defense, growth or frequency of defended types is lagging behind the predator-prey ratio. Here, a phase shift of ~0 days means that the maxima and minima of the two time series occur at the same time (in-phase), a phase shift of ~5 days means that the maxima of one time series co-occurs in time with the minima of the other time series.

For comparing the controls and predation treatment, we used a non-parametric Wilcoxon-Mann-Whitney U test to test for significant differences in mean population density. We calculated the coefficient of variation for prey population densities using the ‘cv’ function as part of the raster package (Hijmans 2015) and a non-parametric Kruskal-Wallis test to test for significant differences between the controls and predation treatment. To test for significant differences between population mean ranks for asexual and sexual offspring fitness (growth rate) at each timepoint in either predation environment (L or H), we used Wilcoxon rank sum tests. This allowed us to test for significant short-term costs and benefits. To test for significant long-term costs and benefits, we used F-tests to compare variances in fitness (growth rate) of sexual and asexual offspring for each timepoint and respective predation environment (L or H). For our defense assay validation, we used linear models to test for the effect of algal phenotype on rotifer growth rate and on the fraction of total cells remaining after one day.  To test for a significant difference in the evolutionary rate of sex (selection for sex) of evolved prey populations compared to the ancestors under control conditions (0 μM nitrogen) and those experienced in the evolution experiment (80 μM nitrogen), we used a Wilcoxon-Mann-Whitney U test.

There is a 9-day gap between the last 2 timepoints of the 62-day experimetal evolution study (last sampled timepoints: days 53 and 62).