Abstract
The ability of a species to colonize newly available habitats is crucial to its overall fitness^{1,2,3}. In general, motility and fast expansion are expected to be beneficial for colonization and hence for the fitness of an organism^{4,5,6,7}. Here we apply an evolution protocol to investigate phenotypical requirements for colonizing habitats of different sizes during range expansion by chemotaxing bacteria^{8}. Contrary to the intuitive expectation that faster is better, we show that there is an optimal expansion speed for a given habitat size. Our analysis showed that this effect arises from interactions among pioneering cells at the front of the expanding population, and revealed a simple, evolutionarily stable strategy for colonizing a habitat of a specific size: to expand at a speed given by the product of the growth rate and the habitat size. These results illustrate stabilitytoinvasion as a powerful principle for the selection of phenotypes in complex ecological processes.
This is a preview of subscription content
Access options
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Rent or Buy article
Get time limited or full article access on ReadCube.
from$8.99
All prices are NET prices.
Data availability
Sequencing data have been deposited to the NCBI Sequence Read Archive (SRA), accession PRJNA559221. Other major experimental data supporting the findings of this study are available within the paper and Supplementary Information. Simulation data can be generated with the custommade code and the parameter sets provided.
Code availability
Custommade simulation code is available via GitHub at https://github.com/jonascremer/chemotaxis_simulation.
References
 1.
Hanski, I. Metapopulation dynamics. Nature 396, 41–49 (1998).
 2.
Skellam, J. G. Random dispersal in theoretical populations. Biometrika 38, 196–218 (1951).
 3.
Andow, D. A., Kareiva, P. M., Levin, S. A. & Okubo, A. Spread of invading organisms. Landsc. Ecol. 4, 177–188 (1990).
 4.
Yi, X. & Dean, A. M. Phenotypic plasticity as an adaptation to a functional tradeoff. eLife 5, e19307 (2016).
 5.
Fraebel, D. T. et al. Environment determines evolutionary trajectory in a constrained phenotypic space. eLife 6, e24669 (2017).
 6.
Ni, B. et al. Evolutionary remodeling of bacterial motility checkpoint control. Cell Rep. 18, 866–877 (2017).
 7.
Shih, H.Y., Mickalide, H., Fraebel, D. T., Goldenfeld, N. & Kuehn, S. Biophysical constraints determine the selection of phenotypic fluctuations during directed evolution. Phys. Biol. 15, 065003 (2018).
 8.
Adler, J. Chemotaxis in bacteria. Science 153, 708–716 (1966).
 9.
Levin, S. A. The problem of pattern and scale in ecology: the Robert H. MacArthur award lecture. Ecology 73, 1943–1967 (1992).
 10.
Hastings, A. et al. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett. 8, 91–101 (2005).
 11.
Hallatschek, O., Hersen, P., Ramanathan, S. & Nelson, D. R. Genetic drift at expanding frontiers promotes gene segregation. Proc. Natl Acad. Sci. USA 104, 19926–19930 (2007).
 12.
Müller, M. J. I., Neugeboren, B. I., Nelson, D. R. & Murray, A. W. Genetic drift opposes mutualism during spatial population expansion. Proc. Natl Acad. Sci. USA 111, 1037–1042 (2014).
 13.
Cao, Y. et al. Collective spacesensing coordinates pattern scaling in engineered bacteria. Cell 165, 620–630 (2016).
 14.
Cremer, J. et al. Chemotaxis as navigation strategy to boost range expansion. Nature https://doi.org/10.1038/s415860191733y (2019).
 15.
Elena, S. F. & Lenski, R. E. Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat. Rev. Genet. 4, 457–469 (2003).
 16.
Bosshard, L. et al. Accumulation of deleterious mutations during bacterial range expansions. Genetics 207, 669–684 (2017).
 17.
Wolfe, A. J. & Berg, H. C. Migration of bacteria in semisolid agar. Proc. Natl Acad. Sci. USA 86, 6973–6977 (1989).
 18.
Deforet, M., CarmonaFontaine, C., Korolev, K. S. & Xavier, J. B. Evolution at the edge of expanding populations. Am. Nat. 194, 291–305 (2019).
 19.
Lenski, R. in Microbial Ecology: Principles, Applications and Methods (eds Levin, M. et al.) 183–198 (McGrawHill, 1992).
 20.
Smith, J. M. Evolution and the Theory of Games (Cambridge Univ. Press, 1982).
 21.
Levin, B. R. Frequencydependent selection in bacterial populations. Phil. Trans. R. Soc. Lond. B 319, 459–472 (1988).
 22.
Alon, U., Surette, M. G., Barkai, N. & Leibler, S. Robustness in bacterial chemotaxis. Nature 397, 168–171 (1999).
 23.
Hansen, C. H., Endres, R. G. & Wingreen, N. S. Chemotaxis in Escherichia coli: a molecular model for robust precise adaptation. PLOS Comput. Biol. 4, e1 (2008).
 24.
Korobkova, E., Emonet, T., Vilar, J. M., Shimizu, T. S. & Cluzel, P. From molecular noise to behavioural variability in a single bacterium. Nature 428, 574–578 (2004).
 25.
Park, H., Guet, C. C., Emonet, T. & Cluzel, P. Finetuning of chemotactic response in E. coli determined by highthroughput capillary assay. Curr. Microbiol. 62, 764–769 (2011).
 26.
Si, G., Wu, T., Ouyang, Q. & Tu, Y. Pathwaybased meanfield model for Escherichia coli chemotaxis. Phys. Rev. Lett. 109, 048101 (2012).
 27.
Tu, Y. Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation. Annu. Rev. Biophys. 42, 337–359 (2013).
 28.
Liu, C. et al. Sequential establishment of stripe patterns in an expanding cell population. Science 334, 238–241 (2011).
 29.
Merrell, D. The Adaptive Seascape: The Mechanism of Evolution (Univ. Minnesota Press, 1994).
 30.
Mustonen, V. & Lässig, M. From fitness landscapes to seascapes: nonequilibrium dynamics of selection and adaptation. Trends Genet. 25, 111–119 (2009).
 31.
Poelwijk, F. J., de Vos, M. G. & Tans, S. J. Tradeoffs and optimality in the evolution of gene regulation. Cell 146, 462–470 (2011).
 32.
Towbin, B. D. et al. Optimality and suboptimality in a bacterial growth law. Nat. Commun. 8, 14123 (2017).
 33.
Hofbauer, J. & Sigmund, K. Evolutionary Games and Population Dynamics (Cambridge Univ. Press, 1998).
 34.
Gore, J., Youk, H. & van Oudenaarden, A. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256 (2009).
 35.
Zheng, H. et al. Interrogating the Escherichia coli cell cycle by cell dimension perturbations. Proc. Natl Acad. Sci. USA 113, 15000–15005 (2016).
 36.
Datta, S., Costantino, N. & Court, D. L. A set of recombineering plasmids for Gramnegative bacteria. Gene 379, 109–115 (2006).
 37.
Jiang, Y. et al. Multigene editing in the Escherichia coli genome via the CRISPR–Cas9 system. Appl. Environ. Microbiol. 81, 2506–2514 (2015).
 38.
Fu, J. L., Kanno, T., Liang, S.C., Matzke, A. J. M. & Matzke, M. GFP lossoffunction mutations in Arabidopsis thaliana. G3 5, 1849–1855 (2015).
 39.
Waite, A. J. et al. Nongenetic diversity modulates population performance. Mol. Syst. Biol. 12, 895 (2016).
 40.
Edelstein, A. D. et al. Advanced methods of microscope control using μManager software. J. Biol. Methods 1, e10 (2014).
 41.
Dufour, Y. S., Gillet, S., Frankel, N. W., Weibel, D. B. & Emonet, T. Direct correlation between motile behavior and protein abundance in single cells. PLOS Comput. Biol. 12, e1005041 (2016).
 42.
Deatherage, D. E. & Barrick, J. E. Identification of mutations in laboratoryevolved microbes from nextgeneration sequencing data using breseq. Methods Mol. Biol. 1151, 165–188 (2014).
 43.
Guyer, J. E., Wheeler, D. & Warren, J. A. FiPy: partial differential equations with Python. Comput. Sci. Eng. 11, 6–15 (2009).
 44.
Fisher, R. The wave of advance of advantageous genes. Ann. Eugen. 7, 355–369 (1937).
 45.
Kolmogorov, A., Petrovsky, I. & Piscounov, N. Étude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique. Mosk. Univ. Bull. Math. 1, 37 (1937).
 46.
Korolev, K. S. Evolution arrests invasions of cooperative populations. Phys. Rev. Lett. 115, 208104 (2015).
 47.
Yang, F., Moss, L. G. & Phillips, G. N. Jr. The molecular structure of green fluorescent protein. Nat. Biotechnol. 14, 1246–1251 (1996).
 48.
Keller, E. F. & Segel, L. A. Model for chemotaxis. J. Theor. Biol. 30, 225–234 (1971).
 49.
Vaknin, A. & Berg, H. C. Physical responses of bacterial chemoreceptors. J. Mol. Biol. 366, 1416–1423 (2007).
 50.
Taylor, J. R. & Stocker, R. Tradeoffs of chemotactic foraging in turbulent water. Science 338, 675–679 (2012).
 51.
Fu, X. et al. Spatial selforganization resolves conflicts between individuality and collective migration. Nat. Commun. 9, 2177 (2018).
Acknowledgements
We thank L. Chao, X. Fu, X. He, J.D. Huang, A. Murray, M. Vergassola, C.I. Wu and G. Zhao for discussions, and Y. Wu and H. Zhou for assistance with bioinformatic analyses. C.L., W.L. and D.L. acknowledge financial support by the Major Research Plan of the National Natural Science Foundation of China (91731302), National Key Research and Development Program of China (2018YFA0902700), Strategic Priority Research Program (XDB29050501), Key Research Program (KFZDSW216) of Chinese Academy of Sciences, and Shenzhen Grants (JCYJ20170818164139781, KQTD2015033117210153, Engineering Laboratory [2016]1194). T.H. and J.C. acknowledge support from the NIH through grant R01GM95903.
Author information
Affiliations
Contributions
C.L. initiated and directed the research. W.L. and D.L. carried out most of the experiments, with contributions from C.L. and J.C. J.C. and T.H. developed the model and carried out the numerical simulations and mathematical analysis. All authors analysed the results. J.C., T.H. and C.L. wrote the manuscript with contributions from W.L.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
Extended Data Fig. 1 Characteristics of cell motility and expansion speed for the ancestor and for mixed populations of evolved cells.
a, Navigated range expansion by bacteria^{14} on a soft agar plate involves a group of pioneering cells (cyan) at the population front that move outwards towards uncolonized space (grey) in response to a signalling cue. During the outward migration, they replicate and leave offspring behind the front. The offspring do not move outwards, but settle wherever they are deposited and grow exponentially (green) until they reach carrying capacity (brown). b, Density profile of a population of ancestor cells containing GFP (Anc_{G}) at various times after inoculation at the centre of a semisolid TB plate with 0.25% agar incubated at 37 °C. The population expanded at a defined expansion speed, covering positions A–E well before selection at 24 h. The density profiles were obtained using confocal microscopy as described previously^{14}. Cellular characteristics were not significantly affected by the expression of GFP (Extended Data Fig. 3). c, The experiment described in Fig. 1a was carried out independently for five distinct selection positions, at distances X_{A}, X_{B}, X_{C}, X_{D} and X_{E}, ranging from 5 to 25 mm from the centre. Three independent lineages of this experiment were propagated in parallel (Extended Data Fig. 2a). Selected samples at cycle n of series S in lineage l are referred to as Sl.n in Supplementary Table 2. d, Front position versus time for the ancestral strain CLM (purple) and populations of evolved cells obtained at the 50th cycle of each of the five evolution series from lineage 1. For each mixed population, collected cells were grown to midexponential phase (OD_{600} = 0.20), and 2 μl of the batch culture was inoculated at the centre of the same TB plate as in b and incubated in the same way. The agar plates were photographed at different times after inoculation and the radius of each expanding population was deduced from the area measured using ImageJ. The expansion speed of each population (Fig. 1b) was obtained as the slope of the linear fit of the data after t = 3 h. e, Growth rates of population samples from each of the five selection series from lineage 1 at various evolution cycles. f, Growth rates of population samples from each of the three selection series (lineage 1) at various cycles of evolution experiments in CAA medium. g, Scatter plot of growth rates and expansion speeds for single strains isolated from frozen samples taken at various cycles of CAA evolution experiments. The evolution cycles are indicated by circles. h, i, Mean run speed (h) and tumble bias (i) of population samples from each of the three selection series from lineage 1 of evolution experiments carried out in TB medium at various evolution cycles. At least 10,000 cells were subjected to singlecell tracking analysis for each experiment (see Methods). Error bars represent s.d. j, Expansion speeds of several strains, each carrying an identified mutation, plotted against the corresponding evolved strains from which the mutations were identified. The mutations are indicated in the legend (see Supplementary Table 1). mutH does not affect motility and serves as a control. Experiments were repeated independently three times for b, d and twice for e, f with similar results. In h, i, data are mean ± s.d. for a single biological replicate, n = 10,000 cells analysed. In j, data are means for n = 3 biological replicates (s.d. error bars are smaller than the symbols).
Extended Data Fig. 2 Experimental evolution of expansion speed in different growth media.
a–d, Expansion speeds of population samples from each selection series are shown at various cycles of evolution experiments carried out in TB (a), LB (b), CAA (c) and glycerolcontaining minimal medium (d). The agar density was 0.25% (w/v). All three lineages of each medium showed similar trajectories. The absence of a chemoattractant in the glycerol minimal medium leads to a very different outcome compared to the experiments carried out in complex media (a–c). In the absence of a chemoattractant, cells follow simple Fisher–Kolmogorov dynamics^{14,44,45}. In line with the absence of a growing population trailing the front, no decrease in swimming behaviour was observed over time. However, as previously investigated, slower swimming behaviour might be selected for when densitydependent growth effects or strong tradeoffs between swimming and cell growth exist^{18,46}. Experiments in a, b, d were carried out once with three biologically independent repeats with similar results; experiment in b was repeated independently twice with similar results.
Extended Data Fig. 3 Fitness effects for cells expressing GFP and its nonfluorescent variant NFP.
a, b, Growth rates (a) and expansion speeds (b) of the ancestor and six mutant clones harbouring GFP, NFP, or no plasmid, along with the corresponding mixed populations from which the six mutant clones were isolated. Plasmids harbouring constitutive GFP or NFP expression (PZA31PtetM2GFP or PZA31PtetM2NFP, respectively) were transformed into the six evolved strains A, B, B′, C, D, and E, to form A_{G}, A_{N} and so on (Supplementary Table 2). c, Fluorescence intensity as a function of cell density. The cell growth of cells expressing GFP was arrested by adding 2 mg ml^{−1} kanamycin, and cells were concentrated to 1.6 × 10^{10} cells per ml. Subsequently, serial dilutions were carried out. For each cell density, cells were vigorously mixed with prewarmed 0.277% (w/v) TB agar containing 2 mg ml^{−1} kanamycin at a ratio of 1:9 and poured into three Petri dishes with 10 ml each. All dishes were allowed to harden at room temperature for 90 min. The cell–agar mixture was subjected to cell counting by fluorescenceactivated cell sorting (FACS) and the fluorescence intensity of the agar plate was detected using a fluorescence microscope. d, The threedimensional structure of GFP with the predicted position of the lossoffunction mutation introduced, Y66C. The shown crystal structure is based on the NFP protein sequence aligned to PDB ID: 1GFL^{47}. The rectangle shows the mutation position. e, Ancestor cells harbouring constitutive GFP or NFP expression (Anc_{G} and Anc_{N} cells, respectively) as viewed by fluorescence microscopy; the images verify the lossoffunction mutation of GFP. f, The relative fitness W of Anc_{G} and Anc_{N} cells at different distances from the agar plate centre. Cells harbouring the GFP or NFP plasmid were equally mixed and inoculated at the centre of 0.25% TB plates. Cell–agar mixtures picked at various distances were subjected to cell counting by FACS 24 h after inoculation and the ratio was reported as W (see Methods). g, h,The genetic circuit and characteristics of the cheZtitratable strains WL1 (g) and WL2 (h). The expression of cheZ is under the control of a P_{tet}tetR (WL1) or P_{lac}lacI (WL2) feedback loop and the native cheZ was seamlessly removed (see Methods). Relative cheZ mRNA levels change around twofold under different concentrations of inducers (aTc for WL1 and IPTG for WL2). Relative cheZ mRNA expression levels, expansion speed, and growth rates of WL1 (g) and WL2 (h) under various concentrations of the respective inducers are shown next to the circuits. Horizontal dashed lines show the corresponding values for the ancestor strain. i–l, Growth rates (i, k) and expansion speeds (j, l) of the aTctitratable cheZ strain WL1 (i, j) and its derivatives expressing GFP (WL7) or NFP (WL8), and the IPTGtitratable cheZ strain WL2 (k, l) and its derivatives expressing GFP (WL9) or NFP (WL10). m, Circles indicate crossover distances between fluorescent derivatives of two cheZtitratable strains (WL1 and WL2), with WL1 strains induced at a fixed concentration of its inducer (0.5 ng ml^{−1} aTc), and WL2 strains induced at various IPTG concentrations (expansion speeds (h, l) shown on the x axis). The background colours again indicate dominance by WL1 (purple) or WL2 (green). In a, b, g–l, data are mean ± s.d. for n = 3 biologically independent repeats (individual data points shown as circles). Error bars in g, h were smaller than the symbols. Experiments shown in c, e, f, m were repeated independently three times with similar results.
Extended Data Fig. 4 Results of twostrain competitions.
a, Expansion speeds and growth rates of ancestor, strain B, and strain D (Supplementary Table 2). b, Calibrated fluorescence intensity profiles of singly grown strains 12 h after inoculation at the centre of 0.25% TB agar plates. Unless noted otherwise, fluorescence intensity is normalized according to the calibration curve shown in Extended Data Fig. 3c; the relative value 1 refers to 5 × 10^{7} cells per ml. c, Relative fluorescence intensities obtained at various times for the fluorescent mutant strains B_{G} (orange) and D_{G} (cyan) and the ancestor Anc_{G} (purple), each grown singly on TB plates. The faster strain has higher fitness everywhere. d, Raw photographs (top) and fluorescence intensity profiles (bottom) before and after merging of a representative twostrain competition between the fluorescent derivatives of the ancestor and strain D 12 h after initial equal inoculation. We used plasmids GFP and NFP in this study to distinguish the two strains from each other in the headtohead competition, as there is no systematic influence caused by the expression of GFP or NFP (Extended Data Fig. 3). Top (from left to right): competition between D_{G} and Anc_{N}, competition between D_{N} and Anc_{G}, and merged photograph in pseudocolour (D_{G}, cyan; Anc_{G}, purple). Bottom: corresponding relative fluorescence intensity profiles. e, As in d, but with strain B instead of strain D. Evolved strains and the ancestor with GFP or NFP were cultured to log phase separately. Two types of the combined mixed strains were prepared and the evolved strain with GFP was mixed with Anc_{N} while the evolved strain with NFP was mixed with Anc_{G} in a 1:1 ratio. Subsequently, 2 μl of the two types of combined mixture was inoculated onto the centre of prewarmed semisolid agar plates separately and allowed to expand at 37 °C for up to 24 h. Photographs and fluorescent intensities of the evolved strains and ancestor with the GFP reporter from these two plates after the expansion were taken at various times and merged (see Methods). f, Relative fitness W_{B} obtained as a ratio of the direct cell count of the fluorescent derivatives of the ancestor and strain B, inoculated at the centre of an agar plate. From left to right: competition between B_{G} and Anc_{N}, competition between B_{N} and Anc_{G}, and averaged curve of both (the relative fitness W_{B}). The competition experiments are the same as in e. The initial and final ratios of the two competitors were measured by cell counting with a flow cytometer (see Methods). g, i, The spatiotemporal development of the bacterial density profiles, indicated by fluorescence intensities, for the competition between the ancestor (purple) and strain D (cyan; g) or strain B (orange; i) at various times during the 24h competition. Beyond the initial period (~6 h), the crossover distance could be clearly defined (vertical dashed line) and was practically timeindependent, with the slower strain gaining advantage in the interior and the faster strain gaining advantage in the exterior. h, j, Relative fitness values W_{D} (h) and W_{B} (j) taken as the ratio of the fluorescence intensities (mutant:ancestor) at various distances for 6 h and beyond. The vertical dashed lines indicate the crossover distance d_{×} where the densities of the two competing strains are equal (W_{i} = 1). Thus, the crossover distance was fixed shortly after the initial period. k–m, Fluorescence intensity profiles (k, l), and relative fitness W (m) of representative twostrain competitions between the ancestor (black solid line) and evolved isolates. The data were taken 12 h (k, m) or 24 h (l) after coinoculation of equal initial mixtures of the two competing strains at the centre of 0.25% TB agar plates, showing that the slower strain spatially outcompetes the faster strain within the crossover distance d_{×} (dashed lines). See Supplementary Table 2 for strain information. Experiments in a–l were repeated independently three times with similar results. In m, the mean ± s.d. of n = 3 biologically independent repeats is shown.
Extended Data Fig. 5 The growth–expansion dynamics of a single strain.
a, We review here the GE model^{14} generated to describe the dynamics of a single strain of E. coli colonizing a soft agar plate. The GE model considers three variables and their dynamics in space x and time t: (i) cell density ρ(x,t); (ii) the concentration a(x,t) of an attractant that cells can sense and move towards; and (iii) the concentration n(x,t) of a nutrient that cells consume to grow. Following the spirit of the classical model introduced by Keller and Segel^{48} (the KS model), cells can move in a random, undirected manner (effective diffusion constant D and diffusion term in equation (S1)) and along the gradient of the signalling molecule (the attractant gradient represented by the convection term in equations (S1) and (S2)), which is generated by cellular consumption (equation (S3)). To account for observed density profiles and their evolution over time, three additional aspects beyond the KS model are important to consider. First, cells can detect and respond to attractant gradients in a scalefree manner only within a limited range between K_{I} and K_{A}, the lower and upper cutoffs describing the molecular limitations of attractant sensing^{26,49} as specified in equation (S2). D_{0} denotes the molecular diffusion of the attractant and nutrient within the agar. The chemotactic parameter c denotes how cells translate the detected attractant gradient into directed movement. Second, cells grow throughout the expansion process, as described by the growth term in equation (S1), with growth rate λ. Third, growth relies not on the presence of the attractant but on the presence of nutrients. We model the latter dependence by the nutrient field n(x,t) and a yield factor Y (equation (S4)). The distinct treatment of the roles of nutrient and attractant is designed to model the dynamics of bacterial cells in complex media, where nonchemotactic components in the complex media are designated as the nutrient, while the (minor) chemotactic components of the complex media are reflected by a low concentration of a single attractant; a detailed discussion and validation of the model, including comparison to other models of chemotactic migration, has been published previously^{14}. b, Emerging density profile (green solid line) and attractant concentration (magenta dashed line) in the GE model. At the front is a density bulge or peak, within which the attractant profile drops steeply, guiding cells to do chemotaxis. Directed movement of cells in the front bulge is coupled to an exponentially rising trailing region. In this region, cells grow and swim randomly, but there is no directed cell movement there owing to the low attractant concentration (a < K_{I}). Note that in cases in which multiple attractants are present in the medium, the population typically exhibits multiple rings, one corresponding to the exhaustion of each attractant. For these cases, the trailing region involving cell growth but no chemotaxis would correspond to the region inside the innermost ring, after the exhaustion of the last attractant. Thus, our model with a single attractant does not attempt to describe the movement of the outer rings, but models the transition of the density profile from the innermost ring to the exponential trailing region. As we describe in Supplementary Analyses 1, 2, the dynamics in this transition region determine strain dominance in multistrain competition processes. c, Stable expansion of the population can be explained by a balance between growth of cells in the front bulge and leakage of cells out of the front due to random movement (cell diffusion). For illustration, consider the green dashed line indicating the density profile at a given earlier time. With only directed movement (along attractant gradient) and cell growth, the front peak would be higher at a later time, as illustrated by the brown dashed line. Instead, growth in the front bulge is compensated by backdiffusion of cells away from the front (green solid line indicating a later time). This compensation mechanism is further shown in the cylinder plots below, illustrating that the observed stable expansion dynamics (constant expansion speed, constant peak density) results from a dynamic balance between growth at the front and backdiffusion. The exponential density profile in the trailing region results from a combination of the steady outward movement of the front and a steady exponential growth of cells leaked out of the front^{14}.
Extended Data Fig. 6 Model of twostrain competitive expansion dynamics.
a, The GE model of bacterial range expansion^{14} (Extended Data Fig. 5) is extended to several strains of bacteria for which the densities are denoted by ρ_{i}(x,t), with i ∈ {1,2} for two strains. The different bacterial strains are assumed to consume the same nutrient (n) and grow at the same rate λ(n), in accordance with Monod’s law. They also sense and consume the same signalling molecule (the attractant (a)) at the same rate μ(a). The random motion of each bacterial strain is described by a diffusion term, with effective diffusion coefficients D_{i} for strain i (equations (S5), (S6)). The spatial attractant profile a(x,t) (resulting from bacterial consumption, equation (S8)) leads to directed motion which is modelled by a drift term in equations (S5) and (S6). The dependence of the drift velocities v_{i} on the attractant profile is given by equation (S7), which describes a range of proportional sensing v ∝ ∂_{x}a/a for K_{I} < a < K_{A}. The magnitude of the chemotactic response is parametrized by the chemotactic coefficient c_{i} for strain i. Finally, the dynamics of the nutrient and attractant are described by equations (S8) and (S9), with D_{0} characterizing molecular diffusion. b, Outcomes of competitive expansion dynamics, showing the density profiles of two strains 24 h after inoculation with equal mixtures at the origin. Competition is run for one strain resembling the ancestor (black line, with expansion speed u_{anc} = 6 mm h^{−1}) and the other strain resembling a mutant, with expansion speed u_{mut} increasing from left to right as indicated above the plots. Experimental data are from Extended Data Fig. 4l. c, Increase and abrupt freezing of the crossover distance over time as observed in simulations (line) and experiments (circles). The expansion speeds of the competing pair are u_{anc} versus u_{mut} = 6 mm h^{−1} versus 7 mm h^{−1}, with the latter mimicking the expansion speed of strain D (Extended Data Fig. 3b). Full temporal evolution of the density profiles is shown in Supplementary Video 2. Experimental data are from Extended Data Fig. 4g. d, The chemotactic competition model in a is repeated to determine the stability distance d* for various u and several smaller growth rate values λ. The results confirm that the linear relation in equation (2) between the stability distance and expansion speed shown in the main text holds for each growth rate simulated. e, The slope u/d*(u) in d is plotted against the growth rate λ. A linear dependence on λ is seen as predicted by Supplementary Analysis 2. f, Effect of differing growth rate (x axis) on the outcome of competition and crossover distance for two strains with different chemotactic coefficients (c_{1}, c_{2}). The main result that the interior is dominated by the slower strain and the exterior by the faster strain holds. The crossover distance that separates the different regions of dominance varies smoothly with growth rate differences, but in opposite ways depending on whether the faster growing strain moves faster (red) or slower (green). The observed changes in crossover distances are minimal when growthrate differences are minimal (for example, <5% as we observed experimentally for TB; Extended Data Fig. 1e) but become substantial when growthrate differences become large (Supplementary Analysis 3). g, The stability distance d* as in d but for an alternative model formulation of chemotactic movement^{50} with different functions describing sensing and the directed movement along gradients. See Supplementary Analysis 4 for more details. Despite the model changes, a linear relation with growthratedependent slopes was still observed. In b–g, the chemotactic coefficient c was varied to modify expansion speeds. The cellular diffusion coefficient D was varied accordingly, with D = c/6.25 remaining constant. Simulations with a fixed diffusion coefficient gave similar results. Experiments in b, c were repeated independently three times with similar results.
Extended Data Fig. 7 Twostrain competition and crossover dynamics.
The results in Fig. 2 show that competition between two strains with the same growth rate and different singlestrain expansion speeds generically gives rise to a distinct spatial structure where the slower strain dominates inside and the faster strain dominates outside, separated by a crossover distance that depends on the two speeds. Here we explain how this feature arises from the underlying GE dynamics. a, Dynamics simulation using equations (S5)–(S9) for two strains differing only in their chemotactic coefficients (c_{slow} = 652 μm^{2} s^{−1} and c_{fast} = 727 μm^{2} s^{−1}, with corresponding singlestrain expansion speeds u_{slow} = 4.95 mm h^{−1} and u_{fast} = 5.17 mm h^{−1}; other parameters, equal for both strains, are provided in Supplementary Table 3. The density profile of the strain with faster singlestrain expansion speed is green, the other purple. Red circles indicate crossover distances, which freeze at a fixed position in the laboratory frame (dashed red line) after the troughs of the two density bulges cross. This feature (see also Extended Data Fig. 4g,i) allows us to define a simple crossover distance d_{×} that is timeinvariant. To understand the crossover dynamics in a, we first provide a qualitative explanation. While the two strains would individually expand at different speeds, the competition dynamics involves a single comigrating population. This is because the two strains chase after the same attractant profile, which can recede only at a single speed^{51}. Because the front is moving faster than the speed at which the slower strain can stably propagate, the slower strain is gradually depleted from the front. As shown by the cylinders, at early times (before the front has reached the crossover distance and where the abundances of the two strains in the front bulge are comparable), the leakage flux of slower cells (purple) at the back exceeds that of the faster cells (green). As the two strains grow at the same rate behind the front, the slower strain dominates there. Because of its faster leakage, the slower strain is preferentially depleted from the front bulge. Subsequently, the leakage flux of the slower strain drops below that of the faster strain despite the faster leakage rate of the slower strain, owing to its reduced abundance. From there on, the back is dominated by the faster strain. At the crossover distance, the two leakage fluxes are the same. b, A simple analysis captures key features of the crossover dynamics and leads to a timeinvariant crossover distance d_{×}(u_{fast}, u_{slow}) for two strains with singlestrain expansion speeds u_{fast} and u_{slow}, and shows how this crossover distance leads to an evolutionarily stable selection distance d* defined in the limit u_{slow}→u_{fast} (Figs. 3, 4). The two density profiles shown in a, ρ_{fast}(x, t) and ρ_{slow}(x, t) in the ‘laboratory frame’ coordinate x, are shown in b in the frame moving with speed u_{fast}. In this moving frame where the spatial coordinate is y = x − u_{fast}t, the density profile of the faster strain (green line) is approximately stationary: \({\rho }_{{\rm{fast}}}(x,t)={\hat{\rho }}_{{\rm{fast}}}(y)\) where \({\hat{\rho }}_{{\rm{fast}}}(y)\) is the stationary solution of the singlestrain dynamics in the reference frame moving at speed u_{fast}. Defining the position of the trough to be y = 0 in the moving frame, \({\hat{\rho }}_{{\rm{fast}}}(y)\) has a bulge for y > 0 and a trailing exponential for y < 0. In this frame, the density profile of the slower strain is generally not stationary, and is written as \({\rho }_{{\rm{slow}}}(x,t)={\hat{\rho }}_{{\rm{slow}}}(y,t)\). As described in Supplementary Analysis 1, because the slower strain expands faster than its singlestrain expansion speed its density bulge, while preserving the spatial structure, is depleted exponentially over time. This is expressed mathematically as \({\hat{\rho }}_{{\rm{s}}{\rm{l}}{\rm{o}}{\rm{w}}}(y,t)={\mathop{\rho }\limits^{ \sim }}_{{\rm{s}}{\rm{l}}{\rm{o}}{\rm{w}}}(y){{\rm{e}}}^{\gamma t}\) where \({\mathop{\rho }\limits^{ \sim }}_{{\rm{slow}}}(y)\) also has a bulge for y > 0 and a trailing exponential for y < 0. Because the density bulges of the two strains are aligned spatially by the common attractant gradient, we denote the position of the bulge peak by y* for both strains. However, in general neither the peak height nor the slope of the trailing exponentials would be the same between the two strains. In b, the density profile of the slower strain, \({\mathop{\rho }\limits^{ \sim }}_{{\rm{s}}{\rm{l}}{\rm{o}}{\rm{w}}}(y){{\rm{e}}}^{\gamma t}\), is illustrated at three times t. Dotted purple line indicates initial time (t = 0) where the density bulge of the slower strain has the same peak value ρ* as the faster strain (green line): \({\mathop{\rho }\limits^{ \sim }}_{{\rm{s}}{\rm{l}}{\rm{o}}{\rm{w}}}({y}^{\ast })={\rho }^{\ast }\). Dashed purple line indicates density profile of the slower strain at time t_{×} where the two density troughs cross: \({\mathop{\rho }\limits^{ \sim }}_{{\rm{s}}{\rm{l}}{\rm{o}}{\rm{w}}}(0){{\rm{e}}}^{\gamma {t}_{\times }}={\rho }_{0}\). The location at which the troughs cross (open red circle) is the crossover distance d_{×} = u_{fast}t_{×} (obtained from y(t_{x}) = 0). For t > t_{×}, the density bulge of the slower strain sinks steadily below that of the faster strain, and its density profile (solid purple line) crosses that of the faster strain (green line) at y_{×} < 0 behind the front and falls steadily backward in the moving frame (filled red circle). Supplementary Analysis 1 shows that for a steady sinking rate γ ∝ (u_{fast} − u_{slow}), the corresponding crossover position in the laboratory frame remains at d_{×} for all t > t_{×} (fixed position of red circle in b). This feature of the crossover dynamics can be understood intuitively: as chemotaxis occurs only within the bulge region, cells at the crossover point (and to its left) experience no net drift. As the two strains grow at the same rate, their densities at this position remain equal to each other for all t > t_{×}, implying that the density crossover is frozen in at d_{×}. The mathematical analysis provides an expression (Supplementary equation (E14); see Supplementary Information) of how the crossover distance d_{×} = u_{fast}t_{×} depends on u_{fast} and u_{slow} and on static properties of the two density profiles, \({\hat{\rho }}_{{\rm{fast}}}(y)\) and \({\tilde{\rho }}_{{\rm{slow}}}(y)\). From this, we can obtain an expression for the stability distance, defined in the limit the two speeds approach each other: d*(u) = ut_{×}(u_{slow}→u). Supplementary Analysis 2 shows that t_{×}(u_{slow}→u) approaches a finite limit that is proportional to 1/λ, the doubling time. The proportionality of d* and u is verified in Extended Data Fig. 6d, and the growth rate dependence of the slope in Extended Data Fig. 6e.
Extended Data Fig. 8 Fitness landscapes for the competitive expansion process at different positions.
The positiondependent fitness landscapes for the competition between two strains (strains 1 and 2, for example, an ancestor and a mutant), with expansion speeds u_{1} and u_{2}, respectively, computed according to the competitive expansion model in Extended Data Fig. 6a. The fitness of the mutant (strain 2) relative to the ancestor (strain 1) is defined as the ratio of its density with respect to the ancestor: \(W(X)\equiv {\mathrm{lim}}_{t\to \infty }\left[\frac{{\rho }_{2}(X,t)}{{\rho }_{1}(X,t)}\right]\). In practice, we take t = 24 h, which is well after the dynamics have halted. We further normalize this fitness by the ratio of the initial inoculant, \({W}_{0}\equiv \frac{{\rho }_{2}(X,0)}{{\rho }_{1}(X,0)}\), to obtain the ‘relative fitness gain’ W/W_{0}. In a, b, green indicates a gain in the fitness of the mutant (W/W_{0} > 1) and blue indicates a loss (W/W_{0} < 1). Top, middle, and bottom rows refer to results at distances X = 20, 15 and 10 mm, respectively. a, Landscape W/W_{0} for equal abundance (50:50) of the two strains at equal inoculation (W_{0} = 1). Cyan and dashed represent W/W_{0} = 1 and u_{1} = u_{2}, respectively. For W_{0} = 1, the cyan line is the crossover distance: X = d_{×}(u, u′) (corresponding to W = 1) (Fig. 4). The middle row (X = 15 mm) is an expanded view of Fig. 4c, except that the phase diagram in Fig. 4c shows only qualitative information on which strain dominates, whereas the fitness landscape here also provides quantitative information on the fitness gain. Top and bottom, corresponding fitness landscapes for X = 20 and 10 mm, respectively. b, Fitness landscape calculations for a mutant:ancestor inoculant ratio of 5:95 (W_{0} ≡ ρ_{2}(X, 0): ρ_{1}(X, 0) = 5:95). The solid brown line still indicates the phase boundary W/W_{0} = 1, but it is no longer directly related to the crossover distance (see below). The topological features of the landscapes in a and b are similar but the details differ, reflecting the frequencydependent nature of the competition process. c, Landscape profiles for three special ancestor expansion speeds, one at each position X indicated by dotted grey lines marked as ① ② and ③ in b, each obtained as the intersection of the phase boundary and the diagonal (solid brown and dashed pink lines, respectively). For X = 15 mm (middle), ancestral speed \({u}_{1}^{\ast }=5.15\,{{\rm{m}}{\rm{m}}{\rm{h}}}^{1}\) marked by grey line ②, the fitness landscape has a maximum at \({u}_{2}={u}_{1}^{\ast }\) (white circle), indicating that ancestral expansion speed \({u}_{1}^{\ast }\) is stable to invasion by mutants with smaller or larger expansion speeds. Similarly, for X = 20 mm (top), the ancestral speed \({u}_{1}^{\ast }=6.77\,{{\rm{m}}{\rm{m}}{\rm{h}}}^{1}\) (①) is the stable expansion speed, whereas for X = 10 mm (bottom), the ancestral speed \({u}_{1}^{\ast }=3.63\,{{\rm{m}}{\rm{m}}{\rm{h}}}^{1}\) (③) is the stable expansion speed. d, The phase boundaries shown as solid cyan and brown lines in a, b intersect the diagonal at the same speed u* for each X. Thus, the stable expansion speeds can be obtained from the crossover distance (cyan line) based on 50:50 inoculant (equation (3)), despite the frequency dependence of the competition process. In other words, even though the overall fitness W depends on the initial frequency W_{0} in a complex way, in the vicinity of the diagonal (u_{1} ≈ u_{2}), W is independent of W_{0} so that the illustration in Fig. 3a of the stability of a strain of speed u at distance d_{×}(u, u) is verified explicitly in c, d. To vary expansion speeds in each panel, the chemotactic coefficients c_{i} were varied within the range 180–1,500 μm^{2} s^{−1}. Diffusion coefficients were changed accordingly (D_{i} = c_{i}/6.25; see Supplementary Model and Supplementary Table 3).
Extended Data Fig. 9 Experimental evolution of expanding bacterial populations in various agar concentrations.
a, Three lineages of evolution trajectories in semisolid TB plates with different agar concentrations. b, c, Illustration of how stability distance is determined for each ancestral expansion speed based on the data shown in a. Using linear regression of the mean evolution trajectories obtained (here for 0.25% agar, ancestral ES = 6 mm h^{−1}, n = 3 biological replicates as in a), we obtained a slope for each series, A–E, as in b. Then we plotted the slope obtained in b against the corresponding selection position to obtain another line (c). The stability distance is estimated as the x intercept of the line. d, Batch culture growth curve of ancestor cells grown in TB (purple) or CAA (brown) medium; the growth rate in CAA is almost 50% slower. Two replicates showed similar results. e, The evolution experiment in CAA medium for five agar densities (same as those shown for TB in a). Three independent lineages were run in parallel as for TB.
Extended Data Fig. 10 Stable expansion speeds in spatial habitats of different sizes.
a, Illustration of a recurring ecological scenario. A mixed community of species with different motility characteristics expands into an enclosed habitat of a certain size. All cells within the habitat have an equal chance of being passaged on to occupy a new habitat of the same size. The competitive expansion model (Extended Data Fig. 6a) was extended to model the coexpansion of five distinct species, chemotaxing on the same attractant gradient (Supplementary Analysis 5). Each species has a distinct expansion speed when expanding alone (between 2 mm h^{−1} and 8 mm h^{−1}), due solely to different motility characteristics that are modelled by different chemotactic coefficients c_{i} in the range {34…274} μm^{2} s^{−1}. A closed (zeroflux) boundary condition is applied so that when cells reach the edge of the habitat, they cannot propagate forward but they can propagate backward by diffusion (backward propagation via chemotaxis along reversed chemoattractant gradients is possible in principle but very limited, because the chemoattractant is mostly consumed when the population expands forward throughout the habitat). b, Average expansion speed of the population changes over different cycles of the simulated expansion–selection process. The average expansion speed increases as the selection proceeds for the largest habitat size (25 mm, pink) and decreases for the smallest habitat size (10 mm, green). The relative abundance of species with different singlestrain expansion speeds is changed after each round of selection. The result after the fifth round is shown on the right for the three habitat sizes shown. The faster species is enriched in the largest habitat (top) and the slower species is enriched in the smallest habitat (bottom). The slower species take more cycles to emerge as the winner since they occupy the habitat interior and receive lower weights owing to twodimensional geometry. c, d, To understand why the faster species were selected against in the smaller habitat, we plotted the radial density function ρ_{i}(r,t) for each species i at t = 10 h and 24 h after expansion during the first selection cycle (before any selection took place). The density profiles show that after the faster species reached the edge of the habitat, they moved backward towards the interior. However, the speed of this backward movement was limited compared to the outward movement. This is attributed to the fact that the chemoattractant behind the front is depleted, so the backward movement cannot rely on chemotaxis but is a result of cell diffusion (via Komolgorov–Fisher dynamics). In d, the cyan and blue lines (expansion speeds 7 and 7.8 mm h^{−1}, respectively) almost overlap and are hard to distinguish visually.
Supplementary information
Supplementary Information
This Supplementary Information contains Supplementary Tables 2–4 and two Supplementary Texts (Supplementary Model and Supplementary Analysis) introducing the competition model and presenting an analysis of the competition dynamics.
SNP mutations in the 50th cycle population genome
Supplementary Table 1 . This Supplementary Table lists the SNP mutations which emerged after 50 cycles of expansion dynamics.
Simulation Parameters
Supplementary Data. Parameter sets to be read in by the custommade Python code to numerically solve the GrowthExpansion model. Simulation code available via GitHub at https://github.com/jonascremer/chemotaxis_simulation. Provided parameter sets cover all simulation variants used in the manuscript, see Supplementary Information, Model 1.5 for details.
Expansion dynamics of the ancestor strain in TB
Video 1 . Growth and expansion of fluorescently labelled ancestor cells (GFP) in TB softagar. The positions A,B,C,D, and E are marked as coloured vertical lines. For the corresponding evolution series, cells were picked at these positions after 24h and subsequently transferred to fresh agar plates. Celldensity was determined using timelapse confocal microscopy. See Supplementary Table 2 for strain details. The experiment was independently repeated 3 times with similar results.
Competition dynamics and temporal shift of crossover position
Video 2 . Simulations showing competition of two strains with singlestrain expansion speeds given by \({u}_{{\rm{a}}{\rm{n}}{\rm{c}}}=6\,{\rm{m}}{\rm{m}}/{\rm{h}}\) (purple) and \({u}_{{\rm{m}}{\rm{u}}{\rm{t}}}=7\,{\rm{m}}{\rm{m}}/{\rm{h}}\) (cyan, comparable to evolved D strain). Starting with abundance ratio 1:1, the faster migrating strain is slowly taking over the front. The crossover distance d × (gray vertical line) moves outwards initially, but becomes frozen once the front is sufficiently dominated by the faster expanding strain; see Extended Data Fig. 6c and Supplementary Analysis 1 for description. The attractant profile is shown as magenta dashed line. Parameters provided in Supplementary Table 3, same as those used to generate Extended Data Fig. 6b.
Source data
Rights and permissions
About this article
Cite this article
Liu, W., Cremer, J., Li, D. et al. An evolutionarily stable strategy to colonize spatially extended habitats. Nature 575, 664–668 (2019). https://doi.org/10.1038/s415860191734x
Received:
Accepted:
Published:
Issue Date:
Further reading

Modeling hostassociating microbes under selection
The ISME Journal (2021)

Dissection of the mutation accumulation process during bacterial range expansions
BMC Genomics (2020)

To boldly go where no cell has gone before
Nature Reviews Microbiology (2020)

Chemotaxis as a navigation strategy to boost range expansion
Nature (2019)

A rule from bacteria to balance growth and expansion
Nature (2019)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.