Mutability of demographic noise in microbial range expansions

Demographic noise, the change in the composition of a population due to random birth and death events, is an important driving force in evolution because it reduces the efficacy of natural selection. Demographic noise is typically thought to be set by the population size and the environment, but recent experiments with microbial range expansions have revealed substantial strain-level differences in demographic noise under the same growth conditions. Many genetic and phenotypic differences exist between strains; to what extent do single mutations change the strength of demographic noise? To investigate this question, we developed a high-throughput method for measuring demographic noise in colonies without the need for genetic manipulation. By applying this method to 191 randomly-selected single gene deletion strains from the E. coli Keio collection, we find that a typical single gene deletion mutation decreases demographic noise by 8% (maximal decrease: 81%). We find that the strength of demographic noise is an emergent trait at the population level that can be predicted by colony-level traits but not cell-level traits. The observed differences in demographic noise from single gene deletions can increase the establishment probability of beneficial mutations by almost an order of magnitude (compared to in the wild type). Our results show that single mutations can substantially alter adaptation through their effects on demographic noise and suggest that demographic noise can be an evolvable trait of a population.

The front roughness was measured using the method described in Ref [7]. A backlight image 30 of the colony was taken at 24 hours, and a custom-written MATLAB code was used to 31 extract the colony boundary using image segmentation. The boundary was fit with a circle 32 and the mean squared displacement in the radial direction was calculated for windows of The growth layer depth was measured using bead displacements between pairs of consecutive 39 timepoints ( Figure S12c). We assumed that the velocity v along the direction of growth has where v y,max is the maximum velocity, y is the position along the direction of growth, y 0 is 42 the front position, and λ is the growth layer depth. The ratio v at two positions within the 43 colony then conveniently takes on the form 44 v = v y (y 2 ) v y (y 1 ) = e −(y 2 −y 1 )/λ .
The average direction of motion of all the trajectories in a field of view was determined, 45 and all subsequent measurements were projected along this axis. For each pair of consecutive 46 timepoints, and for each pair of beads, the ratio of their displacements, vy(y 2 ) vy(y 1 ) , and the calculated. Because noise from inaccurate tracking can dominate the ratio of the velocity 49 when one or both of the bead velocities is small, we restricted our analysis to −50µm < 50 ∆y < 50µm. To minimize the effect of changes in curvature of the front of the colony, we 51 restricted pairs of beads to have ∆x < 50µm. We binned the data (bin size = 10µm, or if 52 range of ∆y < 100µm then we used a bin size of std(∆y)/2) and fit to an exponential decay 53 function using weighted least squares to extract the growth layer depth λ ( Figure S12d). 54 The error in λ represents the error in the fitting parameter. The thickness of the colony was determined using a backlight brightfield image of the whole 57 colony at 23h, where the exposure time was fixed for all colonies. We calibrated the intensity 58 of the transmitted light to the colony thickness using a colony that was grown with fluorescent 59 tracer beads (Methods in main text). We noticed that beads primarily rise to the top of 60 the colony, so we used the height of the fluorescent beads as the ground truth colony height.

61
To measure the average thickness of the whole colony, we calculated the average brightfield 62 intensity within each colony and converted it into a thickness using the calibration curve. 63 We note that this method does not distinguish between changes to intensity due to changes 64 in colony thickness and changes in biomass density, but we make the assumption that any Image segmentation was done with the software Morphometrics [16] and features that were 83 incorrectly segmented were manually rejected. The circularity C is defined as where l is the cell contour length, and A is the cell area. The circularity of a circle is 1 and 85 values larger than 1 correspond to more elongated cells. In order to track the movement of the cells and beads at the single-cell resolution ( Figure   88 1c and Figure S1c), we grew a colony following the procedure in the main experiments

142
The confidence interval was taken from the 2.5% to the 97.5% of the empirically measured beads are at a lower spatial density than cells, we are able to image with an air objective with 162 a lower NA and a lower magnification than if we imaged single cells, and we can also capture 163 a larger field of view. We can also track the beads with images that are taken less frequently 164 because beads move fewer pixels at a lower magnification, and this allows us to image 165 many colonies in parallel. We note that we tried beads with different coatings (unmodified 166 functional group, aminated functional group, Concanavalin A-coated) and the beads with 167 unmodified functional groups yielded the longest bead trajectories. We interpreted this result 168 as the unmodified beads were the best at following the cell trajectories rather than getting lost behind the front. Thus, we chose to use unmodified beads for our experiments.

170
As a first approach to validate the method, we compared the measurement of MSD from 171 sector boundaries, single cell trajectories, and bead trajectories for E. coli DH5α and S.

172
cerevisiae W303 ( Figure S1e), which are two species whose strengths of demographic have 173 previously been compared [7,8]. In all three methods, E. coli has a higher MSD than S.

177
There is a continuous transition from the MSD of E. coli single cells to beads; however,

178
there is a discontinuity between the transition from bead trajectories to sector boundaries.

179
Additionally, there are discontinuities in the transition between the three methods for S.  and Figure S1c). We unfortunately were not able to image for a longer than 1 hour because 189 the images went out of focus. Additionally, we note that adding a coverslip on top of the 190 colony, which was necessary for using a high NA oil immersion objective to get single cell 191 resolution, may change the mechanical behavior of the cells.

192
As a third approach to validate that the method can measure demographic noise, we   µm (see Figure S7). The fit at 25 µm, 50 µm, and 100 µm gave similar reduced chi-squared 217 values, while the fit at 1000 µm had a much higher reduced chi-squared. The poor fit at 218 high window sizes is likely due to the fact that few bead trajectories become that long, 219 and the extrapolation becomes noisy. Note that while the reduced chi-squared values are 220 generally high, we were not able to account for the errors in MSD in the reduced chi-squared 221 calculation; taking those errors into account should reduce the discrepancies in the summary 222 statistic describing the deviation between the fit and the data. 223 We also show the distribution demographic noise for the wild type and knockout strains 224 using the MSD reported at different window sizes in Figure S8. At all window sizes except 225 the largest window size of L = 1000µm (where the MSD also did not agree as well with 226 the fraction of diversity preserved) we see that the wild type and knockout distributions are 227 significantly different from one another to at least a level of 5% (Kolmogorov-Smirnov test). we hypothesize that its discrepancy from its replicates may be due to higher plate moisture.
242 Figure S3b shows that the Pearson correlation coefficient for strains grown in different posi-243 tions on the same plate was 0.55-0.67. Thus, we observe comparable levels of variation due 244 to differences in a strain's position on the plate and differences between plates. 245 We also tested for variation that may result due to the noisiness of averaging a finite 246 number of bead trajectories for each colony. For the selected knockout strains, we randomly 247 split the bead trajectories in each colony's field of view in half, and calculated the MSD 248 separately for each set of trajectories. Figure S3c shows that the Pearson correlation coeffi-249 cient between the two random sets of half of the trajectories is 0.78. Thus, we see that while 250 finite track numbers can lead to variation between replicates, it does not play the largest role 251 in determining variation between replicates, as differences between plates generally leads to 252 more variation.  To better understand if the number of colonies with beneficial sectors is surprising, we 272 can estimate the expected number of de novo beneficial mutations from growth in the colony: where µ b is the beneficial mutation rate per genome per generation, p est is the establishment

295
From this analysis, we did not find any significant GO terms, and we only found one significant KEGG term that was enriched for higher MSD than the weighted median of the 297 entire knockout distribution: ATP-binding cassette transporters. It is possible that expand-  The tail at the lower end of the distributions for all three additionally tested sets of 320 strains does not overlap with the gaussian fit to the wild type distribution nor do they 321 completely overlap the grayed out region from the finite number of wild type measurements.

322
The selected single gene knockout distribution and the mreB point mutants distribution are 323 significantly different from that of the wild type using a two-sample Kolmogorov-Smirnov 324 test (p = 1.9 × 10 −10 for selected knockouts, 6 × 10 −5 for mreB point mutants), while the 325 mrdA point mutant library is less significantly different from the wild type distribution 326 (p = 6.8 × 10 −2 ), most likely due to the fewer number of measurements.

327
The distributions of the strength of demographic noise for the selected single gene knock-328 outs and the mreB point mutants are wider than that of the randomly selected knockouts.

329
Some strains even have MSD close to zero. Thus, we see that the maximally allowed changes 330 to demographic noise are much more extreme than we have sampled with our randomly se-331 lected strains. Sampling more strains will also lead to a more extreme knockout distribution.

332
Similarly to what's observed for the randomly selected knockouts, many more mutants have  it is unclear whether the cell shape could be different in colonies (the condition in this study).

340
To determine whether these datasets would be applied to our experiments, we measured cell  We tested multiple different initial mutant fractions in order to be able to resolve in-364 dividual sectors for different strain backgrounds and fitnesses. Figure S17b shows the es-  Figure S5: Testing for standing variation vs de novo mutations. Probability of seeing at least 1 beneficial sector in colonies grown from (a) the rearrayed glycerol stock or (b) the original Keio collection glycerol stock received from the National BioResource Project. Colonies were grown from liquid culture that was inoculated either directly by scraping off some glycerol stock (blue) or by first streaking the glycerol stock onto a plate, and picking from a single colony (orange). Otherwise, the experimental conditions were the same as that for measuring the distribution of demographic noise. Dots show individual colonies with none (0) or at least 1 beneficial sector (1) and squares show the average across colonies, which represents the probability of seeing a beneficial sector in a given condition. The error bars represent the standard deviation from binomial sampling. In all conditions, beneficial sectors were observed in at least one colony. This supports the hypothesis that beneficial mutations can arise de novo on the timescale of the experiment. Most likely, both beneficial mutations and standing variation from the glycerol stock contribute to the observation of beneficial sectors in the main experiment.  Figure S11: GO and KEGG terms with significantly different MSD from that of the WT at 5% level. The significance of each term was tested by randomly drawing the same number of measurements from the wild type distribution as the number of genotypes associated with that term in our KO subset. We repeated this 106 times and calculated the fraction of times that the weighted median of the random WT subset was more extreme than that of the true values associated with that term. The p value was calculated as twice this fraction (because in determining whether a value is more extreme, we first check whether the weighted median MSD of the term is greater than or less than the weighted median MSD of the full WT distribution). (Left) Weighted median MSD of significant terms and WT weighted median MSD (blue line). Error bars represent the mean absolute deviation. (Right) q-values after applying Benjamini-Hochberg FDR correction to account for multiple testing.  Figure S12 (preceding page): Additional tests of phenotypic traits (a) Colony area is a highly correlated measure of the optical density of the colony resuspended in liquid, which gives the total (alive and dead) biomass. (b) Colony thickness (see Methods) is not correlated with colony area, showing that area is a good measure of growth, rather than just spreading on the surface of the plate. (c) Schematic of measurement growth layer depth using bead displacements from a single colony's field of view. For a consecutive pair of timepoints, ∆y gives the distance between a pair of beads along the average direction of motion, and and v y (y 1 ) and v y (y 2 ) give the bead displacements projected onto the average direction of motion.
(d) Gray points show measurements for all pairs of beads across all pairs of consecutive time points. Black points give binned measurements. Error bars represent standard error of the mean. The binned points roughly followed an exponential decay function between ∆y = -50µm and ∆y = 50µm, and we fit to an exponential decay function within this range, as measurements of bead pairs that are farther away may be dominated by noise. The decay length, d, is extracted from the fit and taken as the growth layer depth. (e) The MSD for bead trajectories from 12-18 hours as a function of the area of the colony at 12 hours (blue points) and the MSD for bead trajectories from 18-24 hours as a function of the area of the colony at 24 hours (orange points). The same colony over time is connected by a gray line. Bead trajectory MSD mostly decreases with increasing colony area over time. Thus, the positive correlation between colony area and MSD seen across genotypes cannot be explained by changes in the bead trajectories over time. (f) Colony front roughness (see Methods) for E. coli DH5α over time. Colony front roughness increases initially but levels off around 19 hours. Since we measure the colony front roughness at 24 hours, we expect to have passed the time when the behavior of the front roughness is transient. Note that this measurement was done on a colony grown in a 10cm diameter petri dish rather than an omniplate, and the colony may have access to a different total amount of nutrients. However, while the growth condition does not exactly match that in the experiment, we used this data to better understand the approximate timescales of when front roughness saturates in time. (e) Relative bead trajectory MSD in cell shape mutants and those predicted by the best fit Lasso regression model from the main text which primarily includes colony-level traits. (f) Partial correlation of cell aspect ratio when controlling for all other phenotypes for the random set of single gene knockouts from the Keio collection. Correlation of demographic noise to single cell aspect ratio in mreB (g) and mrdA (h) single point mutants is higher than that seen for the single gene knockouts, possibly because these strains were enriched for cell shape differences. Strains and data from Ref  Figure S15: Correlation between the bead MSD and (a) front roughness, (b) growth layer depth, or (c) colony area (same as in Figure 3). Below each plot is the partial correlation between the bead MSD and each phenotypic trait (after controlling for correlations with all other traits).