Global epistasis emerges from a generic model of a complex trait

1Despite the idiosyncrasies of epistasis, a number of laboratory microbial evolution experiments show systematic patterns of convergent phenotypic evolution and declining adaptability. A striking example is provided by the Escherichia coli long-term evolution experiment (LTEE) (Figure 1a): 12 replicate populations that adapt in parallel show remarkably similar trajectories of fitness increase over time, despite stochasticity in the identity of fixed mutations and the underlying dynamics of molecular evolution. Similar consistent patterns of fitness evolution characterized by declining adaptability over time have also been observed in parallel yeast populations evolved from different genetic backgrounds and initial fitnesses (; Figure 1b) and in other organisms. Declining adaptability is thought to arise from diminishing-returns epistasis, where a global coupling induced by epistatic interactions systematically reduces the effect size of individual beneficial mutations on fitter backgrounds. Diminishing-returns manifests as a striking linear dependence of the fitness effect of a mutation on background fitness (Figure 1c). While diminishing-returns can be rationalized as the saturation of a trait close to a fitness peak, recent work shows a similar dependence on background fitness even for deleterious mutations, which become more costly on higher fitness backgrounds. This suggests that fitter backgrounds are also less robust to deleterious effects (Figure 1d), a phenomenon that has been termed increasing-costs epistasis. The origin of the global coupling that results in these effects is unknown.

2Put together, these empirical observations suggest that the contributions to the fitness effect, si, of a mutation at a locus in a given genetic background can be written aswhere is the additive effect of the mutation, is its genotype-dependent epistatic contribution independent of the background fitness (i.e., idiosyncratic epistasis), and ci quantifies the magnitude of global epistasis for locus. Equation (1) reflects the observation that the strength of global epistasis depends on the specific mutation and applies independently of whether its additive effect is deleterious (increasing-costs) or beneficial (diminishing-returns). Over the course of adaptation in a fixed environment, global epistatic feedback on mutational effects can lead to a long-term decrease in adaptability. If this feedback dominates, Equation (1) suggests that the dependence of the fitness effect on evolutionary history is summarized entirely by the current fitness, and therefore results in predictable fitness evolution.

3Here, we show that diminishing-returns and increasing-costs epistasis are a simple consequence of widespread epistasis (WE). This is consistent with recent work that proposes a similar argument to explain these phenomena. However, while the core idea is similar, we present here an alternative framework based on the Fourier analysis of fitness landscapes, which leads to new insights and quantitative predictions. In particular, our framework leads to novel predictions for the relationship between the magnitude of global epistasis and the stochastic effects of microscopic epistasis, which we confirm by reanalyzing existing data. Extending this framework, we further quantify how the distribution of fitness effects (DFEs) shifts as the organism adapts and how the fitness effect of a mutation depends on the sequence of mutations that have fixed over the course of adaptation (i.e., historical contingency). While specific historical relationships depend on the genetic architecture, we introduce a novel fitness landscape model with an intuitive architecture for which the entire history is summarized by the current fitness. Using this fitness landscape model, we investigate the long-term dynamics of adaptation and elucidate the architectural features that lead to predictable fitness evolution.

4We begin by examining the most general way to express the relationship between genotype and fitness (i.e., to describe the fitness landscape). A map between a quantitative trait (such as fitness),, and the underlying genotype can be expressed as a sum of combinations of biallelic loci that take on values:where is a constant that sets the overall scale of fitness. The symmetric convention for the two allelic variants is less often used than, but it is an equivalent formulation, which we employ here because it will prove more convenient for our purposes (see for a discussion). The coefficients of terms linear in xi represent the additive contribution of each locus to the fitness (i.e., its fitness effect averaged across genotypes at all other loci), the higher-order terms quantify epistatic interactions of all orders, and is the average fitness across all possible genotypes. Importantly, Equation (2) makes apparent the idiosyncrasies induced by epistasis: a mutation at a locus with interacting partners has an effect composed of contributions.

5To explicitly compute the fitness effect of a mutation at locus on a particular genetic background, we simply flip the sign of xi, keeping all other xj constant, and write down the difference in fitness that results. This fitness effect will generally involve a sum over a large number of terms involving the ’s in Equation (2). While this may suggest that an analysis of fitness effects via Equation (2) is intractable, the analysis in fact simplifies considerably if the locus has a significant number of independent interactions that contribute to the fitness (i.e., provided that the number of independent, nonzero epistatic terms associated to the locus is large). In this case, we show that the fitness effects of individual mutations decrease linearly with background fitness and the fluctuations around this linear trend are normally distributed. In other words, widespread independent idiosyncratic epistatic interactions lead to the observed patterns of diminishing-returns and increasing-costs epistasis.

6We present a derivation of this result in the SI. Here, we explain the key intuition using a heuristic argument. The argument is based on a simple idea: for a well-adapted organism with complex epistatic interactions, a mutation is more likely to disrupt rather than enhance fitness. To be quantitative, consider a highly simplified scenario where some number of the ’s in Equation (2) are ±1 at random and the others are 0. In this case, the fitness of a given genotype is a sum of and interactions that contribute positively and negatively to the trait, respectively, each with unit magnitude, so that. When positive and negative interactions balance, the organism is in a ‘neutrally adapted’ state. By selecting for positive interactions, adaptation generates a bias so that and. If locus involved in a fraction vi of all of interactions is mutated, the effect of the mutation, on average, is to flip the sign of positive interactions and negative interactions. The new fitness is then and thus. The negative linear relation between the background fitness,, and the fitness effect of the mutation, si, is immediately apparent and emerges as a systematic trend simply due to a sampling bias towards positive interactions. Of course, while this relation is true on average, it is possible that locus affects more or less positive interactions due to sampling fluctuations. Provided only that is large and the interactions are independent, these fluctuations are approximately Gaussian with magnitude.

7This basic argument holds beyond the simple model with unit interactions. In the more general case, if the mutation is directed from, we show in the SI that its fitness effect, si, on a background of fitness can be written aswhere,and is a genotype and locus-dependent term that is distributed across genotypes with mean zero and variance expressed in terms of the ’s from Equation (2) (see SI for details). In the following equation and similar ones henceforth, a summation such as is meant to denote a sum over pairs, where each pair appears only once and no pair that includes index appears. Symmetry of the ’s w.r.t. interchanged indices is also assumed (e.g., ). The numerator of in Equation (4) is proportional to the covariance of fitness effects and background fitness, and the denominator is the variance of background fitness across genotypes. A similar equation for the case can be derived. The choice of or is simply a matter of convention. If the convention is reversed, the coefficients of odd-order in Equation (2), that is,, should also switch signs. It can be easily checked that reversing the signs of these quantities in the expression for above leads to the expression for when.

8Note that in general is not guaranteed to be positive and is arbitrary and determined by the genotype-fitness map. However, consistent patterns emerge when locus has a large number of independent, nonzero epistatic terms and the additive effects of its interacting partners are not much larger than the epistatic terms (defined further below), which we call the WE limit. In the WE limit, is normally distributed across genotypes with variance proportional to. This follows from the same reasoning as in our heuristic argument with unit interactions above (see SI for details). In addition, is typically positive, giving rise to a negative linear trend (i.e., diminishing-returns and increasing-costs). We can see this by taking the third- and higher-order terms in Equation (4) to be zero, in which case is positive if. This will typically be true in the WE limit because we expect to scale with the number of interacting partners, while each term in can be positive or negative and thus the sum scales as if the terms are independent. Thus when locus has a large number of interacting partners, is typically positive unless the magnitude of the additive terms is much larger than the magnitude of the epistatic terms,. This argument is easily extended to the case when the third- and higher-order terms are nonzero (see SI); the upshot is that the bias towards positive gets stronger with increasing epistasis.

9The conditions for the WE limit are more likely to hold when the number of loci,, that affect the trait is large. Therefore, we expect to generically observe patterns of diminishing-returns and increasing-costs epistasis for a complex trait involving many loci. Importantly, whether we observe a negative linear trend does not depend on the magnitude of a locus’ epistatic interactions relative to its own additive effect, but rather relative to the additive effects of its interacting partners. If we are not in the WE limit, and instead the additive effects dominate (i.e., ), then Equation (4) suggests that the slope of the linear trend can be either positive or negative. We will show further below that recent experimental data demonstrates that both scenarios can be relevant: some loci have while others have, with the former creating a bias towards the observed negative linear trends that characterize diminishing-returns and increasing-costs epistasis.

10We note that Equation (3) immediately leads to testable quantitative predictions: in the WE limit, the distribution of the residuals,, obtained from regressing si and is entirely determined by the slope of the regression,. Specifically, we predict that these residuals (the deviations of individual genotype fitnesses from the overall diminishing-returns or increasing-costs trend) should be normally distributed with a variance proportional to. However, this condition only applies if diminishing-returns arises from the WE limit. It does not hold if epistasis is negligible, if locus interacts significantly with only a few other dominant loci, or if the epistatic terms are interrelated (e.g., when global epistasis arises from a nonlinearity applied to an unobserved additive trait;;; ). The latter case may still lead to a negative linear trend, but the statistics of the residuals will differ from Equation (3) (see SI for a discussion).

11It is convenient to subsequently work with the symmetric version of Equation (3), where the fitness effects of both and its reversion (whose fitness effect is negative of the former) are included in the regression against their respective background fitness. In this case, the additive term is averaged out, and we show (SI) that in the WE limit,where depends on the genetic background and the locus, and is normally distributed with zero mean and variance, and

12Here, is the total genetic variance due to all loci (i.e., the variance in fitness across all possible genotypes) while is the contribution to the total variance by the ’s involving locus. We therefore refer to vi as the variance fraction (VF) of locus. We show further below that for certain fitness landscapes vi can also be interpreted as the fraction of pathways affected by a locus. For these reasons, we focus on vi, which is half of the negative slope, rather than the slope. Note that the vi’s do not sum to one unless there is no epistasis (with epistasis,, reflecting the fact that the variance contributed by different loci overlap). While the directed mutation case discussed previously is the relevant one when presenting experimental data (e.g., Figure 1c, d), it is conceptually simpler to work with the symmetric case. These two cases coincide and in the WE limit if the additive effect of a locus is small (i.e., ).

13Our results show that the VF vi plays an important role. It determines the slope of the negative relationship between the fitness effect and background fitness. At the same time, it determines the magnitude of the idiosyncratic fluctuations away from this trend. We also note that this slope can be used to experimentally probe the contribution of a locus to the trait (i.e., its VF) taking into account all orders of epistasis, which circumvents the estimation of the individual ’s in Equation (2). The theory additionally predicts that the slope obtained by regressing the sum of fitness effects of two mutations at loci against background fitness is proportional to, where quantifies the magnitude of epistatic interactions of all orders between and (SI).

14Importantly, while the fitness effects of individual mutations (and hence the DFEs) may change over the course of evolution due to epistasis, the distribution of variance fractions (DVF) across loci,, is an invariant measure of the range of effect sizes available to the organism during adaptation. As we will see, this means that the DVF plays an important role in determining long-term adaptability.

15To illustrate our analytical results, we first demonstrate that the effects described above are reproduced in numerical simulations. To do so, we numerically generated a genotype-phenotype map of the form in Equation (2), with loci and an exponential DVF,, where (Materials and methods). This DVF is shown in Figure 2a. Note that corresponds to an epistatic landscape; chosen here thus corresponds to a model within the WE limit (note that in this parameter range). Using this numerical landscape, we measured the fitness effect of mutations at 30 loci across 640 background genotypes with a range of fitnesses (Figure 2b). Our results recapitulate the predicted linear dependence on background fitness (Figure 1c, d), with a negative slope equal to twice the VF predicted from Equation (5). We further simulated the evolution of randomly generated genotypes similar to the experimental procedure used in (Figure 2c), finding that our results reproduce the patterns of declining adaptability observed in experiments (Figure 1b). Note that ∼10 mutations are fixed during this simulated evolution; declining adaptability here is not due to a finite-sites effect.

16As described previously, Equation (5) implies a proportional relationship between the magnitude of global epistasis (quantified by the slope of the relationship between the fitness effect of a mutation and the background fitness) and the magnitude of microscopic epistasis (quantified by the residual variance around this linear trend); see also Figure 3a. We verify this relationship in simulations (Figure 2d). We predict that the slope obtained by regressing the sum of fitness effects of two mutations at loci against background fitness is proportional to. We further assume that (specifically, for the genotype-phenotype map used for numerics). Since vi and vj are typically small for a complex trait, we expect near-additivity and that any deviations are sub-additive, which is confirmed in simulations (Figure 2e, f).

17While testing the latter prediction on double mutants requires further experiments, we can immediately test the relationship between the slope and the distribution of residuals from existing experimental data. To do so, we reanalyzed the data from, which measured the fitness effect of 91 insertion mutants on about 145 backgrounds. These background strains were obtained by crossing two yeast strains that differed by single-nucleotide polymorphisms (SNPs). Of these 40,000 loci, have been identified as causal loci with currently available mapping resolution. In Figure 3, we show the estimated (negative one-half of the slope of the best-fit line) and the VF vi for each of the 91 mutations. These mutations were selected after screening for nonzero effect, and thus the DVF is biased upwards. The mean VF is. The wide range of vi observed in the data implies that the epistatic influence of loci varies greatly across loci, and we will show further below that this is crucial for maintaining a supply of beneficial mutations even when the organism is well-adapted to the environment.

18Our theoretical results imply that we expect the linear relationship between background fitness and fitness effect to be negative if the additive effects of a locus’ interacting partners are not much larger than the epistatic terms. Specifically, we define the additivity of interacting loci (AoIL) for locus aswhich we show can be estimated from data (Materials and methods and SI). If the AoIL is less than half, Equation (4) implies that the linear trend is guaranteed to be negative. If instead the AoIL is greater than 0.5, the trend can be either positive or negative. The data shows a range of AoIL between 0 and 1 across loci. As predicted by our theory, we find that the loci with AoIL lt0.5 always show negative trends and the ones with AoIL gt0.5 show both negative and positive trends (Figure 3c). Importantly, the sign of the trend is determined by the AoIL and not by the additivity of the mutated locus, which we define as

19The additivity across loci also has a wide range. However, small additivity does not necessarily imply a negative trend (Figure 3d).

20We next used the data from to analyze the relationship between the slope of the linear trend and the residual variance around this trend. We find that the experimental data confirms our theoretical prediction that the residual variance is proportional to ) if the AoIL is small (Figure 3e, for loci with AoIL ¡ 0.5 and for all loci). The Gaussian-distributed term in Equation (3) also predicts the shape of the distribution of the residuals given the VFs, which aligns well with the empirical distribution of the residuals (Figure 3f).

21Together, these theoretical results and our reanalysis of experimental data show that linear patterns of global diminishing-returns and increasing-costs epistasis are a simple consequence of widespread epistatic interactions. The DVFs observed in data (Figure 3b) further imply that the epistatic influence of different loci on fitness can vary across a wide range. In what follows, we show that these two observations can be put together to make general predictions about the DFEes, and consequently the long-term dynamics of adaptation. The key ingredient that enables this analysis (including Equation (5)) is that in the WE limit fitness and fitness effects are jointly normal (with respect to a uniform distribution over all possible genotypes), which allows us to quantify complex dependencies between these variables in terms of pairwise covariances.

22Long-term adaptation is determined by the DFEs of possible mutations and the stochastic dynamical processes that lead to fixation. While Equation (5) represents the distribution of the fitness effects of a specific mutation at locus over all genotypes in the population that have fitness, we are instead interested in the DFE, where fitness effects are measured for all the mutations arising in the background of a particular genotype that has fitness. For now we ignore the influence of evolutionary history on the DFE; we expand on that complication in the following section.

23Examining the DFE over loci for a randomly chosen genotype of fitness can be thought of as sampling the fitness effects from the conditional joint distribution, which generally depends on epistasis. If the number of independent, nonzero epistatic terms is large, then is a multivariate normal distribution defined by the means and covariances of the variables, which in turn can be computed in terms of the ’s from Equation (2). In particular, the conditional means and covariances are,, where is the epistatic VF between loci and and. This implies that the conditional correlation between fitness effects is.

24The DFE simplifies considerably if we make certain additional assumptions on the magnitude of epistatic interactions. If we assume the typical VF is small (i.e., ) and also that is, then correlations are and thus negligible. Then, in a particular sample, we can think of each si as being drawn independently with mean and variance. To compute the DFE,, we first sample the VF from the DVF,, and then sample a Gaussian random variable with the aforementioned mean and variance. This leads to the DFEwhere is the standard normal pdf. Curiously, the correlations between si’s vanish when, in which case the above equation is exact and the DFE is determined entirely by the DVF. Further below, we introduce a specific fitness landscape model for which this relation does hold. Diminishing-returns is naturally incorporated in Equation (9): the mean of is, that is, the DFE shifts progressively towards deleterious values with increasing fitness.

25A key unresolved question is the extent to which evolutionary history influences the DFE and the dynamics of adaptation. That is, what does our theory say about historical contingency?

26Suppose a clonal population of fitness y0 accumulates successive mutations resulting in fitnesses. By virtue of arising on the same ancestral background, the fitness gain of a new mutation,, is in general correlated with the full sequence of past fitnesses and the identity of the mutations through its epistatic interactions with them. Based on these correlations, we use well-known properties of conditional normal distributions to writewhere the weights depend on the VF of the new mutation and its epistatic interactions with past mutations. Here, is the normally distributed residual that depends on the initial genotype and the weights (SI). Equation (10) is a generalization to a sequence of mutations of Equation (5), which we can think of as the special case where.

27To gain intuition, it is useful to first analyze Equation (10) when (i.e., to compute the effect of a second mutation conditional on the first). In this case, we show in the SI thatwhere is the fitness effect due to mutation 1. The first term on the right-hand side is the dependence on the fitness of the immediate ancestor, similar to the corresponding term in Equation (5). The second term quantifies the influence of epistasis between loci 1 and 2 on s2. When, dependence on s1 vanishes entirely and s2 depends only on y1. In contrast, if loci 1 and 2 do not interact,, and s2 is, on average, larger if the mutation at 1 is beneficial compared to when it is deleterious. This has an intuitive interpretation: diminishing-returns applies to the overall fitness and the mechanism through which it acts is epistasis. However, if mutations 1 and 2 do not interact, then the increase in fitness corresponding to mutation 1 does not actually reduce the effect of mutation 2 (as expected by diminishing-returns) so the expected effect of mutation 2 is larger. This analysis suggests that during adaptation, since selection favors mutations with stronger fitness effects on the current background, a mutation that interacts less with previous mutations is more likely to be selected.

28To identify the conditions under which history plays a minimal role, we would like to examine when depends only on the current fitness, yk, and is independent of both the past fitnesses and idiosyncratic epistasis. If this were true, then Equation (5) would apply for new mutations that arise through the course of a single evolutionary path (i.e., the fitness effect of a new mutation is ‘memoryless’ and depends only on its VF and the current fitness). Surprisingly, such a condition does exist. We show that this occurs when the magnitude of epistatic interactions between the new mutation and the previous mutations,, satisfies a specific relation:, where is the combined VF of the previous mutations (SI). In general, this condition is not satisfied, implying that there will be historical contingency that can be analyzed using the framework above. Remarkably, it turns out that a fitness landscape model for which the condition is satisfied does exist and arises from certain intuitive assumptions on the organization of biological pathways and cellular processes. This fitness landscape model additionally serves as an example of a landscape where global epistasis can vary substantially across loci. We describe this model below.

29We introduce the ‘connectedness’ model (CN model, for short). In this model, each locus is involved in a fraction of independent ‘pathways,’ where each pathway has epistatic interactions between all loci involved in that pathway (Figure 4a). The probability of an epistatic interaction between three loci is then proportional to since this is the probability that these loci are involved in the same pathway. When the number of loci is large, we show that in this model, and when is small,, where is the average over all loci (SI). The CN model therefore has a specific interpretation: the outsized contribution to the fitness from certain loci (large vi) is due to their involvement in many different complex pathways (large ) and not from an unusually large perturbative effect on a few pathways. The distribution,, across loci determines the DVF.

30Statistical fitness landscapes such as the NK model and the Rough Mt. Fuji model are related to the CN model. Specifically, the CN model is a subclass of the broader class of generalized NK models (see for a review). However, often-studied fitness landscape models have one important difference that distinguishes them and gives qualitatively different dynamics of adaptation (shown further below): in contrast to the CN model, classical fitness landscapes are typically ‘regular.’ That is, the VF of every locus is assumed to be the same (except the star neighborhood model, which has a bimodal DVF; ).

31The CN model is equivalent to a Gaussian fitness landscape with exponentially decaying correlations (SI). The CN model has tunable ruggedness, where the landscape transitions from additivity to maximal epistasis with increasing. Maximal epistasis corresponds to (and hence ) for all. From Equation (5), this implies that the new fitness after a mutation occurs is independent of the previous fitness, consistent with the expectation from a House-of-Cards (HoC) model (where genotypes have uncorrelated fitness). Regular fitness landscape models with exponentially decaying correlations have memoryless fitness effects under the restrictive assumption that every locus is equivalent. We show that the dynamics of adaptation of the more general CN model are also memoryless, that is, the condition detailed in the previous section holds true (SI). Yet, as we show below, the predicted dynamics for the CN model are very different to those from a regular fitness landscape model.

32We emphasize that the well-connectedness assumed for the CN model is not a requirement for Equation (5) to hold. However, how diminishing-returns influences the long-term dynamics of adaptation depends on the specific genetic architecture and the corresponding fitness landscape. Consider, for example, an alternative model of genetic networks organized in a modular structure (Figure 4b). In this model, each locus is part of a single module and interacts epistatically with other loci in that module to determine the fitness of that module; overall fitness is then determined as a function of the module fitnesses. In this case, the variance contributed by a locus is due to its additive contribution and from epistasis between loci restricted to its module. While the argument for diminishing-returns still applies to the fitness as a whole, it follows from the same argument that diminishing-returns should also apply to each module separately. Consequently, the dynamics of adaptation for the modular model are different from the CN model. For simplicity, we analyze the dynamics of adaptation for the CN model and postpone a discussion of how the dynamics differ for different models to subsequent work.

33We now examine the DFE that follows from Equation (5) and what that implies for long-term adaptation under the conditions for memoryless fitness effects. We henceforth assume a large number of loci with sparse epistasis (though the total number of nonzero epistatic terms is still large). This implies that, and; for simplicity, we also assume strong-selection-weak-mutation (SSWM) selection dynamics and, where are fitness effects and is the population size. Under these conditions, a mutation sweeps and fixes in a population before another one arises. The probability of fixation of a beneficial mutation,, is then proportional to its fitness effect.

34It is convenient to rescale fitnesses based on the total variance in fitness across all possible genotypes by defining. Note that is normally distributed with zero mean and unit variance. Here, has an intuitive interpretation as the ‘adaptedness’ of the organism. When the organism is neutrally adapted, positive and negative epistatic contributions to the fitness are balanced and diminishing-returns is negligible. Diminishing-returns is relevant when the organism is well-adapted. Below, we give the intuition behind our analysis, which is presented in full detail in the SI.

35In the neutrally adapted regime, the linear negative feedback in Equation (5) is negligible and the DFE is determined by the distribution of. Loci with large can lead to a DFE with a long tail. If is the typical VF of a locus, the fitness increases as, where ns is the number of substitutions. Since is a measure of overall epistasis, this implies that epistasis speeds adaptation in the neutrally adapted regime by allowing access to more influential beneficial mutations.

36Fitness increases until the effect of the negative feedback cannot be neglected. From Equation (5), this happens when (i.e., when ). Intuitively, fitness begins to plateau when its accumulated benefit from substitutions is comparable to the scale of the total genetic variance and further improvements are due to rare positive fluctuations. In this well-adapted regime, diminishing-returns and increasing-costs epistasis strongly constrain the availability of beneficial mutations, whose effects can be quantified in this model: for a mutation to have a fitness effect, we require from Equation (5) that, which has probability. Beneficial effects of large arise when has a large positive deviation. The most likely that leads to a particular is when is smallest (i.e., at ), in which case, yielding a tail probability. Remarkably, the beneficial DFE in the well-adapted regime is quite generally an exponential distribution independent of the precise form of the DVF (unless it is singular). In particular, we show in the SI that for the DFE,,which depends solely on the adaptedness of the organism. The exponential form arises because of the Gaussianity of, but the argument can be easily extended to with non-Gaussian tails.

37An exponential beneficial DFE has been previously proposed by but arises here due to a qualitatively different argument. Orr’s result instead follows from extreme value theory: suppose the fitness effects of loci are sampled from a DFE and. Then, the probability that a beneficial mutation has at least a certain effect size is, where the latter approximation holds when beneficial mutations are rare (i.e., is small). A well-known result from extreme value theory implies that for a large family of distributions and for, we have (for some constant ) and therefore. This argument is consistent with our results, but does not yield the dependence of on adaptedness and the rate of beneficial mutations without additional information about.

38Under SSWM assumptions, from Equation (12), the typical effect size of a fixed mutation is, which typically has a VF,

39The above relation makes precise the effects of increasing-costs epistasis on adaptation. As adaptation proceeds, the delicate balance of high fitness configurations constrains fixed beneficial mutations to have moderate VFs. A mutation of small VF is likely to confer small benefit and is lost to genetic drift, while one with a large VF is more likely to disrupt an established high fitness configuration.

40This intuition is not captured in regular fitness landscape models, which assume statistically equivalent loci, that is, for all and is singular. From Equation (9), we see that this leads to a Gaussian DFE whose mean decreases linearly with increasing fitness, in contrast to the exponential DFE in our theory. The key difference is the lack of loci with intermediate effect, which drive adaptation in the well-adapted regime. As a consequence, the rate of beneficial mutations declines exponentially and the fitness thus sharply plateaus at. In contrast, our theory predicts a much slower depletion of beneficial mutations, (SI). The rate of adaptation is (since ), which leads to a slow but steady power-law gain in fitness,. The rate of fixation of beneficial mutations is, which gives.

41We verify our analytical results using numerics. As before, we generated a genotype-phenotype map using the CN model with an exponential DVF, and loci. The DFE can be calculated exactly by plugging in this in Equation (9):

42We simulated the evolution of randomly generated genotypes from to and, and the DFE across all loci was measured (we chose so that ). The theoretical prediction for the DFE, Equation (14), closely aligns with the numerical results (Figure 4c).

43Due to computational constraints, it is difficult to simulate evolution deep into the well-adapted regime. To compute the shape of adaptive trajectories and their variability, we instead simulated SSWM dynamics using the DFE directly from Equation (14), beginning from a neutrally adapted fitness. Typical trajectories (Figure 4d) show rapid adaptation to the well-adapted regime beyond which the fitness grows slowly as, as predicted from theory. The predictions for the number of fixed beneficial mutation are also recapitulated (Figure 4e).

44Recent empirical studies have observed consistent patterns of diminishing-returns and increasing-costs epistasis. Our model gives a simple explanation for these observations. In particular, we showed that these patterns are generic consequences of widespread microscopic epistatic interactions. The intuition underlying this result is that a random mutation typically has a larger disruptive effect on the delicate balance of microscopic epistasis that underpins a fitter background. Our model predicts a quantitative relationship between the magnitudes of global epistasis (i.e., the negative slope of diminishing-returns and increasing-costs epistasis) and microscopic epistasis, which we confirmed using existing data (Figure 3).

45A similar explanation for diminishing-returns and increasing-costs epistasis has been recently proposed by. While our core argument for diminishing-returns and increasing-costs epistasis is the same as in that work, our Fourier analysis framework dissects the features of the fitness landscape necessary to observe these phenomena in terms of experimentally measurable average effects (i.e., the ’s in Equation (2)). In particular, we show that the additivity of a locus’ interacting partners critically determines whether the trend is negative or unbiased. In addition, the Fourier analysis framework yields predictions for the DFEs, the historical influence of past mutations on the fitness effect of a newly mutated site, and motivates the proposed ‘connectedness’ fitness landscape model. The analysis of experimental data presented in Lyons et al. complements the experimental data considered here, lending further empirical support for the prevalence of epistasis and its importance in determining long-term adaptability.

46Our model leads to other experimentally testable predictions. The most direct and accessible test of the theory is to measure the fitness for all possible combinations of mutations at 10–15 significant loci and compare (using Equation (6)) the magnitude of global epistasis to the measured fitness coefficients (the ’s). Additionally, we predict that the magnitude of global epistasis of a double mutant should be nearly the sum of magnitudes of the corresponding single mutants, and any deviations should be biased towards sub-additivity. Since the predictions involve measuring residual variance, experimental noise can be an important confounding factor.

47The observation that diminishing-returns occurs as a ‘regression to the mean’ effect on certain fitness landscapes has been noted previously. The theory developed here quantifies precisely when we should expect to observe these patterns. We emphasize that our key result, Equation (5), is a general statistical relation that holds if epistasis is widespread, irrespective of the specific genetic architecture and the corresponding fitness landscape. Weak epistasis with many loci is sufficient to observe noticeable patterns of global epistasis. However, the argument fails if the contribution of a locus is purely additive or when epistasis is limited to one or a handful of other loci. In the latter case, we expect the fitness effect of a mutation to be dominated by the allelic states of its partner loci, and thus take on a few discrete values. A few examples from indeed exhibit this pattern (e.g., cases where the fitness effect of a specific mutation depends primarily on the allelic state at a single other locus).

48We highlight a distinction between global epistasis discussed in this work and another form of global epistasis (also known as ‘nonspecific’ epistasis) typically used in protein evolution to describe nonspecific epistatic interactions due to a nearly additive trait transformed by a nonlinear function. This nonlinear function creates systematic relationships between epistasis terms and breaks the condition of independent epistatic terms required for our arguments to apply. Specific nonlinearities such as an exponential function may indeed lead to a negative linear trend on average, but the structure of the residuals differs from the one in Equation (5) and observed in data.

49A surprising empirical observation is that the negative linear relationship between fitness effect and ancestral fitness characteristic of global epistasis has different slopes for different loci. Our model identifies the negative slope as twice the fraction of variance contributed by a locus to the trait. To explain the wide range of VFs observed in data, we developed the CN model, a framework to think about the organization of cellular processes that can lead to loci of widely varying VFs. In the CN model, loci have a large VF due to their involvement in many different pathways rather than due to a large effect on a single pathway. The CN model can be viewed as a statistical fitness landscape where loci can have a range of VFs, specified by the DVFs. In the special case of every locus having the same VF, the CN model corresponds to a fitness landscape with tunable ruggedness and exponentially decaying correlations.

50Extending our framework to incorporate adaptation, we showed that the DFE depends only on the current fitness, rather than the entire evolutionary history, under the intuitive assumptions behind the CN model. The theory therefore gives a simple explanation for why phenotypic evolution can be predictable, even while the specific mutations that underlie this evolution are highly stochastic.

51Our framework has an implicit notion of ‘adaptedness’ without referencing a Gaussian-shaped phenotypic optimum, often assumed in models of adaptation (e.g., Fisher’s geometric model). Over the course of adaptation, the DFE shifts towards deleterious values, reflecting diminishing-returns, which naturally arises from our basic arguments. For a well-adapted organism, we show that the DFE for beneficial mutations takes on an exponential form and leads to universal adaptive dynamics. While an exponential DFE for beneficial mutations has been proposed previously based on extreme value theory, our result arises due to an entirely different argument: the tail of the beneficial DFE is determined by loci of intermediate size whose disruptive effect due to increasing-costs is small, yet whose effect size is large enough not to be lost due to genetic drift.

52Our theory further predicts declining adaptability, with rapid adaptation in a neutrally adapted regime followed by much slower increases in fitness, resulting in power-law adaptive trajectories when the organism is well-adapted. This is consistent with observations from the E. coli LTEE. Our model predicts a quicker decline in the number of substitutions compared to the near linear trend observed in the LTEE data. However, the dynamics of fixation in the LTEE deviate strongly from SSWM assumptions. This may explain the discrepancy, although we note that existing theory has only analyzed the effects of clonal interference and other breakdowns in SSWM assumptions for a constant DFE and weak epistasis. Further work will be required to understand how these effects interact with global epistasis. For example, we may expect that the effect of a highly beneficial mutation at a segregating locus is more likely to be attenuated due to interference from subsequent deleterious mutations, while a less-fit lineage has a larger pool of beneficial mutations and is thus more likely to ‘leapfrog’ over more-fit lineages.

53The code and data to generate the figures are available at.

54We use a fitness landscape model with loci to generate the genotype-fitness map. Each locus is assigned a sparsity µ from, which is an exponential distribution with mean. Each of independent pathways sample loci with each locus having probability of being selected to a pathway. We choose so that ensures significant epistasis. All loci in a pathway interact with each other, where additive and higher-order coefficient terms of all orders were drawn independently from a standard normal distribution. The total fitness is the sum of contributions from the pathways. We normalize the coefficients so that the sum of squares of all coefficients is 1, that is, the total variance across genotypes is 1. The mean,, is close to zero from our sampling procedure. The above procedure is a simple and efficient way to generate epistatic terms to order 20, beyond which the computational requirements are limited by the exponentially increasing demand. Note that the effects described in the paper were also observed with only pairwise and cubic epistatic terms.

55The VFs shown in Figure 2a can be calculated numerically from the definition. From the theory, given our choice of, these should follow an exponential distribution with mean. There may be deviations since is finite whereas the calculations assume. To generate Figure 2b, in order to get a range of background fitnesses, we first sample 128 random genotypes. These have fitnesses close to zero; in order to obtain a range of fitness values, we simulated the evolution of these 128 genotypes up to under SSWM assumptions to get genotypes at roughly five fitness values. The fitness effect of applying a mutation (i.e., flipping its sign) is measured for 30 randomly chosen loci (which are kept fixed) over each of the 640 genotypes. This is shown for 5 of the 30 and for the mean over the 30 loci in Figure 2b.

56To generate Figure 2c, we sampled 64 random genotypes and 12 replicates of each. The evolution of these 768 genotypes was simulated for a total of 50 generations with a mutation rate of 1 per generation. The mean fitness gain over the 12 replicates is plotted for each of the 64 founders against their initial fitness.

57To generate Figure 2d, the residuals are measured using the same procedure as for the experimental data analysis described below for the initial 128 genotypes at and the 30 loci with the largest VF.

58Double mutants were created by mutating all pairs of the 30 randomly chosen loci on the 640 evolved genotypes. Their mean fitness effect was computed and plotted along with the mean fitness effect for single mutants, shown in Figure 2e. The VF of the pair of loci for the double mutant was estimated as before and compared to the sum of the estimated VFs of the corresponding single mutants. This is shown in Figure 2f.

59To generate the plots in Figure 4c, we simulated the evolution of 128 randomly sampled genotypes to and. The fitness effect of 200 randomly sampled loci was measured and the distribution is plotted.

60The data from consists of the fitness after the addition of 91 insertion mutations on each of 145 background genotypes. The fitness of a particular mutation at locus can be modeled aswhere are the mutant and background fitnesses, respectively, are constants for each locus, and the residual depends on the background genotype.

61We estimate the VF, where the Pearson correlation, where the symbol ⊕ denotes that the mutant and background fitness datasets are concatenated. is estimated as the negative one-half of the slope of the best linear fit of and. The residuals for each of the 145 genotypes for each of the 91 mutations are simplywhere the overline represents an average over the 145 genotypes, which is used as an estimate of the constant term and. In Figure 3b, we plot the distribution of estimated vi and. In Figure 3c, we compute the AoIL for each locus using Equation (7), which we show in the SI to be. In Figure 3d, we compute the additivity using Equation (8). The additive effect is, and gives the sum of squares of the epistatic terms (SI). In Figure 3e, we compute the variance of the residuals across the 145 genotypes for each locus and plot it against the locus’ estimated. In Figure 3f, we plot the distribution of residuals over all genotypes and loci. The prediction is that in the WE limit the distribution of residuals is determined by, where is a Gaussian random variable. We multiply for each locus with 10,000 i.i.d. standard normal random variables, pool the resulting numbers for all loci, and plot the predicted distribution in Figure 3f. The distributions are variance-matched. While Figure 3e shows that the variance of the residuals aligns with the theoretical prediction of being proportional to slope, Figure 3f shows that the data is also consistent with the predicted Gaussianity of the background-genotype-dependent contribution.