Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# An evolutionarily stable strategy to colonize spatially extended habitats

## Abstract

The ability of a species to colonize newly available habitats is crucial to its overall fitness1,2,3. In general, motility and fast expansion are expected to be beneficial for colonization and hence for the fitness of an organism4,5,6,7. Here we apply an evolution protocol to investigate phenotypical requirements for colonizing habitats of different sizes during range expansion by chemotaxing bacteria8. 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 stability-to-invasion as a powerful principle for the selection of phenotypes in complex ecological processes.

This is a preview of subscription content

## Access options

\$32.00

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 custom-made code and the parameter sets provided.

## Code availability

Custom-made simulation code is available via GitHub at https://github.com/jonascremer/chemotaxis_simulation.

## References

1. 1.

Hanski, I. Metapopulation dynamics. Nature 396, 41–49 (1998).

2. 2.

Skellam, J. G. Random dispersal in theoretical populations. Biometrika 38, 196–218 (1951).

3. 3.

Andow, D. A., Kareiva, P. M., Levin, S. A. & Okubo, A. Spread of invading organisms. Landsc. Ecol. 4, 177–188 (1990).

4. 4.

Yi, X. & Dean, A. M. Phenotypic plasticity as an adaptation to a functional trade-off. eLife 5, e19307 (2016).

5. 5.

Fraebel, D. T. et al. Environment determines evolutionary trajectory in a constrained phenotypic space. eLife 6, e24669 (2017).

6. 6.

Ni, B. et al. Evolutionary remodeling of bacterial motility checkpoint control. Cell Rep. 18, 866–877 (2017).

7. 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. 8.

Adler, J. Chemotaxis in bacteria. Science 153, 708–716 (1966).

9. 9.

Levin, S. A. The problem of pattern and scale in ecology: the Robert H. MacArthur award lecture. Ecology 73, 1943–1967 (1992).

10. 10.

Hastings, A. et al. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett. 8, 91–101 (2005).

11. 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. 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. 13.

Cao, Y. et al. Collective space-sensing coordinates pattern scaling in engineered bacteria. Cell 165, 620–630 (2016).

14. 14.

Cremer, J. et al. Chemotaxis as navigation strategy to boost range expansion. Nature https://doi.org/10.1038/s41586-019-1733-y (2019).

15. 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. 16.

Bosshard, L. et al. Accumulation of deleterious mutations during bacterial range expansions. Genetics 207, 669–684 (2017).

17. 17.

Wolfe, A. J. & Berg, H. C. Migration of bacteria in semisolid agar. Proc. Natl Acad. Sci. USA 86, 6973–6977 (1989).

18. 18.

Deforet, M., Carmona-Fontaine, C., Korolev, K. S. & Xavier, J. B. Evolution at the edge of expanding populations. Am. Nat. 194, 291–305 (2019).

19. 19.

Lenski, R. in Microbial Ecology: Principles, Applications and Methods (eds Levin, M. et al.) 183–198 (McGraw-Hill, 1992).

20. 20.

Smith, J. M. Evolution and the Theory of Games (Cambridge Univ. Press, 1982).

21. 21.

Levin, B. R. Frequency-dependent selection in bacterial populations. Phil. Trans. R. Soc. Lond. B 319, 459–472 (1988).

22. 22.

Alon, U., Surette, M. G., Barkai, N. & Leibler, S. Robustness in bacterial chemotaxis. Nature 397, 168–171 (1999).

23. 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. 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. 25.

Park, H., Guet, C. C., Emonet, T. & Cluzel, P. Fine-tuning of chemotactic response in E. coli determined by high-throughput capillary assay. Curr. Microbiol. 62, 764–769 (2011).

26. 26.

Si, G., Wu, T., Ouyang, Q. & Tu, Y. Pathway-based mean-field model for Escherichia coli chemotaxis. Phys. Rev. Lett. 109, 048101 (2012).

27. 27.

Tu, Y. Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation. Annu. Rev. Biophys. 42, 337–359 (2013).

28. 28.

Liu, C. et al. Sequential establishment of stripe patterns in an expanding cell population. Science 334, 238–241 (2011).

29. 29.

Merrell, D. The Adaptive Seascape: The Mechanism of Evolution (Univ. Minnesota Press, 1994).

30. 30.

Mustonen, V. & Lässig, M. From fitness landscapes to seascapes: non-equilibrium dynamics of selection and adaptation. Trends Genet. 25, 111–119 (2009).

31. 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. 32.

Towbin, B. D. et al. Optimality and sub-optimality in a bacterial growth law. Nat. Commun. 8, 14123 (2017).

33. 33.

Hofbauer, J. & Sigmund, K. Evolutionary Games and Population Dynamics (Cambridge Univ. Press, 1998).

34. 34.

Gore, J., Youk, H. & van Oudenaarden, A. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256 (2009).

35. 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. 36.

Datta, S., Costantino, N. & Court, D. L. A set of recombineering plasmids for Gram-negative bacteria. Gene 379, 109–115 (2006).

37. 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. 38.

Fu, J. L., Kanno, T., Liang, S.-C., Matzke, A. J. M. & Matzke, M. GFP loss-of-function mutations in Arabidopsis thaliana. G3 5, 1849–1855 (2015).

39. 39.

Waite, A. J. et al. Non-genetic diversity modulates population performance. Mol. Syst. Biol. 12, 895 (2016).

40. 40.

Edelstein, A. D. et al. Advanced methods of microscope control using μManager software. J. Biol. Methods 1, e10 (2014).

41. 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. 42.

Deatherage, D. E. & Barrick, J. E. Identification of mutations in laboratory-evolved microbes from next-generation sequencing data using breseq. Methods Mol. Biol. 1151, 165–188 (2014).

43. 43.

Guyer, J. E., Wheeler, D. & Warren, J. A. FiPy: partial differential equations with Python. Comput. Sci. Eng. 11, 6–15 (2009).

44. 44.

Fisher, R. The wave of advance of advantageous genes. Ann. Eugen. 7, 355–369 (1937).

45. 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. 46.

Korolev, K. S. Evolution arrests invasions of cooperative populations. Phys. Rev. Lett. 115, 208104 (2015).

47. 47.

Yang, F., Moss, L. G. & Phillips, G. N. Jr. The molecular structure of green fluorescent protein. Nat. Biotechnol. 14, 1246–1251 (1996).

48. 48.

Keller, E. F. & Segel, L. A. Model for chemotaxis. J. Theor. Biol. 30, 225–234 (1971).

49. 49.

Vaknin, A. & Berg, H. C. Physical responses of bacterial chemoreceptors. J. Mol. Biol. 366, 1416–1423 (2007).

50. 50.

Taylor, J. R. & Stocker, R. Trade-offs of chemotactic foraging in turbulent water. Science 338, 675–679 (2012).

51. 51.

Fu, X. et al. Spatial self-organization 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 (KFZD-SW-216) 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

Authors

### 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

Correspondence to Terence Hwa or Chenli Liu.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

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 bacteria14 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 (AncG) at various times after inoculation at the centre of a semi-solid 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 previously14. 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 XA, XB, XC, XD and XE, 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 mid-exponential phase (OD600 = 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 single-cell 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.

ad, 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 glycerol-containing 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 (ac). In the absence of a chemoattractant, cells follow simple Fisher–Kolmogorov dynamics14,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 density-dependent growth effects or strong tradeoffs between swimming and cell growth exist18,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 non-fluorescent 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 (PZA31-Ptet-M2-GFP or PZA31-Ptet-M2-NFP, respectively) were transformed into the six evolved strains A, B, B′, C, D, and E, to form AG, AN 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 × 1010 cells per ml. Subsequently, serial dilutions were carried out. For each cell density, cells were vigorously mixed with pre-warmed 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 fluorescence-activated cell sorting (FACS) and the fluorescence intensity of the agar plate was detected using a fluorescence microscope. d, The three-dimensional structure of GFP with the predicted position of the loss-of-function mutation introduced, Y66C. The shown crystal structure is based on the NFP protein sequence aligned to PDB ID: 1GFL47. The rectangle shows the mutation position. e, Ancestor cells harbouring constitutive GFP or NFP expression (AncG and AncN cells, respectively) as viewed by fluorescence microscopy; the images verify the loss-of-function mutation of GFP. f, The relative fitness W of AncG and AncN 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 cheZ-titratable strains WL1 (g) and WL2 (h). The expression of cheZ is under the control of a Ptet-tetR (WL1) or Plac-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. il, Growth rates (i, k) and expansion speeds (j, l) of the aTc-titratable cheZ strain WL1 (i, j) and its derivatives expressing GFP (WL7) or NFP (WL8), and the IPTG-titratable cheZ strain WL2 (k, l) and its derivatives expressing GFP (WL9) or NFP (WL10). m, Circles indicate crossover distances between fluorescent derivatives of two cheZ-titratable 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, gl, 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 two-strain 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 × 107 cells per ml. c, Relative fluorescence intensities obtained at various times for the fluorescent mutant strains BG (orange) and DG (cyan) and the ancestor AncG (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 two-strain 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 head-to-head 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 DG and AncN, competition between DN and AncG, and merged photograph in pseudocolour (DG, cyan; AncG, 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 AncN while the evolved strain with NFP was mixed with AncG in a 1:1 ratio. Subsequently, 2 μl of the two types of combined mixture was inoculated onto the centre of pre-warmed semi-solid 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 WB 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 BG and AncN, competition between BN and AncG, and averaged curve of both (the relative fitness WB). 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 24-h competition. Beyond the initial period (~6 h), the crossover distance could be clearly defined (vertical dashed line) and was practically time-independent, with the slower strain gaining advantage in the interior and the faster strain gaining advantage in the exterior. h, j, Relative fitness values WD (h) and WB (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 (Wi = 1). Thus, the crossover distance was fixed shortly after the initial period. km, Fluorescence intensity profiles (k, l), and relative fitness W (m) of representative two-strain competitions between the ancestor (black solid line) and evolved isolates. The data were taken 12 h (k, m) or 24 h (l) after co-inoculation 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 al 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 model14 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 Segel48 (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 scale-free manner only within a limited range between KI and KA, the lower and upper cutoffs describing the molecular limitations of attractant sensing26,49 as specified in equation (S2). D0 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 non-chemotactic 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 previously14. 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 < KI). 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 multi-strain 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 back-diffusion 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 back-diffusion. 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 front14.

### Extended Data Fig. 6 Model of two-strain competitive expansion dynamics.

a, The GE model of bacterial range expansion14 (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 Di 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 vi on the attractant profile is given by equation (S7), which describes a range of proportional sensing v xa/a for KI < a < KA. The magnitude of the chemotactic response is parametrized by the chemotactic coefficient ci for strain i. Finally, the dynamics of the nutrient and attractant are described by equations (S8) and (S9), with D0 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 uanc = 6 mm h−1) and the other strain resembling a mutant, with expansion speed umut 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 uanc versus umut = 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 (c1, c2). 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 growth-rate differences are minimal (for example, <5% as we observed experimentally for TB; Extended Data Fig. 1e) but become substantial when growth-rate differences become large (Supplementary Analysis 3). g, The stability distance d* as in d but for an alternative model formulation of chemotactic movement50 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 growth-rate-dependent slopes was still observed. In bg, 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 Two-strain competition and crossover dynamics.

The results in Fig. 2 show that competition between two strains with the same growth rate and different single-strain 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 (cslow = 652 μm2 s−1 and cfast = 727 μm2 s−1, with corresponding single-strain expansion speeds uslow = 4.95 mm h−1 and ufast = 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 single-strain 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 time-invariant. 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 co-migrating population. This is because the two strains chase after the same attractant profile, which can recede only at a single speed51. 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 time-invariant crossover distance d×(ufastuslow) for two strains with single-strain expansion speeds ufast and uslow, and shows how this crossover distance leads to an evolutionarily stable selection distance d* defined in the limit uslowufast (Figs. 3, 4). The two density profiles shown in a, ρfast(xt) and ρslow(xt) in the ‘laboratory frame’ coordinate x, are shown in b in the frame moving with speed ufast. In this moving frame where the spatial coordinate is y = x − ufastt, 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 single-strain dynamics in the reference frame moving at speed ufast. 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 single-strain 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× = ufastt× (obtained from y(tx) = 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 γ (ufast − uslow), 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× = ufastt× depends on ufast and uslow 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×(uslowu). Supplementary Analysis 2 shows that t×(uslowu) 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 position-dependent fitness landscapes for the competition between two strains (strains 1 and 2, for example, an ancestor and a mutant), with expansion speeds u1 and u2, 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/W0. In a, b, green indicates a gain in the fitness of the mutant (W/W0 > 1) and blue indicates a loss (W/W0 < 1). Top, middle, and bottom rows refer to results at distances X = 20, 15 and 10 mm, respectively. a, Landscape W/W0 for equal abundance (50:50) of the two strains at equal inoculation (W0 = 1). Cyan and dashed represent W/W0 = 1 and u1 = u2, respectively. For W0 = 1, the cyan line is the crossover distance: X = d×(uu′) (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 (W0 ≡ ρ2(X, 0): ρ1(X, 0) = 5:95). The solid brown line still indicates the phase boundary W/W0 = 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 frequency-dependent 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 W0 in a complex way, in the vicinity of the diagonal (u1 ≈ u2), W is independent of W0 so that the illustration in Fig. 3a of the stability of a strain of speed u at distance d×(uu) is verified explicitly in c, d. To vary expansion speeds in each panel, the chemotactic coefficients ci were varied within the range 180–1,500 μm2 s−1. Diffusion coefficients were changed accordingly (Di = ci/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 semi-solid 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 co-expansion 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 ci in the range {34…274} μm2 s−1. A closed (zero-flux) 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 single-strain 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 two-dimensional 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 custom-made Python code to numerically solve the Growth-Expansion 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 soft-agar. 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. Cell-density was determined using time-lapse 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 cross-over position

Video 2 . Simulations showing competition of two strains with single-strain 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 cross-over 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.

## Rights and permissions

Reprints and Permissions

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/s41586-019-1734-x

• Accepted:

• Published:

• Issue Date:

• DOI: https://doi.org/10.1038/s41586-019-1734-x

• ### Modeling host-associating microbes under selection

• Florence Bansept
• Nancy Obeng
• Arne Traulsen

The ISME Journal (2021)

• ### Dissection of the mutation accumulation process during bacterial range expansions

• Lars Bosshard
• Stephan Peischl
• Laurent Excoffier

BMC Genomics (2020)

• ### To boldly go where no cell has gone before

• Ashley York

Nature Reviews Microbiology (2020)

• ### Chemotaxis as a navigation strategy to boost range expansion

• Jonas Cremer
• Tomoya Honda
• Terence Hwa

Nature (2019)

• ### A rule from bacteria to balance growth and expansion

• Henry Mattingly
• Thierry Emonet

Nature (2019)