The population genetics of collateral resistance and sensitivity
Abstract
Resistance mutations against one drug can elicit collateral sensitivity against other drugs. Multidrug treatments exploiting such tradeoffs can help slow down the evolution of resistance. However, if mutations with diverse collateral effects are available, a treated population may evolve either collateral sensitivity or collateral resistance. How to design treatments robust to such uncertainty is unclear. We show that many resistance mutations in Escherichia coli against various antibiotics indeed have diverse collateral effects. We propose to characterize such diversity with a joint distribution of fitness effects (JDFE) and develop a theory for describing and predicting collateral evolution based on simple statistics of the JDFE. We show how to robustly rank drug pairs to minimize the risk of collateral resistance and how to estimate JDFEs. In addition to practical applications, these results have implications for our understanding of evolution in variable environments.
Editor's evaluation
When selecting for one particular trait, it is not uncommon for other traits to change. This is due to pleiotropic mutations that affect multiple characters. Ardell and Kryazhimskiy develop a theoretical framework to predict adaptive trajectories observed in environments other than the one selection is operating in. The effects of adaptation across environments have important implication to antibiotic treatments, where resistance evolution to one antibiotic can alter the susceptibility to other antibiotics.
https://doi.org/10.7554/eLife.73250.sa0eLife digest
Drugs known as antibiotics are the main treatment for most serious infections caused by bacteria. However, many bacteria are acquiring genetic mutations that make them resistant to the effects of one or more types of antibiotics, making them harder to eliminate.
One way to tackle drugresistant bacteria is to develop new types of antibiotics; however, in recent years, the rate at which new antibiotics have become available has been dwindling. Using two or more existing drugs, one after another, can also be an effective way to eliminate resistant bacteria. The success of any such ‘multidrug’ treatment lies in being able to predict whether mutations that make the bacteria resistant to one drug simultaneously make it sensitive to another, a phenomenon known as collateral sensitivity.
Different resistance mutations may have different collateral effects: some may increase the bacteria’s sensitivity to the second drug, while others might make the bacteria more resistant. However, it is currently unclear how to design robust multidrug treatments that take this diversity of collateral effects into account.
Here, Ardell and Kryazhimskiy used a concept called JDFE (short for the joint distribution of fitness effects) to describe the diversity of collateral effects in a population of bacteria exposed to a single drug. This information was then used to mathematically model how collateral effects evolved in the population over time.
Ardell and Kryazhimskiy showed that this approach can predict how likely a population is to become collaterally sensitive or collaterally resistant to a second antibiotic. Drug pairs can then be ranked according to the risk of collateral resistance emerging, so long as information on the variety of resistance mutations available to the bacteria are included in the model.
Each year, more than 700,000 people die from infections caused by bacteria that are resistant to one or more antibiotics. The findings of Ardell and Kryazhimskiy may eventually help clinicians design multidrug treatments that effectively eliminate bacterial infections and help to prevent more bacteria from evolving resistance to antibiotics. However, to achieve this goal, more research is needed to fully understand the range collateral effects caused by resistance mutations.
Introduction
The spread of resistance against most antibiotics and the difficulties in developing new ones has sparked considerable interest in using drug combinations and sequential drug treatments to treat bacterial infections, as well as cancers (Pál et al., 2015). Treatments where the drugs are chosen so that resistance against one of them causes the pathogen or cancer population to become sensitive to the other—a phenomenon known as collateral sensitivity—can eliminate the population before multidrug resistance emerges (Pál et al., 2015; Pluchino et al., 2012).
The success of a multidrug treatment hinges on knowing which drugs select for collateral sensitivity against which other drugs. This information is obtained empirically by exposing bacterial and cancercell populations to drugs and observing the evolutionary outcomes (Roemhild et al., 2020; Jensen et al., 1997; Imamovic and Sommer, 2013; Lázár et al., 2018; Maltas and Wood, 2019; Batra et al., 2021; SanzGarcía et al., 2020; Schenk et al., 2015; Lázár et al., 2013; Barbosa et al., 2019; HernandoAmado et al., 2020; Kim et al., 2014; Jahn et al., 2021; Kavanaugh et al., 2020; Laborda et al., 2021; Oz et al., 2014; Munck et al., 2014). Prior studies have largely focused on various empirical questions related to the evolution of collateral sensitivity and resistance, such as identifying their genetic basis (Lázár et al., 2014; Roemhild et al., 2020; Maltas and Wood, 2019; HernandoAmado et al., 2020; Kavanaugh et al., 2020; Laborda et al., 2021), understanding how collateral outcomes depend on treatment design (e.g. sequential versus combination) (Lázár et al., 2014; Munck et al., 2014; Bergstrom et al., 2004; Batra et al., 2021; SanzGarcía et al., 2020; Schenk et al., 2015; Lázár et al., 2013; Kim et al., 2014; Jahn et al., 2021), or testing whether collateral sensitivity is an evolutionarily stable outcome (Barbosa et al., 2019). However, one important feature of these experimental studies has received little attention, namely, the fact that different experiments often produce collateral sensitivity profiles that are inconsistent with each other (e.g. Imamovic and Sommer, 2013; Oz et al., 2014; Barbosa et al., 2017; Maltas and Wood, 2019). Some inconsistencies can be explained by the fact that resistance mutations vary between bacterial strains, drug dosages, etc. (Mira et al., 2015; Barbosa et al., 2017; Das et al., 2020; Pinheiro et al., 2021; Card et al., 2021; Gjini and Wood, 2021). However, wide variation in collateral outcomes is observed even between replicate populations (Oz et al., 2014; Barbosa et al., 2017; Maltas and Wood, 2019; Nichol et al., 2019). This variation suggests that bacteria and cancers have access to multiple resistance mutations with diverse collateral effects and that replicate populations accumulate different resistance mutations due to the intrinsic randomness of the evolutionary process (Jerison et al., 2020; Nichol et al., 2019). However, the diversity of collateral effects among resistance mutations has rarely if ever been systematically measured. Moreover, few existing approaches for designing robust multidrug treatments have modelled this mutational diversity explicitly within the population genetics context (Nichol et al., 2019; Maltas and Wood, 2019). Yet, a theory grounded in population genetics could help us understand how the expected collateral outcomes and the uncertainty around these expectations depend on evolutionary parameters and how these expectations and uncertainties change over time.
Here, we develop such a theory. Collateral sensitivity and resistance are specific examples of the more general evolutionary phenomenon, pleiotropy, which refers to any situation when one mutation affects multiple phenotypes (Wagner and Zhang, 2011; Paaby and Rockman, 2013). In case of drug resistance evolution, the direct effect of resistance mutations is to increase fitness in the presence of one drug (the ‘home’ environment). In addition, they may also provide pleiotropic gains or losses in fitness in the presence of other drugs (the ‘nonhome’ environments) leading to collateral resistance or sensitivity, respectively.
Classical theoretical work on pleiotropy has been done in the field of quantitative genetics (Lande and Arnold, 1983; Rose, 1982; Barton, 1990; Slatkin and Frank, 1990; Jones et al., 2003; Johnson and Barton, 2005). In these models, primarily developed to understand how polygenic traits respond to selection in sexual populations, pleiotropy manifests itself as a correlated temporal change in multiple traits in a given environment. The question of how new strongly beneficial mutations that accumulate in asexual populations evolving in one environment affect its fitness in future environments is outside of the scope of these models.
The pleiotropic consequences of adaptation have also been explored in various ‘fitness landscape’ models (e.g. Connallon and Clark, 2015; Martin and Lenormand, 2015; Harmand et al., 2017; Wang and Dai, 2019; Maltas et al., 2021; Nichol et al., 2019; Tikhonov et al., 2020). In particular, Nichol et al., 2019 specifically addressed the problem of diversity of collateral resistance/sensitivity outcomes in the context of a combinatorially complete fitness landscapes of four mutations in the TEM βlactamase gene. They found that different in silico populations adapting to the same antibiotic often arrive at different fitness peaks which results in different levels of collateral resistance/sensitivity against other drugs. They observed qualitatively similar variability in the collateral outcomes among replicate populations of the bacterium Escherichia coli evolving in the presence of cefotaxime (CTX), although it is unclear whether different populations indeed arrived at different fitness peaks. In general, the fitness landscape approach helps us understand how evolutionary trajectories and outcomes depend on the global structure of the underlying fitness landscape. However, applying this approach in practice is challenging because the global structure of fitness landscapes is unknown and notoriously difficult to estimate, even in controlled laboratory conditions.
Here, we take a different approach which is agnostic with respect to the global structure of the fitness landscape. Instead, we assume only the knowledge of the socalled joint distribution of fitness effects (JDFE), that is, the probability that a new mutation has a certain pair of fitness effects in the home and nonhome environments (Jerison et al., 2014; Martin and Lenormand, 2015; Bono et al., 2017). The JDFE is a natural extension of the DFE, the distribution of fitness effects of new mutations, often used in modeling evolution in a single environment (King, 1972; Ohta, 1987; Orr, 2003; Kassen and Bataillon, 2006; EyreWalker and Keightley, 2007; Martin and Lenormand, 2008; MacLean and Buckling, 2009; Kryazhimskiy et al., 2009; Levy et al., 2015). Like the DFE, the JDFE is a local property of the fitness landscape which means that it can be, at least in principle, estimated by using a variety of modern highthroughput techniques (e.g. Qian et al., 2012; van Opijnen et al., 2009; Chevereau et al., 2015; Levy et al., 2015; Blundell et al., 2019; Aggeli et al., 2021). The downside of this approach is that the JDFE can change over time as the population traverses the fitness landscape (Good et al., 2017; Venkataram et al., 2020; Aggeli et al., 2021). However, in the context of collateral drug resistance and sensitivity, we are primarily interested in short time scales over which the JDFE can be reasonably expected to stay approximately constant.
The rest of the paper is structured as follows. First, we use previously published data to demonstrate that E. coli has access to drug resistance mutations with diverse collateral effects. This implies that, rather than treating collateral effects as deterministic properties of drug pairs, we should think of them probabilistically, in terms of the respective JDFEs. We then show that a naive intuition about how the JDFE determines pleiotropic outcomes of evolution can sometimes fail, and a mathematical model is therefore required. We develop such a model, which reveals two key ‘pleiotropy statistics’ of the JDFE that determine the dynamics of fitness in the nonhome condition. Our theory makes quantitative predictions in a variety of regimes if the population genetic parameters are known. However, we argue that in the case of drug resistance evolution the more important problem is to robustly order drug pairs in terms of their collateral sensitivity profiles even if the population genetic parameters are unknown. We develop a metric that allows us to do so. Finally, we provide some practical guidance for estimating the pleiotropy statistics of empirical JDFEs in the context of ranking drug pairs.
Results
Antibiotic resistance mutations in E. coli have diverse collateral effects
We begin by demonstrating that the JDFE is a useful concept for modeling the evolution of collateral antibiotic resistance and sensitivity. If all resistance mutations against a given drug had identical pleiotropic effects on the fitness of the organism in the presence of another drug, the dynamics of collateral resistance/sensitivity could be understood without the JDFE concept. On the other hand, if different resistance mutations have different pleiotropic fitness effects, predicting the collateral resistance/sensitivity dynamics requires specifying the probabilities with which mutations with various home and nonhome fitness effects arise in the population. The JDFE specifies these probabilities. Therefore, for the JDFE concept to be useful in the context of collateral resistance/sensitivity evolution, we need to show that resistance mutations against common drugs have diverse collateral effects in the presence of other drugs.
To our knowledge, no data sets are currently publicly available that would allow us to systematically explore the diversity of collateral effects among all resistance mutations against any one drug in any organism. Instead, we examined the fitness effects of 3883 gene knockout mutations in the bacterium Escherichia coli, measured in the presence of six antibiotics (Chevereau et al., 2015), as well as the fitness effects of 4997 point mutations in the TEM1 βlactamase gene measured in the presence of two antibiotics (Stiffler et al., 2015).
For four out of six antibiotics used by Chevereau et al., 2015, we find between 12 (0.31 %) and 170 (4.38%) knockout mutations that provide some level of resistance against at least one of the antibiotics (false discovery rate (FDR) ∼25%; Figure 1, Figure 1—source data 1; see Materials and methods for details). Plotting on the $x$axis the fitness effect of each knockout mutation in the presence of the drug assumed to be applied first (i.e. the home environment) against its effect in the presence of another drug assumed to be applied later (i.e. the nonhome environment, $y$axis), we find mutations in all four quadrants of this plane, for all 12 ordered drug pairs (Figure 1, Figure 1—source data 1). Similarly, we find diverse collateral effects among mutations within a single gene (Figure 1—figure supplement 1; see Materials and methods for details).
Since both data sets represent subsets of all resistance mutations, we conclude that E. coli likely have access to resistance mutations with diverse pleiotropic effects, such that a fitness gain in the presence of any one drug can come either with a pleiotropic gain or a pleiotropic loss of fitness in the presence of other drugs. Therefore, the JDFE framework is suitable for modeling the evolution of collateral resistance/sensitivity. In the next section, we formally define a JDFE and probe our intuition for how its shape determines the fitness trajectories in the nonhome environment.
JDFE determines the pleiotropic outcomes of adaptation
For any genotype $g$ that finds itself in one (‘home’) environment and may in the future encounter another ‘nonhome’ environment, we define the JDFE as the probability density ${\mathrm{\Phi}}_{g}(\mathrm{\Delta}x,\mathrm{\Delta}y)$ that a new mutation that arises in this genotype has the selection coefficient $\mathrm{\Delta}x$ in the home environment and the selection coefficient $\mathrm{\Delta}y$ in the nonhome environment (Jerison et al., 2014). For concreteness, we define the fitness of a genotype as its Malthusian parameter (Crow and Kimura, 1972). So, if the home and nonhome fitness of genotype $g$ are $x$ and $y$, respectively, and if this genotype acquires a mutation with selection coefficients $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$, its fitness becomes $x+\mathrm{\Delta}x$ and $y+\mathrm{\Delta}y$. This definition of the JDFE can, of course, be naturally extended to multiple nonhome environments. In principle, the JDFE can vary from one genotype to another. However, to develop a basic intuition for how the JDFE determines pleiotropic outcomes, we assume that all genotypes have the same JDFE. We discuss a possible extension to epistatic JDFEs in Appendix 1.
The JDFE is a complex object. So, we first ask whether some simple and intuitive summary statistics of the JDFE may be sufficient to predict the dynamics of the nonhome fitness of a population that is adapting in the home environment. Intuitively, if there is a tradeoff between home and nonhome fitness, nonhome fitness should decline; if the opposite is true, nonhome fitness should increase. Canonically, a tradeoff occurs when any mutation that improves fitness in one environment decreases it in the other environment and vice versa (Roff and Fairbairn, 2007). Genotypes that experience such ‘hard’ tradeoffs are at the Pareto front (Shoval et al., 2012; Li et al., 2019). For genotypes that are not at the Pareto front, some mutations that are beneficial in the home environment may be beneficial in the nonhome environment and others may be deleterious. In this more general case, tradeoffs are commonly quantified by the degree of negative correlation between the effects of mutations on fitness in the two environments (Roff and Fairbairn, 2007; Tikhonov et al., 2020). Thus, we might expect that evolution on negatively correlated JDFEs would lead to pleiotropic fitness losses and evolution on positively correlated JDFEs would lead to pleiotropic fitness gains.
To test this intuition, we generated a family of Gaussian JDFEs that varied, among other things, by their correlation structure (Figure 2; Materials and methods). We then simulated the evolution of an asexual population on these JDFEs using a standard WrightFisher model (Materials and methods) and tested whether the tradeoff strength, measured by the JDFE’s correlation coefficient, predicts the dynamics of nonhome fitness. Figure 2 shows that our naive expectation is incorrect. Positively correlated JDFEs sometimes lead to pleiotropic fitness losses (Figure 2D,I), and negatively correlated JDFEs sometimes lead to pleiotropic fitness gains (Figure 2B,G). Even if we calculate the correlation coefficient only among mutations that are beneficial in the home environment, the pleiotropic outcomes still do not always conform to the naive expectation, as the sign of the correlation remains the same as for the full JDFEs in all these examples.
There are other properties of the JDFE that we might intuitively expect to be predictive of the pleiotropic outcomes of adaptation. For example, among the JDFEs considered in Figure 2, it is apparent that those with similar relative probability weights in the first and fourth quadrants produce similar pleiotropic outcomes. However, simulations with other JDFE shapes show that even distributions that are similar according to this metric can also result in qualitatively different pleiotropic outcomes (Figure 2—figure supplement 1).
Overall, our simulations show that JDFEs with apparently similar shapes can produce qualitatively different trajectories of pleiotropic fitness changes (e.g. compare Figure 2A, F and B, G or Figure 2D,I and E,J). Conversely, JDFEs with apparently different shapes can result in rather similar pleiotropic outcomes (e.g. compare Figure 2B, G and E, J or Figure 2A,F and D,I). Thus, while the overall shape of the JDFE clearly determines the trajectory of pleiotropic fitness changes, it is not immediately obvious what features of its shape play the most important role, particularly if the JDFE is more complex than a multivariate Gaussian. In other words, even if we have perfect knowledge of the fitness effects of all mutations in multiple environments, converting this knowledge into a qualitative prediction of the expected direction of pleiotropic fitness change (gain or loss) does not appear straightforward. Therefore, we next turn to developing a population genetics model that would allow us to predict not only the direction of pleiotropic fitness change but also the expected rate of this change and the uncertainty around the expectation.
The population genetics of pleiotropy
To systematically investigate which properties of the JDFE determine the pleiotropic fitness changes in the nonhome environment, we consider a population of size $N$ that evolves on a JDFE in the ‘strong selection weak mutation’ (SSWM) regime, also known as the ‘successional mutation’ regime (Orr, 2000; Desai and Fisher, 2007; Kryazhimskiy et al., 2009; Good and Desai, 2015).
We consider an arbitrary JDFE without epistasis, that is a situation when all genotypes have the same JDFE $\mathrm{\Phi}(\mathrm{\Delta}x,\mathrm{\Delta}y)$. We explore an extension to JDFEs with a simple form of epistasis in Appendix 1. We assume that mutations arise at rate $U$ per individual per generation. In the SSWM limit, a mutation that arises in the population either instantaneously fixes or instantaneously dies out. Therefore, the population is essentially monomorphic at all times, such that at any time $t$ we can characterize it by its current pair of fitness values $({X}_{t},{Y}_{t})$. If a new mutation with a pair of selection coefficients $(\mathrm{\Delta}x,\mathrm{\Delta}y)$ arises in the population at time $t$, it fixes with probability $\pi \left(\mathrm{\Delta}x\right)=\frac{1{e}^{2\mathrm{\Delta}x}}{1{e}^{2N\mathrm{\Delta}x}}$(Kimura, 1962) in which case the population’s fitness transitions to a new pair of values $({X}_{t}+\mathrm{\Delta}x,{Y}_{t}+\mathrm{\Delta}y)$. If the mutation dies out, an event that occurs with probability $1\pi \left(\mathrm{\Delta}x\right)$, the population’s fitness does not change. This model specifies a continuoustime twodimensional Markov process.
In general, the dynamics of the probability density $p(x,y,t)$ of observing the random vector $({X}_{t},{Y}_{t})$ at values $(x,y)$ are governed by an integrodifferential forward Kolmogorov equation, which is difficult to solve (Materials and methods). However, if most mutations that contribute to adaptation have small effects, these dynamics are well approximated by a diffusion equation which can be solved exactly (Materials and methods). Then $p(x,y,t)$ is a normal distribution with mean vector
and variancecovariance matrix
where are $r}_{1$ and $r}_{2$, given by Equation 7 and Equation 8 in Materials and methods, are the expected fitness effects in the home and nonhome environments for a mutation fixed in the home environment, and $D}_{11},{D}_{12$, and $D}_{22$, given by Equation 9–Equation 11 in Materials and methods, are the second moments of this distribution. Here, ${U}_{b}=U{\int}_{\mathrm{\infty}}^{\mathrm{\infty}}d\eta \phantom{\rule{thinmathspace}{0ex}}{\int}_{0}^{\mathrm{\infty}}d\xi \mathrm{\Phi}(\xi ,\eta )$ is the total rate of mutations beneficial in the home environment, and $x}_{0$ and $y}_{0$ are the initial values of population’s fitness in the home and nonhome environments.
Equation 1 and Equation 2 show that the distribution of population’s fitness at time $t$ in the nonhome environment is entirely determined by two parameters, $r}_{2$ and $D}_{22$, which we call the ‘pleiotropy statistics’ of the JDFE. The expected rate of fitness change in the nonhome environment depends on the pleiotropy statistic $r}_{2$, which we refer to as the expected pleiotropic effect. Thus, evolution on a JDFE with a positive $r}_{2$ is expected to result in pleiotropic fitness gains and evolution on a JDFE with a negative $r}_{2$ is expected to result in pleiotropic fitness losses. Equation 2 shows that the variance around this expectation is determined by the pleiotropy variance statistic $D}_{22$. Since both the expectation and the variance change linearly with time (provided ${r}_{2}\ne 0$), the change in the nonhome fitness in any replicate population will eventually have the same sign as $r}_{2$, but the time scale of such convergence depends on the ‘collateral resistance risk’ statistic $c={r}_{2}/\sqrt{{D}_{22}}$ (Materials and methods). This observation has important practical implications, and we return to it in the Section ‘Robust ranking of drug pairs’.
These theoretical results suggest a simple explanation for the somewhat counterintuitive observations in Figure 2. We may intuitively believe that evolution on negatively correlated JDFEs should lead to fitness losses in the nonhome environment because on such JDFEs mutations with largest fitness benefits in the home environment typically have negative pleiotropic effects. However, such mutations may be too rare to drive adaptation. At the same time, the more common mutations that do typically drive adaptation may have positive pleiotropic effects, in which case the population would on average gain nonhome fitness, as in Figure 2B. Our theory shows that to predict the direction of nonhome fitness change, the frequency of beneficial mutations with various pleiotropic effects and the strength of these effects need to be weighted by the likelihood that these mutations fix. The expected pleiotropic effect $r}_{2$ accomplishes this weighting.
We tested the validity of Equation 1 and Equation 2 by simulating evolution in the SSWM regime on 125 Gaussian JDFEs with various parameters (Materials and methods) and found excellent agreement (Figure 3A and B). However, many microbes likely evolve in the ‘concurrent mutation’ regime, that is, when multiple beneficial mutations segregate in the population simultaneously (Desai and Fisher, 2007; Lang et al., 2013). As expected, our theory fails to quantitatively predict the pleiotropic fitness trajectories when $N{U}_{b}\S gt;1$ (Figure 3—figure supplement 1). However, the expected rate of change of nonhome fitness and its variances remain surprisingly well correlated with the pleiotropy statistics $r}_{2$ and $D}_{22$ across various JDFEs (Figure 3—figure supplement 1). In other words, we can still use these statistics to correctly predict whether a population would lose or gain fitness in the nonhome environment and to order the nonhome environments according to their expected pleiotropic fitness changes and variances. We will exploit the utility of such ranking in the next section.
We next sought to expand our theory to the concurrent mutation regime. A key characteristic of adaptation in this regime is that mutations whose fitness benefits in the home environment are below a certain ‘effective neutrality’ threshold are usually outcompeted by superior mutations and therefore fix with lower probabilities than predicted by Kimura’s formula (Schiffels et al., 2011; Good et al., 2012). Good et al., 2012 provide an equation for calculating the fixation probability ${\pi}^{*}(\mathrm{\Delta}x)$ for a mutation with home fitness benefit $\mathrm{\Delta}x$ in the concurrent mutation regime (Equation (6) in Good et al., 2012). Thus, by replacing $2\xi $ (the approximate fixation probability in the SSWM regime) in Equation 8 and Equation 11 with ${\pi}^{*}(\xi )$, we obtain the adjusted pleiotropy statistics ${r}_{2}^{*}$ and ${D}_{22}^{*}$ for the concurrent mutation regime (see Materials and methods for details). Note that in contrast to $r}_{2$ and $D}_{22$, the adjusted statistics ${r}_{2}^{*}$ and ${D}_{22}^{*}$ depend on the population genetic parameters $N$ and ${U}_{b}$.
To test how well these statistics predict the dynamics of fitness in the nonhome environment, we simulated evolution on the same 125 JDFEs using the full WrightFisher model with a range of population genetic parameters that span the transition from the successional to the concurrent mutation regimes for 1,000 generations. We find that ${r}_{2}^{*}$ quantitatively predicts the expected rate of nonhome fitness change, with a similar accuracy as Good et al., 2012 predict the rate of fitness change in the home environment, as long as $N{U}_{b}\S gt;1$ (Figure 3C, E and G; compare with Figure 3—figure supplement 1A,C,E). ${D}_{22}^{*}$ also predicts the empirically observed variance in nonhome fitness trajectories much better than D_{22}, although this relationship is more noisy than between mean fitness and ${r}_{2}^{*}$ (Figure 3D, F and H; compare with Figure 3—figure supplement 1B,D,F). Some of this noise can be attributed to sampling, as we estimate both the mean and the variance from 300 replicate simulation runs, and the variance estimation is more noisy. Even in the absence of sampling noise however, we do not expect that ${D}_{22}^{*}$ would predict the nonhome fitness variance perfectly because our theory does not account for the autocorrelation in the fitness trajectories that arise in the concurrent mutation regime but not in the successive mutation regime (see Appendix D in Desai and Fisher, 2007). To our knowledge, a rigorous analytical calculation for ensemble variance in fitness even in the home environment is not yet available.
Overall, our theory allows us to quantitatively predict the dynamics of nonhome fitness in a range of evolutionary regimes if the JDFE and the population genetic parameters $N$ and ${U}_{b}$ are known. However, neither the full JDFE nor the population genetic parameters will likely be known in most practical situations, such as designing a drug treatment for a cancer patient. In the next section, we address the question of how to robustly select drug pairs for a sequential treatment, assuming that the pleiotropy statistics $r}_{2$ and $D}_{22$ are known but the population genetic parameters are not. In the Section ‘Measuring JDFEs’, we provide some guidance on how the JDFE can be measured.
Robust ranking of drug pairs
Consider a hypothetical scenario where a drug treatment is being designed for a patient with a tumor or a bacterial infection. In selecting a drug, it is desirable to take into account not only the standard medical considerations, such as drug availability, toxicity, etc., but also the possibility that the treatment with this drug will fail due to the evolution of resistance. Therefore, it may be prudent to consider a list of drugs pairs (or higherorder combinations), ranked by the propensity of the first drug in the pair to elicit collateral resistance against the second drug in the pair. All else being equal, the drug deployed first should form a highranking pair with at least one other secondary drug. Then, if the treatment with the first drug fails, a second one can be deployed with a minimal risk of collateral resistance. Thus, we set out to develop a metric for ranking drug pairs according to this risk.
Clearly, any drug pair with a negative $r}_{2$ is preferable over any drug pair with a positive $r}_{2$, since the evolution in the presence of the first drug in a pair with ${r}_{2}\S lt;0$ is expected to elicit collateral sensitivity against the second drug in the pair but the opposite is true for drug pairs with ${r}_{2}\S gt;0$. It is also clear that among two drug pairs with negative $r}_{2$, a pair with a more negative $r}_{2$ and lower $D}_{22$ is preferable over a pair with a less negative $r}_{2$ and higher $D}_{22$ because evolution in the presence of the first drug in the former pair will more reliably lead to stronger collateral sensitivity against the second drug in the pair. The difficulty is in how to compare and rank two drug pairs where one pair has a more negative $r}_{2$ but higher $D}_{22$. Our theory shows that the chance of emergence of collateral resistance monotonically increases with the collateral risk statistic $c={r}_{2}/\sqrt{{D}_{22}}$ (see Materials and methods). Thus, we propose to rank drug pairs by $c$ from lowest (most negative and therefore most preferred) to highest (least negative or most positive and therefore least preferred).
To demonstrate the utility of such ranking, consider four hypothetical drug pairs with JDFEs shown in Figure 4A. The similarity between their shapes makes it difficult to predict a priori which one would have the lowest and highest probabilities of collateral resistance. Thus, we rank these JDFEs by their $c$ statistic. To test whether this ranking is accurate with respect to the risk of collateral resistance, we simulate the evolution of a WrightFisher population in the presence of the first drug in each pair for 600 generations and estimate the probability that the evolved population has a positive fitness in the presence of the second drug, that is, the probability that it becomes collaterally resistant (Figure 4B). We find that our a priori ranking corresponds perfectly to the ranking according to this probability, evidenced by the consistently higher collateral resistance risk for JDFEs with higher $c$ (Figure 4B). Interestingly, the top ranked JDFE does not have the lowest expected pleiotropic effect $r}_{2$. Nevertheless, the fact that the pleiotropic variance statistic $D}_{22$ for this JDFE is small ensures that the risk of collateral resistance evolution is the lowest. This 1–1 rank correlation holds more broadly, for all 125 Gaussian JDFEs and all population genetic parameters considered in the previous section (Figure 4C) as well as for the empirical TEM βlactamase JDFEs (Figure 4—figure supplement 1). Overall, we find that we can use the collateral risk statistic $c$ to robustly rank drug pairs according to the risk of collateral resistance evolution, irrespective of the population genetic parameters.
Measuring JDFEs
So far, we assumed that the parameters of the JDFE on which the population evolves are known. In reality, they have to be estimated from data, which opens up at least two practically important questions. The first question is experimental. From what types of data can JDFEs be in principle estimated and how good are different types of data for this purpose? We can imagine, for example, that some properties of JDFEs can be estimated from genome sequencing data (Jerison et al., 2020) or from temporally resolved fitness trajectories (Bakerlee et al., 2021). Here, we focus on the most direct way of estimating JDFE parameters, from the measurements of the home and nonhome fitness effects of individual mutations. The experimental challenge with this approach is to sample those mutations that will most likely contribute to adaptation in the home environment (see ‘Discussion’ for an extended discussion of this problem). Below, we propose two potential strategies for such sampling: the LuriaDelbrück (LD) method and the barcode lineage tracking (BLT) method. The second question is statistical: how many mutants need to be sampled to reliably rank drug pairs according to the risk of collateral resistance? We evaluate both proposed methods with respect to this property.
The idea behind the LD method is to expose the population to a given drug at a concentration above the minimum inhibitory concentration (MIC), so that only resistant mutants survive (Pinheiro et al., 2021). This selection is usually done on agar plates, so that individual resistant mutants form colonies and can be isolated. The LD method is relatively easy to implement experimentally, but it is expected to work only if the drug concentration is high enough to kill almost all nonresistant cells. In reality, resistant mutants may be selected at concentrations much lower than MIC (Andersson and Hughes, 2014). Furthermore, mutants selected at different drug concentrations may be genetically and functionally distinct (Lindsey et al., 2013; Pinheiro et al., 2021) and therefore may have statistically different pleiotropic profiles. As a result, mutants sampled with the LD method may not be most relevant for predicting collateral evolution at low drug concentrations, and other sampling methods may be required for isolating weakly beneficial mutations.
Isolating individual weakly beneficial mutations is more difficult because by the time a mutant reaches a detectable frequency in the population it has accumulated multiple additional driver and passenger mutations (Lang et al., 2013; Nguyen Ba et al., 2019), all of which can potentially have collateral effects. One way to isolate many single beneficial mutations from experimental populations is by using the recently developed barcode lineage tracking (BLT) method (Levy et al., 2015; Venkataram et al., 2016). In a BLT experiment, each cell is initially tagged with a unique DNA barcode. As long as there is no recombination or other DNA exchange, any new mutation is permanently linked to one barcode. A new adaptive mutation causes the frequency of the linked barcode to grow, which can be detected by sequencing. By sampling many random mutants and genotyping them at the barcode locus, one can identify mutants from adapted lineages even if they are rare (Venkataram et al., 2016). As a result, BLT allows one to sample mutants soon after they acquire their first driver mutation, before acquiring secondary mutations.
To evaluate the quality of sampling based on the LD and BLT methods, we consider the following hypothetical experimental setup. $K$ beneficial mutants are sampled from each home environment (with either one of the methods), and their home and nonhome fitness $({X}_{i},{Y}_{i})$ are measured for each mutant $i=1,\mathrm{\dots},K$. Since we are ultimately interested in ranking drug pairs by their risk of collateral resistance, we estimate the collateral risk statistic $\widehat{c}$ from these fitness data for each drug pair and use $\widehat{c}$ to rank them (see Materials and methods for details). We compare such a priori ranking of 125 hypothetical drug pairs with Gaussian JDFEs used in previous sections with their a posteriori ranking based on the risk of collateral resistance observed in simulations.
To model the LD sampling method on a given JDFE, we randomly sample $K$ mutants whose home fitness exceeds a certain cutoff. To model a BLT experiment, we simulate evolution in the home environment and randomly sample $K$ beneficial mutants segregating at generation 250 (see Materials and methods for details). We find that the $\widehat{c}$ranking estimated with either LD or BLT methods captures the a posteriori ranking surprisingly well, even when the number of sampled mutants is as low as 10 per drug pair (Figure 5). Given that the JDFEs with adjacent ranks differ in $c$ by a median of only 0.65%, the strong correlations shown in Figure 5 suggest that even very similar JDFEs can be differentiated with moderate sample sizes. As expected, this correlation is further improved upon increased sampling, and it is insensitive to the specific home fitness threshold that we use in the LD method (Figure 5—figure supplement 1). We conclude that estimating JDFE parameters is in principle feasible with a modest experimental effort, at least for the purpose of ranking drug pairs.
Discussion
We have shown that many resistance mutations against multiple drugs in E. coli exhibit a diversity of collateral effects. If this is true more generally, it implies that there is an unavoidable uncertainty in whether any given population would evolve collateral resistance or sensitivity, which could at least in part explain previously observed inconsistencies among experiments. We quantified the diversity of pleiotropic effects of mutations with a joint distribution of fitness effects (JDFE) and developed a population genetic theory for predicting the expected collateral outcomes of evolution and the uncertainty around these expectations. In the successional mutations regime, our theory shows that the average rate at which fitness in the nonhome environment is gained or lost during adaptation to the home environment is determined by the pleiotropy statistic $r}_{2$ given by Equation 8. How strongly the nonhome fitness in any individual population deviates from this ensemble average is determined by the pleiotropy variance statistic $D}_{22$ given by Equation 11. Importantly, $r}_{2$ and $D}_{22$ are properties of the JDFE alone, that is, they do not depend on the parameters of any specific population. In the concurrent mutations regime, the expected rate of nonhome fitness gain or loss and the associated variance are reasonably well predicted by the adjusted pleiotropy statistics ${r}_{2}^{*}$ and ${D}_{22}^{*}$. Unlike $r}_{2$ and $D}_{22$, the adjusted statistics depend on the population size $N$ and the rate of beneficial mutations ${U}_{b}$.
To quantitatively predict the probability of evolution of collateral drug resistance in practice would require the knowledge of both the JDFE for the focal bacterial or cancercell population in the presence of the specific pair of drugs and its in vivo population genetic parameters. Since estimating the latter parameters is very difficult, it appears unlikely that we would be able to quantitatively predict the dynamics of collateral effects, even if JDFEs were known. A more realistic application of our theory is that it allows us to rank drug pairs according to the risk of collateral resistance even when the population genetic parameters are unknown. Such robust ranking can be computed based on the collateral risk statistic $c={r}_{2}/\sqrt{{D}_{22}}$, a property of the JDFE but not of the evolving population. Drug pairs with positive values of $c$ have a higher chance of eliciting collateral resistance than collateral sensitivity and should be avoided; drug pairs with more negative values of $c$ have a lower risk of collateral resistance evolution than those with less negative values.
We have validated our theory in silico, but how well it would work in vivo (in the clinic) or even in vitro (in the lab) is as of yet unclear. A direct way to validate the theory empirically would be to estimate JDFEs for a model organism, such as E. coli, in a number of drug pairs, rank these pairs according to our collateral rank statistic and then test this ranking by evolving replicate populations and measuring the empirical distributions of collateral resistance/sensitivity outcomes. To the best of our knowledge, the antibiotic resistance JDFEs among genomewide mutations have not yet been measured. One could in principle use existing gene knockout data, such as those obtained by Chevereau et al., 2015 (Figure 1), or the data from deep mutational scanning experiments, such as those obtained by Stiffler et al., 2015 (Figure 1—figure supplement 1), to estimate JDFEs. However, these experiments estimate fitness only for certain subsets of mutations (gene knockouts or point mutations within a single gene, respectively). Since resistance may arise via other types of mutations (Nichol et al., 2019), these data would give us at best an incomplete picture of actual JDFEs. Our results suggest that JDFEs can be reasonably well estimated by sampling resistance mutants at drug concentrations above MIC or by employing the barcode lineage tracking method.
Another obstacle is that, even though many researchers have experimentally evolved various microbes in the presence of drugs, most experiments have maintained too few replicate populations to accurately measure the variation in collateral outcomes of evolution. The study by Nichol et al., 2019, with 60 replicates, is a notable exception. In short, a rigorous test of our theory requires new data on the shapes of wholegenome JDFEs as well as higher throughput evolution experiments.
What the most effective ways of measuring JDFEs are and whether it will be possible to measure JDFE in vivo are open questions. We speculate that the answers will depend on the shapes of the empirical JDFEs because some shapes may be more difficult to estimate than others. For example, if empirical JDFEs resemble multivariate Gaussian distributions, then we can learn all relevant parameters of such JDFE by sampling a handful of random mutants and measuring their fitness in relevant environments. One can also imagine more complex JDFEs where mutations beneficial in the home environment have a dramatically different distribution of nonhome fitness effects than mutations that are deleterious or neutral in the home environment. In this case, very large samples of random mutations would be necessary to correctly predict the pleiotropic outcomes of evolution, so that methods that preferentially sample beneficial mutations may be required. We have considered two such methods, which are experimentally feasible. We have shown that both of them perform extremely well on Gaussian JDFEs in the sense that as few as 10 mutants per drug pair are sufficient to produce largely correct ranking of hypothetical drug pairs. However, it may be difficult to apply these methods in vivo, in which case JDFEs may have to be estimated in the lab, with selection pressures reproducing those in vivo as accurately as possible.
Our model relies on two important simplifications. It describes the evolution of an asexual population where all resistance alleles arise from de novo mutations. In reality, some resistance alleles in bacteria may be transferred horizontally (Sun et al., 2019). Understanding collateral resistance evolution in the presence of horizontal gene transfer events would require incorporating JDFE into other models of evolutionary dynamics (e.g. Neher et al., 2010). Another major simplification is in the assumption that the JDFE stays constant as the population adapts. In reality the JDFE will change over time because of the depletion of the pool of adaptive mutations and because of epistasis (Good et al., 2017; Venkataram et al., 2020). How JDFEs vary among genetic backgrounds is currently unknown. In Appendix 1, we have shown that our main results hold at least in the presence of a simple form of ‘global’ epistasis. Empirically measuring how JDFEs vary across genotypes and theoretically understanding how such variation affects the evolution of pleiotropic outcomes are important open questions.
While we were primarily motivated by the problem of evolution of collateral drug resistance and sensitivity, our theory is applicable more broadly. The shape of JDFE must play a crucial role in determining whether the population evolves toward a generalist or diversifies into multiple specialist ecotypes. Previous literature has viewed this question primarily through the lense of two alternative hypotheses: antagonistic pleiotropy and mutation accumulation (Visher and Boots, 2020). Antagonistic pleiotropy in its strictest sense means that the population is at the Pareto front with respect to the home and nonhome fitness, such that any mutation beneficial in the home environment reduces the fitness in the nonhome environment (Li et al., 2019). The shape of the Pareto front then determines whether selection would favor specialists or generalists (Levins, 1968; Visher and Boots, 2020). Alternatively, a population can evolve to become a homeenvironment specialist even in the absence of tradeoffs, simply by accumulating mutations that are neutral in the home environment but deleterious in the nonhome environment (Kawecki, 1994). More recently, it has been recognized that antagonistic pleiotropy and mutation accumulation are not discrete alternatives but rather extremes of a continuum of models (Bono et al., 2020; Jerison et al., 2014; Jerison et al., 2020). The JDFE provides a mathematical way to describe this continuum. For example, strict antagonistic pleiotropy can be modeled with a JDFE with zero probability weight in the first quadrant and a bulk of probability in the fourth quadrant. A mutation accumulation scenario can be modeled with a ‘+’like JDFE where all mutations beneficial in the home environment are neutral in the nonhome environment (i.e. concentrated on the $x$axis) and all or most mutations neutral in the home environment (i.e. those on the $y$axis) are deleterious in the nonhome environment. Our theory shows that in fact all JDFEs with negative r_{2} lead to loss of fitness in the nonhome environment and therefore can potentially promote specialization. While our theory provides this insight, further work is needed to understand how JDFEs govern adaptation to variable environments. This future theoretical work, together with empirical inquiries into the shapes of JDFEs, will not only advance our ability to predict evolution in practical situations, such as drug resistance, but it will also help us better understand the origins of ecological diversity.
Materials and methods
Analysis of knockout and deep mutational scanning data
Knockout data
Request a detailed protocolChevereau et al., 2015 provide growth rate estimates for 3883 gene knockout mutants of E. coli in the presence of six antibiotics. Our goal is to identify those knockout mutations that provide resistance against one drug and are also collaterally resistant or collaterally sensitive to another drug. However, it is unclear from these original data alone which mutations have statistically significant beneficial and deleterious effects because no measurement noise estimates are provided. To address this problem, we obtained replicate wildtype growth rate measurements in the presence of antibiotics from Guillaume Chevereau and Tobias Bollenbach (available at https://github.com/ardellsarah/JDFEproject; copy archived at swh:1:rev:e91f2940681269511c6bb9fd4560ccd4a7c4d641, Ardell, 2022). In this additional data set, the wildtype E. coli strain is measured on average 476 times in the presence of each drug. We estimate the wildtype growth rate ${r}_{\mathrm{WT}}$ as the mean of these measurements, and we obtain the selection coefficient for all knockout mutants as ${s}_{i}={r}_{i}{r}_{\mathrm{WT}}$. We also obtain the noise distribution ${P}_{\mathrm{noise}}(s)$ from the replicate wildtype measurements (shown in red in the diagonal panels in Figure 1). Modeling ${P}_{\mathrm{noise}}(s)$ as normal distributions, we obtain the Pvalues for each mutation in the presence of each antibiotic.
We then call any knockout mutant as resistant against a given drug if its selection coefficient in the presence of that drug exceeds a critical value ${s}_{\alpha}^{+}\S gt;0$. We choose ${s}_{\alpha}^{+}$ using the BenjaminiHochberg procedure (Benjamini and Hochberg, 1995) so that the false discovery rate (FDR) among the identified resistant mutants is $\alpha \approx 0.25$. We could not find an ${s}_{\alpha}^{+}$ for $\alpha \lesssim 0.25$ for trimethoprim (TMP) and chloramphenicol (CHL), that is, there were not enough knockout mutations with positive selection coefficients to reliably distinguish them from measurement errors.
We apply the same procedure to identify mutations that are collaterally resistant and collaterally sensitive against a second drug among all mutations that are resistant against the first drug, except we aim for FDR $\lesssim 0.10$.
Deep mutational scanning data
Request a detailed protocolStiffler et al., 2015 provide estimates of relative fitness for 4997 point mutations in the TEM1 βlactamase gene in the presence of cefotaxime (CTX) and four concentrations of ampicillin (AMP). We estimate the selection coefficients from the reported relative fitness values by changing the logarithm from log_{10} to natural and dividing it by six, the estimated number of generations that occurred during the 2hr experiment. The latter is based on the fact that Stiffler et al., 2015 chose AMP concentrations which did not significantly alter the E. coli doubling time, which we assumed to be 20 minutes. We used the same number of generations for CTX.
Stiffler et al., 2015 report two replicate measurements per mutant in each concentration of AMP and one measurement per mutant in the presence of CTX. We consider CTX as the home environment and call all mutations with positive measured fitness effects as resistant against CTX. For each such mutation, we use two replicate measurements in each concentration of AMP to estimate its mean fitness effect and the 90% confidence interval around the mean, based on the normal distribution. We call any CTX resistant mutation with the entire confidence interval above (below) zero as collaterally resistant (sensitive) against AMP at that concentration. All remaining CTX resistant mutations are called collaterally neutral.
Theory
Successional mutations regime
Request a detailed protocolWe assume that an asexual population evolves according the WrightFisher model in the strong selection weak mutation (SSWM) limit (Orr, 2000; Kryazhimskiy et al., 2009; Good and Desai, 2015), also known as the ‘successional mutations’ regime (Desai and Fisher, 2007). In this regime, the population remains monomorphic until the arrival of a new mutation that is destined to fix. The waiting time for such new mutation is assumed to be much longer than the time it takes for the mutation to fix, that is, fixation happens almost instantaneously on this time scale, after which point the population is again monomorphic. If the per genome per generation rate of beneficial mutations is $U}_{b$, their typical effect is $s$ and the population size is $N$, the SSWM approximation holds when $N{U}_{b}\ll 1/\mathrm{ln}\left(Ns\right)$ (Desai and Fisher, 2007).
We describe our population by a twodimensional vector of random variables $({X}_{t},{Y}_{t})$, where ${X}_{t}$ and ${Y}_{t}$ are the population’s fitness (growth rate or the Malthusian parameter) in the home and nonhome environments at generation $t$, respectively. We assume that the fitness vector of the population at the initial time point is known and is $({x}_{0},{y}_{0})$. We are interested in characterizing the joint probability density $p(x,y,t)\phantom{\rule{thinmathspace}{0ex}}dx\phantom{\rule{thinmathspace}{0ex}}dy=\mathrm{P}\mathrm{r}\left\{{X}_{t}\in (x,x+dx),{Y}_{t}\in (y,y+dy)\right\}$.
We assume that all genotypes have the same JDFE $\mathrm{\Phi}(\mathrm{\Delta}x,\mathrm{\Delta}y)$, that is, there is no epistasis. In the exponential growth model, the selection coefficient of a mutation is the difference between the mutant and the ancestor growth rates in the home environment, that is, $\mathrm{\Delta}x$. The probability of fixation of the mutant is given by Kimura’s formula, which we approximate by $2\mathrm{\Delta}x$ for $\mathrm{\Delta}x\S gt;0$ and zero otherwise (Crow and Kimura, 1972).
If the total rate of mutations (per genome per generation) is $U$, the rate of mutations beneficial in the home environment is given by $U}_{b}=U{f}_{b$ where ${f}_{b}={\int}_{\mathrm{\infty}}^{\mathrm{\infty}}d\eta {\int}_{0}^{\mathrm{\infty}}d\xi \phantom{\rule{thickmathspace}{0ex}}\mathrm{\Phi}(\xi ,\eta )$ is the fraction of mutations beneficial in the home environment. Once such a mutation arises, its selection coefficients in the home and nonhome environments are drawn from the JDFE of mutations beneficial in the home environment ${\mathrm{\Phi}}_{b}(\mathrm{\Delta}x,\mathrm{\Delta}y)=\mathrm{\Phi}(\mathrm{\Delta}x,\mathrm{\Delta}y)/{f}_{b}$. Then, in the SSWM limit, our population is described by a twodimensional continuoustime continuousspace Markov chain with the transition rate from state $(x,y)$ to state $({x}^{\prime},{y}^{\prime})$ given by
The probability distribution $p(x,y,t)$ satisfies the integrodifferential forward Kolmogorov equation (Van Kampen, 1992)
with the initial condition
When beneficial mutations with large effects are sufficiently rare, Equation 4 can be approximated by the FokkerPlanck equation (Van Kampen, 1992)
where
are the expected fitness effects in the home and nonhome environments for a mutation fixed in the home environment, and
are the second moments of the distribution of the fitness effects of mutations fixed in the home environment. The solution to Equation 6 with the initial condition given by Equation 5 is a multivariate normal distribution with the mean vector $\mathit{m}(t)$ and the variancecovariance matrix ${\sigma}^{2}(t)$ given by Equation 1 and Equation 2.
Concurrent mutations regime
Request a detailed protocolThe theory we developed so far for the successional mutations regime breaks down in the concurrent mutations regime, that is, when multiple adaptive mutations segregate in the population simultaneously (Desai and Fisher, 2007). The main effect of competition between segregating adaptive lineages is that many new beneficial mutations arise in relatively lowfitness genetic backgrounds and have almost no chance of surviving competition (Desai and Fisher, 2007; Schiffels et al., 2011; Good et al., 2012). As a result, the fixation probability of a beneficial mutation with selective effect $\mathrm{\Delta}x$ in the home environment is no longer $2\mathrm{\Delta}x$. Instead, beneficial mutations that provide fitness benefits below a certain threshold $x}_{c$ behave as if they are effectively neutral (i.e. their fixation probability is close to zero), and most adaptation is driven by mutations with benefits above $x}_{c$, where $x}_{c$ depends on the population genetic parameters $N$ and ${U}_{b}$ as well as the shape of the distribution of fitness effects of beneficial mutations. Good et al., 2012 derived equations that allow us to calculate the effective fixation probability ${\pi}^{*}(\mathrm{\Delta}x;N,{U}_{b})$ of a beneficial mutation with the fitness benefit $\mathrm{\Delta}x$ in the home environment in the concurrent mutation regime. Thus, to predict the average rate of nonhome fitness change, we replace the SSWM fixation probability $2\xi $ in Equation 8 with ${\pi}^{*}(\xi ;N,{U}_{b})$ and obtain the adjusted expected pleiotropic effect. We similarly obtain the adjusted pleiotropic variance statistic
although as discussed in Section ‘The population genetics of pleiotropy’, we do not expect ${D}_{22}^{*}$ to capture all of the variation in nonhome fitness trajectories.
To calculate ${\pi}^{*}(\mathrm{\Delta}x;N,{U}_{b})$ for the Gaussian JDFEs shown in Figure 2, we first substitute Equation (20) in Good et al., 2012 with $\beta =2$ into Equation 18, 19 in Good et al., 2012 and then numerically solve these equations for $x}_{c$ and $v$ using the FindRoot numerical method in Mathematica. Note that all our Gaussian JDFEs share the same mean and variance in the home environment, so we need to solve these equations only once for each pair of $N$ and ${U}_{b}$ values. We then substitute the obtained values of $x}_{c$ and $v$ into Equation (4) (9) in Good et al., 2012 and calculate ${\pi}^{*}$ by a numerical integration of Equation (6) in Good et al., 2012 in R (available at https://github.com/ardellsarah/JDFEproject).
Ranking of drug pairs
Request a detailed protocolAccording to Equation 1 and Equation 2, both the expected nonhome fitness and its variance change linearly with time, so that at time $t$ the mean is $Z=c\sqrt{N{U}_{b}t}$ standard deviations above $y}_{0$ (if ${r}_{2}\S gt;0$) or below $y}_{0$ (if ${r}_{2}\S lt;0$), where $c={r}_{2}/\sqrt{{D}_{22}}$. In other words, if ${r}_{2}\S gt;0$, the bulk of the nonhome fitness distribution eventually shifts above $y}_{0$, and if ${r}_{2}\S lt;0$, it shifts below $y}_{0$. All else being equal, a larger value of $c$ implies faster rate of this shift.
The interpretation of these observations in terms of collateral resistance/sensitivity is that adaptation in the presence of the first drug will eventually lead to collateral resistance against the second drug if ${r}_{2}\S gt;0$ and to collateral sensitivity if ${r}_{2}\S lt;0$. Furthermore, all else being equal, collateral sensitivity evolves faster and the chance of evolving collateral resistance is smaller for drug pairs with more negative $c$ (i.e. larger $c$). Thus, we use $c$ to order drug pairs from the most preferred (those with the most negative values of $c$) to least preferred (those with least negative or positive values of $c$).
Generation of JDFEs
Gaussian JDFEs
Request a detailed protocolThe JDFEs in Figure 2 have the following parameters. Mean in the home environment: 0.05. Standard deviation in both home and nonhome environments: 0.1. Means in the nonhome environment: 0.08, 0.145, 0, 0.145, 0.08 in panels A through E, respectively.
The JDFEs in Figure 3 have the following parameters. Mean and standard deviation in the home environment: 0.001 and 0.001, respectively. The nonhome mean varies between 0.0001 and 0.01. The nonhome standard deviation varies between 0.0001 and 0.01. The correlation between home and nonhome fitness varies between 0.9 and 0.9, for a total of 125 JDFEs. All parameter values and the resulting pleiotropy statistics for these JDFEs are given in the Figure 3—source data 1.
JDFEs with equal probabilities of pleiotropically beneficial and deleterious mutations
Request a detailed protocolAll JDFEs in Figure 2—figure supplement 1 are mixtures of two twodimensional uncorrelated Gaussian distributions, which have the following parameters. Mean in the home environment: 0.4. Standard deviation in both home and nonhome environments: 0.1. Means in the nonhome environment: 0.1 and 0.1 in panel A, 0.5 and 0.5 in panel B, 0.17 and 0.5 in panel C, and 0.5 and 0.17 in panel D.
Simulations
We carried out two types of simulations, SSWM model simulations and full WrightFisher model simulations.
Strong selection weak mutation
Request a detailed protocolThe SSWM simulations were carried out using the Gillespie algorithm (Gillespie, 1976), as follows. We initiate the populations with home and nonhome fitness values ${x}_{0}=0$ and ${y}_{0}=0$. At each iteration, we draw the waiting time until the appearance of the next beneficial mutation from the exponential distribution with the rate parameter $N{U}_{b}$ and advance the time by this amount. Then, we draw the selection coefficients $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$ of this mutation in the home and nonhome environment, respectively, from the JDFE (a multivariate normal distribution). With probability $2\mathrm{\Delta}x$, the mutation fixes in the population. If it does, the fitness of the population is updated accordingly.
WrightFisher model
Request a detailed protocolWe simulate evolution in the home environment according to the WrightFisher model with population size $N$ as follows. We initiate the whole population with a single genotype with fitness ${x}_{0}=0$ and ${y}_{0}=0$ in the home and nonhome environments. Suppose that at generation $t$, there are $K(t)$ genotypes, such that genotype $i$ has home and nonhome fitness ${X}_{i}$ and ${Y}_{i}$, respectively, and it is present at frequency ${f}_{i}(t)\S gt;0$ in the population. We generate the genotype frequencies at generation $t+1$ in three steps. In the reproduction step, we draw random numbers ${B}_{i}^{\prime}(t+1)$, $i=1,\mathrm{\dots},K(t)$ from the multinomial distribution with the number of trials $N$ and success probabilities ${p}_{i}(t)={f}_{i}(t)+{f}_{i}(t)\left({X}_{i}(t)\overline{X}(t)\right),$ where $\overline{X}(t)={\sum}_{i=1}^{K(t)}{X}_{i}(t){f}_{i}(t)$ is the mean fitness of the population in the home environment at generation $t$. In the mutation step, we draw a random number $M$ of new mutants from the Poisson distribution with parameter $NU$, where $U$ is the total per individual per generation mutation rate. We randomly determine the ‘parent’ genotypes in which each mutation occurs and turn the appropriate numbers of parent individuals into new mutants. We assume that each new mutation creates a new genotype and has fitness effects $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$ in the home and nonhome environments. $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$ are drawn randomly from the JDFE $\mathrm{\Phi}(\mathrm{\Delta}x,\mathrm{\Delta}y)$. We obtain each mutants fitness by adding these values to the parent genotype’s home and nonhome fitness values. In the final step, all genotypes that are represented by zero individuals are removed and we are left with $K(t+1)$ genotypes with ${B}_{i}(t+1)\S gt;0$ individuals, $i=1,\mathrm{\dots},K(t+1)$. Then we set ${f}_{i}(t+1)={B}_{i}(t+1)/N$.
Sampling beneficial mutants from JDFEs and estimating the collateral risk statistic
Request a detailed protocolWe model the LD sampling method by randomly drawing mutants from the JDFE until the desired number $K$ of mutants whose home fitness exceeds the focal threshold are sampled. We estimate the $c$ statistic from the pairs of home and nonhome fitness effects ${X}_{i}$ and ${Y}_{i}$ of these $i=1,\mathrm{\dots},K$ sampled mutants. To do so, we first estimate $r}_{2$ and $D}_{22$ as ${\widehat{r}}_{2}=1/K{\sum}_{i=1}^{K}{X}_{i}{Y}_{i}$ and ${\widehat{D}}_{22}=1/K{\sum}_{i=1}^{K}{X}_{i}{Y}_{i}^{2}$. We then calculate $\widehat{c}={\widehat{r}}_{2}/\sqrt{{\widehat{D}}_{22}}$.
For the BLT sampling method, we simulate the WrightFisher model as described above for $N={10}^{6}$ and $U={10}^{4}$ for 250 generations. At generation 250, we randomly sample existing beneficial mutants proportional to their frequency in the population without replacement (i.e. the same beneficial mutation is sampled at most once). Sampling more than $\sim 50$ distinct beneficial mutants from a single population becomes difficult because there may simply be not enough such mutants or some of them may be at very low frequencies. Therefore, if the desired number of mutants to sample exceeds 50, we run multiple replicate simulations and sample a maximum of 100 distinct beneficial mutants per replicate until the desired number of mutants is reached. We then estimate the $c$ statistics as with the LD method.
Appendix 1
JDFE with global epistasis
Results in the main text were derived under the assumption that all genotypes have the same JDFE, i.e., in the absence of epistasis. In reality, JDFEs probably vary from one genotype to another, but how they vary is not yet well characterized. Recent studies have found that the fitness effects of many mutations available to a genotype in a given environment depend primarily on the fitness of that genotype in that environment (Khan et al., 2011; Chou et al., 2011; Wiser et al., 2013; Kryazhimskiy et al., 2014; Johnson et al., 2019; Wang et al., 2016; Aggeli et al., 2021; Lukačišinová et al., 2020). This dependence is sometimes referred to as global or fitnessdependent epistasis (Kryazhimskiy et al., 2009; Kryazhimskiy et al., 2014; Reddy and Desai, 2021; Husain and Murugan, 2020). Here, we ask whether our main results would hold if the pathogen population evolves on a JDFE with global epistasis.
Global epistasis can be modeled in our framework by assuming that the JDFE ${\mathrm{\Phi}}_{g}$ of genotype $g$ depends only the fitness of this genotype in the home and nonhome environments, $x(g)$, $y(g)$, i.e. ${\mathrm{\Phi}}_{g}(\mathrm{\Delta}x,\mathrm{\Delta}y)={\mathrm{\Phi}}_{x(g),y(g)}(\mathrm{\Delta}x,\mathrm{\Delta}y)$, which is a twodimensional extension of the model considered by Kryazhimskiy et al., 2009. Thus, in the SSWM regime, the population can still be fully described by its current pair of fitness values in the home and nonhome environments $({X}_{t},{Y}_{t})$. The dynamics of the probability density $p(x,y,t)$ are governed by the same Kolmogorov equation as in the nonepistatic case, which can still be approximated by a diffusion equation (Equation 6). However, while in the nonepistatic case the drift and diffusion coefficients of this equation, r_{1}, r_{2}, D_{11}, D_{12} and D_{22} are constants, in the presence of global epistasis, they become functions of $x$ and $y$. Although this equation cannot be solved analytically in the general case, it can be solved numerically, provided that the functions ${r}_{1}(x,y)$, ${r}_{2}(x,y)$, ${D}_{11}(x,y)$, ${D}_{12}(x,y)$ and ${D}_{22}(x,y)$ are known. Thus, in principle, our theory can predict the trajectories of nonhome fitness in the presence of global epistasis.
To explore the implications of global epistasis for collateral drug resistance evolution, we consider the simplest scenario where the functional form of global epistasis (i.e., how ${\mathrm{\Phi}}_{x,y}$ depends on $x$ and $y$) is the same across different drugs. In this case, we would expect that the ranking of drug pairs according to the risk of collateral resistance would be the same for all genotypes. In particular, the drug pair whose risk of collateral resistance risk is the lowest for the wildtype should also be the pair with the lowest risk for the evolved genotypes.
To test this prediction, we model resistance evolution on Gaussian JDFEs whose mean vector and the correlation coefficient are fixed while the standard deviations and in the home and nonhome environments decrease linearly with the fitness in the respective environment, $\sigma}_{h}(x)=max\{0,{\sigma}_{h,0}{\gamma}_{h}x\$ and $\sigma}_{nh}(y)=max\{0,{\sigma}_{nh,0}{\gamma}_{nh}y\$.
Appendix 1—figure 1A shows how one such JDFE changes along the expected evolutionary trajectory. The corresponding expected home and nonhome fitness trajectories and their variance are shown in Appendix 1—figure 1B. Appendix 1—figure 1C shows how the probability (risk) of collateral resistance changes over time on four different JDFEs with global epistasis. For the ancestral strain (whose fitness we set by convention to $x=y=0$), these four JDFEs are identical to those shown in Figure 4A; as the populations evolve, JDFEs change as specified above with ${\gamma}_{\mathrm{h}}={\gamma}_{\mathrm{nh}}=0.5$. As expected, the ranking of these epistatic JDFEs according to the risk of collateral resistance stays constant over time and can be predicted from estimates of the $c$ parameters for the ancestral strain.
To test this prediction, we model resistance evolution on Gaussian JDFEs whose mean vector and the correlation coefficient are fixed while the standard deviations and in the home and nonhome environments decrease linearly with the fitness in the respective environment, $\sigma}_{h}(x)=max\{0,{\sigma}_{h,0}{\gamma}_{h}x\$ and $\sigma}_{nh}(y)=max\{0,{\sigma}_{nh,0}{\gamma}_{nh}y\$.Appendix 1—figure 1A shows how one such JDFE changes along the expected evolutionary trajectory. The corresponding expected home and nonhome fitness trajectories and their variance are shown in Appendix 1—figure 1B. Appendix 1—figure 1C shows how the probability (risk) of collateral resistance changes over time on four different JDFEs with global epistasis. For the ancestral strain (whose fitness we set by convention to $x=y=0$), these four JDFEs are identical to those shown in Figure 4A; as the populations evolve, JDFEs change as specified above with ${\gamma}_{\mathrm{h}}={\gamma}_{\mathrm{nh}}=0.5$. As expected, the ranking of these epistatic JDFEs according to the risk of collateral resistance stays constant over time and can be predicted from estimates of the $c$ parameters for the ancestral strain.
To test this prediction, we model resistance evolution on Gaussian JDFEs whose mean vector and the correlation coefficient are fixed while the standard deviations and in the home and nonhome environments decrease linearly with the fitness in the respective environment, $\sigma}_{h}(x)=max\{0,{\sigma}_{h,0}{\gamma}_{h}x\$ and $\sigma}_{nh}(y)=max\{0,{\sigma}_{nh,0}{\gamma}_{nh}y\$.Appendix 1—figure 1A shows how one such JDFE changes along the expected evolutionary trajectory. The corresponding expected home and nonhome fitness trajectories and their variance are shown in Appendix 1—figure 1B. Appendix 1—figure 1C shows how the probability (risk) of collateral resistance changes over time on four different JDFEs with global epistasis. For the ancestral strain (whose fitness we set by convention to $x=y=0$), these four JDFEs are identical to those shown in Figure 4A; as the populations evolve, JDFEs change as specified above with ${\gamma}_{\mathrm{h}}={\gamma}_{\mathrm{nh}}=0.5$. As expected, the ranking of these epistatic JDFEs according to the risk of collateral resistance stays constant over time and can be predicted from estimates of the $c$ parameters for the ancestral strain.
Data availability
All code is available on GitHub (https://github.com/ardellsarah/JDFEproject; copy archived at swh:1:rev:e91f2940681269511c6bb9fd4560ccd4a7c4d641). All data are available as Source Data files, included with the manuscript.
References

Microbiological effects of sublethal levels of antibioticsNature Reviews. Microbiology 12:465–478.https://doi.org/10.1038/nrmicro3270

SoftwareJDFEproject, version swh:1:rev:e91f2940681269511c6bb9fd4560ccd4a7c4d641Software Heritage.

Alternative Evolutionary Paths to Bacterial Antibiotic Resistance Cause Distinct Collateral EffectsMolecular Biology and Evolution 34:2229–2244.https://doi.org/10.1093/molbev/msx158

Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple TestingJournal of the Royal Statistical Society 57:289–300.https://doi.org/10.1111/j.25176161.1995.tb02031.x

The dynamics of adaptive genetic diversity during the early stages of clonal evolutionNature Ecology & Evolution 3:293–301.https://doi.org/10.1038/s4155901807581

Evolvability Costs of Niche ExpansionTrends in Genetics 36:14–23.https://doi.org/10.1016/j.tig.2019.10.003

The distribution of fitness effects in an uncertain worldEvolution; International Journal of Organic Evolution 69:1610–1618.https://doi.org/10.1111/evo.12673

The distribution of fitness effects of new mutationsNature Reviews. Genetics 8:610–618.https://doi.org/10.1038/nrg2146

A general method for numerically simulating the stochastic time evolution of coupled chemical reactionsJournal of Computational Physics 22:403–434.https://doi.org/10.1016/00219991(76)900413

Fisher’s geometrical model and the mutational patterns of antibiotic resistance across dose gradientsEvolution; International Journal of Organic Evolution 71:23–37.https://doi.org/10.1111/evo.13111

Physical Constraints on EpistasisMolecular Biology and Evolution 37:2865–2874.https://doi.org/10.1093/molbev/msaa124

Use of collateral sensitivity networks to design drug cycling protocols that avoid resistance developmentScience Translational Medicine 5:204ra132.https://doi.org/10.1126/scitranslmed.3006609

Compatibility of Evolutionary Responses to Constituent Antibiotics Drive Resistance Evolution to Drug PairsMolecular Biology and Evolution 38:2057–2069.https://doi.org/10.1093/molbev/msab006

Chance and necessity in the pleiotropic consequences of adaptation for budding yeastNature Ecology & Evolution 4:601–611.https://doi.org/10.1038/s4155902011283

Theoretical models of selection and mutation on quantitative traitsPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 360:1411–1425.https://doi.org/10.1098/rstb.2005.1667

Stability of the Gmatrix in a population experiencing pleiotropic mutation, stabilizing selection, and genetic driftEvolution; International Journal of Organic Evolution 57:1747–1760.https://doi.org/10.1111/j.00143820.2003.tb00583.x

Reciprocal antibiotic collateral sensitivity in burkholderia multivoransInternational Journal of Antimicrobial Agents 56:105994.https://doi.org/10.1016/j.ijantimicag.2020.105994

Accumulation of Deleterious Mutations and the Evolutionary Cost of Being a GeneralistThe American Naturalist 144:833–838.https://doi.org/10.1086/285709

ConferenceProceedings of the Sixth Berkeley Symposium on Mathematical Statistics and ProbabilityThe role of mutation in evolution.

The measurement of selection on correlated charactersEvolution; International Journal of Organic Evolution 37:1210–1226.https://doi.org/10.1111/j.15585646.1983.tb00236.x

Bacterial evolution of antibiotic hypersensitivityMolecular Systems Biology 9:700.https://doi.org/10.1038/msb.2013.57

BookEvolution in Changing EnvironmentsPrinceton University Press.https://doi.org/10.1515/9780691209418

Single nucleotide mapping of trait space reveals Pareto fronts that constrain adaptationNature Ecology & Evolution 3:1539–1551.https://doi.org/10.1038/s4155901909930

Evolution in alternating environments with tunable interlandscape correlationsEvolution; International Journal of Organic Evolution 75:10–24.https://doi.org/10.1111/evo.14121

The fitness effect of mutations across environments: Fisher’s geometrical model with multiple optimaEvolution; International Journal of Organic Evolution 69:1433–1447.https://doi.org/10.1111/evo.12671

Adaptive Landscapes of Resistance Genes Change as Antibiotic Concentrations ChangeMolecular Biology and Evolution 32:2707–2715.https://doi.org/10.1093/molbev/msv146

Prediction of resistance development against drug combinations by collateral responses to component drugsScience Translational Medicine 6:262ra156.https://doi.org/10.1126/scitranslmed.3009940

Very slightly deleterious mutations and the molecular clockJournal of Molecular Evolution 26:1–6.https://doi.org/10.1007/BF02111276

Strength of selection pressure is an important parameter contributing to the complexity of antibiotic resistance evolutionMolecular Biology and Evolution 31:2387–2401.https://doi.org/10.1093/molbev/msu191

Collateral sensitivity of antibioticresistant microbesTrends in Microbiology 23:401–407.https://doi.org/10.1016/j.tim.2015.02.009

Metabolic fitness landscapes predict the evolution of antibiotic resistanceNature Ecology & Evolution 5:677–687.https://doi.org/10.1038/s41559021013970

Collateral sensitivity as a strategy against cancer multidrug resistanceDrug Resistance Updates 15:98–105.https://doi.org/10.1016/j.drup.2012.03.002

The evolution of tradeoffs: where are we?Journal of Evolutionary Biology 20:433–447.https://doi.org/10.1111/j.14209101.2006.01255.x

Evolutionary landscapes of Pseudomonas aeruginosa towards ribosometargeting antibiotic resistance depend on selection strengthInternational Journal of Antimicrobial Agents 55:105965.https://doi.org/10.1016/j.ijantimicag.2020.105965

Role of pleiotropy during adaptation of TEM1 βlactamase to two novel antibioticsEvolutionary Applications 8:248–260.https://doi.org/10.1111/eva.12200

Editorial: Horizontal Gene Transfer Mediated Bacterial Antibiotic ResistanceFrontiers in Microbiology 10:1933.https://doi.org/10.3389/fmicb.2019.01933

The problem of mediocre generalists: population genetics and ecoevolutionary perspectives on host breadth evolution in pathogensProceedings of the Royal Society B 287:20201230.https://doi.org/10.1098/rspb.2020.1230

The pleiotropic structure of the genotypephenotype map: the evolvability of complex organismsNature Reviews. Genetics 12:204–213.https://doi.org/10.1038/nrg2949

Evolving generalists in switching rugged landscapesPLOS Computational Biology 15:e1007320.https://doi.org/10.1371/journal.pcbi.1007320
Decision letter

Richard A NeherReviewing Editor; University of Basel, Switzerland

Meredith C SchumanSenior Editor; University of Zurich, Switzerland

Joachim KrugReviewer; University of Cologne, Germany
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]
Thank you for submitting your work entitled "The population genetics of pleiotropy, and the evolution of collateral resistance and sensitivity in bacteria" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Richard A Neher as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Joachim Krug (Reviewer #2).
Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife. However, if you feel you can address the comments as outlined below, we encourage you to resubmit.
The reviewers agreed that pleiotropy of mutations and the resulting adaptive trajectories across different environment are important topics that are both of theoretical and applied interest. Your theoretical framework predicts fitness trajectories observed in environments other than the one selection is operating in (home env). These trajectories in nonhome environments are calculated via integrals over the joint fitness effect distribution weighted by the fixation probability in the home environment. However, your framework assumes strong selection and weak mutation (SSWM) and deviations from this assumption seem to have strong effects. We think that these effects need to be at least partially understood. Furthermore, application to the KO library is a useful proofofconcept, but the practical relevance of these patterns for understanding collateral sensitivity/resistance is far from obvious.
In summary, we felt that the manuscript needs to make a more substantiative theoretical advances and/or provide more robust actionable insights into drug resistance evolution to justify publication in eLife. We would be happy to reconsider this manuscript should you be able to make substantive progress towards these issues.
Reviewer #1:
Ardell and Kryazhimskiy use bacterial KO data in multiple conditions to study the structure of pleiotropy, that is the degree to which a genetic perturbation affects multiple phenotypes, and present a theoretical framework to predict and assess fitness trajectories observed in environments other than the one selection is operating in. The work is thoroughly done and has potentially interesting implications for sequential drug therapy.
The central object of their framework is the joint distribution of fitness effects of mutations in multiple environments where the distribution is over all mutations in the genome. The dynamics in the space of fitness in multiple environments is then modeled as a random walk (described by a diffusion equation) assuming that mutations sweep separated in time (SSWM). The model and the calculations necessary to arrive at the predictions are simple and transparent. The results quantitatively predict simulation results with the range of validity of SSWM. Outside this range, the model predicts the qualitative behavior, but is quantitatively wrong.
1) My main disappointment with the paper is the inability to quantitatively describe the dynamics outside the SSWM regime. I would expect that the effects of competing mutations or weak selection could be accounted for at least perturbatively. Alternative, one could determine the distribution of the effects of fixed mutations in the "home" environment in simulations and use this distribution to predict the dynamics in other environments.
2) My other substantial concern is the question whether anything can be learned about drug resistance evolution or collateral sensitivity/resistance from KO experiments. While some drug resistance evolution involves lossoffunction mutations (e.g. porin losses), it often proceeds via point mutations, upregulation, or horizontal acquisition. Furthermore, the statistical treatment here requires many mutations to sample the joint effect distribution to give reliable answers. In clinical resistance evolution, the number of mutations observed is often quite small and their effect distributions are wide. The practical relevance of this is therefore far from clear.
3) While the similarity of this work to similar questions in quantitative genetics is discussed in the introduction, I would like to see an extended discussion whether some limits of the model at hand can be described by the quantitative genetics approach.
Reviewer #2:
The authors present a theoretical framework for analysing pleiotropic effects in populations evolving in different environments based on the concept of a joint distribution of fitness effects (JDFE). Simple correlation measures are derived from the JDFE that allow one to predict the evolutionary outcome in the nonhome environment. Analytic theory is derived in the SSWM regime and complemented by simulations covering the regime of large mutation supply. A proofofconcept application to collateral antibiotic resistance and sensitivity in bacteria based on a published data set for knockout strains is presented. Overall, this is an important, systematic contribution to very timely subject that is well suited for publication in eLife.
1. I do not quite share the authors' surprise at the outcomes shown in Figure 1. In fact, there is a simple heuristic that allows one to predict the direction of the fitness change in the nonhome environment in all cases: Simply look at the ycoordinate of the tail of the JDFE corresponding to the largest beneficial effects along the xaxis.
2. Along the three rows of panels in Figure 2, there appears to be a systematic but in two cases nonmonotonic variation of the slope with the mutation supply NU_b. Do the authors have a (tentative) explanation for this behavior?
Reviewer #3:
The goal of this manuscriptto develop predictive tools for inferring fitness trajectories in new environmentsis an important goal and I appreciate the synthesis of theoretical modeling with parameter estimation from empirical mutation studies.
Reading through the manuscript, however, I found myself repeatedly wondering whether the stated application of the methods developed here doesn't constitute something of a tautology. This could be a misreading on my end, but I'll explain: the authors state that they have the central goal of predicting whether a population adapting to one environment will lose fitness in another "nonhome" environment. Yet the parameter estimation they develop and propose for estimating fitness trajectories requires fitness measurements in both the home and nonhome environments. If one already has fitness measurements for both home and nonhome, how much more information is added by estimating the JDFE? I understand that the authors are estimating the fitness trajectories over time, with the incorporation of population genetic parameters, but again, I was unsure of how much information was added with the JDFE particularly given large discrepancies in the WrightFisher models and the decreasing predictive capacity with time. The bottom row of Figure 1 provided perhaps the most convincing evidence of the usefulness of the JDFE, but the unintuitive result was not adequately explored nor explained (see comment below). Also, perhaps an exploration of how the predictions could be extended to unmeasured environments is possible (as in Kinsler et al. 2020)?
Further specific conceptual comments and suggestions:
1) The authors demonstrate in Figure 1 that JDFEs even with similar shapes produce markedly different fitness trajectories. They argue that the correlation coefficient of the JDFE is not a reliable predictor of fitness trajectories in the home environment. I was struck by this counterintuitive result, and found myself searching for further explanation. Are the authors arguing that the practice of simply looking at the correlation coefficient in tradeoff studies in general is insufficient for predicting the fates of pleiotropic mutations? Either way, it would be helpful to the reader to elaborate on why and under which conditions the discrepancy with the correlation coefficient and fitness trajectories arises.
2) The modeling results throughout the manuscript reveal poor predictive capabilities in WrightFisher simulations. For example, the results in figure 2 show substantial discrepancy between the theoretical predictions and the results of the WrightFisher simulations. The authors address this only briefly stating that outside of the strong selection, weak mutation model (SSWM) the pleiotropy statistics are only "statistical predictors". But the discrepancy was systematic and wide, suggesting rather little insight from the pleiotropy statistics in sequential adaptation scenarios. I could not find discussion of this discrepancy between the SSWM and WrightFisher modeling predictions.
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "The Population Genetics of Collateral Resistance and Sensitivity" for further consideration by eLife. Your revised article has been evaluated by Meredith Schuman (Senior Editor) and a Reviewing Editor.
The manuscript has been improved considerably. The reviewers particularly appreciated the more general theoretical results. However, we would like to see a more thorough discussion of previous literature. In particular, the study by Nichol et al. is important prior work that needs to be discussed in greater depth.
In addition to a more extensive discussion of the literature, we would also like to see, if at all possible, some level of empirical validation of the results beyond the KO data presented so far. The data by Stiffler et al. and Nichol et al. characterize point mutations. Using these data beyond what is currently shown in Figure S1 could be very valuable.
Reviewer #2:
The manuscript has been largely restructured and rewritten, and has improved considerably. The issues that I raised in my previous report have been adequately addressed. As requested in the previous round of reviews, the theory has been generalized beyond the SSWM regime. Moreover, a second data set on resistance mutations in βlactamase has been added, and the discussion of the practical applicability of the approach has been extended significantly. As far as I am concerned, the paper can be accepted in its present form.
Reviewer #3:
The authors thoroughly responded to the reviewers' comments and I found the resubmission to be both clearer and to demonstrate greater prediction accuracy in the Wright Fisher simulations. The addition of the section on estimating JDFE parameters from experimental data was a positive addition to the manuscript in that it provides a bridge for experimentalists to implement the methods developed by the authors.
That being said, as an experimentalist who could potentially implement the proposed modeling in my own work for predicting tradeoffs, I am not yet convinced of the significant advance of the proposed modeling framework for making predictions. Specifically, I found the following two points to present the most significant drawbacks to the manuscript at present:
i) I found the manuscript to lack sufficient discussion of what has been shown before in the field of modeling collateral resistances and how the present manuscript presents a clear advance in light of this work. To the first point, a brief perusal of recent literature on collateral resistance brought me to Nichol et al. 2019 Nature Comm. Ardell et al. reference the Nichol manuscript on line 37 when stating that previous work observes wide variation in collateral outcomes. But Nichol et al. did more than demonstrate variation in collateral outcomes, and instead conducted 60 parallel experimental evolution assays in one antibiotic, measured the probability of collateral resistance/susceptibility and then modeled through SSWM simulations the predicted collateral resistance outcomes for dozens of drug pairs. The present manuscript should explain how its methods/goals/results differ from those of Nichol et al.
My second point is (ii) the manuscript would be significantly strengthened if it could provide proofofconcept validations beyond the KO work and the βlactamase work. If I understand correctly, the authors perform the drugranking experiments with simulated data. I am surprised that the authors cannot find a dataset in which to validate any part of the drugranking predictions. This type of validation would be helpful in convincing the reader of the strength of the proposed methods. As a relevant aside, Beyond Figure S1 I couldn't find where the Βlactamase data was used and the basic conclusion stated in the text for S1 regarding variable resistance pleiotropy is already wellestablished in previous work.
https://doi.org/10.7554/eLife.73250.sa1Author response
[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]
The reviewers agreed that pleiotropy of mutations and the resulting adaptive trajectories across different environment are important topics that are both of theoretical and applied interest. Your theoretical framework predicts fitness trajectories observed in environments other than the one selection is operating in (home env). These trajectories in nonhome environments are calculated via integrals over the joint fitness effect distribution weighted by the fixation probability in the home environment. However, your framework assumes strong selection and weak mutation (SSWM) and deviations from this assumption seem to have strong effects. We think that these effects need to be at least partially understood. Furthermore, application to the KO library is a useful proofofconcept, but the practical relevance of these patterns for understanding collateral sensitivity/resistance is far from obvious.
In summary, we felt that the manuscript needs to make a more substantiative theoretical advances and/or provide more robust actionable insights into drug resistance evolution to justify publication in eLife. We would be happy to reconsider this manuscript should you be able to make substantive progress towards these issues.
Reviewer #1:
Ardell and Kryazhimskiy use bacterial KO data in multiple conditions to study the structure of pleiotropy, that is the degree to which a genetic perturbation affects multiple phenotypes, and present a theoretical framework to predict and assess fitness trajectories observed in environments other than the one selection is operating in. The work is thoroughly done and has potentially interesting implications for sequential drug therapy.
The central object of their framework is the joint distribution of fitness effects of mutations in multiple environments where the distribution is over all mutations in the genome. The dynamics in the space of fitness in multiple environments is then modeled as a random walk (described by a diffusion equation) assuming that mutations sweep separated in time (SSWM). The model and the calculations necessary to arrive at the predictions are simple and transparent. The results quantitatively predict simulation results with the range of validity of SSWM. Outside this range, the model predicts the qualitative behavior, but is quantitatively wrong.
1) My main disappointment with the paper is the inability to quantitatively describe the dynamics outside the SSWM regime. I would expect that the effects of competing mutations or weak selection could be accounted for at least perturbatively. Alternative, one could determine the distribution of the effects of fixed mutations in the "home" environment in simulations and use this distribution to predict the dynamics in other environments.
We agree with both Dr. Neher and reviewer 3 in their assessment of the quantitative disagreement between theory and simulations. Although the development of a full theory of evolution on a JDFE in the concurrent mutations regime goes beyond the scope of this paper, we were able to extend our main results to this regime by using the adjusted fixation probability derived by Good et al. (PNAS 2012). This simple adjustment now produces a good quantitative agreement between theory and simulations.
While developing a quantitative theory of pleiotropy is important, it is also important to keep in mind that quantitative predictions would likely be impossible in many practical situations (such as when decisions need to be made with which drugs to treat a patient) because the population genetic parameters of the treated population will not be known. In such practical situations, a ranking of drug pairs by their risk of collateral resistance evolution may be more useful than a quantitative theory, as long as the ranking is robust with respect to the population genetic parameters. We discuss this problem and our solution to it explicitly in the new section “Robust ranking of drug pairs”.
2) My other substantial concern is the question whether anything can be learned about drug resistance evolution or collateral sensitivity/resistance from KO experiments. While some drug resistance evolution involves lossoffunction mutations (e.g. porin losses), it often proceeds via point mutations, upregulation, or horizontal acquisition. Furthermore, the statistical treatment here requires many mutations to sample the joint effect distribution to give reliable answers. In clinical resistance evolution, the number of mutations observed is often quite small and their effect distributions are wide. The practical relevance of this is therefore far from clear.
Dr. Neher raises two important concerns about the practical relevance of our work. (1) Relevance of the knockout mutations for predicting resistance evolution. (2) Potential difficulties in estimating JDFEs from limited data for the JDFE framework to be useful. Overall, we have reorganized and augmented the manuscript to make it more practically relevant, specifically refocusing it on the evolution of collateral resistance and sensitivity. To address the two specific issues raised by the reviewer, we did the following.
1) First, we note that Chevereau et al. (PLoS Biol, 2015) show that the answer to “the question whether anything can be learned about drug resistance evolution” from KO experiments is “yes”, at least in some cases (see Figure 4 in their paper). This important result notwithstanding, we agree with the reviewer that it is in general unclear what KO mutations alone tell us about full JDFEs. Our presentation in the original manuscript suggested a more direct applicability of koJDFEs than intended. Instead, our central thesis is that the JDFE framework is useful for thinking about collateral resistance/sensitivity evolution, not that actual JDFEs can be estimated from knockout mutations. To make this point clear, we made the following changes.
a) We removed any discussion of the pleiotropy parameters of koJDFE. In fact, we no longer use the term koJDFE.
a) We moved the presentation of the KO data to the beginning of the Results section. We now use these data to demonstrate that resistance mutations are highly diverse in terms of their collateral effects. This fact may seem obvious in retrospect, but, to our knowledge, it has never been made explicitly in the antibiotic resistance literature. We use this argument to justify the JDFE framework for modeling collateral resistance/sensitivity.
c) To further bolster this argument, we now use a second data set, the mutational scanning data in the TEM1 βlactamase gene, obtained by Stiffler et al. (Cell 2015). We show that even point mutations in a single gene exhibit a great degree of variability in their collateral effects.
2) First, we want to clarify that for the JDFE framework to be useful the JDFE need not consist of many mutations. For example, a JDFE can be discrete. But as long as there are multiple classes of mutations with different pleiotropic effects, the pleiotropic outcomes of evolution will be uncertain. Understanding this uncertainty and taking it into account for designing drug treatments requires knowing the probabilities with which different classes of mutations arise, which is the JDFE.
Having said that, we agree with the reviewer that, if JDFEs cannot be reliably estimated, the utility of this framework would be greatly diminished. To address this issue, we added a new section “Measuring JDFEs” where we show how JDFEs can be estimated for practical purposes of drug treatment choice from relatively little data, at least if the JDFEs have relatively simple shapes.
3) While the similarity of this work to similar questions in quantitative genetics is discussed in the introduction, I would like to see an extended discussion whether some limits of the model at hand can be described by the quantitative genetics approach.
We have extensively revised the manuscript, including the Discussion, but we are not quite sure what limits Dr. Neher has in mind here. If he still thinks it is necessary to discuss some additional points, we would appreciate any specific suggestions.
Reviewer #2:
The authors present a theoretical framework for analysing pleiotropic effects in populations evolving in different environments based on the concept of a joint distribution of fitness effects (JDFE). Simple correlation measures are derived from the JDFE that allow one to predict the evolutionary outcome in the nonhome environment. Analytic theory is derived in the SSWM regime and complemented by simulations covering the regime of large mutation supply. A proofofconcept application to collateral antibiotic resistance and sensitivity in bacteria based on a published data set for knockout strains is presented. Overall, this is an important, systematic contribution to very timely subject that is well suited for publication in eLife.
1. I do not quite share the authors' surprise at the outcomes shown in Figure 1. In fact, there is a simple heuristic that allows one to predict the direction of the
fitness change in the nonhome environment in all cases: Simply look at the ycoordinate of the tail of the JDFE corresponding to the largest beneficial effects along the xaxis.
While Dr. Krug did not find the observations in Figure 2 (former Figure 1) surprising, Reviewer 3 did. The degree of surprise naturally varies, but it is quite certain that the majority of eLife readers would not have as well developed intuition about evolution as Dr. Krug. We therefore believe that this figure will help many readers to better appreciate the complexity of the problem at hand. Specifically, one important point that this figure makes is that the correlation between mutational effects—a commonly used metric of tradeoffs—does not predict pleiotropic outcomes of evolution.
As a side note, if we understand the heuristic suggested by Dr. Krug correctly, it may have been an artifact of plotting. After this comment, we noticed that the tip of the outermost contour line in the original version of the figure did indeed predict the direction of collateral evolution. But of course this line is arbitrary, and how to identify “the largest beneficial effects along the xaxis” is unclear. We have now changed the spacing between the contour lines so that there are no more spurious relationships (that we can see) between these lines and the observed collateral outcomes of evolution.
2. Along the three rows of panels in Figure 2, there appears to be a systematic but in two cases nonmonotonic variation of the slope with the mutation supply NU_b. Do the authors have a (tentative) explanation for this behavior?
Thank you for pointing out this issue. An investigating of this issue led us to discover an error in Equation (3). In this equation, the full JDFE had to be replaced with the JDFE of mutations beneficial in the home environment, which we now did. After correcting this error, we find that the SSWM predictions always overestimate the actual rate of fitness change in the home and in the nonhome environments, as they should (see Figure S3). As a consequence, all data points in the first quarter of the plots in Figure S3 (these data were presented in Figure 2 in the previous version of the manuscript) are now below the diagonal and all data points in the third quarter are above the diagonal, and the change in slope is monotonic across panels, as expected.
We observe some nonmonotonicity across panels in the current Figure 3, that is, (^{2}*) underestimates the observed rate of nonhome fitness gains at low and the same issue as above. But there is also no expectation that2* should slightly overestimates it at higher values. This nonmonotonicity is not caused by consistently over or underestimate the observed rates. At this point, we believe that the observed behavior is a property of the approximation derived by Good et al. (PNAS 2015).
Reviewer #3:
The goal of this manuscriptto develop predictive tools for inferring fitness trajectories in new environmentsis an important goal and I appreciate the synthesis of theoretical modeling with parameter estimation from empirical mutation studies.
Reading through the manuscript, however, I found myself repeatedly wondering whether the stated application of the methods developed here doesn't constitute something of a tautology. This could be a misreading on my end, but I'll explain: the authors state that they have the central goal of predicting whether a population adapting to one environment will lose fitness in another "nonhome" environment. Yet the parameter estimation they develop and propose for estimating fitness trajectories requires fitness measurements in both the home and nonhome environments. If one already has fitness measurements for both home and nonhome, how much more information is added by estimating the JDFE? I understand that the authors are estimating the fitness trajectories over time, with the incorporation of population genetic parameters, but again, I was unsure of how much information was added with the JDFE particularly given large discrepancies in the WrightFisher models and the decreasing predictive capacity with time. The bottom row of Figure 1 provided perhaps the most convincing evidence of the usefulness of the JDFE, but the unintuitive result was not adequately explored nor explained (see comment below). Also, perhaps an exploration of how the predictions could be extended to unmeasured environments is possible (as in Kinsler et al. 2020)?
There are several interrelated issues in this comment: (1) Large discrepancies between theory and simulations; (2) Is the theory tautological? (3) If one already has fitness measurements for both home and nonhome, how much more information is added by estimating the JDFE? (4) Unintuitive result in Figure 1 is not explained (5) Can predictions be made about unmeasured environments? We now address these issues one by one.
1) This concern was also raised by Reviewer 1, and we address it in more detail in response to him. The main point is that we have extended our theory and now have good quantitative correspondence between theory and simulations.
2) Our theory is not tautological, which can be asserted on several accounts. Consider Figure 2 (former Figure 1), which is the idealized case when the full JDFE is known. The reviewer admits that even in this idealized case it is not obvious how to make a basic qualitative prediction: will an average population evolve collateral resistance or sensitivity? Prediction is challenging because different mutations contribute differently to adaptation, and a population genetic theory is necessary to properly weigh them. Our theory tells us what summary statistics of the JDFE we need to use to correctly predict the rate of collateral evolution (as well as the uncertainty around this expectation).
While predicting the average qualitative collateral outcome is already challenging, we argue that what really matters in practice is to rank drug pairs by their risk of collateral resistance. Assessing this risk seems impossible without a population genetic model. We have now added a new section where we show how our theory helps us rank drug pairs (section “Robust ranking of drug pairs”).
3) The reviewer alludes here to a more realistic scenario when the full JDFE is not known, but instead “one already has fitness measurements for both home and nonhome” from a finite sample of mutants. The problems are the same as in the idealized case with a fully known JDFE: How do we combine these fitness measurements into an effective estimator of the rate of collateral evolution? If we have fitness measurements for these mutants in the presence of multiple drugs, how do we use these measurements to rank drug pairs according to the risk of collateral resistance? Again, our theory provides answers to these questions. Furthermore, we have now added a new section “Measuring JDFEs” where we discuss the issues related to measurement of JDFE.
4) Thank you for pointing out this omission. Our theory does indeed help us understand the observations made in Figure 1. We have added the relevant discussion (see p. 9, l. 245).
5) The question of how collateral evolution would proceed in unmeasured environments is fascinating and important, but unfortunately it cannot be answered within the current JDFE framework. To answer it, we either need to have some heuristic knowledge of how the effects of mutations vary across environments or we need to know something about the mechanistic basis of resistance and collateral resistance/sensitivity. Fortunately, a lot is known about these mechanisms, at least for some drug pairs, and we are now working in this direction.
Further specific conceptual comments and suggestions:
1) The authors demonstrate in Figure 1 that JDFEs even with similar shapes produce markedly different fitness trajectories. They argue that the correlation coefficient of the JDFE is not a reliable predictor of fitness trajectories in the home environment. I was struck by this counterintuitive result, and found myself searching for further explanation. Are the authors arguing that the practice of simply looking at the correlation coefficient in tradeoff studies in general is insufficient for predicting the fates of pleiotropic mutations? Either way, it would be helpful to the reader to elaborate on why and under which conditions the discrepancy with the correlation coefficient and fitness trajectories arises.
Correct: correlation between the effects of mutations on traits is not sufficient to predict pleiotropic outcomes of evolution. As mentioned above, our theory helps us understand why this is so (see p. 9, l. 245).
2) The modeling results throughout the manuscript reveal poor predictive capabilities in WrightFisher simulations. For example, the results in figure 2 show substantial discrepancy between the theoretical predictions and the results of the WrightFisher simulations. The authors address this only briefly stating that outside of the strong selection, weak mutation model (SSWM) the pleiotropy statistics are only "statistical predictors". But the discrepancy was systematic and wide, suggesting rather little insight from the pleiotropy statistics in sequential adaptation scenarios. I could not find discussion of this discrepancy between the SSWM and WrightFisher modeling predictions.
As mentioned above, this is now resolved.
[Editors’ note: what follows is the authors’ response to the second round of review.]
Essential revisions:
The manuscript has been improved considerably. The reviewers particularly appreciated the more general theoretical results. However, we would like to see a more thorough discussion of previous literature. In particular, the study by Nichol et al. is important prior work that needs to be discussed in greater depth.
We have now expanded our discussion of the prior literature in the Introduction (LL. 25–38). We also added a more extensive discussion of the Nichol et al. paper specifically (LL. 73–79).
In addition to a more extensive discussion of the literature, we would also like to see, if at all possible, some level of empirical validation of the results beyond the KO data presented so far. The data by Stiffler et al. and Nichol et al. characterize point mutations. Using these data beyond what is currently shown in Figure S1 could be very valuable.
We agree that we have not exploited the data by Stiffler et al. to its fullest extent in the previous version of the manuscript. We now utilize it further to show that our theory predicts the correct in silico ranking of drug pairs even for these more realistic nonGaussian JDFEs. This new result, shown in the new Figure 4 —figure supplement 1, bolsters the optimism that our theory will work beyond in silico experiments. However, we realize that this is not an empirical validation and agree that such validation would be highly desirable.
An empirical validation of our theory would require two pieces of data: (1) an estimate of a JDFE for at least one pair of drugs and (2) an evolution experiment in the presence of one of these drugs as well as a fitness measurement after evolution in the presence of the other drug. Notably, the evolution experiment must be carried out with a sufficient number of replicate populations so that we can estimate the probability to evolve collateral resistance, or at least the variance in the collateral outcome. We therefore systematically searched the existing literature for datasets that satisfy these two criteria. We describe this metaanalysis below. The main conclusions of this metaanalysis are the following.
First, there is no single study that provides both pieces of data. Second, there are a handful of studies that can be paired so that one provides a JDFE estimate and the other provides experimental evolution data. Stiffler et al. and Nichol et al. form one such pair. However, none of the pairs are sufficiently aligned in terms of their experimental design that would allow us to test our theory without serious confounding factors. The biggest problem is that fullgenome JDFEs have not yet been reported in the literature. JDFEs can be estimated for single genes from deep mutational scanning data. However, almost all evolution experiments are carried out with whole organisms and resistance mutations arise across several genes.
Thus, it is unfortunately not possible to empirically validate our theory using existing data. We plan to carry out the necessary experiments in our laboratory, but this work goes far beyond the scope of this paper.
Reviewer #3:
The authors thoroughly responded to the reviewers’ comments and I found the resubmission to be both clearer and to demonstrate greater prediction accuracy in the Wright Fisher simulations. The addition of the section on estimating JDFE parameters from experimental data was a positive addition to the manuscript in that it provides a bridge for experimentalists to implement the methods developed by the authors.
That being said, as an experimentalist who could potentially implement the proposed modeling in my own work for predicting tradeoffs, I am not yet convinced of the significant advance of the proposed modellng framework for making predictions. Specifically, I found the following two points to present the most significant drawbacks to the manuscript at present:
i) I found the manuscript to lack sufficient discussion of what has been shown before in the field of modeling collateral resistances and how the present manuscript presents a clear advance in light of this work. To the first point, a brief perusal of recent literature on collateral resistance brought me to Nichol et al. 2019 Nature Comm. Ardell et al. reference the Nichol manuscript on line 37 when stating that previous work observes wide variation in collateral outcomes. But Nichol et al. did more than demonstrate variation in collateral outcomes, and instead conducted 60 parallel experimental evolution assays in one antibiotic, measured the probability of collateral resistance/susceptibility and then modeled through SSWM simulations the predicted collateral resistance outcomes for dozens of drug pairs. The present manuscript should explain how its methods/goals/results differ from those of Nichol et al.
We have expanded the introduction to address this concern, as described in the response to the editor. We agree that the paper by Nichol et al. is important and very relevant to the present study, as we now discuss in LL. 72–79 and mention in several other places in the text. In particular, Nichol et al. are motivated by the same problem as our work. However, our approach to solving it is conceptually different, as we discuss in LL. 79–98. It is also worth noting that Nichol et al. did not quantitatively test their model against their own experimental data. In fact, their experimental design does not allow them to do so, as it suffers from the same problem as we face in testing our theory, i.e., the fact that their theory is based on fitness landscapes consisting of four mutations in a single gene, but resistance mutations in their experiment often arise elsewhere in the genome. In the absence of direct empirical tests, we argue on theoretical grounds that the JDFEbased approach is superior to the approach based on combinatorically complete landscapes (LL. 79–98).
My second point is (ii) the manuscript would be significantly strengthened if it could provide proofofconcept validations beyond the KO work and the βlactamase work. If I understand correctly, the authors perform the drugranking experiments with simulated data. I am surprised that the authors cannot find a dataset in which to validate any part of the drugranking predictions. This type of validation would be helpful in convincing the reader of the strength of the proposed methods. As a relevant aside, Beyond Figure S1 I couldn't find where the Βlactamase data was used and the basic conclusion stated in the text for S1 regarding variable resistance pleiotropy is already wellestablished in previous work.
Please see our detailed response to the editor’s comment above, as well as the results of our systematic literature search described below. In short, we agree that empirical validation is important. However, it is unfortunately impossible to validate our theory with existing data. This may be surprising given the vastness of the literature on antibiotic resistance, but this lack of relevant data can perhaps be explained in hindsight. First, the technical capabilities for estimating genomewide JDFE (e.g., using the barcodelineage tracking method) have appeared only recently and are not yet widely adopted. And second, the importance of stochasticity for collateral outcomes is not yet widely appreciated in the field. We hope that the reviewer appreciates the possibility that publishing this theoretical study may in fact stimulate empiricists to carry out the relevant measurements.
https://doi.org/10.7554/eLife.73250.sa2Article and author information
Author details
Funding
Burroughs Wellcome Fund (1010719.01)
 Sergey Kryazhimskiy
Alfred P. Sloan Foundation (FG20179227)
 Sergey Kryazhimskiy
Hellman Foundation
 Sergey Kryazhimskiy
National Institutes of Health (1R01GM137112)
 Sergey Kryazhimskiy
National Institutes of Health (1T32GM13335101)
 Sarah M Ardell
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Shea Summers and Flora Tang for help and input at the initial stages of the project and Kryazhimskiy and Meyer labs for feedback. We thank Tobias Bollenbach and Guillaume Chevereau for providing wildtype growth rate measurement data. This work was supported by the BWF Career Award at Scientific Interface (Grant 1010719.01), the Alfred P Sloan Foundation (Grant FG2017–9227), the Hellman Foundation and NIH (Grants 1R01GM137112 and 1T32GM13335101).
Senior Editor
 Meredith C Schuman, University of Zurich, Switzerland
Reviewing Editor
 Richard A Neher, University of Basel, Switzerland
Reviewer
 Joachim Krug, University of Cologne, Germany
Publication history
 Received: August 22, 2021
 Accepted: December 6, 2021
 Accepted Manuscript published: December 10, 2021 (version 1)
 Version of Record published: January 18, 2022 (version 2)
Copyright
© 2021, Ardell and Kryazhimskiy
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,652
 Page views

 212
 Downloads

 4
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Ecology
 Evolutionary Biology
Evolutionary theory suggests that individuals should express costly traits at a magnitude that optimizes the trait bearer’s costbenefit difference. Trait expression varies across a species because costs and benefits vary among individuals. For example, if large individuals pay lower costs than small individuals, then larger individuals should reach optimal costbenefit differences at greater trait magnitudes. Using the cavitationshooting weapons found in the big claws of male and female snapping shrimp, we test whether size and sexdependent expenditures explain scaling and sex differences in weapon size. We found that males and females from three snapping shrimp species (Alpheus heterochaelis, Alpheus angulosus, and Alpheus estuariensis) show patterns consistent with tradeoffs between weapon and abdomen size. For male A. heterochaelis, the species for which we had the greatest statistical power, smaller individuals showed steeper tradeoffs. Our extensive dataset in A. heterochaelis also included data about pairing, breeding season, and egg clutch size. Therefore, we could test for reproductive tradeoffs and benefits in this species. Female A. heterochaelis exhibited tradeoffs between weapon size and egg count, average egg volume, and total egg mass volume. For average egg volume, smaller females exhibited steeper tradeoffs. Furthermore, in males but not females, large weapons were positively correlated with the probability of being paired and the relative size of their pair mates. In conclusion, we identified sizedependent tradeoffs that could underlie reliable scaling of costly traits. Furthermore, weapons are especially beneficial to males and burdensome to females, which could explain why males have larger weapons than females.

 Evolutionary Biology
The extinct Steller's sea cow (Hydrodamalis gigas; †1768) was a whalesized marine mammal that manifested profound morphological specializations to exploit the harsh coastal climate of the North Pacific. Yet despite firsthand accounts of their biology, little is known regarding the physiological adjustments underlying their evolution to this environment. Here, the adultexpressed hemoglobin (Hb; a_{2}β/δ_{2}) of this sirenian is shown to harbor a fixed amino acid replacement at an otherwise invariant position (β/δ82Lys→Asn) that alters multiple aspects of Hb function. First, our functional characterization of recombinant sirenian Hb proteins demonstrate that the HbO_{2} affinity of this subArctic species was less affected by temperature than those of living (sub)tropical sea cows. This phenotype presumably safeguarded O_{2} delivery to cool peripheral tissues and largely arises from a reduced intrinsic temperature sensitivity of the H. gigas protein. Additional experiments on H. gigas β/δ82Asn→Lys mutant Hb further reveal this exchange renders Steller's sea cow Hb unresponsive to the potent intraerythrocytic allosteric effector 2,3diphosphoglycerate, a radical modification that is the first documented example of this phenotype among mammals. Notably, β/δ82Lys→Asn moreover underlies the secondary evolution of a reduced bloodO_{2} affinity phenotype that would have promoted heightened tissue and maternal/fetal O_{2} delivery. This conclusion is bolstered by analyses of two Steller's sea cow prenatal Hb proteins (Hb Gower I; z_{2}e_{2} and HbF; a_{2}g_{2}) that suggest an exclusive embryonic stage expression pattern, and reveal uncommon replacements in H. gigas HbF (g38Thr→Ile and g101Glu→Asp) that increased HbO_{2} affinity relative to dugong HbF. Finally, the β/δ82Lys→Asn replacement of the adult/fetal protein is shown to increase protein solubility, which may have elevated red blood cell Hb content within both the adult and fetal circulations and contributed to meeting the elevated metabolic (thermoregulatory) requirements and fetal growth rates associated with this species cold adaptation.