Mutational robustness quantifies the effect of random mutations on fitness. When mutational robustness is high, most mutations do not change fitness or have only a minor effect on it. From the point of view of fitness landscapes, robust genotypes form neutral networks of almost equal fitness. Using deterministic population models it has been shown that selection favors genotypes inside such networks, which results in increased mutational robustness. Here we demonstrate that this effect is massively enhanced by recombination. Our results are based on a detailed analysis of mesa-shaped fitness landscapes, where we derive precise expressions for the dependence of the robustness on the landscape parameters for recombining and non-recombining populations. In addition, we carry out numerical simulations on different types of random holey landscapes as well as on an empirical fitness landscape. We show that the mutational robustness of a genotype generally correlates with its recombination weight, a new measure that quantifies the likelihood for the genotype to arise from recombination. We argue that the favorable effect of recombination on mutational robustness is a highly universal feature that may have played an important role in the emergence and maintenance of mechanisms of genetic exchange.
The reshuffling of genetic material by recombination is a ubiquitous part of the evolutionary process across the entire range of organismal complexity. Starting with viruses as the simplest evolving entities, recombination occurs largely at random during the coinfection of a cell by more than one virus strain [1]. For bacteria the mechanisms involved in recombination are already more elaborate and present themselves in the form of transformation, transduction and conjugation [2,3]. In eukaryotic organisms, sexual reproduction is a nearly universal feature, and recombination is often a necessary condition for the creation of offspring. Although its prevalence in nature is undeniable, the evolution and maintenance of sex is surprising since compared to an asexual population, only half of a sexual population is able to bear offspring and additionally a suitable partner needs to be found [4,5]. Whereas the resulting two-fold cost of sex applies only to organisms with differentiated sexes [6], the fact that genetic reshuffling may break up favorable genetic combinations or introduce harmful variants into the genome poses a problem also to recombining microbes that reproduce asexually [7,8]. Since this dilemma was noticed early on in the development of evolutionary theory, many attempts have been undertaken to identify evolutionary benefits of sex and recombination based on general population genetic principles [9–19]. In this article we approach the evolutionary role of recombination from the perspective of fitness landscapes. The fitness landscape is a mapping from genotype to fitness, which encodes the epistatic interactions between mutations and provides a succinct representation of the possible evolutionary trajectories [20]. Previous computational studies addressing the effect of recombination on populations evolving in epistatic fitness landscapes have revealed a rather complex picture, where evolutionary adaptation can be impeded or facilitated depending on, e. g. , the structure of the landscape, the rate of recombination or the time frame of observation [21–26]. Here we focus specifically on the possible benefit of recombination that derives from its ability to enhance the mutational robustness of the population. A living system is said to be robust if it is able to maintain its function in the presence of perturbations [27–31]. In the case of mutational robustness these perturbations are genetic, and the robustness of a genotype is quantified by the number of mutations that it can tolerate without an appreciable change in fitness. Robust genotypes that are connected by mutations therefore form plateaux in the fitness landscape that are commonly referred to as neutral networks [32–35]. Mutational robustness is known to be abundant at various levels of biological organization, but its origins are not well understood. In particular, it is not clear if mutational robustness should be viewed as an evolutionary adaptation, or rather reflects the intrinsic structural constraints of living systems. Arguments in favor of an adaptive origin of robustness were presented by van Nimwegen et al. [32] and by Bornberg-Bauer and Chan [33], who showed that selection tends to concentrate populations in regions of a neutral network where robustness is higher than average. Whereas this result is widely appreciated, the role of recombination for the evolution of robustness has received much less attention. An early contribution that can be mentioned in this context is due to Boerlijst et al. [36], who discuss the error threshold in a viral quasi-species model with recombination and point out in a side note that “in sequence space recombination is always inwards pointing. ” This observation was picked up by Wilke and Adami [37] in a review on the evolution of mutational robustness, where they conjecture that the enhancement of robustness by selection should be further amplified by recombination, because “recombination alone always creates sequences that are within the boundaries of the current mutant cloud. ” At about the same time, de Visser et al. discussed a mechanism based on the spreading of robustness modifier alleles in recombining populations [27] (see also [38]). In fact indications of a positive effect of recombination on robustness had been reported earlier in computational studies of the evolution of RNA secondary structure [39] and 2D lattice proteins [40] in the presence and absence of recombination. In these systems the native folding structure of a given sequence is determined by its global free energy minimum. Due to the restricted number of attainable folds, most structures are degenerate in the sense that many sequences fold into the same structure. These sequences form neutral networks in sequence space. Xia and Levitt [40] consider two scenarios, in which the evolution of the lattice proteins is dominated by mutation and by recombination, respectively. The results show that in the latter case the concentration of thermodynamically stable protein sequences is enhanced, which is qualitatively explained by the fact that recombination tends to focus the sequences near the center of their respective neutral network. Therefore most often a single mutation does not change the folding structure. More recently, Azevedo et al. [41] used a model of gene regulatory networks to investigate the origin of negative epistasis, which is a requirement for the advantage of recombination according to the mutational deterministic hypothesis [13]. In this study a gene network is encoded by a matrix of interaction coefficients. It is defined to be viable if its dynamics converges to a stable expression pattern and non-viable otherwise. Thus the underlying fitness landscape is again neutral. Based on their simulation results the authors argue that recombination of interaction matrices reduces the recombinational load, which in turn leads to an increase of mutational robustness and induces negative epistasis as a byproduct. In effect, then, recombination selects for conditions that favor its own maintenance. Other studies along similar lines have been reviewed in [42]. Taken together they suggest that the positive effect of recombination on robustness may be largely independent of the precise structure of the space of genotypes or the genotype-phenotype map. Indeed, a related scenario has also been described in the context of computational evolution of linear genetic programs [43]. Finally, in a numerical study that is similar to ours in spirit, Szöllősi and Derényi considered the effect of recombination on the mutational robustness of populations evolving on different types of neutral fitness landscapes [44]. Using neutral networks that were either generated at random or based on RNA secondary structure, they found that recombination generally enhances mutational robustness by a significant amount. Moreover, they showed that this observation holds not only for infinite populations but also for finite populations, as long as these are sufficiently polymorphic. The goal of this article is to explain these scattered observations in a systematic and quantitative way. For this purpose we begin by a detailed examination of the simplest conceivable setting consisting of a haploid two-locus model with three viable and one lethal genotype [35]. We derive explicit expressions for the robustness as a function of the rates of mutation and recombination that demonstrate the basic phenomenon and guide the exploration of more complex situations. The two-locus results are then generalized to mesa landscapes with L diallelic loci, where genotypes carrying up to k mutations are viable and of equal fitness [45–48]. Using a communal recombination scheme and previous results for multilocus mutation-selection models, we arrive at precise asymptotic results for the mutational robustness for large L and small mutation rates. Subsequently two types of random holey landscape models are considered, including a novel class of sea-cliff landscapes in which the fraction of viable genotypes depends on the distance to a reference sequence. For the isotropic percolation landscape analytic upper and lower bounds on the robustness are derived. As a first step towards a unified explanation for the effect of recombination on mutational robustness we introduce the concept of the recombination weight, which is a measure for the likelihood of a genotype to arise from a recombination event. In analogy to the classic fitness landscape concept in the context of selection [20], the recombination weight allows one to identify genotypes that are favored by recombination without referring to any specific evolutionary dynamics. We show that recombination weight correlates with mutational robustness for the landscape structures used in this work, thus providing a mechanistic basis for the enhancement of robustness by recombination. Finally, using an empirical fitness landscape as an example, we quantify the competition between selection and recombination as a function of recombination rate. Throughout we describe the evolutionary dynamics by a deterministic, discrete time model that will be introduced in the next section.
We consider a haploid genome with L loci and the corresponding genotype is represented by a sequence σ = (σ1, σ2, …, σL) of length L. The index i labels genetic loci and each locus carries an allele specified by σi. Here we rely on binary sequences, which means that there are only two different alleles σi ∈ {0,1}. This can be either seen as a simplification in the sense that only two alleles are assumed to exist, or in the sense that the genome consisting of all zeros describes the wild type, and the 1’s in the sequence display mutations for which no further distinctions are made. The resulting genotype space is a hypercube of dimension L, where the 2L genotypes represent vertices, and two genotypes that differ at a single locus and are mutually reachable by a point mutation are connected by an edge. A metric is introduced by the Hamming distance d (σ, κ) =∑i (1−δσiκi), (1) which measures the number of point mutations that separate two genotypes σ and κ. Here and in the following the Kronecker symbol is defined as δxy = 1 if x = y and δxy = 0 otherwise. The genotype σ¯ at maximal distance d (σ, σ¯) =L from a given genotype σ is called its antipodal, and can be defined by σ¯i=1−σi. Finally, in order to generate a fitness landscape, a (Wrightian) fitness value wσ is assigned to each genotype. The forces that drive evolution are selection, mutation and recombination. To model the dynamics we use a deterministic, discrete-time model with non-overlapping generations, which can be viewed as an infinite population limit of the Wright-Fisher model. Demographic stochasticity or genetic drift is thus neglected. Numerical simulations of evolution on neutral networks have shown that the infinite population dynamics is already observable for moderate population sizes, which justifies this approximation [32,44]. We will return to this point in the Discussion. Once the frequency fσ (t) of a genotype σ at generation t is given, the frequency at the next generation is determined in three steps representing selection, mutation, and recombination. After the selection step, the frequency qσ (t) is given as qσ (t) =wσw¯ (t) fσ (t), (2) where w¯≡∑σwσfσ (t) is the mean population fitness at generation t. After the mutation step, the frequency pσ (t) is given as pσ (t) =∑κUσκqκ (t), (3) where Uσκ is the probability that an individual with genotype κ mutates to genotype σ in one generation. Here, we assume that alleles at each locus mutate independently, and the mutation probability μ is the same in both directions (0 → 1 and 1 → 0) and across loci. This leads to the symmetric mutation matrix Uσκ= (1−μ) L−d (σ, κ) μd (σ, κ). (4) In order to incorporate recombination we have to consider the probability that two parents with respective genotypes κ and τ beget a progeny with genotype σ by recombination. This is represented by the following equation: fσ (t+1) =∑κτRσ|κτpκ (t) pτ (t). (5) Descriptively speaking, two genotypes (κ and τ) are taken to recombine with a probability that is equal to their frequency in the population (after selection and mutation). The probability for the offspring genotype σ is then given by Rσ|κτ. These probabilities depend of course on the parent genotypes κ and τ but also on the recombination scheme. Here we consider a uniform and a one-point crossover scheme; see Fig 1 for a graphical representation. These two represent extremes in a spectrum of possible recombination schemes. Nevertheless we will show that both lead to qualitatively similar results in the regimes of interest. In the case of uniform crossover the recombination probabilities are given by Rσ|κτ=r2L (∏iL (δσiκi+δσiτi) ) +1−r2 (δσκ+δστ) (6) and in the case of one point crossover the probabilities can be written as Rσ|κτ=r2 (L−1) ∑n=1L−1[ (∏m=1nδσmκm) (∏m=n+1Lδσmτm) + (∏m=1nδσmτm) (∏m=n+1Lδσmκm) ]+1−r2 (δσκ+δστ). (7) In both equations a variable r ∈ [0,1] appears which describes the recombination rate. For r = 0 no recombination occurs and fσ (t + 1) is the same as pσ (t). For r = 1 recombination is a necessary condition for the creation of offspring (obligate recombination). But also intermediate values of r can be chosen as they occur in nature, e. g. , for bacteria and viruses. In the following we are mostly interested in the equilibrium frequency distribution fσ* of a population, which is determined by the stationarity condition fσ (t+1) =fσ (t) =fσ* (8) for all genotypes σ. From the point of view of fitness landscapes the occurrence of mutational robustness implies that fitness values of neighboring genotypes are degenerate, giving rise to neutral networks in genotype space [29,32–35]. In order to model this situation we use two-level landscapes that only differentiate between genotypes that are viable (wσ = 1) or lethal (wσ = 0). Any selective advantage between viable genotypes is assumed to be negligible. The mutational robustness of a population can then be measured by the average fraction of viable point mutations in an individual, which depends on the population distribution in genotype space [32–34]. It increases if the population mainly adapts to genotypes for which most point mutations are viable. Therefore we define mutational robustness m as the average fraction of viable point mutations of a population, m≡∑σ∈Vmσfσ*withmσ≡nσL. (9) Here the sum is over the set V of all viable genotypes and nσ is the number of viable point mutations of genotype σ. We will refer to mσ as the mutational robustness of the genotype. The expression is normalized by the total number of loci L, since in an optimal setting the entire population has L viable genotypic neighbors and mσ = 1 for all σ ∈ V. The value of m is thus constrained to be in the range [0,1]. We weight the genotypes by their stationary frequencies fσ*, since we want to determine the mutational robustness of populations that are in equilibrium with their environment. In order to elucidate the interplay of recombination and mutational robustness it will prove helpful to introduce a representation of how recombination can transfer genotypes into each other. The number of distinct genotypes that two recombining genotypes are able to create depends on their Hamming distance. In particular, the recombination of two identical genotypes does not create any novelty, whereas a genotype and its antipodal are able to generate all possible genotypes through uniform crossover. Here we introduce a measure which expresses how many pairs of viable genotypes are able to recombine to a specific genotype. It is complementary to the mutational robustness, in the sense that instead of counting the viable mutation neighbors of a genotype, the size of its recombinational neighborhood of viable recombination pairs is determined. The recombinational neighborhood depends on the recombination scheme and the distribution of viable genotypes in the genotype space. For a given recombination scheme the probability for a genotype σ to be the outcome of recombination of two genotypes κ, τ is given by the recombination tensor Rσ|κτ. The recombination weight λσ is therefore obtained by summing the recombination tensor over all ordered pairs of viable genotypes, λσ=12L∑κ∈V, τ∈VRσ|κτ. (10) It can be seen from (5) that λσ = 1 when all genotypes are viable, and hence the normalization by 2L ensures that the recombination weight lies in the range [0,1]. Under this normalization, the recombination weights sum to ∑σ λσ = |V|2/2L, where |V| stands for the number of viable genotypes. In the following the genotype maximizing λσ will be referred to as the recombination center of the landscape. Since neutral landscapes only differentiate between viable (unit fitness) and lethal (zero fitness) genotypes, the recombination weight (10) can alternatively be written as a sum over all ordered pairs of genotypes whereby the recombination tensor is multiplied by the pair’s respective fitness, λσ=12L∑κ, τRσ|κτwκwτ. (11) In this way the concept naturally generalizes to arbitrary fitness landscapes. In the absence of recombination (r = 0) the recombination weight (11) of a genotype is simply proportional to its fitness, λσ=w˜wσ, where w˜=2−L∑σwσ is the unweighted average fitness. Within our recombination schemes, the recombination tensor depends linearly on r and, by definition, so does the recombination weight. Accordingly, for general r the recombination weight interpolates linearly between the limiting values at r = 0 and r = 1. Since λσ for r = 0 is known, the remaining task will be to find λσ for r = 1.
The simplest fitness landscape to study the mutational robustness of a population would be the haploid two-locus model in which all but one genotype are viable [35]; see Fig 2 for a graphical representation of the model. In this setting the population gains mutational robustness if the frequency of the genotype (0,0) for which both point mutations are viable increases relative to the genotypes (0,1) and (1,0). This model has been analyzed previously using a unidirectional mutation scheme where reversions 1 → 0 are suppressed [58,59]. As a consequence, selection cannot contribute to mutational robustness because the genotype (0,0) goes extinct in the absence of recombination. Here we consider the case of bidirectional, symmetric mutations in which both selection and recombination contribute to robustness. A comparison of the two mutation schemes is provided in S1 Appendix. We proceed to solve the equilibrium condition Eq (8). Since the equilibrium genotype frequencies f01* and f10* are the same due to the symmetry of the landscape and the mutation scheme, the recombination step at stationarity reads f00*=p00−ρ (p00p11−p10p01) ⇔f0=p0−ρD, f10/01*=p10/01+ρ (p00p11−p10p01) ⇔f1=p1+2ρD, f11*=p11−ρ (p00p11−p10p01) ⇔f2=p2−ρD, (12) where pσ is the (equilibrium) frequency of genotype σ after the mutation step, fi and pi are the corresponding lumped frequencies [60] of all genotypes with i 1’s, and D≡p00p11−p10p01=p0p2−p12/4 is the linkage disequilibrium after the mutation step. Notice that the one-point and uniform crossover schemes give the same equation form except that the parameter ρ is given by ρ = r in the case of one-point crossover and ρ = r/2 for uniform crossover. However, we would like to emphasize that this is a mere coincidence of the two-locus model which disappears as soon as L is larger than 2. The lumped frequencies qi of all genotypes with i 1’s after the selection step are given by q0=f01−f2, q1=f11−f2, q2=0. (13) Applying the mutation step we obtain p0=q0 (1−μ) 2+μ (1−μ) q1=μ (1−μ) + (1−μ) (1−2μ) q0, p1=q1[ (1−μ) 2+μ2]+2μ (1−μ) q0=1−2μ+2μ2− (1−2μ) 2q0, p2=μ (1−μ) q1+μ2q0=μ (1−μ) −μ (1−2μ) q0, D=p0p2−p12/4=−14 (1−2μ) 2 (1−q0) 2, (14) where we have used the normalization q0 + q1 = 1 to express the right hand sides in terms of q0. Putting everything together, the problem is reduced to solving the following third order polynomial equation for q0,0=q0 (1−f2) −f0=q0 (1−p2+ρD) −p0+ρD=ρ4 (1−2μ) 2 (1−q0) 2 (1+q0) +μ[1−2q0−q02−μ (1−2q0) (1+q0) ], (15) from which we can in principle find exact analytic expressions for fσ*. However, it is difficult to extract useful information from the exact solution. In the following we will therefore provide approximate solutions. If we neglect recombination (ρ = 0), we obtain the following equilibrium genotype frequency distribution: f00* (ρ=0) =1−μ28−16μ+9μ2−12 (2−5μ+3μ2) ≈ (2−1) + (52−22) μ+O (μ2), f01/10* (ρ=0) =14 (4−9μ+6μ2) −1−2μ48−16μ+9μ2≈ (1−12) + (32−94) μ+O (μ2), f11* (ρ=0) =μ2 (4−3μ−8−16μ+9μ2) ≈ (2−2) μ+O (μ2). (16) When ρ = 1, which corresponds to the one-point crossover scheme with r = 1, linkage equilibrium (f00 f11 = f10 f01) is restored after one generation [55]. Accordingly, we can treat each locus independently and get rather simple expressions for fσ* as f00* (ρ=1) =14 (2+μ−μ2+4μ) 2≈1−2μ+2μ+O (μ3/2), f01/10* (ρ=1) =14 (2+μ−μ2+4μ) (μ2+4μ−μ) ≈μ−3μ2+O (μ3/2), f11* (ρ=1) =14 (μ2+4μ−μ) 2≈μ−O (μ3/2). (17) We depict the equilibrium solutions for the above two cases in Fig 3. Now, the mutational robustness m=12 (2f00*+f10*+f01*) =f00*+f10*=f0+12f1 (18) for the above two cases is obtained as m (μ, ρ=0) =14 (μ+8−16μ+9μ2) ≈12− (12−14) μ+O (μ2), (19) m (μ, ρ=1) =12 (2+μ−μ2+4μ) ≈1−μ+μ2+O (μ3/2), (20) which is depicted in Fig 4. These results encapsulate in a simple form the main topic of this paper. Selection alone (ρ = 0) leads to a moderate increase of robustness from the baseline value m=12 corresponding to a random distribution over genotypes, which is attained at μ=12, to m=12 for μ → 0. In contrast, for recombining populations (ρ = 1) robustness is massively enhanced at small mutation rates due to the strong frequency increase of the most robust genotype (0,0) and reaches the maximal value m = 1 at μ = 0. The underlying mechanism is analogous to Kondrashov’s deterministic mutation hypothesis, which posits that recombination makes selection against deleterious mutations more effective when these interact synergistically [13]. In the present case recombination increases the frequency of the double mutant genotype (1,1), which is subsequently purged by selection, and thereby effectively drives the frequency of the allele 1 at both loci to zero. The enhancement of the frequency of the genotype (0,0) by recombination is also reflected in the recombination weights, which take on the values λ00=34+ρ4, λ01=λ10=34−ρ4, λ11=ρ4. (21) Thus the genotype (0,0) is the recombination center of the two-locus landscape. Next we investigate how mutational robustness varies with μ for intermediate recombination rates, assuming that μ is small. As can be seen from Eq (15), the asymptotic behavior of the solution for small ρ and μ depends on which of the two parameters is smaller. We first consider the case ρ ≪ μ ≪ 1. Defining l = ρ/ (4μ) ≪ 1, Eq (15) is approximated by 0=l (1−q0) 2 (1+q0) +1−2q0−q02−μ (1−2q0) (1+q0), (22) where we kept terms up to O (μ), since we have not determined whether l is smaller than μ or not. Since q0=2−1 is the solution of Eq (22) for l = μ = 0, we set q0=2−1+al+bμ and solve the equation to leading order, which gives q0≈2−1+ (3−22) l− (32−2) μ. (23) The mutational robustness then follows as m=f0+f12=12+p0−p22=12+ (1−2μ) q02≈12+3−222l− (12−14) μ, (24) which is consistent with our previous result for ρ = 0; see Eq (19). We note that in this regime it is sufficient for the recombination rate to be of order O (μ2) to compensate the negative effect of mutations on mutational robustness, as the two effects cancel when ρ = ρc with ρc=2 (5+42) μ2≈21. 3×μ2. (25) In the regime ρ ≫ μ, Eq (15) is approximated as (1−4μ) (1−q0) 2 (1+q0) +s (1−2q0−q02) =0, (26) with s = 4μ/ρ. Again we have kept terms up to O (μ) because μ and s are of the same order if ρ = O (1). Since the solution of Eq (26) for μ = s = 0 is q0 = 1, we set q0 = 1 − α with α ≪ 1. Inserting this into Eq (26), we get α≈s. Since α ≫ μ, q0=1−s is the approximate solution to leading order. Hence m=12+ (1−2μ) q02≈1−s4=1−μρ, (27) which is again consistent with our previous result for ρ = 1 in Eq (20). The square root dependence on μ/ρ derives from the corresponding behavior of the genotype frequency f00* and has been noticed previously in the model with unidirectional mutations [58,59]. For arbitrary ρ and μ, we have to use the full Eq (15). Fig 5 illustrates the behaviour of mutational robustness as a function of the recombination rate for different mutation rates and both recombination schemes. For small μ, a low rate of recombination suffices to bring the robustness close to its maximal value m = 1. More precisely, according to Eq (27), a robustness m > 1 − ϵ is reached for recombination rates ρ > μ/ϵ2. To summarize, we have seen that analytic results for the two-locus model are easily attainable. For multi-locus models it is much more challenging to derive analytical results, particularly in the presence of recombination. By way of contrast the dynamics induced only by mutation and selection are easier to understand: While mutations increase the genotype diversity in the population, fitter ones grow in frequency through selection, which reduces diversity. Although one might expect that recombination would increase diversity, a number of studies have shown that recombination is more likely to impede the divergence of populations. Recombining populations tend to cluster on single genotypes or in a limited region of a genotype space and furthermore the waiting times for peak shifts in multipeaked fitness landscapes diverge at a critical recombination rate [22,26,54–56,61]. The results for the two-locus model presented above are consistent with this behaviour, as the genotype heterogeneity of the population decreases with increasing recombination rate (S1 Fig). In the following we will investigate how the focusing effect of recombination enhances the mutational robustness of the population in three different multi-locus models. In the mesa landscape it is assumed that up to a certain number k of mutations all genotypes are functional and have unit fitness, whereas genotypes with more than k mutations are lethal and have fitness zero [48]. Hence the fitness landscape is defined as wσ={1, ifdσ≤k, 0, otherwise, (28) where dσ is the Hamming distance to the wild-type sequence (0,0, …, 0) or, equivalently, the number of loci with allele 1. We will refer to k as the mesa width or as the critical Hamming distance. Such a scenario can for example be observed in the evolution of regulatory motifs, where the fitness depends on the binding affinity of the regulatory proteins and dσ corresponds to the number of mismatches compared to the original binding motif [45,47]. The two-locus model discussed in the preceding section corresponds to the mesa landscape with critical Hamming distance k = 1 and sequence length L = 2. Here we ask to what extent the behavior observed for the two-locus model generalizes to longer sequences and variable k. Numerical simulations suggest that the strong increase of mutational robustness with recombination rate indeed persists in the general setting, and the particular recombination scheme seems to have only a minor influence; see Fig 6. Whereas an analytical treatment for general L, k and intermediate recombination rates appears to be out of reach, accurate approximations are available in the limiting case of strong recombination or of no recombination, assuming mutation rate is small. The full derivations for both cases can be found in S1 Appendix. In the following we summarize the main results. In the percolation landscape genotypes are randomly chosen to be viable (wσ = 1) with probability p and lethal (wσ = 0) with probability 1 − p. An interesting property of the percolation model is the emergence of two different landscape regimes [49,63–65]. Above the percolation threshold pc, viable genotypes connected by single mutational steps form a cluster that extends over the whole landscape, whereas below pc only isolated small clusters appear. Since the percolation threshold depends inversely on the sequence length, pc≈1L, for large L a small fraction of viable genotypes suffices to create large neutral networks. This allows a population to evolve to distant genotypes without going through lethal regions, and correspondingly the percolation model is often used to study speciation [35,49]. A network representation of the percolation model is shown in Fig 8. The algorithm used to generate this visual representation is explained in S1 Appendix. Fig 9 shows three exemplary stationary genotype frequency distributions on the landscape depicted in Fig 8. In the absence of recombination the equilibrium frequency distribution is unique, but in the presence of recombination the non-linearity of the dynamics implies that multiple stationary states may exist [54,55,61]. Fig 9 displays two stationary distributions for r = 1 which are accessed from different initial conditions. It is visually apparent that the recombining populations are concentrated on a small number of highly connected genotypes, leading to a significant increase of mutational robustness. To quantify this effect, the average mutational robustness m¯ is calculated as a function of the recombination rate according to the following numerical protocol: A percolation landscape for given L and p is generated and the initial population is distributed uniformly among all genotypes. The population is evolved in the absence of recombination (r = 0) until the unique equilibrium frequency distribution is reached, for which the mutational robustness m is calculated. Next the recombination rate is increased by predefined increments. After increasing r, the population is again evolved using the stationary state obtained before the increment of r as the initial condition, until it reaches a new stationary state for which the mutational robustness is measured. When the recombination rate has reached r = 1, a new percolation landscape is generated and the process starts all over again. This is done for an adjustable number of runs over which the average is taken. The results of such a computation are shown in Fig 10. Similar to the mesa landscapes, a strong increase of mutational robustness is observed already for small rates of recombination, and the effect is largely independent of the recombination scheme. However, in contrast to the mesa landscape the robustness does not reach its maximal value m = 1 for r = 1 and small μ. This reflects the fact that maximally connected genotypes with mσ = 1 are very rare at this particular value of p. For the purpose of comparison we also determined the average mutational robustness m¯0 of a uniform population distribution for the percolation model. Conditioned on the number v = |V| of viable genotypes and assuming that v ≥ 1, we have m0 (v, L) = n (v, L) /L, where n (v, L) is the average number of viable neighbors of a viable genotype. The latter is given by the expression n (L, v) = (v−1) L2L−1, (38) since for a given viable genotype there are v − 1 remaining genotypes, each of which has the probability L/ (2L − 1) to be a neighboring one. Taking into account that the number of viable genotypes is binomially distributed with parameter p and that the empty hypercube (v = 0) should yield m0 = 0 we obtain m¯0=∑v=12L (v−1) 2L−1 (2Lv) pv (1−p) 2L−v=2Lp−1+ (1−p) 2L2L−1, (39) which simplifies to m¯0=p when 2Lp ≫ 1. Note that the condition 2Lp ≫ 1 is naturally satisfied beyond the percolation threshold. Fig 11 illustrates that the dynamics induced by mutation and selection already increase mutational robustness compared to m¯0 and that the addition of recombination even further increases mutational robustness for all values of p. The figure also displays the expected maximum number of viable neighbors of any genotype in the landscape, m¯max, which provides an upper bound on the robustness. The fact that the numerically determined robustness remains below this bound for all p shows that the ability of recombination to locate the most connected genotype is limited. In S1 Appendix it is shown that limL→∞m¯max=1 for p>12. As outlined above, the algorithm used to generate Figs 10 and 11 computes the mutational robustness of a particular stationary frequency distribution of the recombining population which is smoothly connected to the unique non-recombining stationary state. Although one expects this state to be representative in the sense of being reachable from many initial conditions, for large enough r there can be multiple stationary states that will generally display different robustness (see Fig 9). To illustrate this point, S7 Fig shows the results of a simulation of the percolation model where all stationary states were identified using localized initial conditions, and the mutational robustness was computed separately for each state. Whereas on average the mutational robustness is always enhanced by recombination, there are rare instances when recombination reduces the robustness compared to the non-recombining case. This may happen, for example, if recombination traps the population on a small island of viable genotypes [22,26,55,56]. In this section we introduce a novel class of fitness-landscape models (to be called sea-cliff landscapes) that interpolates between the mesa and percolation landscapes. Similar to the mesa landscape, the fitness values of the sea-cliff model are determined by the distance to a reference genotype κ*. The model differs from the mesa landscape in that it is not assumed that all genotypes have zero fitness beyond a certain number of mutations. Instead, the likelihood for a mutation to be lethal (to “fall off the cliff”) is taken to increase with the Hamming distance from the reference genotype. This is mathematically realized by a Heaviside step function θ (x) that contains an uncorrelated random contribution ησ and the distance measure d (σ, κ*), wσ=θ[ησ−d (σ, κ*) ]={1, ifησ>d (σ, κ*), 0, ifησ<d (σ, κ*). (40) This construction is similar in spirit to the definition of the Rough-Mount-Fuji model [66,67]. The average shape of the landscape can be tuned by the mean c and the standard deviation s of the distribution of the random variables ησ, which we assume to be Gaussian in the following. The average fitness at distance d from the reference sequence is then given by w¯ (d) =Prob (wσ=1) =12[1−erf (d−cs2) ], (41) where erf (x) is the error function. Note that the mesa landscape is reproduced if we take the limit s → 0 for fixed c in the range k < c < k + 1 and the percolation landscape is reproduced if we take a joint limit s, |c| → ∞ with c/s fixed. To fix c and s we introduce two distances d< and and d> such that w¯ (d<) =0. 99 and w¯ (d>) =0. 01, which leads to the relations c=12 (d<+d>) ands≈0. 215 (d>−d<). (42) The model can be generalized to include several predefined reference sequences, w (σ) =θ{∑κ*θ[ησ, κ*−d (σ, κ*) ]}, (43) which allows to create a genotype space with several highly connected clusters. Depending on the Hamming distance between the reference sequences and the variables c and s, clusters can be isolated or connected by viable mutations. Fig 12 shows stationary states in the absence and presence of recombination for two different sea-cliff landscapes with one and two reference genotypes, respectively. Similar to the other landscape models, mutational robustness increases strongly with recombination, due to a population concentration within a neutral cluster. In the presence of two reference genotypes the recombining population should be concentrated within a single cluster. Otherwise lethal genotypes would be predominantly created through recombination of genotypes on different clusters. This observation can also be interpreted in the context of speciation due to genetic incompatibilities [49,61]. Without recombination genotypes on both clusters have a nonvanishing frequency, but still the larger cluster is more populated. In contrast to the percolation landscape, robustness reaches a value close to unity for large r, because highly connected genotypes are abundant close to the reference sequence (S8 Fig). Comparing Figs 6 and 10 and S8 Fig, the dependence of mutational robustness on the recombination rate is seen to be strikingly similar. Despite the very different landscape topographies, in all cases a small amount of recombination gives rise to a massive increase in robustness compared to the non-recombining baseline. For the mesa landscape this effect can be plausibly attributed to the focusing property of recombination, which counteracts the entropic spreading towards the fitness brink and localizes the population inside the plateau of viable genotypes. In the case of the holey landscapes, however, it is not evident that focusing the population towards the center of its genotypic range will on average increase robustness, since viable and lethal genotypes are randomly interspersed. To establish the relation between recombination and mutational robustness on the level of individual genotypes, in Fig 13 we plot the recombination weight of each genotype against its robustness mσ. A clear positive correlation between the two quantities is observed both for percolation and sea-cliff landscapes. Additionally we differentiate between viable and lethal genotypes. In the percolation landscape viable genotypes are uniformly distributed in the genotype space, which implies that lethal and viable genotypes have on average the same number of viable point-mutations. Nevertheless the recombination weight of viable genotypes is larger. The fitness of a genotype influences its own recombination weight, because the genotype itself is a possible parental genotype in the recombination event. In non-neutral fitness landscapes the redistribution of the population through recombination competes with selection responding to fitness differences, and the generalized definition (11) of the recombination weight captures this interplay. To exemplify the relation between recombination weight and mutational robustness in this broader context, we use an empirical fitness landscape for the filamentary fungus Aspergillus niger originally obtained in [68]. In a nutshell, two strains of A. niger (N411 and N890) were fused to a diploid which is unstable and creates two haploids by random chromosome arrangement. Both strains are isogenic to each other, except that N890 has 8 marker mutations on different chromosomes, which were induced by low UV-radiation. Through this process 28 = 256 haploid segregants can theoretically be created of which 186 were isolated in the experiment. As a result of a statistical analysis it was concluded that the missing 70 haploids have zero fitness [69]. In order to illustrate the fitness landscape, a network representation is employed where genotypes are arranged in a plane according to their fitness and their Hamming distance to the wild type, which in this case is the genotype of maximal fitness. In Fig 14A and 14B node sizes are adjusted to the recombination weights and mutational robustness of genotypes, respectively, in order to display the distribution of these quantities. In accordance with the analyses for neutral fitness landscapes, a clear correlation between the recombination weights and mutational robustness is shown in Fig 14C. Since fitness values are not binary we further consider the correlation between the recombination weights and fitness values (Fig 14D). The recombination center is one of the maximally robust genotypes with mσ = 1, but it is not the fittest within this group. The wild type has maximal fitness but, by comparison, lower robustness (mσ = 7/8). Fig 15 highlights how the recombination weights change as a function of the recombination rate and how this affects the stationary state of a population. For small recombination rates the recombination weight of each genotype mainly depends on its own fitness, and therefore the wild type coincides with the recombination center. With increasing recombination rate the connectivity of the surrounding genotype network becomes more important and the recombination center switches to a genotype at Hamming distance d = 2. In contrast to the numerical protocol described previously, in the simulations used to generate Fig 15D–15F the population is reset to a uniform distribution before the recombination rate is increased. Otherwise the population would continue to adapt to the wild type, which has the highest fitness and from which it cannot escape because of peak trapping [22,26]. Starting from an initially uniform distribution the population will adapt to one of three possible final genotypes which depend on the recombination rate. For small and large recombination rates the most abundant genotype coincides with the recombination center (Fig 15D and 15F), whereas for intermediate recombination rates the population chooses another genotype that is also located at Hamming distance d = 2 but has higher fitness (Fig 15E). The recombination center ultimately dominates the population, not only because it is maximally connected (mσ = 1), but also because the genotypes that it is connected to have high fitness. In this sense the sequence of transitions in the most abundant genotype that occur with increasing recombination rate is akin to the scenario described previously in non-recombining populations as the “survival of the flattest” [48,70]. Along this sequence mutational robustness increases monotonically whereas the average fitness of the population actually declines (S9 Fig).
Despite a century of research into the evolutionary bases of recombination, a general mechanism explaining the ubiquity of genetic exchange throughout the domains of life has not been found [17,18]. Even within the idealized scenario of a population evolving in a fixed environment, whether or not recombination speeds up adaptation and leads to higher fitness levels depends in a complicated way on the structure of the fitness landscape and the parameters of the evolutionary dynamics [21–26]. The most important finding of the present work is that, by comparison, the effect of recombination on mutational robustness is much simpler and highly universal. Irrespective of the number of loci, the structure of the fitness landscape or the recombination scheme, recombination leads to a significant increase of robustness that is usually much stronger than the previously identified effect of selection [32–34]. This suggests that the evolution of recombination may be closely linked to the evolution of robustness, and that similar selective benefits are involved in the two cases. Although the relation of robustness to evolutionary fitness is subtle and not fully understood [27], it has been convincingly argued that robustness enhances evolvability and hence becomes adaptive in changing environments [29,31,71,72]. A common perspective on recombination, robustness and evolvability can help to develop novel hypotheses about the evolutionary origins of these phenomena that can be tested in future computational or empirical studies. On a quantitative level, we have shown that robustness generally depends on the ratio of recombination to mutation rates, and that the robustness-enhancing effect saturates when r ≫ μ. This observation highlights the importance of r/μ as an evolutionary parameter. Interestingly, even in bacteria and archaea, which have traditionally been regarded as essentially non-recombining, the majority of species displays values of r/μ that are significantly larger than one [73–75]. Similarly, a recent study of the evolution of Siphoviridae phages revealed a ratio of recombination events to mutational substitutions of about 24 [76]. In eukaryotes this ratio is expected to be considerably higher [40]. This indicates that most organisms maintain a rate of recombination that is sufficient to reap its evolutionary benefits in terms of increased robustness. In order to clarify the mechanism through which recombination enhances robustness, we have introduced the concept of the recombination weight, which is a measure for the likelihood of a genotype to arise from the recombination of two viable parental genotypes. The recombination weight defines a “recombination landscape” over the space of genotypes which is similar in spirit to, but distinct from, previous mathematical approaches to conceptualizing the way in which recombining populations navigate a fitness landscape [77]. It is complementary to the more commonly used notion of a recombination load, which refers to the likelihood for a viable genotype to recombine to a lethal one [41,42]. In many cases the maximum of the recombination weight correctly predicts the most populated genotype in a recombining population at low mutation rate. Moreover, the concept generalizes to non-neutral landscapes and thus permits to address situations where selection and recombination compete. Provided recombination weight is correlated with mutational robustness for the individual genotypes, this explains the positive effect of recombination on the population-level robustness. Whether or not such a correlation exists will generally depend on the structure of the fitness landscapes. For simple neutral landscapes such as the mesa landscape it is an immediate consequence of the focusing property of recombination, but for more complex neutral networks the relationship between the two quantities is nontrivial and needs to be studied on a case-by-case basis. Although a positive correlation was observed numerically both for the holey landscapes and the empirical landscape considered in this work, it is not difficult to construct landscapes where the genotypes with high recombination weight are not highly robust. As a simple but instructive example, in S10 Fig we show results for an ‘atoll’ landscape where a ring of viable genotypes surrounds a central hole of lethals. Throughout this work the effects of genetic drift have been neglected. We expect that our results will be applicable to finite populations as long as the population is sufficiently diverse rather than being monomorphic. This requires the population-wide mutation rate NμL to be much larger than unity [32,44]. If NμL ≪ 1 the population is almost always monomorphic and recombination has no effect. In this regime the population explores the fitness landscape as a random walker and the observed mutational robustness is the uniform robustness m0. In S11 Fig we present the results of finite population simulations on a mesa landscape, which show a sharp transition from the random walk regime to the behavior predicted by the deterministic theory when NμL ∼ 1. Future work should be directed towards extending the present investigation to more realistic genotype-phenotype maps arising, for example, from the secondary structures of biopolymers such as RNA or proteins [39,40,44], or from simple genetic, metabolic or logical networks [29,41,43,78]. There is ample evidence from numerical studies that a favorable effect of recombination on mutational robustness is present also in these more complex systems, but a detailed analysis of the underlying mechanism has not been carried out. This would entail, in particular, the generalization to genotype spaces composed of sequences carrying more than two alleles per site. We expect that at least part of the analysis for the mesa landscapes carries over to this setting, and in fact some results for the non-recombining case have already been obtained [48]. More importantly, the role of the topology of the corresponding neutral networks in shaping the correlation between recombination weight and robustness needs to be explored systematically. Research along these lines will help to corroborate the relationship between recombination and robustness that we have sketched, and to further elucidate the origins of these two pervasive features of biological evolution.
Two long-standing and seemingly unrelated puzzles in evolutionary biology concern the ubiquity of sexual reproduction and the robustness of organisms against genetic perturbations. Using a theoretical approach based on the concept of a fitness landscape, in this article we argue that the two phenomena may in fact be closely related. In our setting the hereditary information of an organism is encoded in its genotype, which determines it to be either viable or non-viable, and robustness is defined as the fraction of mutations that maintain viability. Previous work has demonstrated that the purging of non-viable genotypes from the population by natural selection leads to a moderate increase in robustness. Here we show that genetic recombination acting in combination with selection massively enhances this effect, an observation that is largely independent of how genotypes are connected by mutations. This suggests that the increase of robustness may be a major driver underlying the evolution of sexual recombination and other forms of genetic exchange throughout the living world.
biogeography
organismal evolution
ecology and environmental sciences
classical mechanics
fluid mechanics
population genetics
variant genotypes
microbiology
genetic mapping
mutation
microbial evolution
dna
population biology
geography
fluid dynamics
phylogeography
evolutionary genetics
continuum mechanics
viral evolution
physics
biochemistry
point mutation
percolation
nucleic acids
natural selection
virology
earth sciences
heredity
genetics
biology and life sciences
physical sciences
dna recombination
evolutionary biology
evolutionary processes
Recombination and mutational robustness in neutral fitness landscapes