Stochastic exits from dormancy give rise to heavy‐tailed distributions of descendants in bacterial populations

Variance in reproductive success is a major determinant of the degree of genetic drift in a population. While many plants and animals exhibit high variance in their number of progeny, far less is known about these distributions for microorganisms. Here, we used a strain barcoding approach to quantify variability in offspring number among replicate bacterial populations and developed a Bayesian method to infer the distribution of descendants from this variability. We applied our approach to measure the offspring distributions for five strains of bacteria from the genus Streptomyces after germination and growth in a homogenous laboratory environment. The distributions of descendants were heavy‐tailed, with a few cells effectively ‘winning the jackpot’ to become a disproportionately large fraction of the population. This extreme variability in reproductive success largely traced back to initial populations of spores stochastically exiting dormancy, which provided early‐germinating spores with an exponential advantage. In simulations with multiple dormancy cycles, heavy‐tailed distributions of descendants decreased the effective population size by many orders of magnitude and led to allele dynamics differing substantially from classical population genetics models with matching effective population size. Collectively, these results demonstrate that extreme variability in reproductive success can occur even in growth conditions that are far more homogeneous than the natural environment. Thus, extreme variability in reproductive success might be an important factor shaping microbial population dynamics with implications for predicting the fate of beneficial mutations, interpreting sequence variability within populations and explaining variability in infection outcomes across patients.

In contrast to plants and animals, the offspring distribution is largely unexplored for microorganisms. While it is well known that exponential growth can generate 'jackpots' of mutants at high frequencies (Luria & Delbrück, 1943), far less is known about the distribution of descendants within genetically identical populations.
One reason for this might be that the offspring distribution is seemingly simpler for bacteria undergoing binary fission, since each cell can only leave behind 0 (death), 1 (no doublings) or 2 descendants. However, even clonal populations of bacteria display a distribution of growth rates and lag times, causing them to yield a variable number of offspring after some time (Fridman, Goldberg, Ronin, Shoresh, & Balaban, 2014;Labhsetwar, Cole, Roberts, Price, & Luthey-Schulten, 2013;Wang et al., 2010;Xu & Vetsigian, 2017). In particular, many microorganisms form spores or persister (nongrowing) phenotypes to survive unfavourable environments or disperse (Balaban, 2004;Dworkin & Shah, 2010;Fridman et al., 2014;Vulin, Leimer, Huemer, Ackermann, & Zinkernagel, 2018), and exit from dormancy is often a stochastic process that presumably evolved as a bet-hedging strategy to overcome environmental uncertainty (Sturm & Dworkin, 2015;Villa Martin, Munoz, & Pigolotti, 2019;Xu & Vetsigian, 2017). Importantly, it is unclear how stochastic variability in growth rates and lag times affects genetic drift. One way to study this quantitatively is by examining the distribution of the number of bacteria arising from a single bacterium after a given amount of time τ, where the time τ is substantially longer than the standard doubling time ( Figure 1a). This is a stochastic quantity which can be described by a probability distribution that we term here the 'distribution of descendants'. In a system with seasonality, for example, one might look at this distribution after one season. Defined as such, the distribution of descendants is an important quantity of which little is known for bacteria. Quantifying this distribution and how it varies across species and environments would likely improve our understanding of genetic drift in microbial populations and, ultimately, our ability to correctly interpret the genetic variability observed in sequence data.
Here, we present a scalable methodology for quantifying the distribution of descendants in clonal populations. We used a generalizable barcode tagging approach that enabled us to track descendants from hundreds of subpopulations differing only by a short DNA barcode inserted in their chromosome. We developed analysis methods for determining the distribution of descendants from barcode data and applied these approaches to soil bacteria from the genus Streptomyces. We focused on the Streptomyces because they are ubiquitous saprophytes and have complicated life cycles consisting of germination, mycelial growth (with doubling times of 2-3 hr on rich media) and sporulation (Hopwood, Kieser, Bibb, Buttner, & Chater, 2000). The impact of life cycle stages on the distribution of descendants is particularly poorly understood (Bobek, Šmídová, & Čihák, 2017). Using the variability between replicates, we show that the distribution of descendants is heavytailed-that is some bacteria represent a far greater proportion of the final population than their initial frequency. Furthermore, using microscope time-lapse imaging, we demonstrate that the F I G U R E 1 Measurement of the distribution of descendants arising from a population. (a) Clonal cells, represented by coloured circles, are grown for a period of time (τ) before their relative abundances are measured. (b) The variability in the proportion of descendants between replicate populations of cells is used to determine the distribution of descendants. The distribution of descendants may take on a variety of shapes that have different rates of going to zero in their right tail. A heavy-tailed distribution (solid line) would result in 'jackpots' where individuals have much greater reproductive output than expected based on their initial frequency. (c) In order to track lineages, we constructed a barcoded library of Streptomyces where each spore has a unique 30 base-pair lineage-specific sequence integrated into its chromosome [Colour figure can be viewed at wileyonlinelibrary.com] heavy-tailed nature of the distribution of descendants can, in our case, be largely explained by phenotypic variability in lag time before exponential growth. We then examine the implications of heavily skewed distributions of descendants for the population genetics of microorganisms.

| Construction of barcoded strains of Streptomyces
Oligonucleotides 5′-GATCCACACTCTTTCCCTACACGACGCTCTT CCGATCT-3′ and 5′-S20-N30-AGATCGGAAGAGCGTCGTGTAGG GAAAGAGTGTG/3Phos/ were purchased from Integrated DNA Technologies. The latter oligonucleotide is different for each strain library and contains a unique 20-nucleotide strain barcode (S20), a stretch of 30 random nucleotides that form the set of lineage barcodes (N30), and a 3'-phosphate modification. To permit robust identification of a strain in the presence of sequencing errors, the S20 sequences were designed using EDITTAG (Faircloth & Glenn, 2012 (Hopwood et al., 2000) in which the orientation of EcoRV and BamHI sites in the multiple cloning site is reversed.
To reduce the background of pSRKV004 without inserts after This E. coli culture was used for conjugation into the desired Streptomyces strain according to a standard protocol (Hopwood et al., 2000). Briefly, the transformed conjugation helper strain was mixed with Streptomyces spores, and the bacterial mix was grown on mannitol salt (MS) agar for 16 hr and then overlaid with apramycin (100 µg/ml) and nalidixic acid (50 µg/ml). Strains successfully undergoing conjugation integrate the plasmid at a phage attachment site in their genomic DNA (Sun, Kelemen, Fernández-Abalos, & Bibb, 1999). Barcoded libraries were prepared by scraping spores from exconjugants and selecting against E. coli carry-over by propagating the spores on Streptomyces Isolation Medium (Hopwood et al., 2000) supplemented with 50 µg/ml nalidixic acid and 100 µg/ml apramycin for two propagation cycles. Spores were passed through a 5-µm filter membrane to remove any mycelial fragments and disrupt spore clumping.

| DNA extraction and sequencing
After growth, strains were centrifuged at 376 rcf for 10 min to pellet the cells. A 750 µl volume of supernatant was removed, leaving about 150 µl remaining. Note that about 10% of the original volume was lost to evaporation during growth. The remaining volume containing mycelium was sonicated at 100% amplitude for 3 min using a Model 505 Sonicator with Cup Horn (Qsonica) while the samples were completely enclosed. We did not sequence the initial stock (τ = 0) because we were only able to efficiently lyse mycelium and, therefore, cannot accurately determine the relative abundances of the initial barcoded spore populations. After sonication, the samples were centrifuged, and the supernatant containing DNA was used as template for PCR amplification.
The number of unique dual-index Illumina primers available allowed us to sequence eight replicates of each strain at three different initial concentrations. PCR primers (Table S2) were designed with unique 8-nucleotide i5 and i7 index sequences and Illumina adapters. The random barcode (N30) sequence occurs at the start of the sequencing read to assist with cluster detection on the Illumina platform. Since strains could be distinguished by their sequencespecific barcode (S20), we amplified each replicate using a unique dual-index combination, but used the same set of combinations for all five strains. Hence, the S20 region effectively acted as a third index sequence that allowed the five strains sharing dual-index primers to be correctly demultiplexed (Wright & Vetsigian, 2016b

| DNA sequence analysis
Using the R (R Core Team, 2019) package DECIPHER (Wright, 2016), DNA sequencing reads were filtered at a maximum average error of 0.1% (Q30) to lessen the degree of cross-talk between dual-indexed samples (Wright & Vetsigian, 2016b). Sequences were assigned to the appropriate strain by exact matching the S20 sequence and mapping to the nearest barcode by clustering N30 sequences within an edit distance of 5. To completely eliminate any remaining cross-talk (Wright & Vetsigian, 2016b), we subtracted 0.01% + 5 reads from the count of every barcode by sample. The remaining reads were normalized by dividing by the total number of reads per sample.
The final result of this process was a matrix of read counts for each unique barcode across every sample by strain ( Figure S1).

| Time-lapse imaging of the initial growth
To simultaneously track the growth of many Streptomyces colonies, we inoculated spores onto a device developed as part of another study (Xu & Vetsigian, 2017 Images were processed using in-house MATLAB scripts. First, 25% sized images were aligned between time points by identifying shared features using MATLAB's computer vision toolbox. Regions of the image with remaining misalignment were fixed by local image registration. These transformations were then scaled to larger (50% sized) images used for further analyses. Growing colonies were detected by comparing the difference between subsequent time points, and area was determined through thresholding the image since mycelium is darker than the background. Tracking was terminated at the time point before colonies intersected. Manual validation was applied to remove artefacts that were incorrectly identified by the algorithm as mycelium. Finally, growth curves were removed with extreme jumps between time points or decreases in colony size, which were characteristic of localized failures in thresholding.

| Inferring the distribution of descendants from rare barcodes
The barcode reads that we observe result from three sequential stochastic processes ( Figure 1b):

A common initial pool of barcoded spores is randomly sampled
to inoculate replicates.
2. Individuals in each tube undergo stochastic germination and growth.
3. The final barcode pool in each tube is randomly sampled through qPCR and sequencing to generate the observed read counts.
Our goal is to correctly account for processes 1 and 3 in order to learn about process 2, which is the process of biological interest. Process 2 is quantified by the distribution of descendants, D(x), specifying the probability that an individual in a test tube will leave descendants that constitute a relative fraction x of the final population. Let r ib be the number of observed reads of barcode b in sample i (i = 1…8) and R i = ∑ b r ib be the total number of reads for each sample. We want to estimate the probability density function D(x) based on the data set {r ib }.
To model the initial sampling (process 1), let v b be the unknown initial density of barcode b within the common pool. The units are chosen such that v b is the expected number of cells with barcode b in a tube immediately after inoculation. The actual number of initial cells with barcode b in a sample i, n ib , follows a Poisson distribution: To model the final sampling (process 3), let x ib be the unknown relative frequency of barcode b in sample i at the end of stochastic germination and growth ( ∑ b x ib = 1). The observed reads, r ib , result from stochastically sampling R i reads from the multinomial distribution specified by the frequencies {x ib }. This process can be approximated by assuming independent binomial sampling for each barcode, which is justified if there are many rare barcodes. This leads to: Taking a Bayesian view and considering the unobserved quantities x to be a random variable with an uninformative prior (P(x) = const), we have P(x|;r) ∝ P(r|;x). Substituting the above expression and normalizing the probability density, the factor drops out, and we obtain: Using the above mathematical relationships, we perform an iterative two-step procedure for estimating D(x), starting from an initial guess of the probability density. For example, the initial guess can be set to the singleton distribution smoothed by fitting it to a functional form (e.g. log-normal). The prior distribution is updated by alternating steps 1 and 2, described below, until convergence of D(x). Approximately, five iterations were sufficient for convergence with our data.

| Step 1
Given the observed reads and an estimate of D(x), compute prob- Using Bayes' theorem and the independence of different samples, we get:

Furthermore, summing over all possible intermediate states x ib
for going from v b to r ib , we obtain: and, similarly, summing over all possible n ib : We used P prior (v b ) = const. For experiments in which the initial pool of barcoded cells was diluted 10 times, we also tried setting The above numerical procedure requires us to choose a value for n max , the maximum initial number of cells from a barcode in a tube. If we exclude data for barcodes that are present in all tubes (i.e. m = 8), n max = 20 can be safely used, since it is unlikely that barcodes that are missing from some of the tubes have more than a few initial cells in the tubes in which they were detected.

| Step 2
Given an estimate for the probability distributions of v b , the observed reads, and an estimate of the offspring distribution D prev (x), compute an improved estimate for D(x).
If we know n ib and x ib in a given tube i and barcode b, this would provide information for constructing D(x). If n ib = 1, then x ib constitutes an observation sampled from D(x), which will correspondingly contribute to the probability density of D(x) at x = x ib . If n ib = 2, we cannot immediately determine the contribution of x ib to D(x).
However, this can be done by leveraging our previous best estimate for the distribution of descendants, D prev (x). The hope is that, as we iterate, the algorithm will eventually converge to a self-consistent solution D(x) = D prev (x) and, indeed, it does so in practice. Returning to the case of n ib = 2, we have two initial cells and two contributions, y and x ib − y to the probability density D(x), and each particular value of y contributes a relative weight of . More generally, the contributions to D(y) from a (n ib , x ib ) pair are weighed according to where we define D (0) prev (x) = 1. The factor n ib comes from equivalence between the n ib cells in a tube.
Averaging contributions from all possible (n ib , x ib ) pairs according to their Bayesian probability P(n ib ,x ib |;P(v b ),r ib ), we can compute the overall contribution of a barcode b in sample i to the probability density D(x), and then, we can sum over all i and b to incorporate contributions from all the barcodes and samples: The last factor can be computed as with P(n ib ,x ib │r ib ,v b ) = P(n ib |;v b )P(x ib |;r ib ) (as defined above), and P(v b |;{r ib },D prev (x)) being the output of Step 1.

| Modelling the population genetics consequences of the distribution of descendants
The population genetics simulation consists of discrete time steps.
Each time step captures the dynamics over one ecological cycle.
We assume that each individual (i) starts growing exponentially with growth rate r, after a stochastic lag time t i . Correspondingly, in the beginning of each ecological cycle, N random variables t i are drawn from a lag time distribution, and the corresponding relative abundances of descendants are computed as For models with selection, we set x i → 1 + s i x i and normalized the sum to 1. To complete the cycle, N random individuals are selected from the multinomial distribution specified by x i to start the next cycle. To find the variance effective population size, N e , we computationally determined the distribution of descendants, v, after a full ecological cycle, that is the discrete probability distribution for the number of descendants from one individual after one ecological cycle. The mean of v is 1. We then set N e = N var (v) (Ricky Der, Epstein, & Plotkin, 2011). To perform simulations with a log-normal distribution of descendants, we set t to the normal distribution (with variance 1) and used r = 2 or r = 3 as indicated in figure captions. base-pair random sequence inserted at a fixed site on the chromosome ( Figure 1c). A similar technique has been used previously to tag yeast and Escherichia lineages (Cottinet et al., 2016;Levy et al., 2015). After barcoding, we grew five different strains of Streptomyces in eight separate replicate populations starting from three different initial concentrations. Streptomyces strains first germinate and then grow as interconnected filamentous colonies within liquid medium.

| High-throughput measurement of the distribution of descendants
After 7.5 days of growth, genomic DNA was extracted and the barcoded region was amplified before sequencing (see Section 2). We observed between 211 and 2,534 unique barcoded lineages per strain across all replicates in the experiment. An example of the data collected for one of the five strains is depicted in Figure S1.
Since our analysis methods are based on the variability between replicate populations, they require that the technical variability due to the experimental procedure be far less than the biological F I G U R E 2 Inferring the distribution of descendants from abundant and rare barcodes. (a) Barcodes at the two extremes of relative abundance reflect the shape of the distribution of descendants. (b) Abundant barcodes, those shared by more than 1,000 cells in the initial population, are expected to converge to a normal distribution due to the central limit theorem. However, many of the most abundant barcodes were not normally distributed across replicates, based on their p-values (at right) in the Shapiro-Wilk test. Instead, the abundant barcodes originating from three different strains (colours) were widely scattered in terms of their final proportion of the population (x-axis).
(c) The singletons, those barcodes occurring in only 1 out of 8 replicates, approximate the shape of the distribution of descendants since they likely started from single cells. For the strain with the most singletons at high concentration, Streptomyces G4A3 (808 singletons), we observed many outliers where the relative abundance was much higher than expected (~10 -4 ) on average. (d) Plotting the same histogram on a log-scaled x-axis reveals that the outliers fall beyond the right tail of a fitted log-normal distribution (green curve). These outlying 'jackpots' represent cells that grew to a far higher abundance than the median abundance of singletons, in many case by more than 100-fold. Note that the left side of the distribution is likely truncated because it contains barcodes that fall below the lower detection limit of our method [Colour figure can be viewed at wileyonlinelibrary.com]

Singletons
Abundant barcodes licates. We found that technical replicates had substantially higher correlation than biological replicates ( Figure S2), confirming that most of the variability is biological in nature. This allows the shape of the distribution of descendants to be inferred from fluctuations in the relative abundance of barcodes between biological replicates.
However, it is worth noting that we can only observe the right side of the distribution of descendants, because the lower detection limit of our method is approximately 1 in 10 5 cells based on the finite number of sequencing reads and initial templates in PCR. Therefore, we would not observe the rarest barcodes if they decrease in relative frequency substantially during the course of the experiment.
Nonetheless, we are most interested in the right tail of the distribution of descendants because it might include lineages that increase considerably in relative abundance.

| The distribution of descendants is skewed with a heavy tail
Two extremes of the barcode frequency distribution across replicates reveal characteristics of the distribution of descendants (Figure 2a). At one end, the abundance of a barcode present at high frequency is expected to be normally distributed across replicates.
This is because, for abundant barcodes, each barcode represents a large number of initial cells and the final barcode frequency is a sum of many realizations of the distribution of descendants. Based on the central limit theorem, the variation across replicates will approach normality, so long as the underlying distribution of descendants has a tail that decays sufficiently fast. We tested whether the relative frequencies of the eight replicates belonging to the most abundant barcodes could be normally distributed using the Shapiro-Wilk test.
Each of these barcodes is estimated to be shared by over 1,000 initial cells per replicate based on their final fraction of the population F I G U R E 3 Inferring the distribution of descendants from abundant and rare barcodes. (a) The relative abundance of barcodes (x) is shown according to the number of replicates (m) in which they were detected within the eight total replicates of S. SG4A3 inoculated at medium concentration. Coloured points highlight the proportion of specific barcodes that appeared as outliers in one of the (m) samples in which they were detected. (b) The distribution of descendants is reconstructed using the barcodes that were not observed in all replicates (m < 8). CDF(x) is the probability that the relative abundance of all the descendants of one starting cell is less than x. A log-normal distribution (dashed red) fits the body of the distribution apart from several outliers. (c) Given the median abundances of barcodes observed in all eight replicates (m = 8), it is possible to simulate the expected relative abundances of barcodes according to the fitted log-normal distribution. (d) Repeating the process of simulating barcode abundances 100 times results in a cloud of points (grey) representing expected relative abundances according to a log-normal distribution of descendants (black points correspond to 10 repetitions). Overlaying with the relative abundances of barcodes observed in all eight replicates (coloured points) reveals that the log-normal distribution is insufficient to capture outliers. This provides further support to the idea that the distribution of descendants is heavier-tailed than a log-normal distribution [Colour figure can be viewed at wileyonlinelibrary.com] and the empirically determined initial population size. For four out of nine of these abundant barcodes, the normal distribution was rejected with p-value <.02 (Figure 2b). We calculated a combined p-value of .0006 using Stouffer's method, strongly rejecting the normal distribution. Moreover, the fact that we could repeatedly reject a normal distribution even with a small number of replicates indicates that the deviations from normality are strong. Normality would be expected even for a log-normal distribution of descendants starting from a large enough initial population. Thus, this analysis suggests that the underlying distribution of descendants is heavier-tailed than a log-normal distribution.
At the other extreme, as the initial frequency of a barcode approaches a single cell (Figure 2a), the distribution of final barcode frequencies should reflect the distribution of descendants. We made the approximation that barcodes appearing in only 1 out of 8 replicates of a given initial concentration were sufficiently rare to have originated from a single cell. While we would expect this assumption to be violated in about 10% of cases because a rare barcode had an initial population size of 2, the impact of starting from two cells should be on the order of twofold. The resulting distribution of barcode frequencies for these 'singletons' is heavy-tailed and appeared broader than a log-normal distribution (Figure 2c,d). Surprisingly, for many strains the distribution spanned over three orders of magnitude, meaning that some barcodes were over-represented by more than 1,000-fold that of a typical barcode starting from an identical initial frequency.  Figure 3a). For rare barcodes (m < 8), these outliers cannot be due to exceptionally high initial abundance. Moreover, it is apparent that barcodes that are outliers in some samples can have very low abundance or typical abundance in the other replicate populations. This suggests that outliers are not due to strongly beneficial adaptations that have already established in a barcode prior to inoculation. Furthermore, it is unlikely that outliers are caused by nonrandom sampling or clumping of spores in the inoculum because the initial spore stocks were dispersed by passage through a 5-µm filter membrane. Therefore, the outliers are likely a consequence of a heavy-tailed distribution of descendants.

| The distribution of descendants is heavier than log-normal
We sought to develop a procedure for determining the distribution of descendants that is statistically robust and utilizes more of the data than the singleton (m = 1) barcode approach of Figure 2c,d.
We developed a Bayesian method for recovering the distribution of descendants from relatively rare barcodes (m < 8) given the constraint that we cannot determine initial barcode abundances due to difficulties extracting DNA from spores. Our approach is model-free in the sense that it makes no assumptions about the functional form of the distribution of descendants, and it is essentially an expectation-maximization algorithm with the initial barcode frequencies serving as latent variables. While the method works even without prior information about initial barcode abundances, it allows for such information to be easily incorporated when available. For example, we used the inferred initial barcode abundances from the high concentration experiments as a prior (after dilution) for the lower concentration experiments (see Section 2). We used the raw number of sequencing reads per sample, rather than relative abundances, in order to account for any differences in the depth of sequencing across replicas. Our approach is iterative but reliably converges to almost the same distribution of descendants even when starting from different prior distributions. We validated our approach using data simulated according to different distributions of descendants ( Figure S3). We expect that the robust analysis procedure we developed here will be useful in future studies of the distribution of descendants.
The Bayesian analysis method allowed us to quantitatively reconstructed the distribution of descendants from m < 8 barcodes ( Figure 3b). As shown, a log-normal distribution is an excellent fit to the majority of the data set, but there are several outliers suggesting a heavier tail than the log-normal. We observed these outliers across several strains starting from different initial concentrations ( Figure S4). To corroborate the faster than log-normal tail, we used the log-normal fit from m < 8 to make predictions about the m = 8 data ( Figure 3c). Reassuringly, the prediction based on m < 8 data explains most of the variability in m = 8 data (Figure 3d). But once again, the m = 8 data revealed clear outliers that cannot be explained based on a log-normal distribution of descendants. Taken together, these analyses indicate heavier than log-normal tail for the distribution of descendants.

| Stochastic exits from dormancy largely explain the heavy tail
We reasoned that growth variability would largely result from two sources: differences in lag time before growth (driven by variability in germination times) or variability in growth rates that is auto-correlated across divisions. To determine which source dominated growth variability in our experimental system, we tracked strains under a microscope during their first day of growth on agar medium containing the same nutrients as the liquid experiment. This resulted in large images (Figure 4a) that we aligned between time points to track the growth of each germinated spore (see Section 2). Colony growth was constrained to two dimensions for a long time, which allowed us to estimate the number of genomes present from the area covered by the colonies. This method provides a means of directly assessing the distribution of descendants until the colonies intersect and can no longer be distinguished.
Three of the five strains mostly completed germination during the course of the experiment, while two strains germinated too late to adequately track under the microscope. All three early-germinating strains displayed wide variation in colony size after one day of growth, with the largest colonies being almost 3 orders of magnitude larger in biomass than the smallest (Figures 4b, S5). This likely underestimated the extent of variability, as large colonies can easily overwhelm smaller colonies so that they cannot be identified at later time points and because we sampled only hundreds of spores, thus potentially missing rare instances of extremely early germination.

Nevertheless, it was clear that variation in germination times might
largely account for the approximately log-normal shape of the distribution of descendants (Figure 4c). Considering a simple exponential growth model, descendants = e r * (t−t 0 ) , a constant growth rate (r) and normal distribution of germination times (t 0 ) would result in a lognormal distribution of descendants.
It has been previously shown that growth rates are variable immediately after germination (Xu & Vetsigian, 2017). Therefore, it is possible that differences in the initial growth rate could compound the lag time variability to make the distribution even wider than log-normal. According to the simple growth model above, initial variability in r would amplify variability in t 0 to generate a distribution of descendants that is heavier-tailed than a log-normal distribution. To investigate this effect, we projected the distribution of descendants after 24 hr of growth on solid medium ( Figure S5d-f).
For two of three strains, we observed outliers that were much larger in size than expected for a purely log-normal distribution of descendants. This observation supported our finding that the distribution of descendants is heavier-tailed than a log-normal distribution, and suggests that variability in growth rates may compound variability in lag times. It also implies that processes before stationary phase largely account for the outliers observed in the barcoded replicate populations. However, estimating long-term growth rate variability from microscopy data is difficult because of crowding effects. In this manner, our barcoding approach is preferable to microscopy since it provides data over longer durations and for more cells.

| Selection is an implausible explanation for the observed distribution of descendants
A potential source of growth variation is the existence of genetic differences within the population. One genetic basis for variation is that some barcodes have pre-existing mutations that impart a higher growth rate, resulting in an exponential divergence in relative abundance over time. One way to uncover differences in interbarcode selection coefficients is to look for correlations between the final relative frequencies of rare barcodes. We tested this by plotting the relative frequency of barcodes that were only present in 2 of 8 replicates ( Figure S6). If jackpots within this set are due to selection, we would expect them to manifest in both replicates. In contrast, the correlation between replicate barcodes was extremely low (Pearson's r = .08), indicating that interbarcode selection coefficients are not a major source of the observed variability between replicates. A second way to investigate the role of pre-existing adaptive mutations is to examine whether outliers in the lowest abundance lineages are also outliers in the highest abundance lineages because a beneficial barcode subpopulation at medium concentration would likely also be present in a 10-fold more concentrated inoculum.
However, we found that outliers of rare barcodes at medium concentration are not outliers at higher concentrations, further diminishing the likelihood of pre-adapted barcode lineages as an explanation for the observed jackpots ( Figure S7).
These results do not rule out the possibility that there were rare individuals within a barcode lineage with new or recently acquired beneficial mutations. Such mutants would likely have had to arise after the start of the experiment in order to only be present in a minority of replicates. Given the high number of positively skewed replicates, it is implausible that so many mutants of large effect size could occur and reach high abundance so rapidly. Furthermore, we estimate that most cells only doubled about 10-15 times over the course of the experiment, depending on each strain's initial concentration. Even a large growth rate advantage of 10% would be expected to result in at most a threefold variability in final abundances. Nonetheless, it is well known that mutation is a major cause of fitness variation in populations, and we cannot rule out the fact that some of the variance in the distribution of descendants was attributable to genetic differences. Notwithstanding this limitation, the predominant source of variability, spore germination times, has been shown to be a phenotypic effect (Xu & Vetsigian, 2017).

| The heavy-tailed distribution of descendants yields large deviations from classical population genetics predictions
Previous theoretical work has demonstrated that fat-tailed distributions of descendants with infinite variance can lead to qualitatively different population genetic dynamics (Der et al., 2011;Eldon & Wakeley, 2006;Hallatschek, 2018;Schweinsberg, 2003). While our experimentally determined distributions are fatter than log-normal, it is unlikely that the variance diverges with population size for the particular species and conditions we examined. Importantly, however, heavy-tailed distributions of descendants can still lead to significant deviations from classical population genetics predictions for population sizes (N) and selection coefficients (s) of practical interest, even though the deviations will disappear in the infinite population size limit (N→∞) for fixed values of the product sN. To demonstrate this point and examine the nature of the deviations, we performed population genetics simulations using a log-normal distribution of descendants. We modelled a situation in which we start with N individuals, let them grow exponentially to large numbers following a stochastic exit from dormancy, and then sample N individuals at random to start the next ecological cycle (see Figure 5a, Section 2).
We first determined the distribution of descendants after one ecological cycle, that is the distribution, v(n), for the number of individuals, n, descending from one individual after one cycle. Classical population genetics theory states that the consequences of genetic drift can be captured by a Fisher-Wright model with a (variance) effective population size N e = N/var (v). Figure 5b shows that the heavy-tailed distribution of descendants leads to N e «N, and strikingly, that N e does not scale linearly with N initially, though linearity is ultimately restored given that the assumed log-normal distribution of descendants has a finite variance.
A heavy-tailed distribution of descendants affects the population dynamics beyond reducing N e . We examined how purely neutral dynamics differ from the predictions of a Fisher-Wright model with the same N e . To this end, we computed the distribution of fixation times for a neutral allele starting at 50% abundance ( Figure 5c). We found that neutral mutations take significantly longer to fix with a heavy-tailed distribution of descendants. We next determined the fixation probability of beneficial mutations for different selection coefficients (s) and population sizes (N). We observed (Figure 5d) that as s increases for fixed N, the probability of fixation increases super-linearly and thus more rapidly than expected based on the classical population genetics prediction (for haploid populations) (Kimura, 1962), which states: The classical formula with matching N e agrees with the simulations only for small s since it predicts linear dependence on s for large N (before saturation). Notice that this discrepancy persists at fixed s even as we increase N, despite the fact that the classical limit is recovered at fixed sN (or sN e ) when N → ∞ (Figure 5e). As a control, we verified that this formula agrees well with simulations of the Fisher-Wright model with matching N e across the range of s values ( Figure S8). Therefore, while heavy-tailed distributions of descendants dramatically decrease the effective population size, they increase the efficiency of selection at a given effective population size for mutants with large s.

| D ISCUSS I ON
In this study, we developed and applied a scalable procedure for determining the distribution of descendants arising from a population of bacteria. Surprisingly, the distribution of descendants was heavy-tailed, resulting in a wide range of relative abundances after only a short time (Figure 3). This variation was largely explained by differences in lag time before exponential growth (Figure 4) (Balaban, 2004;Sturm & Dworkin, 2015;Xu & Vetsigian, 2017), it is still unknown how common such stochasticity is among microorganisms. However, it has been argued, for example in the context of desert plants (Gremer & Venable, 2014), that stochastic exit from dormancy is a bet-hedging strategy that increases survival in uncertain environments. Given the generic nature of this argument, it is likely that stochastic exists from dormancy are common across the tree of life. We therefore expect that the findings described here will be relevant to many microbial populations and will stimulate further work on stochastic germination.  The heavy-tailed nature of the distribution of descendants is anticipated to have several effects on bacterial populations. First, extreme stochastic variability can decrease the effective population size dramatically below the census population size (Hedrick, 2005), even when the census size is measured at population bottlenecks within ecological cycles. Moreover, our experimental results supported a population genetics model in which the discrepancy between census and effective population sizes increases with the number of individuals and, therefore, becomes more important for large systems. Such processes can greatly amplify the effects of genetic drift and lead to faster elimination of genetic diversity, larger fluctuations of allele frequencies and an increased lower bound at which weak selective pressure can effectively act. In particular, amplified genetic drift may influence microbial population dynamics on timescales that are important to commercial biotechnologies or bacterial infections. Second, classical population genetics models with matching variance effective population size do not adequately represent dynamics in a population with a heavy-tailed distribution of descendants. We showed that the probability of fixation of beneficial mutations increases faster than linearly with the selection coefficient, and for mutants with large selection coefficients, deviations persist as we increase the population size. Third, since many infections are caused by a small initial number of cells or viruses, wide distributions of descendants may greatly influence the early burden on the host and partly explain the variability in symptoms observed between patients with the same infection. Finally, our results offer support for the notion that true fitness, that is the long-term propensity to have more descendants, is difficult to measure (Mills & Beatty, 1979). Even the largest subpopulations in our experiments exhibited variability in their relative abundance between replicates due to jackpots. Owing to insufficient replication or low initial population size, this variability could easily be interpreted as a long-term heritable fitness difference when potentially none is present.
While, to our knowledge, this is the first measurement of a distribution of descendants for bacteria, it is known that viruses also exhibit large variation in the number of progeny generated from each infected cell. For example, human cells can differ by up to 300-fold in the number of released viruses depending on the stage of the cell cycle in which the infection occurs (Russell, Trapnell, & Bloom, 2018;Timm & Yin, 2012;Zhu, Yongky, & Yin, 2009). The methodology employed here for tracking Streptomyces could be extended to study the distributions of descendants for other species and environments. It would be particularly interesting to determine the distribution of descendants of bacterial populations in their natural environment or as part of the human microbiome, where additional complexities might further broaden the distribution relative to the homogeneous environment explored in this study. We anticipate that heavy-tailed distributions of descendants are common in microbial ecology, as it is possible to envision many sources of variation (e.g. persister states and spatial structuring) that would be amplified by bouts of exponential growth.
Here, we focused on measuring the right tail of the distribution of descendants where ancestral cells produce more progeny than expected. We were unable to observe the extreme left tail of the distribution where cells did not produce progeny or experienced extreme die-offs. It is believed that most major changes in animal population abundances are due to die-offs rather than population explosions (Anderson, Branch, Cooper, & Dulvy, 2017). Since evolution is a function of both birth rates and death rates (Doebeli, Ispolatov, & Simon, 2017), it would be similarly meaningful to measure the left tail of the distribution of descendants. This tail of the distribution is difficult to capture but may shed light on the relative importance of death in evolutionary processes.

ACK N OWLED G EM ENTS
We Ives and Laurence Loewe during the preparation of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing financial or nonfinancial interests.

AUTH O R CO NTR I B UTI O N S
EW performed the experiments. EW and KV designed the study, analysed the data and wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The sequence data are available from the Short Read Archive (SRA) repository under accession number PRJNA353868.