Evolution of population dynamics following invasion by a non‐native predator

Abstract Invasive predatory species are frequently observed to cause evolutionary responses in prey phenotypes, which in turn may lead to evolutionary shifts in the population dynamics of prey. Research has provided a link between rates of predation and the evolution of prey population growth in the lab, but studies from natural populations are rare. Here, we tested for evolutionary changes in population dynamics parameters of zooplankton Daphnia pulicaria following invasion by the predator Bythotrephes longimanus into Lake Kegonsa, Wisconsin, US. We used a resurrection ecological approach, whereby clones from pre‐ and post‐invasive periods were hatched from eggs obtained in sediment cores and were used in a 3‐month growth experiment. Based on these data, we estimated intrinsic population growth rates (r), the shape of density dependence (θ) and carrying capacities (K) using theta‐logistic models. We found that post‐invasion Daphnia maintained a higher r and K under these controlled, predation‐free laboratory conditions. Evidence for changes in θ was weaker. Whereas previous experimental evolution studies of predator–prey interactions have demonstrated that genotypes that have evolved under predation have inferior competitive ability when the predator is absent, this was not the case for the Daphnia. Given that our study was conducted in a laboratory environment and the possibility for genotype‐by‐environment interactions, extrapolating these apparent counterintuitive results to the wild should be done with caution. However, barring such complications, we discuss how selection for reduced predator exposure, either temporally or spatially, may have led to the observed changes. This scenario suggests that complexities in ecological interactions represents a challenge when predicting the evolutionary responses of population dynamics to changes in predation pressure in natural systems.


| INTRODUC TI ON
A direct consequence of globalization is the facilitated movement of species to novel environments (Elton, 2020). Predatory species successfully invading an ecosystem are of particular concern, as they can have devastating effects on native prey species, often causing rapid changes in both population abundance and phenotypic traits (Reznick & Ghalambor, 2001;Strauss et al., 2006;Thompson, 1998).
One obviously relevant aspect in this context is selective predation and resulting evolutionary responses in targeted traits (Åbjörnsson et al., 2004;Melotto et al., 2020). A less obvious aspect is that evolutionary responses in prey may also arise as an indirect byproduct of predator-induced mortality. Such indirect selection may stem from large-scale shifts in community structure and ecosystem function.
For example, a shift in the direction of trophic control in the environment inhabited by the affected prey species may occur, from a state of 'bottom-up control' (i.e. resource limitation) to 'top-down control' (i.e. predator control), or vice versa. The resulting evolutionary responses in traits such as age at maturation, size and number of offspring, and competitive abilities may then be predicted based on density-dependent selection theory (Einum et al., 2008;Mueller et al., 1991;Saether et al., 2016). Indeed, evidence suggests that prey evolutionary responses resulting from such indirect selection can be of similar importance to those resulting from the direct selective mortality (Schmitz et al., 1997;Walsh & Reznick, 2010).
Given that predators may change prey phenotypes it follows logically that this could translate into evolution of their population dynamics. Such changes might include shifts in the intrinsic rate of increase and/or the strength of density dependence. Evolutionary responses in prey population dynamics have been demonstrated by rearing prey under different levels of predation under laboratory conditions (Shertzer et al., 2002;Turcotte et al., 2011;Yoshida et al., 2003). However, complimentary studies of evolved responses in natural populations are scarce (Walsh et al., 2012). Unfortunately, evolutionary response to predation can rarely be addressed by observing such phenomena in the wild, as ecological and evolutionary effects of the predator on population dynamics occur simultaneously and are confounded. For example, this has been considered to be a major constraint in studies of evolutionary responses of fish populations to harvesting (Heino et al., 2015). Ellner et al. (2011) provided an approach to disentangle the contribution from ecological and evolutionary effects in shaping changes in phenotypes of natural populations, but only for populations where pedigree information is available. An alternative approach is to compare the population dynamics of genotypes that have evolved under contrasting predation regimes in a common controlled environment. Bassar et al. (2013) employed such an approach to model population growth rate in guppies (Poecilia reticulata) based on demographic traits, and to compare modeled population growth across densities for populations originating from locations with different predation pressures.
Although this approach has limitations in terms of extrapolating results to the wild (e.g. due to genotype-by-environment interactions), it can at the least demonstrate that a change in predation regimes of wild populations can translate into evolutionary effects on the in- Here, we leverage the invasion of the North European spiny water flea (Bythotrephes longimanus), a predator of herbivorous zooplankton, into Lake Kegonsa, Wisconsin, US (Walsh et al., 2016). Since Bythotrephes was first detected at high abundance in this lake in 2009, the biomass of one of its prey species, the cladoceran Daphnia pulicaria, has been reduced by up to 60% (Walsh et al., 2016). In a previous study of D. pulicaria from this lake, Landy et al. (2020) compared resurrected clones (i.e. genotypes) originating from prior to the invasion of Bythotrephes with those of contemporary (i.e., postinvasion) clones, and provided evidence that the invasion has led to evolutionary change in a suite of life history and behavioral traits.
Specifically, they demonstrated that invasion by Bythotrephes was associated with evolved reductions in size at maturity and fecundity.
In the current study, we follow-up these findings and conduct population growth experiments in a common controlled environment to determine if the proliferation of Bythotrephes has led to evolutionary shifts in D. pulicaria population dynamics characteristics.   (Kilham et al., 1998) and fed non-limiting supply of green algae (Scenedesmus obliquus, ~1.0 mg C L −1 day −1 ). In December 2019, live individuals of each clone were transported to the Norwegian University of Science and Technology and were subsequently kept at 17°C (photoperiod 16 L:8D) in ADaM medium (Klüttgen et al., 1994, SeO 2 concentration reduced by 50%, sea-salt increased to 1.23 g/L) and fed non-limiting supply of Shellfish Diet 1800 (Reed Mariculture Inc.) until the onset of the experiment. The same medium, food, temperature and light regime was used throughout the rest of the study.

| Study animals and husbandry
For each clone, 5-10 adult individuals were randomly chosen from stock cultures and placed in separate 1 L glass beakers (1 clone per beaker) where they were fed Shellfish Diet 1800 (Reed Mariculture Inc.) three times per week at a concentration of 4 ×x 10 5 cells/ml.
When several egg-bearing individuals were identified in each beaker, these were isolated by removing all others. Beakers with egg-bearing individuals were checked for neonates every 24 h and each neonate found within this period was individually transferred to a plastic container containing 100 ml of ADaM. Newborn individuals that died within 6 days were replaced using the same method. In total, 10 newborns were selected from each clone over a span of 7 days, yielding 190 populations that originally consisted of a single individual (10 individuals per clone × 19 clones, 8 pre-invasion, 11 post-invasion). When an individual died after the 7th day, it was reported as dead and not replaced (n = 11). In total, 21 pre-invasion populations (23.3%) and 9 post-invasion (9.0%) populations went extinct during the experiment, and 66.7% of these extinctions occurred within the first 10 days of the experiment. The plastic containers were stored in 2 Memmert Peltier cooled incubator IPP 260plus (Memmert, Germany) climate cabinets at 20.0°C (photoperiod 16 L:8D). Shellfish Diet 1800 (4 × 10 5 cells/ml) was added every second day and medium was changed every 8 days.
Container placement in the climate cabinets was changed haphazardly every 2 days, after feeding. All populations time series were run in parallel during March-May 2020.

| Measuring population growth
To obtain data on population growth, video recordings were made of each population every 8 days (with one exception due to covid-19 regulations) for a period of 3 months, starting 11 days after the start of the population growth experiment, creating a total of 10 censuses.
From these videos, we estimated the number of individuals and the total dry biomass for each population in each census using the R package trackdem (Bruijning et al., 2018). Daily population growth rates, G, were calculated as G = log e (N t + 1 /N t )/d, where N t and N t + 1 is the population abundance (measured either in number of individuals or total biomass) at two consecutive censuses, and d is the duration of time in days between the censuses. For further details on these procedures, see Appendix 1.

| Choice of population dynamics model and calculation of r
Inspection of population growth rate data (both for numerical and biomass growth) revealed strong non-linearity in the density dependence. Thus, the population dynamics is best described by the theta-logistic model G = r(1−[N t /K] θ ) where r is the intrinsic population growth rate, K is the carrying capacity, and θ determines the shape of the density dependence. One concern when fitting this model to data is that different combinations of r and θ can produce model fits of similar likelihood, potentially resulting in ecologically unrealistic estimates of both parameters (Clark et al., 2010). We therefore took advantage of the experimental design, whereby each population was started with a single neonate, which allowed us to obtain direct observations of population growth under low density. We based this calculation of r for each population on the observed population abundance at the start of the experiment (N 0 ) and the second census (N 2 , i.e. r = log e [N 2 /N 0 ]/d). This duration (i.e.

| Statistical analyses
All statistical analyses were conducted in R v.4.1.0 (R Core Team, 2021). We tested for an effect of invasion history on the observed values of r (calculated early in the sequence of population growth and thus under low density, see above) by fitting a linear mixed model to these data using the function lme in the package nlme (Pinheiro et al., 2021), including a random effect of clone ID, and comparing this against a simpler model containing only the random effect using the Akaike information criterion corrected for small sample sizes (AICcmodavg, Mazerolle, 2020). This was done both for numerical and biomass growth rate.
Next, to test for an effect of invasion history on K and θ, we fitted non-linear mixed effect models representing the theta-logistic model to the population growth time series. This was done for the data following the first two censuses, i.e. after the period that had been used to calculate r. Again, this was done separately for numerical and biomass growth rates. In the full model, population growth rate over a given period (between two consecutive censuses) was modeled as a function of the observed value of r for that population and its population size at the start of that period while estimating values of K and θ that depended on invasion history. The model included random effects of clone ID and population nested within clone ID on K. The model was fitted using the function nlme in the package nlme (Pinheiro et al., 2021). This full model was compared with a simpler one where K and/or θ were common to all populations independent of invasion history, again using the Akaike information criterion corrected for small sample sizes. For all analyses, assumptions of normality and homogeneity of residuals were satisfied. All figures were made using the package ggplot2 (Wickham, 2016).

| RE SULTS
The population dynamics of most clones consisted of an initial increase in abundance up to a peak value, followed by a subsequent decline ( Figure A1). According to the model comparison, observed values of intrinsic population growth rate tended to depend on invasion history ( Table 1). The evidence for such an effect in terms of ΔAIC was strongest for biomass (Table 1), but models containing an effect of invasion history suggested higher intrinsic population growth rate in post-invasive clones than in pre-invasive clones for both numerical and biomass growth ( Table 2). The estimated increase in intrinsic population growth rate based on numerical and biomass data were 23% and 15%, respectively ( Table 2).
The difference in population dynamics between pre-and postinvasion clones was also reflected in the comparisons of theta-logistic models of population growth rates. Particularly for numerical data, the model comparison shows strong evidence for an effect of invasion history on K ( Table 1). Evidence for effects of invasion history on θ was weaker for these data, and the best model did not contain such a term ( Table 1). For biomass data, the best model contained effects of invasion history on both K and θ ( Table 1). Models containing such effects on only one of these two parameters received slightly less support, whereas a model containing no effect of invasion history received considerably less support (Table 1). For both these analyses, the best models predicted a higher carrying capacity for post-invasion clones compared with for the pre-invasion clones (Table 3, Figure 1). The estimated increase in carrying capacity based on numerical and biomass data were 27% and 23%, respectively (Table 3).

| DISCUSS ION
We leveraged a well-documented invasive event by a predatory species (Bythotrephes longimanus) in combination with a resurrection ecology approach to investigate evolutionary changes in population dynamics of their main prey species (Daphnia pulicaria). This allowed us test for an association between a change in predation levels and evolution of prey population dynamics in a natural population.
Estimation of population dynamics parameters showed that postinvasion genotypes had an increased intrinsic population growth rate r (strong evidence for biomass data) as well as carrying capacity K (strong evidence for numerical data). For the current study, we had prior information on evolutionary change in individual traits in response to the predator invasion. Specifically, for traits that can be directly used to predict effects on population dynamics, Landy et al. (2020) found that post-invasion clones had a reduced size at maturity and reduced fecundity compared with pre-invasion clones, but found no significant difference in age at maturity. Based on this, one might predict a reduced intrinsic rate of increase in post-invasive clones. Yet, we observed the opposite in our study. One reason for this could be that unmeasured traits that may be genetically correlated with those that were measured may also have influenced  (Shertzer et al., 2002;Yoshida et al., 2003). This system has also been used to test for an evolutionary trade-off between algal population growth and predator defense (Kasada et al., 2014;Meyer et al., 2006;Yoshida et al., 2004). Yoshida et al. (2004) allowed algae to evolve in the presence and absence of the rotifer, after which their growth TA B L E 1 AICc comparisons of candidate models explaining variation in observed intrinsic population growth rate (r) carrying capacity (K) and the shape of density dependence (θ) of experimental populations of Daphnia pulicaria originating from Lake Kegonsa. For r, separate linear mixed effects models are fitted to numerical and biomass population growth rate as dependent variables. Full models include effects of invasion history (whether the population originates from a period before or after invasion by the predatory zooplankton species Bythotrephes longimanus), with clone ID as a random effect. For K and θ, separate theta-logistic models are fitted to numerical population growth rate and biomass population growth rate as dependent variables. Full models include effects of invasion history on K and θ, with clone ID and population (nested within clone ID) included as random effects. having evolved in absence of predation, but only at the most limiting nutrient-level. When nutrients were more abundant no such difference was observed. Thus, in this system it appeared that adaptation to predation primarily caused a reduction in K, with no effect on r.
In the current study, we found no indication of an interaction between predation history and food availability in determining the rate of population growth as observed in the rotifer/algae experiments. If such an interaction had been present, we would expect the relative difference in population growth rate between the two types to differ at low and high population density, which in turn should translate into different effects of type on r and K. Instead, models that contained an effect of invasion history consistently suggested elevated values of both r and K in post-invasion clones compared with pre-invasion clones, independent of measurement type (numerical or biomass population dynamics), meaning that postinvasion clones showed weaker effects of numerical density (and per capita food availability) on per capita growth rate when compared with pre-invasion clones, with slightly higher overall growth rates.
This may indicate that evolutionary responses to invasive predators in natural systems can be complex. One caveat of the result in this study is that they are obtained in a laboratory environment, with different biotic (e.g. competition, parasites, and food quality) and abiotic (e.g. temperature, water chemistry) conditions than what this population experience in the wild. Specifically, in the presence of genotype-by-environment interactions, a given phenotypic difference between two genotypes observed in one environment may disappear or even be reversed in a different environment. This is a general concern with such common environment experiments that aim to describe genetically based differences in phenotype. A second caveat is that whereas all pre-invasion clones were derived from ephippia, and thus represent a sub-set of clones that would hatch at to test for an effect of collection methodology, inspection of clonespecific population dynamics does not indicate a consistent difference in population dynamics between these two groups ( Figure A1).

TA B L E 2
Parameter estimates (obtained using REML) of the best fitting models describing the variation in observed intrinsic population growth rate (r) of Daphnia pulicaria originating from Lake Kegonsa ( Table 1). Post-and pre-invasion populations consist of clones originating from after and before Bythotrephes longimanus invasion, respectively.

TA B L E 3
Parameter estimates (obtained using REML) for K, θ and random effects of K of the best fitting theta-logistic models describing the population dynamics of Daphnia pulicaria originating from Lake Kegonsa ( Table 1). Post-and pre-invasion populations consist of clones originating from after and before Bythotrephes longimanus invasion, respectively. this, such that when reared in a common environment, those from more food-restricted environments actually grow faster (Arendt & Wilson, 1999). Thus, our observation that post-invasive clones outperform pre-invasive clones may be an example of such countergradient variation (Conover & Schultz, 1995), where they have evolved physiological adaptations that increase population growth under a given level of food abundance. In a recent study of D. pulicaria from Lake Mendota, Rani et al. (2022) found that post-invasive clones had a reduced metabolic rate compared with pre-invasion clones, which may represent one such physiological adaptation to a cooler resource deficient environment. Furthermore, Einum et al. (2019) found that variation in somatic growth rate among clones of D. magna was best explained by clone-specific food consumption expressed relative to their rate of energy use. Thus, if post-invasion clones have a reduced rate of metabolism but do not moderate food consumption when reared in a common environment (as in the current study), this could be expected to result in increased somatic growth rate and may thus explain the higher population growth rate.
To conclude, the current study demonstrates an evolutionary shift in the population dynamics of D. pulicaria in parallel with an increase in predation brought about by invasion of the predatory zooplankton Bythotrephes. To our knowledge, this is the first empirical study that directly demonstrates this by comparing genotypes of a single natural population that has experienced a temporal change in exposure to predation. Although the potential F I G U R E 1 Growth rate in terms of (a) number of individuals and (b) total dry mass for pre-and post-invasion populations of Daphnia pulicaria originating from Lake Kegonsa. Pre-and post-invasion populations consist of clones sourced from before and after Bythotrephes longimanus invasion in 2009, respectively. Regression lines give predictions from theta-logistic models with parameter estimates from Table 2 (for r) and Table 3 (for K and θ). for genotype-by-environment interactions prevents strong conclusions regarding fitness of pre-and post-invasion clones in the wild, we suggest that the observed shift in population dynamics may be related to selection for reduced predator exposure, either temporally or spatially. If so, this suggests that complexities in ecological interactions represents a challenge when predicting the evolutionary responses of population dynamics to changes in predation pressure in natural systems. formal analysis (supporting); funding acquisition (equal); methodology (equal); project administration (supporting); resources (equal); software (lead); supervision (equal); writing -review and editing (supporting).

ACK N OWLED G M ENTS
The study was funded through the Research Council of Norway's (RCN) project numbers 230482, 244046 and 223257. We thank V Rani for lab assistance.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

O PEN R E S E A RCH BA D G E
All data used in this study are provided as private-for-peer review at https://figsh are.com/s/7e33e ed4ec f3c36 e56fe. These data will be made publicly and permanently available at Figshare upon acceptance for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used in this study are provided as private-for-peer review at https://figsh are.com/s/7e33e ed4ec f3c36 e56fe. These data will be made publicly and permanently available upon acceptance for publication.