Genetics, Vol. 175, 1813-1822, April 2007, Copyright © 2007
doi:10.1534/genetics.106.066530

This Article Abstract Full Text (PDF) All Versions of this Article: genetics.106.066530v1 175/4/1813    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Gu, X. Search for Related Content PubMed PubMed Citation Articles by Gu, X.

Evolutionary Framework for Protein Sequence Evolution and Gene Pleiotropy

Department of Genetics, Development and Cell Biology, Center for Bioinformatics and Biological Statistics, Iowa State University, Ames, IA 50011

¹ Address for correspondence: Department of Genetics, Developmental and Cell Biology, 536 Science II Hall, Iowa State University, Ames, IA 50011.
E-mail: xgu{at}iastate.edu

	ABSTRACT

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
In this article, we develop an evolutionary model for protein sequence evolution. Gene pleiotropy is characterized by K distinct but correlated components (molecular phenotypes) that affect the organismal fitness. These K molecular phenotypes are under stabilizing selection with microadaptation (SM) due to random optima shifts, the SM model. Random coding mutations generate a correlated distribution of K molecular phenotypes. Under this SM model, we further develop a statistical method to estimate the "effective" number of molecular phenotypes (K_e) of the gene. Therefore, for the first time we can empirically evaluate gene pleiotropy from the protein sequence analysis. Case studies of vertebrate proteins indicate that K_e is typically 6–9. We demonstrate that the newly developed SM model of protein evolution may provide a basis for exploring genomic evolution and correlations.

It is desirable to develop a novel framework that can provide an integrated view of protein evolution, but an unsolved key problem is how to deal with the sophisticated relationship between protein sequence and organismal fitness. We recognized that the substantial genotype–phenotype literature can be used to address this issue (e.g., FISHER 1930; LANDE 1980; G. WAGNER 1989; WAXMAN and PECK 1998; A. WAGNER 2000; POON and OTTO 2000; ZHANG and HILL 2003). As mutations are raw materials for evolution, the distribution of effects of mutations, f(s), has been a central issue in the genotype–phenotype study (LYNCH et al. 1999; BATAILLON 2000; SHAW et al. 2002). Under Fisher's geometric model (e.g., FISHER 1930; HARTL and TAUBES 1998; POON and OTTO 2000) or multivariate model (e.g., LANDE 1980; WELCH and WAXMAN 2003; MARTIN and LENORMAND 2006), f(s) can be predicted on the basis of selective and mutational assumptions. The dimensionality in each framework has been interpreted as the number of loci affecting a quantitative trait (reviewed in KINGSOLVER et al. 2001), or the number of affected phenotypic traits of a mutation, estimated from mutational-accumulation experiments (reviewed in ELENA and LENSKI 2003). On the other hand, an ad hoc gamma distribution for f(s) was introduced in the early days of the neutral theory of molecular evolution (OHTA 1973, 1977; KIMURA 1983), which has been found useful to estimate the proportion of neutral or deleterious mutations (e.g., NIELSEN and YANG 2003; PIGANEAU and EYRE-WALKER 2003; EYRE-WALKER et al. 2006).

Here we extend our recent work (GU 2006). The principle behind our present work in modeling protein sequence evolution is as follows: First, multifunctionality of a gene (pleiotropy) is characterized by K distinct components in the fitness (called molecular phenotypes). Second, these K molecular phenotypes are under stabilizing selection (FISHER 1930; WAXMAN and PECK 1998) with microadaptations due to random shifts of the fitness optimum (GILLESPIE 1991). Third, random mutations in the coding region of the gene generate a mutational distribution for K molecular phenotypes. Correlated structures among K molecular phenotypes are taken into account at both selectional and mutational levels. More importantly, we develop a statistical method to estimate the "effective" number of molecular phenotypes (K_e) of the gene. Case studies of vertebrate proteins demonstrate that for the first time we can empirically evaluate gene pleiotropy only from phylogenetic sequence analysis.

	THE MODEL

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Molecular phenotypes and gene pleiotropy:
Pleiotropy of a gene (protein) function can be modeled by multiple (K) molecular phenotypes, denoted by (y₁, ..., y_K). Each molecular phenotype y_i represents a nontrivial component of genetic variation, contrasting with organismal fitness as a result of a specific (yet unknown) biological process (Figure 1). At one extreme, molecular phenotypes could correspond to subcomponents of protein function, regardless of the biological processes. At another extreme, molecular phenotypes could be determined mainly by distinct physiological processes. Since these underlying biological processes are usually intractable, the concept of molecular phenotypes based on fitness effects avoids this difficulty.

View larger version (19K): [in this window] [in a new window] [Download PPT slide]   FIGURE 1.— Schematic presentation for the concept of molecular phenotypes that affect the fitness of an organism.

(1)

where µ is the optimum of the fitness function, and _w is a (positive definite) symmetric matrix characterizing correlated stabilizing selection on K molecular phenotypes. The ith diagonal element measures the strength of stabilizing selection on the ith molecular phenotype, while the ijth nondiagonal element _w,ij measures the correlated stabilizing selection on y_i and y_j.

Let y₀ be the population mean of molecular phenotypes of a gene, which may not necessarily be at the optimum (µ). By definition (FISHER 1930), the coefficient of selection for the molecular phenotype y is given by

(2)

Classical stabilizing selection (TURELLI 1985; WAXMAN and PECK 1998) assumes that the population mean is always at the fitness optimum, i.e., y₀ = µ, leading to purifying selection; that is, (y | y₀ = µ) 0 always holds.

Next, we consider the effects of mutations on molecular phenotypes. A random mutation in the coding region of a gene may affect molecular phenotypes in a correlated fashion (Figure 1). Such a pleiotropic mutational effect can be described by a multivariate normal distribution

(3)

where y_m is the mean mutational effect on molecular phenotypes, and the covariance matrix _m measures the correlated mutational effects. The coefficient of selection (y | y₀, µ) and the distribution of mutational effects p(y) provide a connection between organismal fitness and protein sequence evolution through molecular phenotypes y (Figure 1). Below we discuss the classical stabilizing selection model and the newly developed stabilizing selection with microadaptation (SM) model. See Figure 2 for illustration.

View larger version (9K): [in this window] [in a new window] [Download PPT slide]   FIGURE 2.— (A) Stabilizing selection model with microadaptation in the case of a single molecular phenotype. (B) Distribution of mutational effects on the molecular phenotype. In both cases, we set K = 1 for simplicity.

For mathematical convenience, we assume that the mutational distribution, p(y) in Equation 3, is centered at the population mean, namely, y_m = y₀ = 0. In the following we consider mainly the selection intensity S(y) = 4N_e0(y) = –2N_e(), where N_e is the effective population size. Given Equation 3 of p(y), the distribution of S, f(S), can be uniquely determined by the probability theory (see the APPENDIX). Although the analytical form of f(S) is available only for a few special cases, the moments of S can be obtained readily. For instance, the mean of S can be computed by

(4)

which is always negative.

Stabilizing selection with microadaptation model:
Consider the case when the fitness optimum (µ) of a gene is no longer fixed during the evolution. Rather, µ can be shifted either by environmental changes or by internal physiological perturbations. Each shift of µ results in a microadaptation toward a new fitness optimum (GILLESPIE 1991). Since the direction and strength of the µ-shift are unobservable, we treat µ as a (K-dimensional) random variable. Similar to the above, we assume that the population mean and the mutation mean satisfy y₀ = y_m = 0. Then, under the SM model, the coefficient of correlation defined by Equation 2 can be specified as follows:

(5)

After assuming that µ follows a multivariate normal distribution (µ) with mean 0 and covariance matrix _µ, one can integrate out the unobservable µ by . Therefore, one can show that (y), the coefficient of selection of y averaged over µ-shifts during the evolution, is given by

(6)

where the matrix characterizes the correlated nature of molecular phenotypes in fitness after stabilizing selection and microadaptation. It follows that the selection intensity defined by S(y) = 4N_e(y) is given by

(7)

Apparently, the stabilizing selection model is a special case when _µ = 0 so that U = _w. Moreover, unlike the stabilizing selection model that S(y) is always negative, the selection intensity S(y) under the SM model can be positive when the µ-shift becomes huge. See ANALYTICAL RESULTS for more details.

Evolutionary rate of protein sequence:
Suppose that a protein carrying one (nonsynonymous) mutation has the molecular phenotype vector y. The selection intensity for this mutation, S(y), is given by Equation 7. Hence, the well-known formula for the evolutionary rate (KIMURA 1983) depends on the molecular phenotypes y; that is,

(8)

where v is the mutation rate. At the gene level, molecular phenotypes y are random variables, resulting from random mutations at nucleotide sites (or amino acid residues). Therefore, the evolutionary rate is a random variable that varies among sites. Under the SM model, the mean evolutionary rate of a protein sequence is given by

(9)

given the distribution of mutational effects on y, p(y).

	ANALYTICAL RESULTS

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Here we present analytical results under the SM model, which are useful to study the underlying mechanisms governing protein sequence evolution and to develop a statistical method for estimating important parameters such as K, the number of molecular phenotypes of a gene.

Single molecular phenotype (K = 1):
When K = 1, the fitness function in Equation 1 can be reduced as and the selection coefficient in Equation 6 as , where measures the strength of single stabilizing selection. Similarly, the shift of µ follows a single normal distribution (µ) with mean 0 and . Together, Equation 8 can be simplified as

(10)

Given the distribution of mutational effects , one can show the mean of selection intensity, , as follows:

(11)

Several interesting results need to be mentioned. First, is determined by the mutational variance () and the variance of microadaptation () relative to the coefficient of stabilizing selection (). Hence, the scale of molecular phenotype is not important. Second, if < , and vice versa. In particular, the microadaptation model is reduced to the stabilizing selection model when = 0 (GU 2007). Third, microadaptation reduces the sequence conservation by decreasing the absolute values of , although may remain negative. In other words, microadaptation may increase the evolutionary rate of protein sequence. And fourth, from Equation 10, one can show that S follows a negative gamma distribution with the shape parameter . For instance, when < , S follows a negative gamma distribution

(12)

K molecular phenotypes:
For K 2, the distribution of selection intensity S, f(S), is generally not analytically tractable. Nevertheless, one can derive the analytical forms of moments on the basis of the theory of multiple normal distribution. These results are crucial for estimating K from the sequence data.

Mean and variance of S—general results:
Denote matrix A = U^–1_m = Z – , where and . Briefly, under stabilizing selection, matrix Z describes correlated mutational effects, while describes correlated microadaptations. Together, matrix A characterizes the net effects of correlated mutations on fitness under the SM model. In the APPENDIX (also see MARTIN and LENORMAND 2006 for phenotypic data), we show that the joint effects of all selective, microadaptive, and mutational covariances can be reduced to K eigenvalues of matrix A, ₁, ..., _K. Consequently, the mean and variance of S are given by

(13)

and

(14)

respectively. In short, each _i corresponds to an independent molecular phenotypic direction, on which mutation, microadaptation, and stabilizing selection act with an average effect –2N_ei on the mean selection intensity () of the gene.

Canonical representation of :
Although Equation 13 is mathematically elegant, we wish to find an equivalent representation of that is biologically more interpretable. Due to the arbitrary nature of the original K molecular phenotypes, without loss of generality one may choose a convenient coordinate system. Here we adopt a canonical form of K molecular phenotypes, to define independent phenotypes experiencing stabilizing selection. This leads to a diagonal matrix _w; each diagonal element measures the independent stabilizing selection on the ith (canonical) molecular phenotype. Let be the ith diagonal element of matrix _w or the mutational variance for the ith (canonical) molecular phenotype. We have shown that the canonical presentation for the mean of selection intensity is given by

(15)

where the parameter _i measures the net effect of microadaptation on the ith (canonical) molecular phenotype. In particular, when the mutational effect or microadaptation is independent among the canonical molecular phenotypes, we have . It should be noted that the pleiotropic nature of mutations affects the variance of S, Var(S) (see the APPENDIX), while is not affected, as shown in Equation 15.

Model classification:
As the stabilizing selection acts against nucleotide substitutions (negative selection), microadaptation may provide an opposite force (positive selection) to increase the evolutionary rate. The mean selection intensity given by Equations 13 and 15 can be used to classify the pattern of sequence conservation under the SM model:

1. Strong stabilizing selection with weak microadaptation (SM_w): In this case, the magnitude of µ-shift of molecular phenotypes during the evolution is small, relative to the strength of stabilizing selection. It indicates a dominant purifying selection in the protein sequence evolution. The SM_w model can be defined as having matrix A positive definite, such that each eigenvalue _i > 0 and . In the canonical form of Equation 15, it means that _i < 1 for any i = 1, ... , K. The pure stabilizing selection is a special case of _i = 0 (no microadaptation). Moreover, implies that the ratio of nonsynonymous to synonymous substitutions (d_N/d_S) is <1. Therefore, the SM_w model may describe the general pattern in protein sequence evolution, evidenced by genomewide analyses that have shown d_N/d_S < 1 for most genes.
   
2. Episodic microadaptation under strong stabilizing selection (SM_E): Some genes may have experienced episodic adaptive processes in a few molecular phenotypes, resulting in negative eigenvalues of matrix A (i.e., _i < 0 for some i's, but the mean selection intensity remains negative, ). That is, matrix A is no longer positive definite but a positive trace remains. In the canonical form, some of _i are >1, although . Apparently, the d_N/d_S ratio under the SM_E model may be increased by episodic microadaptations, but d_N/d_S < 1 remains. A special case is the neutral-like behavior of protein sequence evolution when , resulting in . That is, the observed neutral-like evolution may reflect the zero net effect on fitness, canceling each other between episodic microadaptation and stabilizing selection.
   
3. Strong microadaptation under stabilizing selection (SM_s): In a few genes, positive selection and adaptation may dominate the evolution in many molecular phenotypes. In these cases, results in and the ratio d_N/d_S > 1. An extreme case occurs if _i < 0 or _i > 1 (the canonical form) for all i = 1, ..., K. In this case, matrix A is negative definite, indicating the overwhelming adaptation forces for virtually all of the substitutions.
   

The W-model for S-gamma distribution:
Under some strong assumptions, GU (2007) showed that S is distributed as a (negative) gamma. A more general condition for S-gamma distribution exists when stabilizing selection, mutations, and microadaptation have very similar correlation structures among molecular phenotypes. Indeed, when the three corresponding matrices can be written as , , and , respectively, where W is a (positive-definite) symmetric matrix, the eigenvalues of matrix A become identical; i.e., ₁, ..., _K = , as given by

(16)

Then, one can easily show that when > 0 (the SM_w model), f(S) follows a negative gamma distribution

(17)

where the parameter h = 1/(4N_e). Although it is rare, when < 0 (the SM_s model), f(S) follows a (positive) gamma distribution. In both situations, the gamma shape parameter of S is the (half) number of molecular phenotypes (K/2) (GU 2007).

	ESTIMATION AND EXAMPLES

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
The newly developed SM model is useful in practice only if the underlying parameters can be estimated. But this is difficult to accomplish using conventional statistical approaches, since the number of unknown parameters may be huge. Instead, we attempt to focus on two important parameters, K, critical for understanding gene pleiotropy, and , the measure for overall sequence conservation. Our purpose here is to estimate K and , without specifying other more detailed parameters.

Moments of evolutionary rate:
Under the SM model, the mean evolutionary rate of a protein sequence () is given by Equation 9. Therefore, any kth moment of the evolutionary rate is given by

(18)

To apply Equation 18 to estimate unknown parameters, we use the approximation S/(1 – e^–S) e^–|S|(1 + c|S|), where c 0.5772 is Euler's constant. As shown by Figure 3, this approximation is satisfactory for S 0. Thus, under the SM_w model (stabilizing selection with weak microadaptation), which assures S(y) = –2N_ey'U^–1y < 0, the kth moment of the evolutionary rate can be approximated as follows:

(19)

View larger version (15K): [in this window] [in a new window] [Download PPT slide]   FIGURE 3.— (A) The /v – S (the rate ratio of evolution to mutation vs. the selection intensity) relationship, for the accurate formula and the approximation. (B) The gK K plotting for estimating the effective number of molecular phenotypes.

(20)

and the second moment of (k = 2) by

(21)

Effective number (K_e) of molecular phenotypes:
In general, the mean evolutionary rate depends on K distinct eigenvalue parameters (N_ei's), while the second moment depends on N_ei's and , and so forth. To develop a simple method for estimating K on the basis of Equations 20 and 21, we invoke the assumption of N_ei >> 1, leading to

(22)

and

(23)

respectively. Notably, the first and second moments of evolutionary rate depend only on two parameters, K and the product . Further, we found that the ratio of the second moment to the first moment (), denoted by , is given by

(24)

which is only K dependent. As shown in Figure 3, g_K decreases when K increases; g_K = 1 when K = 0, and when .

Equations 22–24 provide the theoretical foundation for estimating K when N_ei >> 1. In this case, the parameter K is interpreted as the number of molecular phenotypes that have experienced strong stabilizing selection. Thereafter we refer to K_e as the effective number of molecular phenotypes, which is less than the true number of molecular phenotypes. If the ratio g_K is known, K_e can be estimated according to Equation 24.

A widely used measure for is d_N/d_S, the ratio of nonsynonymous to synonymous substitutions, under the assumption that synonymous substitutions are almost neutral. We have realized that the second moment of rate is related to the H-measure (GU et al. 1995) for the rate variation among sites (UZZEL and CORBIN 1971), defined by . Ranging from 0 to 1, a higher value of H means a greater rate variation among sites and vice versa. One can easily verify

(25)

Therefore, the effective number of molecular phenotypes (K_e) is the solution of the following equation:

(26)

Effective selection intensity (S_e):
Equation 13 indicates that the mean selection intensity can be written as , where is the (arithmetic) average of selection intensity over molecular phenotypes. However, without knowing each N_ei, it is difficult to estimate . After realizing that the product is actually estimable, as implied by Equation 22, we define the effective selection intensity as , where is the geometric average of selection intensity per molecular phenotype; that is,

(27)

In general, , and when all _i's are roughly the same. At any rate, both geometric and arithmetric means of selection intensity have virtually the same effect on the rate of protein evolution. Replacing by d_N/d_S in Equation 22 and K by K_e, the effective selection intensity can be estimated by

(28)

Estimation procedure:
We propose an analytical pipeline to estimate the effective number of molecular phenotypes (K_e) and the effective selection intensity (S_e). Although it needs to be refined both statistically and computationally, our procedure consists of the following steps:

1. Infer the phylogenetic tree from a multiple alignment of homologous protein sequences. Although there is no methodological preference, we require that the inferred tree topology should be roughly the same over methods.
   
2. Estimate the nonsynonymous to synonymous ratio (d_N/d_S) from closely related coding sequences that satisfy d_N/d_S < 1.
   
3. Estimate the H-measure for rate variation among sites: Use the method of GU and ZHANG (1997) to infer the (corrected) number of changes at each site, under the inferred phylogeny.
   

Let and V(x) be the mean and variance of number of changes over sites, respectively. Assuming a Poisson process at each site, we obtain the mean evolutionary , where T is the total evolutionary time along the tree. Similarly, the variance of evolutionary rate among sites is given by . Then, H can be estimated by

(29)

And finally, one can estimate K_e and S_e according to Equations 26 and 28, respectively.

Examples:
For case studies, we compiled 20 vertebrate protein families. Each family includes eight complete sequences from human, mouse, dog, cow, chicken, Xenopus, fugu, and zebrafish, respectively, which were downloaded from Ensembl EnsMar (http://www.ensembl.org/Multi/martview) (KASPRZYK et al. 2004). The synonymous distance (d_S) and the nonsynonymous distance (d_N) between mammalian orthologs were estimated by codeml of the PAML package (YANG 1997). We used the same phylogeny (Figure 4) to estimate the H-index of rate variation among sites; using an alternative tree gave a very similar result (not shown).

View larger version (9K): [in this window] [in a new window] [Download PPT slide]   FIGURE 4.— Phylogenetic tree of vertebrates used in our data analysis.

	DISCUSSION

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

Our SM model provides an approach to understanding the role of gene pleiotropy in a gene network. Although a detailed analysis will be published elsewhere, we test the hypothesis that gene pleiotropy underlies genomic correlations related to the evolutionary rate of protein sequence (HE and ZHANG 2006; PAL et al. 2006; WOLF et al. 2006). That is, (i) evolutionary rate is inversely related to gene pleiotropy (K), (ii) gene pleiotropy determines gene dispensability, and (iii) gene pleiotropy underlies "hub" genes. In short, estimation of gene pleiotropy (K_e) provides a model-based data-analysis approach that avoids the interpretative difficulty in principle-component or partial-correlation analysis.

In the future, we shall refine the analytical pipeline for estimating K_e. We shall examine how statistical properties of d_N/d_S and H would affect the estimation efficiency of K_e. In addition to the delta method, another method to estimate the sampling variance of the estimated K_e is nonparametric bootstrapping. Under the W-model when the selection intensity S follows a (negative) gamma distribution, a likelihood approach can be developed. Since the method of NIELSEN and YANG (2003) is applicable only for closely related sequences, we shall develop a protein sequence-based likelihood approach.

Gene pleiotropy and allelic pleiotropy:
It should be noted that the effective number of molecular phenotypes (K_e) measures the functional pleiotropy of a gene, rather than a single mutation. For instance, K_e = 8.23 for trans-aldolase is the estimated number of distinct fitness dimensions implied by all random mutations in the coding sequence. Although the estimated K_e indicates that most genes may be highly pleiotropic, what it exactly means is very much gene specific.

The number of fitness dimensions affected by a single mutation () could be much less than K_e, say, 2–4 (DUDLEY et al. 2005). WAXMAN and PECK (1998) argued that a single optimal gene sequence may become common when three or more characters (as molecular phenotypes here) are affected by each mutation, (i.e., 3). The problem of estimating effectively from sequence data remains a question that requires further study.

Several studies (NIELSEN and YANG 2003; PIGANEAU and EYRE-WALKER 2003; EYRE-WALKER et al. 2006) estimated the shape parameter of gamma distribution for S. NIELSEN and YANG (2003) obtained the shape parameter of 3.22 from the primate mitochondrial genome sequences. Under the SM framework, it can be interpreted as K_e/2; see the above section for the W-model. That is, K_e = 2 x 3.22 = 6.44 for primate mitochondria proteins, consistent with our results (Table 1). In contrast, EYRE-WALKER et al. (2006) obtained a much smaller value (0.23) of the shape parameter from the human SNP data. We speculate that the estimate from population genetics data perhaps reflects the allelic pleiotropy (). Hence, EYRE-WALKER et al. (2006) may provide a rough estimate for = 2 x 0.23 = 0.46.

Assumptions and alternative models:
The SM model involves several simplifying assumptions, since an exhaustive consideration of all possible mutational effects on molecular phenotypes is intractable. Many are shown with previous phenotype–genotype models. One may refer to MARTIN and LENORMAND (2006) for a recent summary about rationales and criticisms. The first assumption is the single fitness optimum. Under this condition we assume that disruptive selection with alternative optima following an adaptive process (reviewed in KINGSOLVER et al. 2001 and ELENA and LENSKI 2003) occurs rarely at the molecular level. Moreover, we used the Gaussian fitness function (see Equation 1) to measure the fitness consequences of a deviation from the optimum. LANDE (1980) showed that when the population mean is close to the optimum, a Gaussian-like function can be a local kernel approximation for many arbitrary fitness functions.

Second, the effect of mutations on molecular phenotypes is assumed to be continuous and symmetric. The continuous property is consistent with some empirical evidence from random mutations of genes (KEIGHTLEY 1994; KEIGHTLEY and CABALLERO 1997; IMHOF and SCHLOTTERER 2001). The symmetric assumption implies that even if mutations bias toward higher or lower values of molecular phenotypes, the trend should be small compared to the mutational variance. In practice we used a multiple Gaussian distribution for mutational effects. In the non-Gaussian case, it is possible to reduce kurtosis after an appropriate transformation such that mutational effects on a "transformed" molecular phenotype are roughly Gaussian-like. In short, we conclude that in the absence of practically useful alternatives, single optimum, Gaussian fitness function, and continuous symmetric effect of mutations remain useful working assumptions. Further work will examine the effects of these assumptions using computer simulations.

A more important issue is whether the estimation procedure for the effective number of molecular phenotypes (K_e) is robust under various alternative models. The first alternative model is Fisher's geometric model, in which adaptive change is represented by stepwise movement of a point toward the center of the hypersphere. In this case, the number of molecular phenotypes (K) is Fisher's geometric dimension. Several studies (HARTL and TAUBES 1996, 1998; SELLA and HIRSH 2005) showed that at equilibrium, Fisher's geometric model can be approximated by allowing the current population to vary randomly around the fixed optimum. A second alternative model is to define the fitness to be a function of flux through a metabolic pathway (KACSER and BURNS 1979; HARTL and TAUBES 1996), rather than a Gaussian-like function. Preliminary results have indicated that in both cases, estimation of K_e varies only slightly. It seems that K_e is a robust measure of gene pleiotropy, which may not be very sensitive to the specific evolutionary model. We shall address this issue more extensively in future studies.

	APPENDIX

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
Mean and variances of S—derivation of Equations 13 and 14:
For the distribution of y that follows a multi(K)-variate normal distribution,

(A1)

we note that for any quadratic function y'U^–1y, the first and second moments can be expressed as

(A2)

respectively, where the operator tr[.] means the trace of a square matrix. From Equation 7 it is straightforward that and . Denote the matrix A = U^–1_m. After noting that the trace of the matrix is the sum of the eigenvalues, i.e., and , we obtain Equations 13 and 14.

Canonical representation of S—derivation of Equation 15:
The canonical representation will choose molecular phenotypes such that _w = diag(_,1, ..., ). Recall that A = U^–1_m = Z – , where and . Then, one can verify , since the trace of a matrix is the sum of its diagonal elements. Next, matrix can be written as . It follows that , where is the ith eigenvalue of matrix _µm. From the relation tr[A] = tr[Z] – tr[], together we have

(A3)

leading to Equation 15, where /[].

The canonical form for the variance of S can be obtained similarly, but the algebra is somewhat tedious, largely due to the detail of A² = Z² + ² – 2Z. Let . One thus can show the trace and . Noting that Z = diag(_,1, ..., _,K)_mµm, we obtain , where is the ith eigenvalue of matrix _mµm. Therefore, from the fact tr[A²] = tr[Z²] + tr[²] – 2 tr[Z], we have

(A4)

Condition of S-gamma distribution:
Given p(y) characterized by the covariance _w, the distribution of selection intensity, f(S), can be uniquely determined by the probability theory that for any S* and y* that satisfy S* = –2N_e(y*'y*) in Equation 7, the relationship holds. However, the analytical form of f(S) is not available, except for the case when the matrix U in the quadratic form and the covariance in p(y), _m, satisfies A = _mU^–1 = I. In other words, the matrix A is diagonal, having the same eigenvalues .

Moments of evolutionary rate—derivation of Equations 20 and 21:
Here we show several analytical integral results. First we consider , which is given by

(A5)

where . It follows that

(A6)

In the same manner, we have given by

(A7)

Next we consider . From the derivation of I₁, one can verify

(A8)

From Equation A2, we obtain

(A9)

In the same manner, we show that as follows:

(A10)

Finally, we consider the integral . Similar to the derivations of I₁–I₅, one can write

(A11)

where matrix B* = ( + 8N_eU^–1)^–1. Since U^–1B* = (U + 8N_eI)^–1 = (A^–1 + 8N_eI)^–1, and |B*/_m| = |I + 8N_eA|, I₅ can be written as follows:

(A12)

From Equation 19, one can show that and , both of which are determined a single matrix A. Hence, we have derived Equations 20 and 21 on the basis of the matrix theory that, for a square matrix A, the trace is the sum of the eigenvalues, the determinant is the production of the eigenvalues, and any eigenvalue of the matrix that is a function of A is the function of the corresponding eigenvalue of A.

	ACKNOWLEDGEMENTS

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 
The author is grateful to Adam Eyre-Walker and the reviewing editor Antony Long for constructive comments on the manuscript. Assistance in data analysis from Dongping Xu and Zhixi Su is appreciated.

	LITERATURE CITED

TOP ABSTRACT THE MODEL ANALYTICAL RESULTS ESTIMATION AND EXAMPLES DISCUSSION APPENDIX ACKNOWLEDGEMENTS LITERATURE CITED

 

BARTON, N., 1990 Pleiotropic models of quantitative variation. Genetics 124: 773–782.[Abstract]

BATAILLON, T., 2000 Estimation of spontaneous genome-wide mutation rate parameters: Whither beneficial mutations? Heredity 84: 497–501.[CrossRef][Medline]

DRUMMOND, D. A., J. D. BLOOM, C. ADAMI, C. O. WILKE and F. H. ARNOLD, 2005 Why highly expressed proteins evolve slowly. Proc. Natl. Acad. Sci. USA 102: 14338–14343.[Abstract/Free Full Text]

DUDLEY, A., D. JANSE, A. TANAY, R. SHAMIR and G. CHURCH, 2005 A global view of pleiotropy and phenotypically derived gene function in yeast. Mol. Syst. Biol. 1: 2005.0001.

ELENA, S. F., and R. E. LENSKI, 2003 Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat. Rev. Genet. 4: 457–469.[Medline]

EYRE-WALKER, A., M. WOOLFIT and T. PHELPS, 2006 The distribution of fitness effects of new deleterious amino acid mutations in humans. Genetics 173: 891–900.[Abstract/Free Full Text]

FISHER, R. A., 1930 The genetical theory of natural selection. Oxford University Press, Oxford.

GILLESPIE, J. H., 1991 The Causes of Molecular Evolution. Oxford University Press, Oxford.

GU, X., and J. ZHANG, 1997 A simple method for estimating the parameter of substitution rate variation among sites. Mol. Biol. Evol. 14: 1106–1113.[Abstract]

GU, X., 2007 Stabilizing selection of protein function and distribution of selection coefficient among sites. Genetica (in press).

GU, X., Y. X. FU and W. H. LI, 1995 Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites. Mol. Biol. Evol. 12: 546–557.[Abstract]

HARTL, D. L., and C. H. TAUBES, 1996 Compensatory nearly neutral mutations: selection without adaptation. J. Theor. Biol. 182: 303–309.[CrossRef][Medline]

HARTL, D. L., and C. H. TAUBES, 1998 Towards a theory of evolutionary adaptation. Genetica 102/103: 525–533.

HE, X., and J. ZHANG, 2006 Toward a molecular understanding of pleiotropy. Genetics 173: 1885–1891.[Abstract/Free Full Text]

IMHOF, M, and C. SCHLOTTERER, 2001 Fitness effects of advantageous mutations in evolving Escherichia coli populations. Proc. Natl. Acad. Sci. USA 98: 1113–1117.[Abstract/Free Full Text]

KACSER, H., and J. A. BURNS, 1979 Molecular democracy: Who shares the controls? Biochem. Soc. Trans. 7: 1149–1161.[Medline]

KASPRZYK, A., D. KEEFE, D. SMEDLEY, D. LONDON, W. SPOONER et al., 2004 EnsMart: a generic system for fast and flexible access to biological data. Genome Res. 14: 160–169.[Abstract/Free Full Text]

KEIGHTLEY, P. D., 1994 The distribution of mutation effects of viability in Drosophila melanogaster. Genetics 138: 1315–1322.[Abstract]

KEIGHTLEY, P. D., and A. CABALLERO, 1997 Genomic mutation rates for lifetime reproductive output and lifespan in Caenorhabditis elegans. Proc. Natl. Acad. Sci. USA 94: 3823–3827.[Abstract/Free Full Text]

KIMURA, M., 1979 Model of effectively neutral mutations in which selective constraint is incorporated. Proc. Natl. Acad. Sci. USA 75: 1934–1937.[CrossRef]

KIMURA, M., 1983 The Neutral Theory of Molecular Evolution. Cambridge University Press, Cambridge, UK.

KINGSOLVER, J. G., H. E. HOEKSTRA, J. M. HOEKSTRA, D. BERRIGAN, S. N. VIGNIERI et al., 2001 The strength of phenotypic selection in natural populations. Am. Nat. 157: 245–261.[CrossRef]

LANDE, R., 1980 The genetic covariance between characters maintained by pleiotropic mutations. Genetics 94: 203–215.[Abstract/Free Full Text]

LYNCH, M., J. BLANCHARD, D. HOULE, T. KIBOTA, S. SCHULTZ et al., 1999 Perspective: spontaneous deleterious mutation. Evolution 53: 645–663.[CrossRef]

MACLEAN, R. C., G. BELL and P. B. RAINEY, 2004 The evolution of a pleiotropic fitness tradeoff in Pseudomonas fluorescens. Proc. Natl. Acad. Sci. USA 101: 8072–8077.[Abstract/Free Full Text]

MARTIN, G., and T. LENORMAND, 2006 A general multivariate extension of Fisher's geometrical model and the distribution of mutation fitness effects across species. Evolution 60: 893–907.[CrossRef][Medline]

NIELSEN, R., and Z. YANG, 2003 Estimating the distribution of selection coefficients from phylogenetic data with applications to mitochondrial and viral DNA. Mol. Biol. Evol. 20: 1231–1239.[Abstract/Free Full Text]

OHTA, T., 1973 Slightly deleterious mutant substitutions in evolution. Nature 246: 96–98.[CrossRef][Medline]

OHTA, T., 1977 Extension to the neutral mutation random drift hypothesis, pp. 148–176 in Molecular Evolution and Polymorphism, edited by M. KIMURA. National Institute of Genetics, Mishima, Japan.

OTTO, S. P., 2004 Two steps forward, one step back: the pleiotropic effects of favoured alleles. Proc. Biol. Sci. 271: 705–714.[CrossRef][Medline]

PAL, C., B. PAPP and M. J. LERCHER, 2006 An integrated view of protein evolution. Nat. Rev. Genet. 7: 337–348.[CrossRef][Medline]

PIGANEAU, G., and A. EYRE-WALKER, 2003 Estimating the distribution of fitness effects from DNA sequence data: implications for the molecular clock. Proc. Natl. Acad. Sci. USA 100: 10335–10340.[Abstract/Free Full Text]

POON, A., and S. P. OTTO, 2000 Compensating for our load of mutants: freezing the meltdown of small populations. Evolution 54: 1467–1479.[CrossRef][Medline]

SHAW, F. H., C. J. GEYER and R. G. SHAW, 2002 A comprehensive model of mutations affecting fitness and inferences for Arabidopsis thaliana. Evolution 56: 453–463.[CrossRef][Medline]

SELLA, G., and A. E. HIRSH, 2005 The application of statistical physics to evolutionary biology. Proc. Natl. Acad. Sci. USA 102: 9541–9546.[Abstract/Free Full Text]

TURELLI, M., 1985 Effects of pleiotropy on predictions concerning mutation-selection balance for polygenic traits. Genetics 111: 165–195.[Abstract/Free Full Text]

UZZEL, T., and K. W. CORBIN, 1971 Fitting discrete probability distribution to evolutionary events. Science 172: 1089–1096.[Abstract/Free Full Text]

WAGNER, A, 2000 The role of population size, pleiotropy and fitness effects of mutations in the evolution of overlapping gene functions. Genetics 154: 1389–1401.[Abstract/Free Full Text]

WAGNER, G. P., 1989 Multivariate mutation-selection balance with constrained pleiotropic effects. Genetics 122: 223–234.[Abstract/Free Full Text]

WALL, D. P., A. E. HIRSH, H. B. FRASER, J. KUMM, G. GIAEVER et al., 2005 Functional genomic analysis of the rates of protein evolution. Proc. Natl. Acad. Sci. USA 102: 5483–5488.[Abstract/Free Full Text]

WAXMAN, D., and J. R. PECK, 1998 Pleiotropy and the preservation of perfection. Science 279: 1210–1213.[CrossRef][Medline]

WELCH, J. J., and D. WAXMAN, 2003 Modularity and the cost of complexity. Evolution 57: 1723–1734.[CrossRef][Medline]

WOLF, Y. I., L. CARMEL and E. V. KOONON, 2006 Unifying measures of gene function and evolution. Proc. Biol. Sci. 273: 1507–1515.[CrossRef][Medline]

WRIGHT, S., 1968 Evolution and the Genetics of Populations. University of Chicago Press, Chicago.

YANG, Z., 1997 PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. 13: 555–556.[Free Full Text]

ZHANG, X. S., and W. G. HILL, 2003 Multivariate stabilizing selection and pleiotropy in the maintenance of quantitative genetic variation. Evolution 57: 1761–1775.[CrossRef][Medline]