Evolution of adaptive phenotypic variation patterns by direct selection for evolvability

A basic assumption of the Darwinian theory of evolution is that heritable variation arises randomly. In this context, randomness means that mutations arise irrespective of the current adaptive needs imposed by the environment. It is broadly accepted, however, that phenotypic variation is not uniformly distributed among phenotypic traits, some traits tend to covary, while others vary independently, and again others barely vary at all. Furthermore, it is well established that patterns of trait variation differ among species. Specifically, traits that serve different functions tend to be less correlated, as for instance forelimbs and hind limbs in bats and humans, compared with the limbs of quadrupedal mammals. Recently, a novel class of genetic elements has been identified in mouse gene-mapping studies that modify correlations among quantitative traits. These loci are called relationship loci, or relationship Quantitative Trait Loci (rQTL), and affect trait correlations by changing the expression of the existing genetic variation through gene interaction. Here, we present a population genetic model of how natural selection acts on rQTL. Contrary to the usual neo-Darwinian theory, in this model, new heritable phenotypic variation is produced along the selected dimension in response to directional selection. The results predict that selection on rQTL leads to higher correlations among traits that are simultaneously under directional selection. On the other hand, traits that are not simultaneously under directional selection are predicted to evolve lower correlations. These results and the previously demonstrated existence of rQTL variation, show a mechanism by which natural selection can directly enhance the evolvability of complex organisms along lines of adaptive change.


INTRODUCTION
It has long been recognized that phenotypic variation is highly structured [1 -4]. The best documented pattern is called morphological integration as identified by Olson & Miller [5], where functionally related traits tend to be more strongly correlated than functionally less-related traits. Morphological integration has been documented in many systems, from plants to vertebrate animals (e.g. [6][7][8][9][10]) and also confirmed at the level of quantitative trait locus (QTL) effects rather than phenotypic or genetic correlations (e.g. [11 -14]). Furthermore, it has been shown that patterns of integration change in evolution in a predictable way. For example, in quadrupedal animals, the corresponding elements in the forelimbs and hindlimbs are more correlated than in animals with specialized limbs. For instance, in bats, the forelimbs and hindlimbs perform different functions and have lower correlations than those of other mammals [15]. Moreover, the lower correlation among forelimbs and hindlimbs in hominoid primates is associated with a greater disparity in forelimb and hindlimb proportions than in monkeys [16]. All these data suggest that correlation patterns among traits evolve in a predictable way and influence the evolvability of the affected traits. These observations raise the question of what the evolutionary mechanisms that shape the pattern of phenotypic variation are. How much is genetic variation affecting the covariation among phenotypic characters? How does natural selection mould these variation patterns? The first question has been answered by QTL mapping and mutation accumulation studies which show that there is ample genetic variation affecting the covariation among characters (see below). In this paper we address the second question by modelling how selection acts on that form of genetic variation.
A variety of methods have been used to determine the amount of genetic variation affecting trait correlations. For instance, a recent quantitative genetic study compared correlations between traits that vary across families in a population [17]. The results show that a substantial proportion (about 40%) of this variation is genetic. Evidence for such variation has also been found in previous experimental work. For example, Lenski [18,19] documented genetic variation in Escherichia coli of the association between resistances to two different phages. Templeton and colleagues [20,21] found that the expression and association of the name-giving phenotype abnormal abdomen with other phenotypic traits that are affected by the same locus, varies across fly populations. Recently, Kim et al. [22] documented genetic variation in correlation of phenotypic traits in yeast, causing varying degrees of rust coloration being associated with different degrees of drug resistance.
To map the loci associated with genetic variation in trait correlation, Cheverud et al. [23] developed a QTL mapping method. By associating marker genotypes with the regression coefficient between two phenotypic traits, it was found that there are many polymorphic loci that affect the statistical relationships among phenotypic traits ( [17,23,24]; figure 1b). These loci have been called 'relationship Quantitative Trait Loci' (rQTL), because they change the variational relationships between characters. Genotypes can differ in the degree of canalization of variance both owing to gene-by-gene and gene-by-environment interaction. Two recent papers ( [17], see also [24]) established that the effect of the rQTL is because of differential epistatic interactions, i.e. epistatic interactions are trait-specific. Epistasis was also found responsible for variable intertrait associations in studies of single loci [18 -22]. In brief, the change in genetic background affects the traits of a pleiotropic domain of a locus in different ways, leading to differences in correlations among traits across genotypes at the focal locus (figure 1b). When multiple traits are affected by epistasis simultaneously, differential canalization of traits changes intertrait correlations through differentially modifying the expression of pre-existing genetic variation. Interestingly, rQTL mapping studies have shown that a substantial proportion of these rQTL have no detectable additive effect. These 'neutral' QTL (i.e. neutral with respect to the trait mean) constitute 30 -84% of the rQTL loci detected in mice (75% in Cheverud et al. [23]; 30% in Pavlicev et al. [17]; 84% in Leamy et al. [24]). This fraction of 'neutral' rQTL is particularly interesting because selection can act on their variational effects without conflicting with selection on their additive effects. Here, we present a phenomenological model of selection on these neutral rQTL, i.e. genes that change the variance and covariance among traits but do not directly affect the mean phenotype (see scheme in figure 1b). We show that this model can lead to the alignment of high levels of phenotypic variation with the direction of selection (figure 1c).

THE MODEL
To investigate how directional selection affects rQTL variation, we construct a simple haploid model of a population with two alternative genotypes at an rQT locus. We model the life cycle in three steps: in the first step, we model the evolution of the quantitative traits within each rQTL genotype class using Lande's phenotypic models [25,26]. Because each rQTL allele causes a different genetic covariance matrix, the selection response is different in the two rQTL genotype classes. This will lead to different average fitness of the two rQTL genotypes. In the second step, we model selection among the two rQTL classes using a deterministic one-locus, twoallele model. Finally, we model sexual reproduction as an incomplete mixture process between the genotypic values associated with each rQTL allele.
Specifically we consider a phenotype with two traits z 1 , z 2 , with sexual reproduction and a haploid genotype.
Modelling the haploid case aids transparency. However, the same scenario can be modelled for a diploid genome, although with considerably greater bookkeeping costs. The rQTL has two alternative alleles, A i and A j . We assume random mating, absence of gene-by-environment interaction and a constant G within each rQTL genotype class. The two alternative genotypes at the haploid rQTL differ in genotypic (G i and G j ) and phenotypic (P i and P j ) covariance matrices, but initially not in their average phenotypic values (figure 1b). The mean fitness values of the two rQTL genotype classes are denoted as w 1; ; w 2 and rQ TL ge no typ es t r a i t the difference between them determines the selection coefficient of the rQTL alleles, A i and A j .
The genotypic values of the next generation are then derived in three steps. In the first step, the trait means of the selected population are derived according to Lande's D z ¼ Gr ln w. The model is based on the selection regime (see below) and the respective genetic covariance matrices, with In our case with two traits, we write for rQTL allele A i : In most of our scenarios we assume constant directional selection on both traits, with the fitness function w(z 1 , z 2 ) ¼ w 0 exp[a 1 z 1 þ a 2 z 2 ], where w 0 is an arbitrary scaling coefficient [27]. The coefficients a 1 and a 2 determine the strength of directional selection on traits z 1 and z 2 , respectively. We use the multiplicative form of selection because a linear selection function causes a decrease in selection gradient with increasing population mean, while the selection gradient remains constant with multiplicative fitness. Using the exponential fitness function w(z) ¼ w 0 exp[a 1 z 1 þ a 2 z 2 ] and the bivariate normal distribution for the distribution of phenotypic values, the mean fitness reads as Differentiating the log mean fitness with respect to mean values gives respectively. Thus, following the selection of the breeding population, the change in traits affected by the two genotypes is In the second step, the fitness differences between the rQTL classes are used to determine the frequency change of the rQTL alleles according to Wright's selection equation: for A i with the frequency p, and analogously for A j with frequency q.
In the third step, the recombination between the rQTL and the genetic background affecting trait means is implemented as a mixture process of genotypes. This step reflects recombination during the transition to the next generation. Thus, the mean phenotypic value of the genotype i in the next generation (z 0 i ) depends on the phenotypic value of both genotypes in the previous generation (z S i ,z S j ), to the extent of recombination: where p, q are the allele frequencies of the two alleles at rQTL, and L is a measure of linkage between the rQTL and other loci affecting the mean, L [ [0.75, 1], with L ¼ 1 representing full linkage and L ¼ 0.75, free recombination. Linkage acts oppositely to recombination, maintaining the linkage disequilibrium (LD) between the advantageous alleles that build up owing to selection. Thus, instead of modelling recombination between pairs of loci explicitly, we model its composite effect on the phenotypic value of the particular rQTL genotype using the term L (see appendix A for derivation and an alternative method). The rationale for the range of linkage values is the following: if L would be 0.5, each mating among genotypes with different rQTL would lead to offspring with complete linkage equilibrium. Since free recombination at most leads to one-half of complete linkage equilibrium, the lowest level of linkage in this model is given by the mean between L ¼ 0.5 and L ¼ 1.
The resulting genotypic trait means of the next generation are used to calculate the corresponding mean fitness of the two rQTL genotypes, w 1 ; w 2 , in the next generation: Implementing recombination with the mean genotypic values generates trait means in the next generation. Owing to the different genetic covariances in the two rQTL genotype classes, the traits respond differently to selection within the rQTL classes. The genotype class with more genetic variation in the direction of selection will have a greater selection response and thus also acquire higher mean fitness than the alternative genotype class. Finally, the fitness differences between rQTL genotype classes determine the change in frequency of the rQTL alleles. The difference in mean fitness of the two rQTL genotype classes can be calculated as: where D w is the fitness difference of rQTL genotypes in the previous generation, and g kli are the elements of the genetic variance -covariance matrix for the two traits. The first two indices determine the coefficient in the G matrix and i or j stand for either allele at the rQTL. We simplify equation (2.1) by replacing the constant terms with S j ¼ a 2 1 g 11j þ 2a 1 a 2 g 12j þ a 2 2 g 22j and S i ¼ a 2 1 g 11i þ 2a 1 a 2 g 12i þ a 2 2 g 22i and DS ¼ S i 2 S j resulting in the pleasantly simple recursion equation When D w 0 ði; jÞ = 0, selection occurs among the rQTL alleles. Because the genotypes at the rQTL are associated with different G matrices, selection on the rQTL changes the covariance structure of the traits.
Several conclusions can be drawn from equation (2.2). First, we can see that the fitness difference between the rQTL alleles will increase or decrease depending on the sign of DS. This term gives the selection response of the two genotype classes to directional selection, depending on the rQTL allele-specific genetic covariance matrices and the direction and strength of selection. This change in mean fitness diminishes owing to recombination depending on the degree of linkage L. Note that even at the lowest level of linkage, the fitness difference does not completely disappear, because the LD does not relax within one generation.
Next, the bivariate genotypic values and the difference in mean fitness can be analytically derived in any generation. Given that the two alleles at rQTL have equal fitness in the first generation, iterating equation Eventually, the mean fitness difference reaches a stable fixed point if linkage is not perfect, L , 1. The stationary fitness difference is After the fitness difference has equilibrated, selection at the rQT locus proceeds with constant selection coefficients. This simplifies the analysis of selection dynamics.
We also explored the situation where one trait is under directional selection while the second trait is under stabilizing selection, a fitness function known as the 'corridor model' [28,29]. We applied an analogous logic as above but a different combination of selective pressures to select for dissociation of traits. This model is described by the following fitness function: where a is a coefficient of directional selection acting on trait z 1 , and b is a coefficient of stabilizing selection acting on trait z 2 .
The average fitness for the subpopulation with allele i at the rQTL is then and differentiation yields the following selection gradients: 22 :

SIMULATION RESULTS
The model was numerically simulated to illustrate the effects of rQTL selection on the evolution of trait correlations. Figure 2a -c shows the change in the trait correlation with different degrees of linkage L (different curves in each plot) and different directions and kinds of selection. For transparency, we show only symmetrical cases, where the two rQTL genotype classes differ in covariance while the genetic variances are held equal, and where the selection on the two traits is either equal or of equal magnitude but different signs, or where only one trait is under directional selection and the other under stabilizing selection. It is easy to see from equation (2.1), however, that differences in variances combine with covariance differences in determining rQTL selection. In other words, the genotype with higher variance along the direction of selection, be it univariate or multivariate (i.e. the selection along the linear combination of trait values and correspondingly the covariance of traits), has a selective advantage.
In figure 2a, we show a scenario where the two traits are co-selected, i.e. both traits are under directional selection of equal strength. The initially most frequent rQTL allele has zero correlation, and the alternative allele causes a correlation of 0.9. As can be seen, the correlation among the traits increases, because of the selection of the rQTL allele that causes higher genetic covariance. The reason is intuitively easy to see. The rQTL genotype class with higher genetic covariance has more genetic variance in the direction selection. This means that the mean fitness of this genotype class increases faster than that of the rQTL genotype class with lower genetic covariance. Thus, the rQTL causing higher genetic correlation is selected, even in the face of moderate to high recombination. High levels of recombination decrease the rate of selection on the rQTL and thus also diminish the rate of increase in trait correlation.
With antagonistic selection, that is selection increasing the mean of one character while decreasing the mean of the other character, the rQTL allele gets selected that reduces the correlation or can even lead to negative correlation (figure 2b). When directional selection acts only on one character, while the other character is neutral, however, rQTL selection is not possible in this model (not shown). The reason is simply that if one character is neutral, then neither the mean value of that character nor the presence or absence of genetic correlation affects the fitness of the rQTL allele. The assumption that one of the two characters is totally neutral does not seem to be plausible, however. It is more likely that if one trait is under directional selection then the second trait is under stabilizing selection, a fitness function known as the 'corridor model' [28,29]. Figure 2c shows the results of simulating rQTL selection in the corridor model. Directional selection on a single character, while the other is under stabilizing selection, favours zero correlation among the characters. In this case, a genetic correlation is causing a correlated selection response in the trait under stabilizing selection, which causes it to deviate from its optimal value. Hence, the correlated selection response leads to lower fitness compared with genotypes that have no genetic correlation, and thus genotypes with lower correlation are selected. rQTL selection is also effective even under very weak directional selection. In figure 2a,b, the coefficients a 1 and a 2 are +0.05, while in figure 2c, a ¼ 0.002.
In the scenarios considered so far, directional selection continued over many generations in the same direction.
While this might be plausible in the case of long-term climatic trends or artificial selection, this assumption seems to be extreme. We therefore considered scenarios where directional selection is fluctuating, meaning that a pair of co-selected characters experience simultaneous selection for increasing or decreasing the trait means, as in situations where fluctuating temperatures may lead to alternating selection for greater or lesser body size. From the symmetry properties of equations (2.1) and (2.2), it can be seen that the same rQTL allele would be selected regardless of whether the traits are co-selected for larger or smaller means. It follows that fluctuating co-selection should also lead to increased absolute correlations among co-selected traits even when the direction of selection is fluctuating. We confirmed this intuition with simulations shown in figure 3.
The simulations presented in figure 3 are analogous to the previous ones in figure 2a,b, except that the signs of both selection coefficients change simultaneously every 100 generations. We have implemented two scenarios: in figure 3a, the population is allowed to reach linkage equilibrium between episodes of directional selection. This scenario simulates decay of positive linkage equilibrium during a short period of stabilizing selection between epochs of directional selection. In the second scenario in figure 3b, selection is reversed immediately following the previous epoch of directional selection. starting with a population that contains 1% of an allele with high correlation between traits (r i ¼ 0.9), and 99% of an allele causing zero correlation. Selection leads to an increase in genetic correlation between the two traits until the allele causing high correlation is fixed. (b) A scenario in which selection is acting on the traits in opposite directions, i.e. antagonistic selection, causing increase in one and decrease in another trait (a 1 ¼ 0.05, a 2 ¼ 20.05). At the beginning, the most frequent allele is the one that causes zero genetic correlation, and the alternative allele causes a correlation of 20.9. Antagonistic selection favours the allele causing negative correlation. (c) Evolution of genetic correlation in a corridor model, i.e. directional selection on one trait, while the other trait is under stabilizing selection. The initially most frequent allele causes a positive correlation of 0.9, and the alternative allele leads to zero correlation. The genetic correlation decreases as the allele causing the lower correlation between traits is selected (r 1 ¼ 0, r 2 ¼ 0.9, directional selection coefficient ¼ 0.001, stabilizing selection coefficient ¼ 0.0025). Note that rQTL selection in the corridor model is much more effective than with directional selection alone. The strength of directional selection in this simulation is only 0.001, while in plots (b) and (c) it is 0.05 in order to reach the same level of correlation change in the similar amount of time.
The difference between these scenarios is important. LD develops between the rQTL alleles and the rest of the genes owing to selection. The rQTL allele causing faster evolution is associated with genes that cause trait means more advanced in the direction of selection. When the direction of selection changes, the advantageous rQTL allele is associated with the 'wrong' genes and is temporarily selected against until the LD is reversed. Hence, if selection changes direction without allowing linkage equilibrium to decay between epochs, the evolution of trait correlation takes longer. These episodes of counter-selection at the beginning of each epoch, owing to unresolved 'wrong' LD, are reflected in a characteristic saw-tooth pattern in figure 3b. If however the linkage equilibrium is established between opposing epochs of directional selection, the correlation between the traits evolves at about the same rate. This can be seen from the comparison of figures 2a and 3a, both showing positive co-selection of both traits. Figure 3a also shows slight 'wrinkles' in the curve. These are caused by the fact that at the beginning of each epoch, the LD necessary for selecting the rQTL allele has to be re-established, leading to a short period of slow selection. Figure 3c shows the same comparison but for the case of  (a) In this scenario, the linkage equilibrium between the rQTL and the rest of the genome is restored between each epoch of directional selection. This is equivalent to the assumption that between each epoch there is a period of stabilizing selection during which linkage equilibrium is restored (not simulated, but linkage equilibrium is enforced by averaging genotypic values). The rate of evolution of genetic correlations is similar to the scenario of sustained directional selection, only slightly delayed at the beginning of each epoch. The reasons for these delays are two. First, at the beginning of each epoch, LD has to be rebuilt that fuels the selection of the rQTL allele (see the plot of the rQTL allele frequency in the electronic supplementary material, figure S1a, which does not show a drop in frequency but only a slowdown in increase). Second, the relaxation of LD between epochs leads to a decrease of the population-level genetic correlation, which explains the instantaneous drop of correlation visible in the L ¼ 0.9 curve at generation 100 and somewhat less so in the other curves and epochs. (b) Same scenario as in (a) except that the LD is not relaxed between epochs of directional selection. This scenario leads to a pronounced sawtooth pattern at the beginning of each new epoch. The reason is that initially, at the switch from one epoch to the next, the favoured rQTL allele is associated with genotypic values of the traits that cause lower fitness in the new selection regime. This 'wrong' LD is dissipated over a few generations and then selection for the favoured rQTL allele resumes. Overall, these effects do not significantly delay the evolution of the favoured form of genetic correlation. (c) Fluctuating selection in the corridor model, with equilibration between the epochs. The direction of the directional selection changes, while the stabilizing selection is constant. There is no effect of linkage relaxation in this model (not shown).
fluctuating directional selection on one trait and constant stabilizing selection on the second trait (i.e. corridor model). Other than the delays caused by the need to invert LD after each epoch of directional selection, this type of fluctuating selection is as effective as continued directional selection in shaping the correlation pattern of traits. Hence, rQTL selection under fluctuating environmental conditions requires some degree of recombination, because the LD needs to be inverted at the beginning of each epoch of directional selection. With full-linkage rQTL, selection is impossible under fluctuating environments (see the electronic supplementary material, figure S2). Note that the analytical solution in equation (2.2) and its consequences refer only to the persistent directional selection on both traits, but not to the corridor model or fluctuating selection.

DISCUSSION
Modelling the evolution of variational properties requires modelling multi-loci dynamics with gene interaction and multiple phenotypic characters under directional selection including the effect of LD. In most mathematical models most, if not all, of these complications are ignored in order to make the models tractable. The recent pioneering work by Jones et al. [30 -32] focuses on stabilizing selection. In this project, we chose to make a number of alternative simplifications that preserve the presence of gene interaction and LD but use other abstractions. Specifically, we model the evolution of genes affecting two characters in the form of a phenotypic Lande equation and represent LD by an association between rQTL alleles and the genotypic values of the quantitative traits. To our knowledge, this model is the first to accommodate gene interaction, LD and directional selection and is still mathematically tractable. Certainly, the model is an abstraction that involves several simplifications. In particular, the G-matrix associated with a particular rQTL is kept constant and does not account for G-matrix dynamics caused by allele frequency changes. Furthermore, we assume that the rQTL alleles do not develop an additive effect. Epistatic rQTL can change their additive effect owing to background change [33]. We currently lack a general theory predicting the effect of epistasis on the multitrait phenotype. Nevertheless, the parametrization of an analogous epistatic effect on univariate characters has been developed in the framework of the multi-linear model [33,34]. Within this framework, the overall effect of epistasis on the locus' additive genetic effect on a single trait can be assessed. The parameter of the multilinear model that measures the change in additive effect of the epistatic locus is the directionality of epistasis. Directionality describes how the additive effect of this locus on a phenotypic trait is modified on average owing to genetic interactions with the background loci. Negative values of directionality indicate diminished additive effect at the locus owing to epistasis, while the positive directionality indicates an increase of the additive effect owing to epistasis. A value close to zero implies that the effects of epistasis on the additive genetic effect of a locus cancel out in total across background loci. In this terminology, the assumption of the constancy of genetic effect across generations implies that the directionality of epistasis of the rQTL is zero for both phenotypic traits. A recent empirical study documented low net effect of epistasis on additive variance in mice intercross [35]. Additional studies of directionality are required to show how frequently this assumption holds.
Hence, at this stage we do not include all the complexities of the G-matrix dynamics. It is known that even with an additive model, the evolutionary dynamics of the G matrix is very difficult to model mathematically. The difficulties are further increased greatly if epistasis is taken into account. The benefits of doing so, however, could be quite limited. It is well understood that under sustained directional selection, the genetic variance approaches quasi-equilibrium between novel mutations and fixation of alleles [36]. This equilibrium is equal to the neutral mutation selection equilibrium, and thus is proportional to the mutational variance of the character. Furthermore, artificial directional selection has relatively small effects on genetic variance [37]. Another potential expansion of a current model is by including multiple rQTL, as in all documented cases, intertrait relationship is affected by more than a single rQTL. While our model manifestly is an abstraction, we suggest it is one that can be highly productive until a stronger mathematical approach is found that can deal with the complexities of multiple characters, polygenic inheritance with gene interaction and LD.
In conclusion, the model presented here shows that, given the presence of 'neutral' rQTL, directional selection can cause an adaptive change in the genetic covariance structure. The result is that the direction of highest genetic variance aligns with the direction of selection ( [38]; see figure 1c). This mechanism requires the build-up of LD between the rQTL and the rest of the genome. Our simulations show that even strong recombination can sustain rQTL selection because recombination does not instantly destroy LD. Selection on rQTL can decrease variational constraints caused by trait covariance. If directional selection acts on only one character, and the other character is under stabilizing selection, lower covariance, i.e. the dissociation of traits, is favoured. Such increased variational autonomy of traits facilitates independent adaptation and differentiation of traits and can contribute to an increase in organismal complexity, i.e. the individualization of body parts. As mentioned in §1, this is what happened in the stem lineage of anthropoid primates, probably facilitating the limb shape diversity of higher primates [16]. The model furthermore also accommodates fluctuating directional selection, which seems more plausible than long-term directional selection. In fluctuating selection, recombination plays an important role in facilitating the evolution of trait correlations.
In this model natural selection acts directly on evolvability by selecting genotypes that increase the rate of response to directional selection. The idea that selection could directly act on evolvability is controversial [39], but in recent years computational models have shown that selection for evolvability is probably a generic property of all but the most simplistic genetic models ([40-45]; see also [46] for an empirical example). The present model, together with the demonstration of the existence of quasi-'neutral' rQTL variation [17,23,24], shows that selection may be able to modify the genetic architecture to increase evolvability and, specifically, increase the response rate to directional selection. In essence, we show that through differential epistatic interactions, selection can produce the 'right kind' of heritable phenotypic variation to accomplish adaptive change from the pre-existing genetic, but phenotypically cryptic, variation. Finally, what are the possible physiological mechanisms causing this kind of variation? Complete account of the possible mechanisms of variation in pleiotropy is beyond the scope of the current study, yet we include an example here. He & Zhang [47] found in yeast that most of pleiotropy is owing to a single gene product involved in different spatial or temporal contexts, each context resulting in a different phenotype. They furthermore found that the degree of pleiotropy correlates positively with the number of protein -protein interactions of the product [47]. It is plausible then that the variation in the relationship between the end-phenotypes of these pathways arises because genes interacting with the focal product in one pathway do not overlap with the genes involved in other pathways, leading to somewhat independent evolvability of the multiple functions at the pleiotropic locus. This was exemplified in a study of yeast by Kim et al. [22], in which the authors analysed the molecular basis of a pleiotropic mutation at CYS4 (cystathionine beta-synthase) that causes rust coloration of colonies growing on copper sulphate and differences in drug sensitivity [48,49]. The single nucleotide mutation was shown to lower the levels of cysteine, an essential amino acid involved both in a sulphur assimilation pathway and in a detoxication pathway. Low levels of cysteine cause increased assimilation of sulphur and therefore rust coloration, and sensitivity to multiple drugs. Different gene products interact with the focal product, cysteine, in different pathways. Some of the interacting genes involved are specific to individual pathways, hence mutations in these genes (i.e. different backgrounds) affect the two phenotypes separately, causing variation in the strength of the relationship between them, in spite of a common gene product involved in both. The rQTL found in this case are the epistatically interacting loci affecting only one of the functions (pathways) of the pleiotropic locus, and causing different levels of correlation in different genotypes. Selection on a particular intertrait relationship according to our model would therefore fix the genotype at the epistatic locus. More focus on the molecular mechanisms involved will reveal further possibilities promoting variation in pleiotropy.
We thank Bruce Walsh for the comments on the manuscript. M.P. acknowledges support from Erwin Schrö dinger Postdoctoral Fellowship of Austrian Science Foundation (FWF). J.M.C. acknowledges funding from NSF (no. BCS-0725068) and NIH (nos AR053224, DK055736 and GM065509). Research in the GPW laboratory is funded by the John Templeton Foundation. The opinions expressed in this article are not those of the John Templeton Foundation.

APPENDIX A
We present here two conceptually distinct ways to model the effect of recombination on the mean trait value. The model we used in the study is a phenomenological model. We also explain a classical population-genetic alternative below.
In the phenomenological model, we do not consider an explicit mechanism of how the genotypes of the next generation are created, but assume that the mean phenotype of the offspring genotypes of a mating of genotypes with different average phenotype is a simple mixture of the average parental phenotypes. For instance, let us assume that the average phenotype of the genotypes with R i is z i and that of R j is z j . Now let us calculate the average phenotype of the R i genotype in the next generation owing to recombination.
In the offspring generation, R i genotypes can arise from two types of matings: (i) matings of type R i Â R i , and (ii) matings of R i Â R j . Under random mating, the first kind of pairing, R i Â R i occurs with frequency p 2 . In this case, all offspring will be of R i type and will have the average phenotype z i , where p is the frequency of R i genotypes and q ¼ 1 2 p is the frequency of the R j genotypes. In the second, matings of R i Â R j , case pairings will occur with a frequency 2pq. Of the offspring with this pairing, one-half will be of R i type and these will have the average phenotype L z i þ ð1 À LÞ z j ; where L is a phenomenological parameter of linkage. If L ¼ 1, there would be no recombination and the R i genotypes would retain their average parental phenotype. Assuming L ¼ 0.5 would mean that the R i genotypes, coming from a R i Â R j mating, would have the same average phenotype as the R j genotypes coming from the same mating. This would mean that, within one generation, the LD between R i and R j genotypes would be completely relaxed. We know that LD does not relax within one generation and thus restrict our phenomenological linkage parameter to the range [0.75, 1]. If L ¼ 0.75, the resulting average phenotype will be half way between the average phenotypes of the R i parents and the linkage equilibrium phenotype ð z i þ z j Þ=2. Hence, the assumption is that the highest level of recombination will decrease the LD by a factor of 0.5. Taking the average phenotype resulting from the R i Â R i and the R i Â R j matings together and normalizing it by the frequency of R i genotypes, p, we obtain z 0 i ¼ ð p þ qLÞ z i þ qð1 À LÞ z j ; and analogously for the average phenotype of the R j genotypes we obtain z 0 j ¼ ð pL þ qÞ z j þ pð1 À LÞ z i : An alternative method of deriving an analogous result is to consider the mean effect loci as a block, with no recombination within, only recombination relative to the rQTL. This way the whole block becomes associated with an allele at the rQTL, and the model becomes a haploid two-locus model. This can be modelled in the conventional way, using the rate of recombination c, with parameter space 0 -0.5. Consider a haploid organism with two alternative gametes, i and j. We treat mean-effect-loci as a composite locus M, and thus we consider two loci on each gamete: R (¼rQTL) and M. Given the rate c of recombination between R and M, the two-locus genotypes in the offspring generation have the following frequencies and phenotypes: parents frequency (mating) gametes frequency (gametes) phenotype (1 2 c)/2 z j Selection for evolvability M. Pavlicev et al. 1911