## Abstract

Although pleiotropy, the capability of a gene to affect multiple phenotypes, has been well known as one of the common gene properties, a quantitative estimation remains a great challenge, simply because of the phenotype complexity. Not surprisingly, it is hard for general readers to understand how, without counting phenotypes, gene pleiotropy can be effectively estimated from the genetics data. In this article we extensively discuss the Gu-2007 method that estimated pleiotropy from the protein sequence analysis. We show that this method is actually to estimate the rank (*K*) of genotype–phenotype mapping that can be concisely written as *K* = min(*r*, *P*_{min}), where *P*_{min} is the minimum pleiotropy among all legitimate measures including the fitness components, and *r* is the rank of mutational effects of an amino acid site. Together, the *effective gene pleiotropy* (*K*_{e}) estimated by the Gu-2007 method has the following meanings: (i) *K*_{e} is an estimate of *K =* min(*r*, *P*_{min}), the rank of a genotype–phenotype map; (ii) *K*_{e} is an estimate for the minimum pleiotropy *P*_{min} only if *P*_{min} < *r*; (iii) the Gu-2007 method attempted to estimate the pleiotropy of amino acid sites, a conserved proxy to the true gene pleiotropy; (iv) with a sufficiently large phylogeny such that the rank of mutational effects at an amino acid site is *r* → 19, one can estimate *P*_{min} between 1 and 19; and (v) *K*_{e} is a conserved estimate of *K* because those slightly affected components in fitness have been effectively removed by the estimation procedure. In addition, we conclude that mutational pleiotropy (number of traits affected by a single mutation) cannot be estimated without knowing the phenotypes.

PLEIOTROPY, or the capability of a gene to affect multiple phenotypes (Fisher 1930; Wright 1968), has played a central role in genetics, development, and evolution (*e.g.*, Lande 1980; Turelli 1985; Wagner 1989; Barton 1990; Keightley 1994; Hartl and Taubes 1996, 1998; Waxman and Peck 1998; Lynch *et al.* 1999; Bataillon 2000; Orr 2000; Poon and Otto 2000; Wagner 2000; Elena and Lenski 2003; Welch and Waxman 2003; Wingreen *et al.* 2003; Zhang and Hill 2003; MacLean *et al.* 2004; Otto 2004; Eyre-Walker *et al.* 2006; Pigliucci 2008; Wagner *et al.* 2008). Although functional genomics has brought high-throughput data to bear on the nature and extent of pleiotropy (Dudley *et al.* 2005; Ohya *et al.* 2005; Pal *et al.* 2006; Cooper *et al.* 2007), this issue remains highly controversial, largely because of the phenotypic complexity; see Wagner and Zhang (2011), Hill and Zhang (2012), and Paaby and Rockman (2013) for recent reviews and comments.

Meanwhile, a new research line has emerged in the past decade, aiming at the estimation of gene pleiotropy from genetics or sequence data, rather than the measurement of affected phenotypes (Martin and Lenormand 2006; Gu 2007a,b; Chevin *et al.* 2010; Su *et al.* 2010; Zeng and Gu 2010; Razeto-Barry *et al.* 2012; Chen *et al.* 2013). For instance, Gu (2007a) developed a statistical method to estimate the “effective pleiotropy” (*K*_{e}) of a gene from the multiple-sequence alignment (MSA) of protein sequences. Most genes have *K*_{e} in the range between 1 and 20 (Su *et al.* 2010), with the medium *K*_{e} = 6.5 of these estimates that is comparable to some empirical pleiotropy measures (Wagner and Zhang 2011; Paaby and Rockman 2013). However, we have to acknowledge that this approach is not very easy to be understood; that is, How can we know the amount of multifunctionality (pleiotropy) of a gene without biologically knowing each functional component?

In this article we try to answer the following questions: What quantity was actually estimated by the method proposed by Gu (2007a)? And under which condition can this estimate be interpreted as effective gene pleiotropy? To answer these questions, we introduce the concept of minimum pleiotropy (*P*_{min}), which has an intrinsic relationship with the rank of the genotype–phenotype map. Therefore, we propose that the method of Gu (2007a) is to estimate the rank of the genotype–phenotype map that can be regarded as effective estimation of *P*_{min}. Simulations and empirical studies are carried out to test several predictions from the new theory.

## Materials and Methods

### Data sets

Random single-copy gene samples were from three sets of model organisms, respectively. MSAs of vertebrate genes were from the data set of Su *et al.* (2010), while those of five yeast homologous genes (*Saccharomyces cerevisiae*, *Candida glabrata*, *Kluyveromyces lactis*, *Debaryomyces hansenii*, and *Yarrowia lipolytica*) were from the Génolevures database (http://cbi.labri.fr/Genolevures), and those of 12 *Drosophila* genes were from the Consortium of *Drosophila* genome projects.

### Estimation of *K*_{e} from protein sequences

We calculated the *d*_{N}*/d*_{S} (the ratio of nonsynonymous to synonymous rates) of each gene, based on the orthologous sequences from two closely related species: human and mouse for the vertebrate data set, *S. cerevisiae* and *S. paradoxus* for the yeast data set, and *D*. *melanogaster* and *D. sechellia* for the *Drosophila* data set. After the gene phylogeny was reconstructed by the neighbor-joining method, we estimated *H*, a normalized index for the rate variation among sites (Gu 2007a). It follows that the effective gene pleiotropy (*K*_{e}) was estimated by solving the equation [*d*_{N}*/d*_{S}]/(1 − *H*) = 2* ^{−Kj}*[1 +

*ϕ*(

*K*)], where

_{j}*ϕ*(

*K*

_{e}) = 0.0208

*K*

_{e}(

*K*

_{e}+ 2)/(1 + 0.289

*K*

_{e}). The technical procedures for estimating

*K*

_{e}are described by Gu (2007a) and Su

*et al.*(2010) in detail.

### Estimation of *K*_{e} from nucleotide sequences

Vertebrate genes used in this analysis were from the data set of Su *et al.* (2010). For each gene, the MSA of nucleotides was transformed exactly from the MSA of amino acids. We first calculated the quantity *g* = (*d/d*_{S})/(1 − *H*) for each gene. Here *d* is the evolutionary distance at the first and second nucleotide positions (*d*_{12}) or at all three positions (*d*_{123}), and *d*_{S} is the synonymous distance. All these distances were estimated between human and mouse, by the Kimura two-parameter method. When the phylogeny is given, the normalization index of rate variation among nucleotide sites (*H*) can be calculated by the algorithm of Gu and Zhang (1997) with some modifications. In practice, when the method of Gu (2007a) is applied for the nucleotide sequences, we may face a technical problem: Due to the saturated synonymous substitutions at many nucleotide sites, estimation of the quantity *g* would be subject to a large sampling variance, resulting in a high number of not-applicable (NA) cases. To avoid this problem, we used the following procedure: First, calculate the mean of *g* over all genes, as well as the confidence interval (25% and 75% quantiles, respectively); and, second, use the mean of *g* to estimate the mean of *K*_{e}, as well as the 25% and 75% quantiles, respectively.

### Computer simulation

In the simulation procedure, we first generate the matrix **A** = **Σ _{w}**

^{−1}

**Σ**as follows: Given any fixed

_{m}*n*(

*n*= 6, 10, 14, and 18), we generated

*n*real random variables between 1 and 10 and assigned these variables randomly as diagonal elements of

**Σ**; hence it is of full rank. However, matrix

_{w}**Σ**is partitioned into two components: a full,

_{m}*r*-rank matrix

**Σ***that was generated similar to procedure (

_{m}*r*is less than

*n*) and then a diagonal matrix with

*n − r*small but the same positive values of diagonal elements such as

*B*

_{low}= 0.001. Together, matrix

**A**is given by

**A**= C

**Σ**

_{w}^{−1}

**Σ**, where the constant

_{m}*C*is chosen such that the mean of first

*r*eigenvalues of

**A**equals a specified value

*B*

_{0}. The simulation study was carried out in Matlab. Note that Matlab script (KeEstimation.m.txt) for the computer simulation and the software (GenePleio_win-v15.rar, codes included) for the estimation of

*K*

_{e}are available in supporting information, File S1. The real data sets used in this study are deposited in DryAd (http://datadryad.org/; doi: 10.5061/dryad.s84fm).

## Results and Discussion

### Many-face problem of pleiotropy under a genotype–phenotype map

Like many biological terminologies, the concept of pleiotropy is easy to understand but difficult to measure. Practically, biologists address this issue by experimentally identifying distinct phenotypes, *in vitro* or *in vivo*, caused by the same mutation (mutational pleiotropy) or mutations from the same gene (gene pleiotropy). Although these studies largely focused on the nature of phenotypes, the results did provide a “minimum” estimate about the extent of pleiotropy. Thus, one may compare whether a gene is more pleiotropic than others, as long as the same technological platform is used, *e.g.*, the mouse knockout experimentation.

However, different pleiotropy measures based on different technological platforms certainly represent biologically significant differences that may even lead to contradictory interpretations (Wagner and Zhang 2011; Paaby and Rockman 2013). To illustrate this problem, we invoke the basic system formulated by Paaby and Rockman (2013), which has three steps:

1. Molecular pleiotropy (

*f*) refers to the number of protein functions, most likely*in vitro*. For instance, an enzyme is called pleiotropic if it can catalyze many distinct substrates.2. Developmental pleiotropy (

*d*) refers to the number of developmental pathway autonomies that underlie diverse phenotypes, including diseases. Developmental pleiotropy provides concrete examples to demonstrate how a single gene or mutation may affect distinct phenotypes*in vivo*.3. Selectional pleiotropy (

*n*) refers to the number of molecular phenotypes (Gu 2007a), each of which corresponds to a single nontrivial fitness component. In essence, selectional pleiotropy is the dimension of Fisher’s geometric model, which underlies most statistical studies of the genotype–phenotype map (Lande 1980; Turelli 1985; Wagner 1989; Barton 1990).

While the difference between *f* and *d* illustrates the many-face problem of pleiotropy, the difference between *n* and *f* (or *d*) represents the conceptual gap between empirical pleiotropy and Fisher’s geometric model. Although the relationship between *f*, *d*, and *n* could be astonishingly complex, we assume that, from the genotype (*G*) to the phenotype, the map direction can be represented as *G* → *f* → *d* → *n*. Meanwhile, under Fisher’s model, the map direction is as simple as *G* → *n*. While one may view this one-step map in Fisher’s model as merely a theoretical abstract of many cryptic steps, we argue that the dimension of Fisher’s model could be misled by this oversimplification; it equals *n* only when *n < f* and *n < d*; otherwise, it is <*n*. Actually, it is the minimum pleiotropy *P*_{min} = min(*f*, *d*, *n*) that determines the dimension of Fisher’s model. In other words, *P*_{min} is the rank of phenotype complexity, providing a link between the real genotype–phenotype map and Fisher’s model.

The concept of *P*_{min} can be naturally extended to any general genotype–phenotype map with *N* steps. Without loss of generality, we suppose that the lowest level is the molecular functions (*f*) related to protein activities and the highest level is the molecular phenotypes (*n*) related to fitness components. At each *i*th step node, the pleiotropy level can be symbolically denoted by *P _{i}*,

*i*= 1, … ,

*N*, where

*P*

_{1}=

*f*and

*P*=

_{N}*n*. Thus, the minimum pleiotropy can be generally defined as follows: (1)While minimum pleiotropy may provide one solution to the many-face problem, the question that remains is under what condition we are able to estimate

*P*

_{min}. To this end, we should address the rank of the genotype–phenotype map first.

### Rank of the genotype–phenotype map

Without loss of generality, the genotype–phenotype map can be schematically represented as *G* → *f* → (*P*_{2}, … , *P _{N−}*

_{1}) →

*n*, where (

*P*

_{2}, … ,

*P*

_{N−}_{1}) represents the black box to indicate the phenotype complexity. We then define the

*rank of the genotype–phenotype map*(

*K*) as the canonical number of mapping paths from the genotype (

*G*) to the fitness components (

*n*) through the entire phenotypic space. In other words,

*K*is the degrees of freedom for a particular genotype–phenotype mapping, which highly depends on the nature of the genotype. In the case of single mutations treated as genotypes, the

*mutational pleiotropy*(a single mutation can affect multiple phenotypes) means that those corresponding mapping paths start from the same mutational event. Indeed, these affected phenotypes are expected to be highly correlated, without many degrees of freedom, simply because the same biochemical–structural property of the mutation underlies most phenotypic consequences. By contrast, a pleiotropic gene may contain many mutations, each of which affects only one specific phenotype. As a result, a gene may contain a bunch of related but distinct mapping paths that have high degrees of freedom. This issue can be well addressed by the

*rank of mutational effects*(

*r*), defined by the canonical number of mutations within a genotype unit that cause distinct phenotypes. Intuitively,

*r*is the number of starting points of map paths; obviously,

*r*= 1 for single mutations.

Putting this together, the rank (*K*) of the genotype–phenotype map is the smaller one of the rank of phenotypic complexity (minimum pleiotropy) and the rank of mutational effects; that is, (2)As is discussed later, this geometric structure of the genotype–phenotype map can be rigorously derived under the linear transformation model of mapping. Obviously, to what extent *K* can be interpreted as the minimum pleiotropy (*P*_{min}) depends on the unknown rank (*r*) of mutational effects; *K* = *P*_{min} if *P*_{min} < *r*. In particular, when *r* = 1, *K* always = 1, regardless of how large *P*_{min} is. Nevertheless, the up limit of *r* for given genotype data is controlled by the number of generic states (*M*). For instance, *M =* 2 for random nucleotides, *M =* 4 for nucleotide sites, and *M =* 20 for amino acid sites. In fact, the maximum rank (*r*_{max}) for a given genotype is the number of generic states minus one (the wild type); *i.e.*, *r*_{max} = *M* − 1. Moreover, the principle of “microergodicity” claims that the mutational process can reach over all generic states for a sufficiently long time; in this case we have *r* → *r*_{max}. In practice, the “real” rank of mutational effects (*r*) may be <<*r*_{max} if the underlying sampling process has no sufficient time for mutations to reach all states.

### Rank of the genotype–phenotype map estimated by the Gu (2007a) method

In this section, we show that the statistical method developed by Gu (2007a) (Gu-2007 method) actually was to estimate the rank (*K*) of the genotype–phenotype map from the protein sequences. Since the genotype–phenotype map can convey the geometric information, the rank *K*, between the stabilizing selection and the protein sequence evolution, we are able to, in retrospect, estimate *K* based on the analysis of protein sequence data without the help of observed phenotypes. Since *K* = min(*r*, *P*_{min}), we thus conclude that the minimum pleiotropy (*P*_{min}) can be estimated without counting the number of phenotypes, if it is in the range of [1, *r*].

#### Derivation of *K* = min(*r*, *P*_{min}) under a linear mapping model:

Consider the case where the genotype–phenotype map is schematically represented as a linear mapping model *G = r* → *P*_{1} → , … , → *P _{N} = n*, which can be characterized by three matrices (Wagner 1989): an

*r × r*variance–covariance matrix (

**V**) for the correlated mutational effects of genotypes (

*G*), an

*r × n*transformation matrix

**T**that links between

*r*number of mutational effects and

*n*number of fitness components (molecular phenotypes), and an

*n × n*stabilizing selection matrix

**Σ**that is full rank of

_{w}*n*. The linear mapping model claims that the mutation matrix

**V**and transformation matrix

**T**together determine how mutations ultimately affect

*n*molecular phenotypes, as characterized by an (

*n × n*) matrix defined by

**Σ**. However, the rank of

_{m}= T′VT**Σ**is not

_{m}*n*except for some special cases, because

**T**is a multiplication of those consecutive step-transformation matrices; that is,

**T**

*=*

*T*

_{1}

*T*

_{2},

*…*,

*T*

_{N}, where

*T*

_{1}is the matrix with

*r*rows and

*P*

_{1}columns,

*T*

_{2}is the one with

*P*

_{1}rows and

*P*

_{2}columns, and so forth. Thus, the linear algebra theory states that the rank of

**Σ**is the minimum of (

_{m}*r*,

*P*

_{1},

*…*,

*P*

_{N−}_{1},

*P*). Noting that the minimum pleiotropy

_{N}*P*

_{min}= min(

*P*

_{1},

*…*,

*P*

_{N−}_{1},

*P*), we have the rank of

_{N}**Σ**given by min(

_{m}*r*,

*P*

_{min}).

Under this transformation model, the rank (*K*) of the genotype–phenotype map turns out to be the rank of matrix **A = Σ _{w}^{−1}Σ_{m}**. As

**Σ**is a full rank of

_{w}*n*, the rank of

**A**is the same as that of

**Σ**. Thus, we have shown

_{m}*K*= min(

*r*,

*P*

_{min}) ≤ min (

*r*,

*n*). The equals sign holds when

*n = P*is the smallest among

_{N}*P*

_{1}, … ,

*P*or the one-step (

_{N}*N*= 1) genotype–phenotype map (Chen

*et al.*2013).

#### Estimation of *K* = min(*r*, *P*_{min}) by the Gu-2007 method:

Gu (2007a) formulated the statistical model of protein sequence under *n*-dimensional stabilizing (purifying) selection. The results showed that, in addition to mutation rate *v*, the first (mean) and the second moments of the evolutionary rate are dominated by the *n* eigenvalues (*α _{i}*,

*i*= 1, … ,

*n*) of matrix

**A**=

**Σ**, which are nonnegative,

_{w}^{−1}Σ_{m}*i.e.*, either positive or zero. Without loss of generality, one may set

*α*

_{1}≥ α

_{2}… ≥

*α*≥ 0. For instance, the mean evolutionary rate can be written as (3)where

_{n}*c*= 0.5772,

*B*= 4

_{i}*N*

_{e}

*α*is the

_{i}*i*th component of selection intensity, and

*N*

_{e}is the effective population size. A notable property of Equation 3 is that only those nonzero

*B*’s affect the evolutionary rate. Since some eigenvalues (

_{i}*α*’s) of matrix

_{i}**A**may be zero, the number of eigenvalues (

*n*) in Equation 3 should be replaced by that of positive eigenvalues. We formulate this argument as follows. Since the rank of matrix

**A**is

*K*= min(

*r*,

*P*

_{min}) ≤ min(

*r*,

*n*) ≤

*n*, matrix

**A**has the first

*K*nonzero positive eigenvalues; that is,

*α*> 0 for any

_{i}*i*≤

*K*, and

*α*= 0 for all

_{i}*i*>

*K*. It immediately follows that

*B*

_{1}≥

*B*

_{2}≥, … , ≥

*B*> 0, whereas

_{K}*B*

_{K+}_{1}=, … , =

*B*= 0, suggesting that Equation 3 is equivalent to (4)As similar results can be derived for any

_{n}*k*th moment of evolutionary rate (not shown), we conclude that the Gu-2007 method is to estimate the rank of the genotype–phenotype map

*K*= min(

*r*,

*P*

_{min}), rather than the gene pleiotropy directly.

### From the rank of the genotype–phenotype map to effective gene pleiotropy

Putting this together, we are able to explain why gene pleiotropy can be effectively estimated without counting the phenotypes. Precisely, let *P*_{min} be the minimum pleiotropy over all legitimate measures. Under the nearly neutral model (Kimura 1983), the Gu-2007 method attempts to estimate *K* = min(*r*, *P*_{min}), the rank of the genotype–phenotype map, from a given MSA of protein sequences. Moreover, the effective estimate of *K* (denoted by *K*_{e}) can be interpreted as the effective gene pleiotropy (*K = P*_{min}) when *P*_{min} *< r*, or *K = r* as a minimum estimate of *P*_{min} when *P*_{min} *> r*. We also suggest, to be practically meaningful, that (i) amino acid sequence alignment (*r*_{max} = 19) is preferred, rather than nucleotide sequence alignment (*r*_{max} = 3), and that (ii) sequence sampling should be from a large phylogeny so that the real rank of mutational effects (*r*) can be reasonably close to *r*_{max}. For amino acid sequence alignment, the principle of microergodicity ensures that, when the phylogeny is sufficiently large, the rank (*r*) of mutational effects is close to the maximum value *r*_{max} = 19, leading to *K* → min(19, *P*_{min}).

In the literature, the term “mutation size” (*s*) was used as the overall mutational effect on fitness. At the molecular level, it is virtually the same as the coefficient of selection. One criticism (Razeto-Barry *et al.* 2011, 2012; Razeto-Barry and Maldonado 2011) was that the Gu-2007 method assumed that gene pleiotropy and the overall mutation size (*s*) were correlated. We try to clarify this issue in the following two arguments. First, the overall mutation size can be written as *s* = *s*_{1} + … + *s _{n}*, where

*s*is the mutation size for the

_{i}*i*th fitness component. If

*s*remains a rough constant, we have

_{i}*s*≈

*s*

_{0}

*n*, implying a strong correlation between

*s*and

*n*. By contrast, if there is only a very small number (

*n**) of large

*s*’s, but the others are very small, we expect that

_{i}*s*correlates only with

*n** rather than with

*n*. Indeed, the Gu-2007 method estimates the number of fitness components with strong mutation size. Second, we have shown that the rank of the genotype–phenotype map

*K*= min(

*r*,

*P*

_{min}), rather than the pleiotropy

*P*

_{min}, was involved in the statistical model that underlies the Gu-2007 method. Hence, the question turns out to be whether

*K*and

*s*are correlated. As there seems no correlation between the rank (

*r*) of mutational effects and the mutation size (

*s*), we claim that

*K*is not correlated with

*s*when

*r*<

*P*

_{min}. When

*r*>

*P*

_{min}, the estimated

*K*is for the minimum pleiotropy

*P*

_{min}. In practice, whether the overall mutation size (

*s*) is correlated with

*P*

_{min}may affect the interpretation of the Gu-2007 method, but not the procedure of estimation.

### Computer simulation and data analysis

#### Simulation study for *K*_{e} estimation:

We carried out computer simulations to illustrate our main conclusion. To be concise, the simulation procedure was designed under the simple one-step mapping model so that *K* = min(*r*, *n*). For each of the simulation data, the Gu-2007 method was implemented to estimate the effective gene pleiotropy (*K*_{e}). Table 1 shows the results. In case *A*, we set *r* = *n*/2 and so *K* = min(*r*, *n*) = *n*/2. It appears that the estimate (*K*_{e}) is, on average, less than but close to *r = n*/2, because the estimation of the Gu-2007 method is conserved. In case *B*, we set *r* = 15 but allow *K* to vary from 2 to 20. Although *K*_{e} is on average <*n* in all cases, it increases almost linearly with the increasing of *n* until it reaches the point of *n = r*. As expected, when *n > r*, the estimated *K*_{e}, on average, becomes quickly saturated toward the line of *r* = 15, regardless of the increasing of *n*.

#### Range of *K*_{e} estimated from protein sequences:

One impressive prediction from our theory is that, based on the protein sequence alignment, the Gu-2007 method can estimate the degree of gene pleiotropy effectively only between 1 and 19, whereas this estimate would converge to 19 when the real gene pleiotropy becomes higher. We have examined three data sets: 321 vertebrate genes (Su *et al.* 2010), 580 single-copy genes from five yeast genomes, and 437 single-copy genes from 12 *Drosophila* genomes. Although each data set shows a great range of gene pleiotropy, almost all estimates (99%) are <19 (Table 2), consistent with the theoretical prediction.

When the total evolutionary time of the phylogeny is not large enough to meet the criterion of microergodicity, the real rank (*r*) of mutational effects could be <<*r*_{max}, resulting in a lower rank (*K*) of the genotype–phenotype map. We tested this prediction by the vertebrate gene data set (Su *et al.* 2010). For each gene, we estimated *K*_{e} from four species (mammals), five species (mammals and chicken), six species (mammals, chicken, and frog), seven species (mammals, chicken, frog, and fugu), and eight species (mammals, chicken, frog, fugu, and zebrafish), respectively. We observed the average of *K*_{e} = 3.1, 4.2, 6.0, 6.3, and 6.5 for the data sets of four species to eight species, respectively. Indeed, for the same gene set, the estimation of *K*_{e} tends to be small when the protein phylogeny is small. One may wonder the ultimate value that *K*_{e} may approach when the phylogeny becomes very large. Tentatively, we used a simple extrapolation model after calculating the total branch length of the phylogeny in each case and estimated and found that, roughly, the average of *K*_{e} approaches 7.17.

*K*_{e} estimated from nucleotide sequences reflecting the rank (*r*) of mutational effects:

While the Gu-2007 method is applied to the nucleotide sequence alignment, we have *K =* min(*r*, *P*_{min}) < 3 because of the maximum rank *r*_{max} = 3 at nucleotide sites. We used the vertebrate gene set (Su *et al.* 2010) to test this prediction. Two types of nucleotide sets were analyzed: The first data set includes the first and second codon positions (data set 1st2nd), and the second one includes all sites (data set all). Since many saturated mutations may have occurred at synonymous sites, the estimation procedure in the Gu-2007 method would be statistically unreliable. We propose a modified procedure to avoid this problem; see *Materials and Methods*. We analyzed 321 vertebrate single-copy genes and found that the average of *K*_{12}, the effective gene pleiotropy based on the first and second codon positions, is 2.91, with the 25–75% sampling quantile (1.75–5.59). This finding can be reasonably explained as the result of *r*_{max} = 3 such that *K*_{12} *=* min(3, *P*_{min}) ≈ 3 despite that the real gene pleiotropy *P*_{min} may be much higher. Moreover, when all nucleotide positions are used, the average of effective gene pleiotropy (*K*_{123}) is as low as 0.48, with the 25–75% sampling quantile (0.00–1.23). Indeed, synonymous substitutions mainly at the third codon position have little functional effect (Table 3).

#### Distribution of coefficient of selection:

Although the distribution of *s* (the coefficient of selection or, broadly, mutation size), *f*(*s*), can be approximately modeled as a (negative) gamma distribution, a long-standing controversy was about how to determine the shape parameter *β*. One class of models determined that *β* is a universal constant between 0.5 and 1.0, based on the argument that gene pleiotropy and mutation size are not correlated (Orr 2000; Razeto-Barry *et al.* 2011, 2012; Razeto-Barry and Maldonado 2011). By contrast, another class of models suggested that mutation size is the sum of effects from all (*n*) fitness components, leading to the relationship *β* ≈ *n*/2 (Poon and Otto 2000; Welch and Waxman 2003; Gu 2007a,b). We compromise between these two extremes as follows. Rephrasing the derivation of Gu (2007b) with the interpretation *K* = min(*r*, *P*_{min}), we have *β* ≈ *K*/2 = min(*r*, *P*_{min})/2, which includes two special cases: (i) *β* ≈ *r*/2 if *r* < *P*_{min}, predicting that *β* could be any value between 0.5 and 2.0, and (ii) *β* ≈ *P*_{min}/2 if *r* > *P*_{min}. Our model explains most empirical distributions of *s* with *β* = 0.5–1 simply because they were from nucleotide data (Keightley 1994). In the following we briefly discuss two examples.

Piganeau and Eyre-Walker (2003) integrated the distribution *f*(*s*) into a population genetics model. They analyzed polymorphisms of 13 genes encoded in the mouse mitochondrial genome and estimated, on average, *β* = 0.86. It means that the rank of the genotype–phenotype map is *K* = 2*β* = 1.72, consistent with our prediction of *K* = min(3, *P*_{min}) < 3 for nucleotide sites. Since the short-term coalescent time of the mouse population allows virtually no multiple (two or more) mutations to appear at the same nucleotide site, the rank (*r*) of mutational effects is likely to be <2.

In the same manner, the low-pleiotropy estimation by Martin and Lenormand (2006) was not very surprising. Since their study was based on six mutation accumulation (MA) experiments, the rank of mutational effects is virtually *r =* 1 so that *K =* 1, regardless of the real mutation pleiotropy. Our explanation is consistent with the follow-up analyses from the same research group (Martin *et al.* 2007).

### Concluding remarks

The exact meaning of the term *effective* designed for *K*_{e} in the original article (Gu 2007a) was somewhat ambiguous, which is clarified as follows. First, *K*_{e} estimated by the Gu-2007 method is an estimate for the minimum pleiotropy *P*_{min}. Second, *K*_{e} is an estimate of *K =* min(*r*, *P*_{min}), the rank of the genotype–phenotype map. Third, what the Gu-2007 method attempted to estimate is the pleiotropy of amino acid sites, a conserved proxy to the true gene pleiotropy. Finally, *K*_{e} is a conserved estimate of *K* because those fitness components that are slightly affected by gene mutations have been effectively removed by the estimation procedure. One may wonder whether the Gu-2007 method can be expanded beyond single amino acid sites so that more generic states are allowable. For instance, we can develop an evolutionary model that allows coevolution of amino acid pairs or hyplotypes in the population genetics data. While it is certainly worthwhile to investigate in a future study, we suspect that the problem of microergodicity may become severe. Consequently, a huge number of sequences (or any genotype samples) or a large phylogeny is required.

We may acknowledge that our main result *K* = min(*r*, *P*_{min}) should be a principal limitation of any method to estimate pleiotropy from genetic data instead of phenotype counting. One impressive claim is that we cannot estimate the mutational pleiotropy directly from any genetic data because of *r* = 1 virtually. Nevertheless, some important implications of our study may be interesting for the general community of evolutionary geneticists. For instance, under the framework of the genotype–phenotype map, minimum pleiotropy is a measure of a gene-specific canalization, the reduced sensitivity of a phenotype to genetic changes, recognized by observing that most genetic changes leave the phenotypes invariant. Moreover, the rank of the genotype–phenotype map *K* may be the key to determining under which condition the cost of complexity may occur. Indeed, our work is only a first attempt to understand the evolution of the genotype–phenotype map.

## Acknowledgments

The author is grateful to David Waxman, Zhixi Su, Yanwu Zeng, and Wenhai Chen for constructive discussions and criticisms in the early version of the manuscript and also grateful to the reviewing editor and two anonymous reviewers for constructive comments that have helped to address the pleiotropy problem clearly.

## Footnotes

*Communicating editor: I. Hoeschele*

- Received March 28, 2014.
- Accepted May 28, 2014.

- Copyright © 2014 by the Genetics Society of America