Survival of the simplest in microbial evolution

The evolution of microbial and viral organisms often generates clonal interference, a mode of competition between genetic clades within a population. Here we show how interference impacts systems biology by constraining genetic and phenotypic complexity. Our analysis uses biophysically grounded evolutionary models for molecular phenotypes, such as fold stability and enzymatic activity of genes. We find a generic mode of phenotypic interference that couples the function of individual genes and the population’s global evolutionary dynamics. Biological implications of phenotypic interference include rapid collateral system degradation in adaptation experiments and long-term selection against genome complexity: each additional gene carries a cost proportional to the total number of genes. Recombination above a threshold rate can eliminate this cost, which establishes a universal, biophysically grounded scenario for the evolution of sex. In a broader context, our analysis suggests that the systems biology of microbes is strongly intertwined with their mode of evolution.

In this paper, the authors build a theory of "phenotypic interference". In short, this is a theory for the evolution of clonal populations in which (even housekeeping) selection for many additive traits (e.g. thermodynamic stability of proteins) results in a load that scales quadratically with the number of traits due to clonal interference. This mode of evolution has been previously unobserved, because (i) of neglecting biophysical realism due to nonlinear behavior of individual traits; (ii) of nontrivial coupling between the local (single trait) and global (fitness wave) levels of evolution.
I believe that this paper makes a fundamental contribution to the field and recommend that it is published.
My major recommendation is that the clarity of initial exposition, as I describe below, could be increased to make the paper parseable more widely than to the specialized audience dealing with clonal interference. Several concrete suggestions: 1) For someone not familiar with the concept of "fitness wave" in detail it is hard to parse the intro para "Over a wide range of model parameters..." What I would find more important is to highlight first the feedback between local and global selection and explain what this is in words. Here, or perhaps later (in "theory of phenotypic interference" initial para) it would really help to highlight how such feedback emerges: through the fact that the \tilde{sigma} (related to eff. pop size) is set by the evolutionary dynamics itself. This is interesting to wider field as by default eff. pop size is often just assumed to be some constant, which would be wrong in this context. This is, if I understand properly, the gist of point (ii) above; 2) I would also highlight more point (i) above, namely the necessity of treating the essential biophysical nonlinearity (figure 1); right now, this is said in caption to Fig 1 (where you compare biophys model in red and discrete fitness model in blue), but I suggest you highlight this when you call out Figure 1B.
3) As above, for "Theory of phenotypic interference" first para, prepare the reader how the local and global fitness dynamics will be coupled, i.e., by mentioning specifically that the local part will assume some \tilde{sigma}, which will be then set by the global fitness wave part. This I understood only after reading through the whole manuscript first and being confused. 4) Minor, in the subsection "Evolution of a quantitative trait under interference selection" --that subsection is just a general section on the evolution of a quantitative trait, as the interference doesn't come in yet, correct? Perhaps you can amend the title. Also: can you make clear what parts of this subsection are novel (e.g., relations in Eq (1) and (2) appear straightforward, has this been done before)? 5) I had the most trouble parsing the crucial section "Housekeeping evolution of multiple traits". First, it would be good to be more clear what are the necessary conditions for the fitness wave regime to hold in this case. Make it more explicit that the fitness wave theory for your case contributes the crucial equation that relates the fitness variance and the coalescence rate, and perhaps give a few sentence intuition behind it. Moreover, currently you assert c0~100, which requires at least a few words about where the magic constant emerges from. Second, you refer to c as the "slowly varying parameter". In what sense is this slow, can you be more precise? Third, can you comment briefly about how Eq 5 is consistent with the sigma^2/\tilde{sigma}^2 = c0 log(N\sigma) from the fitness wave theory? Fourth, you say "As sown in Methods, this parameter [c] has a simple interpretation:..." I think a few sentences more here in the main text would help, c is the crucial quantity. Fifth, you say "The calibration between theory and data involves the 1. The paper's main novelty seems to be studying quantitative traits in the regime of clonal interference, thereby introducing the concept of "phenotypic interference." It was not clear to me what is gained from this approach beyond the existing body of work on clonal interference. The authors imply that "phenotypic interference" is a fundamentally new concept, but to me it just seems like clonal interference applied to quantitative traits. I furthermore found the emphasis on biophysical traits like protein stability misleading, because it didn't seem important to the authors' results: they didn't apply their results to any protein evolution data, nor did they make the argument that destabilizing mutations to proteins should be of dominant importance. Is the point just to use any nonlinear, diminishing-returns fitness function? If so, I think it would be clearer to consider abstract traits with such properties, making the argument that these fitness functions should be common among molecular traits (with protein stability being just one example).
2. The main result of this paper seems to be predicting how clonal interference of quantitative traits constrains genome size. This is of course an important problem, but there is a long history of studies addressing it, and I did not feel that literature was adequately discussed and compared to the present result. This made it hard for me to tell how significantly different the authors' results are. In particular, constraints on genome size resulting from accumulation of mutations destabilizing protein folding was discussed by Zeldovich et al. (PNAS 2007). As far as I can tell, this paper is cited in the Methods, but not discussed in the context of the results.
3. Overall the paper is quite technical and difficult to read. In particular, the main text contains quite a lot of notation (some of which I think is nonstandard and possibly unnecessary ---see specific comments below). The authors should consider significant simplications to the presentation of their results if they hope for them to be accessible to an audience beyond theoretical population geneticists.
1. It would be helpful for the authors to include line numbers in their manuscript ---this would make it easier to refer to specific items in the text. I will do my best below to indicate where in the text my comments apply.
2. In my opinion Fig. 1b comes too early in the text ---the introduction, at which point I was definitely not prepared to understand it. Shouldn't this plot be part of a main figure in the Results section, where it can be fully explained?
3. The colors in Fig. 2 are called "blue" and "orange," although for me they appeared to be more of a teal and red. I'm not sure if this is a difference in software/hardware or my own vision, but the authors might double-check these colors in case readers get confused.
4. The main text introduces a lot of new notation, some of which appears to be nonstandard and unnecessary. For example, they use $G$ for the stability of a protein rather than the more standard $\Delta G$, which is compounded by defining $\Delta_G$ as the variance. Also, why use introduce $\Gamma$ instead of $\bar{G}$, and $\tilde{\sigma}$ instead of $N_e^{-1}$? To me this made many of the results confusing, because these expressions fundamentally involve classical combinations of selection coefficients, mutation rates, and effective population sizes, but their notation obscures this intuition.
5. The authors make a lot of mathematical statements without adequate explanation. For example, the authors need to provide conceptual explanations for why $\Delta_G f'(\Gamma)$ is the rate of stability increase by selection, and where the mutation-selection balanace condition ($\Delta_G f'(\Gamma) = u\epsilon_G$) comes from. Full derivations in the Methods/SI were also not cited for these.
6. In the section called "Evolution of a quantitative trait under interference selection," where is the effect of linkage? I couldn't tell the difference between these results and what one would expect for an unlinked trait evolving independently. 7. "The strength of selection on genetic variants is not fixed a priori, but is an emergent property of the global evolutionary process. A faster pace of evolution, i.e., an increase in coalescence rate $\tilde{\sigma}$, reduces the efficacy of selection." This sounds like restating a very standard concept in population genetics, that absolute selection strength doesn't matter so much as its value relative to competing processes like genetic drift (coalescence) and mutation. If the authors are implying something beyond this, then they should clarify.
8. "The non-adaptive, stationary scenario of housekeeping evolution..." I think calling that scenario "non-adaptive" is misleading: it is a steady state, so the mean fitness doesn't change, but there is still selection. To me "non-adaptive evolution" is devoid of any selection, e.g., purely genetic drift, mutation, recombination. 9. Above Eq. 3, the parameter $c_0$ appears and is given a value but with no explanation. Is it related to $c$ in Eq. 5?
10. In general Eqs. 3, 4, 5 are presented without adequate explanation. I gather they define properties of the fitness wave, but the authors need to introduce fitness waves more systematically and explain how they are parameterized. Figure 3 purports to show the two scaling regimes derived in the text between independentlyevolving genes and the linkage-dominated case. However, the plots are not convincing to me that their analytical results are particularly accurate ---the data does not match either regime all that well, nor does it show a clear crossover. This is especially concerning given that the data is from simulations, so whatever approximations the authors use in their analytical results may be problematic.

11.
12. The text leading up to Eq. 6 should precisely define genetic load before introducing results on it.
13. At the bottom of page 6, the authors start mentioning the effect of all sorts of other scenarios on the genetic load result (Eq. 6), but this is way too terse to understand, and there was no reference to the SI which I later realized fully explains these cases.
14. I was surprised to see no mention of Neher and Shraiman (PNAS 2009) in the section on "The pathway to sexual evolution." I think that paper, which also shows the transition between linkagedominated and recombination-dominated cases, is an important comparison to this manuscript's results.
15. "The spectrum of site selection coefficients is not a fixed input, but a dynamical output of the evolutionary process. That is the main difference of our approach to previous population-genetic models of asexual evolution [24][25][26][27][28][29][30]." As I understand it, this statement does not seem to be supported by the references or the preceding text. The citations mentioned here also have nonlinear fitness functions, where the strength of selection on a trait depends on the current value of the trait. The main diffrence, I would say, is a lack of clonal interference.
16. The parameterization of the fitness function in Eqs. 10 and 11 does not seem to allow for essential proteins or those that are otherwise lethal when unfolded. As long as $f_0$ is finite, then even $G \to -\infty$ for a single protein will not result in zero total fitness (Eq. 11).
In this paper, the authors build a theory of "phenotypic interference". In short, this is a theory for the evolution of clonal populations in which (even housekeeping) selection for many additive traits (e.g. thermodynamic stability of proteins) results in a load that scales quadratically with the number of traits due to clonal interference. This mode of evolution has been previously unobserved, because (i) of neglecting biophysical realism due to nonlinear behavior of individual traits; (ii) of nontrivial coupling between the local (single trait) and global (fitness wave) levels of evolution.
I believe that this paper makes a fundamental contribution to the field and recommend that it is published. My major recommendation is that the clarity of initial exposition, as I describe below, could be increased to make the paper parseable more widely than to the specialized audience dealing with clonal interference. Several concrete suggestions: We thank the referee for his general positive feedback and for his concrete criticism that helped us to significantly improve the manuscript. We revised the manuscript extensively to make it more comprehensible for a general audience, e.g., with Box 1. Moreover, we went over all derivations again to streamline them. You find specific points in the detailed response below. 1) For someone not familiar with the concept of "fitness wave" in detail it is hard to parse the intro para "Over a wide range of model parameters..." What I would find more important is to highlight first the feedback between local and global selection and explain what this is in words. Here, or perhaps later (in "theory of phenotypic interference" initial para) it would really help to highlight how such feedback emerges: through the fact that theσ (related to eff. pop size) is set by the evolutionary dynamics itself. This is interesting to wider field as by default eff. pop size is often just assumed to be some constant, which would be wrong in this context. This is, if I understand properly, the gist of point (ii) above; Thanks for these suggestions. We have now completely reorganized the introduction to highlight the feedback between global and local selection as the central new concept of the paper. The relevant paragraph now reads Selection against genome complexity can be traced to a specific evolutionary mode of phenotypic interference, which is shown to occur in a parameter regime appropriate for typical microbial systems. The key characteristics of this mode are summarized in Box 1. First, local quantitative traits of a given gene are in an evolutionary equilibrium where the long-term average of the trait value, and its position of on the fitness landscape, are determined by the uphill force of selection and the downhill force of mutations (Box 1, top). Second, global genome evolution takes place in a so-called fitness wave; that is, genetic and phenotypic variants in multiple genes co-exist in a population and generate a broad distribution of fitness values [Desai and Fisher, 2007, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008 (Box 1, center). In a fitness wave, selection rather than genetic drift determines the coalescence rate of a population, which is the inverse of the effective population size , Rice et al., 2015. Third, the defining feature of phenotypic interference is that global and local evolution are coupled in a specific way: the global coalescence rate sets the average selection coefficient of local trait changes. This feedback and its consequences are not covered by the existing models of fitness waves [Desai and Fisher, 2007, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008 and of quantitative trait evolution [Berg et al., 2004, Chi and Liberles, 2016, Gerland and Hwa, 2002, Goldstein, 2011, Manhart and Morozov, 2015, Nourmohammad et al., 2013, Serohijos and Shakhnovich, 2014.
In addition, we have added a new Box 1 that summarizes the selective feedback and the key scales of selection and coalescence in a non-technical way. Box 1 also serves as a shortcut through the theory section for readers interested primarily in the biological results.
2) I would also highlight more point (i) above, namely the necessity of treating the essential biophysical nonlinearity (figure 1); right now, this is said in caption to Fig 1 (where you compare biophys model in red and discrete fitness model in blue), but I suggest you highlight this when you call out Figure 1B.
The importance of nonlinearities in gene fitness landscapes is now discussed in the introductory paragraph on biophysical models: Selection on a gene is described by a thermodynamic fitness landscape f (G), which is a sigmoid function with a high-fitness plateau corresponding to stable proteins and a low-fitness plateau corresponding to unfolded proteins (Box 1, top). We also discuss a stability-affinity protein model with a two-dimensional fitness landscape f (G, E); this model includes enzymatic or regulatory functions of genes, specifically the protein binding affinity E to a molecular target. From the perspective of molecular evolution, these landscapes provide a biophysically grounded model of local fitness epistasis, which couples all sequence sites contributing to a stability or affinity trait in the same gene. Importantly, genome-wide local epistasis in protein-coding sequence operates independently on any assumptions of fitness interactions across genes. Beyond proteins, local epistasis occurs ubiquitously in quantitative molecular traits associated with binding interactions. This form of epistasis will be an important building block of our dynamical analysis but is not covered by the standard theory of asexual evolution [Desai and Fisher, 2007, Gerrish and Lenski, 1998, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008, Tsimring et al., 1996.
The role of nonlinearities and fitness epistasis is again taken up in the Discussion. 3) As above, for "Theory of phenotypic interference" first para, prepare the reader how the local and global fitness dynamics will be coupled, i.e., by mentioning specifically that the local part will assume someσ, which will be then set by the global fitness wave part. This I understood only after reading through the whole manuscript first and being confused.
To take this point into account, we have reorganized the entire theory section to make the logical structure clearer. The introductory theory paragraph now reads The solution of the minimal model depends on a self-consistent closure of the local and global evolutionary dynamics, which are illustrated in Box 1. We will discuss these components in order. First, we study the evolutionary equilibrium of an individual biophysical trait in a process with given coalescence rateσ. Such traits evolve in a low-mutation regime characterized by strongly peaked trait distributions within a population and strong fluctuations over time and across genes (Supplementary Figure S1a). We show that the single-trait equilibrium depends in a simple way on the global parameterσ. Second, for the steady state of housekeeping evolution, the fitness variances of all traits combine to a global fitness wave (Supplementary Figure S1b). In tune with fitness wave theory, the wave pattern is of approximately Gaussian form and reflects broadly distributed, standing fitness variation maintained by mutations, selection, and coalescence [Desai and Fisher, 2007, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008. The global fitness wave and, hence, the coalescence rateσ, are determined by the ensemble of traits evolving under genetic linkage. Conversely,σ determines the single-trait equilibria, leading to a closure of the derivation. 4) Minor, in the subsection "Evolution of a quantitative trait under interference selection" -that subsection is just a general section on the evolution of a quantitative trait, as the interference doesn't come in yet, correct? Perhaps you can amend the title.
Also: can you make clear what parts of this subsection are novel (e.g., relations in Eq (1) and (2) appear straightforward, has this been done before)?
We agree that the scaling relations of an individual trait as a function of the coalescence rateσ = N −1 e hold regardless if this quantity is set by genetic drift (as is the case for recombining populations) or by selection (which applies to fitness waves). We now state this explicitly: In a fitness wave, the parameterσ couples each individual trait to the global evolutionary dynamics of all genetically linked genes (Box 1). In contrast, an independently evolving trait depends on an effective population size N e set by genetic drift.
The specifics of interference selection enter through the derivation of the diversity formula (6), which is given in Methods. It turns out that the "effectively neutral" form of the diversity applies throughout the phenotypic interference regime even for traits with many sequence sites, while it would be violated in large populations governed by genetic drift. We also agree that the calculations in this section are quite straightforward and extend previous results for neutral sequence diversity [Gillespie, 2000, Good et al., 2014, Kimura, 1983, Rice et al., 2015 and for quantitative traits under genetic drift [Charlesworth, 2013, Keightley and Hill, 1988, Lynch and Hill, 1986, Nourmohammad et al., 2013, as stated in the text. The new part is the derivation of eqs. (1), (3) for traits under selection in the phenotypic interference regime; see eqs. (16) -(17). Because interference selection is quite different from genetic drift, this is not clear a priori, and it is an important ingredient for the entire self-consistent solution. We have rewritten this section to make the flow of arguments clearer. 5) I had the most trouble parsing the crucial section "Housekeeping evolution of multiple traits". First, it would be good to be more clear what are the necessary conditions for the fitness wave regime to hold in this case. Make it more explicit that the fitness wave theory for your case contributes the crucial equation that relates the fitness variance and the coalescence rate, and perhaps give a few sentence intuition behind it.
We have rewritten the text according to your suggestion: In the housekeeping state, the total fitness variance σ 2 is simply the sum of the fitness variances Δ f of the individual genes (Supplementary Figure S4). Using Equation (3), this sum rule takes the form σ 2 = gΔ f = 2ugσ, which relates the scales of global selection and coalescence, σ and σ. Given a sufficient supply of non-neutral mutations, global evolution proceeds in a fitness wave (the condition for wave occurrence will be made precise below). General fitness wave theory then provides a second relation between global selection and coalescence, where N is the population size and c 0 is a model-dependent prefactor Hallatschek, 2013, Neher et al., 2013] (Methods). Combining these relations, we obtain the global fitness wave of phenotypic interference, ...
We also added a new Methods section "Housekeeping equilibrium and fitness waves" that provides some background on fitness waves including short derivations, to make the text more self-contained. Moreover, currently you assert c 0 ∼ 100, which requires at least a few words about where the magic constant emerges from. The prefactor c 0 in eq. (4) is model-dependent and known only in asymptotic cases, which is now explained in Methods, "Housekeeping equilibrium and fitness waves". For this reason, we treat c 0 as a fit parameter; this has been clarified in the main text and in Fig. 1. Second, you refer to c as the "slowly varying parameter". In what sense is this slow, can you be more precise?
This refers to the logarithmic dependence of c on evolutionary parameters. We have made this explicit in eq. (7) and in Methods, "Housekeeping equilibrium and fitness waves". Third, can you comment briefly about how Eq 5 is consistent with the σ 2 /σ 2 = c0 log(Nσ) from the fitness wave theory?
The step is now made explicit in the main text, from eq. (4) to eq. (7), and is further explained in Methods. It involves a series expansion. From (4), we have

This equation can be solved analytically
Fourth, you say "As sown in Methods, this parameter [c] has a simple interpretation:..." I think a few sentences more here in the main text would help, c is the crucial quantity. Thanks, we have added this to the main text below eq. (7), to Box 1, and to Methods, "Housekeeping equilibrium and fitness waves". The parameter c is indeed ubiquitous: in any fitness wave, it relates both the coalescence time and the total fitness span to the fitness variance (Box 1 and Methods). Under phenotypic interference, it is related in a specific way to the wave complexity (Box 1 and eq. (19)). In a nutshell, c measures how many mutational steps (of size ∼σ) are required to cover one fitness spanσ. Fifth, you say "The calibration between theory and data involves the fitting of a single global amplitude c 0 "; going back to the above assertions about c 0 , can you give some intuition of what is the meaning of this parameter being fitted?
As discussed above, c 0 is a model-dependent prefactor in the general fitness wave relation (4). Therefore, we treat it as a global fit parameter, which is now explicitly quoted in Fig. 1. 6) Can you comment on the sharpness of the transition in Fig 5a? The derivation in the section "The pathway to sexual evolution", which predicts a first order transition, is based on infinite system size. In a finite system, first order transitions are always rounded. In the present system, there are two sources of finite-size scaling: the number of genes g, and the population size N . As expected, the simulation results are rounded consistently for different gene numbers. We now comment on this directly in Fig. 4a. 7) There is an interesting analogy that you might want to highlight between population genetic interference discussed in this paper and crosstalk interference (also leading to superlinear degradation of performance in signaling networks) discussed in Friedlander, Prizak et al, Nat Comms 7 (2016); in both cases, it is crucial to self-consistently consider the functioning / evolution of one component in a background of others that are structurally similar. A followup paper by the same authors (Nat Comms 8, 2017) is a further example of nontrivial evolution in sigmoidal biophysical fitness landscapes (refs 24-30 in your paper) of the same kind as the stability/affinity in your paper.
Thanks for pointing us to these papers. Indeed the analogy to [Friedlander et al., 2016] is very interesting and we mention the paper in the discussion. Moreover, we list the fitness landscape of [Friedlander et al., 2017] as another example for the model of the fitness landscape.

Reviewer #2
This manuscript studies clonal interference dynamics in a model of quantitative biophysical traits, focusing particularly on how such dynamics may constrain the size of genomes. I consider the general problem interesting and their approach appealing. However, I did not gain many new insights from their approach, but rather found their results mainly technical. Therefore as it stands, I suspect the paper will be accessible only to specialists in theoretical population genetics who will be interested in the model. Please see below for lists of my general and specific comments.
We thank the reviewer for assessing our manuscript and giving valuable feedback. In response, we have rewritten the paper, and we hope the revised version highlights the biological results in a better way. In the Discussion, we recap four novel points: the new picture of interference selection at the phenotypic level, the systems-biological implications for evolution experiments and for macro-evolution, and the specific pathway to recombination discussed in the paper. Throughout the paper, the role of the underlying biophysics for these results is explained better. General comments: 1. The paper's main novelty seems to be studying quantitative traits in the regime of clonal interference, thereby introducing the concept of "phenotypic interference." It was not clear to me what is gained from this approach beyond the existing body of work on clonal interference. The authors imply that "phenotypic interference" is a fundamentally new concept, but to me it just seems like clonal interference applied to quantitative traits.
We now address this point comprehensively in several places of the manuscript. This paper is the first to combine fitness wave theory and the theory of quantitative traits. While asexual evolution in global fitness waves and evolution on biophysical fitness landscapes (e.g., for protein stability) are individually well-studied subjects, the merger of these into a consistent model produces qualitatively new results. The main reason is the distinct selective forces that drive the fitness wave in our model: we use mechanistically grounded biophysical fitness landscapes instead of ad hoc assumptions on the sprectrum of site selection coefficients. This is now made explicit already in the introduction: From the perspective of molecular evolution, these landscapes provide a biophysically grounded model of local fitness epistasis, which couples all sequence sites contributing to a stability or affinity trait in the same gene. Importantly, genome-wide local epistasis in protein-coding sequence operates independently on any assumptions of fitness interactions across genes. Beyond proteins, local epistasis occurs ubiquitously in quantitative molecular traits associated with binding interactions. This form of epistasis will be an important building block of our dynamical analysis but is not covered by the standard theory of asexual evolution [Desai and Fisher, 2007, Gerrish and Lenski, 1998, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008, Tsimring et al., 1996. .... Selection against genome complexity can be traced to a specific evolutionary mode of phenotypic interference, which is shown to occur in a parameter regime appropriate for typical microbial systems. The key characteristics of this mode are summarized in Box 1. First, local quantitative traits of a given gene are in an evolutionary equilibrium where the long-term average of the trait value, and its position on the fitness landscape, are determined by the uphill force of selection and the downhill force of mutations (Box 1, top). Second, global genome evolution takes place in a so-called fitness wave; that is, genetic and phenotypic variants in multiple genes co-exist in a population and generate a broad distribution of fitness values [Desai and Fisher, 2007, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008 (Box 1, center). In a fitness wave, selection rather than genetic drift determines the coalescence rate of a population, which is the inverse of the effective population size , Rice et al., 2015. Third, the defining feature of phenotypic interference is that global and local evolution are coupled in a specific way: the global coalescence rate sets the average selection coefficient of local trait changes. This feedback and its consequences are not covered by the existing models of fitness waves [Desai and Fisher, 2007, Good et al., 2012, Hallatschek, 2011, Neher and Hallatschek, 2013, Rouzine et al., 2008 and of quantitative trait evolution [Berg et al., 2004, Chi and Liberles, 2016, Gerland and Hwa, 2002, Goldstein, 2011, Manhart and Morozov, 2015, Nourmohammad et al., 2013, Serohijos and Shakhnovich, 2014.
In addition, Box 1 summarizes the salient features of phenotypic interference in a non-technical way.
The new results are now succintly captured in the discussion. In particular, the nonlinear scaling of the load and its implications, as well as the specific pathway for the evolution of sex, have not been reported before. This pathway builds again on mechanistically grounded epistasis; its differences to previous models based on fitness waves , Weissman and Barton, 2012 and to models based on global epistasis are now clearly described.
I furthermore found the emphasis on biophysical traits like protein stability misleading, because it didn't seem important to the authors' results: they didn't apply their results to any protein evolution data, nor did they make the argument that destabilizing mutations to proteins should be of dominant importance. Is the point just to use any nonlinear, diminishing-returns fitness function? If so, I think it would be clearer to consider abstract traits with such properties, making the argument that these fitness functions should be common among molecular traits (with protein stability being just one example).
Thank you for this comment, which has been instrumental for the revision. We now explain the role of biophysics in a better way throughout the manuscript. In a nutshell, this role is not just to give examples of molecular quantitative traits, but also to provide mechanistically grounded models of epistatic selection. A cell does have about 5000 proteins, so the genome-wide selection map does include about 5000 local fitness landscapes of sigmoid form for fold stability alone. As we show in this paper, the ensemble of these landscapes leads to different evolutionary regimes than previous models of asexual evolution based on abstract assumptions on selection. Without the biophysical underpinning, our epistatic model would be just as good as any other assumption.
Regarding the specifics of protein evolution, we use experimental data in our model, largely informed by refs. , Serohijos and Shakhnovich, 2014. Evidence for the pervasive role of protein stability in molecular evolution has come from a number of recent experimental works. We discuss in detail one example, the analysis of Lenski's long-term experiment, in light of our model [Couce et al., 2017, Good et al., 2017, Leiby and Marx, 2014, Tenaillon et al., 2016. But the thrust of the macro-evolutionary part of the paper is indeed somewhat different: instead of arguing for the importance of protein evolution compared to other molecular units, we show that already a minimal model with only housekeeping evolution of proteins has strong implications for genome complexity. More comprehensive models with more molecules and with ecological adaptations, which doubtlessly are relevant in parallel to protein stability evolution, would only increase these effects. But even the minimal model needs a biophysical underpinning for credibility. 2. The main result of this paper seems to be predicting how clonal interference of quantitative traits constrains genome size. This is of course an important problem, but there is a long history of studies addressing it, and I did not feel that literature was adequately discussed and compared to the present result. This made it hard for me to tell how significantly different the authors' results are. In particular, constraints on genome size resulting from accumulation of mutations destabilizing protein folding was discussed by Zeldovich et al. (PNAS 2007). As far as I can tell, this paper is cited in the Methods, but not discussed in the context of the results.
The paper by  and follow-up papers are seminal works that informed much of our thinking about protein evolution. We believe we have carefully cited this group of papers throughout the manuscript; the relation between  and our work is now discussed at several places of the main text. In a nutshell, we use exactly the same "mutational wind" as  but there are two main differences between  and our work 1 : (i) We treat the evolution of genetically linked proteins under interference selection, while  assume factorization across genes, which amounts to free recombination: "Because mutations in all proteins are independent, and boundary conditions are the same in each dimension (i.e., for each gene), one can write P (E) = Γ i=1 p(E i ), separate the variables, and reduce the problem to a product of Γ 1D diffusion problems for each gene: [...]" . Their paper mentions possible deviations from factorization at higher mutation rates but does not treat interference selection between these mutations.
(ii) We show that in the phenotypic interference regime, individual traits evolve in a monomorphic lowmutation regime, while  assume a highly polymorphic quasispecies (the solution of an imaginary-time Schrödinger equation) for each trait (the regime Nm > 1 in their notation, see also reply to point 6).
These differences translate into qualitatively and quantitatively different results. Linked evolution of quantitative traits produces a superlinear load, L ∼ ug 2 (data and red line in Fig. 2);  find a linear load, L ∼ ug, which we recover for recombining populations (brown line in Fig. 2). The difference between these can be orders of magnitude. 3. Overall the paper is quite technical and difficult to read. In particular, the main text contains quite a lot of notation (some of which I think is nonstandard and possibly unnecessary -see specific comments below).
The authors should consider significant simplifications to the presentation of their results if they hope for them to be accessible to an audience beyond theoretical population geneticists.
We addressed this comment throughout the paper. In particular, Box 1 provides a non-technical summary of phenotypic interference aimed at a broad readership. Throughout the main text, we try to keep technicalities at a minimum, while the Methods section now provides more background to make the paper self-contained. The biological implications section can now be read with just Box 1 as background. Specific comments: 1. It would be helpful for the authors to include line numbers in their manuscript -this would make it easier to refer to specific items in the text. I will do my best below to indicate where in the text my comments apply.
We are sorry for the inconvenience and added line numbers to the new version.
2. In my opinion Fig. 1b comes too early in the text -the introduction, at which point I was definitely not prepared to understand it. Shouldn't this plot be part of a main figure in the Results section, where it can be fully explained?
In the revised version, this figure has been trimmed and is now a schematic within Box 1. We think it is important to give this take-home message at an early point in the manuscript. All quantitative analysis is in the later figures, which come in parallel with the detailed explanation. 3. The colors in Fig. 2 are called "blue" and "orange," although for me they appeared to be more of a teal and red. I'm not sure if this is a difference in software/hardware or my own vision, but the authors might double-check these colors in case readers get confused.
Thanks, we will make sure that the final version of figures has unambiguous colors. 4. The main text introduces a lot of new notation, some of which appears to be nonstandard and unnecessary. For example, they use G for the stability of a protein rather than the more standard ΔG, which is compounded by defining Δ G as the variance. Also, why use introduce Γ instead ofḠ, andσ instead of N −1 e ? To me this made many of the results confusing, because these expressions fundamentally involve classical combinations of selection coefficients, mutation rates, and effective population sizes, but their notation obscures this intuition.
We are well aware of other notations and their conflicts, especially with the community of protein evolution. However, there is a trade-off between the notations between scientific communities: ΔG would collide with the Δ widely used for diversities in population genetics. For example, Δ ΔG from population genetics collides with ΔΔG for mutational effects in protein evolution. Denoting the mean trait Γ byḠ can again cause confusion, for example in the probability distribution Q(Γ), which is the cross-species distribution of the within-species mean trait. Similarly,σ = N −1 e is the inverse of the effective population size, but it is fundamentally a selection scale and not related to genetic drift as the classical notation N e may suggest. To emphasize the links between selection scales, e.g. in Box 1, we prefer to use a notation as introduced in [Schiffels et al., Genetics 2011;Rice et al., Genetics 2015]. However, we added more explanatory text and alternative notations throughout the text. 5. The authors make a lot of mathematical statements without adequate explanation. For example, the authors need to provide conceptual explanations for why Δ G f (Γ) is the rate of stability increase by selection, and where the mutation-selection balance condition (Δ G f (Γ) = u G ) comes from. Full derivations in the Methods/SI were also not cited for these.
There is a trade-off between readability and level of detail in the main text. To address your comment, we added several more detailed derivations to Methods, including a careful derivation of the mutation-selection balance condition (Δ G f (Γ) = u G ) in the section "Evolutionary model". We also linked the explanations in Methods to the main text in a better way. 6. In the section called "Evolution of a quantitative trait under interference selection," where is the effect of linkage? I couldn't tell the difference between these results and what one would expect for an unlinked trait evolving independently.
We agree that the relations (1) -(3) would formally take the same form in standard population genetics, if we use the identification N e =σ −1 . This is now clearly stated in the text: This form extends previous results on neutral sequence diversity [Gillespie, 2000, Good et al., 2014, Kimura, 1983, Nourmohammad et al., 2013, Rice et al., 2015 and on traits under genetic drift [Keightley and Hill, 1988, Lynch and Hill, 1986, Nourmohammad et al., 2013.
The question is whether this form is actually valid for traits under diminishing-return selection and genetic linkage. This is far from obvious and has, to our knowledge, not been established before. In Methods, we show that the form (1) for the diversity is self-consistently valid in the phenotypic interference regime. Here is where the specifics ofσ as a selection scale independent of genetic drift come in. In contrast, Δ G = 2N e u 2 G would often cease to be valid for complex traits evolving independently. This is the case, for example, in the regime of , where Δ G is shaped by selection (in the notation of , the trait variance is the width of the distribution p(ΔG) in their eq. [6a]). Our main text states this: In Methods, we derive Equation (1 for quantitative traits in a fitness landscape f (G) by showing that stabilizing selection on Δ G can be neglected throughout the phenotypic interference regime; this scaling is confirmed by simulations (Supplementary Figure S2a). In a fitness wave, the parameterσ couples each individual trait to the global evolutionary dynamics of all genetically linked genes (Box 1). In contrast, an independently evolving trait depends on an effective population size N e set by genetic drift.
7. "The strength of selection on genetic variants is not fixed a priori, but is an emergent property of the global evolutionary process. A faster pace of evolution, i.e., an increase in coalescence rateσ, reduces the efficacy of selection." This sounds like restating a very standard concept in population genetics, that absolute selection strength doesn't matter so much as its value relative to competing processes like genetic drift (coalescence) and mutation. If the authors are implying something beyond this, then they should clarify.
The statement does say more than standard population genetics. The standard concept you refer to is that selection coefficients enter fixation probabilities via the product with the effective population size, 2N e s = s/σ. Eq. (1) states that this product is tuned to a fixed average value, 2N e s = 2.
The text below eq. (1) gives a non-technical explanation of this result. We have rewritten this text to make it clear that the standard concept is just one part of the argument. 8. "The non-adaptive, stationary scenario of housekeeping evolution..." I think calling that scenario "nonadaptive" is misleading: it is a steady state, so the mean fitness doesn't change, but there is still selection. To me "non-adaptive evolution" is devoid of any selection, e.g., purely genetic drift, mutation, recombination.
There are different definitions of adaptation in the literature. We clarified in revised the text that we use a standard definition of evolutionary and population biology; see e.g., Kreitman, 1991, Orr, 2005]. Here adaptation is defined by a net surplus of beneficial mutations and is clearly distinguished from evolutionary equilbrium. According to that definition, the housekeeping state has no adaptation because the average protein stability remains constant in time. 9. Above Eq. 3, the parameter c 0 appears and is given a value but with no explanation. Is it related to c in Eq. 5?
c 0 is a model-dependent prefactor in the fitness wave relation (4). We have clarified the main text of this section; in Methods "Housekeeping equilibrium and fitness waves", we provide more background on eq. (4). 10. In general Eqs. 3, 4, 5 are presented without adequate explanation. I gather they define properties of the fitness wave, but the authors need to introduce fitness waves more systematically and explain how they are parameterized.
Thanks for this suggestion. In response, the section "Housekeeping evolution of multiple traits" has been completely rewritten to improve clarity. In addition, the key fitness wave relation (4) is explained in Methods, "Housekeeping equilibrium and fitness waves". Box 1 gives a non-technical summary of all this. 11. Figure 3 purports to show the two scaling regimes derived in the text between independently-evolving genes and the linkage-dominated case. However, the plots are not convincing to me that their analytical results are particularly accurate -the data does not match either regime all that well, nor does it show a clear crossover. This is especially concerning given that the data is from simulations, so whatever approximations the authors use in their analytical results may be problematic. Figure 1 (which was Fig. 3 before) shows the data (dots) and two theoretical scaling regimes: independently evolving genes (brown) and genes under phenotypic interference (red). We are not sure if we understand your comment since the data do interpolate between these theoretical curves. In particular, there is little doubt that the data show strong deviations from the brown line of independently evolving genes. The same is true for the genomic load in Fig. 2, which deviates by 2 orders of magnitude from the linear expectation (brown line) and is well approximated by a quadratic form. So the main point of the paper, the superlinear form of the genetic load, is unambiguosly supported by Fig. 1 and Fig. 2.
Moreover, there is a number of important consistency checks: (a) In line with theory, the data decouple well from the selection amplitude f 0 of the gene. This has been tested for 4 orders of magnitude of selection coefficients (color code). (b) The crossover point g 0 ∼ c is shown consistently in all panels (red dashed line). (c) Even relative amplitudes are reproduced consistently over all panels of Fig. 1 and Fig. 2, there is only one global fit parameter c 0 (see caption of Fig. 1). 12. The text leading up to Eq. 6 should precisely define genetic load before introducing results on it.
Thanks, we now carefully define the genetic load: The evolutionary cost of deleterious mutations can be quantified by the genetic load, which is defined as the mean fitness of a population compared to the fitness maximum. In the biophysical fitness landscape f (G) of the minimal model, the load of a given gene takes the approximate form f 0 − f (Γ), where Γ denotes the population mean stability and f 0 is the fitness of a fully stable gene (G 0); see Box 1 (top) and Equation (13) in Methods.
13. At the bottom of page 6, the authors start mentioning the effect of all sorts of other scenarios on the genetic load result (Eq. 6), but this is way too terse to understand, and there was no reference to the SI which I later realized fully explains these cases. We have extended this paragraph and improved the connection to the full explanations in SI: In Supplementary Notes and Supplementary Figure S5, we discuss phenotypic interference in a number of extended biophysical models. These include active protein degradation at the cellular level, a ubiquitous process that drives the thermodynamics of folding out of equilibrium [Hochstrasser, 1996]. We show that thermodynamic non-equilibrium influences the fitness landscape but does not affect our qualitative and quantitative conclusions. Another example is the stability-affinity model, which has two quantitative traits per gene that evolve in a twodimensional sigmoid fitness landscape f (G, E) with both traits under (correlated) stabilizing selection [Chéron et al., 2016, Manhart andMorozov, 2015]. We show that under reasonable biophysical assumptions, evolution in a stability-affinity model produces a 2-fold higher interference load than the minimal model, L int (g) ≈ 8ug 2 /c. Alternative models with a quadratic single-peak fitness landscape describe, for example, gene expression levels under stabilizing selection [Nourmohammad et al., 2017]. We show that such landscapes generate an even stronger nonlinearity of the load, L int (g) ∼ g 3 . In contrast, a discrete model with a fitness effect f 0 of each gene shows a linear load up to a characteristic gene number g m = (f 0 /u) log(Nf 0 ) associated with the onset of mutational meltdown by Muller's ratchet [Gordo and Charlesworth, 2000, Muller, 1964, Rouzine et al., 2008. These examples suggest that superlinear scaling of the genetic load holds under quite general conditions, given a sufficient number of quantitative traits evolving under genetic linkage and in fitness landscapes with negative epistasis. This type of fitness landscape is ubiquitous in biophysical models.
14. I was surprised to see no mention of Neher and Shraiman (PNAS 2009) in the section on "The pathway to sexual evolution." I think that paper, which also shows the transition between linkage-dominated and recombination-dominated cases, is an important comparison to this manuscript's results.
Thank you for pointing out this reference. It gives another example for a crossover from selection on genotypes towards selection on genes; related crossovers are treated in the cited refs. . Model and results of [Neher and Shraiman, 2009] clearly differ from ours. While we discuss a long-term evolutionary mode with the competition of trait mutations, [Neher and Shraiman, 2009] discuss the loss of diversity from an initial standing variation. Therefore, no new variants arise there and they describe short-term dynamics. This scenario cannot describe the evolution of recombination. Moreover, they show for the short-term fixation that small recombination rates may brings a small benefit, but the critical transition towards selection on genes instead of genotypes comes with a cost, see Fig. 4b of [Neher and Shraiman, 2009]. This is a property of heterogeneous global epistasis, which acts very differently from the local epistasis of biophysical fitness landscapes. We added a reference to [Neher and Shraiman, 2009] and a comparative discussion to the section "Comparison of models for the evolution of recombination": In another model with broad, heterogeneous pairwise epistasis, recombination generates a transition from genotype selection to gene selection [Neher and Shraiman, 2009]. This causes a shortterm benefit, but the transition causes a net fitness loss by breaking linkage between epistatic sites [Neher and Shraiman, 2009]. Our model replaces the assumption of global epistasis across the genome by local diminishing-return epistasis between loci affecting the same protein stability or affinity trait. This form of epistasis follows directly from the underlying thermodynamic nonlinearities, and its evolutionary effects differ from global epistasis.
15. "The spectrum of site selection coefficients is not a fixed input, but a dynamical output of the evolutionary process. That is the main difference of our approach to previous population-genetic models of asexual evolution [24][25][26][27][28][29][30]." As I understand it, this statement does not seem to be supported by the references or the preceding text. The citations mentioned here also have nonlinear fitness functions, where the strength of selection on a trait depends on the current value of the trait. The main difference, I would say, is a lack of clonal interference.
16. The parameterization of the fitness function in Eqs. 10 and 11 does not seem to allow for essential proteins or those that are otherwise lethal when unfolded. As long as f 0 is finite, then even G → −∞ for a single protein will not result in zero total fitness (Eq. 11).
We added a comment to Methods on how our approach takes into account essential genes. As introduced in Shakhnovich, 2009, Zeldovich et al., 2007], essential genes (for which lack of folding is lethal) can be modeled by a sharp cutoff of the fitness landscape. That is, our Malthusian fitness landscape f (G) acquires a singularity in the regime of marginal folding (f → −∞ for G → G 0 ∼ 0). Importantly, however, the landscape f (G) retains the form (13) in the regime of stable folding (G/k B T 1), which implies that our conclusions remain unaffected. This is also consistent with our results on gene loss, Equation (10) and Fig. 3. We show that loss affects weakly selected genes, while essential genes are protected from loss by strong selection.
The coalescence rate is the inverse of the average coalescence time, i.e. the time from the common ancestor of the population. It is the timescale that defines the scale of effective selection . We improved the explanation of this term in the Introduction and in Fig. 1. 8. I did not understand the derivation of Eq. 8. The authors say that for a stable wave, mutation rate must exceed average selection, so ug > s. But if I substitute in Eq. 6 (s = 4ug/c) as the authors say, I get ug > 4ug/c, which leads to a condition on c (c > 4) not g, since the g drops out.
Thanks, we corrected this inconsistency and added a straightforward derivation of the onset threshold in Methods (lines 580 -585). 9. Lines 234-235: The previous sentences indicate the mutational load L for an individual gene can be approximated as f 0 − f (Γ) ≈ f 0 e −Γ/kT . But then this sentence says I should use Eq. 6 -how? That equation does not have any relation for Γ.
We improved the text (now lines 239 -242). The key point is that the product s = 0 f has been evaluated in equation (6). 10. Lines 284-285 assert that the physiological mutational load is similar to the interference load using E. coli parameters, but don't the numbers show the interference load is much larger? The estimate for the interference load in microbes in 3e-2 (line 267), but the physiological load estimate seems to be 1/5000 2e-4 (lines 283-284).
As we specify in line 266, the cost in line 267 is for a 10 % more genes, i.e. for 500 additional genes. Therefore, the costs per gene are of the same order of magnitude. 11. Lines 294-295 mention the "stability condition f (G)" -shouldn't there be an inequality here?
We added a reference to the detailed balance condition, Equation (2). Moreover, we specified the discussion of (in-)stability regimes below with respect to f 0 , which clarifies the condition and the equation (10). 12. Lines 307-308 and Eq. 11: What is the logic from "each additional deleterious trait changes are only marginally selected" to "the relaxation time is of the order of the inverse mutation rate"?
The logic is, marginally selected mutations evolve in a a near-neutral way. Therefore, the fixation rate of these mutations becomes similar to the mutation rate. We specified the intermediate argument in the text. 13. Figure 3 caption: I think the two parts of panel (a) should be marked "top" and "bottom," not "left" and "right." Thanks, we fixed this. 14. Line 347 includes the relation Rξ/gσ (ξ). Is this a definition of ξ? Or an independent constraint on R and ξ?
This relation is discussed in detail in refs. , Weissman and Barton, 2012. It is a self-consistent relation between block size and recombination following from simple time-scale arguments, which is now specified in the text. 15. Line 354 mentions that there is a first-order phase transition from genotypic selection at low recombination rates to allelic selection at high recombination rates, but I didn't see any explanation or evidence for this being a sharp transition.
The evidence comes from analytical considerations in the section "The pathway to sexual evolution" predicting the first order transition. We improved the discussion of these analytical results. In particular, in a finite system, first order transitions are always rounded. In the present system, there are two sources of finite-size scaling: the number of genes g, and the population size N . As expected, the simulation results are rounded consistently for different gene numbers. We now comment on this directly in Fig. 5a. 16. I'm still skeptical about the mathematical formulation of the fitness function of protein stabilities in Eqs. 13 and 14, especially as it applies to essential genes. The authors say it can account for essential proteins by taking f 0 → ∞, but that will make such a protein always lethal for any finite stability. (Also, it seems like an undesirable parameterization of the model to make essential proteins require some sort of infinite limit.) We improved the text in the Methods section "Biophysical fitness models" to make it clear that we do not model essential genes by f 0 → ∞. Rather, zero growth (Wright fitness w = 0, lethality) occurring at a finite 4 threshold free energy G 0 translates into a singularity of the Malthusian (log) fitness, f (G) = log w(G) → −∞ for G → G 0 at finite f 0 . This is similar to the modeling of essential genes in Shakhnovich, 2009, Zeldovich et al., 2007]. Our simulations with essential genes are done with this well-defined fitness. I think this problem is related to their choice of adding fitness functions for each protein to get the total fitness of the genotype, since I would argue sigmoidal fitness functions are more naturally multiplicative. I think a better formulation would be the following. Let f (G) = (1 + (1 − s)e −G/kT )/(1 + e −G/kT ), so a perfectly-stable protein has maximum fitness 1, and a perfectly-unstable protein has fitness 1 − s. That is, s is the fitness cost of the unfolded protein, which ranges from 0 (neutral protein) to 1 (essential protein). The fitness for the entire genotype is the product of these fitness functions for each protein. Therefore if an essential protein becomes completely unstable, the fitness for the entire genotype becomes zero.
Thanks for pointing to this interesting model of selection on essential genes, which we included in the section "Biophysical fitness models", translating it into the language of Malthusian (log) fitness. We then recall a salient point of our analysis: the mutation-selection balance drives the average load per gene to a fixed value 2σ independently of the gene selection amplitude f 0 . That is, the equilibrium load of essential genes is no different from others, the larger f 0 is compensated by a larger equilibrium G. Therefore, the model assumption (14) of multiplicative fitness across genes, i.e., additive Malthusian (log) fitness, is consistent in this regime.