## Abstract

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.

ALTHOUGH the wide availability of high-throughput data has greatly facilitated the study of evolutionary genomics, the foundation of molecular evolution has been relatively “static” (Kimura 1983). Basically, the evolutionary fate of a mutant is largely determined by the prespecified fitness value (the coefficient of selection) and the effective population size. Although this classical treatment has successfully removed the complexity of genotype–phenotype association, problems may appear when the association itself has been one of the research highlights. The recent debate about genomic correlations among evolutionary rate, interaction, expression, and dispensability has indicated the limitation of current evolutionary models (reviewed in Pal *et al*. 2006). Alternatively, researchers have to use the partial-correlation or principle-component methods to analyze major determinants governing protein evolution (*e.g*., Drummond *et al*. 2005; Wall *et al*. 2005; Wolf *et al*. 2006), because these approaches are not model dependent.

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

#### Molecular phenotypes and gene pleiotropy:

Pleiotropy of a gene (protein) function can be modeled by multiple (*K*) molecular phenotypes, denoted by (*y*_{1}, …, *y _{K}*). Each molecular phenotype

*y*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.

_{i}Following the widely used procedure (*e.g*., Lande 1980; Turelli 1985; Waxman and Peck 1998), we assume that the fitness function of molecular phenotypes **y** = (*y*_{1}, …, *y _{K}*)′ is Gaussian-like;

*i.e*.,(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

*i*th diagonal element measures the strength of stabilizing selection on the

*i*th molecular phenotype, while the

*ij*th nondiagonal element σ

_{w,ij}measures the correlated stabilizing selection on

*y*and

_{i}*y*.

_{j}Let **y**_{0} 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 _{0}** = μ, leading to purifying selection; that is, ρ(

**y**|

**y**= μ) ≤ 0 always holds.

_{0}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

_{0}*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.

#### Stabilizing selection model:

Under the classical stabilizing model, the population mean of molecular phenotypes (**y**_{0}) is always fixed at the optimum (μ). Without loss of generality, we assume that **y**_{0} = μ = **0**. From Equation 2 one can show that , reflecting the stabilizing (purifying) selection against deleterious mutations, which necessarily results in a deviation from the optimum. A consequence of this model is that sequence evolution is dominated by the fixation of very slightly deleterious mutations (Ohta 1973; Kimura 1979, 1983).

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} = **0**. In the following we consider mainly the selection intensity *S*(**y**) = 4*N*_{e}ρ_{0}(**y**) = −2*N*_{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 _{0}** =

**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**) = 4*N*_{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

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**, α_{1}, …, α_{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 −2*N*_{e}α_{i} on the mean selection intensity () of the gene.

##### Canonical representation of S̄:

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 *i*th (canonical) molecular phenotype. Let be the *i*th diagonal element of matrix **Σ**_{w} or the mutational variance for the *i*th (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 *i*th (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:

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.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.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*., α_{1}, …, α_{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/(4*N*_{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

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 *k*th 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**) = −2*N*_{e}**y′U**^{−1}**y** < 0, the *k*th moment of the evolutionary rate λ can be approximated as follows:(19)

As shown in the appendix, the mean evolutionary rate (*k* = 1) is given by(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*_{e}α_{i}'s), while the second moment depends on *N*_{e}α_{i}'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*_{e}α_{i} ≫ 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*= 1 when

_{K}*K*= 0, and when .

Equations 22–24 provide the theoretical foundation for estimating *K* when *N*_{e}α_{i} ≫ 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*_{e}α_{i}, 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:

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.

Estimate the nonsynonymous to synonymous ratio (

*d*_{N}/*d*_{S}) from closely related coding sequences that satisfy*d*_{N}/*d*_{S}< 1.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).

Table 1 shows the results. The *d*_{N}/*d*_{S} ratio varies substantially among genes (from 0.001 to 0.274), and the *H*-index ranges from 0.286 to 0.886. Consequently, the estimated effective number of molecular phenotypes *K*_{e} ranges from 2.3 to 20.4, with an average of 8.77 (the median is 8.23). Our analysis suggests that most genes are highly pleiotropic. That is, there are on average eight distinct components in the fitness that may be related to a gene. Roughly, ∼73% of the among-gene variation in the *d*_{N}/*d*_{S} ratio can be explained by the variation in *K*_{e} among genes. In other words, the level of gene pleiotropy associated with a gene may be the dominant factor for predicting sequence conservation, as predicted by Fisher (1930). The effective mean selection intensity (*S*_{e}) ranges from 5.39 to 23.63, with a median of 12.6. Note that *d*_{N}/*d*_{S} ratios in Table 1 were estimated from the human and mouse orthologous genes; using other mammals or taking an average would lead to ∼5% difference in the *K*_{e} estimation.

## DISCUSSION

#### Estimation of gene pleiotropy:

Although the theory of gene pleiotropy, the capacity of a gene to affect multiple phenotypic characters, has been used to explain many biological phenomena (Fisher 1930; Wright 1968; Barton 1990; Waxman and Peck 1998; Otto 2004; MacLean *et al*. 2004; Dudley *et al*. 2005), the characteristic level of pleiotropy remains largely unknown. Many phenotype–genotype models (*e.g*., Fisher 1930; Wright 1968; Hartl and Taubes 1996, 1998; Waxman and Peck 1998; Poon and Otto 2000) have linked gene pleiotropy to the dimensionality of the fitness function of the organism. The major contribution of this article is to formulate a statistical framework for estimating the capacity of a gene to significantly affect distinct components in the fitness, the effective number of molecular phenotypes (*K*_{e}). Thus, for the first time we have an empirical measure of gene pleiotropy. Examples in Table 1 show that *K*_{e} is typically ∼6–9, the number of distinct fitness components a gene may affect, which has been further confirmed by a large-scale data analysis (not shown).

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 × 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 × 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

#### 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**^{−1}**y**, 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 *i*th 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**^{2} = **Z**^{2} + **Ω**^{2} − 2**ZΩ**. Let . One thus can show the trace and . Noting that **ZΩ** = diag(_{,1}, …, _{,K})**Σ**_{m}**Σ**_{μ}**Σ**_{m}, we obtain , where is the *i*th eigenvalue of matrix **Σ**_{m}**Σ**_{μ}**Σ**_{m}. Therefore, from the fact tr[**A**^{2}] = tr[**Z**^{2}] + tr[**Ω**^{2}] − 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** = −2*N*_{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** = **Σ**_{m}**U**^{−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*_{1}, 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*_{1}–*I*_{5}, one can write(A11)where matrix **B*** = ( + 8*N*_{e}**U**^{−1})^{−1}. Since **U**^{−1}**B*** = (**U** + 8*N*_{e}**I**)^{−1} = (**A**^{−1} + 8*N*_{e}**I**)^{−1}, and |**B***/**Σ**_{m}| = |**I** + 8*N*_{e}**A**|, *I*_{5} 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**.

## Acknowledgments

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.

## Footnotes

Communicating editor: A. D. Long

- Received October 5, 2006.
- Accepted January 8, 2007.

- Copyright © 2007 by the Genetics Society of America