## Abstract

For different fitness mutational models, with epistasis introduced, we simulated the consequences of drift (*D* scenario) or mutation, selection, and drift (*MSD* scenario) in populations at the *MSD* balance subsequently subjected to bottlenecks of size *N* = 2, 10, 50 during 100 generations. No “conversion” of nonadditive into additive variance was observed, all components of the fitness genetic variance initially increasing with the inbreeding coefficient *F* and subsequently decreasing to zero (*D*) or to an equilibrium value (*MSD*). In the *D* scenario, epistasis had no appreciable effect on inbreeding depression and that on the temporal change of variance components was relevant only for high rates of strong epistatic mutation. In parallel, between-line differentiation in mean fitness accelerated with *F* and that in additive variance reached a maximum at *F* ∼ 0.6–0.7, both processes being intensified by strong epistasis. In the *MSD* scenario, however, the increase in additive variance was smaller, as it was used by selection to purge inbreeding depression (*N* ≥ 10), and selection prevented between-line differentiation. Epistasis, either synergistic or antagonistic (this leading to multiple adaptive peaks), had no appreciable effect on *MSD* results nor, therefore, on the evolutionary rate of fitness change.

THE roles of genetic drift and natural selection in shaping the genetic variation of fitness due to segregation at epistatic loci have often been discussed since Wright's (1931) pioneering treatment of the subject. In general, the pertinent analyses have been usually elaborated within an analytical framework where changes in the mean and the components of the genetic variance exclusively due to drift were first considered, this being followed by an examination of the conditions that may subsequently allow for a more rapid selection response and/or facilitate the movement of populations to new adaptive peaks.

Theoretically, it is well known that the contribution of neutral additive loci to the additive genetic variance of metric traits in populations decreases linearly as the inbreeding coefficient *F* increases, until it ultimately vanishes when fixation is attained (Wright 1951). For neutral nonadditive loci, however, that contribution may initially increase until a critical *F* value is reached and then subsequently decline to zero. This is the case of simple dominant loci (Robertson 1952; Willis and Orr 1993), and it also applies to two-locus models showing either additive × additive epistasis (Cockerham and Tachida 1988; Goodnight 1988) or more complex epistasis involving dominance at the single-locus level (Cheverud and Routman 1996; López-Fanjul *et al*. 1999, 2000; Goodnight 2000). Furthermore, those models have been extended to cover multiple additive × additive epistatic systems (Barton and Turelli 2004, López-Fanjul *et al*. 2006).

In parallel, laboratory experiments have also studied the impact of population bottlenecks on the additive variance of metric traits (see reviews by López-Fanjul *et al*. 2003 and Van Buskirk and Willi 2006). For morphological traits not strongly correlated with fitness, a decrease in their additive variance together with little or no inbreeding depression was often observed, both results being compatible with the corresponding additive expectations and suggesting that the standing variation of those traits is mainly controlled by quasi-neutral additive alleles. Using typical estimates of mutational parameters, Zhang *et al*. (2004) showed that these experimental results can be explained by assuming a model of pleiotropic and real stabilizing selection acting on the pertinent trait. On the other hand, life-history traits closely connected to fitness usually show strong inbreeding depression and a dramatic increase in additive variance after a brief period of inbreeding or bottlenecking, indicating that much of that variance should be due to deleterious recessive alleles segregating at low frequencies. However, it should be kept in mind that experimental results cannot discern between simple dominance and dominance with additional epistasis as causes of inbreeding-induced changes in the additive variance.

In their discussion of the shifting-balance theory (Wright 1931), Wade and Goodnight emphasized the evolutionary importance of the “conversion” of epistatic variance into additive variance, proposing that drift-induced excesses in the additive variance for fitness available to selection could enhance the potential for local adaptation, a phenomenon that was not discussed in the original formulation of Wright's theory (Wade and Goodnight 1998; Goodnight and Wade 2000; but see Coyne *et al*. 1997, 2000). However, the additive variance is inflated only under restrictive conditions that often involve low-frequency deleterious recessive alleles (Robertson 1952; López-Fanjul *et al*. 2002), so that a drift-induced excess in the additive variance of fitness will be associated with inbreeding depression and, therefore, it is unlikely to produce a net increase in the adaptive potential of populations. In addition, previous considerations were based on the theoretical analysis of the behavior of neutral genetic variation after bottlenecks, and the role of selection acting on epistatic systems controlling fitness has not been studied.

In this article we used analytical and simulation methods to investigate the contribution of epistatic systems to the change in the mean and the genetic components of variance of fitness during bottlenecking, due to the joint action of mutation, natural selection, and genetic drift (*MSD*). To develop a biologically reasonable model, we assumed that mutations show a distribution of homozygous and heterozygous effects close to those experimentally observed in *Drosophila melanogaster*, and we imposed different types of epistasis on this basic system. The pattern and strength of epistatic effects on fitness is largely unknown, but synergism between homozygous deleterious mutations at different loci has often been reported in Drosophila mutation-accumulation experiments (Mukai 1969; Ávila *et al*. 2006). Therefore, we studied the consequences of synergistic epistasis in pairs of loci by increasing the deleterious effect of the double homozygote above that expected from the deleterious effects of the homozygotes at both loci involved. However, to explore the consequences of bottlenecking in a multiple-peak adaptive surface, we also considered cases of antagonistic epistasis where, at each pair of loci, the fitness of the double homozygote for the deleterious alleles was larger than expected. Of course, other epistatic models could also be considered, including those showing higher-order interaction effects, but the severe shortage of relevant empirical data makes the choice highly subjective and, consequently, we restricted our analysis to the simplest case. On the other hand, our procedure has the practical advantage of allowing the definition of epistasis by the addition of a single parameter to those describing the properties of individual loci.

Our aim was to describe and analyze drift-induced changes in the components of the genetic variance of fitness, where neutral predictions will be reliable only during extreme and brief bottlenecks. For moderate bottleneck sizes or long-term inbreeding, it becomes necessary to consider the concurrent effects of natural selection both on the standing variation and on that arisen by new mutation. Moreover, the nature of the genetic variability of fitness in the base population, arisen by mutation and shaped by natural selection and drift, is critical for the assessment of the consequences of subsequent bottlenecks. For nonepistatic models, the genetic properties of the trait can be theoretically inferred from the pertinent mutational parameters and effective population sizes by assuming a balance between mutation, selection, and drift. This can be numerically achieved using diffusion theory, and reliable approximations can be easily calculated by analytical methods (García-Dorado 2007). Notwithstanding, the analytical study of the contribution of epistasis to the genetic properties of fitness at the *MSD* balance becomes particularly difficult and it must be complemented with computer simulation.

## MODEL AND METHODS

#### The neutral model:

We consider an extension of the model developed by López-Fanjul *et al*. (2002), with two neutral independent diallelic loci (*i* = 1, 2) at Hardy–Weinberg equilibrium segregating with frequencies *p _{i}* and

*q*(= 1 −

_{i}*p*). Genotypic values are given in Table 1, using the standard notation for fitness models, although the following analysis is concerned only with neutral predictions. At the

_{i}*i*th locus, the homozygous effect is

*s*and the degree of dominance in the absence of epistasis is

_{i}*h*(

_{i}*h*= 0, 0.5, or 1 for complete recessive, additive, or complete dominance gene action, respectively). This basic gene action can be viewed as that shown by single loci segregating against a fixed genetic background. Epistasis has been imposed on that basic system and it is represented by a factor

_{i}*k*affecting the genotypic value of the double homozygote for the allele decreasing the metric trait at each locus (

*k*< 1,

*k*> 1, or

*k*= 1 for antagonistic, synergistic, or no epistasis, respectively). At the

*i*th locus, the marginal average effect of gene substitution (α

*), the marginal genotypic value of the heterozygote (δ*

_{i}*, expressed as deviation from the midhomozygote value), and the marginal degree of dominance (γ*

_{i}*) are given by(1)(2)and(3)where*

_{i}*c =*(

*k*− 1) . Thus, epistasis modifies the basic properties of the

*i*th locus, as α

*, δ*

_{i}*, and γ*

_{i}*become dependent on the allele frequencies at the other locus (*

_{i}*q*). That is, they are contingent on the genetic background. For a given

_{j}*k*value, the basic (

*h*) and the marginal (γ

_{i}*) degrees of dominance become closer as*

_{i}*q*decreases, and γ

_{j}*approaches zero (complete recessive) as*

_{i}*k*and

*q*increase.

_{j}In an infinitely large panmictic population (base population), the mean (ancestral mean *M*) and the components of the genetic variance (additive, *V*_{A}, dominance, *V*_{D}, and epistatic, *V*_{I}, ancestral components of variance) for the two-locus system considered are given by(4)(5)(6)and(7)Thus, the epistatic component of variance *V*_{I} is independent of the degree of dominance of the loci involved.

Assume now that a number of replicates of size *N* are drawn at random from the base population, the *N* individuals in each replicate mating in panmixia to produce a large sample of offspring in which Hardy–Weinberg equilibrium is restored (ignoring transient linkage disequilibrium). As Equations 4–7 are polynomial functions of (*i* = 1, 2, *m* = 1–4), explicit expressions for their expected values at the Hardy–Weinberg equilibrium after *t* consecutive bottlenecks of *N* parents (*i.e*., derived mean *M _{t}* and variance components

*V*

_{At},

*V*

_{Dt}, and

*V*

_{It}) can be readily obtained by substituting by the corresponding exact

*m*th moment of the allelic frequency distribution with binomial sampling given by Crow and Kimura (1970, p. 335). As can also be expressed in terms of the inbreeding coefficient after

*t*generations (

*F*), the equations for the derived mean and variance components also apply when bottleneck sizes are not constant from generation to generation and, thus, a single parameter describes the outcome of an arbitrary bout of random drift. Notwithstanding, equations for the derived variance components can be managed only numerically.

_{t}The derived values of the marginal average effect of gene substitution at locus *i* (α* _{it}*) and the marginal genotypic value of the heterozygote (δ

*) can be obtained by taking averages in Equations 1 and 2, givingandThus, temporal changes in α*

_{it}*and δ*

_{i}*will occur only in the presence of epistasis (*

_{i}*c*≠ 0). Due to sampling, there will also be variance in α

*and δ*

_{i}*among replicated bottleneck lines at a given generation (descriptive of the spatial variation in a metapopulation structure), and taking variances in Equations 1 and 2 givesandThus, spatial changes in δ*

_{i}*will occur only in the presence of epistasis (*

_{i}*c*≠ 0) but those for α

*require only dominance (*

_{i}*h*≠ ), with or without additional epistasis.

_{i}From Equations 5 and 6, the derived additive and dominance components of variance after *t* bottlenecks can be expressed asandwhere *H _{i}* is the heterozygosity at the

*i*th locus (= 2

*p*). These equations explicitly show how the derived components of variance depend on the spatial and temporal changes in α

_{i}q_{i}*and δ*

_{i}*and also on a covariance factor. Although cov(δ*

_{i}

_{it}^{2},

*H*

_{i}^{2}) does exist only with linkage disequilibrium, cov(α

_{i}^{2},

*H*) depends on the type of basic gene action and requires only dominance (with or without epistasis).

_{i}Taking expectations in Equation 4, the change in mean after *t* bottlenecks is given by(8)where is the dominance × dominance epistatic component of variance (Crow and Kimura 1970, p. 80). Thus, Equation 8 indicates that antagonistic epistasis (*k* < 1) and/or basic dominant gene action (*h _{i}* > 0.5;

*i.e*., δ

*< 0) may imply an unrealistic enhancement of the mean with inbreeding.*

_{i}Finally, the between-line variance *V*_{B} after *t* bottlenecks can be obtained by taking variances in Equation 4. Again, the corresponding analytical expression is intractable but numerical solutions can be calculated for any combination of allele frequencies.

Considering the whole genome as consisting of a number of pairs of loci, each defined as in the above model, neutral predictions for the change in mean Δ*M _{t}* and the

*V*

_{At},

*V*

_{Dt},

*V*

_{It}, and

*V*

_{Bt}variance components after

*t*bottlenecks were obtained by adding up the corresponding contributions from each pair, ignoring linkage disequilibrium.

#### The simulation model:

We used simulation to study the consequences of mutation and selection in the bottleneck process described above. The trait studied was fitness, which is usually assigned to causes acting multiplicatively so that the trait value is always positive. Therefore, we simulated the variation in fitness due to segregation at several independent pairs of loci in a multiplicative way. Genotypic values for a pair of loci are given in Table 1, and those for individuals were obtained as the product of the corresponding contributions from each pair of loci involved. This ensures that, in the absence of epistasis (*k* = 1), relative deleterious effects at each locus are independent of the genetic background. For each base population or bottlenecked line, the contributions of each pair of loci to the within-line *V*_{A}, *V*_{D}, and *V*_{I} variance components were computed from expressions (5)–(7), and these were added up for all pairs of loci involved. Therefore, under neutrality (*D* scenario, see below), estimates from simulation data, computed in the real scale, were directly comparable to their corresponding analytical predictions. However, due to the multiplicative model used in the simulations, the temporal decay in mean fitness was smaller than that corresponding to the nonmultiplicative model, and the between-line variance for fitness became dependent on the fitness average. Mean fitnesses from simulation results are given in the real scale, since this provides a more intuitive picture from which satisfactory predictions can be obtained under the multiplicative model (García-Dorado 2008). Nevertheless, regarding the between-line variance, a logarithmic transformation was used to smooth the comparison between simulation and analytical results; *i.e*., simulation estimates were computed from log-transformed fitness values, and the corresponding predictions were calculated using mutational effects for log-scaled fitness, where detrimental effects are larger than those corresponding to the real scale, particularly in the case of severely deleterious alleles. Thus, for the *i*th neutral locus, predictions for the between-line variance were obtained usingand

##### Base population and mutational parameters:

Individual-based simulations were used to obtain a base population under mutation–selection–drift balance (*MSD*), by allowing the occurrence of nonrecurrent mutations in a population of 10^{3} individuals during 10^{4} generations. Potential parents for each individual were sampled with replacement from the corresponding parental population with a probability proportional to the parents' fitnesses. Thus, a parent was chosen if a random number was lower than its individual fitness and disregarded otherwise.

Mutations were assumed to occur at a rate λ per haploid genome and generation, and their homozygous effects were sampled from a gamma distribution *f*(*s*) with shape parameter β and mean effect . Three different models covering a wide range of parameters (first four columns in Table 2 and Figure 1) were used in the simulations. Model I assumed few mutations of large average effect. Model III considered a 10-fold mutation rate and accounted for about the same rate of severely deleterious mutations (say those with *s* > 0.3), but it also allowed for a higher rate of mutations with tiny to moderate effect (say those with *s* < 0.3), so that the average deleterious effect was much smaller. Model II represented an intermediate situation, similar to model I regarding the rate of moderate to severe deleterious mutation, but accounting for a larger rate of tiny to mild deleterious mutations. The shape parameters of the gamma distributions were chosen such that the mutational variance for each model was 0.002, a value roughly close to that experimentally observed for fitness in *D. melanogaster* (García-Dorado *et al*. 1999, 2004). The dominance coefficient of mutations, *h*, was obtained from a uniform distribution ranging between 0 and exp(−*Ks*). This model also has experimental support (Caballero and Keightley 1994; García-Dorado and Caballero 2000), and a value of *K* = 7.5 was used in all cases to maintain the constant relationship between and *s*, as suggested by García-Dorado (2003).

The total number of loci available for mutation (fifth column in Table 2) was set up to allow a proportion of segregating loci in the base population, which was ∼70–80% for each model (sixth column in Table 2), so that new mutations could always be assigned to nonsegregating sites as in a nonrecurrent model, and ∼55–60% of those loci were segregating in pairs (seventh column in Table 2). If a mutation got fixed, its effect in homozygosis was accumulated to the basal fitness and the corresponding locus was made available for new mutation to economize computational resources. For each combination of mutational models and epistatic effects we obtained a set of 50 base populations (replicates).

Scenarios were run (i) without epistasis (*k* = 1); (ii) with synergistic epistasis (*k* = 2 or 4), *i.e*., in the theoretical model the deleterious effect of the double homozygote was twofold or fourfold that without epistasis; or (iii) with antagonistic epistasis (*k* = 0 or 0.5), *i.e*., the double homozygote did not show a deleterious effect or this was half that without epistasis, respectively. Note that the case *k* = 0 produces a fitness surface with two equal adaptive peaks, corresponding to the fixation of the wild or the mutant alleles at both loci. Moreover, for the particular case where both loci had the same (*s*, *h*) values, *k* < *h* always produces a fitness surface with two adaptive peaks.

##### Simulated bottlenecks:

Starting from the base populations at the *MSD* balance, we simulated lines that were maintained with a number *N* of breeding individuals per generation (equivalent to consecutive bottlenecks) under two different scenarios. The first one (*D*) considered drift as the only agent of genetic change, to contrast the results with the analytical neutral predictions. The second scenario (*MSD*) allowed for natural selection and mutation during the bottleneck process. From each base population, 100 individuals were sampled without replacement to estimate the mean and the variance components at generation zero, thus simulating an experimental situation. After evaluation, a set of *N* parents was chosen at random in the *D* scenario or, as explained above, checked for survival in the *MSD* scenario. These parents were randomly paired, and 100 offspring were used to estimate the new mean and variance components. This process was repeated during 100 generations and was replicated 1000 times for each base population. Different bottleneck sizes of *N* = 2, 10, and 50 individuals were considered for each mutational model.

## RESULTS

#### The base populations:

In the base populations, the allele frequency distributions were “L”-shaped and, as expected, their medians and variances were both inversely related to the average strength of selection (*i.e*., they increased from model I to III, Figure 2). Within models, different *k* values did not result in appreciable differences between those distributions and only results for the nonepistatic case (*k* = 1) are shown in Figure 2.

Table 3 gives averages over replicates for the mean fitness in the base populations, the corresponding genotypic variance *V*_{G} and its components (*V*_{A}, *V*_{D}, and *V*_{I}, expressed as percentages of *V*_{G}), and the gametic disequilibrium averaged over all pairs of segregating loci, which was measured by *D* [the difference *g*(*A*_{1}*B*_{1})*g*(*A*_{2}*B*_{2}) *− g*(*A*_{1}*B*_{2})*g*(*A*_{2}*B*_{1}), where *g* stands for gametic frequency] or *Z* [the ratio *g*(*A*_{1}*B*_{1})*g*(*A*_{2}*B*_{2})/*g*(*A*_{1}*B*_{2})*g*(*A*_{2}*B*_{1})]. In the absence of epistasis, the expected mean at the *MSD* balance can be computed using diffusion theory, where fitness is relative to that of a genotype carrying none of the segregating deleterious alleles, giving 0.916, 0.842, and 0.465 for mutational models I, II, and III, respectively (see computing procedure in García-Dorado 2003). These values were in good agreement with the corresponding mean fitnesses of the simulated base populations relative to the genotype carrying none of the segregating deleterious alleles. They also were in reasonable agreement with their deterministic expectations for an infinite population (exp[−2λ], *i.e*., 0.905, 0.819, and 0.368 for mutational models I, II, and III, respectively), implying that natural selection was the main force shaping the standing fitness genetic variance. Table 3 also shows the mean fitness relative to a genotype carrying none of the deleterious alleles occurred during the 10^{4} generations, indicating that only the mutational model III produced an appreciable load from deleterious fixation.

The genetic variance for fitness was small (averaging 6.4 × 10^{−3} for models I and II and 12.5 × 10^{−3} for model III), the fraction corresponding to *V*_{D} being about double that for *V*_{A} for models I and II (*V*_{D}/*V*_{G} = 65.4% *vs. V*_{A}/*V*_{G} = 33.73%, on average) and the reverse for model III (*V*_{D}/*V*_{G} = 35.0% *vs. V*_{A}/*V*_{G} = 63.3%, on average). On the whole, a considerably lower mean fitness as well as a greater *V*_{G} and *V*_{A}/*V*_{G} ratio were obtained for model III, characterized by the highest rate of tiny to mild deleterious mutation. Epistasis, either synergistic or antagonistic, had little influence on the above results. For all models and epistatic systems, the *V*_{I}/*V*_{G} ratio obviously increased with the strength of epistasis, although its value never exceeded 5% and it was unrelated to those of the *V*_{A}/*V*_{G} and *V*_{D}/*V*_{G} ratios (as shown by Equation 7).

Linkage disequilibrium, as measured by *D*, was negligible (of the order of 10^{−6}), and using a sequential Bonferroni test, the only *D* value significantly different from zero was negative, as expected from the action of natural selection since it corresponded to *k* > 1. To study the consequences of selection in a deterministic model, the *Z* measure of disequilibrium can be more informative since, for deleterious allele frequencies approaching zero, it remains at its quasi-equilibrium value (*Z* = 1, *Z* < 1, or *Z* > 1 for a nonepistatic, synergistic, or antagonistic multiplicative fitness model, respectively). However, random fluctuation in gametic frequencies from mutation and drift will increase the average *Z* value over its deterministic expectation. On the other hand, linkage disequilibrium can be underestimated, since *Z* cannot be computed in the absence of a repulsion gamete, and a logarithmic transformation will be of no use in this case. In agreement with this, the *Z* values reported in Table 3 had large standard errors (ranging from 0.02 to 0.25) and oscillated erratically around 1 (overall *Z* = 1.03). Only two *Z* values significantly departed from 1 after using a sequential Bonferroni test, but they did so in the opposite direction to that expected from the action of natural selection, evidencing that random processes and sampling bias were more important than natural selection in determining the *Z* values.

#### The consequences of bottlenecking: synergistic *vs*. nonepistatic systems:

Estimates of the additive variance, averaged over replicates, are shown in Table 4, both for the base population (*N* = 1000) and after 100 consecutive bottlenecks (*N* = 2, 10, 50). Without epistasis (*k* = 1), expectations at the *MSD* balance computed from diffusion theory and from a shortcut method (García-Dorado 2007) are also given in Table 4 and they were in satisfactory agreement with each other for *N* ≥ 10. For mild bottlenecks (*N* = 50), the final inbreeding coefficient was *F* = 0.63 and the population should still be approaching the corresponding *MSD* balance. In this situation, the derived variance after bottlenecks was larger than expected at the new equilibrium, so that some further reduction would follow as *F* approaches its final value of 1. For strong bottlenecks (*N* = 2, 10), however, the final inbreeding coefficient was *F* ≈ 1 and, therefore, a new *MSD* balance should have been attained at this stage. In agreement with this, *MSD* expectations fitted the pertinent simulation results (*k* = 1) for the base population (*N* = 1000) and, after 100 generations, for *N* = 10. The fit was still acceptable for *N* = 2 although, in this particular case, a continuous theoretical treatment of the process is scarcely appropriate. In general, the effect of epistasis on the value of the *MSD* additive variance obtained from simulation was significant only for strong epistasis (*k* = 4) and extreme bottlenecks (*N* = 2).

The decay in mean fitness after *t* successive bottlenecks (up to 100) of size *N* (2, 10, or 50) is shown in Figure 3 for the three mutational models, the two scenarios, and the different types of gene action. In the *D* scenario, simulation results fitted very precisely the analytical neutral predictions, as shown in Table 5 where only values for model I and *N* = 50 are given. The decay was larger for model III, characterized by the highest rate of mildly deleterious mutations, and it was practically unaffected by epistasis. In the *MSD* scenario, the decay was even less affected by epistasis, so that, as shown in Figure 3, it reasonably agreed with the analytical predictions for the *MSD* nonepistatic case (García-Dorado 2008). For *N* = 2, these predictions were slightly above the true mean in the long term, due to the bias introduced by assuming additive gene action in the prediction of the deleterious fixation rate at the new *MSD* balance. In some cases, fitness showed a slight fitness rebound after an initial inbreeding depression. This was appreciable only for models I and II with *N* = 10 and cannot be accounted for by García-Dorado's approach. For *N* ≥ 10, the *MSD* decline (both predicted and simulated) was much smaller than that obtained in the *D* scenario, an outcome that must be ascribed to natural selection (partially) purging recessive deleterious alleles. For *N* = 2, however, selection was quite inefficient, implying early fitness declines almost as large as those obtained in the *D* scenario and substantial later declines from new deleterious mutation.

Figure 4 shows the additive (*V*_{A}), dominance (*V*_{D}), and epistatic (*V*_{I}) components of the genetic variance after *t* consecutive bottlenecks of size *N* = 2, 10, 50, plotted against the corresponding inbreeding coefficient *F*. In each model, variance components are given for the *D* and *MSD* scenarios and the three types of gene action considered per scenario. In the *D* scenario, the magnitude of each variance component for given *F* and *k* values was of course the same for all bottleneck sizes (black lines), and the simulation results fitted very precisely the corresponding analytical predictions for the additive variance (Table 5). In the *MSD* scenario, however, the magnitude of a variance component for the same *F* and *k* values depends on the pertinent bottleneck size. In this case, predictions are available only for the additive variance contributed by nonepistatic loci (García-Dorado 2007) but they are not shown since they underestimate the transient increase in additive variance after bottlenecks when the degree of dominance was very small (García-Dorado 2008).

In all instances, the additive, dominance, and epistatic components of the genetic variance behaved in a similar manner, initially increasing with *F* up to a maximum value and subsequently declining to zero in the *D* scenario or to an equilibrium value in the *MSD* scenario, where continued mutation prevented the final exhaustion of those variance components (in the case of the epistatic variance, the equilibrium value was so small that it cannot be appreciated in Figure 4). In the *MSD* scenario, selection eroded the variance components after bottlenecks with respect to their corresponding values in the *D* scenario, slightly for *N* = 2 but substantially for *N* ≥ 10. In general, the highest peaks were reached in the *D* scenario (or in the *MSD* scenario for *N* = 2), at a value of *F* ∼ 0.5 (*V*_{A}), 0.3 (*V*_{D}), or 0.4 (*V*_{I}), while peaks for *N* ≥ 10 in the *MSD* scenario were much lower and were attained at smaller *F* values as *N* increased. For all models and scenarios, only strong synergistic epistasis induced a significant departure from the pattern of change observed for the additive and dominance variance without epistasis. This departure consisted of an increase in the pertinent variance value that, in the *MSD* scenario, became negligible for *N* ≥ 10, particularly in the case of the additive component.

Table 5 shows that the value of the between-line variance *V*_{B} observed in the *D* scenario for the log-scaled trait was in good agreement with the corresponding neutral predictions computed from appropriately transformed effects (see model and methods). The evolution of the between-line differentiation after bottlenecks (log-scaled), plotted against *F*, is shown in Figure 5 for *N* = 10. As expected from theory, the increase in the between-line variance in the *D* scenario accelerated with inbreeding, at a rate that was larger for epistatic systems than for simple dominants, particularly for model III and strong epistasis. On the contrary, the increase in *V*_{B} in the *MSD* scenario was very small and for models I and II, where most deleterious mutations had important effects that are efficiently selected against, the corresponding plot visually overlaps with the abscissa axis in Figure 5. Summarizing, uniform purifying selection impeded differentiation among the lines for all *k* values, despite the continuous input of new nonrecurrent mutations.

The change in the variance of the additive variances of the lines *V*(*V*_{A}) after bottlenecks plotted against the pertinent *F* values is shown in Figure 6 for model I, the *D* and *MSD* (*N* = 10) scenarios, and all *k* values. In the *D* scenario, *V*(*V*_{A}) initially increased until an *F* value ∼0.6–0.7 was reached and subsequently declined to zero, the excess for strong synergistic epistasis being much larger than those for moderate epistasis or simple dominance. In the *MSD* scenario, however, the transient increase in *V*(*V*_{A}) for strong epistasis was much smaller and it was practically negligible in the remaining cases.

#### The consequences of bottlenecking: antagonistic *vs*. nonepistatic systems:

Given the qualitative agreement between the results obtained from the three mutational models reported above, only the intermediate model II was explored for the effects of antagonistic epistasis (*k* = 0 or 0.5) and, as the results for the two *k* values considered were very similar, only those for *k* = 0 are given. The evolution of the overall mean, the additive variance, and the between-line variance after bottlenecks is presented in Figure 7, both for the epistatic (*k* = 0) and the nonepistatic case (*k* = 1). In the *D* scenario, the behavior of the above parameters was very close for those *k* values and also for synergistic epistasis (*k* = 2), although the mean fitness depression with antagonistic (synergistic) epistasis was slightly smaller (larger) than that in the absence of epistasis, and the opposite was observed for the increases in the components of the genetic variance and the between-line variance. In the *MSD* scenario, the discrepancies between the results obtained with (*k* = 0) or without epistasis (*k* = 1) were even smaller than in the *D* scenario. Thus, in practice, the efficiency of natural selection to promote new favorable epistatic combinations was not enhanced by bottlenecking. Furthermore, selection prevented between-line differentiation irrespective of the presence of epistasis. Therefore, the antagonistic epistasis systems studied did not qualitatively modify the evolutionary pattern for synergistic epistasis described above, in spite of the two-peak adaptive surface produced by the former system.

## DISCUSSION

#### The genetic model:

Computer simulations were carried out to produce a biologically plausible representation of the genetic variation for fitness in a population at the *MSD* balance (base population) with or without epistasis. To do that, we used three different models to describe the deleterious properties of spontaneous mutation covering the typical range of empirical estimates (García-Dorado *et al*. 1999, 2004). Thus, all models accounted for similar inputs of mutational variance and mutation rates for moderate to severe deleterious effects, and they also assumed the same negative exponential relationship between the expected degree of dominance and the homozygous deleterious effect, but they differed in the mutation rates for tiny to mild deleterious effects.

Epistasis embraces a huge variety of types but reports of significant epistatic effects on metric traits usually come from crosses between highly divergent lines, these results being scarcely relevant to the description of the genetic architecture of outbred populations. For morphological traits, epistasis may be ascribed to interactions between the average effects of quasi-additive loci, which will mainly contribute to the additive × additive epistatic component of variance. For fitness, however, epistasis is generally assumed to be pervasive, although the type and magnitude of the interactions between the nonadditive genes involved have been poorly explored (see reviews by Barton and Keightley 2002 and Carlborg and Haley 2004). Notwithstanding, synergism between the deleterious effects of homozygous mutations at different loci has often been reported in *D. melanogaster* (Mukai 1969; Whitlock and Bourguet 2000; Ávila *et al*. 2006) and *Saccharomyces cerevisiae* (Dickinson 2008). Furthermore, a recent analysis suggests that the average degree of epistasis for fitness correlates to genomic complexity, interactions tending to be antagonistic in unicellular eukaryotes and synergistic in higher eukaryotes (Sanjuán and Elena 2006). Thus, of the countless ways in which epistasis for fitness could be modeled, we chose a simple one where, for each pair of loci, the genotypic effect of the double homozygote for the alleles that decrease fitness was increased (synergistic epistasis) or diminished (antagonistic epistasis), therefore contributing to all possible epistatic components of variance (additive × additive, additive × dominant, and dominant × dominant). Of course, higher-order interactions could also be considered but they will not be prevalent, as blocks composed of large numbers of epistatic loci are unlikely to segregate at the *MSD* balance. Even in our case, where the genome consisted of pairs of epistatic loci, only about half of the segregating loci at the base population showed epistasis (Table 2). Notwithstanding, a substantial proportion of the genome still remains involved in epistatic interactions, thus providing conservative results with respect to the impact of epistasis on the change of the mean and the genetic variance after consecutive bottlenecks.

The genetic variation of fitness at the base population described above was subsequently exposed to the joint action of continued mutation, selection, and drift (*MSD* scenario) or to drift alone (*D* scenario). The pertinent changes in the mean and the components of the within-line genetic variance, as well as those in the between-line variance, were compared among scenarios. In the neutral case, a close fit was found between the simulation results and those computed by analytical methods for all types of gene action considered, but natural selection completely transformed the whole picture.

#### The base populations:

For mutational models I and II, a close agreement was found between the mean fitness of the base populations in the absence of epistasis, relative to that of the original deleterious-free genotype, and their expectation at the *MSD* balance computed from diffusion theory. This implies that deleterious fixation caused no substantial load after 10^{4} generations, when the base population was evaluated. Furthermore, the observed mean fitness values were in rough agreement with their deterministic predictions for an infinite population, implying that, for *N* ≥ 10^{3}, the main factor limiting the segregation load was the deleterious mutation rate, rather than the form of the distribution of deleterious effects or the type of gene action. Thus, for those mutational models, natural selection was the main force shaping the genetic architecture of fitness in populations of effective size ≥10^{3}. For model III, however, with a larger rate of slightly deleterious mutations, the load from deleterious fixation was appreciable, although small, and the segregating load was somewhat smaller than its deterministic prediction, indicating that the role of genetic drift in determining the genetic architecture of fitness was a mild one. In addition, the population mean fitness was virtually unaffected by epistasis. This is to be expected, due to the scarce influence of epistasis on the ancestral fitness additive variance, as is discussed below.

With loose linkage, Fisher's theorem of natural selection will hold at the quasi-linkage equilibrium even in the presence of epistasis, because the fitness change due to natural selection acting upon the additive × additive component of the epistatic variance will be compensated by that caused by the erosion of linkage disequilibrium through recombination (Crow and Kimura 1970, p. 217). Therefore, the rate of fitness change will be equal to the genic variance for fitness, which includes the effects of gametic disequilibrium (Kimura 1965). Strictly, these considerations apply to large populations where the evolution of gene frequencies is driven by natural selection and the consequences of new mutation are disregarded. Our results show that, in populations at the *MSD* balance, the epistatic component of the genetic variance for fitness was negligible, so that the gametic disequilibrium induced by epistasis between unlinked deleterious mutations is expected to be irrelevant. In fact, when linkage disequilibrium was computed as *Z* (the parameter to which the quasi-equilibrium predictions refer), the values obtained were more affected by new mutation and drift than by selection acting on epistatic effects. In addition, we carried out some simulations where, after running the base population, all loci were randomized before bottlenecking to eliminate gametic disequilibrium, but we found no change in the results.

Epistasis did not appreciably constrain the frequency distribution or the number of segregating deleterious alleles. This is due to those frequencies being very low at the *MSD* balance, making it very unlikely that pertinent mutations in each pair of loci will be present in the same individual and, therefore, their epistatic effect will not be expressed. As a consequence, the magnitude of the additive variance was also practically independent of epistasis. It should be noted that this additive variance is that corresponding to a panmictic population at linkage equilibrium, obtained by adding up the theoretical contributions of all pairs of loci computed from Equations 1 and 5. Nevertheless, since linkage disequilibrium became irrelevant, this additive variance should be practically equal both to the genic variance for fitness and to the rate of fitness increase from natural selection. In general, the situation was characterized by a small genetic variance with a relatively important additive component, a description that might seem inconsistent with most of the genetic variance being due to the segregation of deleterious (partially) recessive mutations at low frequencies, which will mostly contribute to the dominance variance. However, this result is in qualitative agreement both with the empirical observations and with the theoretical expectations. The reason is that, at the *MSD* balance, the frequency distribution of alleles with different deleterious effects and dominance coefficients is spread by drift, so that those loci generating the smallest *V*_{A}/*V*_{G} ratios (those segregating for rare severely deleterious alleles with important recessive action) will make a minor contribution to the overall fitness genetic variance, as discussed by Hill *et al*. (2008).

#### Temporal changes after continued bottlenecking:

##### Changes in mean fitness:

In the *D* scenario, a dramatic decay in mean fitness was observed for the three mutational models, but the effect of epistasis was very slight unless epistasis was strong (*k* = 4) and the rate of tiny to mild deleterious mutation was high (model III). On the contrary, in the *MSD* scenario, substantial inbreeding depression was attained only after extreme bottlenecks, which also produced an additional decline in mean from the accumulation of new deleterious mutations that became important in the long term. On the whole, the simulation showed that the selective purge of partially or totally recessive deleterious alleles due to their increased effect during inbreeding was very efficient, except in the case of extreme bottlenecks. It also showed that the synergistic effects of epistatic deleterious alleles segregating in the base population would scarcely contribute to depression even in the absence of selection. In all cases where natural selection efficiently purged inbreeding depression it also removed any further decline in mean due to synergy. However, in the case of antagonistic epistasis, selection was unable to promote the transition to the alternative adaptive peak corresponding to fixation of the initially deleterious mutations at the loci involved.

##### Changes in the components of the within-line genetic variance:

In both scenarios, all components of the genetic variance increased with *F* until a maximum value was reached, subsequently decreasing to zero (*D* scenario) or to a final equilibrium value (*MSD* scenario). Purging selection, however, eliminated a large fraction of the excess in variance components, and the remaining part of that excess was attained at lower *F* values. Only for extreme bottlenecks (*N* = 2) drift dominated over selection so that similar results were obtained in both scenarios. This is relevant to most experimental results that refer to lines that passed through one generation of a bottleneck of two individuals, where a considerable increase in the within-line additive variance together with a large inbreeding depression was detected for *D. melanogaster* or *Tribolium castaneum* viability (see reviews by López-Fanjul *et al*. 2003, Zhang *et al*. 2004, and Van Buskirk and Willi 2006). In Drosophila, that increase continued after three bottlenecks (up to *F* = 0.50) and was followed by a decrease after six bottlenecks (*F* = 0.73), also in agreement with theory (García *et al*. 1994). However, from the evolutionary point of view, such extreme and continued bottlenecks are unlikely to be relevant. Our results illustrate that bottleneck-induced excesses in additive variance are associated with drastic inbreeding depression. Therefore, after population expansion, natural selection acting upon such excesses will remove just the fraction of the inbreeding depression due to deleterious alleles still segregating at the population. Obviously, for traits that are not strongly related to fitness, no association between increased additive variance and inbreeding depression is expected. In this instance, any excess in additive variance after bottlenecking could be used by natural selection if the trait becomes adaptive, but it should be noted that such traits will not usually show the nonadditive (dominance) gene action required to produce an excess in additive variance with inbreeding.

In the particular case of multiple neutral loci showing additive × additive epistasis at different levels, it has been analytically shown that an increase in the additive component of variance after bottlenecks is invariably associated with a parallel erosion of the epistatic components of variance as inbreeding progresses, starting from the reduction of that corresponding to the highest-order interaction and subsequently followed by the decrease of the lower-order components. This result, first obtained by Goodnight (1988) for two-locus systems and later extended to multiple-locus systems by Barton and Turelli (2004) and López-Fanjul *et al*. (2006), suggested a conversion process where the additive variance inflates over time at the expense of a loss in nonadditive variance. However, for neutral nonadditive loci (with or without epistasis), it has been theoretically shown that an excess in additive variance after bottlenecks will not be generally compensated by a corresponding decrease in nonadditive variance (López-Fanjul *et al*. 2002; Barton and Turelli 2004), except in specific epistatic situations such as those involving marginal overdominance (underdominance) for basic loci segregating at equilibrium (intermediate) frequencies. In agreement with this prediction, our results clearly show that the values of all variance components were inflated after bottlenecks, both under drift alone and when selection was acting.

On the whole, theoretical analyses indicate that dominance can be generally considered as the primary cause of an increase in the components of the additive variance after bottlenecks, this being merely modulated by epistasis (López-Fanjul *et al*. 2002; Barton and Turelli 2004). Only for specific models involving high levels of multilocus interaction, the role of dominance was subordinated to that of epistasis (Naciri-Graven and Goudet 2003) but, as discussed by Turelli and Barton (2006), these models are biologically implausible as, contrary to empirical observations, they imply a large contribution of the epistatic variance to the total genetic variance. Our simulation results for simple dominant nonepistatic loci indicate that the changes in the additive and dominance components of the genetic variance after bottlenecks, both in the *D* and in the *MSD* scenarios, were practically indistinguishable from those obtained with additional moderate epistasis (synergistic or antagonistic) and significant departures from this pattern were observed only in the case of strong epistasis, although they were generally small. As expected when most genetic variance can be ascribed to the segregation of rare (partially) recessive alleles, the values of the additive and dominance components of variance after an equal number of bottlenecks were of the same order, with *V*_{A} *> V*_{D} in the *D* scenario and *V*_{A} *< V*_{D} in the *MSD* scenario. Furthermore, the value of the epistatic component of variance was always very small, about one order of magnitude below those of the additive or the dominant components in the case of strong epistasis and much lower for moderate epistasis.

The effect of the mutational model used was generally small unless the level of epistasis was strong. However, higher mutation rates and smaller average deleterious effects (model III) resulted in larger changes in the magnitude of the variance components after bottlenecks. After many bottlenecks, the values attained with dominance or moderate epistasis were quite similar and significant departures were detected only in those cases combining strong epistasis (*k* = 4), extreme bottlenecks (*N* = 2), and the highest mutation rate for tiny to mild deleterious effects (model III).

#### Spatial changes after continued bottlenecking:

At low levels of inbreeding, it has been shown that the behavior of the between-line variance due to neutral loci is not greatly affected by their type of gene action (López-Fanjul *et al*. 2002; Barton and Turelli 2004). For high levels of inbreeding, however, the between-line variance contributed by rare recessives will greatly exceed the additive expectation and this excess will be larger in the presence of epistasis. These theoretical predictions have also been substantiated by our simulations where, in the *D* scenario, between-line differentiation in mean fitness showed an accelerated increase with *F* that can be intensified (diminished) with the strength of synergistic (antagonistic) epistasis. In the *MSD* scenario, however, the between-line variance practically vanished for all types of gene action considered.

After bottlenecks, most of the spatial variation of the genetic variance generated by neutral additive loci among replicated lines [*V*(*V*_{A})] is expected to be due to the buildup of transient linkage disequilibrium between pairs of loci caused by sampling (Avery and Hill 1977). However, for neutral nonadditive loci, our simulation estimates of the variance components, which do not include the effect of gametic disequilibrium as they were computed from allele frequencies, showed that *V*(*V*_{A}) increased up to a maximum value that was attained at *F* ∼ 0.6–0.7 and subsequently declined to zero, as *V*_{A} dissipated in all replicates as fixation progressed. With strong epistasis, the magnitude of this temporal increment in *V*(*V*_{A}) induced by drift was much higher than that observed with dominance only. With selection, however, the values of *V*(*V*_{A}) were practically indistinguishable from zero unless strong epistasis was considered. Thus, in both scenarios, drift-induced changes in *V*_{A} and *V*(*V*_{A}) showed similar patterns.

The effect of epistasis on diversification in an adaptive surface was also considered. This was achieved by studying a model with antagonistic epistasis that determines two adaptive peaks of equal height, corresponding to the double homozygotes for the wild or the mutant allele at both loci, respectively. In this case, the increases in additive variance and between-line variance after bottlenecks were very small, unless extreme bottlenecking drastically reduced the effectiveness of selection. Therefore, natural selection was unable to push the population to the new peak. Of course, a peak shift might occur if the population happens to drift toward the new peak during extreme bottlenecking, followed by an expansion allowing for efficient selection, but the probability of this case should be small. The situation may be different if combinations of new antagonistic mutations produce a new adaptive peak higher than the original one, as selection will be stronger in the neighborhood of this peak. Although this pattern of gene action is probably very uncommon, it could be worthwhile to explore its consequences under brief bottlenecks followed by expansion and to assay how likely is natural selection to promote genetic diversification by driving populations to different adaptive peaks. Genetic diversification from spatially varying selective pressures, implying genotype × environment interactions for fitness, could also be affected by epistasis, but this exceeds the scope of our work that is concerned only with uniform and constant natural selection.

Summarizing, drift, acting as the only evolutionary agent, can induce excesses both in the additive variance of a trait and in its between-line differentiation in mean and additive variance, particularly when strong epistasis is present. For fitness, however, the deleterious properties of spontaneous mutations primarily imply purifying selection. Our results show that the balance between drift and purifying selection is such that the actual increase in additive variance was very limited and the between-line differentiation was efficiently halted. Thus, nonadditive gene action, including epistasis, does not usually increase evolutionary rates, even for systems generating two-peak adaptive surfaces.

## Acknowledgments

We thank J. M. Alvarez-Castro for useful discussions and Charles Goodnight and two other anonymous referees for helpful comments. A.P.-F. is currently funded by an Ángeles Alvariño research fellowship from Xunta de Galicia. This work was supported by grants CPE03-004-C2 (Plan Estratégico del Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria, Spain), CGL2006-13445-C02/BOS, and CGL2008-02343/BOS (Ministerio de Educación y Ciencia, Spain, and Fondos Feder) and by Xunta de Galicia.

## Footnotes

Communicating editor: J. Wakeley

- Received May 8, 2009.
- Accepted July 5, 2009.

- Copyright © 2009 by the Genetics Society of America