Genetic Paths to Evolutionary Rescue and the Distribution of Fitness Effects Along Them

Abstract

Novel environments can cause strong selection and rapid adaptation. The genetic basis of such rapid adaptation tends to be composed of few loci of large effect. Current theory qualitatively agrees but largely neglects the demographic...

The past century has seen substantial theoretical and empirical progress on the genetic basis of adaptation. Over this same period, a pressing need to prevent the evolution of drug resistance has uncovered much about the potential genetic basis of persistence in declining populations. However, we have little theory to predict and generalize how persistence—by sufficiently rapid adaptation—might be realized in this explicitly demographic scenario. Here, we use Fisher’s geometric model with absolute fitness to begin a line of theoretical inquiry into the genetic basis of evolutionary rescue, focusing here on asexual populations that adapt through de novo mutations. We show how the dominant genetic path to rescue switches from a single mutation to multiple as mutation rates and the severity of the environmental change increase. In multi-step rescue, intermediate genotypes that themselves go extinct provide a “springboard” to rescue genotypes. Comparing to a scenario where persistence is assured, our approach allows us to quantify how a race between evolution and extinction leads to a genetic basis of adaptation that is composed of fewer loci of larger effect. We hope this work brings awareness to the impact of demography on the genetic basis of adaptation.

Antimicrobial drug resistance, Evolutionary escape, Fisher’s geometric model, Genetic basis of adaptation, Mathematical theory

OUR understanding of the genetic basis of adaptation is improving rapidly due to the now widespread use of genomic sequencing (see examples in Bell 2009; Stapley et al. 2010; Dettman et al. 2012; Schlötterer et al. 2015). A recurrent observation, especially in experimental evolution with asexual microbes, is that the more novel the environment and the stronger the selection pressure, the more likely it is that adaptation primarily proceeds by fewer mutations of larger effect (i.e., that adaptation is oligogenic sensu  Bell 2009). An extreme case is the evolution of drug resistance, which is often achieved by just one or two mutations (e.g., Bataillon et al. 2011; Pennings et al. 2014).

However, drugs, and other sufficiently novel environments, will often induce not only strong selection but also population decline. Such declines hinder both the production and maintenance of adaptive genetic variation (Otto and Whitlock 1997), thus impeding evolution and threatening extinction. Drug resistance evolution is a particular instance of the more general phenomenon of evolutionary rescue (Gomulkiewicz and Holt 1995; Bell 2017), where persistence requires sufficiently fast adaptive evolution.

Most theory on the genetics of adaptation (reviewed in Orr 2005) assumes constant population size, and, therefore, does not capture the characteristic “race” between adaptation and extinction that occurs during evolutionary rescue. Many models have been created to describe this race (reviewed in Alexander et al. 2014), but, so far, largely focus on two extreme genetic bases, both already introduced in Gomulkiewicz and Holt (1995): rescue is either caused by minute changes in allele frequencies across many loci in sexuals (i.e., the infinitesimal model; Fisher 1918) or by the substitution of a single large effect “resistance” mutation (e.g., one locus, two allele models). We therefore largely lack a theoretical framework for the genetic basis of evolutionary rescue that captures the arguably more realistic situation where an intermediate number of mutations are at play (but see exceptions below). The near absence of such a framework prevents us from predicting the number of mutations that evolutionary rescue will take and the distribution of their effect sizes. The existence of a more complete framework could therefore provide valuable information for those investigating the genetic basis of drug resistance (e.g., the expected number and effect sizes of mutations) and would extend our understanding of the genetic basis of adaptation to cases of nonequilibrial demography (i.e., rapid evolution and “eco-evo” dynamics).

Despite these gaps in the theory on the genetic basis of evolutionary rescue, there is a wealth of data. For example, the genetic basis of resistance to a variety of drugs is known in many species of bacteria (reviewed in MacLean et al. 2010), fungi (reviewed in Robbins et al. 2017), and viruses (reviewed in Yilmaz et al. 2016). This abundance of data reflects both the applied need to prevent drug resistance and the relative ease of isolating the genotypes that survive (hereafter “rescue genotypes”), e.g., in a Luria-Delbrück fluctuation assay (reviewed in Bataillon and Bailey 2014). Assaying fitness in the environment used to isolate mutants (e.g., in the drug) then provides the distribution of fitness effects of potential rescue genotypes. Additional data on the genetic basis of drug resistance arise from the construction of mutant libraries (e.g., Weinreich et al. 2006) and the sequencing of natural populations (e.g., Pennings et al. 2014). Together, the data show that resistance often appears to arise by a single mutation (e.g., MacLean and Buckling 2009; Gerstein et al. 2012; Lindsey et al. 2013) but not always (e.g., Bataillon et al. 2011; Pennings et al. 2014; Gerstein et al. 2015; Williams and Pennings 2019). The data also indicate that the fitness effect of rescue genotypes is more often large than small, creating a hump-shaped distribution of selection coefficients (e.g., Kassen and Bataillon 2006; MacLean and Buckling 2009; Gerstein et al. 2012; Lindsey et al. 2013; Gerstein et al. 2015) that is similar in shape to that proposed by Kimura (1983) [see Orr (1998), for more discussion], but with a lower bound that is often much greater than zero.

Theory on evolutionary rescue (reviewed in Alexander et al. 2014) has focused primarily on the probability of rescue rather than its genetic basis. However, a few studies have varied the potential genetic basis enough to make some inference about how evolutionary rescue is likely to happen. For instance, in the context of pathogen host-switching, Antia et al. (2003) numerically explored the probability of rescue starting from a single ancestral individual when k sequential mutations are required for a positive growth rate, each mutation occurring from the previous genotype with the same probability and all intermediate genotypes being selectively neutral. The authors found that rescue became less likely as the number of intermediate mutations increased, suggesting that rescue will generally proceed by the fewest possible mutations. This framework was expanded greatly by Iwasa et al. (2004a), who allowed for arbitrary mutational networks (i.e., different mutation rates between any two genotypes) and standing genetic variation in the ancestral population. Assuming the probability of mutation between any two genotypes is of the same order, they showed that genetic paths with fewer mutational steps contributed more to the total probability of rescue, again suggesting rescue will occur by the fewest possible mutations. Iwasa et al. (2004a) also found that multiple simultaneous mutations (i.e., arising in the same meiosis) can contribute more to rescue than paths that gain these same mutations sequentially (i.e., over multiple generations) when the growth rates of the intermediate mutations are small enough, suggesting that rare large mutations can be the most likely path to rescue when the population is very maladapted or there is a fitness valley separating the wild-type and rescue genotype. This point was also demonstrated by Alexander and Day (2010), who emphasized that multiple simultaneous mutations become the dominant path to rescue in the most challenging environments. As a counterpoint, Uecker and Hermisson (2016) explored a greater range of fitness values in a two-locus two-allele model, showing that, with standing genetic variation, rescue by sequential mutations at two loci (two mutational steps) can be more likely than rescue by mutation at a single locus (one simultaneous mutational step), particularly when the wild type is very maladapted, where the single mutants can act as a buffer in the face of environmental change. In summary, current theory indicates that the genetic basis of rescue hinges on the chosen set of genotypes, their fitnesses, and the mutation rates between them. So far these choices have been in large part arbitrary or chosen for mathematical convenience.

Here, we follow the lead of Anciaux et al. (2018) in allowing the genotypes that contribute to rescue, as well as their fitnesses and the mutational distribution, to arise from an empirically justified fitness-landscape model (Tenaillon 2014). In particular, we use Fisher’s geometric model to describe adaptation following an abrupt environmental change that instigates population decline. There are two key differences between this approach and earlier models using Fisher’s geometric model (e.g., Orr 1998): here (1) the dynamics of each genotype depends on their absolute fitness (instead of only on their relative fitness), and (2) multiple mutations can segregate simultaneously (instead of assuming only sequential fixation), allowing multiple mutations to fix—and in our case, rescue—the population together as a single haplotype (i.e., stochastic tunnelling, Iwasa et al. 2004b). In this nonequilibrium scenario, variation in absolute fitness, which allows population size to vary, can create feedbacks between demography and evolution, which could strongly impact the genetic basis of adaptation relative to the constant population size case. In contrast to Anciaux et al. (2018), our focus here is on the genetic basis of evolutionary rescue and we also explore the possibility of rescue by mutant haplotypes containing more than one mutation. In particular, we ask: (1) how many mutational steps is evolutionary rescue likely to take, and (2) what is the expected distribution of fitness effects of the surviving genotypes and their component mutations?

We first introduce the modeling framework before summarizing our main results. We then present the mathematical analyses we have used to understand these results, and end with a discussion of our key findings.

Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article and Supplementary Material. Supplementary figures are provided in File S1. Code used to derive analytical and numerical results and produce figures (Mathematica, version 9.0; Wolfram Research Inc. 2012) is provided as File S2. These files, as well as code used to run individual-based simulations (Python, version 3.5; Python Software Foundation), have all been deposited at figshare: https://doi.org/10.25386/genetics.11343704. All these files, as well as simulation data and freely accessible versions of File S2 (CDF and PDF), are also available at https://github.com/mmosmond/GeneticBasisOfRescue.

Model

Frequently used notation in our mathematical analysis are listed in Table 1.

Frequently used notation

Table 1

Frequently used notation

Symbol	Meaning
n	Number of (scaled) phenotypic dimensions
λ	Variance in mutant phenotypes along each dimension
$m_{max}$	Maximum growth rate
$f\left(m\prime |m\right)$	Distribution of growth rates among mutants from a genotype with growth rate m (eq. 1)
U	Per genome mutation probability
$N_0$	Initial number of wild-type individuals
$m_0$	Wild-type growth rate
$p_0$	Probability a wild-type individual has descendants that rescue the population
P	Probability of rescue (eq. 2)
$p\left(m,\text{Λ}\left(m\right)\right)$	Probability a genotype with growth rate m, itself fated for extinction, has descendants that rescue the population (eq. 3)
$p_{est}\left(m\right)$	Probability a genotype with growth rate m establishes, i.e., rescues the population (eq. 4)
$\text{Λ}\left(m\right)$	Probability that an individual with growth rate m produces a mutant that has descendants that rescue the population
${\text{Λ}}_i\left(m\right)$	Probability that an individual with growth rate m produces a mutant that has descendants with $i-1$ additional mutations that rescue the population
${\text{Λ}}_2^i\left(m\right)$	Probability that an individual with growth rate m produces sufficiently subcritical $\left(i=”-”\right)$⁠, critical $\left(i=0\right)$⁠, or supercritical $\left(i=”+”\right)$ first-step mutants that eventually lead to two-step rescue (eq. 8)
ψ	$2\left(1-\sqrt{1-m/m_{max}}\right)$
${\psi }_0$	$2\left(1-\sqrt{1-m_0/m_{max}}\right)$
${\rho }_{max}$	$m_{max}/\lambda$
α	${\rho }_{max}{\psi }_0^2/4$

Fisher’s geometric model

We map genotype to phenotype to fitness using Fisher’s geometric model, originally introduced by Fisher (1930, pp.38–41) and reviewed by Tenaillon (2014). In this model, each genotype is characterized by a point in n-dimensional phenotypic space, $\overrightarrow{z}$⁠. We ignore environmental effects, and, thus, the phenotype is the breeding value. At any given time, there is a phenotype, $\overrightarrow{o}$⁠, that has maximum fitness, and fitness declines monotonically as phenotypes depart from $\overrightarrow{o}$⁠. We assume that n phenotypic axes can be chosen and scaled such that fitness is described by a multivariate Gaussian function, with variance one in each dimension, no covariance, and height $W_{\mathit{max}}$ (which can always be done when considering genotypes close enough to an nondegenerate optimum; Martin 2014). Thus the fitness of phenotype $\overrightarrow{z}$ is $W\left(\overrightarrow{z}\right)=W_{\mathit{max}}\hspace{0.17em}\text{exp}\left(-{\Vert \overrightarrow{z}-\overrightarrow{o}\Vert }^2/2\right)$⁠, where $\Vert \overrightarrow{z}-\overrightarrow{o}\Vert =\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(z_i-o_i\right)}^2}}$ is the Euclidean distance of $\overrightarrow{z}$ from the optimum, $\overrightarrow{o}$⁠. Here we are interested in absolute fitness; we take $\text{ln}\left[W\left(\overrightarrow{z}\right)\right]=m\left(\overrightarrow{z}\right)=m_{\mathit{max}}-{\Vert \overrightarrow{z}-\overrightarrow{o}\Vert }^2/2$ to be the continuous-time growth rate (m is for Malthusian fitness) of phenotype $\overrightarrow{z}$⁠. We ignore density- and frequency-dependence in $m\left(\overrightarrow{z}\right)$ for simplicity. The fitness effect, i.e., selection coefficient, of phenotype $z\prime$ relative to z in a continuous-time model is exactly $s=\text{ln}\left[W\left(z\prime \right)/W\left(z\right)\right]=m\left(z\prime \right)-m\left(z\right)$ (Martin and Lenormand 2015). This is approximately equal to the selection coefficient in discrete time $\left(W\left(z\prime \right)/W\left(z\right)-1\right)$ when selection is weak $\left(W\left(z\prime \right)-W\left(z\right)\ll 1\right)$⁠.

To make analytical progress, we use the isotropic version of Fisher’s geometric model, where mutations (in addition to selection) are assumed to be uncorrelated across the scaled traits. Universal pleiotropy is also assumed, so that each mutation affects all scaled phenotypes. In particular we use the “classic” form of Fisher’s geometric model (Harmand et al. 2017), where the probability density function of a mutant phenotype is multivariate normal, centered on the current phenotype, with variance λ in each dimension and no covariance. Using a probability density function of mutant phenotypes implies a continuum-of-alleles (Kimura 1965), i.e., phenotype is continuous and each mutation is unique. Mutations are assumed to be additive in phenotype, which induces epistasis in fitness (as well as dominance under diploid selection), as fitness is a nonlinear function of phenotype. We assume asexual reproduction, i.e., no recombination, which is appropriate for many cases of antimicrobial drug resistance and experimental evolution, while recognizing the value of expanding this work to sexual populations.

An obvious, and important, extension would be to relax the simplifying assumptions of isotropy and universal pleiotropy, which we leave for future work. Note that mild anisotropy yields the same bulk distribution of fitness effects as an isotropic model with fewer dimensions (Martin and Lenormand 2006), but this does not extend to the tails of the distribution. Therefore, whether anisotropy can be reduced to isotropy with fewer dimensions in the case of evolutionary rescue, where the tails are essential, is unknown. In the Discussion we briefly explore the effects of non-Gaussian distributions of mutant phenotypes.

Given this phenotype-to-fitness mapping and phenotypic distribution of new mutations, the distribution of fitness effects (and therefore growth rates) of new mutations can be derived exactly. Let m be the growth rate of some particular focal genotype, and $m\prime$ the growth rate of a mutant immediately derived from it. Then, let $s_o=m_{\mathit{max}}-m$ be the selective effect of a mutant with the optimum genotype and $s=m\prime -m$ the selective effect of the mutant with growth rate $m\prime$⁠. The probability density function of the selective effects of new mutations, s, is then given by Equation 3 in Martin and Lenormand (2015). Converting fitness effects to growth rate $\left(m\prime =s+m\right),$ the probability density function for mutant growth rate $m\prime$ from an ancestor with growth rate m is (cf. Equation 2 in Anciaux et al. 2018)

$$f\left(m\prime |m\right)=\frac{2}{\lambda }\hspace{0.17em}{f}_{{\chi }_n^2}\left(\frac{2\left(m_{\mathit{max}}-m\prime \right)}{\lambda },\frac{2\left(m_{\mathit{max}}-m\right)}{\lambda }\right),$$

(1)

where $f_{{\chi }_n^2}\left(x,c\right)$ is the probability density function over positive real numbers x of ${\chi }_n^2\left(c\right)$⁠, a noncentral chi-square deviate with n degrees of freedom and noncentrality $c>0$ (Equation 26.4.25 in Abramowitz and Stegun 1972).

Lifecycle

We are envisioning a scenario where $N_0$ wild-type individuals, each of which have phenotype ${\overrightarrow{z}}_0$⁠, experience an environmental change causing population decline, $m_0\equiv m\left({\overrightarrow{z}}_0\right)<0$⁠. Each generation, an individual with phenotype $\overrightarrow{z}$ produces a Poisson number of offspring, with mean $\text{exp}\left[m\left(\overrightarrow{z}\right)\right]$⁠, and dies. This process implicitly assumes no interaction between individuals, i.e., a branching process with density- and frequency-independent growth and fitness and no clonal interference. Each offspring mutates with probability U (we ignore the possibility of multiple simultaneous mutations within a single genome), and mutations are distributed as described above (see Fisher’s geometric model).

Simulation procedure

We ran individual-based simulations of the above process to compare with our numeric and analytic results. Except where noted, populations were considered rescued when there were $\ge 100$ individuals with positive growth rates (all other replicates went extinct). The most common genotype at the time of rescue was considered the rescue genotype, and the number of mutational steps to rescue was set as the number of mutations in that genotype.

Probability of rescue

Let $p_0$ be the probability that a given wild-type individual is “successful”, i.e., has descendants that rescue the population. The probability of rescue is then 1 minus the probability that none of the initial wild-type individuals are successful,

$$P=1-{\left(1-p_0\right)}^{N_0}\approx 1-\text{exp}\left(-N_0{p}_0\right),$$

(2)

where the approximation assumes small $p_0$ and large $N_0$⁠. What remains is to find $p_0$⁠.

Results

We start with a heuristic explanation of our main results before turning to more detailed derivations in the next section.

Rescue by multiple mutations

A characteristic pattern of evolutionary rescue is a “U”-shaped population size trajectory (e.g., Orr and Unckless 2014). This is the result of an exponentially declining wild-type genotype being replaced by an exponentially increasing mutant genotype. On a log scale, this population size trajectory becomes “V”-shaped (we denote it a “V-shaped log-trajectory”). On this scale, the population declines at a constant rate (producing a line with slope $m_0<0$⁠) until the growing mutant subpopulation becomes relatively common, at which point the population begins growing at a constant rate (a line with slope $m_1>0$⁠). This characteristic V-shaped log-trajectory is observed in many of our simulations where evolutionary rescue occurs (Figure 1A). Alternatively, when the wild type declines faster, and the mutation rate is larger, we sometimes see “U-shaped log-trajectories” (e.g., the red and blue replicates in Figure 2A). Here there are three phases instead of two; the initial rate of decline (a line with slope $m_0<0$⁠) is first reduced (transitioning to a line with slope $m_1<0$⁠) before the population begins growing (a line with slope $m_2>0$⁠).

Figure 1

Typical dynamics with a relatively slow wild-type decline and a small mutation rate (⁠$m_0=-0.1$⁠, $U={10}^{-4}$⁠). (A) Population size trajectories on a log scale. Each line is a unique replicate simulation (100 replicates). Replicates that went extinct are gray, replicates that were rescued are in color (and are roughly V-shaped). (B) The number of individuals with a given derived allele, again on a log scale, for the yellow replicate in (A). The number of individuals without any derived alleles (wild type) is shown in gray, the rescue mutation is shown in yellow, and all other mutations are shown in black. Other parameters: $n=4$⁠, $\lambda =0.005$⁠, $m_{max}=0.5$⁠.

Figure 2

Typical dynamics with a relatively fast wild-type decline and a large mutation rate (⁠$m_0=-0.3$⁠, $U={10}^{-2}$⁠). (A) Population size trajectories on a log scale. Each line is a unique replicate simulation (500 replicates). Replicates that went extinct are gray, replicates that were rescued are in color. Note that the blue and red replicates are cases of two-step rescue (and roughly U-shaped), while the yellow replicate is one-step rescue (and therefore V-shaped). (B) The number of individuals with a given derived allele, again on a log scale, for the red replicate in (A). The number of individuals without any derived alleles (wild type) is shown in gray, the rescue mutations are shown in red, and all other mutations in black. Here, a single mutant with growth rate <0 arises early and outlives the wild type (solid red). A second mutation then arises on that background (dashed red), making a double mutant with a growth rate >0 that rescues the population. Other parameters: $n=4$⁠, $\lambda =0.005$⁠, $m_{max}=0.5$⁠.

As expected, V-shaped log-trajectories are the result of a single mutation creating a genotype with a positive growth rate that escapes loss when rare and rescues the population (Figure 1B), i.e., one-step rescue. U-shaped log-trajectories, on the other hand, occur when a single mutation creates a genotype with a negative (or potentially very small positive) growth rate, itself doomed to extinction, which outpersists the wild type, and gives rise to a double-mutant genotype that rescues the population (Figure 2B), i.e., two-step rescue. These two types of rescue comprise the overwhelming majority of rescue events observed in our simulations, across a wide range of wild-type decline rates (e.g., Figure 3).

Figure 3

The probability of evolutionary rescue as a function of initial maladaptation. Shown are the probabilities of 1-, 2-, 3-, and 4-step rescue (using Equations 2–7), as well as the probability of rescue by up to four mutational steps (“total”, using $\text{Λ}\left(m_0\right)={\displaystyle {\sum }_{i=1}^4{\text{Λ}}_i\left(m_0\right)}$⁠). Circles are individual-based simulation results (ranging from ${10}^5$ to ${10}^6$ replicates per wild-type growth rate). Populations were considered rescued when there were ≥ 1000 individuals with positive growth rate. Open circles denote the fraction of simulations where the rescue genotype (see Simulation procedure) had a given number of mutations, and closed circles are the sum of these fractions. Parameters: $N_0={10}^4$⁠, $U=2\times {10}^{-3}$⁠, $n=4$⁠, $\lambda =0.005$⁠, $m_{max}=0.5$⁠.

In the text, we focus on low to moderate mutation rates affecting growth rate. With sufficiently high mutation rates, rescue by three or more mutations comes to dominate (Figure S1). It has recently been suggested that when the mutation rate, U, is substantially less than a critical value, $U_C=\lambda n^2/4$⁠, we are in a “strong selection, weak mutation” regime, where selection is strong enough relative to mutation that essentially all mutations arise on a wild-type background (Martin and Roques 2016), consistent with the House of Cards approximation (Turelli 1984, 1985). Thus, in this regime, rescue tends to occur by a single mutation of large effect (Anciaux et al. 2018). In the other extreme, when $U\gg U_C$⁠, we are in a “weak selection, strong mutation” regime, where selection is weak enough relative to mutation that many cosegregating mutations are present within each genome, creating a multivariate normal phenotypic distribution (Martin and Roques 2016), consistent with the Gaussian approximation (Kimura 1965; Lande 1980). Thus, in this regime, rescue tends to occur by many mutations of small effect (Anciaux et al. 2019). As shown in Figure 3 (where $U=U_C/10$⁠) and Figure S1 (where $U_C=0.02$⁠), rescue by a small number of mutations (but >1) can become commonplace in the transition zone (where U is neither much smaller or much larger than $U_C$⁠), where there are often a considerable number of cosegregating mutations (e.g., Figure 2B, where $U=U_C/2$⁠).

The probability of k-step rescue

Approximations for the probability of one-step rescue under the strong selection, weak mutation regime were derived by Anciaux et al. (2018). Here, we extend this study by exploring the contribution of k-step rescue, deriving approximations for the probability of such events, as well as dissecting the genetic basis of both one- and two-step rescue in terms of the distribution of fitness effects of rescue genotypes and their component mutations.

Although requiring a sufficiently beneficial mutation to arise on a rare mutant genotype doomed to extinction, multistep evolutionary rescue can be the dominant form of rescue when the wild type is sufficiently maladapted (Figure 3 and Figure S1). Indeed, on this fitness landscape, the probability of producing a rescue genotype in one mutational step mutant drops very sharply with maladaptation (Anciaux et al. 2018); the probability of multi-step rescue declines more slowly as mutants with intermediate growth rates can be a “springboard”—albeit not always a very bouncy one—from which rescue mutants are produced. These intermediates contribute more as mutation rates and the decline rate of the wild-type increase (Figure 3 and Figure S1), the former because double mutants become more likely and the latter because the springboard becomes more necessary. With a large enough number of wild-type individuals or a high enough mutation rate (Figure S1), multistep rescue can not only be more likely than one-step, but also very likely in an absolute sense.

Classifying two-step rescue regimes

Two-step rescue can occur through first-step mutants with a wide range of growth rates. As shown below (see Approximating the probability of two-step rescue), these first-step mutants can be divided into three regimes: “sufficiently subcritical”, “sufficiently critical”, and “sufficiently supercritical” (we will often drop “sufficiently” for brevity; Figure 4). Sufficiently critical first-step mutants are defined by having growth rates close enough to zero that the most likely way for such a mutation to lead to two-step rescue is for it to persist for such an unusually long period of time, and accordingly grow to such an unusually large subpopulation size, that it will almost certainly produce successful double mutants. Sufficiently subcritical first-step mutants are then defined by having growth rates that are negative enough to almost certainly prevent such long persistence times. Instead, these mutations tend to persist for an expected number of generations, proportional to the inverse of their growth rate $\left(1/\left|m\right|\right)$⁠, while maintaining relatively small subpopulation sizes (on the order of one individual per generation). Mutations conferring a positive growth rate can also go extinct, and thus can also act as springboards to rescue. Conditioned on extinction, supercritical mutations behave like subcritical mutations with a growth rate of the same absolute value (Maruyama and Kimura 1974). Sufficiently supercritical first-step mutants are therefore defined analogously to subcritical first-step mutants, having positive (rather than negative) growth rates that are large enough to prevent sufficiently long persistence times once conditioned on extinction. Despite having similar extinction trajectories as subcritical mutations, “doomed” supercritical mutations arise less frequently by mutation from the wild type but mutate to rescue genotypes at a higher rate. Overall, they too can contribute substantially to rescue. Note that supercritical two-step rescue is not one-step rescue with subsequent adaptation, as we condition on the first-step mutation going extinct in the absence of the second mutation. However, empirically, it will be impossible to tell if the first-step mutation was indeed doomed to extinction if it is found to have a positive growth rate in the selective environment.

Figure 4

One- and two-step genetic paths to evolutionary rescue. Here we show an $n=2$ dimensional phenotypic landscape. Continuous-time (Malthusian) growth rate (m) declines quadratically from the center, becoming negative outside the thick black line. The gray zone indicates where growth rates are “sufficiently critical” (see text for details). Blue circles show wild-type phenotypes, red circles show the phenotypes of intermediate first-step mutations, and yellow circles show the phenotypes of rescue genotypes.

The relative contribution of each regime changes with both the initial degree of maladaptation and the mutation rate (Figure 5 and Figure S2). When the wild type is very maladapted (relative to mutational variance), most two-step rescue events occur through subcritical first-step mutants (Figure 5A), which arise at a higher rate than critical or supercritical mutants, and, yet, persist longer than the wild type. When the wild type is less maladapted, however, critical and supercritical mutations become increasingly likely to arise and contribute to two-step rescue, both due to their closer proximity to the wild type in phenotypic space as well as the slower decline of the wild type increasing the cumulative number of mutations that occur. The mutation rate also plays an interesting role in determining the relative contributions of each regime (Figure 5B and Figure S2). When mutations are rare, only first-step mutations that are very nearly neutral $\left(m\sim 0\right)$ will persist long enough to give rise to a two-step rescue mutation. As the mutation rate increases, however, the range of first-step mutant growth rates that can persist long enough to lead to two-step rescue widens because fewer individuals carrying the first-step mutation are needed before a successful double mutant arises.

Figure 5

The relative contribution of sufficiently subcritical, critical, and supercritical single mutants to two-step rescue. The curves are drawn using Equations 10–14 (Equation 12 is used for $m_0<0.2$ while Equation 13 is used for $m_0>0.2$⁠). The dots are numerical evaluations of Equation 8. Parameters: $n=4$⁠, $\lambda =0.005$⁠, $m_{\text{max}}=0.5$⁠, (A) $U={10}^{-3}$⁠, (B) $m_0=-0.1$⁠.

The distribution of fitness effects among rescue mutations

Mutants causing one-step rescue have growth rates that cluster around small positive values (⁠$m\mathrm{\gtrsim }0$⁠; blue curves in Figure 6). Consequently, the distribution of fitness effects (DFE) among these rescue mutants is shifted to the right relative to mutations that establish in a population of constant size (compare solid blue and gray curves in Figure 6), with a DFE beginning at $s=m-m_0\ge -m_0>0$ rather than $s=0$ (Kimura 1983). As a result of this increased threshold, the one-step rescue DFE has a smaller variance than both the DFE of random mutations, and the DFE of mutations that establish in a constant population (compare blue and gray curves in Figure 6). Further, while the variance in the DFE of random mutations and of those that establish in a population of constant size increases slightly with initial maladaptation (due to the curvature of the phenotype-to-fitness function), the variance in the one-step rescue DFE decreases substantially (compare panels in Figure 6), as rescue becomes restricted to a rapidly decreasing proportion of the available mutants.

Figure 6

The distribution of growth rates among rescue genotypes under one-step (blue; Equation 15 solid and 16 dashed) and two-step (red; Equation 17) rescue for three different levels of initial maladaptation. For comparison, the distribution of random mutations (dashed gray; Equation 1) and the distribution of beneficial mutations that establish in a population of constant size (solid gray; Equation 1 times Equation 4 and normalized) are shown. Intervals (horizontal lines) indicate the size of the most common fitness effect $\left(s=m_0-m\right)$ in a population of constant size (gray) and in one-step rescue (blue). The histograms show the distribution of growth rates among rescue genotypes observed across (A) ${10}^4$⁠, (B) ${10}^5$⁠, and (C) ${10}^6$ simulated replicates. Other parameters: $N_0={10}^4$⁠, $U=2\times {10}^{-3}$⁠, $n=4$⁠, $\lambda =0.005$⁠, $m_{\text{max}}=0.5$⁠.

The DFE of genotypes that cause two-step rescue (the combined effect of two mutations) is also clustered at small positive growth rates, but it has a variance that is less affected by the rate of wild-type decline (red curves in Figure 6). This is because double-mutant rescue genotypes are created via first-step mutant genotypes that have larger growth rates than the wild type (i.e., are closer to the optimum), allowing them to create double mutants with a larger range of positive growth rates.

Finally, we can also look at the distribution of growth rates among first-step mutations that lead to two-step rescue, i.e., “springboard mutants” (Figure 7 and Figure S2). Here there are two main factors to consider: (1) the probability that a mutation with a given growth rate arises on the wild-type background but does not by itself rescue the population, and (2) the probability that such a mutation persists long enough for a sufficiently beneficial second mutation to arise on that same background and together rescue the population. Subcritical mutations conferring growth rates closer to zero persist longer but are less likely to arise from the wild type, creating a trade-off between mutational input and the probability of rescue that can lead to a wide distribution of contributing subcritical growth rates (blue shading in Figure 7). In contrast, supercritical mutations with growth rates nearer to zero are more likely arise by mutation, to go extinct in the absence of further mutation, and to persist for longer once conditioned on extinction, together creating a relatively narrow distribution of contributing supercritical growth rates (yellow shading in Figure 7). As explained above, increasing the rate of wild-type decline (or decreasing the rate of mutation) increases the contribution of subcritical first-step mutants, and the importance of mutational input, lowering the mode and increasing the variance of the first-step DFE (compare panels in Figure 7).

Figure 7

The distribution of growth rates among first-step mutations that lead to two-step rescue (black; Equation 18) for three different levels of initial maladaptation. Shading represents our sufficiently subcritical approximation (blue; replacing $p\left(m,{\text{Λ}}_1\left(m\right)\right)$ with ${\text{Λ}}_1\left(m\right)/\left|m\right|$ in the numerator of Equation 18), our sufficiently critical approximation (red; using $Uf\left(0|m_0\right)\sqrt{2{\text{Λ}}_1\left(0\right)}$ as the numerator in Equation 18), and our sufficiently supercritical approximation (yellow; replacing $p\left(m,{\text{Λ}}_1\left(m\right)\right)$ with ${\text{Λ}}_1\left(m\right)/\left|m\right|$ in the numerator of Equation 18). The histograms show the distribution of growth rates among first-step mutations in rescue genotypes with two mutations observed across (A and B) ${10}^5$ or (C) ${10}^6$ simulated replicates. We hypothesize that the overabundance of supercriticals (especially in panel (A)) is likely due to us sampling only the most common rescue genotype in each replicate, which is not necessarily the first genotype that rescues. See Figure 6 for additional details.

Note that, given two-step rescue, the growth rate of both the first-step and second-step mutation may be negative when considered by themselves in the wild-type background. This potentially obscures empirical detection of the individual mutations involved in evolutionary rescue.

Mathematical Analysis

The probability of k-step rescue

Generic expressions for the probability of one- and two-step rescue were given by Martin et al. (2013), using a diffusion approximation of the underlying demographics. The key result that we will use is the probability that a single copy of a genotype with growth rate m, itself fated for extinction, but which produces rescue mutants at rate $\text{Λ}\left(m\right)$⁠, rescues the population (equation S1.5 in Martin et al. 2013). With our lifecycle this is (cf. equation A.3 in Iwasa et al. 2004a)

$$p\left(m,\text{Λ}\left(m\right)\right)=1-\text{exp}\left[\left|m\right|\left(1-\sqrt{1+\frac{2\text{Λ}\left(m\right)}{m^2}}\right)\right].$$

(3)

We can therefore use $p_0=p\left(m_0,\text{Λ}\left(m_0\right)\right)$ as the probability that a wild-type individual has descendants that rescue the population and what remains in calculating the total probability of rescue (Equation 2) is $\text{Λ}\left(m_0\right)$⁠. We break this down by letting ${\text{Λ}}_i\left(m\right)$ be the rate at which rescue genotypes with i mutations are created; the total probability of rescue is then given by Equation 2 with $p_0=p\left(m_0,{\displaystyle {\sum }_{i=1}^{\mathrm{\infty }}{\text{Λ}}_i\left(m_0\right)}\right)$⁠.

In one-step rescue, ${\text{Λ}}_1\left(m_0\right)$ is just the rate of production of rescue mutants directly from a wild-type genotype. This is the probability that a wild type gives rise to a mutant with growth rate m (given by $Uf\left(m|m_0\right)$⁠) times the probability that a genotype with growth rate m establishes. Again approximating our discrete time process with a diffusion process, the probability that a lineage with growth rate $m\ll 1$ establishes, ignoring further mutation, is (e.g., Martin et al. 2013)

$$p_{est}\left(m\right)\approx \{\begin{array}{ll}0\hfill & m\le 0\hfill \\ 1-\text{exp}\left(-2m\right)\hfill & m>0\hfill \end{array}.$$

(4)

This reduces to the $2\left(s+m_0\right)$ result in Otto and Whitlock (1997), when $m=s+m_0$ is small, which further reduces to $2s$ in a population of constant size, where $m_0=0$ (Haldane 1927). Using this, the rate of one-step rescue is

$${\text{Λ}}_1\left(m_0\right)=U{\displaystyle {\int }_0^{m_{max}}f\left(m|m_0\right)p_{est}\left(m\right)\text{d}m.}$$

(5)

Taking the first order approximation of $p\left(m_0,{\text{Λ}}_1\left(m_0\right)\right)$ with ${\text{Λ}}_1\left(m_0\right)/m_0^2$ small gives the probability of one-step rescue (Equation 5 of Anciaux et al. 2018), which effectively assumes deterministic wild-type decline. For completeness we rederive their closed-form approximation in File S2 (and give the results in the Appendix, see Approximating the probability of one-step rescue).

The probability of two-step rescue is only slightly more complicated. Here ${\text{Λ}}_2\left(m_0\right)$ is the probability that a mutation arising on the wild-type background creates a genotype that is also fated for extinction but persists long enough for a second mutation to arise on this mutant background, creating a double-mutant genotype that rescues the population. We therefore have

$${\text{Λ}}_2\left(m_0\right)=U{\displaystyle {\int }_{-\mathrm{\infty }}^{m_{max}}f\left(m|m_0\right)\left[1-p_{\text{est}}\left(m\right)\right]p\left(m,{\text{Λ}}_1\left(m\right)\right)\text{d}m.}$$

(6)

Following this logic, we can retrieve the probability of k-step rescue, for arbitrary $k\ge 2,$ using the recursion

$${\text{Λ}}_k\left(m_0\right)=U{\displaystyle {\int }_{-\mathrm{\infty }}^{m_{max}}f\left(m|m_0\right)\left[1-p_{est}\left(m\right)\right]}$$

(7)

$$p\left(m,{\text{Λ}}_{k-1}\left(m\right)\right)\text{d}m,$$

with the initial condition given by Equation 5. The probability of 1-, 2-, 3-, and 4-step rescue is plotted in Figure 3.

Approximating the probability of two-step rescue

The probability of two-step rescue is given by Equation 2, with $p_0=p\left(m_0,{\text{Λ}}_2\left(m_0\right)\right)$ (Equations 3–6). We next develop some intuition by approximating this for different classes of single mutants.

First, note that, when the growth rate of a first-step mutation is close enough to zero such that $m^2\ll {\text{Λ}}_1\left(m\right)$⁠, we can approximate the probability that such a genotype leads to rescue before itself going extinct, $p\left(m,{\text{Λ}}_1\left(m\right)\right)$⁠, using a Taylor series, as $\sqrt{2{\text{Λ}}_1\left(m\right)}$ [cf. equation A.4b in Iwasa et al. (2004a), see also File S2]. We can also derive this result heuristically by considering the probability that a lineage will persist long enough that it will incur a secondary rescue mutation. As shown in the Appendix (see Mutant lineage dynamics), while $t<1/\left|m\right|,$ a mutant lineage with growth rate m that is destined for extinction persists for t generations with probability $\sim 2/t$ (Equation 21), and in generation t since it has arisen has $\sim t/2$ individuals (Equation 22). Thus, while $T<1/\left|m\right|,$ a mutant lineage that persists for T generations will have produced a cumulative number $\sim T^2/4$ individuals. Such lineages will then lead to two-step rescue with probability $\sim {\text{Λ}}_1\left(m\right)T^2/4$ until this approaches 1, near $T=2/\sqrt{{\text{Λ}}_1\left(m\right)}$⁠. Since the probability of rescue increases like $T^2$ while the probability of persisting to time T declines only like $1/T$⁠, most rescue events will be the result of rare long-lived single mutant genotypes. Considering only the most long-lived genotypes, the probability that a first-step mutation leads to rescue is then the probability that it survives long enough to almost surely rescue, i.e., for $T\sim 2/\sqrt{{\text{Λ}}_1\left(m\right)}$ generations. Since the probability of such a long-lived lineage is $2/T\sim \sqrt{{\text{Λ}}_1\left(m\right)}$⁠, this heuristic result agrees with our Taylor series approximation of Equation 5. Thus, for first-step mutants with growth rates satisfying $2/\sqrt{{\text{Λ}}_1\left(m\right)}<1/\left|m\right|$⁠, implying $m^2\ll {\text{Λ}}_1\left(m\right)$⁠, with probability $\sim \sqrt{{\text{Λ}}_1\left(m\right)}$ persistence is long enough to almost certainly ensure rescue. This same reasoning has been used to explain why the probability that a neutral mutation segregates long enough to produce a second mutation is $\sim \sqrt{U}$ in a population of constant size (Weissman et al. 2009).

At the other extreme, when the growth rate of a first-step mutation is far enough from zero such that $m^2\gg {\text{Λ}}_1\left(m\right)$⁠, we can approximate $p\left(m,{\text{Λ}}_1\left(m\right)\right)$⁠, again using a Taylor series, with ${\text{Λ}}_1\left(m\right)/\left|m\right|$ [cf. equation A.4c in Iwasa et al. (2004a), see also File S2]. Conditioned on extinction such genotypes cannot persist long enough to almost surely lead to two-step rescue. Instead, we expect such mutations to persist for, at most, $\sim 1/\left|m\right|$ generations (Equation 21) with a lineage size of $\sim 1$ individual per generation (Equation 22), and thus produce a cumulative total of $\sim 1/\left|m\right|$ individuals. The probability of two-step rescue from such a first-step mutation is therefore ${\text{Λ}}_1\left(m\right)/\left|m\right|$⁠, and, again, this heuristic argument matches our Taylor series approach. This same reasoning explains why a rare mutant genotype with selection coefficient $\left|s\right|\gg 0$ in a constant population size model is expected to have a cumulative number of $\sim 1/\left|s\right|$ descendants, given it eventually goes extinct (Weissman et al. 2009).

The transitions between these two regimes occur when ${\text{Λ}}_1\left(m\right)/\left|m\right|=\sqrt{2{\text{Λ}}_1\left(m\right)}$⁠, i.e., when $\left|m\right|=\sqrt{{\text{Λ}}_1\left(m\right)/2}$⁠. We call single mutants with growth rates $m<-\sqrt{{\text{Λ}}_1\left(m\right)/2}$ “sufficiently subcritical”, those with $\left|m\right|<\sqrt{{\text{Λ}}_1\left(m\right)/2}$ “sufficiently critical”, and those with $m>\sqrt{{\text{Λ}}_1\left(m\right)/2}$ “sufficiently supercritical”. Given that U, and, thus, ${\text{Λ}}_1\left(m\right)$ will generally be small, m will also be small at these transition points, meaning we can approximate the transition points as $m^{\text{*}}=\sqrt{{\text{Λ}}_1\left(0\right)/2}$ and $-m^{\text{*}}$⁠. We then have an approximation for the rate of two-step rescue,

$$\begin{array}{c}{\text{Λ}}_2\left(m_0\right)={\text{Λ}}_2^{(-)}\left(m_0\right)+{\text{Λ}}_2^{\left(0\right)}\left(m_0\right)+{\text{Λ}}_2^{(+)}\left(m_0\right)\\ {\text{Λ}}_2^{(-)}\left(m_0\right)=U{\displaystyle {\int }_{-\mathrm{\infty }}^{-m^{\text{*}}}f(m|m_0)\frac{{\text{Λ}}_1\left(m\right)}{\left|m\right|}\text{d}m}\\ {\text{Λ}}_2^{\left(0\right)}\left(m_0\right)=U{\displaystyle {\int }_{-m^{\text{*}}}^{m^{\text{*}}}f\left(m|m_0\right)\left[1-p_{est}\left(m\right)\right]\sqrt{2{\text{Λ}}_1\left(m\right)}\text{d}m}\\ {\text{Λ}}_2^{(+)}\left(m_0\right)=U{\displaystyle {\int }_{m^{\text{*}}}^{m_{max}}f\left(m|m_0\right)\left[1-p_{est}\left(m\right)\right]\frac{{\text{Λ}}_1\left(m\right)}{\left|m\right|}\text{d}m}\end{array}$$

(8)

where ${\text{Λ}}_2^{\left(i\right)}\left(m_0\right)$ is the rate of two-step rescue through sufficiently subcritical first-step mutants $\left(i=“-”\right)$⁠, sufficiently critical first-step mutants $\left(i=0\right)$⁠, or sufficiently supercritical first-step mutants $\left(i=“+”\right)$⁠. A schematic depicting the one- and two-step genetic paths to rescue is given in Figure 4.

Closed-form approximation for critical two-step rescue

When U is small $m^{\text{*}}$ is also small, allowing us to use $m=0$ within the integrand of ${\text{Λ}}_2^{\left(0\right)}\left(m_0\right)$⁠, which spans a range, $\left[-m^{\text{*}},m^{\text{*}}\right]$⁠, of width $2m^{\text{*}}\approx \sqrt{2{\text{Λ}}_1\left(0\right)}$⁠, giving

$$\begin{array}{l}{\text{Λ}}_2^{\left(0\right)}\left(m_0\right)\approx Uf\left(0|m_0\right)\sqrt{2{\text{Λ}}_1\left(0\right)}2m^{\text{*}}\\ =2Uf\left(0|m_0\right){\text{Λ}}_1\left(0\right).\end{array}$$

(9)

We can then approximate ${\text{Λ}}_1\left(m\right)$ with ${\stackrel{\sim }{\text{Λ}}}_1\left(m\right)$ (Equation 19), and take $m\to 0$ (Equation 20), giving a closed form approximation for the rate of two-step rescue through critical single mutants in Fisher’s geometric model,

$${\text{Λ}}_2^{\left(0\right)}\left(m_0\right)\approx 4U^{2}f\left(0|m_0\right)\sqrt{m_{max}\lambda /\pi }.$$

(10)

This well approximates numerical integration of ${\text{Λ}}_2^{\left(0\right)}\left(m_0\right)$ (Equation 8; see Figure 5 and File S2). In general, it will perform better when the critical zone, and thus $U\sqrt{m_{max}\lambda }$⁠, becomes smaller.

To get a better understanding of how the rate of two-step critical rescue depends on the underlying parameters of Fisher’s geometric model, we approximate $f\left(m|m_0\right)$⁠, assuming that the distance from the wild type to the optimal phenotype is large relative to the distribution of mutations (i.e., ${\rho }_{\mathit{max}}=m_{\mathit{max}}/\lambda$ is large), and convert this to a distribution over $\psi =2\left(1-\sqrt{1-m/m_{\mathit{max}}}\right)$⁠, a convenient rescaling (for details see File S2 and Anciaux et al. 2018). Evaluating this at $m=0$ gives

$${\text{Λ}}_2^{\left(0\right)}\left(m_0\right)\approx U^2{\left(1-{\psi }_0/2\right)}^{\left(1-n\right)/2}{e}^{-\alpha }\frac{2}{\pi },$$

(11)

where ${\psi }_0=2\left(1-\sqrt{1-m_0/m_{\mathit{max}}}\right)<0$ and $\alpha ={\rho }_{\mathit{max}}{\psi }_0^2/4$⁠.

Closed-form approximations for noncritical two-step rescue

We can also approximate ${\text{Λ}}_1\left(m\right)$ in ${\text{Λ}}_2^{(-)}\left(m_0\right)$ and ${\text{Λ}}_2^{(+)}\left(m_0\right)$ with ${\stackrel{\sim }{\text{Λ}}}_1\left(m\right)$ (Equation 19), leaving us with just one integral over the growth rates of the first-step mutations. We then replace $f\left(m|m_0\right)$ with its approximate distribution over ψ as above.

In the case of subcritical rescue, we can then make two contrasting approximations (see File S2 for details). First, when the ψ (and thus m) that contribute most are close enough to zero (meaning maladaptation is not too large relative to mutational variance), and we ignore mutations that are less fit than the wild type, we find the rate of subcritical two-step rescue is roughly

$${\text{Λ}}_2^{(-)}\left(m_0\right)\approx U^2\frac{{\left(1-{\psi }_0/2\right)}^{1-n}}{1-{\psi }_0/4}{e}^{-\alpha }\frac{\text{log}\left({\psi }_0/{\psi }_{-}^{\text{*}}\right)}{\pi },$$

(12)

where ${\psi }_{-}^{\text{*}}=2\left(1-\sqrt{1+{\stackrel{\sim }{m}}^{\text{*}}/m_{\mathit{max}}}\right)<0$ and ${\stackrel{\sim }{m}}^{\text{*}}=\sqrt{{\stackrel{\sim }{\text{Λ}}}_1\left(0\right)/2}$ (Equation 20). Second, when the mutational variance, λ, is very small relative to maladaptation, implying that mutants far from $m=0$ substantially contribute, we find the rate of subcritical two-step rescue to be nearly

$${\text{Λ}}_2^{(-)}\left(m_0\right)\approx -U^2\frac{{\left(1-{\psi }_0/2\right)}^{1-n}}{1-{\psi }_0/4}{\left(e^{-\alpha }\frac{1}{{\left(\alpha /2\right)}^3\pi }\right)}^{1/2}.$$

(13)

These two approximations do well compared with numerical integration of ${\text{Λ}}_2^{(-)}\left(m_0\right)$ (Equation 8; see Figure 5 and File S2). As expected, we find that Equation 13 does better under fast wild-type decline, while Equation 12 does better when the wild type is declining more slowly.

For supercritical two-step rescue, only first-step mutants with growth rates near $m^{\text{*}}$ will contribute (larger m will rescue themselves and are also less likely to arise by mutation), and so we can capture the entire distribution with a small m approximation (following the same approach that led to Equation 12). As shown in File S2, this approximation works well for sufficiently small first-step mutant growth rates, $\psi <\sqrt{2/{\rho }_{\mathit{max}}}$⁠, beyond which the rate of two-step rescue through such first-step mutants falls off very quickly due to a lack of mutational input. Thus, considering only supercritical single mutants with scaled growth rate $m<\sqrt{2/{\rho }_{\mathit{max}}}$⁠, our approximation is

$${\text{Λ}}_2^{(+)}\left(m_0\right)\approx U^2\frac{{\left(1-{\psi }_0/2\right)}^{1-n}}{1-{\psi }_0/4}{e}^{-\alpha }\frac{\text{log}\left({\psi }_{\mathit{max}}/{\psi }_{+}^{\text{*}}\right)}{\pi },$$

(14)

with ${\psi }_{+}^{\text{*}}=2\left(1-\sqrt{1-{\stackrel{\sim }{m}}^{\text{*}}/m_{\mathit{max}}}\right)$ and ${\psi }_{\mathit{max}}=\sqrt{2/{\rho }_{\mathit{max}}}$⁠. This approximation tends to provide a slight overestimate of ${\text{Λ}}_2^{(+)}\left(m_0\right)$ (Equation 8; see Figure 5 and File S2).

Comparing two-step regimes:

These rough but simple closed-form approximations (Equations 11–14) show that, while the contribution of critical mutants to two-step rescue scales with $U^2$⁠, the contribution of noncritical single mutants scales at a rate less than $U^2$ (Figure 5B), due to a decrease in ${\psi }_{-}^{\text{*}}$ (decreasing the range of subcritical mutants) and an increase in ${\psi }_{+}^{\text{*}}$ (decreasing the range of supercritical mutants) with U. This difference in scaling with U is stronger when the wild type is not very maladapted relative to the mutational variance, i.e., when Equation 12 is the better approximation for subcritical rescue. The approximations also show that when initial maladaptation is small, the ratio of supercritical to subcritical contributions (Equation 12 divided by Equation 14) depends primarily on the range of growth rates included in each regime, while with larger initial maladaptation this ratio (Equation 13 divided by Equation 14) begins to depend more strongly on initial maladaptation and mutational variance (α). The effect of maladaptation and mutation rate on the relative contributions of each regime is shown in Figure 5.

The distribution of growth rates among rescue genotypes

We next explore the distribution of growth rates among rescue genotypes, i.e., the distribution of growth rates that we expect to observe among the survivors across many replicates.

We begin with one-step rescue. The rate of one-step rescue by genotypes with growth rate m is simply $Uf\left(m|m_0\right)p_{\text{est}}\left(m\right)$⁠. Dividing this by the rate of one-step rescue through any m (Equation 5) gives the distribution of growth rates among the survivors

$$g_1\left(m\right)=\frac{Uf\left(m|m_0\right)p_{\text{est}}\left(m\right)}{{\text{Λ}}_1\left(m_0\right)},$$

(15)

where the mutation rate, U, cancels out. This distribution is shown in blue in Figure 6. The distribution has a mode at small but positive m as a result of two conflicting processes: smaller growth rates are more likely to arise from a declining wild type but larger growth rates are more likely to establish given they arise. As the rate of wild-type decline increases, the former process exerts more influence, causing the mode to move toward zero and reducing the variance.

We can also give a simple, nearly closed-form approximation here using the same approach taken to reach Equation 19. On the ψ scale, the distribution of effects among one-step rescue mutations is

$${\stackrel{\sim }{g}}_1\left(\psi \right)=\frac{\text{exp}\left(\alpha \right)\sqrt{\alpha {\rho }_{\mathit{max}}}}{\left[\text{exp}\left(\alpha \right)\sqrt{\pi \alpha }\text{Erfc}\left(\sqrt{\alpha }\right)-1\right]{\psi }_0}{e}^{-{\rho }_{\mathit{max}}{\left(\psi -{\psi }_0\right)}^2/4}\psi ,$$

(16)

implying the ψ are distributed like a normal truncated below $\psi =0$ and weighted by ψ. This often provides a very good approximation (see dashed blue curves in Figure 6).

In two-step rescue, the rate of rescue by double mutants with growth rate $m_2$ is given by Equation 6, with ${\text{Λ}}_1\left(m\right)$ replaced by $Uf\left(m_2|m\right)p_{est}\left(m_2\right)$⁠. Normalizing gives the distribution of growth rates among the double-mutant genotypes that rescue the population

$$g_2\left(m_2\right)\approx \frac{A\left(m_2\right)}{{\displaystyle {\int }_0^{m_{max}}A\left(m_2\right)\text{d}{m}_2}}$$

(17)

$$A\left(m_2\right)={\displaystyle {\int }_{-\mathrm{\infty }}^{m_{max}}f\left(m|m_0\right)\left[1-p_{est}\left(m\right)\right]}$$

$$p\left(m,Uf\left(m_2|m\right)p_{est}\left(m_2\right)\right)\text{d}m.$$

This distribution, $g_2\left(m\right)$⁠, is shown in red in Figure 6. Because the first-step mutants contributing to two-step rescue tend to be nearer the optimum than the wild type, this allows them to produce double-mutant rescue genotypes with higher growth rates than in one-step rescue (as seen by comparing the mode between blue and red curves in Figure 6). The fact that these first-step mutants are closer to the optimum also allows for a greater variance in the growth rates of rescue genotypes than in one-step rescue. Thus, the two-step distribution maintains a more similar mode and variance across wild-type decline rates than the one-step distribution. Note that because $g_2\left(m_2\right)$ depends on U, the buffering effect of first-step mutants depends on the mutation rate (see The distribution of growth rates among rescue intermediates below for more discussion).

The distribution of growth rates among rescue intermediates

Finally, our analyses above readily allow us to explore the distribution of first-step mutant growth rates that contribute to two-step rescue. Analogously to Equation 15, we drop the integral in ${\text{Λ}}_2\left(m_0\right)$ (Equation 6) and normalize, giving

$$h\left(m\right)=\frac{Uf\left(m|m_0\right)\left[1-p_{\text{est}}\left(m\right)\right]p\left(m,{\text{Λ}}_1\left(m\right)\right)}{{\text{Λ}}_2\left(m_0\right)},$$

(18)

where the first U cancels but the U within ${\text{Λ}}_1\left(m\right)$ does not. This distribution is shown in black in Figure 7. At slow wild-type decline rates, the overwhelming majority of two-step rescue events arise from first-step mutants with growth rates near 0. As indicated by Equation 8, the contribution of first-step mutants with growth rate m declines as $\sim 1/\left|m\right|$ as m departs from zero, due to shorter persistence times given eventual extinction. As wild-type growth rate declines, the relative importance of mutational input, $f\left(m|m_0\right)$⁠, grows, causing the distribution to flatten, and first-step mutants with substantially negative growth rates begin to contribute (compare panels in Figure 7; see also Figure 5A). Decreasing the mutation rate disproportionately increases the contribution of first-step mutants with growth rates near zero (while simultaneously shrinking the range of growth rates that are sufficiently critical; Figure 5B), making the distribution of first-step mutant growth rates contributing to two-step rescue more sharply peaked around $m=0$ (Figure S2). Correspondingly, with a higher mutation rate, a greater proportion of the contributing single mutants have substantially negative growth rates.

Discussion

Here, we have explored the probability and genetic basis of evolutionary rescue by multiple mutations on a simple fitness landscape. We find that rescue by multiple mutations can be the most likely path to persistence under high mutation rates, or when the population is initially very maladapted. Under these scenarios, intermediate genotypes that are declining less quickly provide a “springboard” from which rescue genotypes emerge. In two-step rescue, these springboard single mutants come from one of three regimes: those that have growth rates near enough to zero (“sufficiently critical”) that rescue is most likely when a mutation persists for an unusually long period of time, and grows to an unusually large subpopulation size; and those with growth rates that are either negative or positive enough (“sufficiently subcritical” or “sufficiently supercritical,” respectively) to restrict persistence times and subpopulation sizes, conditioned upon the loss of the first mutation in the absence of a second, rescuing mutation. The relative contribution of each regime shifts with initial maladaptation and mutation rate; rare mutations that can occasionally reach unusually large subpopulation sizes play a larger role when the population is not severely maladapted (e.g., Figure 7A) or mutation rate is high (e.g., Figure S2C). In contrast, when populations are initially very maladapted (e.g., Figure 7C), most first-step mutations are themselves also very maladapted, and, thus, restricted in the subpopulation sizes they are expected to reach before being lost. All three regimes help to maintain the variance in the distribution of fitness effects among rescue genotypes as initial maladaptation increases; meanwhile, in one-step rescue, the variance declines due to ever more extreme sampling of the tail of the mutational distribution (compare blue and red curves in Figure 6).

Our prediction, that rescue by more de novo mutations can be more likely than rescue by fewer, is novel. In previous models (e.g., Antia et al. 2003; Iwasa et al. 2004a; Alexander and Day 2010), the general conclusion has been that, since the probability of rescue scales with $U^k$ (where U is the mutation rate and k is the minimum number of mutations required for rescue), the probability of rescue declines with the number of mutations. This assumes, however, that the probability of a mutation occurring, U, is the limiting factor. Here, we have shown that, when the probability of a beneficial mutation arising declines with its selective advantage, the probability of sampling once from the extreme tail of the DFE can be lower than sampling multiple mutations closer to the bulk of the DFE, so that rescue via multiple mutations can become the dominant path. Rescue by multiple mutations may also be more likely with standing genetic variation, as small-effect intermediate mutations may segregate at higher frequencies than large-effect rescue mutations before the environmental change (and also decline less quickly than the wild type following environmental change); this is especially true with recombination, where rescue genotypes can arise from segregating intermediate mutations without mutation (Uecker and Hermisson 2016).

How often rescue arises as a result of multiple mutations is an open question. It is clear that more than one mutation can contribute to adaptation to near-lethal stress, but experiments are often designed to avoid extinction (reviewed in Cowen et al. 2002), and, therefore, greatly expand the scope for multiple mutations to arise on a single genotype. A few exceptions provide some insight. For example, populations of Saccharomyces cerevisiae that survived high concentrations of copper acquired multiple mutations (Gerstein et al. 2015)—in fact, the authors argue for the “springboard effect” discussed above, where first-step mutations prolong persistence, and, thereby, allow further mutations to arise. In Pseudomonas fluorescens, fluctuation tests with nalidixic acid showed that nearly one-third of the most resistant surviving strains were double mutants (Bataillon et al. 2011), which were able to tolerate 10 times higher drug concentrations than single mutants, suggesting two-step rescue might dominate at high drug concentrations. While suggestive, it is unclear if our prediction—that rescue takes more mutational steps with greater initial maladaptation—holds true generally. Verification will require more experiments that allow extinction and uncover the genetic basis of adaptation at different severities of environmental change (e.g., drug concentration).

In describing the genetic basis of adaptation in populations of constant size, Orr (1998) showed that the mean phenotypic displacement toward the optimum scales roughly linearly with initial displacement. Converting phenotype to fitness, this implies that the mean fitness effect of fixed mutations $\left(s=m-m_0\right)$ increases exponentially as initial Malthusian fitness $\left(m_0\right)$ declines (i.e., $s\sim \text{exp}\left(-m_0\right)$⁠), which is a roughly linear increase when initial fitness is small $\left(\left|m_0\right|\ll 1\right)$⁠. Here, we see that, under one-step rescue, the mean fitness effect also increases roughly linearly as the initial growth rate declines (see horizontal blue lines in Figure 6). However, the rate of this linear increase in fitness effect is much larger under rescue than in a population of constant size (compare blue and gray horizontal lines in Figure 6), where declines in wild-type fitness not only allow larger mutations to be beneficial, but also require larger mutations for persistence. Thus the race between extinction and adaptation during evolutionary rescue is expected to produce a genetic basis of adaptation with fewer mutations of larger effect.

While, under one-step rescue, the fitness effect of the first mutation increases roughly linearly as wild-type fitness declines, most rescue events will be two-step for wild-type fitnesses below some value (e.g., at $m_0\approx -0.25$ in Figure 3; this threshold value of $m_0$ increases with mutation rate, Figure S1). At this junction, the effect size of the first mutation will no longer increase as quickly (and potentially even decrease), as it switches from a rescue mutant to an intermediate mutant whose expected fitness begins to decline substantially with the fitness of the wild type (Figure 7). Thus, as rescue switches from dominantly k-step to dominantly $\left(k+1\right)$-step the genetic basis of adaptation becomes more diffuse, with each mutation having a smaller individual fitness effect as the contributing fitness effects spread over more loci. In the limit of large k (due to large initial maladaptation or high mutation rates), the genetic basis of adaptation should, at some point, converge to many loci with small effect, as would also be expected in a population of constant size. Indeed, at very high mutation rates the rate of adaptation (the change in mean fitness) is the same under rescue as it is in populations of constant size (Anciaux et al. 2019), implying that the genetic basis of adaptation no longer depends on demography. It is therefore at intermediate levels of initial maladaptation and low mutation rates, where rescue primarily occurs from a few large effect mutations, that the race between adaptation and persistence is predicted to have the largest effect on the genetic basis of adaptation.

Fluctuation tests (Luria and Delbrück 1943) provide a means to generate random mutations, and then isolate potential rescue genotypes (typically assumed to be one-step only), whose growth rates can be measured under the selective conditions. These experiments are designed such that there is substantial standing genetic variation at the time of exposure to the selective conditions, which should increase the contributions of mutations with small growth rates (Orr and Betancourt 2001), although these could be outcompeted by mutations with higher growth rates and/or be undersampled. Regardless, consistent with our theory (Figure 6), the resulting growth rate distributions in both bacteria and yeast often find modes that are substantially greater than zero (as opposed to, say, an exponential distribution; Kassen and Bataillon 2006; MacLean and Buckling 2009; Gerstein et al. 2012; Lindsey et al. 2013; Gerstein et al. 2015). A number of these conform even more closely to our expected shape (Kassen and Bataillon 2006; Gerstein et al. 2015), while others appear to be substantially more clumped around the mode, perhaps due to a very restricted number of possible rescue mutations in any one circumstance, the size of the experiment, or the way in which growth rates are measured. Finally, Gerstein et al. (2015) not only provide the distribution of growth rates among rescue genotypes, but also the growth rates of individual mutations that compose multi-step rescue genotypes. In four lines where multiple mutations were detected and a segregation analysis performed, one mutation in each line was inferred to have a minor effect, and the other mutation was an amplification of the copper metallothionein gene (CUP1) with a major fitness effect. These results are consistent with the minor effect mutations being subcritical mutations that provided a springboard for the larger CUP1 mutations.

Pinpointing the mutations responsible for adaptation is hampered by genetic hitchhiking, as beneficial alleles elevate the frequency of linked neutral and mildly deleterious alleles (Barton 2000). The problem is particularly severe under strong selection and low recombination, and therefore reaches an extreme in the case of evolutionary rescue in asexuals, especially if many neutral and deleterious mutations are segregating at the time of environmental change. To circumvent this, mutations that have risen to high frequency in multiple replicates are often introduced in a wild-type background, in isolation, and sometimes also in combination with a small number of other common high-frequency mutations, and grown under selective conditions (e.g., Jochumsen et al. 2016; Ono et al. 2017). As we have demonstrated above (e.g., Figure 7C), however, under multi-step rescue, there may be no one mutation that individually confers growth in the selective conditions. Thus, a mutation that was essential for rescue may go undetected, or be mistaken as a hitchhiker, if the appropriate multiple-mutation genotypes are not tested. Unfortunately, reverse engineering all combinations of mutations quickly becomes unwieldy as the number of mutations grows, and thus this approach will not be practical under severe initial maladaptation and high mutation rates, where we predict rescue to occur by many mutations. Interestingly, our simulations show that the population dynamics themselves may help differentiate how many mutations contribute to rescue (e.g., V- vs. U-shaped log-trajectories; Figures 1 and 2), and fitting models of k-step rescue could produce estimates for growth rates of the k genotypes.

Environmental change often selects for mutator alleles, which elevate the rate at which beneficial alleles arise, and, subsequently, increase in frequency with them (Tenaillon et al. 2001). When beneficial alleles are required for persistence, as in evolutionary rescue, mutator alleles can reach very high frequencies or rapidly fix (e.g., Mao et al. 1997). Consistent with this, mutator alleles are often associated with antibiotic resistance in clinical isolates (see examples in Bell 2017). Further, the more beneficial mutations available, the larger the advantage of a mutator allele; for a mutator that increases the mutation rate m-fold, its relative contribution to the production of n beneficial mutations scales as $m^n$ (Tenaillon et al. 1999). Thus, conditions that cause multi-step rescue to be more likely than one-step rescue should also impose stronger selection for mutator alleles. There are a number of examples where lineages with higher mutation rates acquired multiple mutations and persisted at higher doses of antibiotics (Couce et al. 2015; San Millan et al. 2016). The number of mutations required for persistence is, however, often unknown, making it difficult to compare situations where rescue requires different numbers of mutations. Experiments with a combination of drugs may provide a glimpse; for instance, Escherichia coli populations only evolved resistance to a combination of two drugs (presumably through the well-known mutations specific to each drug) when mutators were present, despite the fact that mutators were not required for resistance to either drug in isolation (Gifford et al. 2019). In cases where we have less information on the genetic basis of resistance, our model suggests that mutators will be more advantageous when initial maladaptation is severe (e.g., higher drug concentrations or a larger number of drugs), as rescue will then be dominated by genetic paths with more mutational steps.

Here, we have investigated the genetic basis of evolutionary rescue in an asexual population that is initially genetically uniform. Extending this work to allow for recombination and standing genetic variation at the time of environmental change—as expected for many natural populations—would be valuable. The effect of standing genetic variation on the probability of one-step rescue is relatively straightforward to incorporate, depending only on the expected number of rescue mutations initially present and their mean establishment probability (Martin et al. 2013). In the case of the fluctuation tests discussed above, where mutations accumulated in the short interval before the onset of selection are assumed to be relatively neutral, the effect of standing genetic variance on one-step rescue might be incorporated by a simple rescaling of $N_0$⁠, to account for the additional mutants present in the standing variation. When considering longer periods of time in populations that are not rapidly expanding, mutation-selection balance may be reached before the onset of selection. In this case, the probability of one-step rescue from standing genetic variance in Fisher’s geometric model was given by Anciaux et al. (2018), whose Equation 3 and Equation 5 immediately give the distribution of fitness effects among those that rescue. Allowing these standing genetic variants to be springboards to multi-step rescue will help clarify the role of standing genetic variation on the genetic basis of rescue more generally. Recombination can help combine such springboard mutations into rescue genotypes but will also break these combinations apart, as demonstrated in a two-locus two-allele model of rescue (Uecker and Hermisson 2016). How recombination affects the genetic basis of evolutionary rescue when more loci can potentially contribute remains to be seen. Also left unexplored is the effect of density-dependent fitness; for example, competition may reduce mutant growth rates, and, thereby, increase the size of mutations that are required for rescue, especially when the wild type declines slowly. Combining density-dependence and standing genetic variance is known to create complex dynamics in a one-locus two-allele model of rescue (Uecker et al. 2014), and adding more potential genotypes is sure to add yet more complexity.

Many of our simple closed-form results rely upon knowing the distribution of mutant growth rates (Equation 1), which arises from the assumption that mutant phenotypes are distributed normally about their ancestor, and that Malthusian fitness is a quadratic, on some scaled phenotypic axes. It is clear that deviations from these assumptions will affect our results, at least quantitatively. For instance, mutant phenotype distributions with truncated or fat tails are likely to lead to smaller or larger mutational steps, respectively, with downstream effects on the probability of rescue, the number of contributing mutations, and the resulting DFEs. As a preliminary investigation of this prediction, we performed simulations with mutant phenotype distributions having the same expectation and covariances as assumed above under normality, but with truncated (platykurtic) or fat (leptokurtic) tails (Figure S3A). While our qualitative results above hold, the probability of rescue declines slower with wild-type maladaptation when the mutational distribution has fatter tails (compare dotted and solid black in Figure S3C). Fatter tails also reduce the number of mutations contributing to rescue (e.g., one-step rescue dominates for all wild-type decline rates in Figure S3C). Finally, fatter tails cause the distributions of rescue genotype growth rates following one- and two-step rescue to have more variance and become more similar to one another (Figure S4B), and also tend to increase the contribution of supercritical single mutants in two-step rescue (Figure S5). All told, the genetic basis of rescue is expected to consist of fewer mutations of larger effect, with less consistent effect sizes across replicate populations, as the tails of the mutant phenotype distribution become fatter.

In the numerical examples above, we have not varied the number of scaled phenotypic axes, n, i.e., the dimensionality of the phenotypic landscape (although the analytical results apply for arbitrary n). Because increasing the number of dimensions changes the distribution of fitness effects, and, in particular, decreases the proportion of mutations that are beneficial (Fisher 1930), this may have cascading influences on our results. As shown in Anciaux et al. (2018), the probability of one-step rescue by de novo mutation declines with dimensionality, and is only weakly dependent on dimensionality when initial maladaptation is small [such that ${\text{Λ}}_1\left(m_0\right)\approx -m_{0}Ug\left(\alpha \right)$⁠, Equation 19]. Here, we show that the distribution of fitness effects among one-step rescue mutants is nearly independent of dimensionality for any degree of initial maladaptation (Equation 16 and the blue curves in Figure S6B). Further, as seen by comparing Equations 11–14 to Equation 19, the probability of two-step rescue depends on dimensionality much like one-step rescue does, suggesting that, while increasing dimensionality may decrease the probability of rescue, it may have little effect on the number of steps rescue tends to take. This is demonstrated more generally in Figure S6A, where an order of magnitude increase in the number of dimensions decreases the probability of rescue by roughly an order of magnitude, but has little effect on the relative rates of 1-, 2-, 3-, and 4-step rescue. Finally, Figure S6, B and C shows that dimensionality has very little effect on the distribution of fitness effects among two-step rescue genotypes (Equation 17), and among first-step mutants leading to two-step rescue (Equation 18). To conclude, while the probability of rescue declines with the complexity of the organism and its environment, the genetic basis of rescue is expected to be relatively invariant across complexity, as with the genetic basis of adaptation in populations of constant size (Orr 1998, see also gray curves in Figure S6, B and C).

In the numerical examples above, we have also focused on a particular value of mutational variance, λ. Clearly, since rescue relies on mutations of large effect, decreasing λ should decrease the probability of rescue, much like decreasing the mutation rate, U, does (Figure S1). While our analysis (Equation 11–14 and Equation 19) and numerical results (see File S2) show that this is true, we find that λ and U have very different effects on the genetic basis of rescue (File S2). In particular, given a similar effect on the total probability of rescue, decreasing U generally restricts rescue to fewer mutational steps while decreasing λ forces rescue to occur by more mutations. Further, the distribution of fitness effects of mutations contributing to rescue is nearly independent of U but a decrease in λ strongly reduces the mode of the DFE. This demonstrates that populations with similar probabilities of rescue can vary greatly in the way they achieve it genetically.

Appendix

Approximating the probability of one-step rescue

The probability of one-step rescue in this model has been derived by Anciaux et al. (2018). As replicated in File S2 and given by their Equation 7, when ${\rho }_{\mathit{max}}=m_{\mathit{max}}/\lambda$ is large a simple, nearly closed-form approximation is

$${\text{Λ}}_1\left(m_0\right)\approx {\stackrel{\sim }{\text{Λ}}}_1\left(m_0\right)\equiv -m_{0}U\frac{{\left(1-{\psi }_0/2\right)}^{\left(1-n\right)/2}}{1-{\psi }_0/4}g\left(\alpha \right),$$

(19)

where ${\psi }_0=2\left(1-\sqrt{1-m_0/m_{\mathit{max}}}\right)$⁠, $g\left(\alpha \right)=\text{exp}\left(-\alpha \right)/\sqrt{\pi \alpha }-\text{erfc}\left(\sqrt{\alpha }\right)$⁠, and $\alpha ={\rho }_{\mathit{max}}{\psi }_0^2/4$⁠, with $\text{erfc}(.)$ the complimentary error function. When the wild-type declines slowly $m_0$ and thus ${\psi }_0$ is small and ${\text{Λ}}_1\left(m_0\right)\approx Ug\left(\alpha \right)$⁠. In the limit $m_0\to 0$⁠, Equation 19 becomes

$${\stackrel{\sim }{\text{Λ}}}_1\left(0\right)\equiv \underset{m_0\to 0}{\text{lim}}{\stackrel{\sim }{\text{Λ}}}_1\left(m_0\right)=2U\sqrt{m_{\mathit{max}}\lambda /\pi }.$$

(20)

Mutant lineage dynamics

Here, we follow the lead of Weissman et al. (2010) and Uecker and Hermisson (2016) in approximating our discrete-time process with a continuous-time branching process (see Chapter 6 in Allen, 2010). Consider a birth–death process, where individuals give birth at rate b and die at rate d. One can then obtain the probability generating function for the number of individuals at a given time, $n\left(t\right)$⁠, given the initial number, $n\left(0\right)$⁠. We are primarily interested in new mutant lineages, $n\left(0\right)=1$⁠. The generating function then allows us to calculate the probability that a lineage persists at least until time t, and the distribution of $n\left(t\right)$ given it does so (see below).

To convert between birth and death rates and our compound Malthusian parameter we follow Uecker and Hermisson (2016) in equally distributing the growth rate m between birth and death, $b=\left(1+m\right)/2$ and $d=\left(1-m\right)/2$⁠, such that $m=b-d$ and the continuous-time process exhibits the same amount of drift as the discrete time process (and matches discrete-time simulations well; Uecker et al., 2014). We can now report the necessary results in terms of m (assuming $\left|m\right|<1$⁠).

Denoting the extinction time as T, the probability a mutant with growth rate m persists until time t is approximately (see File S2 for derivation)

$$P\left(T>t\right)\approx \{\begin{array}{cc}2/t& t\ll \left|1/m\right|\\ -2m\text{\hspace{0.05em}\hspace{0.17em}exp}\left(mt\right)& t\gg -1/m>0\end{array}$$

(21)

As pointed out in Weissman et al. (2010) (whose equation A2 differs from Equation 21 by a factor of two because they have $b+d=2$⁠), the distribution of persistence times has a long tail (like $1/t$⁠) until being cut off (declining exponentially) at $t=-1/m$⁠.

Given that a lineage persists until t, the distribution of $n\left(t\right)$ is roughly (see File S2 for derivation)

$$P\left(n\left(t\right)=n|n\left(t\right)>0\right)\approx \{\begin{array}{cc}2\left(1/t\right){\left(1+2/t\right)}^{-n}& t\ll \left|1/m\right|\\ -2m{\left(1+m\right)}^{n-1}& t\gg -1/m>0\end{array}$$

(22)

As pointed out in Weissman et al. (2010) (whose equation A3 only differs from Equation 22 by constants), the distribution of $n\left(t\right)$ is approximately geometric for small or large t, implying $n\left(t\right)$ is very unlikely to be greater than the minimum of t and $-1/m$⁠.

Acknowledgments

We would like to thank the Otto and Doebeli laboratories for helpful feedback at various stages, Ophélie Ronce and Thomas Lenormand for their hospitality and valuable input at the beginning of this project, and Mike Whitlock, Amy Angert, Luis-Miguel Chevin, and Joachim Hermisson for constructive criticism on previous versions of the manuscript. Funding provided by the National Science and Engineering Research Council (CGS-D 6564 to M.M.O., RGPIN-2016-03711 to S.P.O.), the University of British Columbia, Banting, and the University of California Davis (fellowships to M.M.O.), the National Institute of General Medical Sciences of the National Institutes of Health (NIH R01 G.M.108779 to Graham Coop), the Agence Nationale de la Recherche (ANR-18-CE45-0019 “RESISTE” to G.M.), and the Centre Méditerranéen Environment et Biodiversité (“BACTPHI” to G.M.).

Footnotes

Supplemental material available at figshare: https://doi.org/10.25386/genetics.11343704.

Communicating editor: J. Hermisson

Author notes

Present address: Center for Population Biology, University of California, Davis, CA 95616.