According to a central principle of molecular evolution, the likelihood that a given mutation occurs is independent of the mutation's phenotypic consequences. Organisms cannot choose specific mutations. This tenet was challenged by [1], who observed that under a certain selective pressure, E. coli cells appeared to acquire an excess of beneficial mutations. The idea that cells can somehow 'direct' evolution was thought provoking, and stimulated many investigations (for reviews see [2–6]). While the notion that cells can directly decide in which genomic regions to increase their mutation rate has been mostly abandoned [4, 7], the original observations by [1] have been corroborated (see above reviews).

If mutations arise independently of their phenotypic consequences, then how can adaptations occur that require multiple amino acid mutations and for which the intermediate stages are either selectively neutral or disadvantageous? Large populations can climb multiple fitness peaks, even with disadvantageous intermediate alleles [8, 9]. Although no new mechanisms are therefore required to explain the evolution of complex proteins [10], we propose that errors in transcription and translation (phenotypic mutations) allow the selection of the intermediate mutations of a multiple-mutation requiring trait, and can thus speed up the evolution of complex traits.

Studies on the phenotypic mutation rate indicate that it is orders of magnitude larger than the genotypic mutation rate [11, 12]: the misreading error rate during protein synthesis is estimated to be between 10^-3 to 10^-4 misreadings per codon [13], compared with a genotypic mutation rate of between ~10^-7 to 10^-11 [14]. Consequently, for a protein of 300 residues, on average more than 1 in 10 copies of the protein will contain a mutation. Using mutation rates derived from the literature and conservative biological assumptions, we show via mathematical modeling and simulations that phenotypic mutations allow evolution to select for neutral intermediate alleles of a multi-mutation trait, actually selecting for proteins whose exact DNA sequence is not in the organism under selection. Evolution is then able to look ahead for evolutionary jackpots in sequence space.

Our theory is based on the following hypothetical scenario. A protein can increase the fitness of an individual if it evolved a specific trait. This trait requires two mutations, for example a disulphide bridge between two cysteine residues. A modification of only two residues can result in large structural changes [15]. Having only one of the required mutations is either selectively neutral or deleterious, however when an individual has only one mutation, small amounts of the protein with both mutations will be produced due to phenotypic mutations. If the presence of both mutations at low concentrations provides even a small fitness improvement then the allele with one mutation will spread though the population. As the frequency of the intermediate allele increases, there is a greater probability that if the second mutation occurs, it will be in the presence of the first mutation, and thus provide the full fitness benefit.

Our hypothesis is similar to an effect proposed in 1896 by J.M. Baldwin [16], known now as the "Baldwin Effect" or "Organic Selection". The core idea is that the probability of a trait occuring can be selected for, not just the trait itself. If the phenotypic plasticity of an organism allows it to learn a trait, then during the course of evolution, the organism's descendants may get better at aquiring the trait, and may ultimately acquire genes that code for the trait directly. Using a genetic-algorithm based model, Hinton and Nowlan [17] tested this idea, and concluded that organisms with phenotypic plasticity evolved faster, even though the learned traits did not become hereditary in their study. Subsequent studies futher explored the Baldwin effect in different fitness landscapes in the context of machine learning (see [18] for a review). This work clarified two aspects of the Baldwin effect: that lifetime learning can accelerate evolution in certain contexts, and that this learning usually comes with a cost. Organisms that have to acquire a beneficial trait by learning will generally have a selective disadvantage over organisms that genetically encode the trait. That evolution can select for the ability to learn has been demonstrated experimentally in fruit flies [19]. Our model differs in a subtle but important way from the Baldwin effect. The Baldwin effect describes learning from the perspective of the individual, meaning that an organism starts its life without the trait, but then later has a chance of acquiring the trait. In our model, phenotypic mutations occur at a given rate, although some individuals are more predisposed to the highly beneficial phenotypic mutations than others. The traits are not learned, because individuals that have the neutral, intermediate mutation will always express the beneficial phenotypic mutation at a low rate. The organisms in our model are not exactly phenotypically variable, rather phenotypically diffuse. They do not learn, but rather possess a small part of the many phenotypes close to their genotype. We call this effect the "look-ahead effect" as opposed to the Baldwin effect to highlight that no learning takes place.

Our work is more closely related to works on phenotypic plasticity [20] and random or noisy phenotypes [21]. The aim of this article is to derive explicit analytic expressions for the fixation process of genes whose fitness is modulated by phenotypic mutations, and to show that adaptive phenotypic mutations can undergo positive selection under biologically plausible conditions.

We model the scenario of a protein evolving a trait that requires two mutations. The model is based on a population-genetics framework where a single gene can evolve into different alleles. We do not consider duplication and divergence of genes. In addition, the process described here will likely only occur for proteins with sufficiently long half-lives, as the protein must persist for some time to exert a phenotypic effect. As we model only a single gene, we expect our results to be more relevant for single-celled organisms and viruses than for multicellular organisms, which tend to have larger genomes and smaller effective population sizes than microorganisms.

The model consists of the evolution of three non-recombining haploid genotypes, where each genotype contains one of the three alleles shown in Figure 1. The three different alleles are named according to number of relevant mutations, corresponding to zero mutations (allele 0), a single mutation (allele 1), and both mutations (allele 2) required for the adaptive feature. Having both mutations of the adaptive feature provides a selective advantage s. We assume that the intermediate allele (allele 1) is selectively neutral if transcribed and translated without error. We specifically take into consideration errors in transcription and translation, that is, phenotypic mutations.

In the model, the population initially consists of one individual carrying allele 1 and N - 1 individuals carrying allele 0. So long as allele 1 is present, allele 2 can be generated by mutations. The population evolves for a fixed time period, during which allele 2 can be generated by mutation and go to fixation. In each generation, selection increases the frequency of the alleles according to their corresponding fitness values. Allele 0 has a fitness of 1. Allele 2 has a fitness of 1 + s, where s is the selection coefficient provided by the adaptive feature. The fitness of allele 1, the intermediate allele with only a single mutation, depends on the phenotypic error rate. Most phenotypic errors will be neutral or deleterious, however some will be beneficial. For simplicity, we assume that the length of the protein and the expression level are both constant. In addition, we do not explicitly model deleterious phenotypic mutations. As long as the spectrum of deleterious phenotypic mutations does not change substantially among alleles 0, 1, and 2, we can treat its fitness effect as a common factor which we divide out of all fitness values. This assumption becomes invalid if, for example, phenotypic mutations for allele 1 are significantly more deleterious than those for either allele 0 or allele 2. What happens when we relax these assumptions will be the subject of future work.

If there are no phenotypic mutations, allele 1 has the same fitness as allele 0. However, if phenotypic mutations occur, allele 1 can produce a small number of allele 2 proteins due to phenotypic errors. The fitness of allele 1 is therefore dependent on the number of such errors. We assume, for the sake of simplicity, that fitness is a linear sum of individual proteins, meaning that if some phenotypic variants of a protein have a higher fitness, then the overall fitness of that allele is proportionally increased.

We let r be the number of residues that can potentially complement the first mutation to provide the full two-residue adaptive feature. These r residues represent, e.g., the sites at which the second cysteine of a cysteine bridge could arise; other similar two-residue mutations that significantly improve functionality can be proposed. Two residues that comprise an adaptive trait are likely to co-evolve, because if a mutation occurs in one of the residues, selection strongly favors a compensatory mutation in the other. Based on a large data set, [22] found that co-evolving residues are spatially near. Co-evolving residues were, on average, 98.6 amino-acids apart along the sequence, but had a mean spatial distance of 6.9 Å. This spatial distance can be compared to the width of the van der Waals volume of an amino-acid (5–6 Å), showing that most co-evolving residues are effectively in contact proximity. Therefore, r is mostly independent of the size of the protein, as long as the protein is of sufficient length. [23] calculated the mean contact density (the mean number of residues in contact with a given residue) for 194 yeast proteins, and found that most residues have a mean contact density of seven to eight residues. In this work we use r = 8. Given r possible positions for the second residue, and assuming that each position requires a specific residue, the fraction of proteins of allele 1 containing the second (now highly beneficial) mutation is $\beta =\frac{r}{19}\lambda$, where λ is the per codon non-synonymous phenotypic mutation rate. In this model, we use λ = 4.5 × 10^-4 mistranslations per codon [24, 25]. The fraction β of allele-2 proteins contribute to the fitness, giving allele 1 a fitness of 1 + sβ.

When considering genetic (i.e. inherited) mutations, for simplicity we neglect back mutations (e.g. from allele 1 to allele 0), and assume there are no recurrent mutations of allele 1 from allele 0 (the model starts with a single copy of allele 1). Allele 2 arises via a mutation from allele 1. We ignore the possibility of a double mutation directly from allele 0 to allele 2, as this probability is extremely small in the parameter range we are interested in. The genetic mutation rate for allele 1 mutating into allele 2 is derived as follows: For microbes, the rate of mutations per nucleotide per generation is between ~10^-7 to 10^-11 [14]. Here we use 10^-8 as the non-synonymous mutation rate per codon per generation. The resulting mutation rate for changing allele 1 into allele 2 is $U=\frac{r}{19}{10}^{-8}=\frac{8}{19}{10}^{-8}$.

Genes can also acquire null mutations, rendering the gene non-functional and therefore eliminating the organism. The null mutation rate for protein-encoding genes is on the order of 10^-6 per generation [14]. However, this rate will depend on the length (L) of the protein. Assuming an average protein length of 300 residues, the per-residue null mutation rate is given by 10^-6/300 = ~3.3^-9. For a protein of length L, the null mutation rate is given by μ = 3.3^-9L.

We have described a model demonstrating the consequences of positive phenotypic mutations on the evolution of a single gene. We have compared numerical simulations with the analytical approximations and found them to be in good agreement. When phenotypic mutations exert an effect on fitness, selection can operate on the intermediate allele of a complex trait, which otherwise (without phenotypic mutations) would be neutral. We refer to selection for the intermediate allele as the look-ahead effect, because this effect allows evolution to select for sequences not yet in the genome.

The approximation for small sβ, Eq.(6), shows most clearly the relationship between the parameters. The look-ahead effect is proportional to N, s, and β, and sets in when N is on the order of 1/(sβ). For large Nsβ, the look-ahead effect saturates. The asymptotic value of ξ is approximately N for NU ≪ 1.

Therefore, large populations have two advantages over small populations in terms of the look-ahead effect: the effect sets in for smaller values of s, and saturates at a larger asymptotic value ξ. Of course, even in the absence of the look-ahead effect, larger populations can more easily traverse multiple local fitness peaks [9]. Because the selection coefficient s depends on the environment, a valid question is how often does s reach sufficiently high levels so that the look-ahead effect can operate. For microbial species such as bacteria, sufficiently large s should be reasonably common. Many bacteria experience highly fluctuating [27] and structured [28] environments, where growth is limited by the lack of a key trait. An obvious and extreme example is antibiotic resistance. Evolving a defense against an antibiotic molecule can involve only a few amino acids [29], and those individuals that can generate an enzyme capable of degrading the antibiotic, even if briefly or weakly, will have a very large fitness increase. In fact, if the efficacy of the antibiotic is 100% on susceptible genotypes, a mutation providing only moderate resistance has an infinite selective advantage. And even for very small antibiotic concentrations, mutants diffiering by only two amino acids at a single β-lactamase gene can be selected effectively [30, 31]. Thus, bacteria may frequently experience environments in which a large fitness increase (large s) is only a few mutations away. Similarly, in bacteriophages, selective coefficients s of 10 or more are not uncommon, even for individual mutations [32]. Our work is entirely theoretical, but we expect that it will be possible to experimentally verify our predictions in future work. For experimentally observing the look-ahead effect, we would need a system where s and N are both large, while β (the phenotypic mutation rate) can be modified. The values of both N and s used in this work are well within biologically realistic ranges achievable in a microbiological laboratory. Conditions for large s may be created with e.g. antibiotic resistance, which is a common laboratory workhorse. Unfortunately, many antibiotics function by reducing translation fidelity [33], and thus would conflate s and β. Changing β could involve a mutated ribosome. Ribosomes appear to be optimized for accurate and efficient translation of mRNA [34], and several examples of altered ribosome fidelity exist, both increasing [35] and decreasing fidelity [36]. Specific regions of the ribosome rRNA sequence have been identified as influencing fidelity [37], and various agents can reduce fidelity, e.g., streptomycin, magnesium, and ethanol [36]. Few mutations may be sufficient to alter the fidelity of a ribosome, for example, a single mutation in the S5 ribosomal protein in E. coli increases frameshifting and nonsense mutations [38]. In yeast, mutations in the 18S RNA have been found that both increase and decrease translational fidelity [39].

In fact, a similar system to what we propose was already used to estimate the effects of tRNA competition on misreading error rates [13]. Here, several firefly luciferases were constructed with inactivating point mutations at an essential active-site lysine residue. Mistranslations of the transcript occur, rescuing the mutant and restoring the wild-type function. This system demonstrated a possible evolutionary constraint for the system presented in this work, in that some codons are more likely to be misread than others, depending on the relative amounts of tRNAs.

In this work, we have calculated the look-ahead effect from a comparison between the two cases of β > 0 and β = 0. The latter may not be experimentally possible; any experiment will likely compare two different positive values of β. Nevertheless, Figure 5 shows that a larger look-ahead effect can be achieved with a higher β, where increasing β by one order of magnitude both increases the look-ahead effect by an order of magnitude and lowers the smallest s where an effect is observed. Of course, our model does not take into account the loss of fitness or other confounding effects from a higher phenotypic mutation rate. Thus, a balance must be found in having two different values of β that are different enough to measure, while at the same time minimizing the confounding effects. The most obvious consequence of increasing the phenotypic mutation rate is that overall fitness may be reduced, for example in E. coli, where a higher translational error rate activates stress responses [40], or in mouse, where such errors are implicated in neurodegeneration caused by misfolded proteins that aggregate [41]. Increasing translational fidelity may not come without fitness cost either. The hyperaccurate mutations in the 18S RNA in yeast [39] cause an increase in oxidative stress. This observation suggests that cells consume more energy to achieve hyperaccuracy. It may also partially explain why the phenotypic error rate is much higher than the genotypic error rate, as there is possibly a direct disadvantage in reducing the phenotypic error rate, rather than only reducing the selective advantage that occurs if the phenotypic error rate is reduced, as discussed in [42].

Buerger et. al [42] asked whether evolution has selected for the current phenotypic error rate, which does not differ significantly between eukaryotes and prokaryotes [24, 25] even though the source of errors is different. They suggested that the increase in fitness becomes incrementally smaller for improvements to transcription and translation fidelity. We would like to speculate that the phenotypic error rate is on the border between minimal costs (of e.g. misfolded proteins) and maximum payoff (via the look-ahead effect). The goal of our analysis was to demonstrate that the look-ahead effect is theoretically possible, and as such, we intentionally excluded confounding factors for the sake of clarity. There are several aspects not considered in our model that may play important roles. For example, in this work we did not consider the expression level. For low expressed genes, the mutation from allele 1 to allele 2 will occur less frequently compared to highly expressed genes. However, if allele 2 is produced it will be at a higher concentration (of allele 2 mutant proteins in a population of allele 1 proteins), as the overall copy number of allele 1 is low. This difference in expression levels is likely reduced in a large population, where beneficial mutations occur with sufficient frequency. Another factor related to the expression level is translational robustness. It has been proposed that highly expressed genes are under selection to properly fold despite phenotypic mutations, and consequently evolve slower [43, 44]. If a gene is robust to translational errors, then it can tolerate a larger variety of mutations, of which some may be intermediates to a new adaptive multi-residue trait. Thus, translational robustness may increase the sequences available for experimentation at the phenotypic level. However, if the intermediate allele is itself not robust to errors in translation, then it will not be neutral, and may be selected against. The location of the protein trait will also influence the viability of the intermediate allele: mutations near the surface of the protein are less likely to disrupt the protein compared to mutations in the core [45].

In the presence of noise, phenotypic mutations may also help purge negative mutations [46]. If we have a system similar to the one described in this work but the final two-mutation trait is deleterious, then the phenotypic errors will lead to a selective disadvantage of the intermediate genotype. To give a concrete example, consider the case of prions, where an intermediate mutation favouring the formation of prions would be expressed at a small rate and would increase the liklihood of forming the misfolded proteins [47]. Since the majority of mutations are deleterious, the negative look-ahead effect is probably more common than the positive look-ahead effect on which we focused here.

In this work we use a single fitness optimum, and do not take into consideration multiple local optima as done by Borenstein et al. [21], who studied the effect of learning and of noisy phenotypes on evolution. Borenstein et al. considered varying learning rates, and showed that there is a trade-off between the amount of phenotypic plasticity and both the speed of reaching a local optimum and the genetic stability of the evolving individuals. It would be of interest to see if the same conclusions apply using our model, replacing phenotypic plasticity by a diffuse phenotype and learning rate by the phenotypic mutation rate. Such an analysis will require, however, that we explicitly model the deleterious spectrum of phenotypic mutations, and allow for different distributions of phenotypic mutations for allele 0, 1, and 2.

In conclusion, we propose that organisms can experiment with protein sequences that are mutationally close to the current sequence, but not yet in the genome. This effect allows selection for intermediates of complex traits, opening up a more direct route to the trait and thus reducing the time needed for fixation in the population.

The numerical simulations were written in Java using the Colt scientific library [48] for the generation of random numbers. The analytic expressions were evaluated using both Mathematica and Python, the latter in conjunction with the SciPy package [49]. Source code for the numerical simulations is available on request from DJW.

The population in each simulation is represented by three numbers, corresponding to the abundance of each of the three alleles. As described, the initial abundances are N - 1, 1, 0 for alleles 0, 1, 2, respectively. The simulation runs for a specified number of generations T. We used T = 5 × 10⁵ throughout this work. Strictly speaking, T is the number of generations in which allele 1 can mutate into allele 2; for later generations this possibility of mutation is disabled. If allele 2 is present at time T, then the simulation is continued until allele 2 is either lost or has reached fixation. Generations are discrete, with mutations, selection, and drift occurring at each generation. During each generation we perform the following steps. First we check if either allele 0 or allele 2 has reached fixation; if so, we stop the simulation, as both cases are absorbing states. Next, for each allele we check for null mutations by drawing a random number from the Poisson distribution where the expected number of events is the null mutation rate μ multiplied by the total number of individuals with the given allele. Mutations from allele 1 to allele 2 are computed in a similar manner, where the expected number of events is U multiplied by the number of allele 1 individuals. Then, after the possible production of the mutant allele 2, selection acts on the fitness of the alleles, where the frequency of each allele is multiplied by its corresponding fitness, [1, 1 + sβ, 1 + s] for alleles [0, 1, 2], giving the new number of alleles in a possibly larger population. Finally, the next population of N individuals is chosen by recursively sampling from the binomial distribution, representing random genetic drift. Allele 0 is first sampled with the mean = (frequency of allele 0), and the (number of trials) = N. Allele 1 is then sampled from the combined allele 1 and 2 individuals. The number of simulations where allele 2 becomes fixed is divided by the total number of simulations, giving an estimate of the fixation probability. The number of simulations for each parameter set was between 10⁸ and 10⁹.

Here, we present the details of our analytic derivations.