Giacomo Zilio,Sascha Krenek,Claire Gougat-Barbera,Emanuel A Fronhofer,Oliver Kaltz
Predicting evolution in experimental range expansions of an aquatic model system
dispersal evolution,reaction-diffusion model,predictability,eco-evolutionary dynamics,invasions
Abstract Predicting range expansion dynamics is an important goal of both fundamental and applied research in conservation and global change biology. However, this is challenging if ecological and evolutionary processes occur on the same time scale. Using the freshwater ciliate Paramecium caudatum, we combined experimental evolution and mathematical modeling to assess the predictability of evolutionary change during range expansions. In the experiment, we followed ecological dynamics and trait evolution in independently replicated microcosm populations in range core and front treatments, where episodes of natural dispersal alternated with periods of population growth. These eco-evolutionary conditions were recreated in a predictive mathematical model, parametrized with dispersal and growth data of the 20 founder strains in the experiment. We found that short-term evolution was driven by selection for increased dispersal in the front treatment and general selection for higher growth rates in all treatments. There was a good quantitative match between predicted and observed trait changes. Phenotypic divergence was further mirrored by genetic divergence between range core and front treatments. In each treatment, we found the repeated fixation of the same cytochrome c oxidase I (COI) marker genotype, carried by strains that also were the most likely winners in our model. Long-term evolution in the experimental range front lines resulted in the emergence of a dispersal syndrome, namely a competition—colonization trade-off. Altogether, both model and experiment highlight the potential importance of dispersal evolution as a driver of range expansions. Thus, evolution at range fronts may follow predictable trajectories, at least for simple scenarios, and predicting these dynamics may be possible from knowledge of few key parameters.
Predicting evolution in experimental range expansions of an aquatic model system

Predicting range expansion dynamics is an important goal of both fundamental and applied research in conservation and global change biology. However, this is challenging if ecological and evolutionary processes occur on the same time scale. Using the freshwater ciliate Paramecium caudatum, we combined experimental evolution and mathematical modeling to assess the predictability of evolutionary change during range expansions. In the experiment, we followed ecological dynamics and trait evolution in independently replicated microcosm populations in range core and front treatments, where episodes of natural dispersal alternated with periods of population growth. These eco-evolutionary conditions were recreated in a predictive mathematical model, parametrized with dispersal and growth data of the 20 founder strains in the experiment. We found that short-term evolution was driven by selection for increased dispersal in the front treatment and general selection for higher growth rates in all treatments. There was a good quantitative match between predicted and observed trait changes. Phenotypic divergence was further mirrored by genetic divergence between range core and front treatments. In each treatment, we found the repeated fixation of the same cytochrome c oxidase I (COI) marker genotype, carried by strains that also were the most likely winners in our model. Long-term evolution in the experimental range front lines resulted in the emergence of a dispersal syndrome, namely a competition—colonization trade-off. Altogether, both model and experiment highlight the potential importance of dispersal evolution as a driver of range expansions. Thus, evolution at range fronts may follow predictable trajectories, at least for simple scenarios, and predicting these dynamics may be possible from knowledge of few key parameters.

Predicting ecological dynamics and species’ range shifts has become a major goal for conservation and management strategies in times of global climate and environmental change (Petchey et al., 2015). Indeed, whether the outcomes of range expansions or biological invasions can be predicted at all remains highly debated in ecology even in simple settings, due to the intrinsic stochasticity of these phenomena (Giometto et al., 2014; Melbourne & Hastings, 2009). Moreover, evolutionary processes occur at the same time scale as ecological dynamics during range expansions (Perkins et al., 2013; Williams et al., 2016), potentially exacerbating the uncertainty of outcomes (Williams et al., 2019). Thus, assessing the predictability of evolutionary trait change during range expansion is fundamental to our general understanding of the phenomenon, but can also guide the management of such eco-evolutionary dynamics. Theory shows that range expansions can involve the concurrent evolution of dispersal and other traits (Kubisch et al., 2014; Perkins et al., 2013) and lead to the emergence of dispersal syndromes (Clobert et al., 2012; Cote et al., 2017). Individuals with greater dispersal propensity are the first to reach the range front, and they will reproduce with conspecifics that have the same fast spreader characteristics (Hughes et al., 2007; Thomas et al., 2001). Consequently, high dispersal ability and correlated life-history traits evolve in the range front populations due to spatial selection and spatially assortative mating (Phillips et al., 2008; Shine et al., 2011). Since expansion speeds are mainly influenced by dispersal and reproduction (Fisher, 1937; Kolmogorov et al., 1937), the two traits can be rapidly selected and evolve simultaneously. However, if dispersal is costly (Bonte et al., 2012) there may be trade-offs with other traits. Higher reproduction at the range front may come at the expense of lower competitive ability (Burton et al., 2010), recalling the competition–colonization trade-off in classic species coexistence models (Calcagno et al., 2006). Fast evolution in range front populations can produce eco-evolutionary feedbacks that may speed up the expansion process (Chuang & Peterson, 2016; Miller et al., 2020; Ochocki et al., 2019; Shine et al., 2011; Williams et al., 2019). In the emblematic example of the cane toad (Rhinella marina) expansion in Australia, increased dispersal at the range front coincided with evolutionary change in behavioral, morphological, and demographic traits, accelerating the speed of the toad expansion (Perkins et al., 2013; Phillips et al., 2006). Growing empirical evidence from other natural populations and biological systems (Lombaert et al., 2014; Simmons & Thomas, 2004) suggest that dispersal evolution at range fronts is a common phenomenon. Recently, experimental evolution and microcosm landscapes have been used to test fundamental predictions and mimic range expansions in the laboratory. Experiments with ciliates (Fronhofer & Altermatt, 2015; Zilio et al., 2023), arthropods (Ochocki & Miller, 2017; Petegem et al., 2018; Szűcs et al., 2017; Weiss-Lehman et al., 2017), or plants (Williams et al., 2016) showed the rapid evolution of dispersal and other dispersal-related traits during the experimental range expansions, such as growth, size, locomotory capacity, or exploratory behavior. However, whether we can accurately predict these processes and the accompanying trait evolution from prior information on the genetic or phenotypic characteristics of the expanding populations remains an open question (Angert et al., 2020). Coupling microcosm experiments with mathematical modeling and genetic analyses provides a possible way forward to assess the predictability of range expansions, when evolution is at play (Nosil et al., 2020). In micro/mesocosm landscapes, we can study the repeatability of outcomes through independent replicates under controlled conditions. Using specifically tailored and parameterized mathematical models, we can formalize putative processes of range expansion dynamics and confront predicted with observed results. Genetic analysis can further characterize the degree of similarity among experimental replicates and link phenotypic trait change to genetic change. Here, we employed such a combined approach to assess the predictability of evolutionary outcomes of range expansions in an aquatic model organism, the freshwater protozoan Paramecium caudatum. Following previous studies (Fronhofer & Altermatt, 2015; Nørgaard et al., 2021), we used interconnected two-patch systems to mimic populations at an expanding range front, where recurrent episodes of natural dispersal alternated with periods of population growth (range front treatment; Figure 1A). In the contrasting range core treatment, only the non-dispersing individuals were maintained, while in a third control treatment the dispersing individuals were always put back together with the non-dispersing ones (Figure 1A). We recreated these experimental treatments in a predictive mathematical model, parameterized for dispersal and growth characteristics of the 20 Paramecium strains that were used to assemble the founder populations (total of 15 lines) in the evolutionary experiment. Thus, starting from standing genetic variation, we compared predicted and observed short-term trait evolution in front and core populations, and assessed the repeatability of evolutionary outcomes at the genotypic level after 30 dispersal/growth cycles. Long-term evolution was tracked over the course of three years (160 cycles). Our main finding is that short-term evolutionary outcomes were highly predictable and essentially depend on two parameters (genetic variation in dispersal and growth). In the long run, dispersal syndrome emerged at range fronts, possibly due to de novo evolution.

Paramecium caudatum is a freshwater ciliate with a world-wide distribution, feeding on bacteria and detritus. Asexual reproduction occurs by mitotic division and represents the main mode of population growth. Swimming is accomplished through the coordinated movement of ciliary bands on the cell surface (Wichterman, 1986). Previous work on P. caudatum indicated a genetic basis of dispersal propensity (Zilio et al., 2021). Here, we used 20 P. caudatum strains (i.e., clonal cultures derived from a single individual) from various geographic origins (Weiler et al., 2020; Zilio et al., 2021) and representing different groups of cytochrome c oxidase I (COI) genotypes (Supplementary Table S1). COI genotype here refers to haplotypes of the mitochondrial COI gene, which can be assigned to different haplogroups or clades (i.e., all founder strains belong to clade A or B; Johri et al., 2017). All cultures were reared under standard laboratory conditions in lettuce medium with the food bacterium Serratia marcescens at 23°C, allowing up to three asexual doublings per day (Nidelet & Kaltz, 2007).

Prior to the start of the long-term experiment, we assayed the 20 founder strains for dispersal and population growth characteristics (Supplementary Table S1). For the dispersal assay, we placed aliquots of 8 mL of culture (at equilibrium density, 145 individuals/mL on average) in a two-patch system (additional details below) and let the Paramecium disperse for 3 hr. Once connections were blocked, we estimated the number of residents and dispersers by taking 150–600 µL samples from the two tubes and counting the number of individuals under a dissecting microscope. Dispersal was taken as the proportion of dispersers of the total number of individuals in the system. We tested four replicates per strain, over two experimental blocks. Dispersal rates were then estimated from a generalized linear mixed model (frequentist approach with binomial error distribution; glmer in “lme4” R package with bobyqa optimizer (Bates et al., 2015)) with strain as fixed factor and experimental block and observation (to account for overdispersion) as random effects. The posteriors obtained with this approach were used to parametrize the mathematical model for dispersal (see below). For the growth assay, we placed ca. 200 individuals (from cultures at equilibrium density) in 20 mL of fresh medium. Over the course of 6 days, we tracked population density by counting the number of individuals in daily samples of 100–200 µL. We tested three replicates per strain. Using a Bayesian approach (Rosenbaum et al., 2019), we estimated the intrinsic population growth rate (r0) and equilibrium density () for each replicate by fitting a r-K population growth model to the time series data. The posterior distributions of r0 and were used to parametrize the mathematical model as for dispersal, details of the Bayesian fitting are given in Supplementary Information (Supplementary Appendix).

The evolutionary experiment comprised a sequence of cycles, where dispersal events alternated with periods of population growth. The founder population was created by mixing the 20 strains at equal proportions in a single culture, which was then divided up into 15 replicate lines, assigned to the following three treatments. First, in the range front treatment (six lines), we placed the Paramecium in one of the two tubes in two-patch dispersal systems (interconnected 15-mL tubes, Figure 1A). Connections were opened for 3 hr, during which time individuals were allowed to swim to the other tube. We then collected the dispersers and cultured them for 1 week under permissive conditions in 20 mL of fresh medium (in 50-mL plastic tubes), until we initiated a new round of dispersal, again only retaining the dispersers and culturing them for 1 week, and so on. Second, the range core treatment (six lines) followed the same cycles of dispersal and growth, but only the non-dispersing residents were retained after each dispersal episode. Third, in the control treatment (three lines), residents and dispersers were mixed after each dispersal event and then cultured for 1 week, as in the other treatments. In corollary, the range front treatment mimics the advancing cohort of a spatially expanding population, whereas populations from the core treatment remain in place and constantly lose emigrants. The control treatment is similar to the core treatment, except for the loss of emigrants. A total of 161 cycles were accomplished. Prior to each dispersal event, ca. 1,800 individuals (median; 25%/75% quantile range: 1,400/2,700) were placed in the dispersal systems. After dispersal, the number of individuals starting the 1-week growth period were matched between treatments. Because dispersal rates were low at the beginning, these starting numbers were initially set to 200 individuals (placed in a total volume of 20 mL of fresh medium). During the following 1-week growth period, stable population sizes were typically reached within 3–4 days, with densities of ca. 240 individuals per ml (median; 25%/75% quantiles: 180/360). After cycle 32, when dispersal had already reached higher levels (see Results), we adjusted the starting numbers to ca. 1,500 (median; 25%/75% quantiles: 1,100/2,000).

For each line, dispersal was measured at each dispersal event and equilibrium densities (), taken at the end of the 1-week growth period at each cycle. Furthermore, growth rate (r0) was determined in assays conducted at cycle 21 (year 1), 78 (year 2), and 160 (year 3), as described above, with two to three replicates per line and year. Bayesian model fitting was used to estimate r0 (Supplementary Figure S1).

All founder strains were genotyped for the mitochondrial COI gene, based on DNA extracted from 10 cells per strain, using the Chelex100 (Bio-Rad Laboratories GmbH, Germany) method (Barth et al., 2006). Sequence data were compared with the non-redundant sequence database using NCBI-BLAST and an in-house database to infer COI genotype affiliation. All sequence data are deposited in GenBank (Barth et al., 2006; Weiler et al., 2020). At cycle 30 in the evolutionary experiment, we tested for the absence or presence of founder strains in the evolved lines. To this end, DNA from mixes of 50 cells from each line was extracted and analyzed for (multiple) COI marker signals, using a restriction fragment length polymorphism method, established previously (Killeen et al., 2017). This method characterizes the line for the most frequent COI genotype, and it has a resolution threshold of c. 5%, i.e., it can detect two to three cells of a minority genotype in the sample of 50 cells (Killeen et al., 2017). Details of the protocols are provided in Supplementary S9. Note that, although being a convenient and widely used marker for intraspecific genetic variation in Paramecium (e.g. Greczek-Stachura et al., 2021; Przyboś et al., 2019; Tarcz et al., 2012) the COI gene is only moderately variable. In our study, 12 strains shared the same COI genotype (Supplementary Table S1) and could therefore not be distinguished in the founder mix. Data for more informative markers (such as SNP genotypes) were not available.

Our model is designed to capture the specificities of the evolutionary experiment and the characteristics of the strains in the founder population. Thus, we model the population dynamics of Paramecium strains, assuming logistic growth following the Verhulst equation (Henle et al., 2004) expanded to include both intra- and inter-strain competition: where is the population size of strain i, is its intrinsic rate of increase and as the competition coefficients. The model is parametrized with the posteriors extracted from growth curve fitting (r0, ), as described above. We make the simplifying assumption that intra- and interspecific competition is of equal strength. We further assume a quasi-extinction threshold of 0.7 (we have tested the effect of different quasi-extinction thresholds ranging from 0.0001 to 0.9), which implies that strains experience an extinction if they exhibit densities below this value. We model the community dynamics of the strains for 7 days, followed by a 3-hr dispersal phase in a two-patch metapopulation. During this 3 hr dispersal phase, all strains can disperse from their patch of origin to their destination patch according to the dispersal rates estimated from dispersal assay; the model is parametrized with posteriors extracted from the frequentist statistical analysis described above. After the dispersal phase, we follow the patch of origin (residents in the range core treatment), the destination patch (dispersers in the front treatment) or the combined patches (dispersers and residents mixed in the control treatment). We repeat this procedure for a total of 10 iterations. As in the experiment, we control for densities between rounds of iteration by selecting the equivalent of 10 mL samples. This approach allows us to predict, based only on measurements of growth parameters and dispersal rates of the founder strains, which strains should predominate in each of the three treatments at the end of the experiment. It is important to keep in mind that the underlying model is deterministic. However, since we parametrize the model with draws from posteriors, our approach takes into account the uncertainty associated with the data and yields a distribution of likely outcomes, given these uncertainties. Note that our model depicts a scenario of selection from standing genetic variation; it does not include mutational change.

All statistical analyses were performed in R (ver. 4.2.0) and JMP 14 (SAS Institute, Inc., 2018) using a frequentist approach except for the pairwise trait correlations, where a Bayesian approach was preferred due to its higher flexibility. We analyzed dispersal (proportion of dispersers), using generalized linear models (GLM) with binomial error distribution. We considered evolutionary treatment (core, front, control), experimental cycle and line (nested within treatment) as explanatory variables. We analyzed variation in intrinsic population growth rate (r0) and equilibrium density (; averages per line and year) using GLMs, with evolutionary treatment, year and line as explanatory variables. To illustrate how selection acts on standing genetic variation in our model, we associated the winning probability of each of the 20 founder strains (i.e., the fixation probability among the 20 strains in 10,000 model runs) with their respective median values of dispersal, r0, and from the distributions used by the model. For each treatment, we then performed multiple regressions, with winning probability as response variable and the three traits as explanatory variables. To investigate associations between dispersal, r0 and , we constructed a data matrix based on trait means per year and line (3 years × 15 lines, n = 45), after centering and scaling trait distributions. One range-front line was lost in year 3, leading to n = 44. We also performed a principal component analysis (PCA), considering the joint variation in all three traits. Finally, we used a Bayesian approach (“rstan” package version 2.21.2 (Carpenter et al., 2017)) to estimate pairwise trait correlations (Supplementary Figure S3).

Over the first 25 cycles of the evolutionary experiment, we observed a strong increase in dispersal in the range front treatment (Figure 1B). Dispersal reached 22.3% (± 0.012 standard error [SE], averaged over cycles 15–25) at the front, compared to only 4.4% (± 0.004 SE) in the core treatment and 6% (± 0.012 SE) in the control treatment (effect of selection treatment: χ22 = 119.7; p < .001). Increased dispersal at the front established within only a few cycles, and was formally significant for the first time at cycle 8 (cycle-by-cycle analysis: p < .001). Our parametrized model captured this rapid increase of dispersal in the range front treatment (Figure 1B), and there was a quantitative match between the distribution of endpoint levels of dispersal in the model and observed values in the experiments (Figure 2A). The model further predicted general increases in growth rate (r0) and equilibrium density () in all treatments. For r0, the values ranged from 0.07 in the ancestral mix to 0.08 in core and front end-point populations. This is consistent with results from the growth assay conducted at cycle 21, where estimates of r0 for the 15 lines are well within the central range of predicted values in the model (Figure 2B). As predicted, the evolutionary treatments did not significantly differ in r0 (treatment: F2,12 = 1.2; p = .354). Unlike in the model (Figure 2C), range front lines produced nearly 20% higher equilibrium densities than did range core and control lines (treatment: F2,12 = 11.1; p = .003).

Our model finds strong variation in the fixation probability among the 20 strains, and different treatments have different most likely winners (range: 0.7%–16.8%; Supplementary Figure S7). For the range front treatment, multiple regression analysis (Supplementary Table S2) shows that both dispersal and r0 are positively associated with strain winning probability, and this with equal strength (standardized beta [β] regression coefficients: +0.55 and +0.64, respectively; Figure 3). Thus, selection is predicted to favor strains that both disperse more and grow faster. In contrast, in the core and control treatments, strain winning probability is mainly associated with high growth rate (β > +0.96), accompanied by weak selection against clones with higher dispersal (β ≤ −0.27) or equilibrium density (β ≤ −0.33). Molecular analysis of the 15 lines based on COI genotype indicates complete genetic divergence between selection treatments. For all nine range core and control lines, only the b05 COI genotype was detected. The two strains in the founder population that carry this genotype (Supplementary Table S1) have very high growth rates and very low dispersal, the trait combination favored in the model. Indeed, the candidate strain AMF_11_1A has the highest growth rate overall and is the most likely winner in core and control treatments according to our model (Figure 3). In contrast, all six range front lines appear to be fixed for the b07 COI genotype. This genotype is shared by 12 founder strains (Supplementary Table S1), which may thus have gone to fixation in groups or individually. Among these candidate strains is the most likely winner (goe_14) predicted by the model: it has the highest growth rate and the third-highest dispersal, in line with the prediction of the two traits being under joint positive selection in this treatment. As shown in Figure 2, trait values of the most likely front and core winner strains (goe_14 vs. AMF_11_1A; strain posterior distributions on the right) show a good match with both the predicted model outcomes (left distributions; Figure 2) and the experimental data.

In addition to the short-term evolution, we also observed a long-term increase in dispersal in the range front treatment over the entire time span of the 3 years of the experiment (cycle × treatment interaction: χ22 = 88.8; p < .001; Figure 1B). This trend is significant, even when omitting the first 50 cycles (χ22 = 51.7; p < .001). We found little evidence for a dispersal difference between range core and control lines, neither overall (contrast core vs. control: p > .68) nor when considering individual cycles (11 cycle-by-cycle contrasts with .0078 < p < .09, none significant after correction for multiple testing). While no significant treatment effects were detected in the first growth assay (cycle 21, see above), range front lines had nearly twofold lower values of r0 than range core lines in assays conducted in year 2 and 3 (year × treatment: F4 = 6.66; p < .001; Supplementary Figure S1A). Furthermore, while beginning to grow more slowly, range front lines continued to produce up to twofold higher than range core and control lines (treatment: F2 = 34.21; p < .001; Supplementary Figure S1B). Figure 4A–C illustrates short- and long-term trends in pairwise trait associations, in relation to the model predictions. For dispersal and r0 (Figure 4A), there was no clear relationship between the two traits after short-term selection (year 1). However, in year 2 and 3, observed data points tend to fall outside the main predicted ranges, and a negative relationship between dispersal and r0 emerged (Figure 4A). This negative association is highly significant over all lines and years combined (r = −0.627, 95% credible interval [CI] [−0.771; −0.434]), but also holds for year 2 and 3 separately (Supplementary Figure S3). The positive relationship between dispersal and , already observed as a short-term trend, further consolidated in year 2 and 3 (Figure 4B), again with values mostly falling outside the main predicted short-term ranges. The correlation is significant overall (r = 0.599, 95% CI [0.347; 0.725]), as well as for each year separately (Supplementary Figure S3). Furthermore, diverging trends in core and front lines lead to a negative association between r0 and (Figure 4C). The negative correlation is of intermediate effect size overall (r = −0.325, 95% CI [−0.575; −0.031]), and is significant in all three years separately (Supplementary Figure S3). It should be noted that all of these main trends of divergence hold, when we correct for year effects, by expressing front and core line data relative to the control treatment in each year (Supplementary Figure S5). PCA (Figure 4D) summarizes the patterns of phenotypic divergence. Demography-related traits and dispersal are pulling in approximately equal strength on PC axis 1, but in opposite directions (PC1 loadings: r0 = −0.53; = +0.57; dispersal = +0.62). Thus, range front lines are characterized by a combination of higher equilibrium density and dispersal, but lower intrinsic population growth rate relative to range-core and control lines (MANOVA: F2,37 = 10.85, p < .001). The separation of clouds indicates the progressive divergence through time, with a maximum in year 3. There is little differentiation between range core and control treatments.

Predicting range expansions with ecology and evolution occurring on the same timescale is a challenging task. While previous ecological range expansion studies (Giometto et al., 2014; Melbourne & Hastings, 2009) investigated the predictability of range expansion speed, we focused on the predictability of evolutionary change during replicated range expansions. We included short-term evolution from standing genetic variation in a simple model parameterized for our laboratory system and confronted predicted evolutionary outcomes with results from experimental range expansions. Both model and experiment show rapid divergence between range core and front treatments, with selection for higher dispersal at the front. Over longer time scales, experimental range core and front populations continued to diverge, indicating de novo evolution and resulting in the emergence of dispersal syndromes.

In the context of reaction-diffusion models, dispersal (diffusion) and population growth at low densities are the two key traits for understanding and predicting range expansion dynamics (Fisher, 1937; Kolmogorov et al., 1937). Consistent with this view and previous studies (Phillips et al., 2010; Shine et al., 2011), dispersal and population growth rate were here identified as main targets of selection. Higher dispersal was immediately selected from standing genetic variation at the range front and weakly selected against in the range core in the model as well as in the experiment, where range front populations showed increased dispersal already after the first few cycles. Such strong and fast selection on dispersal in the vanguard front populations has been found in similar experiments (Fronhofer & Altermatt, 2015; Ochocki & Miller, 2017; Petegem et al., 2018; Szűcs et al., 2017; Weiss-Lehman et al., 2017; Williams et al., 2016), but also in natural populations (Perkins et al., 2013; Phillips et al., 2006). Dispersal evolution might therefore accelerate the speed of range expansion already over very short time scales (Miller et al., 2020). Contrary to more standard views of range expansion with r- and K-selection (Burton et al., 2010; Miller et al., 2020), growth rate was under positive short-term selection in both range core and front treatments. This can be explained by the fact that populations in all treatments experienced regular bottlenecks, thus imposing general selection for increased growth rate, a trait for which there was ample variation among founder strains (Supplementary Figure S4). Importantly, however, our model shows that dispersal and growth rate can be simultaneously selected in the range front treatment (Figure 3). Whether one or the other trait has more weight depends on the stochasticity introduced by the quasi-extinction threshold. With small values, even weak dispersers make it into the new patch and can subsequently regrow to high density without experiencing an extinction. Indeed, additional model scenarios show that, when we decrease the quasi-extinction threshold, selection for growth rate overrides selection for dispersal and the strain with the highest growth rate becomes fixed in all treatments (Supplementary Figures S8.1–S8.3). However, the model scenario that fits the observed data indicates a large enough extinction threshold in our experiment, putting equal selective weight on dispersal and growth rate (Figure 2) and allowing selection to pick the best possible disperser strain that still has a high growth rate. This depicts a combination of pushed and pulled wave scenarios (Miller et al., 2020), where range expansion strongly depends on processes during relatively long growth period just behind the front (typical of pushed waves), but also on dispersal bottlenecks (typical of pulled waves). Unlike dispersal and population growth rate, equilibrium density did not appear to be under direct short-term selection in any of the treatments (see Figure 3 and multiple regression analysis), even though it reached higher levels in the range front treatment (Figure 4). This suggests that the increase in equilibrium density might be an emergent property. Indeed, previous work on another ciliate (Tetrahymena) shows that evolution in the presence of relatively slowly growing food bacteria can promote more prudent consumption rates (Fronhofer & Altermatt, 2015), thus associating lower maximum population growth rate (selected trait) with higher equilibrium density (emergent trait; see also Zilio et al., 2023).

As previously shown for ecological models (Giometto et al., 2014), realistic predictions can be made about the traveling speed of range expansions, at least in controlled laboratory settings. Here, we find that the predictability of certain aspects of range expansions may still hold when evolution is at play. Indeed, our model accurately predicted the observed divergence in dispersal and/or growth characteristics in core versus front treatments. Thus, range core and control treatments and the range front treatment were fixed for different COI genotypes (b05 vs. b07, Supplementary Table S1), while within treatments all populations had the same COI genotype. This finding confirms the complete divergence between treatments and suggests that selection from standing genetic variation is repeatable for independent experimental replicates of the same treatment. However, because certain strains share the same COI marker genotype (Supplementary Table S1), we cannot conclude with certainty that the same strains were fixed across replicates. Nonetheless, despite the lack of resolution, the two fixed COI genotypes correspond to those carried by strains most likely to go to fixation in the model. This is particular evident for the range core and control treatments (fixed for the b05 genotype), where the only two b05 strains in the founder population also were the two most likely winners in the model, due to their high growth rate. For the range front treatment, there is more uncertainty, as 12 strains carry the b07 genotype fixed in this treatment. The most likely winning strain (goe_14) in the model indeed carries the b07 genotype, and among the 12 candidates it is the only strain with both high dispersal and high growth rate, the trait combination favored in the model (Figure 3; Supplementary Table S1). Additional genome sequencing would be required to determine whether these lines are fixed for the same or different (combinations of) strains. One inconsistency between predicted and observed outcomes is that, although correctly identifying the most likely winner strains, our model nonetheless predicts the frequent fixation of strains with alternative COI genotypes (Figure 3). According to the model, our exclusive findings of b05 strains in all nine range core and control lines and b07 strains in all six front lines are highly unlikely (p = .289 < .0001 and .486 = .01, respectively). Possibly, when we determined dispersal and growth of the founder strains, measurement error added to biologically relevant variation (Supplementary Figure S4). This additional noise then cascades through the model, resulting in model predictions with larger error. Alternatively, our model may be missing other factors, such as direct strain–strain interactions or behavioral complexities like density-dependent dispersal, potentially amplifying among-strain variation in performance. Clearly, our experimental design employed a minimalist range expansion scenario that kept the advancing front isolated from the core population. Nonetheless, we believe that both model and data demonstrate that we captured relevant eco-evolutionary dynamics that is more complex than, for example, a simple selection experiment for increased or decreased dispersal. Namely, we allowed for the free interplay between reproductive traits and dispersal, and outcomes are in line with general predictions for real-world range expansions. Our experimental results are also consistent with previous work in another ciliate species, using two-patch systems (Fronhofer & Altermatt, 2015) as well as interconnected multi-patch systems (Fronhofer et al., 2017). One next step toward realism would be to repeat our model-experiment confrontation for more complex experimental landscapes and to test whether we can predict the ecological dynamics (e.g., range expansion speed) or the spatial distribution of genotypes.

Experimental evolution studies show that adaptation to novel conditions may reduce performance in other environments (Kassen, 2014). The emergence of such trade-offs depends on underlying biochemical and life-history constraints (Walsh & Blows, 2009), but also on historical contingency, determining the composition and genetic architecture of the ancestral population, and thus the available trait space for selection. Generally, natural populations present high levels of standing genetic variation, on which selection can act and produce rapid adaptation (Bitter et al., 2019; Chaturvedi et al., 2021). Our case mimicked this situation and generated short-term responses to selection, but no clear signs of trade-offs. In the long run, however, range front and core populations continued to diverge in multiple traits (Figure 4D), and the increase in dispersal in the front treatment was associated with a decrease in growth rate (Figure 4A). Such coupled responses in dispersal and life-history traits are referred to as dispersal syndrome (Clobert et al., 2012; Cote et al., 2017). Typically, they involve the emergence of a competition–colonization trade-off, where dispersal evolution coincides with selection for opportunistic growth strategies (r-selection). Theoretical and empirical studies have demonstrated the importance of dispersal syndromes in generating eco-evolutionary feedbacks and accelerating the pace of range expansions and biological invasions (Burton et al., 2010; Miller et al., 2020; Ochocki et al., 2019; Perkins et al., 2013). Dispersal–growth trade-offs were previously reported for this (Zilio et al., 2023) and another ciliate species (Fronhofer & Altermatt, 2015). In these systems, growth rate is a good indicator of competitive ability, and the trade-off with dispersal likely reflects a true life-history constraint, mediated through energy costs of foraging activity (Fronhofer & Altermatt, 2015). This might be also related to concomitant changes in different aspects of swimming behavior (for details, see Supplementary S2). The evolved differences between core and front lines are stable, even after switching core and front treatments for multiple cycles (Supplementary Figure S6.1). Moreover, mixes of core and front lines readily respond to dispersal selection (Supplementary Figure S6.2), making new evolutionary experiments possible, where phenotypic measurements change can easily be combined with the tracking of COI genotype frequencies.

Our study considers an asexually reproducing organism, where the fixation of (advantageous) allele combinations is not affected by reshuffling through sex and recombination (Lehtonen et al., 2012). Similarly, the plant Arabidopsis thaliana was used in a range expansion experiment, where the fastest-dispersing clonal genotype became predominant in multiple replicate lines, all starting from the same initial mix of clones (Williams et al., 2016). Thus, asexual reproduction narrows down the variability in the range expansion outcomes and, as we show here, makes predictions possible with relatively simple models. Recombination may make predictions more difficult and likely increase the variability of range expansions outcomes, as shown in studies using sexually reproducing organisms (Ochocki & Miller, 2017; Petegem et al., 2018; Weiss-Lehman et al., 2017). For example, recombination may slow down range expansions in the short term, but speed up longer-term responses by creating novel trait associations previously unavailable. Effects of sex in our paramecia may not be readily predictable, due to the nuclear dimorphism typical of ciliates. While creating novel genetic variants (in the germline micronucleus), sexual reproduction also involves the recreation of a new somatic macronucleus, thereby erasing any (somatic) adaptation acquired during asexual life (Verdonck et al., 2021).

Predicting evolution is arduous because of the intrinsic tension between determinism and contingency (Blount et al., 2018), and it demands an adequate theoretical representation of the eco-evolutionary processes in the biological system in question and reliable information on the genetic variation in the relevant traits (Nosil et al., 2020), as we describe in this work. At least in simple settings such as ours, accurate predictions of the short-term evolutionary trait changes during range expansions require surprisingly few parameters, and independent biological realizations can be highly repeatable. Future studies will need to consider, for example, different landscape scenarios and interactions with other species occurring during range expansion. This would imply a more systems-biology approach, with simulations calibrated on the empirical knowledge of the specific ecological scenario and biological players (Duputié et al., 2012; Papp et al., 2011). More generally, increasing our capacities to make reliable quantitative predictions of invasive eco-evolutionary processes is critical to a variety of issues, from conservation and biocontrol strategies to antibiotic development and disease management.

Click here for additional data file.