Rapid evolution leads to differential population dynamics and top‐down control in resurrected Daphnia populations

Abstract There is growing evidence of rapid genetic adaptation of natural populations to environmental change, opening the perspective that evolutionary trait change may subsequently impact ecological processes such as population dynamics, community composition, and ecosystem functioning. To study such eco‐evolutionary feedbacks in natural populations, however, requires samples across time. Here, we capitalize on a resurrection ecology study that documented rapid and adaptive evolution in a natural population of the water flea Daphnia magna in response to strong changes in predation pressure by fish, and carry out a follow‐up mesocosm experiment to test whether the observed genetic changes influence population dynamics and top‐down control of phytoplankton. We inoculated populations of the water flea D. magna derived from three time periods of the same natural population known to have genetically adapted to changes in predation pressure in replicate mesocosms and monitored both Daphnia population densities and phytoplankton biomass in the presence and absence of fish. Our results revealed differences in population dynamics and top‐down control of algae between mesocosms harboring populations from the time period before, during, and after a peak in fish predation pressure caused by human fish stocking. The differences, however, deviated from our a priori expectations. An S‐map approach on time series revealed that the interactions between adults and juveniles strongly impacted the dynamics of populations and their top‐down control on algae in the mesocosms, and that the strength of these interactions was modulated by rapid evolution as it occurred in nature. Our study provides an example of an evolutionary response that fundamentally alters the processes structuring population dynamics and impacts ecosystem features.


Clonal identities and contamination
The 36 clones used in the present experiment were the same as used in Stoks et al. (2016), however, we detected that some contamination had occurred when we screened the clones for their genotype using 9 microsatellite markers (see above). For the pre-fish population our multi-locus genotype identification suggests that one clone (B7) was accidently replaced by another clone from the same population (B9). For this population, the contamination only occurred within the population and resulted in 11 instead of 12 clones being used in the experiment ( Figure SI1 a & b). For the high-fish population, two clones (M10 and M12) were accidently replaced by another clone from the same population (M11), one clone (M3) was accidently replaced by a clone with the same multilocus genotype as clone B6 of the pre-fish population, and one clone (M2) was accidently replaced by a clone with the same multilocus genotype as clone T12 from the reduced-fish population. So for this population, 8 clones were inoculated in addition to two clones that was derived from another population (Figure SI1 c & d). In the reduced-fish population, multi-locus genotype identification suggests that one clone (T3) was accidently replaced by a clone from the high-fish population (M7).
For this population, 11 clones were used in the experiment in addition to one clone from another population ( Figure SI1 e & f). As we screened the lineages for their genotypes after inoculating the mesocosms, we could not prevent the contamination. Although the contamination was unfortunate, only the among-population contaminations (one clone in the reduced-fish population and two clones in the high-fish population) can potentially interfere with the interpretation of our results, which are based on comparisons among populations. Moreover, to the extent that the interpretation of our results depends on the detection of among-population differences, the contamination does not induce false positives, but rather results in a conservative assessment of the impact of evolutionary change. The average detected relative abundance of the pre-fish clone in the high-fish population treatment at the end of the experiment was 0% in both Control and Predation treatments ( Figure SI1 c & d), which translates into at most very low abundances of this clone. The average observed relative abundance of the reduced-fish population clone in the high-fish population treatment at the end of the experiment was 1% in the Control mesocosms and 26% in the Predation mesocosms (Figure SI1 c & d). The average observed relative abundance of the high-fish population clone in the reduced-fish population treatment at the end of the experiment was 38.1 % in the Control mesocosms and 28.7% in the Predation mesocosms ( Figure SI1 e & f).

SI2. Data analysis Empirical dynamic modeling
Empirical dynamic modeling (EDM), uses time-series to reconstruct the attractor manifold (see further) and allows for the exploration of the mechanisms underlying the dynamics of the system (Deyle, May, Munch, & Sugihara, 2016;Sugihara, 1994;Sugihara et al., 2012;Sugihara & May, 1990).
Simplex projections are a forecasting technique from the EDM framework (Sugihara & May, 1990).
The forecast skill of simplex projections using one group of time-series as a library (i.e. learning set) to make forecasts for data points in another time-series can be used to assess the similarity in the attractor manifold of those time-series. The S-map method is a technique from the EDM framework that can be used to estimate interaction strengths. S-maps do so by recovering the Jacobians (i.e. partial derivatives) at each time-point (Deyle, May, et al., 2016;Sugihara, 1994). Another technique from the EDM framework is convergent cross mapping (CCM) (Sugihara et al., 2012). CCMs can be used to identify causal links in the system. A brief explanation of EDM and the techniques that we used are given below, together with details of our implementation of them. For more in-depth explanations and further examples, we refer the reader to the rEDM user guide Ye, Clark, Deyele, Keyes, and Sugihara (2016) and empirical dynamic modeling for beginners (Chang, Ushio, & Hsieh, 2017).
Dynamic systems are often described as a set of multiple equations, in which each equation describes how the dynamics of a certain variable depends on itself and other variables. When different states of a system are very similar, the state of the system will over the short run evolve very similarly. Representing all the relevant variables of a system as a set of Cartesian coordinates in state space and the observations in the time-series of these variables as coordinates visited by that system, results in a collection of trajectories forming a geometric object called an attractor manifold.
The attractor manifold is the product of the specific rules and equations that describe the interactions between variables of the system, and thus is an empirical description of the dynamics of the system. An animation explaining the reconstruction of the attractor manifold by its variables is given in Sugihara et al. (2012; https://youtu.be/8DikuwwPWsY). As similar states evolve similarly over the short run, so do nearby states in state space. Thus, when time-series of the relevant variables are available, short-term forecasts for a given state can be made based on the predicted short-term future of nearby states in state space. However, when time-series for some variables in the system are not available, trajectories cross and nearby states will not go in exactly the same direction. In reality there might be countless variables influencing every system, but often the majority of the changes over time in a certain state variable are caused by only a few other variables.
Hence, relatively skillful forecasts can be made based on these few relevant variables. It is, however, not always feasible to obtain measurements of or know all relevant variables of the system. Takens (1981) addressed this problem by using the fact that in a dynamic system, time-series of a variable that is influenced by other variables also contain information on these variables. An everyday used example of this principle is our ability to estimate the future location of moving objects by using consecutive snapshots of these objects, locations, without directly observing the momentum of the object and the forces (e.g. gravity) acting on it. Thus, instead of representing the state of the system as a vector (i.e. a multidimensional point) in state space with as coordinates the relevant state variables, one can use time-lagged observations (i.e. snapshots) of one variable as the coordinates. This is called a time-lagged embedding. There are a minimum number of time-lagged observations needed to capture all the necessary information and thereby prevent trajectories from crossing in the time-lagged embedding (i.e. for the embedding to be diffeomorphic). Before Takens' theorem, it was not clear if the number of time-lags needed would be too high for any practical usage. Takens (1981) found a connection between the number of required lags, i.e. the embedding dimension E, and the number of state variables, i.e. the number of dimensions D. He demonstrated that a timelagged embedding using just E = 2 * D + 1 lags is the maximum needed to obtain a diffeomorphism of the original attractor manifold of a dynamic system (i.e. to prevent lines crossing in the embedding).
Thus, if the relevant variables of a system are the two variables X and Y, no more than five lags are needed in the embedding (i.e. {Y(t), Y(t-1), Y(t-2), Y(t-3), Y(t-4)}). This means one can obtain a shadow version (i.e. a time-lagged embedding) of the original attractor manifold by using only a few lags of one variable. Although this shadow manifold is a globally distorted (e.g. stretched or bent) version of the original manifold, this distortion is a smooth invertible change in coordinates. A short animation by Sugihara et al. (2012) explaining Takens' theorem can be found here: https://youtu.be/QQwtrWBwxQg. The same points in time that are close on the original attractor manifold are also close on the shadow manifold. Thus, as the state of the system changes over time and visits different neighborhoods on the attractor manifold, it will pass by neighborhoods on the manifold that it has visited before and the time points in history close to each other on the original attractor are also close on the shadow manifold. In (univariate) simplex projections, this fact is used to make forecasts using only one variable (Sugihara & May, 1990). As differences in dynamics between different populations can be more pronounced in some variables (i.e. dimensions in state space) than others, we decided to test for significant differences between populations independently for different variables. For this test we were thus able to use univariate simplex projections.
An interesting consequence from Takens' theorem is that when a variable X influences another variable Y, then X will leave its mark on the dynamics of Y. Time-points close on the time-lagged embedding of X will also be close on the time-lagged embedding of Y. Thus, if nearby points on the shadow manifold of Y are also close on the shadow manifold of X, than X likely influences variable Y. This is the basis for convergent cross mapping (Sugihara et al., 2012). An animation explaining this principle by Sugihara et al. (2012) can be found here: https://youtu.be/NrFdIz-D2yM.

Simplex projections
We used univariate simplex projections to compare the similarity in dynamics of individual variables between populations. In simplex projections, for each state in the time-lagged embedding of the prediction set (i.e. for each target state x(t*)), the E+1 closest points on the shadow manifold of X, which is reconstructed using only the library set, are taken and a weight is calculated based on their Euclidian distance from x(t*). To make forecasts with a specified time step, t p , into the future, the states that the E+1 nearest states have t p into the future are multiplied with their respective weights and the average of these products is used as the forecast (Sugihara & May, 1990). If the time-series used in the library can be used to make skillful forecasts of the time-series in the prediction set, then the dynamics underlying the time-series in the prediction set must be similar to those in the timeseries used in the library. We expressed forecast skill in MAE (mean absolute error). For each population, to determine whether the dynamics of replicates of the same population are more similar than between different populations, we used a one sided Mann-Whitney test comparing the skill of within population forecasts to forecasts from other populations. We made forecasts for all possible combinations of (replicate) time-series as library and prediction sets, with 3 time-series as library predicting one other time-series (i.e. excluding combinations where the same time-series occurred in both library and prediction set). This leads to four MAEs for each set of within population forecasts and 16 for the between population forecasts. For simplex projections the number of lags used in the reconstructed state space has to be specified (i.e. the embedding dimension E). We here each time tried all embedding dimensions below 7 and used the E that resulted in the highest value.
Another parameter that has to be set in simplex predictions is the forecast time step t p . The standard choice for this is to use the smallest step (i.e. 3 days for Daphnia time-series and 1 day for phytoplankton). However, when dynamics are relatively slow, small time steps might not be sufficiently challenging to distinguish in forecast skill. In contrast, as is characteristic to non-linear systems, large forecast time steps result in lower predictability and forecast skill decreases in all models. Therefore, in testing the difference between populations, we tested different forecast time-steps (3, 6, 9, 12 and 15 days). Although within population simplex projections might not be significantly more skillful (i.e. lower MAE) with very small time-steps or too large time-steps, if they significantly differ from another population, than at least within a certain range of forecast time steps, they would do significantly better when using libraries from the same population than when using libraries from a different population. In Simplex projections, unlike the other EDM methods we applied, each time-series is analyzed separately. Hence, we were able to make use of the finer resolution of the phytoplankton time-series and use a time lag of one day between the time lags in the embedding. The P-values of these tests against the forecast time-step used are shown in figure   SI8 and the number of significant p-values (<0.05) among the tests using the five different time steps, are shown in Table SI1 in the main text.

Cross mapping
For each time point t p , the E+1 closest points on the shadow manifold of Y are taken and a weight is calculated based on their distance from t p . Multiplying the E+1 nearest states with their respective weights and averaging them gives the value of Y at t p . The same E+1 time points, but from the timeseries of X are then multiplied with these weights to make a prediction. Note that these do not have to be the E+1 closest points on X as well, as long as they are close enough, the prediction will be reasonably skillful. The Pearson correlation coefficient between the predicted values of X based on the manifold of Y at each time point and the true value of X at those time points is the prediction skill of the cross map (ρ ccm ). In the presence of a causal link between the considered variables, the prediction skill will increase until they converge to a certain ρ ccm when more time-points are used to make the shadow manifold of the library variable (in our case Y), (Sugihara et al., 2012).
Given that each selection of days from the time-series to use in the library will result in a different cross map skill (ρ ccm ), the days of the time-series used at each library size is drawn multiple times (in our case 100 times) resulting in a distribution of cross map skills (ρ ccm ). In all our analysis, we combined time-series data of the four replicates of each population within each treatment, but in such a way that no one vector (i.e. set of time lagged observations) contained data points from different replicates (Clark et al., 2015;Hsieh, Anderson, & Sugihara, 2007). To avoid problems with overfitting, we first performed CCMs predicting the state of X three days before based, on the shadow manifold of Y (see Deyle, Maher, Hernandez, Basu, & Sugihara, 2016;Deyle, May, et al., 2016) and then used the embedding dimension resulting in the highest ρ ccm from this as embedding dimension for the actual CCMs. In the actual CCMs, the shadow manifold of Y was used to predict the state of X on the same day (see Figure SI10). False signs of positive cross maps were further eliminated by testing against randomly generated surrogate time-series that were used as a null distribution (Deyle, Maher, et al., 2016;Deyle, May, et al., 2016). We generated 500 surrogate timeseries in which the order of the days was randomly permutated. If the cross map skill is significantly greater in the original time-series than in the surrogate time-series based null distribution, the properties that were incorporated in the surrogate time-series are not enough to explain the size of the cross map skill. Importantly, we used the same permutation of days for each of the 4 replicates in generating null distributions to also consider potential false signs of causal influences resulting from synchrony of variables with external forces that acted upon the four replicates simultaneously. We determined the embedding dimensions for both the surrogate and the regular time-series similarly to the regular CCMs (i.e. based on the ρ ccm in 3 day backward predictions). When the original timeseries performed better than 95% of the surrogates using ρ ccm as a criterion, we considered them to be significant (i.e. p<0.05). Results of this analysis are given in Figure SI9. All cases with significant surrogate tests ( Figure SI9) showed convergence in the CCM plots ( Figure SI10). An overview of significant interactions between variables is given in Figure 5 in the main text.

S-maps
The closer one zooms into a small neighborhood on the attractor manifold, the more linear the manifold becomes. The S-map method is a locally weighted linear regression scheme (Sugihara, 1994). It approximates the best local linear model at each measured state by giving more weight to states on the attractor manifold that are more nearby that state. Similar to a multivariate linear regression, S-maps average out noise by using all data points, rather than just a few neighboring points in state space. In contrast to multivariate linear regression, S-maps allow points more closely located on the manifold to the target point x(t) to be given a higher weight in the forecast, thereby accounting for potential state-dependent-differences in interaction strengths over time, which is typical for non-linear dynamic systems. S-maps contain one variable, theta (θ), which sets the degree  S-maps can be used on both univariate embeddings of one variable or using a multivariate embedding (Deyle & Sugihara, 2011). When using a multivariate embedding, the S-map coefficients, i.e. the regression coefficients of each locally weighted linear regression (which are equivalent to the partial derivatives on the manifold or the Jacobian) give the strengths and signs of the interactions between variables (Deyle, May, et al., 2016).
The theta resulting in the best forecast skill ρ was used for each population, treatment and forecasted variable. The forecast skill ρ was found to be always significantly (<0.0001) better than zero using Fisher's z-transformation (Table SI2). The resulting S-map estimates of interaction strengths (i.e. the elements of C) are plotted against time in Figure SI6 and SI7.
The effect of X on the future (3 days later) of Y is given by But in the non-linear forecasts by the S-maps, for each forecast different parameters are estimated because of the local weighting, and the ! ! , ! ! and ! ! become the elements of the Jacobian matrix: Non-linear systems, as analyzed by S-maps, can thus be interpreted as linear systems, but with changing parameters depending on the position in state space of the system. By interpreting interactions this way, the straightforward intuition we have with linear systems can be extended to explore non-linear systems, the difference being that the interaction strengths change over time depending on the states. So we see a cloud of interaction strengths in Figure 3, Figure 4, Figure SI4 and Figure SI5 that shows how interaction strengths vary with state. By plotting the interaction strengths estimated at each time point against measured observations at the same time points of an important variable (e.g. phytoplankton in Figure 3, Figure 4a-d, Figure SI4 and Figure SI5) we are able to see how the interaction strengths from the locally weighted linear regressions predicted by the model depend on this variable. To test for the significance of relations between estimated interaction strengths and variables, we performed several linear regressions (Figure 4). The statistics of these linear regressions are shown in Table SI2. All three variables capture some aspects of juvenile energy content, but should be considered loose approximations rather than precise estimates. We plot neonate body size because of its link to body mass. We plot average number of offspring assuming that investment per individual offspring would become lower as their numbers increase. This assumes an equal amount of available energy of all mothers. The third index (average size at maturity divided by average fecundity) tries to take variation in energy content of mothers into consideration by scaling the number of offspring to the size of the mother.  Table SI1. Degree of non-linearity (theta) resulting in optimal S-map forecast skill (i.e. highest ρ) and probability that the ρ is greater than zero using Fisher's z-transformation.  Table SI2. Results of the simple linear regressions in Figure 4. The S-map estimates of the effect of Daphnia on phytoplankton were regressed on chlorophyll a concentration (log(chla)) and (only for juveniles) the ratio of Adult Daphnia : chlorophyll a (Adult/log(chla)) for pre-fish, high-fish and reduced-fish in the Control and Predation treatments. Linear regressions were only performed for all population x treatment combinations when the interactions were significant (p<0.05) in the CCM test.  Predation treatment at the end of the experiment. Clones codes that start with a letter "T" belong to the reduced-fish population ("Top"), codes starting with "M" indicate clones that belong to the high-fish population ("Middle") and "B" indicates clones from the pre-fish population ("Bottom"). In the prefish population clone B7 was accidently replaced by clone B9 from the same population. In the high-fish population clones M8 and M9 are the same multi-locus genotype (cannot be discriminated with the marker set used), clones M10 and M12 were accidently replaced by clone M11 from the same population, clone M3 was accidently replaced by a clone with the same multilocus genotype as B6 of the pre-fish population and clone M2 was accidently replaced by a clone with the same multilocus genotype as clone T12 from the reduced-fish population. In the reduced-fish population clone T3 was replaced by clone M7 from the high-fish population. Error bars indicate standard error.        Figure SI10. Convergent cross-maps for each population in both treatments. Each plot shows the cross map between two variables in both directions against the length of the library set used to make the cross map. At each library length 100 random samples were used to make cross maps. The solid line shows the average cross map skill (ρ ccm ) and the dashed line shows the standard deviation. If variable X influences variable Y, the skill of the cross map from Y to X should initially increase and then converge to an upper limit, as the length of the library set used is increased.