We experimentally and numerically investigate the evolutionary dynamics of four competing strains of E. coli with differing expansion velocities in radially expanding colonies. We compare experimental measurements of the average fraction, correlation functions between strains, and the relative rates of genetic domain wall annihilations and coalescences to simulations modeling the population as a one-dimensional ring of annihilating and coalescing random walkers with deterministic biases due to selection. The simulations reveal that the evolutionary dynamics can be collapsed onto master curves governed by three essential parameters: (1) an expansion length beyond which selection dominates over genetic drift; (2) a characteristic angular correlation describing the size of genetic domains; and (3) a dimensionless constant quantifying the interplay between a colony’s curvature at the frontier and its selection length scale. We measure these parameters with a new technique that precisely measures small selective differences between spatially competing strains and show that our simulations accurately predict the dynamics without additional fitting. Our results suggest that the random walk model can act as a useful predictive tool for describing the evolutionary dynamics of range expansions composed of an arbitrary number of genotypes with different fitnesses.
A competition between stochastic and deterministic effects underlies evolution. In a well-mixed system such as a shaken culture of the yeast microorganism Saccharomyces cerevisiae, stochastic competition between individuals, mutations, and selection dictate the dynamics of the population [1]. In spatially structured environments, active or passive dispersal of individuals also plays an important role. The local “well-mixed” dynamics must be coupled to the motion of individuals, leading to strikingly different evolutionary dynamics, even in the absence of selection [2–7]. A model laboratory system that can be used to explore the coupling between local “well-mixed” effects and spatial deterministic and stochastic dynamics is a microbial range expansion [8], in which a population expands into an unoccupied region of a hard agar Petri dish. Non-motile microbes expand outwards from their initial position due to a combination of growth coupled with random pushing by neighboring cells and leave behind a record of their genetic competition as they cannot move and cease reproducing once the population becomes too dense [8]. A frozen genetic pattern of four competing strains of E. coli marked by different fluorescent colors can be seen in Fig 1. Spatial structure is present in the frozen genetic patterns because the microbes at the expanding frontier produce daughter cells of the same color that migrate only a small fraction of the front circumference within a generation. Hallatschek et al. [8] identified the key role of genetic drift in producing these sectored patterns; the small population size at the front of an expanding population [9,10] enhances number fluctuations (i. e. genetic drift), eventually leading to the local fixation of one strain past a critical expansion radius R0. The decrease in genetic diversity as the small number of individuals at the frontier expands is referred to as the “Founder effect” [11]. Outside of the laboratory, range expansions occur naturally during the spread of invasive species such as the bank vole in Ireland [12] or the cane toad in Australia [13], and played a role in the evolutionary history of humans when migrating out of Africa [14]. In these natural expansions, populations may have many competing genotypes, or alleles, each instilling a different fitness. Even if a population is originally clonal, mutations may create new alleles that compete with one another to proliferate, a phenomenon known as clonal interference [15]. An allele’s fitness is often determined by its corresponding expansion velocity. Faster expanding individuals will colonize more territory and will block slower strains from expanding, resulting in the increased abundance of ‘faster’ alleles at the frontier [13,16,17]. If the curvature of a microbial colony can be neglected and its front is sufficiently smooth, it has been shown both theoretically and experimentally that the domain wall of a faster expanding strain will displace a slower expanding strain at a constant rate per length expanded after an initial transient, resulting in a characteristic triangular shape [17] as shown on the right side of Fig 1. If the curvature of the expansion is not negligible, the sector boundaries will trace logarithmic spirals [17]. Even in the most simple scenario when de-novo mutations and mutualistic or antagonistic interactions are ignored, the dynamics of many competing alleles with varying fitnesses at the front of a range expansion have neither been quantified theoretically nor explored in laboratory experiments. Prior laboratory experiments focused on the dynamics of a single sector of a more fit strain (representing a competing alelle) of yeast sweeping through a less fit strain [17] in regimes where stochastic wandering of genetic boundaries was not expected to be important. Recent experimental work studied how fast a single more fit strain swept through a less fit strain in a range expansion and compared the dynamics to the same strains in a well mixed test tube [9]. In this paper, we experimentally and numerically investigate the dynamics of four competing strains (alleles) of E. coli with varying selective advantages initially distributed randomly at the front of a radial range expansion. The eCFP (blue) and eYFP-labeled (yellow) strains expanded the fastest, followed by the non-fluorescent (black) strain, and finally the mCherry-labeled (red) strain. The differences in expansion speeds are reflected in Fig 1 as follows: the yellow/blue bulges at the front of the expansion are larger than the black bulges which are larger than the red bulges. The significant random undulations at the frontier, however, significantly mask the selection-induced bulges. As is evident from Fig 1, the size and location of a monoclonal sector can be described by the locations of its boundaries. When two boundaries collide, they either annihilate if the neighbors to the left and right of the collision are the same or coalesce if the neighbors are different, as illustrated by the A and C respectively on the left side of Fig 1. We therefore describe our expansions as a one-dimensional line of annihilating and coalescing random walkers, a description that has been used extensively in previous work (see Ref. [2] for a review). To account for the radial geometry of our colonies, we allow the frontier to inflate, corresponding to the increasing perimeter of the colony as its radius increases. Past the radius R0 where genetic domains originally form, we describe the random motion of genetic domains by a diffusion constant per length expanded Dw (see Fig 1) [18]. If dx characterizes the displacement of a domain wall perpendicular to the expansion direction and dL is the distance the colony has expanded (the radius that the colony has grown) as illustrated on the right side of Fig 1, where we neglect the circumferential curvature in this small region, we define the diffusion constant per length expanded as 2Dw = dVar (x) /dL where Var (x) ≡ 〈x2〉 − 〈x〉2 is the variance and the brackets indicate an average over many domain walls. Note that Dw has dimensions of length. Similarly, differences in expansion velocities between neighboring strains will lead to the deterministic displacement of domain walls per length expanded as the faster expanding strain will reach the contested point on the front before a slower growing strain as mentioned above [17]; we characterize this deterministic motion by a dimensionless “wall velocity, ” [18] v w i j = d 〈 x 〉 / d L, where i is the strain to the left of the domain wall and j is the strain to the right. Note that v w i j = − v w j i. The dynamics of an arbitrary number of neutral competing strains in an expansion (i. e. v w i j = 0 for all domain walls) is well understood as the dynamics can be described as a one-dimensional system of annihilating and coalescing random walkers [19–21] which is equivalent to a one-dimensional q-state Potts model [22,23] governed by zero-temperature Glauber dynamics [24] or a q-opinion Voter model [25,26]. Many theoretical predictions and analyses of this system exist; of particular relevance to this paper are the relative annihilation and coalescence rates per collision as q is varied [27–29] and the calculation of spatial correlation functions [28]. To map standard linear results onto an inflating ring (i. e. including R0 in the models), one can use a conformal time transformation [30–32]. Fewer results are available in the presence of selection, i. e. when domain walls have deterministic biases (nonzero v w i j) [33]. Analytical results are rare because the moment hierarchy of this model does not close [2] as discussed in S1 Appendix. In this paper, we measure and predict three quantities relevant to the evolutionary dynamics of our four competing strains of E. coli in radial range expansions: the average fraction of each of our four strains, the two-point correlation functions between our strains, and the relative annihilation and coalescence probabilities per domain wall collision (see Fig 1), a quantity that has received theoretical attention [27–29] but has neither been explored experimentally nor investigated in the presence of selection. We measure these three quantities using an image analysis toolkit (available on GitHub, complete with examples of how to use it [34]) that extends experimental techniques for two-color (two-allele) range expansions [8,9, 17,18,35,36] to an arbitrary number of competing strains. We next use an efficient radial simulation (also on GitHub [34]) of annihilating and coalescing random walkers with deterministic wall velocities to determine what sets the scale of the dynamics and to synthesize our experimental and theoretical results. We show that three key combinations of R0, Dw and v w i j control the dynamics of our four strains. We conclude with suggestions for future studies. The details of our experimental, theoretical, and simulation methods are given in the last section.
We begin by reporting our measurements of the average fraction of each strain, the two-point correlation functions between strains, and the relative rates of annihilations and coalescences as a function of length expanded for our four competing strains of E. coli. As discussed in the Materials and Methods, we found that our eCFP and eYFP strains had the fastest expansion velocities followed by the black strain and finally the mCherry strain (see Table 1). We expected that our experimental measurements would reflect this hierarchy of speeds; faster expanding strains should have a larger fitness than slower expanding ones. To illustrate the presence of selection, we used neutral theory (discussed in detail in S1 Appendix) as a null expectation; selection caused deviations from the neutral predictions. To calibrate neutral theory to our experiments we fit R0 and Dw, two model parameters illustrated in Fig 1, following the procedures discussed in the Materials and Methods. The fit values of R0 and Dw can be seen in Table 2. In later sections, we show how to predict the average fraction, two-point correlation functions, and relative rates of annihilation and coalescences using our random-walk model and simulation. In this section, we introduce three key combinations of our random walk model’s input parameters R0, Dw, and v w i j (see Fig 1) that control the evolutionary dynamics of our four competing E. coli strains. Using simulation, we show that we can utilize these key combinations to collapse the simulated evolutionary dynamics (focusing on the experimental quantities we measured above: the average fraction, two-point correlation function, and annihilation asymmetry) of an arbitrary number of competing strains in a range expansion. A major goal of this paper is to test if the annihilating and coalescing random-walk model can predict the experimental evolutionary dynamics of our four competing strains (alleles) with different fitnesses (radial expansion velocities). To the best of our knowledge, analytical results for the random-walk model are unavailable (as discussed in S1 Appendix); we consequently used our simulations to predict the dynamics. In this section we quantify the three key parameter combinations for our experimental expansions and then use them to predict the evolutionary dynamics of all four of our competing E. coli strains in an independent experiment. In the last section, we found that our simulation dynamics could be collapsed onto master curves for a fixed set of κ i j = R 0 / L s i j by rescaling the length expanded L by any single L s i j = D w / (v w i j) 2 and by rescaling ϕ by ϕ c = 8 D w / R 0. These simulated master curves were invariant to the alteration of simulation parameters provided that the set of κij remained the same. This insight allowed us to develop a novel method of characterizing the experimental dynamics. Namely, we could experimentally determine L s i j, κij, and ϕc, collapse the experimental data the same way as the simulations (i. e. L / L s i j, ϕ/ϕc), and compare the two to predict the dynamics of many competing alleles in a range expansion. As discussed below, this technique ultimately allowed for accurate predictions of the evolutionary dynamics of the four competing strains and, surprisingly, allowed us to make much more precise measurements of selective differences between strains. As mentioned above in the Experimental Results section, using the procedures outlined in the Materials and Methods, we had previously determined R0 = 3. 50 ± 0. 05 mm and Dw = 0. 100 ± 0. 005 mm (Table 2). In order to fit L s i j = D w / (v w i j) 2 and κ i j = R 0 / L s i j, however, we needed to measure v w i j. By tracking the growth of a more fit sector sweeping through a less fit strain (see the Materials and methods), we found that each strain swept through mCherry with a wall velocity of v w i R = 0. 06 ± 0. 02 (as seen in Table 2); we could not detect the wall velocity of the eYFP and eCFP sweeping through the black strain. In principle, the measured values of R0, Dw, and v w i j should have allowed us to totally calibrate the three key parameter combinations. For example, ϕ c = 8 D w / R 0 = 0. 48 ± 0. 01. The value of κ i j = R 0 / L s i j followed from the measurement of L s i j using the known value of R0. Unfortunately, the final parameter L s i j = D w / (v w i j) 2 was more difficult to calibrate. Using v w i R = 0. 06 ± 0. 02, we found that L s i R = 30 ± 20 m m; the error on this value was too large for it to be predictive in our simulations. Furthermore, as we were unable to accurately measure the wall velocity of the eCFP/eYFP strains sweeping through the black strain, we could not calculate the corresponding selection length scale. We therefore needed a new technique to determine L s i j. As our eCFP and eYFP strains were neutral within error, we treated our system as composed of one neutral (N) eCFP/eYFP strain, a red (R) mCherry strain, and a black (B) strain (q = 3 colors). As the eCFP/eYFP expanded faster than the black followed by the mCherry strain, we needed to determine the values of L s N R, L s N B, and L s B R. To fit L s i j more precisely than that from our direct measurement of wall velocity, we competed pairwise combinations of strains in range expansions (i. e. the eCFP/eYFP strain vs. mCherry) and calculated the two-point correlation functions Fij (L, ϕ) at the maximum length expanded of L = 6. 5 mm. As there were only two competing strains, there was only one Ls. To fit the value of Ls, we began by rescaling the experimental length expanded L by Ls and ϕ by ϕc (Table 2) and calculated the resulting κ = R 0 / L s. Note that Ls simultaneously rescales the length expanded L by Ls and sets the value of κ = R 0 / L s, changing the shape of the collapsed correlation function. We then ran a simulation at the set value of κ (the chosen simulated values of Ls and ϕc did not matter due to the collapse) and then compared the collapsed experimental dynamics to our simulation. Fig 7 illustrates the fitting procedure by displaying the experimentally rescaled two-point correlation function FNR (the solid red line) at a length expanded of L = 6. 5 mm between our eCFP/eYFP strain (N) and our mCherry strain (R) (inoculated at fractions of 2/3 and 1/3 respectively) and simulated universal correlation functions corresponding to different values of Ls (dashed lines). To determine the best-fitting value of Ls, we calculated the sum of the squared displacements weighted by the inverse of the experimental standard error squared between experiment and simulation. The best-fitting Ls was determined by finding the value which minimized the weighted sum of squares. To estimate the error in our fit, we assigned each potential value of Ls a probability proportional to the inverse of the weighted sum of squares, normalized the probability distribution, and set the error in our fit of Ls to the confidence intervals of the probability distribution. Our fit values of L s i j and κij using this technique are listed in Table 3; the values of κij are also plotted in Fig 5. Although this technique was about a factor of 5 more precise than using the measured wall velocities v w i R to determine L s i j, the upper bounds of the 95% confidence intervals were still very large as seen in Table 3; the potential values of Ls had a very large tail. To test that the resulting Ls and κ could accurately predict the experimental dynamics at all L and not just the L where the correlation functions were fit, we plotted the experimental average fraction and correlation functions (solid lines, Fig 8) as we varied L and compared their values to those predicted by simulation (dashed lines, Fig 8). Fig 8 uses the same set of experimental data as that from Fig 7. The simulation using the fit parameters always closely tracked the experimental values at all L, suggesting that our fitting technique was robust and could be used to describe the dynamics of our strains. Having determined L s i j and κij from pairwise competitions between strains, we tested if we could predict the average fraction, the two-point correlation functions, and the annihilation asymmetry when the four E. coli strains were grown together (treating the eYFP and eCFP strains as neutral, so q = 3) in an independent experiment. We inoculated the four strains in equal proportions. Fig 9 shows experimental measurements of the average fractions and two-point correlation functions (solid lines) together with simulated predictions (dashed lines) that used the independently fit values; no additional fitting parameters were used. The predicted average fractions and correlation functions closely tracked the dynamics for L ≳ 3 mm. We attribute the deviations for L ≲ 3 mm to image analysis artifacts resulting from the presence of the black strain (see the Image Analysis section in the Materials and methods). At the largest length expanded of L = 6. 5 mm where artifacts were minimal, the experiments matched the predictions within error. All average correlation functions at this length expanded were successfully predicted by the simulations; we only display FNR for simplicity. In addition to predicting the average fractions and correlation functions, the simulation with our fit L s i j and κij predicted that the annihilation asymmetry would deviate only slightly from neutrality (at most a change of 0. 1) over the length expanded by our colonies in every experiment, agreeing with our findings (Fig 4). This can be readily observed in Fig 6 which displays a simulation of two neutral strains (our eCFP and eYFP strains) and a less fit strain (our mCherry strain) inoculated in equal fractions. If we rescale the maximum distance expanded by our colonies of Lmax = 6. 5 mm by the smallest selection length (this results in the largest possible change of ΔP) of L s N R = 13 mm, Lmax/Ls ∼ 0. 5 and ΔP increases from 0 (neutrality) to at most 0. 1. This small deviation from neutrality is within the uncertainty of our experimental measurement of ΔP. Evidently, certain quantities, like the average fraction and correlation functions, show signs of selection before others (in this case, the annihilation asymmetry). The quantitative agreement between our model and our experiments suggests that the one-dimensional annihilating-coalescing random walk model can indeed be used to predict the dynamics of many competing strains with different fitnesses in a range expansion.
We investigated the evolutionary dynamics of four competing strains of E. coli in radial range expansions with differing selective advantages. We measured the average fraction Fi of each strain, the two-point correlation functions Fij between strains, and the annihilation asymmetry ΔP with our image analysis toolkit [34]. Our simulations, which model the expansions as a one-dimensional line of random walkers subject to deterministic drift on an inflating ring, showed us that these three quantities could be collapsed onto universal curves for a fixed set κ i j = R 0 / L s i j when the length expanded by the colony L was rescaled by any of the lengths L s i j = D w / (v w i j) 2 and the angular distance between strains ϕ was rescaled by ϕ c = 8 D w / R 0. To test if the random walk model could predict experimental dynamics, we independently calculated experimental values of L s i j, κij, and ϕc and compared the dynamics between the two. The simulations accurately predicted the dynamics of the average fraction, correlation functions, and annihilation asymmetry when all four of our strains were present with no additional fitting parameters. The annihilation asymmetry ΔP is a quantity unique to range expansions with three or more strains and, to the best of our knowledge, has not been studied previously. Our results illustrate the importance of considering domain wall annihilation and coalescence when more than two strains are present and suggest that the annihilating-coalescing random-walk model can act as a useful predictive tool when describing the evolutionary dynamics of range expansions composed of an arbitrary number of competing alleles with different fitnesses. Along the way, we introduced a new technique that compared universal simulated correlation functions to experimental correlations to fit L s i j. The resulting values of L s i j were about a factor of 5 more precise than directly evaluating L s i j = D w / (v w i j) 2 with the wall velocities extracted from the growth of sectors. Given our fit L s i j, we evaluated v w i j using L s i j = D w / (v w i j) 2 and the known Dw; we compare these values to those extracted from single sectors in Table 4. The wall velocities from both measurements agreed within error, but the wall velocities obtained from our correlation method were at least a factor of two more precise than tracking single sectors. Although the correlation technique dramatically increased the precision in our evaluation of L s i j, the resulting precision increase for v w i j was less pronounced as v w i j ∝ 1 / L s i j. Nevertheless, it is clear that the correlation technique can be used to precisely extract small differences in fitness between spatially competing strains. Our work illustrates that the annihilating-coalescing random walk model can predict the experimental dynamics of an arbitrary number of competing alleles with different fitnesses in microbial range expansions. It is possible that this model could predict the dynamics of range expansions occurring outside of the laboratory, especially if the expanding organisms’ underlying motion did not completely smear out the population’s spatial structure; the organismal motion could potentially be accounted for by increasing the domain wall diffusion coefficient Dw. To predict the dynamics of expansions, however, the annihilating-coalescing walk model relies on a key set of parameters: the set of L s i j, the set of κij, and ϕc. We found that the set of L s i j could not be predicted from the independent radial expansion velocities of our strains; standard techniques [17] using the relative ratio of expansion velocities to predict v w i j, and thus L s i j, yielded inconsistent results (see S2 Appendix where we quantify the discrepancy and postulate why it occurred). As the set of L s i j is so fundamental to the evolutionary dynamics of range expansions, future work should investigate why relative radial expansion velocities could not be used to accurately predict v w i j and thus L s i j and whether this phenomenon is specific to E. coli range expansions or our specific strains. It would also be interesting to incorporate the reported super-diffusive motion of domain walls [8,9] into our simplified simulations and theoretical analysis. The random walk model’s ability to successfully predict the evolutionary dynamics of our experiments suggests that annihilating and coalescing genetic domain walls subject to diffusion and selection-induced displacement provide a useful conceptual framework from which to understand range expansion dynamics.
We used four E. coli strains (labelled BW001, BW002, BW003, and BW012) with a DH5α background and plasmids whose sequences coded for spectrally distinguishable fluorescent proteins. The unique colors were obtained by using the plasmid vector pTrc99a [39] and the open reading frame for the respective fluorescent proteins. Strains BW001, BW002, and BW003 expressed eCFP (cyan/blue), Venus YFP (yellow), and mCherry (red) respectively, and were identical to the E. coli strains eWM282, eWM284, and eWM40 used in Ref. [40]. Note that these three strains were isogenic and differed only by the open reading frames corresponding to their respective fluorescent proteins. The final strain, BW012, was a mutated descendant of strain BW002 (yellow) that fluoresced at a decreased intensity, appearing black, while retaining its ampicillin resistance from the pTrc99a vector. Throughout this work, no additional mutations were introduced or observed. We therefore consider that these four strains correspond to four different alleles. Throughout the paper, we refer to the strains as eCFP, eYFP, mCherry, and black. To prepare saturated cultures, strains were inoculated in 10mL of 2xYT media and were shaken for approximately 16 hours at 37°C. After vortexing each saturated culture and obtaining their concentration via optical density (OD-600) measurements, appropriate volumes (e. g. , 1: 1: 1 mixtures of three strains) were added to an Eppendorf tube with a final volume of 1mL. The Eppendorf tube was then vortexed to uniformly mix the strains. A volume of 2 μL was taken from the vortexed tube and placed on center of a 100 mm diameter Petri dish containing 35 mL of lysogeny broth (LB), ampicillin at a concentration of 100 μg/mL, and 1. 25% w/v bacto-agar. The carrier fluid in the resulting circular drop evaporated within 2-3 minutes, depositing a circular “homeland” of well-mixed bacteria onto the plate. After inoculation, plates were stored for 8 days upside down (to avoid condensation) in a Rubbermaid 7J77 box at 37°C with a beaker filled with water; the water acted as a humidifier and prevented the plates from drying out. The plates were occasionally removed from the box and imaged (at roughly 24 hour intervals) using the brightfield channel to determine the radius of the colony as a function of time. On the eighth day, the plates were imaged in both fluorescent and brightfield channels. The number of replicate plates used are stated next to the respective experimental results. If we noticed that a mutation had occurred during an expansion (mutations usually presented themselves as unexpected large bulges at the front of a colony or as distortions in fluorescent intensity), we discounted the colony. We imaged our range expansions with a Zeiss SteREO Lumar. V12 stereoscope in four channels: eCFP, eYFP, mCherry (fluorescent channels), and brightfield. In order to analyze a colony with a maximum radius of approximately 10 mm using a single image, we stitched four images together with an overlap of 20% using AxioVision 4. 8. 2, the software accompanying the microscope. We blended the overlapping areas of the images to lessen the impact of background inhomogeneities. An example of a stitched image can be seen on the left side of Fig 10. Stitching introduced small artifacts such as vertical lines near the center of our expansions; we verified that these did not affect our results. To extract the local fraction of each strain per pixel, we first created binary masks for each fluorescence channel indicating if the corresponding E. coli strain was present. We utilized the “Enhance Local Contrast” (CLAHE) algorithm [41] in Fiji [42], an open-source image analysis platform, to help correct for inhomogeneities in background illumination. After applying the CLAHE algorithm, a combination of automatic thresholding and manual tracing yielded a binary mask of each channel, an example of which is shown in Fig 10; the image on the left is an overlay of an experimental range expansion’s fluorescent channels and the image on the right is the overlay of the corresponding binary masks. A small amount of manual tracing was required near the edges of our colonies because our fluorescent lamp provided uneven illumination; resulting dark regions could barely be identified above background noise. As we mainly used manual tracing near the edge of the colonies where the monoclonal sectors were well defined, we found that our procedure was very reproducible. To alleviate this problem, future work could utilize brighter strains or a more advanced imaging setup. We mapped the binary images to the local fraction of each E. coli strain in the following way: if N binary masks (corresponding to N colors) were “on” at a pixel, the local fraction of their corresponding channels was assigned to be 1/N. Although this assignment produces inaccuracies (i. e. , if one strain occupied 90% of a pixel and the other occupied 10%, our algorithm would register both as 50%), domain boundaries were the only areas besides the homeland and the early stages of the range expansions where multiple strains were colocalized. The black strain was defined to be present at pixels reached by the range expansion in which no other strains were present. Although this definition introduced errors at radii close to the homeland with significant color overlap, the error became negligible at large radii as quantified in Supplementary S5 Fig. Once we determined the fraction of each strain at each pixel, we were able to extract quantities such as the total fraction of each strain in the colony and spatial correlations between strains at a given expansion radius. The mask in Fig 10 highlights that sector boundaries can be used to determine local strain abundance. Although it is possible to extract the position of every domain wall from each strains’ local fraction, it is challenging to actually track a single wall due to collisions between walls. To address this problem, we created a binary mask of the edges in our images and labelled the edges of each domain. Annihilations and coalescences were counted manually within Fiji [42]; automated measures were not accurate enough. It is worth pointing out that in this paper, we ignore the three-dimensional structure of our colonies and describe them by our two-dimensional images taken with the stereoscope. We justify this approximation because the initial diameter of our colonies is at least a factor of 10 larger than their height (less than 1 mm as judged by a ruler), so they are effectively two-dimensional, and because the strain composition of our colonies does not vary with height inside the colony. We confirmed that strain composition does not vary with height by using a confocal microscope to probe the internal structure and also by taking a pipette tip, scratching it through a sector, growing the cells touched by the tip in overnight culture, and verifying that plated single colonies from the culture were the same color as the sector. We used the average expansion velocity of each strain for radii R > R0 as a proxy for selective advantage, similar to previous work [17,35]. In three independent sets of experiments using different batches of agar plates (the main source of variability in our experiments), we measured the diameter of 12 expansions of each strain approximately every 24 hours following the protocol for range expansions with two or more strains. To account for biological variance, sets of four of the 12 colonies were created from independent single colonies; no statistical difference was seen between biological replicates. The diameters were determined by manually fitting a circle to a brightfield image of the expansion three times and averaging the measured diameters. Fig 11 shows the average radius increasing with time for each strain from one of our experiments. In every experiment, the eCFP and eYFP strains had the fastest expansion velocities (the respective datapoints overlap in Fig 11), followed by the black strain, and then finally the mCherry strain. The expansion velocity slowly decreased as a function of time; we attribute this to nutrient depletion in the plates. The radial expansion velocity of each strain was obtained by using linear regression to fit the radius versus time for radii greater than R0. We calculated the average radial expansion velocity between the three sets of plates and reported its error as the standard error of the mean; see Table 1. Additionally, we quantified the dimensionless selective advantage of each strain relative to the slowest growing mCherry strain following [17] via siR = ui/uR − 1 where the R indicates the mCherry strain (red) in each experiment. The selective advantages were consistent, within error, when we calculated the velocities ui and uR over different time intervals. We averaged siR across our three experiments and reported its error as the standard error of the mean as seen in Table 1. The eCFP and eYFP strains had an average selective advantage of 9%, similar to the experiments of Weber et al. [35] which found, despite the fact that they used different E. coli strains and plasmids, that the expression of mCherry decreased the expansion velocity of their strains by approximately 15% in certain “fast growth” environmental conditions. Our black strain had an approximately 6% enhancement over the mCherry strain. Differences in radial expansion velocities of this magnitude have been used to study yeast S. cerevisiae and E. coli range expansions in the past [9,17]. To investigate the source of this fitness defect, we took the plasmids from our original strains, inserted them into a different set of clonal DH5α cells, and inoculated the new eCFP, eYFP, and mCherry strains in equal proportions in a range expansion. We saw that the average mCherry fraction decreased by 10% at a radius expanded of R = 10 mm, matching the results of Fig 2, suggesting that the presence of the plasmids was responsible for the fitness defect. From Table 1, it is clear that the variance in siR was large between different sets of agar plates. Although siR varied significantly, the order of expansion velocities between the strains was consistent; the eCFP and eYFP strains always expanded faster than the black strain which expanded faster than the red strain. Importantly and in stark contrast to siR, the demixing radius R0, wall velocities v w i j, and diffusion coefficient Dw were very consistent between sets of plates (measured below), resulting in consistent evolutionary dynamics between our competing strains. To test if the radial expansion velocity differences were related to the basal growth rates of our strains in liquid-culture, we competed all of our strains against mCherry in 10mL flasks of 2xYT growth media. We created three independent replicates of each pairwise competition (9 tubes in total) in the flasks by inoculating 60% of mCherry and 40% of the other strain in mid-log phase. We passaged saturated E. coli samples into new media every 24 hours; we determined the population composition of each flask using a BD LSR Fortessa FACS machine when the cells were passaged. We competed the strains for 72 hours, corresponding to approximately 45 generations (doubling times). We found that every strain grew faster than mCherry in every replicate in liquid culture as judged by a decrease in fraction of mCherry over time. Following previous work [43], we used the decrease in mCherry fraction to determine the dimensionless “well-mixed” selective advantage s i R w m = g i / g R − 1 of the strain competing against it, where gi is the growth rate of strain i and gR is the growth rate of mCherry. We list the measured values of s i R w m in Table 1. The radial expansion velocity fitness siR did not agree with the liquid-culture fitness s i R w m within error, in contrast to previous experiments with the yeast Saccharomyces cerevisiae [17]. However, every strain in liquid culture still grew faster than mCherry. Interestingly, in well-mixed culture, the black strain had the largest growth rate followed by eCFP and eYFP (they had the same growth rate) and then mCherry, disagreeing with the order of radial expansion velocities (where black expanded slower than eCFP and eYFP). As the growth rate differences were small, it is possible that additional factors allowed the eCFP and eYFP strains to expand faster than black on agar. Our E. coli switched from log to stationary phase in the 24 hour cycle; the changing environment may have resulted in a different order of fitnesses compared to the agar plates as well. Future work should investigate how the eCFP and eYFP strains expanded faster than the black strain despite a smaller basal growth rate and should also investigate how such small growth rate differences in liquid culture resulted in such large differences in radial expansion velocities on solid agar. When calibrating our model to experiment, the precise value of R0 did not matter as long as each strain’s local fraction could be accurately measured at that radius. Therefore, to maximize the length over which we could quantify range expansion growth, we defined the local fixation radius R0 as the minimum radius where our image analysis package became accurate. For R < R0, our package predicted equal fractions of each strain due to the overlap of each channel in the homeland (see Fig 10). Therefore, to determine R0, we inoculated radial expansions with three strains in unequal proportions; we used 10% of two strains and 80% of another. The minimum radius where the fractions agreed with their inoculated values was R0 = 3. 50 ± 0. 05 mm as seen in Supplementary S6 Fig. We found that this value of R0 worked for all colonies. Past work has found that E. coli colony domain walls fluctuate diffusively in certain conditions [18] and super-diffusively in others [8]. In our expansions, the domain walls appeared to fluctuate super-diffusively (as judged by tracking the position of domain walls and determining their variance vs. length expanded), but we were able to successfully fit the evolutionary dynamics using a diffusive theory. Creating a super-diffusive theory to describe the evolutionary dynamics of our system is beyond the scope of this paper. To obtain an effective diffusion constant Dw and to test if the diffusive approximation adequately described our experimental dynamics, we fit the neutral Voter model’s prediction of heterozygosity. The heterozygosity is the probability that two points separated by an angle of ϕ at a length expanded of L = R − R0 are occupied by different strains and is thus a measure of spatial genetic diversity. The neutral Voter model’s prediction of heterozygosity can be given in terms of the two-point correlations used in the main text or can be explicitly written as (see S1 Appendix) H (ϕ, L) = ∑ i ∑ j ≠ i F i j (ϕ, L) = H 0 erf [ (1 + R 0 / L) | ϕ / ϕ c | ]. (11) H0 = H (ϕ, L = 0) is the heterozygosity when L = 0 and ϕ c = 8 D w / R 0 is a characteristic angular correlation length (one of the key combinations of model parameters from the main text). For q colors inoculated in equal fractions, H0 = 1 − 1/q. We fit H (ϕ, L) to our experimentally measured heterozygosity of two neutral strains (eCFP and eYFP) on three independent sets of agar plates each with 14 range expansions. We averaged the heterozygosity at each L as can be seen in Fig 12 (error bars were omitted for readability; the same figure with error bars can be found as Supplementary S7 Fig). As we had previously measured R0 = 3. 50 ± 0. 05 mm, and H0 = 1/2 for two neutral strains inoculated at equal fractions, Dw is the single free parameter in eq (11). We consequently fit Dw at each L with non-linear least-squares, averaged the Dw from the three independent experiments, and found Dw = 0. 100 ± 0. 005 mm; the reported error is the standard error of the mean between the experiments. The value of the diffusion constant is on the same order of magnitude as that from previous work [18]. Fig 12 shows the Voter model’s fit (dashed lines) together with the experimental heterozygosity (solid lines) for one set of plates using our values of Dw and R0. The fit closely matches the experimental heterozygosity suggesting that a diffusive description of E. coli domain motion is justified. We use this value of Dw for all strains. In principle, Dw may depend on ij, the particular domain wall type. However, we checked that the measured value of Dw did not vary for our all ij (all strain) combinations by examining the variance in domain wall position versus length expanded; the variances agreed within error and were thus consistent with a constant Dw. The two-point correlation functions in the main-text were well fit by a constant Dw as well. Unlike the Voter model and our simulations, the experimental heterozygosity at zero separation H (L, ϕ = 0) fails to vanish due to overlap between strains at domain boundaries; this effect is less pronounced at large radii because the effective angular width of boundaries decreased. The discrepancy between the theoretical and experimental heterozygosity is larger at small lengths expanded because the overlap between strains is larger; our image analysis is consequently less accurate. We used image analysis to directly quantify v w i j from the angular growth of more-fit sectors. Characteristic single sectors of each strain sweeping through the mCherry strain can be seen on the left side of Fig 13. In radial expansions, more fit strains should, on average, sweep logarithmic spirals through less fit strains at large lengths expanded, as verified in yeast expansions [17]. It can be shown that the average angular width of a sector of strain i sweeping through strain j is given by (see S1 Appendix for more details) ⟨ ϕ - ϕ 0 ⟩ = 2 v w i j ln (R / R 0) (12) where ϕ is the angular width at radius R and ϕ0 is the initial angular width of the domain at R0. 2 v w i j can thus be extracted from the slope of a linear regression fit of 〈ϕ − ϕ0〉 vs. ln (R/R0) as seen on the right side of Fig 13. By tracking domain walls directly, we found that more fit strains (eCFP, eYFP, black) swept through the less fit mCherry strain with a wall velocity of v w i R = 0. 06 ± 0. 02. We could not accurately measure the wall velocity of the eCFP and eYFP strains sweeping through the black strain. The wall velocity was significantly smaller than expected from the basal independent expansion velocities of our strains (Table 1); potential explanations for this phenomenon are discussed in S2 Appendix. The magnitude of the velocities were consistent between experiments (using 40 single sectors on three sets of plates) but were too imprecise to be predictive in our models. Lattice simulations of range expansions, especially radial ones, can suffer from artifacts arising from the preferred directions of the lattice. It is possible to use an amorphous Bennett model lattice [44] to mitigate some of these effects [32]. Instead, we developed a simple off-lattice method that treats the domain walls as annihilating and coalescing random walkers moving along the edge of an inflating ring. The basic idea of the simulation is illustrated in Fig 14. We incorporate both the random, diffusive motion of the domain walls as well as deterministic movement due to selection. The radial expansion procedure is most easily understood by first considering a linear range expansion simulation for which the simulation steps are as follows: Our simulation’s diffusion coefficient per length expanded, characterizing the random motion of the domain walls, can be shown to be Dw = a/2 when rij is small while its wall velocity per length expanded, characterizing the deterministic displacement of domain walls due to selection, can be shown to be v w i j = r i j ≤ 1. Our algorithm, thus far simulating only linear expansions, can easily be extended to simulate radial geometries. To incorporate the radially inflating perimeter, we note that a domain wall at a radius R will jump an angular distance of δϕ = a/R, as shown in Fig 14. As the radius of our experimental expansions increases approximately linearly with generation time, we describe its radius as R = R0 + at. We thus account for inflation by using a time-varying angular jump length of δ ϕ (t) = a R 0 + a t. (13) If there are N0 individuals at the frontier, R0 is given by R0 = N0a/ (2π). This modification of the domain wall step size δϕ is the only difference between the radial and linear cases! In contrast to algorithms that follow the position and state of every organism at the front of a colony, our algorithm only tracks the positions of domain walls and is consequently much faster per generation as the sectors coarsen, allowing for simulations of larger colonies. Fig 15 displays a radial and linear simulation with three neutral colors and a fourth red color with a selective disadvantage comparable to our experiments. We check that our simulation correctly reproduces the behavior of a single more fit domain wall sweeping through a less fit strain as we vary simulation parameters in Supplementary S8 Fig. Our implementation of this algorithm and examples of how to use it are available on GitHub [34].
Population expansions occur naturally during the spread of invasive species and have played a role in our evolutionary history when humans migrated out of Africa. We use a colony of non-motile bacteria expanding into unoccupied, nutrient-rich territory on an agar plate as a model system to explore how an expanding population’s spatial structure impacts its evolutionary dynamics. Spatial structure is present in expanding microbial colonies because daughter cells migrate only a small distance away from their mothers each generation. Generally, the constituents of expansions occurring in nature and in the lab have different genetic compositions (genotypes, or alleles if a single gene differs), each instilling different fitnesses, which compete to proliferate at the frontier. Here, we show that a random-walk model can accurately predict the dynamics of four expanding strains of E. coli with different fitnesses; each strain represents a competing allele. Our results can be extended to describe any number of competing genotypes with different fitnesses in a naturally occurring expansion as long as the underlying motility of the organisms does not cause our model to break down. Our model can also be used to precisely measure small selective differences between spatially competing genotypes in controlled laboratory settings.
bacteriology
heterozygosity
organismal evolution
population genetics
radii
microbiology
mathematical models
geometry
simulation and modeling
luminescent proteins
mathematics
yellow fluorescent protein
microbial evolution
population biology
research and analysis methods
random walk
imaging techniques
proteins
mathematical and statistical techniques
image analysis
genetic drift
biochemistry
heredity
genetics
biology and life sciences
physical sciences
evolutionary biology
bacterial evolution
evolutionary processes
Genetic drift and selection in many-allele range expansions