THE G MATRIX UNDER FLUCTUATING CORRELATIONAL MUTATION AND SELECTION (original) (raw)

Abstract

Theoretical quantitative genetics provides a framework for reconstructing past selection and predicting future patterns of phenotypic differentiation. However, the usefulness of the equations of quantitative genetics for evolutionary inference relies on the evolutionary stability of the additive genetic variance–covariance matrix (G matrix). A fruitful new approach for exploring the evolutionary dynamics of G involves the use of individual-based computer simulations. Previous studies have focused on the evolution of the eigenstructure of G. An alternative approach employed in this paper uses the multivariate response-to-selection equation to evaluate the stability of G. In this approach, I measure similarity by the correlation between response-to-selection vectors due to random selection gradients. I analyze the dynamics of G under several conditions of correlational mutation and selection. As found in a previous study, the eigenstructure of G is stabilized by correlational mutation and selection. However, over broad conditions, instability of G did not result in a decreased consistency of the response to selection. I also analyze the stability of G when the correlation coefficients of correlational mutation and selection and the effective population size change through time. To my knowledge, no prior study has used computer simulations to investigate the stability of G when correlational mutation and selection fluctuate. Under these conditions, the eigenstructure of G is unstable under some simulation conditions. Different results are obtained if G matrix stability is assessed by eigenanalysis or by the response to random selection gradients. In this case, the response to selection is most consistent when certain aspects of the eigenstructure of G are least stable and vice versa.

The additive genetic variance–covariance matrix, or G matrix, has become an important tool in predicting the course of phenotypic evolution under natural selection and genetic drift (Lande and Arnold 1983; Arnold 1992; Arnold et al. 2001). G is a symmetrical matrix in which the diagonal elements are additive genetic variances for traits and the off-diagonal elements are additive genetic covariances (Lynch and Walsh 1998). Quantitative genetic theory comprises, at its base, a set of equations that use this matrix to describe the course and pattern of phenotypic evolution (Lande 1979; Lande and Arnold 1983; Arnold et al. 2001). Central among these equations is the multivariate version of the classic breeder's equation (Lande 1979), which provides the single-generation response to selection,

1

in which is a vector of single-generation changes in population means for traits and β is a vector of the partial regression coefficients of fitness on phenotype (i.e., the selection gradient; Lande 1979; Lande and Arnold 1983). As the response to selection depends both on the genetic variation for single traits and on the genetic correlations between them, the G matrix plays a central role in determining the course of adaptive evolution.

Lande (1979) also proposed that the net selection gradient (β_net_) responsible for phenotypic change over many generations could be obtained by summing equation (1) from time 0 to time t, such that

2

in which is a vector of trait means at time t. In addition, if the populations or species had phenotypically diverged under drift, then the covariance matrix for species mean trait values, the D matrix, has an expectation proportional to the time average of G,

3

in which Ne is the effective population size. An implicit assumption of equation (2) and a practical necessity of equation (3) is that G is stable over the evolutionary time scale of interest. As such, successful retrospective selection analysis (eq. 2) and prospective prediction of diversification (eq. 3) depend on the evolutionary dynamics of G (Turelli 1988; Arnold 1992; Arnold et al. 2001).

Theoretical considerations provide no concrete predictions about the long-term dynamics of G and consequently the stability of G is fundamentally an empirical question (Barton and Turelli 1987; Turelli 1988; Shaw et al. 1995; Roff 2000). Although there is accumulating evidence that G may be conserved, particularly at low taxonomic levels (Roff 1997; Roff and Mousseau 1999; Steppan et al. 2002; Bégin and Roff 2004; Revell et al. 2007), theoretical investigations are needed to explore the evolutionary dynamics of G. Because the dynamics of G are not particularly amenable to analytic exploration (see Jones et al. 2003), several authors have made successful use of individual-based, stochastic computer simulations to examine the shape and stability of G under various conditions of selection, mutation, and genetic drift (Bürger et al. 1989; Bürger and Lande 1994; Reeve 2000; Jones et al. 2003, 2004).

In the most thorough such study to date, Jones et al. (2003) showed that the stability of some aspects of the eigenstructure of G is increased by correlational natural selection and correlated effects of pleiotropic mutation. In particular, the stability of the orientation of the primary eigenvector of G is increased by correlational mutation or selection or both. They also note, however, that under several mutation–selection parameter combinations (henceforward “mutation–selection scenarios”), instability in the orientation of the primary eigenvector of G is associated with low matrix eccentricity and thus is probably inconsequential to evolution by natural selection (Jones et al. 2003; p. 1758). This is because evolution is not constrained by a matrix that is not eccentric and thus matrix orientation is irrelevant to evolution. They also show that stability in the magnitude of the elements of G is lower for small effective population size (Jones et al. 2003).

In this paper, I use eigenanalysis and the response to random selection vectors to measure stability and to explore the conditions under which G can be considered to be stable. The latter approach, which is adapted from the random skewers method of Cheverud et al. (1983; Cheverud 1996), measures the evolutionary dynamics of G in terms of the correlation between the response to random selection vectors. This measure of the dynamic of G differs from the eigenanalysis in which the stability of G is measured directly. As such, I will henceforward refer to the random skewers approach as measuring the evolutionary consistency of G. This is because random skewers measures the consistency of the expected response to selection given the same selection gradient imposed on two different G matrices. As this method uses the quantitative genetic equation (1), it may be useful in diagnosing differences between G matrices that are likely to influence the response to natural selection.

I also simulate previously unexplored conditions under which G may be more severely destabilized. In particular, I explore the degree to which the stability and evolutionary consistency of G are compromised by fluctuating coefficients of correlational selection and pleiotropic mutation, and by fluctuating effective population size. To my knowledge, no prior simulation study has investigated the consequences of fluctuating correlational selection and mutation on the evolution of the G matrix.

In summary, therefore, this study uses different methods to: (1) evaluate the evolutionary dynamics of G under various conditions when correlational selection, correlational pleiotropic mutation, and effective population size are constant over time; and (2) evaluate the evolutionary dynamics of G when correlational selection and mutation, and effective population size fluctuate over time.

Methods

SIMULATING THE EVOLUTION OF G

The simulation model used in this study is a Monte Carlo simulation with all individuals modeled and is similar to that employed by Jones et al. (2003, 2004), and in a univariate context by Bürger et al. (1989) and Bürger and Lande (1994). In this model, two traits were determined by _n_= 50 unlinked pleiotropic loci. A mutation at a locus produced a new allele with new effects on both traits. The new effects were drawn from a bivariate normal distribution with mean [0, 0], variances α12=α22= 0.05, and correlation _r_μ. The new state at the locus was then determined by adding new effects to the old. This corresponds to a continuum of alleles mutation model (Crow and Kimura 1964; Kimura 1965). I determined individual phenotypes by summing across all loci and adding to both traits random normal environmental deviates with mean [0, 0], variances _E_11=_E_22= 1.0, and correlation _re_= 0.

I simulated a diploid, hermaphroditic, sexually reproducing population with discrete generations and population size N, which, for most simulations, was held constant at _N_= 1000. I used a slightly different life cycle than that employed by Jones et al. (2003). The life cycle was as follows: (1) I calculated phenotype and expected fitness for each adult; (2) I randomly mated the adults with probability equivalent to their expected fitness—each mating produced one offspring and having participated in one mating did not diminish the probability of mating again. This was repeated until the population of progeny was equal to N. (3) I produced offspring genotypes by independently assorting one random gene copy from each parent at each locus into each offspring, which corresponds to a model of no physical linkage among loci, and simulated mutations as described above with a probability per allele generation of μ= 0.0002.

The main difference between the life cycle simulated in this study and that of Jones et al. (2003) is that Jones et al. impose mortality selection prior to reproduction, whereas I mated adults with a probability determined by their fitness. For a subset of the simulations presented in this study I also conducted the simulation using prereproductive mortality selection (not shown). I found no difference between the results in those cases. This should not be surprising given that the probability of any mating pair is the product of the fitness parameter of each mate (calculated as described below). This makes the selection model used in this study equivalent to multiplicative fecundity selection, which has been shown elsewhere to be exactly analogous to the two-sex viability selection model used in Jones et al. (2003; Bodmer 1965; Hartl and Clark 1989).

In addition, in Jones et al. (2003), the effective population size Ne was larger than the population size N, whereas in my simulations Ne and N were equivalent. This circumvents the necessity of complicated calculations to evaluate Ne given N.

I assigned expected fitness, W(z), to each individual according to the bivariate Gaussian fitness function,

4

where z is a column vector of individual trait values (the superscript T indicating that a transpose was taken), and ω is a symmetrical two by two matrix describing the curvature and orientation of the selection surface (Lande 1979). Following Jones et al. (2003, 2004), I set ω11=ω22= 49, which corresponds to a scenario of weak stabilizing selection. The strength of correlational selection, ω12, was set equal to _r_ω(ω11ω22)1/2, where _r_ω is the correlation coefficient of the fitness surface.

Each simulation run involved the initiation of the population at genetic uniformity, followed by 10,000 generations of stabilizing selection, mutation, and drift, during which the starting population attained equilibrium conditions. This “burn-in” was then followed by 2000 experimental generations during which values for the elements of the additive genetic variance–covariance matrix, three aspects of its eigenstructure (the size, Σ; a measure of the eccentricity, ɛ; and the orientation, ϕ), as well as the phenotypic means for the traits were recorded in each generation. Size, eccentricity, and orientation of G were calculated following the methods described in Jones et al. (2003). Their specific calculation is also described below.

I conducted two sets of experimental manipulations using this design. The first set involved repeating the simulation conditions of Jones et al. (2003), in which mutation, selection, and drift were simulated for constant correlational selection and mutation and constant population size. I repeated their simulations to analyze the results using a response-to-selection vector correlation method to assess evolutionary consistency (described below), and to consider the results obtained using this method in the context of those calculated when stability is evaluated directly by decomposing the eigenstructure of the matrix. For this analysis, I used six different mutation–selection scenarios (following Jones et al. 2003). Although in this text I assign each scenario a letter indicated in parentheses (e.g., scenario [_a_]), for clarity when the scenarios are subsequently referenced I will provide both the letter index and description of the scenario wherever possible (so that the reader need not memorize the specific parameter combinations of each scenario for which I conducted the simulations). The six mutation–selection scenarios were as follows: (a) no correlational selection or mutation (_r_ω= 0, _r_μ= 0); (b) correlational selection with uncorrelated mutation (_r_ω= 0.75, _r_μ= 0); (c) correlational mutation with uncorrelated selection (_r_ω= 0, _r_μ= 0.5); (d) correlational selection and mutation (_r_ω= 0.75, _r_μ= 0.5); (e) correlational selection and mutation with high values of the correlation coefficients (_r_ω= 0.9, _r_μ= 0.9; Jones et al. 2003); and (f) correlational mutation and selection in which the correlation coefficients have opposite signs (_r_ω= 0.75, _r_μ=−0.5). This final scenario corresponds to one of antagonistic pleiotropic mutation, in which a mutation with a beneficial or only mildly deleterious effect on one trait has a severely deleterious effect on the other trait. These simulation conditions explore a range of parameter values for _r_ω and _r_μ, both of which are poorly known empirically. I simulated each set of conditions with 20 runs.

The second set of manipulations involved simulating conditions under which G might be severely destabilized. In these simulations I manipulated the constancy of _r_ω, _r_μ, and N. I initiated _r_ω at 0 and allowed _r_ω to change over time as a random walk, where _r_ω(_t_+ 1) was equal to _r_ω(t) plus a deviate drawn from a random normal distribution with mean 0 and variances (g) σ2(r) = 0.0001, (h) σ2(r) = 0.001, (i) σ2(r) = 0.01, and (j) σ2(r) = 0.1. These range from very slowly (g) to very rapidly (j) changing _r_ω. In a separate set of analyses, I initiated _r_μ at 0 and allowed _r_μ to change over time, in which _r_μ(_t_+ 1) was equal to _r_μ(t) plus a deviate drawn from a random normal distribution with mean 0 and variances σ2(r) as described above for _r_ω (mutation-selection scenarios [_k_] through [_n_]). Necessarily, −1 ≤_r_≤ 1, so if a correlation hit the boundary it was reflected by an amount equivalent to that by which it would have otherwise exceeded the bound. Example time series for r are provided in Appendix A. Finally, I manipulated population size N. In these simulations, run for conditions of no mutational and selectional correlation (as in mutation-selection scenario [_a_], above), N(0) was initiated at _N_= 1000, after which N(_t_+ 1) was equal to N(t) plus a deviate drawn from a random normal distribution with mean 0 and variances (o) σ2(N) = 10, (p) σ2(N) = 100, (q) σ2(N) = 1000, or (r) σ2(N) = 10,000. The boundary conditions of 20 ≤_N_≤ 2000 were set to avoid population extinction as well as the complete loss of genetic variation, and to increase computational expediency. Reflection off the bound was performed as it was for r in the previous simulations. Example time series for N are provided in Appendix B. Although _r_ω, _r_μ, or N were varied, other parameters were held constant. In particular, _r_ω (except in simulations [_g_] through [_j_]) and _r_μ (except in simulations [_k_] through [_n_]) were held fixed at 0, and (except in simulations [_o_] through [_r_]) N was held fixed at 1000, as in simulations (a) through (f). Each set of simulation conditions was explored with 20 replicate runs.

MEASURING G-MATRIX STABILITY

I measured the evolutionary dynamics of G using two different approaches. The first was the approach of eigenanalysis (e.g., Jones et al. 2003, 2004). I analyzed stability by evaluating the constancy of various aspects of the eigenstructure of G: the matrix size (Σ), calculated as the sum of the eigenvalues of G; an inverse measure of the eccentricity (ɛ; henceforward for brevity referred to as the eccentricity), calculated as the ratio of the smaller to the larger eigenvalue of the matrix; and the orientation (ϕ) of G, calculated as the angle in degrees between the leading eigenvector and the first trait axis. To evaluate the stability of G, I calculated the mean cross-generational change in the elements of G, as well as in the total size (ΔΣ), the eccentricity (Δɛ), and the orientation (Δϕ; Jones et al. 2003). These were used as measures of the stability of G and of important aspects of the eigenstructure of G. This approach allowed direct measurement of the temporal stability of G and its eigenstructure.

The second approach that I used to measure the evolutionary dynamics of G uses an adapted version of the random skewers method (Cheverud et al. 1983; Cheverud 1996; Marroig and Cheverud 2001). The random skewers method measures the similarity between the two variance–covariance matrices by first multiplying each by a single randomly generated selection gradient. For each matrix, the product is a vector that is equivalent to the expected response to selection given the covariance matrix and the random selection gradient. The vector correlation between the response-to-selection vectors from each of the two matrices, averaged across many random selection gradients, is a measure of the similarity of the matrices (Cheverud 1996).

Specifically, for matrices G1 and G2, and m randomly generated selection vectors, β1 through βm, the random skewers correlation (rs) is calculated as

5

where corr denotes the calculation of a vector correlation

6

Elements of the random vectors were drawn from a random normal distribution with mean 0.0 and variance 1.0. Generating random vectors in this way results in a set of vectors whose orientations are random and uniformly distributed (i.e., they are random vectors from a unit circle). As before, superscriptT indicates that a transpose is used. This calculation should be familiar as a standardized vector dot product.

The random skewers correlation is a useful metric of similarity because it uses the multivariate response-to-selection equation (eq. 1), and thus provides a measure of G matrix resemblance in terms of the expected similarity of the multivariate response to natural selection. As the random skewers matrix correlation is a measure of the consistency of the expected evolutionary responses of two matrices to natural selection, and to distinguish it from the direct measures of matrix stability, I henceforward refer to random skewers as providing a measure of the evolutionary consistency of G.

My modification of this approach, which I call “serial random skewers,” measures the temporal evolutionary consistency of G using the mathematics of the random skewers method. It involves calculating the average arccosine of the random skewers vector correlation across all possible time steps in each 2000-generation run, and then averaging this result across runs. The arccosine provides the angle between the response-to-selection vectors (henceforward θ is used to indicate the angle, in degrees, between response-to-selection vectors, with subscripts or hats to differentiate specific contrasts and averages; Fig. 1). θ has the convenient interpretation of representing the angle, in degrees, over which the response to selection is “bent” by substituting matrix 2 for matrix 1 or vice versa. Its average across all possible time steps of length t is (e.g., for 2000 generations between _t_= 0 and t_= 2000 there are 2001 −Δ_t time steps of length t). I averaged for Δ_t_= 1 to Δ_t_= 2000 across 20 replicate runs of each simulation scenario. _m_= 100 random selection vectors were used at each time step (see eq. 5). was thus evaluated as follows:

7

Obviously is calculated from both independent and nonindependent observations, where nonindependence results both from time-step overlap and from temporal autocorrelation of G. As such, was used as a quantitative measure of evolutionary consistency but not to assess whether significant differences in consistency were observed among different mutation–selection scenarios. For that purpose a different measure of matrix consistency was used (see below).

1

An illustration of the random skewers procedure. G matrices are represented by ellipses GA and GB, for which orientation is determined by the primary eigenvector of G and length and width are proportional to the square roots of the eigenvalues of G. Each Gi is multiplied by the random selection vector, β. A measure of matrix similarity is provided by the angle, θ, between the response-to-selection vectors, RA and RB. Note how RA and RB are “bent” toward the primary eigenvectors of GA and GB, respectively. Matrix similarity is evaluated as the mean value of the angle, θ, from many random selection vectors.

Interpretation of is straightforward. If is relatively small compared to a null model, it indicates that the mean deviation from collinearity of the response to selection because of the differences among G matrices sampled Δ_t_ generations apart in time is small, and thus the evolutionary consistency of G is high. In other words, if is small, then swapping G matrices over the time interval Δ_t_ would be of little consequence to the direction of the response to selection. By contrast, relatively large implies that evolutionarily important differences between the matrices exist (Fig. 1).

I also evaluated the significance of the differences in evolutionary consistency among mutation–selection scenarios. Admittedly, hypothesis testing in this context has some artificiality associated with it—because we know a priori that our null hypothesis is false (the data were generated under different simulation conditions). Nonetheless, evaluating the significance of the difference between simulation conditions can provide a reasonable indication of whether discernibly different evolutionary consequences of the different simulation models are implied.

Because the means at each Δ_t_ (denoted as ) in the preceding analysis were calculated from both independent (nonoverlapping and between-run) and nonindependent (overlapping, within-run) observations (see eq. 7, above), standard errors cannot be estimated in the typical way. Thus, to facilitate hypothesis tests of G matrix consistency, I also calculated θ_t_,0, the angle, in degrees, between the response vectors to random selection gradients at times t and t_= 0, for each value of t > 0 and averaged across the 20 replicates of each simulation scenario. Variance in θ_t,0 was calculated in the standard way across runs at each t > 0, which is possible because at each time step only independent (between-run) observations are used. I evaluate the significance of differences in G matrix consistency among mutation–selection scenarios by calculating the _t_-statistic at all time intervals as

8

where time is indicated in the subscript and the genetic scenarios are shown in parentheses. tt(i,j) is expected to be distributed as a _t_-statistic with _df_= 38. Genetic scenario (a: _r_ω=_r_μ=σ2(_r_ω) =σ2(_r_μ) =σ2(N) = 0, and _N_= 1000) was used as the null scenario i in all analyses such that a positive value for the _t_-statistic indicated that the alternative scenario promoted the evolutionary consistency of G relative to its stability with constant correlational mutation and selection equal to zero, and constant population size, whereas a negative result indicated that the alternative mutation–selection scenario decreased the consistency of G.

Although each value of the _t_-statistic was calculated from a set of independent (between-run) observations, sequential values of t are highly serially autocorrelated, making the calculation of specific _P_-values for the test very difficult. Consequently, I present a temporal profile of the _t_-statistic for each comparison. Only profiles that consistently exceed the significance threshold of the two-tailed critical value from the _t_-distribution, _t_0.05(df = 38) = 2.02, are considered to indicate a significant difference in consistency with respect to null scenario (a).

I also measured similarity in the relative magnitude of the response to selection. Magnitudinal similarity, M, was measured as the ratio of the lengths of the response-to-selection vectors, the larger vector over the smaller. This ratio can also be expressed in terms of the absolute value of a difference in logarithms,

9

Here, is the _j_th of n elements of the response-to-selection vector associated with the _i_th matrix. This metric measures the standardized mean difference in the magnitude of the response to selection produced by G because of the same random selection vector. It was averaged across many random selection vectors and many pairs of G matrices for each time interval. I calculated means and conducted hypothesis testing of M as described for θ, above.

θ and M were used to assess the evolutionary dynamics of G because, given two G matrices: (1) θ provides a measure of the expected difference in the direction of the response to selection induced by differences between matrices (Fig. 1); and (2) M provides a measure of the expected difference in the magnitude of the response to selection induced by differences between matrices. Thus θ and M facilitate evaluating the potential evolutionary importance of differences in G. We might consider differences in G to be of greater potential evolutionary consequence, specifically in terms of the response-to-selection equation (eq. 1), if they induce a large θ or M. That is not to say that other changes in G are biologically unimportant or uninteresting, but rather that those associated with a large θ or M are more likely to affect the effectiveness of the evolutionary inferences based on the multivariate response-to-selection equation (eq. 1).

Results

STABILITY OF G

The mean total size (Σ) and eccentricity (ɛ) of G differed substantially among simulation conditions (Table 1). Σ was maximized when correlational mutation and selection was aligned and high [scenario (_e_: _r_ω= 0.9, _r_μ= 0.9); Table 1], and minimized when the correlation coefficients of mutation and selection had opposite sign (f: _r_ω= 0.75, _r_μ=−0.5). Mean ɛ was maximized (i.e., the eccentricity was minimized) near ɛ= 0.75 when mutation and selection were both uncorrelated ([a: _r_ω= 0, _r_μ= 0]; Table 1) and when mutation and selection were both correlated but with opposite sign ([f: _r_ω= 0.75, _r_μ=−0.5]; Table 1), and minimized (eccentricity maximized) for high _r_ω and _r_μ (e).

1

The influence of correlational natural selection (rω) and pleiotropic mutation (rμ) on the stability of the eigenstructure of G. Following Jones et al. (2003), the following symbols are used: rg : genetic correlation; λ1 : primary eigenvalue of G; λ2 : secondary eigenvalue of G; Σ : total size of G; ɛ : eccentricity of G; ϕ : orientation of the primary eigenvector of G. Δ_x_ is the mean single generation change in parameter x, divided by the value of x, except for in the case of ΔG12 and Δϕ, for which the raw value is given. Parameter values not shown in the table are: _N_= 1000, ω11=ω22= 49, α12=α22= 0.05, μ= 0.0002, and _E_11=_E_22= 1.0.

1

The influence of correlational natural selection (rω) and pleiotropic mutation (rμ) on the stability of the eigenstructure of G. Following Jones et al. (2003), the following symbols are used: rg : genetic correlation; λ1 : primary eigenvalue of G; λ2 : secondary eigenvalue of G; Σ : total size of G; ɛ : eccentricity of G; ϕ : orientation of the primary eigenvector of G. Δ_x_ is the mean single generation change in parameter x, divided by the value of x, except for in the case of ΔG12 and Δϕ, for which the raw value is given. Parameter values not shown in the table are: _N_= 1000, ω11=ω22= 49, α12=α22= 0.05, μ= 0.0002, and _E_11=_E_22= 1.0.

The stability of Σ and ɛ (ΔΣ and Δɛ) was similar for all mutation–selection scenarios, whereas stability of ϕ (Δϕ) differed considerably among simulation conditions. Stability of ϕ was low under scenarios (a: _r_ω= 0, _r_μ= 0) and (f: _r_ω= 0.75, _r_μ=−0.5), but was improved by increased positive _r_ω and _r_μ (Table 1). Aside from scenarios (a) and (f), in which ϕ fluctuates wildly and thus for which the mean orientation is meaningless, the mean of ϕ varied little among simulation conditions (Table 1).

When matrix evolutionary consistency was evaluated in terms of the predictability of the response to selection by serial random skewers, slightly different findings were revealed (Fig. 2A). Evolutionary consistency was not persistently significantly different from the mutation–selection scenario (a: _r_ω= 0, _r_μ= 0; constant, uncorrelated selection and mutation) in mutation–selection scenarios (b: _r_ω= 0.75, _r_μ= 0), (c: _r_ω= 0, _r_μ= 0.5), and (f: _r_ω= 0.75, _r_μ=−0.5) (Fig. 2B). Consistency of the response to selection was mildly but significantly improved under the mutation–selection scenario (d: _r_ω= 0.75, _r_μ= 0.5), and even more greatly increased under strong correlational mutation and selection (i.e., is decreased; [e: _r_ω= 0.9, _r_μ= 0.9]; Figs. 2A, B).

2

Evolutionary consistency of G for constant mutational and selectional correlation (rμ and rω). (A) Mean angle in degrees, , between the response vectors to random selection vectors (random skewers) applied to all pairs of G matrices over all time intervals of each length, Δ_t_. The rate at which increases with Δ_t_ is a measure of the consistency of G. Larger implies less consistent G. (B) The t_-statistic for θ_t,0, the angle between the response vectors to random selection gradients at times t and _t_= 0, for each value of t > 0, across 20 replicate runs. Significantly greater consistency than the null mutation-selection scenario (a) is indicated if t consistently exceeds the horizontal dashed line at t_0.05(df = 38) = 2.02, and significant inconsistency relative to the null scenario (a) is indicated when t < −_t_0.05(df = 38) =−2.02. (C) Mean magnitudinal similarity, M, measured as the ratio of the lengths of the response-to-selection vectors from random skewers, the larger vector over the smaller, between all pairs of G matrices over all time intervals, Δ_t. (D) The _t_-statistic for _M_t,0. N was invariant at _N_= 1000, ω11=ω22= 49, α12=α22= 0.05, μ= 0.0002, and _E_11=_E_22= 1.0 for all simulations. Other parameters are listed in the figure legend. See text for additional explanation.

Matrix consistency evaluated by similarity in the magnitude of the response to selection, M, was unaffected by the mutation–selection scenario used in the simulation (Figs. 2C, D). An equilibrium mean value of M for a large Δ_t_ of approximately 1.2 was approached asymptotically at similar rates by the different mutation–selection scenarios (Fig. 2C).

In Figures 2A and C, it is important to note that for large values of Δ_t_ many fewer comparisons were used to estimate and M than for small Δ_t_ (only 20 for the largest Δ_t_). Consequently, less weight should be given to the extreme right sides of Figures 2A and C; however this does not affect Figures 2B and D for which 20 comparisons (only one per run) are calculated at each time-step.

INSTABILITY OF G

Under a fluctuating coefficient of correlational selection, an increase in the rate of change of _r_ω, σ2(_r_ω), was associated with a decrease in the magnitude of the diagonal elements of G and thus of the overall size of G (Table 2). The mean value of ɛ was low (i.e., matrix eccentricity was high) for a slow rate of change of _r_ω, and high for high σ2(_r_ω) and σ2(_r_ω) = 0. Mean ϕ was relatively unaffected (Table 2; compare (a) and (f) vs. other scenarios in Table 1). The stability of the total size (ΔΣ) and eccentricity (Δɛ) were unaffected by fluctuating _r_ω; however, the stability of G (measured by Δϕ) was highest under slowly changing _r_ω and lowest under unchanging or rapidly changing _r_ω (Table 2).

2

The influence of changing values of the correlation coefficients of correlated natural selection (_r_ω) and pleiotropic mutation (_r_μ) on the stability of the eigenstructure of G. σ2(r) is the intergenerational variance in r. All other column heads are the same as in Table 1. All parameter values not shown are as in Table 1.

2

The influence of changing values of the correlation coefficients of correlated natural selection (_r_ω) and pleiotropic mutation (_r_μ) on the stability of the eigenstructure of G. σ2(r) is the intergenerational variance in r. All other column heads are the same as in Table 1. All parameter values not shown are as in Table 1.

For fluctuating correlational mutation, neither the diagonal elements of G nor Σ were substantially affected. The mean value of ɛ was lowest for slowly changing _r_μ and highest for rapidly changing _r_μ or unchanging _r_μ= 0. The mean orientation, ϕ, was relatively unaffected by the rate of change of _r_μ (Table 2; see comment above). As with fluctuating _r_ω, the stability of the size (ΔΣ) and eccentricity (Δɛ) were unaffected by fluctuating _r_μ; however, the stability of the orientation of G (measured by Δϕ) was highest under slowly changing _r_μ and lowest under unchanging or rapidly changing _r_μ (Table 2).

Under fluctuating population size with _r_ω=_r_μ= 0 and fixed, the mean value of Σ decreased with increasing σ2(N), whereas mean ɛ and ϕ were largely unaffected (Table 3). The stability of the eigenvalues and Σ were lowest for high σ2(N) (Table 3).

3

The influence of changing values of the effective population size (N) on the stability of the eigenstructure of G. σ2(N) is the intergenerational variance in N. All other column heads are the same as in Table 1. _r_ω=_r_μ= 0. All other parameter values not shown are as in Table 1.

3

The influence of changing values of the effective population size (N) on the stability of the eigenstructure of G. σ2(N) is the intergenerational variance in N. All other column heads are the same as in Table 1. _r_ω=_r_μ= 0. All other parameter values not shown are as in Table 1.

These results suggest that the eigenstructure of G is most stable when _r_ω and _r_μ change at a low rate and least stable under unchanging or rapidly changing _r_ω and _r_μ (Table 2). The stability of the total size and eigenvalues of G is weakly compromised by increasing σ2(N) (Table 3).

When I analyzed the simulation results using serial random skewers, G was also fairly consistent in terms of the predictability of the response to selection when _r_ω and _r_μ changed at a low rate. However, the response of G to natural selection was even more consistent if _r_ω and _r_μ were constant and zero, or if they changed rapidly (Fig. 3A: scenario [j: σ2(_r_μ) = 0.1]; Fig. 4A: scenario [n: σ2(_r_ω) = 0.1]). Evolutionary consistency of G was also significantly decreased relative to the scenario of constant correlational selection and mutation only under conditions of _r_ω and _r_μ changing at an intermediate rate (σ2(r) = 0.001; scenarios (h) and (l); Figs. 3B, 4B), and not under other conditions.

3

Evolutionary consistency of G under scenarios involving a changing coefficient of correlational selection, rω. rω changes as a random walk with variance σ2(rω). (A) Mean angle in degrees, , between the response vectors to random skewers as in Figure 1. (B) The t_-statistic for θ_t,0 for each value of t > 0. Significant consistency and inconsistency is indicated by the horizontal dashed lines as in Figure 1. (C) Mean magnitudinal similarity, M, as in Figure 1. d) The _t_-statistic for M t,0. Selection scenarios are listed in the figure legend. rμ is held fixed at _rμ_= 0. All other parameters are as in Figure 1. See text for additional explanation.

4

Evolutionary consistency of G under scenarios involving a changing coefficient of correlational mutation, rμ. rμ changes as a random walk with variance σ2(rμ). (A) Mean angle in degrees, , between the response vectors to random skewers as in Figure 1. (B) The t_-statistic for θ_t,0 for each value of t > 0. Significant consistency and inconsistency is indicated by the horizontal dashed lines as in Figure 1. (C) Mean magnitudinal similarity, M, as in Figure 1. (D) The _t_-statistic for M t,0. Mutation scenarios are listed in the figure legend. rω is held fixed at _rω_= 0. All other parameters are as in Figure 1. See text for additional explanation.

Also in contrast to the result derived under constant correlational mutation and selection, for conditions of changing _r_ω and _r_μ the magnitude of the response to selection is severely affected under some circumstances. In particular, M is highest and most significantly inconsistent when compared to mutation–selection scenario (a: _r_ω=_r_μ=σ2(_r_ω) =σ2(_r_μ) = 0; constant correlational selection and mutation) under conditions in which _r_ω and _r_μ change at an intermediate rate [σ2(_r_) = 0.001; scenarios (_h_) and (_l_); Figs. 3C, D, 4C, D].

When N changes over time, the evolutionary consistency of G also depends on the rate of change in N. In particular, fluctuating population size, in which σ2(N) = 1000, severely increased inconsistency of M (Figs. 5C, D). Higher or lower rates of change of N induced less instability in M (Figs. 5C, D). There is also little evidence for a simple monotonic relationship between evolutionary consistency and the effective population size (Ne)—calculated as the harmonic mean of N ( Appendix B).

5

Evolutionary consistency of G under scenarios of changing effective population size, N. N changes as a random walk with variance σ2(N). (A) Mean angle in degrees, , between the response vectors to random skewers as in Figure 1. (B) The t_-statistic for θ_t,0 for each value of t > 0. Significant consistency and inconsistency is indicated by the horizontal dashed lines as in Figure 1. (C) Mean magnitudinal similarity, M, as in Figure 1. (D) The _t_-statistic for _M_t,0. rμ and rω were held fixed at _rμ_=_rω_= 0. All other parameters are as in Figure 1. See text for additional explanation.

The same disclaimer applying to the far right Figures 2A and C, which is that many fewer comparisons are used to estimate and M for large Δ_t_ than small, also applies to Figures 3A and C, Figures 4A and C, and Figures 5A and C. Consequently, less weight should be accorded to the far right of each panel figure.

Discussion

IS G STABILIZED BY CORRELATED SELECTION AND MUTATION?

The G matrix plays a central role in modern evolutionary theory because of its importance in determining the multivariate response to selection (Lande 1979; Lande and Arnold 1983). However, attempts to extend the equations of quantitative genetics to predict multivariate evolution over macroevolutionary time scales rely on a highly stable or predictably evolving G matrix (Arnold et al. 2001). There are no theoretical reasons to presuppose that G is constant: theory predicts it to be stable in some circumstances and unstable in others (Turelli 1988; Roff 2000). Empirical findings that G may be stable over some time scales (reviewed in Roff 1997), suggest that the reconstruction of past selection and the prediction of morphological diversification based on quantitative genetic models may not be hopeless ventures (Lande 1979; Arnold et al. 2001). Still, there is a paucity of empirical studies in which G from more than a few related taxa are compared (but see Bégin and Roff 2004), and there is little theory to suggest the factors likely to result in high stability or instability of G.

In this paper, I provide insight into the impact of selection, mutation, and effective population size on the evolutionary dynamics of G. In addition to the eigenstructure comparison method (previously employed; see Jones et al. 2003), I also examined the evolutionary consistency of G using the vector correlation approach of serial random skewers (Cheverud et al. 1983; Cheverud 1996). In this latter method, G matrix similarity is assessed explicitly in terms of the predictability of the response to selection given random selection vectors, herein termed “evolutionary consistency.”

I found that the evolutionary consistency of the response to selection was fairly high across a range of conditions and even under circumstances in which the eigenstructure of G was unstable (Jones et al. 2003). This does not in any way negate the importance of the G matrix for evolutionary inference—to the contrary, it implies that G may be sufficiently consistent for evolutionary inference so long as the conditions of correlational selection and mutation are stable over time, irrespective of the particular value of the coefficients of correlational selection and mutation. This conclusion is quite compatible with the findings of Jones et al. (2004) in which the reconstruction of past selection is effective when correlational selection and mutation are constant over time—regardless of the actual value of the correlation coefficient of correlational selection or mutation.

Evolutionary consistency is suggested by the small mean deviation from collinearity ( always < 7° even for large Δ_t_— which corresponds to a very high vector correlation of _r_= cos() ≈ 0.99) between the response-to-selection vectors among all G matrices evolving under conditions of constant correlational mutation and selection (regardless of the actual values of _r_ω and _r_μ). This result implies that the evolutionary importance of the differences among matrices may be small (Fig. 2A). Furthermore, I found negligible and statistically insignificant differences in the evolutionary consistency of G among several of the mutation–selection scenarios (Figs. 2B, D). Nonetheless, although I find that neither moderately correlated selection nor moderately correlated mutation in isolation significantly improves the evolutionary consistency of G as assessed by serial random skewers, the increased stability of G induced by correlational selection and mutation enhances the consistency of G when both are high (Table 1; Figs. 2A, B).

A previous study identified factors that influence the stability of the eigenstructure of G (Jones et al. 2003). I found some instances in which although aspects of the eigenstructure of G were destabilized, G nonetheless produced a consistent response to selection (was evolutionarily consistent). This result is predicted by Jones et al. (2003) who point out that the same conditions that result in highly variable G matrix orientation also result in low eccentricity (ɛ close to 1.0; p. 1758). This superficial discordance between eigenstructure stability and evolutionary consistency will generally be found when no eigenvector of G dominates, as is the case when matrix eccentricity is low (corresponding to high values for ɛ). If no eigenvector dominates G, then the direction of the response to selection is not constrained by G, and instability in the orientation of G has little bearing on the response to natural selection.

INSTABILITY OF G: AN “INTERMEDIATE DISTURBANCE” HYPOTHESIS

I also found that although G may be evolutionarily consistent over a broad range of conditions, it can be rendered evolutionarily inconsistent by coefficients of correlational mutation and selection that change over time, particularly if the rate of change is intermediate. If the correlation coefficient changes too slowly, it will not change substantially during the course of the simulation and thus will fail to induce any inconsistency in G. By contrast, if r changes sufficiently rapidly, the covariance of r from one generation to the next is so low that the system evolves as if mutation or selection were uncorrelated.

This result would not be obtained by evaluating the eigenstructure of G. In particular, when evaluated as the cross-generational change in the orientation of G, stability is highest for the slowest changing correlation coefficient of correlational selection or mutation and lowest for the most rapidly changing and unchanging correlations. This is contrary to the result from serial random skewers in which an intermediate rate of change renders G least evolutionarily consistent. Because serial random skewers evaluates the evolutionary dynamics of G explicitly in terms of the response to natural selection, my findings suggest that, holding all else constant, evolutionary analyses based on the equations of theoretical quantitative genetics will be most severely affected by the correlation coefficients of correlational mutation and selection that change at an intermediate rate.

I also found a similarly nonmonotonic relationship between the rate of fluctuation in population size and the evolutionary consistency of G. When population size did not change at all, G was most evolutionarily consistent. G was least consistent when effective size changed at an intermediate rate. Unlike conditions of fluctuating correlational selection and mutation, in which both the direction and magnitude of the response to selection are affected, under fluctuating effective population size only the magnitude of the response to selection is inconsistent over time. This is because effective population size primarily affects the overall size of G, rather than its shape (see Jones et al. 2003).

WHAT IS MEANT BY AN INTERMEDIATE RATE OF CHANGE?

One central result of this study is that the evolutionary consistency of G evaluated by the response to selection does not increase or decrease monotonically with the rate of change of the coefficients of correlational selection or mutation, or the effective population size. However, it is important to keep in mind that the specific conditions of this study that were responsible for substantial instability in G cannot be easily translated literally to the empirical realm.

Two factors in particular prohibit this kind of extrapolation. (1) The mutation rate used in this study (μ= 0.0002 mutations/allele generation) was used following Jones et al. (2003, 2004) and to facilitate computational expediency. This rate may be unrealistically high for a single quantitative trait locus, but may be realistic if considered to be the mutation rate for several physically linked loci affecting the same trait (see Jones et al. 2003[p. 1753]). It is possible that similar matrix inconsistency would be induced by lower σ2(_r_μ) if the value of μ were also decreased and vice versa. (2) The time scale over which the evolution of G is studied in these analyses is quite short by macroevolutionary standards, yet is longer than almost any ecological study (2000 generations). It is possible that lower or higher σ2(r) and σ2(N) would induce comparable matrix dynamics over longer and shorter evolutionary time scales respectively. Figure 6 compares the evolutionary consistency of G when σ2(_r_μ) = 0.001 and μ= 0.0002 (model [_l_], above) to that when both parameters are decreased by a factor of 10 and , respectively (which is the same as decreasing σ and μ by the same factor, : σ2(_r_μ) = 0.0001 and μ≈ 0.000063). They exhibit almost identical evolutionary dynamics and neither is similar to σ2(_r_μ) = 0.0001 and μ= 0.0002 (Fig. 6). This shows that simultaneously changing some parameters of the mutation–selection scenario can sometimes have compensatory effects on the consistency of G.

6

Evolutionary consistency of G under several conditions of fluctuating, rμ. rμ changes as a random walk with variance σ2(rμ). σ2(rμ) = 0.001 with mutation rate μ= 0.0002 is compared to σ2(rμ) = 0.0001 with mutation rate μ= 0.0002 × 10−1/2, which is the same as decreasing σ and μ by the same factor (). σ2(rμ) = 0.0001 and μ= 0.0002 are also shown. All other parameters are the same as in Figure 4.

Although the nature of the compensatory effect is an empirical observation of this study, its discovery is not surprising. This is because for a given starting value of _r_μ, the number of pleiotropic mutations that arise with expected correlation _r_μ is determined by the relationship between μ and σ. If μ is high relative to σ, then many pleiotropic mutations will take place before there has been substantive change in _r_μ. By contrast, if μ is low relative to σ, then _r_μ is likely to change too rapidly for there to be substantial accumulation of pleiotropic mutations with correlation _r_μ.

In spite of the limitations on the quantitative predictions that can be made from this study, it may be useful to consider whether empirical values for σ2(rω), σ2(_r_μ), and σ2(N) vary on a range equal to or greater than that explored herein. If so, then the intermediate minimum in consistency will be experienced by empiricists. If not, then the effective empirical relationship between σ2(rω), σ2(_r_μ), and σ2(N) may be monotonic and not hump-shaped as my results imply.

Compared to σ2(rω) and σ2(_r_μ), σ2(N) is very easy to address. Empirical studies have shown both persistent population stability and wide population fluctuation over ecological time scales (e.g., Hanski et al. 1991; Hörnfeldt et al. 2005; Venkateswaran and Parthasarathy 2005; Elias et al. 2006; Pellet et al. 2006). These fluctuations can exist across several orders of magnitude (e.g., Hörnfeldt et al. 2005), which exceeds the range of N for this study. Furthermore, the rate of population fluctuation can vary both among populations of a single species (e.g., Hanski et al. 1991) and over time within a single population (e.g., Hörnfeldt et al. 2005). In addition, genetic studies of mtDNA have shown that present-day census sizes frequently exceed long-term population genetic effective sizes (e.g., Avise 2000), an observation consistent with populations that fluctuate in size over time (Hartl and Clark 1989).

Empirical values for σ2(rω) and σ2(_r_μ) are much more difficult to find. A recent review (Kingsolver et al. 2001) observed that many studies of natural selection have insufficient power to detect quadratic selection (of which correlational selection is a component) or neglect to even estimate it. Fluctuating linear natural selection is well known, however, in some cases involving natural selection that reverses sign (e.g., Grant and Grant 1989; Losos et al. 2006). Fluctuating correlational selection may be as common but is poorly documented.

Similarly, the temporal stability of _r_μ is poorly known. One study designed to experimentally estimate M, the mutational variance–covariance matrix of which _r_μ is a component, showed that although M in different experimental lines were correlated, there were also significant differences in M among lines (Camara and Pigliucci 1999; Camara et al. 2000). How this result applies to the temporal stability of _r_μ is unclear; however, if M differs among experimental isogenic lines, it is also highly plausible that M and _r_μ vary over time.

In all, it is clear than in some species fluctuating N is sufficiently large to drive significant inconsistency in G. However, in this study variable N affected primarily the evolutionary consistency of the magnitude of the response to selection, impacting little the evolutionary consistency of the direction of the response to selection (Fig. 5). By contrast, fluctuations in rω and _r_μ severely influenced the evolutionary consistency of G, but it is unclear over what range rω and _r_μ naturally vary.

EXTENSION TO A MULTIVARIATE G MATRIX

Following previous authors (e.g., Jones et al. 2003, 2004) I use simple bivariate quantitative genetic simulations in this study. Most empirical quantitative genetic experiments estimate G for three or more traits. It is possible that the results of this study apply equally to the evolution of higher dimensional G matrices; however, this cannot be assumed.

Multivariate quantitative genetic simulations create many difficulties that do not affect the bivariate case (discussed in Appendix C). However, it can be shown in a very limited set of numerical simulations that the results for four traits are consistent with those obtained herein for two traits ( Appendix C). In particular, evolutionary consistency is lowest when the matrix of correlational pleiotropic mutation is unstable.

Although my multivariate simulations are quite limited in scope, to my knowledge this is the first study presenting individual-based multivariate quantitative genetic simulations of the evolution of the G matrix for more than two phenotypic traits. Future similar studies might consider a more expansive set of simulations of the evolution of the multivariate G matrix.

CONCLUSIONS

In this study the evolutionary dynamics of G in computer simulations were evaluated using eigenstructure analysis and an approach utilizing the response to random selection gradients. This “evolutionary consistency” provides a measure of the usefulness of G for evolutionary inferences that rely on a stable covariance structure.

Evolutionary consistency was high under some simulation conditions in which the eigenstructure of G was unstable. This was because of the fact that under many of the simulation conditions explored in this and previous studies (e.g., Jones et al. 2003), the circumstances in which the orientation of G is least stable coincide with those in which matrix eccentricity is low. Consequently, instability of the matrix orientation in these simulations does not affect the consistency of G.

Under conditions of fluctuating correlational selection and mutation, and fluctuating population size, G was most inconsistent when fluctuations occurred at an intermediate rate—that is, matrix inconsistency does not increase monotonically with increasing rates of fluctuation in the coefficients of correlational selection and mutation, and the population size. Population fluctuation over ecological time scales is well known. However, at the present time, estimates of the temporal stability of the coefficients of correlational selection and mutation from natural populations are lacking. Consequently, although it is highly plausible that the rates of fluctuation in these parameters vary on the range explored by this study, concrete empirical recommendations await estimates of the temporal stability of correlational selection and mutation.

Associate Editor: S. Steppan

ACKNOWLEDGMENTS

This work was supported in part by a grant from the National Science Foundation (DEB- 0519777). I thank L. Harmon, B. Langerhans, J. Losos, S. Steppan, and several anonymous reviewers for insightful comments on earlier versions of the paper.

LITERATURE CITED

Arnold

,

S. J.

1992

.

Constraints on phenotypic evolution

.

Am. Nat

.

140

:

S85

–

S107

.

Arnold

,

S. J.

,

M. E.

Pfrender

, and

A. G.

Jones

.

2001

.

The adaptive landscape as a conceptual bridge between micro- and macroevolution

.

Genetica

112-113

:

9

–

32

.

Avise

,

J. C.

2000

.

Phylogeography: the history and formation of species

.

Harvard Univ. Press

, Cambridge MA.

Barton

,

N. H.

, and

M.

Turelli

.

1987

.

Adaptive landscapes, genetic distance and the evolution of quantitative characters

.

Genet. Res

.

49

:

157

–

173

.

Bégin

,

M.

, and

D. A.

Roff

.

2004

.

From micro- to macroevolution through quantitative genetic variation: positive evidence from field crickets

.

Evolution

58

:

2287

–

2304

.

Bodmer

,

W. F.

1965

.

Differential fertility in population genetics models

.

Genetics

51

:

411

–

424

.

Bürger

,

R.

, and

R.

Lande

.

1994

.

On the distribution of the mean and variance of a quantitative trait under mutation-selection-drift balance

.

Genetics

138

:

901

–

912

.

Bürger

,

R.

,

G. P.

Wagner

, and

F.

Stettinger

.

1989

.

How much heritable variation can be maintained in finite populations by mutation-selection balance?

Evolution

43

:

1748

–

1766

.

Camara

,

M. D.

, and

M.

Pigliucci

.

1999

.

Mutational contributions to genetic variance-covariance matrices: an experimental approach using induced mutations in Arabidopsis thaliana

.

Evolution

53

:

1692

–

1703

.

Camara

,

M. D.

,

C. A.

Ancell

, and

M.

Pigliucci

.

2000

.

Induced mutations: a novel tool to study phenotypic integration and evolutionary constraints in Arabidopsis thaliana

.

Evol. Ecol. Res

.

2

:

1009

–

1029

.

Cheverud

,

J. M.

1996

.

Quantitative genetic analysis of cranial morphology in the cotton-top (Saguinus oedipus) and saddle-back (S. fuscicollis) tamarins

.

J. Evol. Biol

.

9

:

5

–

42

.

Cheverud

,

J. M.

,

J. J.

Rutledge

, and

W. R.

Atchley

.

1983

.

Quantitative genetics of development: genetic correlations among age-specific trait values and the evolution of ontogeny

.

Evolution

37

:

895

–

905

.

Crow

,

J. F.

, and

M.

Kimura

.

1964

. The theory of genetic loads. Pp.

495

–

505

in

S. J.

Geerts

, ed.

Proceedings of the XI international congress of genetics

.

Pergamon

, Oxford , U.K.

Elias

,

S. P.

,

J. W.

Witham

, and

M. L.

Hunter

,Jr.

2006

.

A cyclic red-backed vole (Clethrionomys gapperi) population and seedfall over 22 years in Maine

.

J. Mammal

.

87

:

440

–

445

.

Grant

,

B. R.

, and

P. R.

Grant

.

1989

.

Natural selection in a population of Darwin's finches

.

Am. Nat

.

133

:

377

–

393

.

Hanski

,

I.

,

L.

Hansson

, and

H.

Henttonen

.

1991

.

Specialist predators, generalist predators, and the microtine rodent cycle

.

J. Anim. Ecol

.

60

:

353

–

367

.

Hartl

,

D. L.

, and

A. G.

Clark

.

1989

.

Principles of population genetics

. 2nd ed.

Sinauer

, Sunderland, MA.

Hörnfeldt

,

B.

,

T.

Hipkiss

, and

U.

Eklund

.

2005

.

Fading out of vole and predator cycles?

Proc. R. Soc. Lond. B

272

:

2045

–

2049

.

Jones

,

A. G.

,

S. J.

Arnold

, and

R.

Bürger

.

2003

.

Stability of the G-matrix in a population experiencing pleiotropic mutation, stabilizing selection, and genetic drift

.

Evolution

57

:

1747

–

1760

.

Jones

,

A. G.

,

S. J.

Arnold

, and

R.

Bürger

.

2004

.

Evolution and stability of the G-matrix on a landscape with a moving optimum

.

Evolution

58

:

1639

–

1654

.

Kimura

,

M.

1965

.

A stochastic model concerning the maintenance of genetic variability in quantitative characters

.

Proc. Nat. Acad. Sci

.

54

:

731

–

736

.

Kingsolver

,

J. G.

,

H. E.

Hoekstra

,

J. M.

Hoekstra

,

D.

Berrigan

,

S. N.

Vignieri

,

C. E.

Hill

,

A.

Hoang

,

P.

Gibert

, and

P.

Beerli

.

2001

.

The strength of phenotypic selection in natural populations

.

Am. Nat

.

157

:

245

–

261

.

Lande

,

R.

1979

.

Quantitative genetic analysis of multivariate evolution, applied to brain:body size allometry

.

Evolution

33

:

402

–

416

.

Lande

,

R.

, and

S. J.

Arnold

.

1983

.

The measurement of selection on correlated characters

.

Evolution

37

:

1210

–

1226

.

Losos

,

J. B.

,

T. W.

Schoener

,

R. B.

Langerhans

, and

D. A.

Spiller

.

2006

.

Rapid temporal reversal in predator-driven natural selection

.

Science

314

:

1111

.

Lynch

,

M.

, and

B.

Walsh

.

1998

.

Genetics and analysis of quantitative traits

.

Sinauer Associates

, Sunderland, MA.

Marroig

,

G.

, and

J. M.

Cheverud

.

2001

.

A comparison of phenotypic variation and covariation patterns and the role of phylogeny, ecology, and ontogeny during cranial evolution of new world monkeys

.

Evolution

55

:

2576

–

2600

.

Pellet

,

J.

,

B. R.

Schmidt

,

F.

Fivaz

,

N.

Perrin

, and

K.

Grossenbacher

.

2006

.

Density, climate and varying return points: an analysis of long-term population fluctuations in the threatened European tree frog

.

Oecol

.

149

:

65

–

71

.

Reeve

,

J. P.

2000

.

Predicting long-term response to selection

.

Genet. Res

.

75

:

83

–

94

.

Revell

,

L. J.

,

L. J.

Harmon

,

R. B.

Langerhans

, and

J. J.

Kolbe

.

2007

.

A phylogenetic approach to determining the importance of constraint on phenotypic evolution in the neotropical lizard Anolis cristatellus

.

Evol. Ecol. Res

.

9

:

261

–

282

.

Roff

,

D. A.

1997

.

Evolutionary quantitative genetics

.

Chapman & Hall

, New York.

Roff

,

D. A.

.

2000

.

The evolution of the G matrix: selection or drift?

Heredity

84

:

135

–

142

.

Roff

,

D. A.

, and

T. A.

Mousseau

.

1999

.

Does natural selection alter genetic architecture? An analysis of quantitative genetic variation among populations of Allenomobius socius and A. fasciatus

.

J. Evol. Biol

.

12

:

361

–

369

.

Shaw

,

F. H.

,

R. G.

Shaw

,

G. S.

Wilkinson

, and

M.

Turelli

.

1995

.

Changes in genetic variances and covariances: G whiz!

Evolution

49

:

1260

–

1267

.

Steppan

,

S. J.

,

P. C.

Phillips

, and

D.

Houle

.

2002

.

Comparative quantitative genetics: evolution of the G matrix

.

Trends Ecol. Evol

.

17

:

320

–

327

.

Turelli

,

M.

1988

.

Phenotypic evolution, constant covariances, and the maintenance of additive variance

.

Evolution

42

:

1342

–

1347

.

Venkateswaran

,

R.

, and

N.

Parthasarathy

.

2005

.

Tree population changes in a tropical dry evergreen forest of south India over a decade (1992-2002)

.

Biodivers. Conserv

.

14

:

1335

–

1344

.

Appendix

Appendix A. Fluctuating coefficients of correlational mutation and selection and the evolutionary dynamics of G.

To study the evolutionary dynamics of G when the coefficients of correlational selection and mutation varied over time, I simulated changing _r_ω and _r_μ under a Brownian motion process with bounds −1.0 and 1.0, and variance σ2(r).

Figure A1 provides example time series for _r_ω or _r_μ over 2000 generations using the same values for σ2(r) as in the main text.

A1

Example time series for the correlation coefficient, r, of correlational mutation or natural selection. σ2(r) is the expected variance across generations. r was evolved according to a Brownian motion process with bounds −1.0 ≤_r_≤ 1.0.

Appendix B. Effective population size and the evolutionary dynamics of G when the population size varies over time.

To determine the relationship between the effective population size (_N_e) and the evolutionary dynamics of G, I first calculated _N_e from each simulation in which the population size N was varied over time. I calculated _N_e as the harmonic mean of N over the course of the simulation (Hartl and Clark 1989). In these simulations, as discussed in the main text, I varied N according to a Brownian motion process with N(0) = 1000, bounds 20 ≤_N_≤ 2000, and variance σ2(N).

Table B1 lists the mean effective population sizes for each value of σ2(N), along with the standard deviation of ,. I obtained by calculating the arithmetic mean across runs of the harmonic mean of N within runs for a given simulation scenario.

B1

Mean effective population size for various rates of change in N. Mean effective size for each simulation condition was calculated as the arithmetic mean across runs of the harmonic mean within runs of N for each of 20 runs. was calculated as the square root of the variance among runs in Ne.

B1

Mean effective population size for various rates of change in N. Mean effective size for each simulation condition was calculated as the arithmetic mean across runs of the harmonic mean within runs of N for each of 20 runs. was calculated as the square root of the variance among runs in Ne.

To provide examples of the sort of population size fluctuation that results from different σ2(N), I simulated N over time for 2000 generations under each of five σ2(N) [including the trivial case of σ2(N) = 0]. Quantitative genetic simulations were performed over 12,000 generations, see main text, although evolutionary dynamics were only assessed in the final 2000 generations. These values of σ2(N) correspond to mutation-selection scenarios (a) and (o) through (r) in the main text. Figure B1 plots N over time under each scenario.

B1

Example time series for the population size, N, under various conditions of fluctuating population size. σ2(N) is the expected variance across generations. N was changed according to a Brownian motion process with bounds 20 ≤_N_≤ 2000.

Appendix C. Evolutionary dynamics of G matrices containing more than two traits.

To explore the generality of the findings presented in this study in which, following Jones et al. (2003, 2004), I investigate only the evolutionary dynamics of a bivariate G matrix, I simulated the evolution of G containing four traits.

Multivariate simulations are subject to several limitations that do not affect the bivariate case. Firstly, multivariate simulations are much more computationally intensive. Second, the mutation-selection parameter space dimensionality increases in proportion to the square of the number of traits, _m_[in particular the number of parameters of correlational selection and mutation is equal to m(m + 1)]. Third, the coefficients of correlational selection or mutation cannot be arbitrarily specified and cannot be arbitrarily changed over time. This is because randomly chosen elements in a correlation/covariance matrix will probably not satisfy the requirement of positive definiteness—nor will a positive definite matrix whose elements are randomly changed remain positive definite.

Nonetheless, I performed some limited simulations of the evolution of G in four dimensions.

For these analyses, I used a highly similar population genetic model to that used in the bivariate simulations—with some exceptions made for computational reasons. In particular, whereas the number of loci, n, was set to 50 in the bivariate simulations, it was adjusted to 20 for the four-dimensional simulations; in contrast to all bivariate simulations, no selection was simulated; and the mutation rate, μ, was set to 0.0025 at all loci.

I used two different correlation matrices for the mutational effects of pleiotropic mutations. The first, Rμ(1), represents the very simple case in which there are no correlations between the pleiotropic effects of mutations on the four traits. I also used this matrix as a starting point for the case in which Rμ was varied over time:

In the second case, I specified Rμ(2) so as to contain the range of correlation coefficients used for mutation (and selection) in the bivariate simulations. Obviously, the specification of Rμ(2) was constrained by the requirement of positive definiteness:

I also set the mutational variances for the four traits as follows:

I then ran three sets of 20 individual-based, numerical simulations under otherwise similar conditions to those used in the bivariate simulations. In the first two sets of simulations I initialized the correlation matrix of correlational pleiotropic mutation with Rμ(1) and Rμ(2), respectively, and held the mutation matrix constant throughout the simulation. In the third set, I initialized the correlation matrix with Rμ(1) and randomly perturbed it by small steps over time, to simulate fluctuating correlational pleiotropic mutation [henceforward, the conditions of this simulation are referred to as Rμ(3)].

As discussed in a preceding paragraph of this appendix (above), it is not possible to directly, randomly change the elements of a correlation matrix and expect the matrix to retain its necessary property of positive definiteness. However, it is possible to evolve the elements of the Cholesky decomposition of the correlation matrix and then each generation recover the correlation matrix from its Cholesky decomposition matrix, where Rμ=CTC, in which C is the upper diagonal Cholesky decomposition matrix—whose columns have been standardized to have a sum of squares equal to 1.0—and T indicates that a matrix transpose was calculated. This approach is affected by two shortcomings apparent to the author. In particular (1) although the elements of the Cholesky decomposition matrix might be evolved by some model, say Brownian motion, it is not clear under what model the correlation matrix is evolving; and (2) the elements of Rμ evolve at different rates, even though the elements of C may be perturbed at a constant rate. That said, I can empirically verify that if the Cholesky decomposition matrix is perturbed by small steps, then the correlation matrix is perturbed by small steps, and if the Cholesky decomposition matrix is evolved by large steps, then the correlation matrix will also evolve by large steps. In these simulations, I set the variance of the Brownian motion process by which the Cholesky decomposition matrix of the correlation matrix of the effects of pleiotropic mutations was evolved to σ2= 0.001.

I then performed serial random skewers, as described in the main body of this article. Figure C1 shows that G is more consistent when the covariance matrix of correlational pleiotropic mutation is stable over time and less consistent when the mutation matrix changes over time. The value of Rμ also affects the evolutionary consistency of G. G was more consistent if off-diagonal elements of G are nonzero [**Rμ**(2)].

C1

The evolutionary consistency of G for four traits under different conditions. Conditions are: constant correlation matrix of the effects pleiotropic mutations, no correlation among traits [Rμ(1)]; constant correlation matrix of the effects of pleiotropic mutations, correlations among traits [Rμ(2)]; and a correlation matrix of the effects of pleiotropic mutation that was initiated as in Rμ(1), but then evolved over time [Rμ(3)]. Specific parameter values for the simulations are provided in the text.

Because these results are not incompatible with those obtained by way of bivariate simulations, it is highly possible that the findings of my study may extend to the more realistic situation of a G matrix containing three or more traits. However, given the limited nature of the multivariate simulations presented herein, this should be the subject of future more detailed study.

© 2007, Society for the Study of Evolution