Genotyping
All founder strains were genotyped for the mitochondrial cytochrome c oxidase I (COI) gene, by extracting DNA from 10 cells per strain and applying a PCR amplification protocol and sequence analysis (56). At cycle 30 in the evolutionary experiment, mixes of 50 cells from each line were processed in the same way and analysed for (multiple) COI marker signals. This method characterises the line for the most frequent COI genotype, and it has a resolution threshold of c. 5%, i.e., it can detect 2-3 cells of a minority genotype in the sample (56). The sequences used are deposited in GenBank (50, 57).
Range expansion modelOur 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 expanded to include both intra- and inter-strain competition:
\begin{equation} \frac{\text{dN}_{\text{i\ }}}{\text{dt}}=(r_{0,i}-\sum_{j}{(\alpha_{\text{ij}}N_{j}}))\ N_{i}\nonumber \\ \end{equation}
where \(N_{i}\) is the population size of strain i, \(r_{0,i}\) is its intrinsic rate of increase and \(\alpha_{\text{ij}}\) as the competition coefficients. The model is parametrized with the posteriors extracted from growth curve fitting (r0, \(\overline{N}\)), 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 (the bottleneck occurring during dispersal, 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 3h dispersal phase in a 2-patch metapopulation. During this 3h 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 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.
Statistical analysisStatistical analyses were performed in R (ver. 4.1.2) and JMP 14 (SAS Institute Inc. 2018). We analysed dispersal (proportion of dispersers), using generalised 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 analysed variation in intrinsic population growth rate (r0) and equilibrium density (\(\overline{N}\); 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 with their respective median values of dispersal, r0, and \(\overline{N}\)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 \(\overline{N}\), we constructed a data matrix based on trait means per year and line (3 years x 15 lines, n = 45), after centering and scaling trait distributions. One range-front line was lost in year 3, leading to n = 44. We used a Bayesian approach with the “rstan” package version 2.21.2 (58) to estimate pairwise trait correlations (SI Appendix, Fig. S3). We also performed a principal component analysis (PCA), considering the joint variation in all three traits.
AcknowledgementsGZ was supported by a grant from the Agence Nationale de Recherche (n° ANR-20-CE02-0023-01) to OK. This is publication ISEM-XXXX-XXX of the Institut des Sciences de l’Evolution. Florent Deshors and Thomas Teissedre assisted with movement and growth assays. Alison Duncan, Flore Zelé and Alexey Potekhin gave helpful comments for the experimental set up and the interpretation of the results.