Populations exposed to a high mutation rate harbor abundant deleterious genetic variation, leading to depressed mean fitness. This reduction in mean fitness presents an opportunity for selection to restore fitness through the evolution of mutational robustness. In extreme cases, selection for mutational robustness can lead to flat genotypes (with low fitness but high robustness) outcompeting fit genotypes (with high fitness but low robustness)—a phenomenon known as survival of the flattest. While this effect was previously explored using the digital evolution system Avida, a complete analysis of the local fitness landscapes of fit and flat genotypes has been lacking, leading to uncertainty about the genetic basis of the survival-of-the-flattest effect. Here, we repeated the survival-of-the-flattest study and analyzed the mutational neighborhoods of fit and flat genotypes. We found that the flat genotypes, compared to the fit genotypes, had a reduced likelihood of deleterious mutations as well as an increased likelihood of neutral and, surprisingly, of lethal mutations. This trend holds for mutants one to four substitutions away from the wild-type sequence. We also found that flat genotypes have, on average, no epistasis between mutations, while fit genotypes have, on average, positive epistasis. Our results demonstrate that the genetic causes of mutational robustness on complex fitness landscapes are multifaceted. While the traditional idea of the survival of the flattest emphasized the evolution of increased neutrality, others have argued for increased mutational sensitivity in response to strong mutational loads. Our results show that both increased neutrality and increased lethality can lead to the evolution of mutational robustness. Furthermore, strong negative epistasis is not required for mutational sensitivity to lead to mutational robustness. Overall, these results suggest that mutational robustness is achieved by minimizing heritable deleterious variation.
All well-adapted populations receive an influx of deleterious variation every generation due to de novo mutations. This influx reduces the population's average fitness and causes a mutational load [2, 11]. For populations with high mutation rates, deleterious variation can alter the dynamics of natural selection and lead to selection acting not just on individual genotypes, but on genotypes and their likely mutant offspring [5, 6, 31, 42]. This combination of a wild-type sequence and its mutants is known as a quasispecies and is thought to be a unit of selection in organisms such as the first RNA replicators [14, 15] and modern-day RNA viruses .
Since the quasispecies concept was first proposed [14, 15], many studies have analyzed the conditions and consequences of evolutionary dynamics in the strong mutational regime (e.g., [4, 23, 42]). One consequence of evolution under high mutation rates is that populations are expected to evolve mutational robustness by fixing mutations that decrease the likelihood (and/or decrease the average effect) of deleterious mutations [12, 44]. In the extreme case, this trend can lead to genotypes with lower fitness outcompeting genotypes with higher fitness if the lower-fitness genotypes are more mutationally robust . This effect has been termed the survival-of-the-flattest effect, where flat genotypes, that is, those with lower fitness but greater mutational robustness, outcompete fit genotypes, meaning those with higher fitness but decreased mutational robustness . This effect was first observed in the digital evolution system Avida; it was later demonstrated with experimental evolution of RNA viruses [8, 28, 36].
Even though the evolutionary dynamics of the survival-of-the-flattest effect are well understood, the exact genetic causes of robustness are unknown. Indeed, the differences in the distribution of fitness effects  (that is, the local fitness landscapes) between fit genotypes and flat genotypes have not been explored. While it is generally assumed that mutational robustness evolves due to the evolution of an increased number of neutral mutational neighbors [18, 39, 41], other authors have argued that increased mutational sensitivity may evolve [3, 23, 26, 32, 46]. Additionally, recent work has indicated that populations evolving under strong genetic drift are expected to evolve drift-robust genetic architectures, giving rise to genotypes with a decreased likelihood of small-effect deleterious mutations and an increased likelihood of neutral and/or strongly deleterious mutations . A mapping of the local fitness landscapes of fit genotypes and flat genotypes is needed to compare the survival-of-the-flattest effect with these other proposed evolutionary mechanisms.
To explore the local fitness landscapes of fit genotypes and flat genotypes, we first repeated the original survival-of-the-flattest experiment. We then took genotypes adapted to either low or high mutation rates and obtained the distribution of fitness effects for these genotypes by examining mutants up to four substitutions away from the wild-type sequence. We found that flat genotypes had a greater likelihood of beneficial and neutral mutations (as expected for mutationally robust genotypes), while fit genotypes had a greater likelihood of deleterious mutations, as expected from mutationally sensitive genotypes. Contrary to expectations, we also found that flat genotypes had a greater likelihood of lethal mutations, suggesting that the survival-of-the-flattest effect relies on a reduction of heritable deleterious mutation, which can be achieved both by an increase in neutral variations and by an increase of lethal mutations. These results illustrate the complexity of mutationally robust genome architectures on nontrivial fitness landscapes.
2 Materials and Methods
Here we describe the relevant details of Avida (version 2.14; our version, with added code for data analysis, is available at https://github.com/joshf93/MappingThePeaks) for the current experiments (see  for a full overview of the Avida software). In Avida, a finite population of computer programs (Avidians) compete for the resources (memory, space, and processing time) required for reproduction. Each program consists of a circular genome of computer instructions that encode the ability to self-replicate. During this replication process, instructions may be copied inaccurately, resulting in mutations that are passed onto an Avidian's offspring. As Avidians with faster replication speeds enter the population, they reproduce faster than slower-replicating Avidians, leading to the spread of faster replicators and the extinction of slower replicators. In other words, there is selection for faster replicators. Therefore, because there is heritable variation and selection in Avida, populations of Avidians undergo evolution by natural selection . Avida has been used to test many evolutionary concepts difficult to test in biological systems [1, 10, 19, 20, 24, 29, 38].
The Avida world consists of a toroidal grid of N cells; each cell can be occupied by at most one Avidian, and thus N is the maximum population size. During reproduction, each offspring Avidian is placed in one of the nine adjacent cells of its parent, including the cell that contains the parent. If some of these cells are unoccupied, a random empty cell is chosen. If all neighboring cells are occupied, a random cell is chosen, its occupant is removed, and the new Avidian then occupies the cell. This random selection of replacement adds an element of genetic drift to Avida.
The limited number of cells in the Avida environment is one of the limiting resources over which Avidians compete. The other limiting resource is the opportunity to execute the instructions in an Avidian's genome. The resource required to execute one instruction is called a single instruction processing unit (SIP). During each unit of time in Avida (called an update), 30N SIPs are available to the population for execution. These SIPs are probabilistically distributed to Avidians according to a figure of merit that is related to their ability to perform certain computations. These computations (binary logic operation) are the traits that Avidians display. In a monoclonal population, all Avidians express the same traits and thus have the same merit and thus receive on average 30 SIPs (i.e., they execute 30 genome instructions during that update). However, Avidian lineages can evolve the ability to increase merit, and thus increase the number of instructions that can be executed per update.
Fitness in Avida is implicit and can only be estimated by executing the computer instructions within an Avidian's genome. In this regard, Avida differs from classical evolutionary simulations where genotypes are assigned a fitness and their frequency increases or decreases proportionally to this assigned value. Here, faster replicators outcompete slower replicators and fitness emerges from this dynamic. There are two means by which an Avidian lineage can evolve increased replication speed. First, genotypes that require fewer instruction executions to undergo reproduction will reproduce faster than Avidians that require a greater number of instruction executions. In environments where genome size can vary, this is often achieved by decreasing the genome size. In fixed-genome-size environments, this trajectory of fitness gain still occurs, often through the fixation of instructions that copy more genome information from parent to offspring within a single execution of the copy loop. The second method of fitness increase in Avida is through the evolution of complex adaptations: the ability to perform one- and two-input Boolean logic calculations (the traits discussed earlier). In the setup used here, nine of these functions are potentially rewarded with SIPs; the number of SIPs awarded is proportional to the complexity of the evolved trait .1 In other words, the more calculations an Avidian can perform, the more SIPs it receives per update, and the larger the portion of its genome it can then execute per unit time. This increases the Avidian's replication speed and leads to selection for the ability to perform these calculations. When Avida estimates a genotype's fitness, it is estimated as the genotype's merit divided by the number of instruction executions required for replication (otherwise known as its gestation time).
Avidian populations evolve these complex adaptations by fixing a sequence of mutations. This genetic variation is introduced into a population when an Avidian copies one instruction inaccurately into its daughter's genome. The rate at which these mutations occur is set by the experimenter and is usually set as a rate per instruction copied (a copy-mutation rate). When mutations occur during this copying process, multiple mutations can occur per replication event. Additionally, in experiments where genome size can evolve, insertion and deletion mutations occur upon division (not during the replication process) at a prespecified rate. These indel mutations insert or delete one instruction. Finally, large-scale genome changes can occur, but these changes are caused by the specific replication algorithm encoded by an Avidian's genome; these mutations are inherent in the Avida genetic code (implicit mutations) and are not under the control of the experimenter.
2.2 Experimental Design
Here, we repeated the experimental protocol of the original survival-of-the-flattest study . We differed from the original study in one aspect: the placement of offspring. Here, we used Avida's current default setting, where a new offspring randomly replaces one of the nine neighbors of its parent (including possibly the parent). The original article used mass-action reproduction, where new offspring could randomly replace any offspring in the population. We first generated a collection of genotypes during the initial adaptation step. We evolved 201 populations of 3,600 individuals for 5 × 104 updates.2
All populations initially started with the default Avida ancestor with a 100-instruction genome. Point mutations occurred at a rate of 0.0075 mutations per instruction copied, and insertion/deletion mutations occurred at a rate of 0.05 mutations per division each. At the end of these populations' evolution, we extracted the most abundant genotype.
For the mutation-rate adaptation phase of the experiment, we duplicated each of the 201 genomes and used one as progenitor of an experiment at a fixed low genomic mutation rate of 0.5 mutations per genome per generation, and one at a fixed high genomic rate of 2.0 mutations per genome per generation. To achieve this, we adjusted the per-site mutation rate of each genome so that the targeted genomic rate was achieved, and disallowed sequence length changes. In this manner, the initial adaptation created genetic variation used as input to the adaptation step at a fixed genomic rate. From this stage onward each genotype was limited to the tasks it could complete at the end of the initial adaptation phase—that is, novel traits were not rewarded (and as a consequence could not become fixed in populations of this size). However, they could lose traits they already possessed, and could re-evolve them subsequently. Each population adapting to a high or low mutation rate consisted of 3,600 individuals and evolved for 104 generations. After evolution we again selected the most abundant genotype from each population and classified genotypes from the low-mutation environment as fit genotypes, and genotypes from the high-mutation environment as flat genotypes.
Next, we performed competition experiments between each pair of fit and flat genotypes (each pair has the same ancestor). We only performed these competitions with genotype pairs where the fit genotype was 1.5 times fitter than the flat genotype. For the competitions, we seeded populations of 3,600 individuals with 1,800 individuals of each genotype and ran the simulation for 200 generations. We repeated these competitions across a range of genomic mutation rates (0.5 to 3.0 mutations per genome per generation in increments of 0.5) and performed five replicates per mutation rate. We tracked the abundance of the fit genotypes over time to determine the outcome of the competition. We classified the flat genotype as the winner if it reached at least 95% of the population in three of five replicates at a given mutation rate.
2.3 Data Analysis
We performed mutational analyses to determine the fitness landscapes of the fit and flat genotypes used in the competitions, using Avida's analyze mode. In analyze mode, Avidians run through their life cycle in isolation (rather than in a population), and characteristics such as their fitness can be estimated. We first estimated the base fitness of all Avidian genotypes used in the experiment. Then, we took the genotypes used in the competition and generated all possible single mutants. For each genotype, there are 25L single mutants, where L is the genome size and 25 is the number of possible alternative alleles at each position in the genome, as there are 26 instructions in the Avida instruction set. We then repeated this procedure for all double mutants. For triple and quadruple mutants, we randomly sampled 108 mutants, as the number of all possible mutants made generating every mutant computationally prohibitive.3
We first evolved 201 populations and isolated the most abundant genotype from each population. We then evolved these genotypes in both a high-mutation-rate environment and a low-mutation-rate environment and isolated the most abundant genotype from each environment. The genotypes from the low-mutation-rate environment did not significantly change in fitness compared to their ancestor (Figure 1(a); Wilcoxon rank sum test, p = 0.411). The genotypes from the high mutation rate environment significantly decreased in fitness compared to their ancestors (Figure 1(a); p = 2.216 × 10−7). By tracking the fitness of the high-mutation-rate populations through time, we observe that these genotypes first severely decrease in fitness upon transfer, but then recover fitness as they search out new, and presumably more mutationally robust, fitness peaks (Figure 1(b)).
Of our 201 pairs of genotypes adapted to low mutation rates and high mutation rates, 166 pairs evolved so that the low mutation rate genotype (the fit genotype) had a fitness greater than 1.5 times the high-mutation-rate genotype (the flat genotype), which is the threshold used in the original study . We took these pairs and performed competition experiments between the two genotypes across a range of mutation rates. In 76 (45.8%) of these competition experiments we noted that the eventual winner of the competition depended on the mutation rate, and that a switch between winners occurred at a critical mutation rate, just as was observed in the initial study . Fit genotypes win the competition at low mutation rates, but flat genotypes win the competition at high mutation rates even with the large fitness difference (Figure 2). It is worth noting that our definition of winning excluded many genotype pairs that showed the effect to a lesser degree.
Next, we analyzed the mutational neighborhood (the local fitness landscape) of all 166 pairs of fit and flat genotypes. We first generated all genotypes containing one point mutation and measured their fitness. There are significant differences in particular classes of mutations. Flat genotypes have a greater likelihood of beneficial (Figure 3(a); Wilcoxon rank sum test, p = 1.35 × 10−20), neutral (Figure 3(b); p = 2.66 × 10−35), and lethal (Figure 3(d); p = 0.0168) mutations. Fit organisms have a greater likelihood of deleterious mutations (Figure 3(c); p = 1.31 × 10−44). Additionally, fit genotypes have lower mean relative fitness than flat genotypes (Figure 3(e); p = 2.24 × 10−27).
Furthermore, we analyzed the mutational neighborhood of these genotypes for sequences that are two, three, and four mutations away from the wild-type sequence in order to estimate differences in epistasis between fit and flat genotypes. We found that the fraction of lethal mutants increased with mutational distance for both the fit and the flat genotypes, with the flat genotypes always showing significantly more lethals than the fit types (see Figure 4). The fraction of deleterious mutations decreases with increasing mutational distance, with the fit types showing a higher fraction of deleterious mutations than the flat types for all distances we tested.
As expected from the data on the relative fitness of one-step mutants (Figure 3(e)), fit genotypes have greater α values (Figure 5(a); p = 1.06 × 10−27), that is, they are less robust to mutations. Flat genotypes have greater β values (Figure 5(b); p = 5.44 × 10−34), indicating that flat genotypes have somewhat negative epistasis (β > 1) compared to the fit genotypes, whose average β is solidly in the positive epistasis territory (β < 1). However, it should be noted that flat genotypes have, on average, almost no epistasis (β very close to 1). These data demonstrate how the selection pressure for mutational robustness drives populations to alternative areas of the fitness landscape.
We used digital experimental evolution to map the local fitness landscape of high-fitness, but mutationally fragile, fit genotypes adapted under low mutation rates, and low-fitness, but mutationally robust, flat genotypes adapted under high mutation rates. After repeating the original survival-of-the-flattest effect , we calculated the distribution of fitness effects for de novo point mutations up to four substitutions from the wild-type genotypes. We found that, as expected, flat genotypes were more mutationally robust than fit genotypes. Flat genotypes had a lower likelihood of deleterious mutations and a greater likelihood of neutral and beneficial mutations than fit genotypes. Surprisingly, flat genotypes also had a greater likelihood of lethal mutations and more negative epistasis than fit genotypes, suggesting that the factors responsible for mutational robustness involve only limiting heritable deleterious variation.
Around the time when the original survival-of-the-flattest article first appeared, mutational robustness was largely conceptualized as a decrease in the likelihood of a deleterious mutation [39, 45]. And while these changes in genetic architecture can lead to mutational robustness, work since that original publication has suggested that it is increases in the effect of deleterious mutations, not decreases, that lead to robustness [3, 32]. Those authors argue that increased severity of deleterious mutations can lead to population-level robustness because purifying selection is more effective against strongly deleterious mutations. It has also been proposed that increased negative (or synergistic) epistasis can evolve as a mechanism to increase the strength of purifying selection and hence robustness [26, 46].
Our results suggest that all of these genetic mechanisms (increased neutrality, increased severity of deleterious mutations, and increased negative epistasis) can contribute to mutational robustness (although in this study the “increased negative epistasis” really is more of an absence of positive epistasis). These data support a previously proposed relationship between increased neutrality (low α) and negative epistasis (β > 1), based on data from Avida, RNA folding, and neutral-network-fitness landscapes . In other words, a population can only maximize immediate mutational robustness at the cost of increased deleterious severity of mutants farther away in the fitness landscape. We should also note that mutations on a flat genetic background are not, on average, negatively epistatic. Flat genotypes evolve to have, on average, no epistasis between mutations, as previously shown to occur in high-mutation-rate Avida populations . However, they are more negatively epistatic than fit genotypes, which have positive epistasis between mutations. The lack of negative epistatic interactions in flat genotypes is likely related to the increased likelihood of lethal mutations. Additional mutations on a genetic background with a lethal mutation cannot, by definition, lead to a greater decline in fitness (negative epistasis).
Theoretical arguments suggest that increasing the severity of mutations not only protects against high mutation rates but also protects against fitness loss due to genetic drift , although this depends on the specific distribution of fitness effects . Recently, LaBar and Adami proposed the concept of drift robustness: the idea that small populations tend to populate drift-robust fitness peaks with a low likelihood of slightly deleterious mutations (and a high likelihood of both neutral and strongly deleterious mutations), while large populations would adapt to drift-fragile fitness peaks with a high likelihood of slightly deleterious mutations . Both flat genotypes and drift-robust genotypes have some noticeable similarities in their distribution of fitness effects (increased neutrality, in particular), so it is worth asking whether drift robustness and survival of the flattest are really the same phenomenon.
The main difference between drift robustness and survival of the flattest is the role of natural selection in the two phenomena. The evolution of drift robustness does not require competition between drift-robust and drift-fragile genotypes . Instead, small populations tend to evolve towards drift robustness because they can only maintain fitness on drift-robust genetic backgrounds. Survival of the flattest, however, occurs when a flat genotype outcompetes a fit genotype, and is driven by selection for mutational robustness. It has also been established in survival of the flattest that the critical high mutation rate above which this selective advantage occurs is independent of population size , while drift robustness (almost by definition) is strongly dependent on population size. Furthermore, our previous work demonstrated that drift-robust genotypes are not more mutationally robust than drift-fragile genotypes, if one measures mutational robustness as the mean relative fitness of all single mutants . Single mutants of flat genotypes are more fit, on average, than those of fit genotypes, again illustrating a difference between drift robustness and survival of the flattest. Thus, while similarities exist between drift-robust genetic architecture and flat genetic architecture, it appears that they represent distinct evolutionary responses to different population-genetic environments.
It would be worthwhile to determine if the results observed to hold here can be replicated with RNA viruses, the biological model for the survival-of-the-flattest effect [8, 28, 36]. It would be illuminating to test if the likelihood of lethal mutations does increase in these evolved viruses, especially in light of studies arguing that RNA viruses have high mutational fragility, not high mutational robustness . Such a study would lead to a better understanding of the mutational mechanisms that lead to robust organisms in nature.
In high-mutation-rate environments, populations will evolve alternative genetic architectures that increase their robustness to deleterious mutations. Classically, this is viewed as a survival-of-the-flattest effect , where populations adapt to fitness peaks with increased neutrality . Here, we repeated the original experiments to finely map the fitness landscapes of fit and flat genotypes. Flat, mutationally robust genotypes adapt to fitness peaks with an increased likelihood of both neutral genotypes and lethal genotypes, a decreased likelihood of deleterious mutations, and greater negative epistasis compared to fit genotypes. These results demonstrate the genetic complexity of mutationally robust genomes that minimize the impact of heritable deleterious genetic variation.
T.L. acknowledges support from a Michigan State University Distinguished Fellowship, a BEACON fellowship, and the Russell B. DuVall award. This work was supported in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research. This material is based in part upon work supported by the National Science Foundation under Cooperative Agreement No. DBI-0939454. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.
The performance of a calculation increases an Avidian's merit as a power of the number of operations needed to perform the calculation, with increasing powers for calculations of increasing difficulty. These powers are the default Avida setting. The exact values appear in the config file that is available on the github repository.
The number of possible n-mutants of a length L genome encoded in an alphabet of size D is (D − 1)n, which translates to over 1.5 × 1012 possible mutants for L = 100, n = 4, and D = 26.
List of Abbreviations
RNA: Ribonucleic acids
SIP: Single Instruction Processing
Availability of Data and Material
All Avida configuration files and data analysis scripts to recreate these experiments, analyze the data, and recreate the figures are available at https://github.com/joshf93/MappingThePeaks.
The authors declare that they have no competing interests.
J.F. performed the experiments and analyzed the data. J.F., T.L., and C.A. designed the study and wrote the manuscript.