J Mol Evol (2006) 63:635–653
DOI: 10.1007/s00239-005-0233-x

Codon Usage and Selection on Proteins
Joshua B. Plotkin,1 Jonathan Dushoﬀ,2 Michael M. Desai,3 Hunter B. Fraser4
1

Department of Biology, University of Pennsylvania, Philadelphia, PA 19104, USA
Department of Ecology and Evolutionary Biology, Princeton University, Princeton, NJ, USA
3
Department of Molecular and Cellular Biology and Department of Physics, Harvard University, Cambridge, MA 02138, USA
4
Department of Molecular and Cellular Biology, University of California, Berkeley, Berkeley, CA 94707, USA2

Received: 3 October 2005 / Accepted: 12 June 2006 [Reviewing Editor: Dr. Lauren Meyers]

Abstract. Selection pressures on proteins are usually measured by comparing homologous nucleotide
sequences (Zuckerkandl and Pauling 1965). Recently
we introduced a novel method, termed volatility, toestimate selection pressures on proteins on the basisof their synonymous codon usage (Plotkin andDushoﬀ 2003; Plotkin et al. 2004). Here we provide atheoretical foundation for this approach. Under the
Fisher-Wright model, we derive the expected frequencies of synonymous codons as a function of thestrength of selection on amino acids, the mutation
rate, and the eﬀective population size. We analyze theconditions under which we can expect to draw
inferences from biased codon usage, and we estimate
the time scales required to establish and maintainsuch a signal. We ﬁnd that synonymous codon usage
can reliably distinguish between negative selection
and neutrality only for organisms, such as some microbes, that experience large eﬀective population sizes
or periods of elevated mutation rates. The power ofvolatility to detect positive selection is also modest—requiring approximately 100 selected sites—but
it depends less strongly on population size. We showthat phenomena such as transient hyper-mutators canimprove the power of volatility to detect selection,even when the neutral site heterozygosity is low. We
also discuss several confounding factors, neglected by
the Fisher-Wright model, that may limit the applicability of volatility in practice.

Correspondence to: Joshua B. Plotkin; email: jplotkin@sas.upenn.
edu

Key words: Codon usage — Selection — FisherWright — Diﬀusion — Volatility — Protein evolution

Introduction
Background
Nucleotide coding sequences of many organisms exhibit signiﬁcant codon bias—that is, unequal usage ofsynonymous codons. Codon bias has been attributed
both to neutral processes, such as asymmetric mutation rates, and to selection acting on the synonymouscodons themselves. The most common selective
explanation of codon bias posits that synonymouscodons diﬀer in their ﬁtness according to the relative
abundances of iso-accepting tRNAs; a codon corresponding to a more abundant tRNA would be usedpreferentially so as increase translational eﬃciency(Ikemura 1981; Debry and Marzluﬀ 1994; Sorensen
et al. 1989). To a large extent, this hypothesis has
successfully explained interspeciﬁc variation in codon
usage for organisms ranging from Escherichia coli to
Drosophila melanogaster (Akashi 2001).
Recently, we have noted that codon bias in a protein sequence can also result from selection at theamino acid level, even in the absence of direct selection on synonymous codons themselves (Plotkin andDushoﬀ 2003; Plotkin et al. 2004). Codon bias arises
from selection at the amino acid level because ofasymmetries in the structure of the standard genetic

636

code. Proteins that experience diﬀerent selective
regimes should exhibit diﬀerent synonymous codon
usage. Following from this observation, we haveintroduced methods to screen a single genome
sequence for estimates of the relative selection pressures on its proteins by comparing their synonymouscodon usage, controlling for their amino acid content(Plotkin et al. 2004).
In this paper, we provide a theoretical analysis ofcodon usage biases that result from selection at theamino acid level. Our analysis, based on the FisherWright model of population genetics, provides a theoretical grounding for techniques of estimating selection pressures on proteins using signals gathered fromtheir synonymous codon usage. Throughout most of
this paper, we will ignore any source of direct selection
on synonymous codons and focus, instead, on codon
biases that result from selection at the amino acid level. To the extent that any other sources of codon bias
apply equally across the genome, we have devised arandomization method to control for these externalsources of codon bias when estimating relative selection pressures on proteins (Plotkin et al. (2004). In thediscussion, we describe a range of confounding factorsthat may vary across the genome and limit the applicability of methods to detect selection from synonymous codon usage. In Supplementary Information we
respond to criticisms of volatility.
Codon Volatility
Codon usage biases can arise from selection on proteins because synonymous codons may diﬀer in theirvolatility—deﬁned, loosely, as the proportion of acodonÕs point mutations that result in an amino acidsubstitution. Although there are several possible
deﬁnitions of volatility, which can all be informative,
we have recently used the following formal deﬁnition(Plotkin et al. 2004).
We index the 61 sense codons in an arbitrary orderi = 1,...,61. We use the notation aa(i) to denote theamino acid encoded by codon i. For each codon i, let
B(i) denote the set of sense codons that diﬀer fromcodon i by a single point mutation. We deﬁne thevolatility of codon i by
1 X
m ðiÞ ¼
D½aaðiÞ; aaðiÞð1Þ#BðiÞ j2BðiÞ
where D denotes the Hamming metric, which is zero if
two amino acids are identical, and one otherwise. The
deﬁnition in Eq. 1 applies when all nucleotide mutations occur at the same rate. When diﬀerential nucleotide mutation rates are known (e.g., a transition/transversion bias), these rates can be incorporated
into the deﬁnition of codon volatility by appropriately
weighting the ancestor codons (Plotkin et al. 2004).

Minor variants of Eq. 1 yield related deﬁnitions ofcodon volatility. For some applications, one may
want to allow termination codons in the deﬁnition of
B(i). It is also natural to consider alternatives to theHamming metric, D, that weight substitutions
between amino acids depending upon the diﬀerencesin their stereochemical properties (Miyata et al. 1979;
Plotkin and Dushoﬀ 2003). A variety of other metrics(Tang et al. 2004; Yampolsky and Stoltzfus 2004)
that reﬂect the eﬀects of diﬀerent amino acid substitutions on protein structure may likewise be incorporated into the deﬁnition of codon volatility. In thispaper, however, we will focus on the most basic
deﬁnition of codon volatility (Eq. 1, using the Hamming metric), because variant deﬁnitions are based on
the same underlying principle and produce similar
results in practice (Plotkin and Dushoﬀ 2003).
The volatility of a codon is strongly inﬂuenced by
the amino acid which it encodes. Therefore, in orderto estimate the relative selection pressures on proteins
across a genome, we have introduced a simple randomization technique to produce ‘‘volatility p-values’’
that control for the amino acid sequence of each gene
(Plotkin et al. 2004). Under the most basic deﬁnitionof volatility, there are four amino acids (glycine, leucine, arginine, and serine) whose synonymous codons
diﬀer in their volatility. As a result, when controllingfor amino acid content, we will obtain a volatility
signal from only those sites that contain one of thesefour amino acids—which amounts to roughly 30% of
the sites in a typical gene. (If one uses stereochemical
metrics for D in the deﬁnition of volatility 75% of
the sites in a gene contain a volatility signal). Although 30% may seem like a small proportion of sites
from which to obtain a signal of selective pressures, it
is larger than the proportion of sites often used todetect selection via sequence comparison of recently
diverged species (Fleischmann et al. 2002; Clark
2003). For example, fewer than 4% of neutral sites
exhibit substitutions when comparing human andchimpanzee sequences (Clark et al. 2003).
In the following sections we analyze the consequences of selection on proteins for codon usage in
general, as well as for the volatility measure in particular. Under the Fisher-Wright model, we demonstrate that the expected codon usage at a site, as wellas its temporal dynamics, depends upon the strengthof positive or negative selection on the amino acidsequence. We ﬁrst analyze the dynamics of negative
selection in inﬁnite and ﬁnite populations, and then
we analyze positive selection. Our analysis is usuallyconﬁned to the patterns of codon usage at a single site
under selection at the amino acid level. However, we
also discuss codon usage over many sites within agene or genome, and we analyze how many sites are
required in principle to detect a reliable signal of
selection by inspecting synonymous codon usage.

637

Results
Negative Selection in an Inﬁnite Population
Most nonsynonymous mutations in a protein coding
sequence presumably reduce the ﬁtness of an organism. For a large proportion of sites, therefore, natural
selection opposes any change in the amino acid. We
refer to this type of selection as negative selection. For
the purposes of exploring the eﬀect of negative
selection on codon usage, we assume that selection
cannot directly discriminate between the synonymouscodons for the favored amino acid at a site. However,mutations are more likely to be nonsynonymous, andhence deleterious, if the codon at that site has highvolatility. As we will show, this fact results in aneﬀective preference for the less volatile codons,among those codons that code for the favored aminoacid at the site. We emphasize that this preference fora codon of low volatility at a site under negative
selection is not caused by a direct ﬁtness diﬀerence
between synonyms. Rather, more volatile codons will
occur less frequently as a second-order consequenceof negative selection at the amino acid level, and thestructure of the genetic code.
Proteins containing a larger number of sites undernegative selection will exhibit a statistical bias towards less volatile codons, after controlling for theiramino acid content. In this section we calculate theexpected magnitude of the codon bias as a function of
the mutation rate, the strength of negative selection,and the population size. We also analyze the conditions under which we can expect to detect and draw
inferences from this bias, and we estimate the timescales needed to establish and maintain such a signal.
A simpliﬁed genetic code
In an inﬁnite population, we can describe thedynamics of codon usage at an individual site byusing the standard multi-allele model, based on thework of Haldane and used throughout the literature(e.g., Nagylaki 1992, Eq. 2.25; Higgs 1994). Thismodel describes a single site which can assume any of
K states. In order to investigate codon usage, we
consider K = 64 states, corresponding to each of the
64 possible codons. In continuous time, the proportion xi of individuals with codon i evolves according
to
K
X
dxi
¼
xj ðtÞwj Mij  xi WðtÞ
dt
j¼1

ð2Þ

where wjPis the Malthusian ﬁtness of codon j,WðtÞ 
j wj xj ðtÞ is the mean ﬁtness of the population, and Mij is probability that codonPj will produce codon i upon replication, withiMij = 1.

Although Eq. 2 is non-linear, the equilibrium frequencies of the codons i = 1,2,...K are given by theprincipal right eigenvector of the matrix wj Mij(Thompson and McBride 1973). These frequencies
determine the expected equilibrium codon usage at asite. A very similar model, and method of solution,
has been used to study the evolution of mutational
robustness (van Nimwegen et al. 1999; Wilke 2001).
For the purposes of this paper, alternative formulations of the K-allele model that treat the processes of
selection and mutation separately (e.g, Crow andKimura 1970, Eq. 6.4.1) yield the exact same results.
Our modeling framework is similar to that ofGolding and Strobeck (1982), who also formulated amutation-selection model for the full genetic codeand pointed out that mutations to deleterious acids
will have an eﬀect on synonymous codon usage.
The equilibrium solution to Eq. 2 for the real genetic code does not lend itself to intuitive understanding. Transient dynamics are also diﬃcult to
calculate in this high-dimensional system. Therefore,
in order to highlight the essential points of ouranalysis, we ﬁrst consider a ‘‘toy’’ genetic code thatretains those features of the true genetic code relevantto the study of synonymous codon usage under negative selection. As we will demonstrate, the solution
for the simpliﬁed genetic code yields a complete
understanding for the full genetic code as well.
We imagine a simpliﬁed genetic system with only
three possible codons, a1, a2, and b. Codons a1 and a2
code for amino acid A, which is favored, and codon b
encodes amino acid B, which has selective disadvantage r. We assume that mutations occur at rate u
between these codons according to the structure
a1 ¿ a2 ¿ b so that of the two synonymous codons,
a2 is more volatile.
According to the standard multi-allele model (Eq.
2), the relative frequencies of codons a1, a2 and b are
described by the equation0
1 0
1
a1 ðtÞ
1u
u
0
d B
C B
C
1  2uu ð1  rÞ A
@ a2 ðtÞ A ¼ @ u
dt
b ðtÞ
0
u
ð1  uÞ ð1  rÞ
0
1
01
a1 ðtÞ
a1 ðtÞ
B
C
B
C
 @ a2 ðtÞ A  W ðtÞ @ a2 ðtÞ A
b ðtÞ
b ðtÞ
ð3Þ
where W(t) = a1 (t) + a2 (t) + (1)r)b(t).
The equilibrium frequencies of codons are given by
the principal right eigenvector of the matrix in Eq. 3.
A simple perturbation analysis of this eigenvector
shows that the equilibrium proportion of a1 depends
monotonically on r, and it exhibits a sharp transition between two regimes: the weak selection regimer  u and the strong selection regime r  u. In the638

weak selection regime, the equilibrium relative frequencies of synonyms are given by the expansion
 2
1
1 r
r
a^1
þO
¼ þ
ð4Þ
12 u
u2
a^1 þ a^2 2
And in the strong selection regime, the equilibrium
relative frequencies are given by
pﬃﬃﬃ
pﬃﬃﬃ
 2
5  1 ð5  2 5Þð1  rÞ u
u
a^1

þO 2
¼
2
5
r
r
a^1 þ a^2
ð5ÞIn the absence of selection (r = 0) all three codons occur with equal frequency, as we would expect.
In particular, the relative proportion of the two
synonymous codons a1 and a2 equals 12, regardless ofthe mutation rate. For weak selection (r  u), thisresult is still approximately true, according to theperturbation expansion above. In the case of strong
negative selection (r  u), the relative proportionof the two synonymous codons is given approximately by the reciprocal of the golden meanpﬃﬃﬃ
ð 5  1Þ=2  0:62:
We note that the results above, in the limits ofstrong and weak selection, are fairly straightforward
when formulated in the language of mutational
robustness: codons coding for the favored amino acid(under strong selection) or for any amino acid (under
weak selection) form a neutral network, and the steadystate concentrations are given, to ﬁrst order, by theprincipal eigenvector of the connection matrix of thisnetwork (van Nimwegen et al. 1999; Wilke 2001)
whose perturbation expansions we have given above.
The sharp transition between the weak and thestrong selection regimes deﬁnes r = u as a criticalvalue for negative selection. For r  u negative
selection is ineﬀective at favoring the less volatile codon, and the site is eﬀectively neutral. But when r  u,negative selection favors the less volatile codon, andthe magnitude of this eﬀect depends only weakly on thevalue of r. This is an essential point. In the strong
selection regime, the magnitude of negative selection isrelatively unimportant; volatile codons are disfavoredat all sites where r  u. The transition between the
weak and strong selection regimes is shown in Fig. 1.
The eﬀective disadvantage of a volatile codon
The critical value of r discussed above can beunderstood intuitively by considering the eﬀective
selective disadvantage of the more volatile codon a2that results indirectly from its volatility. We will use
the notion of an eﬀective selective disadvantage to aid
in our analysis of codon usage at a site under negative
selection. But we emphasize that our model (Eq. 2)
does not assume any direct ﬁtness diﬀerence between
synonymous codons.

Fig. 1. The relationship between selection at the amino acid level
and resulting synonymous codon usage. The graph shows relative
a1 þ a^2 Þ as aequilibrium frequency of synonymous codons, a^1 =ð^function of the strength of negative selection, r. The relative frequency of codon a1 is approximately ½ in the weak selection repﬃﬃﬃ
gime (r  u) and approximatelyð 5  1Þ= 2 in the strong selection
regime (r  u). In this ﬁgure u = 10)5.

When the disfavored amino acid B is lethal to theorganism, then the eﬀective selective disadvantage ofcodon a2 is particularly simple to understand. In thiscase, individuals with codon a2 are removed from the
population at rate u because they mutate to the lethal
codon b, but receive no back-mutations. Hence theeﬀective selective disadvantage, denoted s, of codon
a2 versus codon a1 is given by s = u. The eﬀective
selective disadvantage of a2 does not arise from aﬁtness diﬀerence between synonyms but, rather, fromselection at the level of amino acids and the structureof the genetic code.
When amino acid B is not lethal the situation isslightly more complicated. Nevertheless, for r  u,mutations from a2 to b typically die due to negative
selection before they mutate back from b to a2. As aresult, the eﬀective selective disadvantage will still bes = u in the regime of strong selection. We can makethis argument concrete by considering the mutationselection balance between codon b and codon a2.
According to the standard mutation-selection balance, the equilibrium proportion of codon b relative
to codon a2 equals u/r in the regime r u. Thus foreach mutant from a2 to b, there are at most of orderu/r mutations from b to a2. The net mutation ratefrom a2 to b is therefore u(1 ) (u/r)). This is the rate
at which individuals of type a2 are lost from the
population due to the fact that a2 is more volatile
than a1. Thus the eﬀective selective disadvantage ofcodon a2 relative to a1 is given by s = u (1 ) (u/r)).
In the strong selection regime we neglect u/r compared to 1, and the eﬀective selective disadvantage ofcodon a2 is simply s = u.
A similar argument holds for the real genetic code.
In this case, the favored amino acid may correspondto several synonymous codons, each with a potentially diﬀerent volatility. However, the eﬀective

639

selective disadvantage, s, of a more volatile codon
relative to a less volatile synonym is simply the difference in the number of mutations leading to a disfavored amino acid ((r  u) times u/3, where u is thenucleotide mutation rate. (Note that u/3 is the rate of
mutation between any two particular nucleotides.)For example, when considering the relative frequencies of codons AGA and CGG at a site under negative selection for arginine, AGA has selective
disadvantage s ¼ 2u=3 compared to CGG, since
AGA has two more disfavored neighbors than CGG.
An analogous argument can be used to calculate
the eﬀective selective disadvantage of codon a2 in theregime of weak selection (r  u). In this regime, therelative equilibrium proportion of codon b versus
codon a2 equals 12  r=ð8uÞ. As a result, the eﬀective
selective disadvantage of a2 versus a1 is approximately s = 0, plus a small correction of order r. In
other words, when r  u selection between a1 and a2
is eﬀectively neutral; it cannot generate codon bias.
We therefore refer to the regime r  u as the ‘‘almost-neutral regime.’’ This result holds both for thesimpliﬁed three-codon model and for the real geneticcode.
It is also important to calculate the amount of timerequired to reach equilibrium codon usage in thepresence of strong negative selection. Explicit solution of Eq. 3, assuming r  u, indicates that the efold relaxation time is of order 1/u (the selection
coeﬃcient is s  u, and so the time scale for population sizes to change under selection is of order 1/s 
1/u). In other words, starting from any initial frequencies a1(0) and a2(0), these frequencies will become e-fold closer to their equilibrium values after aduration of order 1/u generations. The same timescale holds for almost-neutral sites (r  u) and forthe real genetic code. In practice, u will be quite small,and equilibrium volatility is approached very slowly.
We will revisit this point when we discuss ﬁnite
populations and, again, when we discuss positive
selection.
A speciﬁc example of negative selection
In this section we consider a simple example that
demonstrates how our analysis applies to the real
genetic code. We use Eq. 2 to model the dynamics of
K = 64 alleles corresponding to the 64 codons, indexed in an arbitrary order. For our example, we
consider a single site under negative selection for an
arginine codon. In this case we deﬁne8< 1  3u if i ¼ j
Mij ¼ u=3
if i and j differ by a point mutation:0
otherwise
ð6Þ

where u is the nucleotide mutation rate. We deﬁne
8
1
if i encodes arginine>><
1  r if i encodes a non-arginine amino
wi ¼
acid
>
>
:
1  c if i encodes stop
ð7Þ
so that a codon encoding an amino acid other thanarginine has ﬁtness 1 ) r, and a termination codon
has ﬁtness 1 ) c. We analyze this model numericallyby calculating the principal right eigenvector of thematrix wjMij, which yields the equilibrium proportions of all 64 codons.
In the case of no selection (r = c = 0), we ﬁndthat all codons occur with the same equilibrium frequency, independent of mutation rate, as we would
expect. For almost-neutral selection (r  c  u),codon usage is still approximately uniform. In the
opposite case when arginine is favored and all otheramino acids (or termination codons) are strongly
disfavored (i.e., r  c  u), the arginine codons
CGA, CGG, CGC, CGT, AGA, and AGG occur
withequilibrium
relative
proportions

0.214:0.214:0.191:0.191:0.095:0.095. As expected,
under negative selection the more volatile argininecodons occur with lower relative frequency in equilibrium.
The equilibrium frequencies of arginine codons
determine the expected volatility at a single arginine
site under negative selection. Assuming free recombination (Sawyer and Hartl 1992), an individual geneconsists of many such sites randomly assembled; themean and standard deviation in the volatility (persite) of a randomly sampled gene are shown in Fig. 2,
as a function of the strength of a negative selection r.
Note that the stronger the negative selection, thelower the expected equilibrium volatility. The expected volatility exhibits a sharp transition from highto low values when the strength of negative selectionr reaches the mutation rate u, as discussed above. On
either side of this transition, the volatility is insensitive to r. The standard deviations plotted in Fig. 2
correspond to a gene comprised of L = 200 arginine
sites, each modeled independently by the multi-allele
equation.
According to Fig. 2, L = 200 independent arginine sites that each experience neutrality (r  u) canbe distinguished on the basis of their volatility from L= 200 sites that experience negative selection (r u). The diﬀerence in the expected volatility between
these two regimes is greater than four standard
deviations of the volatility within either regime.
Similar results hold for serine and leucine, but glycine
is less informative (Table 1).
In reality, the selective constraint r will vary
greatly across the sites of a given protein. In this case,640
Table 1. Equilibrium codon usage under neutrality versus
negative selection

Fig. 2. The relationship between selective strength and volatility
for a gene comprised of L = 200 freely recombining sites underselection for arginine. The graph shows expected volatility per site
in the gene (±1 standard deviation, dashed) as a function of thestrength of negative selection, r. The nucleotide mutation rate isu = 10)5. The expected volatility is signiﬁcantly depressed in theregime of strong negative selection, r  u. (For this ﬁgure weassume c = 1; virtually identical results hold for c = r).

disregarding the possibility of positive selection, thevolatility of a gene (after controlling for its aminoacid sequence) essentially reﬂects the relative numberof informative sites that experience negative selection
versus neutrality. For example, the volatility of geneX that contains L = 200 informative sites undernegative selection and an equal number of neutral
sites will be signiﬁcantly greater (with a Z-score ofabout three) than the volatility of gene Y that consistsof 2L informative sites all under negative selection. A
more thorough discussion of variable selection pressures across genes is described in a subsequent section(Eq. 18), below.
Table 1 shows the equilibrium relative frequencies of synonymous codons for each of the informative amino acids (G, L, R, and S) under
neutrality versus negative selection. In Table 1 weassume, as we do throughout this paper, that volatility is measured using the Hamming metric andthat there is no transition/transversion bias. Corresponding values for diﬀerent metrics or including
a mutational bias may be calculated using the same
approach. As shown in Table 1, the diﬀerence inthe expected volatility between selective regimes isleast extreme (indeed, barely informative) for glycine sites. The volatility diﬀerence is most extreme
for serine sites: the highly volatile codons AGT and
AGC are not expected to occur at a site undernegative selection, but they are expected occur at asite under neutrality or, especially, under positive
selection (see below). This extreme case results fromthe fact that codons AGT and AGC are not connected by synonymous point mutations to the otherserine codons. As a result, the relaxation time toequilibrium for codons AGT and AGC is longer
than 1/u generations, and these codons may be lessinformative for negative selection because they takeLeucine
cta
ctc
ctg
ctt
tta
ttg
E [m]
r [m]
Arginine
aga
agg
cga
cgc
cgg
cgt
E [m]
r [m]
Serine
agc
agt
tca
tcc
tcg
tct
E [m]
r [m]
Glycine
gga
ggc
ggg
ggt
E [m]
r [m]

NeutralNeutral*Negativem

0.16667
0.16667
0.16667
0.16667
0.16667
0.16667
0.65146
0.073620.17300
0.185800.178900.185800.129900.14650
0.645900.07259

0.21353
0.19098
0.21353
0.19098
0.09549
0.09549
0.63172
0.07022

5/9
6/9
5/9
6/9
5/7
6/8

0.16667
0.16667
0.16667
0.16667
0.16667
0.16667
0.65278
0.098540.152100.17050
0.152100.177400.17050
0.177400.65400
0.09660

0.09549
0.09549
0.21353
0.19098
0.21353
0.19098
0.62592
0.09354

6/8
7/9
4/8
6/9
5/9
6/9

0.16667
0.16667
0.16667
0.16667
0.16667
0.16667
0.71792
0.12504

0.185100.185100.134400.171900.151620.171900.729810.125610.00000
0.00000
0.25000
0.25000
0.25000
0.25000
0.632430.03913

8/9
8/9
4/7
6/9
5/8
6/9

0.25000
0.25000
0.25000
0.25000
0.656250.01804

0.22460
0.261800.251700.261800.657240.018590.25000
0.25000
0.25000
0.25000
0.656250.01804

5/8
6/9
6/9
6/9

Note. In each regime, we report the equilibrium relative abundanceof codons, and the resulting mean and standard deviation involatility per site. The second column corresponds to neutrality(r = c  u): the third column corresponds to neutrality but with
disfavored termination codons (r  u, c = 1); the fourth column
corresponds to negative selection in an inﬁnite population (r  u, c u). The ﬁnal column gives the volatility of each codon, assuming
no transition/transversion bias (Plotkin et al. 2004).

too long to be eﬀectively selected against. Thissituation does not necessarily imply that codons
AGT and AGC should be treated separately fromthe other serine codons. In fact, when treated as anentire group, the serine codons are particularly
informative for discriminating between positive
selection and neutrality (Table 2).
Negative Selection in a Finite Population
The models presented in the previous section describe
the processes of mutation and negative selection in an
inﬁnite population. In ﬁnite populations, however,genetic drift also aﬀects allelic frequencies. In thissection, we study the combined eﬀects of mutation,

641
Table 2. The mean and standard deviation of volatility at a site under various models of positive selection

Arg
E [m]
r [m]
Leu
E [m]
r [m]
Ser
E [m]
r [m]
Gly
E [m]
r [m]

Pos1

Pos2

h=1

h = 0.1

h = 0.01

0.66830.09410.6607
0.09270.66840.08190.66830.08890.66750.08960.65540.07220.64840.06980.6541
0.06360.65320.06910.65380.06840.74680.12710.73370.12680.73720.12750.73960.13010.74220.13010.65760.01720.65760.01720.65760.01630.65760.01710.65760.0171

Columns 2 and 3 summarize the expected volatility after a positively selected sweep according to heuristics 1 and 2 described in the maintext. Columns 4 through 6 report simulations of the full Fisher-Wright model in a population of N = 1000 individuals. For each simulated
selective sweep, the population is initially ﬁxed for a random, disfavored codon. Codons for the favored amino acid are assigned ﬁtness 1,
and all others ﬁtness 1 ) r (here r = 0.1 but results do not depend on r provided r  1/N and r  u). After the favored amino acid has
swept the population (proportion > 99%), the average volatility at the site is recorded. Results reported for each amino acid and h value
represent the average and standard deviation obtained in >2000 independent runs. As discussed in the main text, the appropriate eﬀective
population size that determines the value of h does not necessarily equal the average neutral site heterozygosity.

negative selection, and drift, which we analyze using
diﬀusion equations. These equations can be very
complex. A full treatment of even the simpliﬁed
three-codon genetic code requires a two-dimensional
diﬀusion process, and the real genetic code involves a63-dimensional process. To make this problem tractable, we use the notion of the ‘‘eﬀective selective
disadvantage’’ of more volatile codons, discussed
above. This allows us to consider the dynamics onlyat the favored codons, thereby reducing the dimensionality of the diﬀusion process.
The neutral (r = 0) or almost-neutral (r  u)regimes are straightforward: here all synonymouscodons for the favored amino acid have the sameeﬀective ﬁtness. In this regime, each synonymouscodon occurs with the same probability in steadystate, independent of population size.
For the remainder of this section, we analyze thecase of strong negative selection (r  u) at a single
site. We assume that Nr  1, so that selection at theamino acid level is eﬀective. (This assumption is not
restrictive, because when Nr[1, then either h = 2Nuor r will be too small to support a volatility signatureof selection.) We consider a diﬀusion approximation
to the process of mutation, selection, and drift
operating only on the synonymous codons, to each ofwhich we assign an eﬀective selective coeﬃcient. For
the simpliﬁed three-codon genetic system, the morevolatile codon a2 has an eﬀective selective disadvantage of s = u compared to codon a1. For the real
genetic code, more volatile codons will have a selective disadvantage of this order, but the precise valueof s will depend on the speciﬁc codon in question. In
the following analysis, we consider the case of thesimpliﬁed three-codon system. However, we do not
explicitly make the substitution s =u, so that our

results can also be applied (with a slightly diﬀerent
value of s) to the real genetic code.
The time-dependent proportion f(x,t) of allele a1
relative to allele a2 can be described in the diﬀusion
limit of the Fisher-Wright model by the Komolgorov
forward equation (Kimura and Crow 1964),@fðx; tÞ@1 @2
¼  faðxÞfðx; tÞg þfbðxÞfðx; tÞg
@t
@x
2 @x2
ð8Þ
where the instantaneous mean and variance in thechange of allelic frequency are given byaðxÞ ¼ sxð1  xÞ  ux þ uð1  xÞbðxÞ ¼ xð1  xÞ=N^
The stationary distribution of allele frequencies fðxÞsatisﬁes the equation@^^fbðxÞfðxÞg
¼ 2aðxÞfðxÞ
@x

ð9Þ

which has the solution (Wright 1931)^ ¼ Cxh1 ð1  xÞh1 eSx
fðxÞ

ð10Þ

h = 2Nu, S = 2Ns, and C is chosen so thatRwhere1^fðxÞdx¼ 1. Since s  u (and thus S  h), the shape0^ falls into two catof the stationary distribution fðxÞegories: a bell-shaped distribution in the regime h >1
and a U-shaped distribution in the regime h < 1. Inother words, for h > 1 the steady-state population istypically polymorphic at the locus, much like theinﬁnite population mutation-selection balance. In
contrast, for h < 1 the steady-state population isusually near-monomorphic at the locus, occasionallyswitching between alleles a1 and a2, with a bias
(whose strength is determined by S) toward allele a1.

642

In stationary state, the expected proportion ofallele a1 is given byZ 1
1 Bðh þ 1=2; S=2Þ^ð11ÞxfðxÞdx¼ þ
Mðh; SÞ ¼
2 2Bðh  1=2; S=2Þ0experience a direct ﬁtness diﬀerence. In the limit ofsmall h, we ﬁnd that

where Bðx; yÞ is the modiﬁed Bessel function of theﬁrst kind. Similarly, the variance in the frequency ofallele a1 is given byZ 1^Vðh; SÞ ¼
 MðhÞ2x2 fðxÞdx
ð12Þ

which agrees with BulmerÕs result (his Eq. 6). In otherwords, BulmerÕs approximation applies only in thelimit of small mutation rates or population sizes. Our
results also agree with prior work on mutational
robustness: for h  1, the population assumes theinﬁnite-population solution, and for h  1 all nodes
are populated with equal frequencies (van Nimwegen
et al. 1999; Wilke 2001). Equation 11 extends earlier
work on mutational robustness by providing ananalytic solution for all intermediate values of h.
We can again use the standard Taylor expansionof the Bessel function to obtain a simple expression
for the variance in the stationary frequency of allele
a1,

0

1
þ
4 þ 8hf½2h B ðh  1=2; S=2Þ B ðh þ 3=2; S=2Þwhere; Vðh; SÞ ¼

ð1 þ 2hÞ B ðh þ 1=2; S=2Þ2 g
f½ð4 þ 8hÞ B ðh  1=2; S=2Þ2 g

ð13Þ

We use the standard Taylor series expansion of
Bðx; yÞ to obtain a simple approximation for themean stationary frequency of allele a1:Mðh; SÞ ¼

1 S
þ þ Oðh2 Þ
2 4

ð14Þ

valid for h  S  1. Thus the diﬀerence in expectedvolatility at a site under neutral versus negative
selection is of order S, when h  1.
When h = S = 1, the mean stationary frequencyof allele a1 assumes the value 1/(e )1)  0.58. For h 
S  1, the mean frequency quickly approaches
the
pﬃﬃﬃ

asymptotic value limh ﬁ ¥ M (h, h) = ð 5  1Þ 2,
in agreement with our earlier result for an inﬁnite
population.
The results in this section generalize our analysisof an inﬁnite population. For an inﬁnite population,
we found that the expected relative frequency of codon a1 equals
pﬃﬃﬃ 1=2
 in the almost neutral regime, and it
equals ð 5  1Þ 2 in the strong selection regime. In aﬁnite population with h  1, the same results hold. Ina ﬁnite population with h  1, the expected relative
frequency of the more volatile codon equals 12 in theneutral regime, and it equals 1=2 + (Ns/2) in thestrong selection regime. For any population size, therelative frequency of codon a1 depends monotonically on the strength of selection at the amino acidlevel r (Eq. 11), and it exhibits a sharp transition at
the critical value r = u. Note that the eﬀective
selective disadvantage s of codon a2 versus a1 is due
solely to the lower ﬁtness of amino acid B, and thus it
depends monotonically on r.
It is worth noting that our exact expression (Eq.
11) for the mean stationary frequency of allele a1
generalizes earlier work by Bulmer (1991) on therelative frequency of two synonymous codons that1Bð1=2; S=2Þ1¼limMðh; SÞ ¼ þ
2 2Bð1=2; S=2Þ 1 þ eS

h!0

Vðh; SÞ 

ð3 þ 2hÞð4 þ 8hÞ  3S2
16ð3 þ 2hÞð1 þ 2hÞ2

ð15Þ

ð16Þ

which is a highly accurate approximation for all h,provided as usual that S is of order h or smaller. Notethat when h  1 the variance is approximated by1
4  ðh=2Þ, and when h  1 the variance is of order 1/h.
Inferring negative selection in a ﬁnite population
Our exact (Eq. 11) or approximate (Eq. 14) expressions for the stationary mean frequency of codon a1
allow us to determine the minimum number of sites
required for codon volatility to distinguish between
neutrality and negative selection. When sites aremodeled independently (equivalent to the assumptionof linkage equilibrium [Sawyer and Hartl 1992]),
under neutrality (r  u; s = 0) the relative frequency
of codon a1 versus codon a2 across a gene of length L
is binomially distributed with mean 12 and variance 1/
4L. If, on the other hand, the gene experiences negative selection (r  u; s = u), then the relative frequency of codon a1 is binomially distributed, withmean M(h, S) and variance M(h, S)[1 ) M(h, S)]/L.
Therefore, in order to reject neutrality at the 95%
conﬁdence level with 50% power, we requirerﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1
1=4 þ Mðh; SÞ½1  Mðh; SÞMðh; SÞ  > 1:96
2
L
ð17Þ
Using this equation, Fig. 3 shows the minimum
number of sites required to distinguish between negative selection and neutrality on the basis of codon
volatility, under our simpliﬁed ‘‘genetic code’’ consisting of three codons.

643

Fig. 3. The relationship between the scaled population size,
h = 2Nu, and the minimum number of sites required to distinguishnegative selection from neutrality, at the 95% conﬁdence level. Sites
are assumed to be unlinked. It is important to note that theappropriate eﬀective population size that determines the value of hdoes not necessarily equal the average neutral site heterozygosity(see text).

Nevertheless, provided J  L, the power to discriminate between gene X and gene Y is decreased by afactor of two at most (compared to J = 0), and so theminimum number of sites required to detect negative
selection (Fig. 3) would change by a factor of two at
most.
The results in this section were derived for a simpliﬁed genetic code, whereas the real genetic code
would require analysis of multi-dimensional diﬀusion
equations. The scaling behavior of the one-dimensional case is expected to hold in higher dimensions—
i.e., for the real genetic code when comparing neutrality to negative selection, for h  1 the expected
diﬀerence in volatility per site will be of order h, andfor h  1 the expected diﬀerence in volatility can becalculated from the inﬁnite population model (Eq. 2
and Table 1).
Relaxation toward steady state

Under the real genetic code, about one-third of thesites in a typical gene are informative, to varying
degrees. Thus, the y-axis in Fig. 3 would typically be
increased by roughly a factor of three, depending
upon the amino acid content. Therefore, within anorder of magnitude, Fig. 3 reﬂects the power of thevolatility technique for detecting negative selection.
Figure 3 indicates that, according to the FisherWright model, synonymous codon usage will notprovide a robust signal of negative selection on
individual genes except for micro-organisms thatexperience relatively large eﬀective population sizes(or mutator episodes, see below).
Equation 17 applies when comparing a collectionof neutral sites against a collection of sites undernegative selection. In most situations, however, theselective constraint r will vary across the sites of aprotein. For example, consider gene X, with L + J
sites under negative selection, compared to gene Y,
with L neutral sites and J sites under negative selection. In this case, the expected frequency of codon a1
in gene Y is (L/2 + JM(h, S))/(L + J). Therefore, inorder to infer that gene X experiences more negative
selection than gene Y, at the 95% conﬁdence level werequireL=2 þ JMðh; SÞðL þ JÞ
sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
L=4 þ ðL þ 2JÞMðh; SÞ½1  Mðh; SÞ> 1:96
ðL þ JÞ2

Mðh; SÞ 

ð18Þ
As Eq. 18 indicates, the power to discriminate between two genes is decreased when both genes contain many sites under the same selective regimes andonly a few sites under diﬀerent selective regimes.

Although Eq. 11 predicts the steady-state relative
frequencies of codons a1 and a2 in the selected regime(r  u), we have not yet discussed how long it takes,
on average, to reach this steady state. In the case of avery large population, h  1, we know from theinﬁnite population model (Eq. 2) that the e-foldrelaxation time to equilibrium is of order 1/u generations. In this section, we demonstrate that the same
result applies to the time scale of relaxation toward
steady state in the regime h  1.
Again, we consider a single site under negative
selection. In the regime h  1, we have seen that thesteady-state population will spend most of the timein a nearly monomorphic state, with a preference (oforder S) for the less volatile codon, a1. Therefore, inorder to calculate the time scale of relaxation toward steady state, we may simply calculate theamount of time required such that, starting with apopulation ﬁxed for allele a2, the probability of the
population remaining ﬁxed for allele a2 has been
reduced e-fold.
Given a population initially ﬁxed for codon a2,
there are Nu mutations to codon a1 generated per
generation. Each of these mutations has an eﬀective
selective advantage s=u over allele a2 and will
therefore ﬁx with probability 2s/(1 ) e)2Ns) (Crowand Kimura 1970). Hence the rate of production of amutation that will eventually ﬁx is given by
Pfix ¼

2Nus u;
1  e2Nsð19Þ

assuming h  1. According to this calculation, themean time until ﬁxation of codon a1 is of order 1/u
generations, which gives the time scale of relaxation
to the steady-state codon usage in a ﬁnite population
under negative selection.

644

About Population Sizes
As discussed above, the strength of the signal of
negative selection depends upon the parameter h. h is
twice the product of the eﬀective population size andthe (per-site) mutation rate: h = 2Nu. What is theappropriate value of h in practice?
Unfortunately, this question is far easier asked
than answered. Population geneticists have struggled
to reconcile estimates of h deduced from polymorphism data with direct measurements of N and u
across broad taxonomic ranges. The eﬀective population sizes of micro-organisms, in particular, aretopics of active debate. Estimates of h are usuallyobtained by comparing polymorphism data at neutral(or presumably neutral) sites against the expected site
diversity or the expected number of segregating sites
under a neutral model (Ewens 2004). The assumptionof synonymous site neutrality will typically lead tounderestimates of h. In one of the few cases where
neutrality was not assumed, h was estimated as 0.18
per site in E. coli (Hartl and Sawyer 1994). In a recent
survey where neutrality was assumed (Lynch andConery 2003) the authors reported an average valueof h  0.15 among the prokaryotes studied. But
estimates of h for a microbial species can vary by fourorders of magnitude, and they depend strongly on
assumptions about population structure (Berg 1996).
To complicate matters further, heterogeneity inmutation rates leads to substantial underestimates ofh (Tajima 1996).
Aside from uncertainty in its estimation, the valueof h deduced from neutral SNP data (Lynch andConery 2003) may not be relevant to questions of
selection and volatility. Monomorphism observed atneutral sites may result from non-neutral processes,
such as background selection (Charlesworth et al.
1993) or hitchhiking on periodically sweeping sites(Maynard Smith and Haigh 1974). As a result, thevariance eﬀective population size estimated from
SNP data may not be relevant to other aspects ofevolution, such as substitutions at linked weakly selected sites (Gillespie 2001).
One particularly striking example of a discrepancy
in the appropriate eﬀective population sizes arises
from the consideration of mutator phenotypes. Populations of microbial species periodically experience atransient increase in the mutation rate, usually 100–
1000 times greater than that of a non-mutator strain(Bjedov et al. 2003). Between 2 and 36% of bacterial
populations isolated in the wild at any given time
exhibit a mutator phenotype (Giraud et al. 2001;
Oliver et al. 2000; LeClerc et al. 1996). The mutatorphase can be induced in several ways. A defective
DNA repair gene (e.g., mismatch repair gene) may
arise and sweep to ﬁxation by hitchhiking on a positively selected mutation. The entire population thenexperiences an elevated mutation rate until a nonmutator allele sweeps and replaces the mutator(Notley-McRobb et al. 2001; Denamur et al. 2000). A
second, perhaps more common, mechanism is stressinduced mutagenesis; natural isolates of E. coli often
experience an increase in their mutation rate in response to stress (Bjedov et al. 2003). As a result of
these and other observations, researchers have argued
that bacterial populations evolve primarily by periodic acquisition of mutator phenotypes followed byadaptive sweeps and subsequent loss of the mutator(Giraud et al. 2001; Denamur et al. 2000; NotleyMcRobb et al. 2001). For example, Denamur et al.
(2000) concluded that ‘‘a signiﬁcant fraction of
genomic diversity of E. coli natural isolates was
generated during the MMR-deﬁcient intermezzos of
their evolutionary past’’; and Tenallion et al. (2004)
concluded that ‘‘mutations produced by stress-induced mutator alleles can represent a large fraction of
the mutations produced by a clone during its evolution.’’ As we shall see, the eﬀect of mutators on
synonymous codon usage is dramatic: the expectedsite diversity is driven by the value of h in the wild
type regime (hw = 2Nuw) but the pattern of synonymous codon usage at a site under negative selection isdriven by the value of h in the mutator regime(hm = 2Num  hw).
As a simple example of this phenomenon, we
have simulated a Fisher-Wright model of a single
locus in a population of constant size N = 1000.
The simulated site is subject to recurrent mutation
between ‘‘alleles’’ a1 and a2 at wildtype rate
uw = 10)5. As above, the alleles a1 and a2 diﬀer inﬁtness by s, where s equals the mutation rate. In
each generation, we record the frequency, x, of allele
a1. Periodically, we model the ﬁxation of a mutator
allele (or, equivalently, the stress-induced mutagenesis across the entire population) by exogenously
increasing the mutation rate to um = 103 · uw for
100 generations; thereafter we enforce a selective
sweep at the site, followed by reversion to thewildtype mutation rate. Overall, the population
experiences the mutator regime for 5% of the time,
consistent with observed frequencies of mutator
phenotypes in the wild (Giraud et al. 2001; Oliver et
al. 2000; LeClerc et al. 1996). According to oursimulations, the average site diversity, 2x(1 ) x), ata randomly chosen time equals 0.028, which is close
to its expected value assuming that h is given by hw:
E [2x(1 ) x)]  hw= 0.02. But the average frequency of allele a1 equals 0.611, which is close to its
expectation assuming that h is given by hm:
E[x] = M(hm, hm) = 0.616 (Eq. 11). In other words,
the average frequency of the less volatile codon a1 isdominated by the mutator periods, but the averagesite heterozygosity (and any estimate of h based onit) is dominated by the non-mutator periods.

645

There is a simple, intuitive explanation for thisresult. The average heterozygosity at the site is low at
virtually all times (except during the brief mutator
periods) because selective sweeps cause monomorphism, followed by long periods of low h. Therefore,
the eﬀective h for SNP diversity is small, i.e., close to
Nuw. But the site converges quickly toward the lessvolatile codon during the mutator periods, since therate of convergence is determined by s=um. And thesite is essentially frozen during the non-mutator
periods, since the decay rate of volatility is only uw.
Therefore the expected frequency of a1 at a random
time is primarily determined by the frequency reachedduring the mutator regime. As is clear from this
explanation, the expected frequency of codon a1 will,in general, depend on the stochastic scheduling of
mutator periods. For example, the site will converge
toward M(hm, hm) provided the population experiences at least one mutator phase of duration of order
1/um generations, within every 1/uw generations. In
fact, even if the mutator phases are very brief andinfrequent, the average frequency of allele a1 cangreatly exceed the value predicted by h estimatedfrom the neutral site heterozygosity. This result
indicates that, for organisms that experience mutator
periods, the value of h relevant to volatility may be
100- to 1000-fold greater than the value of h generally
reported in the literature, because the latter value isusually estimated from neutral site heterozygosity.
Although the simple model used in this section
does not describe any but the most phenomenological features of mutator alleles, it does reveal animportant general observation: the value of h estimated from presumably neutral polymorphism data
does not in general equal the eﬀective value of h that
determines synonymous codon usage at a site undernegative selection. This result is of utmost importanceto any discussion of the relationship between h andthe power of volatility to detect selection.
Positive Selection
In the sections above, we have considered selectionthat opposes a change to the amino acid at a site.
Under the Fisher-Wright model, we have seen thatnegative selection induces a bias toward the lessvolatile codons for the favored amino acid at a site.
However, selection sometimes favors a change in theamino acid at a particular site. As we will demonstrate, after a positively selected sweep has occurred,a site is more likely to be occupied by a codon ofgreater than average volatility (controlling for theamino acid now present at the site).
A variety of mechanisms are known to cause
positive selection, such as exogenous changes in theenvironment or frequency dependent ﬁtnesses. We do
not here model all of the various types of positive

selection but, rather, focus on the essential aspectshared by these mechanisms: pressure for an aminoacid substitution at a site. One common model ofpositive selection (Goldman and Yang 1994; Yanget al. 2000) stipulates continual pressure to change
the amino acid at any and all sites in a gene as quickly
as possible: as soon as an amino acid substitution has
occurred at a site, there is immediately pressure foranother substitution. Over the time scales of ourinterest, we view this scenario as unlikely. Instead, we
will study a more straightforward conception of
positive selection: we analyze the dynamics at a single
site that has, for a period of time, experienced selection for amino acid A, and that subsequently experiences selection for a diﬀerent amino acid, B (forwhatever reason). We refer to the change in thefavored amino acid at the site as a positive selection
event.
Prior to the onset of positive selection, amino acidA is assigned ﬁtness 1 and all other amino acids ﬁtness 1 ) r; subsequently, amino acid B is assigned
ﬁtness 1 and all others ﬁtness 1 ) r. We assume thatNr  1 (otherwise, the site is eﬀectively neutral at theamino acid level) and that r  u (otherwise, the expected codon frequencies are uniform). At somepoint after the change in selective regimes, a codon
for amino acid B will eventually arise and quickly
sweep close to ﬁxation. We are interested in calculating the relative probabilities of the various B
codons after the positively selected sweep has occurred. Whatever B codon arises and ﬁxes, we can be
sure that it will have arisen through a non-synonymous mutation, since we have assumed that the
population was initially dominated by amino acidA „ B. Intuition would therefore suggest that ahighly volatile B codon (which has more non-synonymous neighbors) is more likely to occur immediately after the sweep than a less volatile B codon(which as has fewer non-synonymous neighbors). In
other words, intuition would suggest that the expected volatility at the B-site is elevated immediately
following the selective sweep. Subsequent to thesweep, since B is now favored, negative selection will
reduce the volatility at the site. However, this process
takes time. Thus, for some time after the positive
selection event, we expect to ﬁnd a bias toward elevated volatility at the site, which gradually decays. In
this section, we analyze this process under the FisherWright model.

Positive selection under a simpliﬁed genetic code
Similar to previous sections, we initially consider asimpliﬁed genetic code consisting of four codons, a1,
a2, b1, and b2, the ﬁrst two of which encode aminoacid A, and the latter two amino acid B. Mutations

646

can only occur between codons a1 and a2, a2 and b2
and b2 and b1, creating the mutation structure
a1 ¿a2 ¿b2 ¿b1
In this simpliﬁed genetic code, codons a2 and b2 are
the more volatile codons for their respective aminoacids.
After the change in selection from amino acid A to
B, a mutation to codon b2 that survives stochastic
drift will eventually arise. Thus, at least initially, themore volatile codon b2 is more prevalent than the lessvolatile codon b1. During this period, we can detect
the signature of the positively selected sweep becauseof the elevated volatility at the site. However, negative selection for amino acid B will eventually favorcodon b1. Therefore, the volatility signature of thepositive selection event will be present provided thatthe time scale of decay toward codon b1 is longer than
the interval since the positive selection event.
Fortunately, the time scale of decay toward b1 isquite long. For h  1, we can use the inﬁnite population model to ﬁnd this time scale: As discussed
above, the time required to reduce the volatility e-fold
is of order 1/u. For h  1, we must use a ﬁnite
population size calculation. In this regime, the population is nearly monomorphic at almost all times.
Following the selective sweep, the site will be monomorphic for b2 with almost unit probability. We areinterested in the duration of time required such thatthe probability of being monomorphic for b2 (asopposed to b1) has been reduced e-fold. The probability of switching between b2 and b1, however, is oforder u per unit time (even before b2 has ﬁnishedoutcompeting a2), according to Eq. 19. Thus, the timescale of decay in a ﬁnite population is also 1/u.
According to this analysis, a selective sweep will
result in the presence of a more volatile codon fororder 1/u generations—a very long time indeed. (Inthe case of E. coli, for example, there are estimatedonly to be only 100–300 generations per year [Gibbons and Kapsimalis 1967; Ochman et al. 1999].)
Equivalently, repeated sweeps for amino acid changes
at a site will result in the presence of more volatilecodons at almost all times, provided that new sweeps
occur more often than every 1/u generations. It is
important to recall, however, that an episode of elevated mutation rates following a selective sweep will
induce a faster decay rate of elevated volatility.

Positive selection under the real genetic code
The analysis above for a simpliﬁed genetic code may
be generalized to the real genetic code, for arbitrary
amino acids A and B. We do not assume that A and B
are separated by a single point mutation. After achange in selective regimes from A to B, the popu-

lation will eventually acquire a B codon that survives
drift and ﬁxes in the population. As above, we areinterested in the relative probabilities of each of the Bcodons—and thus the expected volatility at thesite—immediately after the selective sweep has occurred. Under the Fisher-Wright model applied to
the real genetic code, the relative probabilities of theB codons following the selective sweep are diﬃcult to
calculate analytically. Instead, we will introduce twoheuristic approximations, and compare theseapproximations to simulations of the Fisher-Wrightmodel.
Our ﬁrst heuristic assumes that the relative probability of occurrence of a particular B codon after aselective sweep is proportional to the number of its
non-synonymous non-termination neighbors, compared to the other B codons. For example, codon
AGG (Arg) has 7 non-synonymous non-terminationneighbors; whereas codons AGA, CGA, CGC, CGG,
and CGT have 6, 4, 6, 5, and 6 such neighbors,
respectively. Thus, the probability that AGG will
dominate the population immediately after a selective
sweep for arginine is given by 7/(7 + 6 + 4 + 6 + 5+ 6) according to heuristic 1. The intuition behindheuristic 1 is simple: we know that the mutation
giving rise to the selective sweep was non-synonymous, and there are more non-synonymous routes
into certain arginine codons than into others.
We also consider a second heuristic that makes adiﬀerent simplifying assumption. Under heuristic 2,
we assume that each of the possible non-synonymouscodons within a point mutation of amino acid B isequally likely to have given rise to the B codon thateventually sweeps the population. In addition, if apotential ancestor codon has more than one way ofmutating into B, we assume that each of thesemutations was equally likely. For example, there are
27 non-termination, non-arginine codons within apoint mutation of arginine. According to heuristic 2,
we assume that each of these possible ancestor codons was equally likely to precede the arginine sweep.
Thus the probability that codon AGA (Arg) will
sweep is given by 1/27 · (1 + 1 + 1 + 1 + 1/3 + 1/3), where the various terms represent the probabilities
of a mutation to AGA from ATA, ACA, AAA,GGA, ACT, and AGC respectively. The fractions inthis expression arise because codons AnA and GGA
each have only one way of mutating to arginine,
whereas codons AGT and AGC each have threeways, which we have assumed to be equally likely.
Heuristic 2 can likewise be applied to calculate theprobability of each arginine codon after the sweep,
and thus the expected volatility at the site (Table 2).
The diﬀerence between the two heuristics is subtle,
but important. Heuristic 1 assumes that each possible
mutational route into the new amino acid is equallyprobable, whereas heuristic 2 assumes that each

647

possible ancestor codon is equally probable. Neither
heuristic yields a completely accurate description of
the Fisher-Wright model. In a small population,
however, we might expect heuristic 2 to be a moreaccurate approximation, because the population willusually be almost monomorphic before the sweep, as
it drifts between non-B codons. As soon as a substantial number of individuals drift to a B-adjacentcodon, however, one of the adjacent B codons will ﬁx
in the population long before the population drifts
into other non-B codons. Thus each ancestral codon
is equally probable. If a non-B ancestor has more B
neighbors, this does not make it a more likely
ancestor to the B codon: it simply chooses one ofthese paths with equal probability. On the other
hand, in a large population we expect heuristic 1 to bea more accurate approximation of the Fisher-Wrightmodel. In this case, the population will tend to bepolymorphic for non-B codons before the sweep. All
mutational routes into B will be used. Thus all possible mutations yielding a B codon are equallyprobable. If a non-B ancestor has more B neighbors,
it will produce more mutations to B codons.
Neither heuristic, however, yields a completely
accurate description of the expected volatility at a site
after a positively selected sweep (see Table 2). Bothheuristics neglect second- and higher-order eﬀects
that determine volatility after a sweep under the fullFisher-Wright model. For instance, if amino acid B is
two mutations away from amino acid A, then therelative probabilities of various B codons after the
sweep will depend not only on the number of their
non-B neighbors but also on the number of theirneighbors’ non-B neighbors, and so on. Therefore,
the heuristics are at best ﬁrst-order approximations to
the Fisher-Wright model.
Table 2 shows the expected volatility at a positively selected site as predicted by our two heuristics,
compared to simulations of the full Fisher-Wrightmodel. According to either heuristic and according to
Fisher-Wright simulations over a range of h values,
the expected volatility at a site after a positively
selected sweep is always greater than the expectedvolatility of a neutral or negatively selected site
(Table 1). As discussed above, following a positively
selected sweep, the elevated volatility at the site willsubsequently decay on a time scale of order of 1/u
generations.
There is an important distinction between thevolatility signature of positive selection versus that ofnegative selection. The depressed volatility at a site
under negative selection is caused by a mutationselection-drift balance. When the eﬀective population
size is small, a large number of sites are required to
distinguish negative selection from neutrality. By
contrast, the volatility signature of positive selection
is not an equilibrium property, and it is not as sen-sitive to population size. As shown in Table 2, thevolatility of a site after a positively selected sweep
depends only weakly on h.
It is worth noting that the elevated volatility after apositively selected sweep for serine will decay evenmore slowly than for other amino acids, because thehighly volatile codons AGC and AGT are not connected by synonymous mutations to other serine codons. A mutator period subsequent to a selective
sweep will accelerate the decay toward lower volatility.
Inferring positive selection
As we have seen, a gene that contains many sites
under positive selection will exhibit a greater volatility (controlling for its amino acid composition) than agene under mostly neutral or, especially, negative
selection. How many positively selected sites are required in order to detect a reliable signal? In the caseof arginine, for example, at h = 1 the expectedvolatility of a site that has recently experienced
a positively selected sweep is approximately
0.6684 ± 0.0819, whereas a neutral arginine site has
expected volatility 0.6528 ± 0.0985, and an arginine
site under negative selection has expected volatility
0.6404 ± 0.0693 (obtained by Fisher-Wright simulation at h = 1; cf. Table 1). Therefore, the volatilityof 56 positively selected arginine sites will be signiﬁcantly greater than that of 56 negatively selected sites,at the 95% conﬁdence level. Similarly, the volatility of
260 arginine sites under positive selection will besigniﬁcantly greater than that of 260 neutral sites. For
h = 0.1, 190 sites are required to reliably distinguishpositive selection from negative selection; and 280
sites are required to distinguish positive selection
from neutrality. Similar results hold for leucine andserine; glycine is less informative.

Discussion
In the preceding sections we have analyzed the consequences of selection on proteins for synonymouscodon usage, and for volatility in particular. For allparameter values, we have demonstrated that negative
selection at the amino acid level depresses volatility
and positive selection increases volatility, comparedto neutrality. We have elucidated the conditions underwhich we can expect to draw inferences about selection on proteins from synonymous codon usage. Allof these results have been derived under the FisherWright model of population genetics—the samemodel that has been used to quantify the power ofmost other methods to detect selection (McDonald
and Kreitman 1991; Sawyer and Hartl 1992; Goldman
and Yang 1994; Ewens 2004).

648

As others have noted (Dagan and Graur 2004;
Hahn et al. 2005), the idea of estimating selection ona protein using a single nucleotide sequence challenges the essentiality of the comparative method in
evolutionary research—a paradigm that has dominated the ﬁeld for four decades (Zuckerkandl andPauling 1965). Although we have demonstrated thereliability of this new approach in certain parameter
regimes, there is the potential for confounding factors
ignored by the Fisher-Wright model. In the discussion below we compare methods of detecting selection on proteins, and we discuss several important,practical limitations of the volatility method.
Codon Volatility Versus Comparative Sequence
Analysis
Selection pressures on proteins are usually estimatedby comparing homologous nucleotide sequences(Zuckerkandl and Pauling 1965). Orthologous genes
are identiﬁed in diﬀerent organisms and sequenced;
their sequences are then aligned, and the changes thathave accumulated since divergence are used to infer
the selection pressures that have been acting (Goldman and Yang 1994). When available, sequencevariation sampled from individuals within a species
can be compared with variation across species to
produce an elegant test for adaptive evolution at a
locus (McDonald and Kreitman 1991; Sawyer andHartl 1992). In addition, there is a variety of statistical tests designed to detect a departure from neutrality in the site frequency spectrum sampled within
a single species (see Kreitman 2000 and references
therein). In many cases, the complete distribution ofthese statistics under the neutral null model are difﬁcult to derive, but they have been studied throughcomputer simulation (Simonsen et al. 1995).
Techniques for estimating selective constraints via
comparison of divergent sequences are typically applied to one or several genes at a time. When extensive intra- or inter-speciﬁc sequence data are available
at a locus of interest, such techniques have provenenormously useful for measuring selection, and it isunlikely that they will be signiﬁcantly improved byincorporating information about synonymous codon
usage. But the accurate estimation of selective constraints requires a large number (approximately six ormore [Anisimova et al. 2001]) of orthologous sequences for each gene of interest. At the genome-wide
scale, comparative data (i.e., orthologous gene sequences) will not be available for all genes, andmethods to estimate selective constraints based onsequence comparison will often be inapplicable.
Furthermore, the genes under positive selection areoften of particular interest, but such genes are evenless likely to have identiﬁable orthologs in related
species due to their rapid sequence divergence (Plot-kin et al. 2004). Even in the lineage of the Saccharomyces genus, which is currently the best-casescenario for comparative genomics, the genomes offour species have been fully sequenced and only twothirds of the genes in S. cerevisiae have unambiguously identiﬁable orthologs in related species (Kellis
et al. 2003). Unlike comparative techniques, theanalysis of synonymous codon usage oﬀers a computational tool to screen for selection pressures on all
genes in a sequenced genome. Genome-wide screens
based on analyzing synonymous codon usage may
prove useful in identifying important classes of genes
under strong selection, such as the antigens ofpathogens (Plotkin et al. 2004). Nevertheless, at thescale of individual genes, screens for selection basedon codon usage will probably not have enough statistical power except for some micro-organisms,
which experience very large eﬀective population sizes
or periods of hyper mutation.
Unlike most statistics that test for a departurefrom neutrality based on comparative data, estimatesof selection based on ‘‘volatility p-values’’ (Plotkin
et al. 2004) are not estimators in a rigorous statisticalsense—i.e., statistics whose sampling properties canbe derived from a null model, and which can be usedin likelihood ratio tests of a null hypothesis (Yanget al. 2000; Clark et al. 2003). Given the expected
relative frequencies of codons that we have derived
for each of the three regimes (neutral, negative, andpositive selection; Tables 1 and 2), it would be possible to develop maximum-likelihood methods that
estimate the number of sites of a gene in each regime.
This approach will be complicated, however, by theneed to control for other sources of codon bias (seebelow). It is critical to note that volatility p-values,
which control for each geneÕs amino acid content,measure only the relative selection pressures on theproteins in a genome, and therefore cannot distinguish between positive or relaxed negative selection.
This point has been emphasized previously (Plotkin
et al. 2004, 2005) but misunderstood by several
authors (see Supplementary Information).
Aside from the diﬀerent situations in which they
are applicable, estimates of selection based on codon
volatility diﬀer in a fundamental way from mostestimates based on sequence comparison. Homologous sequence comparison between species is often
used to assess, by either maximum likelihood(Goldman and Yang 1994) or maximum parsimony(Li 1993), the rates of synonymous and non-synonymous substitutions in a coding sequence. The ratioof these rates, dN/dS, is used as a measure of theselective constraints that have been acting on a protein since the divergence of the species being compared. An alternative approach, based on a Poisson
Random Field (PRF) model of polymorphism frequencies, uses the site frequency spectrum at a locus

649

sampled from a population to deduce the average
selective pressure for or against amino acid changes
in a gene (Sawyer and Hartl 1992). PRF models canalso be used to construct likelihood ratio tests ofdeparture from neutrality (Bustamante et al. 2001).
Like most comparative methods, however, both ofthese models typically assume that all the sites withina gene experience the same selective pressure against
amino acid substitutions (but see the site-by-site
likelihood tests of Yang et al. [2000]). Using the PRF
method, for example, authors have estimated a verysmall ‘‘average’’ selection pressure against amino acidchanges in E. coli genes: r  10)8 (Hartl and Sawyer
1994). This value does not represent the arithmetic
average of the true r values across sites, but rather
the best-ﬁt constant value of r that would make the
PRF model consistent with observed sequence variation at polymorphic sites.
When evolutionary rates are estimated at individual residues (Yang 2000; Yang et al. 2000), however,
we ﬁnd great variation across sites. Moreover, direct
experimental measurements of the ﬁtness consequences of amino acid substitutions in micro-organisms reveals huge variation in selection pressuresacross the residues of an individual protein: a substantial proportion of substitutions is lethal and asubstantial proportion has undetectable eﬀect(Wertman et al. 1992; Wlocha et al. 2001; Zyle and
De Visser 2001; Sanjuan et al. 2004). Therefore, it is
not entirely clear how best to interpret the value ofr  10)8 estimated for E. coli genes using the PRF
model, which assumes constant pressure at eachresidue.
Compared to dN/dS or r estimated by the PRF
model, codon volatility quantiﬁes selection pressuresin a very diﬀerent, coarser manner. As discussed
above, volatility essentially measures the number ofsites in a gene that experience negative (r  u) versus
neutral (r  u) versus positive selection. Given thatamino acid changes at many sites in a protein sequence are lethal while changes at other sites have noeﬀect, it is reasonable and meaningful to estimate thenumber of sites in the selected versus neutral regimes.
But volatility is not sensitive to variation in selective
pressures within either of these regimes. Hence, thevolatility measure is in some ways a coarser description of selective pressure than PRF or dN/dS. One
should not necessarily expect that volatility will correlate very strongly with dN/dS or PRF estimates,
because the latter measures represent some sort ofaverage r over the entire gene, and are thus presumably sensitive to the full range of variation in r. A
measure of selective constraint based on codon volatility is therefore diﬀerent from and complementary
to measures based on dN/dS or the PRF model.
In our analysis of synonymous signatures of
selection on proteins, we have treated sites as inde-

pendent and freely recombining. The same independence assumption is made under the PRF model, andunder most incarnations of dN/dS. In reality, sites
within a gene are often tightly linked. Linkage increases the variance in estimates of selection pressuresinferred by the PRF model (Akashi and Schaeﬀer1997), and it is likely to have similar eﬀects on
inferences made from synonymous codon usage. The
precise eﬀects of linkage on estimates of selection
remains a topic for future research. We also note that
temporal variations in the selection pressure at a site
may aﬀect synonymous codon usage (Myers et al.
2005).
It is important to note that the most common
model used to estimate dN/dS from orthologous
nucleotide sequences does not itself reﬂect the true
relationship between selection and volatility. dN/dS isoften estimated by ﬁtting maximum-likelihood
parameters to a simpliﬁed Markov-chain model ofsequence evolution that ignores population variability (Goldman and Yang 1994). Models that ignore
population variability are perfectly reasonableapproximations when comparing the sequences ofdivergent lineages; but such models fail to detect theeﬀect of amino acid selection on synonymous codon
usage. Such models consider only a single sequencethat is assumed to represent the dominant genotype
in the population at any time. Mutation and selection
are modeled simultaneously by adjusting the transition rates between codon states in the sequence(Goldman and Yang 1994). As a result, in equilibrium, the number of transitions into a state per unit
time must equal the number of transitions out of thatstate; and so equilibrium synonymous codon usage
does not depend the strength of selection in these
simpliﬁed models (Goldman and Yang 1994). In fact,
such models typically require as parameters the speciﬁcation of the equilibrium codon usage (Goldman
and Yang 1994), and so they clearly cannot be used to
predict equilibrium codon usage. Therefore, simulations of sequence evolution based on these simpliﬁed
models (such as the non-frequency-dependent simulations of Zhang [2004] or those of Nielsen and Hubisz [2005] or Dagan and Graur [2004]) will fail todetect the relationship between selection and codon
usage. By contrast, more realistic simulations that
account for population variability (such as the frequency-dependent simulations of Zhang [2004] andthe Fisher-Wright simulations in this work) will
properly reﬂect the relationship between selection andcodon usage in a replicating population.

Other Sources of Codon Bias
Our analysis of volatility, much like the theoreticalunderpinning of the PRF method and dN/dS, is based

650on the Fisher-Wright model of population genetics,
and it speciﬁcally ignores sources of direct selection onsynonymous codons. In this section, however, we will
discuss other processes that induce codon bias andthat may conﬂate estimates of selection on proteins.
Although it came as a surprise to early theorists(King and Jukes 1969), it is now clear that several
processes induce unequal usage of synonymous codons. In micro-organisms, some sources of bias, suchas biased nucleotide content, apply roughly equally to
all genes in a genome. To the extent that othersources of codon bias apply equally across a genome,
there is a straightforward randomization procedure
to control for these biases when comparing the volatilities of genes in a genome (Plotkin et al. 2004). If,
however, other sources of codon bias diﬀer from geneto gene within a genome, they may (if not properly
controlled for) introduce errors into estimates ofrelative selection pressures on proteins inferred fromcodon volatility (Plotkin et al. 2004). By the same
token, selection on synonymous codons—particularly
selection that varies from gene to gene—will also
conﬂate estimates of selection obtained by classicaltechniques such as dN/dS (Sharp and Li 1987; Hirsh
et al. 2005; Chamary et al. 2006).
One process that may vary across the genome isthe transition/transversion mutation bias. Results on
S. cerevisiae, whose genome exhibits variation in thetr/tv bias, suggest that this source of variable codon
bias will not distort estimates of selection based onvolatility: whether or not one accounts for the variation in the tr/tv bias across the genome of S. cerevisiae one obtains virtually the same rankings ofvolatility p-values (r > 0.99). The same result holds
for the E. coli genome as well.
Aside from mutational biases, there are othersources of codon bias that vary from gene to gene insome organisms. In the yeast S. cerevisiae, researchers
have observed that synonymous codon usage, measured by the Codon Adaptation Index (CAI) (Sharpand Li 1987), is correlated with a geneÕs expression
level in laboratory conditions (Coghlan and Wolfe
2000). This correlation is thought to be caused by
selection for translational eﬃciency and/or accuracy:a codon corresponding to a more abundant tRNA is
expected to be translated more quickly (due to thehigher probability per unit time that the appropriatetRNA will ‘‘ﬁnd’’ the codon) and more accurately(since the correct tRNA will likely have the greatest
chance of pairing if it is the most abundant). Considering this alternative source of biased codon usage,
two questions should be asked: Do other sources ofcodon bias distort estimates of selection? And how
can we control for these confounding factors? We do
not have a truly satisfactory answer for either of thesequestions, but the discussion below may shed somelight on the issues involved.

The degree to which other sources of codon bias
distort estimates of selection on proteins will dependon the organism being studied. In some species, suchas humans, codon frequencies exhibit a much weakercorrespondence with tRNA abundances than in otherspecies. In a species with a strong correspondence
between codon usage and tRNA abundances, theextent to which variation in this source of codon bias
aﬀects volatility will depend on whether volatile codons are (un)preferred: if there is no correlation between volatility and tRNA abundances, then theother sources of codon bias will introduce only random error into volatility estimates, making the conclusions based on volatility less powerful but still
reliable. If instead the preferred codons tend to have
either high or low volatility, then this eﬀect will
introduce systematic errors into volatility estimates of
selection (Stoletzki et al. 2005). In the latter case, inorder to quantify how much codon usage bias iscaused by volatility as opposed to other factors, one
would require a method to quantify for individual
genes the amount of codon bias due to these otherfactors. Unfortunately we do not have near the necessary level of predictive power for other sources ofcodon bias in any organism. Although expression
levels are somewhat predictive of codon bias,
expression levels do not explain most of the variationin codon bias in any genome studied thus far (Akashi
2001; Coghlan and Wolfe 2000). Until the various
sources of biased codon usage can be reliably disentangled, we cannot reliably quantify the eﬀects ofthese biases on estimates of selection obtained via
volatility. Nor can we quantify the eﬀects of thesebiases on dN/dS, which is also conﬂated by translational selection on synonymous sites.
The second question, how to control for othersources of codon bias, is also diﬃcult to answer atpresent. As discussed above, a truly satisfactory
method to control for other sources of bias would
require us to quantify other sources of codon bias in a
predictive manner for each gene. While this degree ofprecision is not currently possible, one approach is to
assume that the codon bias measured by CAI is entirely independent of volatility, and to then control
for CAI using partial correlations. For several reasons, we expect this approach to be conservative, as
we illustrate using the yeast S. cerevisiae (we use thisspecies as an example because it shows a strong
preference for codons that match abundant tRNAs,
and because we have reliable dN/dS values for almost
two-thirds of its genes, calculated from multiplealignments of closely related species [Hirsh et al.
2005]). First, we note that dN/dS is itself strongly
correlated with both CAI and gene expression levels
in yeast (Pal et al. 2001), and it is likely therefore thatany measure of selective constraint consistent with
dN/dS will itself be strongly correlated with CAI and651

expression. Second, it is possible that the codon bias
measured by CAI is in part caused by volatility (i.e.,highly expressed genes tend to experience stronger
purifying selection and therefore exhibit codon usage
biased towards lower volatility), and so controllingfor CAI would be inappropriate. Despite several
biological hypotheses (Drummond et al 2005), there
is no accepted mechanistic explanation for the correlation between CAI and dN/dS in yeast (Pal et al.
2001; Akashi 2001), and so it is unclear whether
controlling for CAI is entirely appropriate.
Nevertheless, we have examined the relationshipbetween volatility and dN/dS while controlling for
CAI. We ﬁnd that even when controlling for CAI
there remains a highly signiﬁcant partial rank correlation between volatility p-values and dN/dS in yeast(p < 10)22) (Plotkin et al. 2005). We interpret thisresult as evidence that volatility is measuring selective
constraints above and beyond any signal inherent in
CAI Partial correlations are unreliable, however,when controlling for a noisy variable (Drummond etal. 2006). Although CAI is measured without noise,
CAI is a noisy proxy for translation rates. An alternative approach proposed by Drummond et al.
(2006), principal component regression, cannot beapplied in the case of only two predictor variables,
since the PCA will choose axes weighted equally byeach predictor variable. Therefore, we performed aprincipal component rank regression on three predictor variables, and we ﬁnd that volatility p-values
explain roughly the same amount of independent
variation in dN/dS as is explained by CAI or mRNA
expression levels (8.2%, 11.1%, and 9.3%, respectively). This analysis further supports our conclusionthat volatility is correlated with selective constraints
on yeast proteins, above and beyond any signal
determined by CAI or expression levels.
Finally, we note that there may be direct selection
on synonymous codons in order to evade mis-translation (Konopka 1985). Since mis-translation is most
likely to occur between a codon and an anticodon
that diﬀer by a single nucleotide, the deﬁnition ofvolatility (Eq. 1) is appropriate for measuring theselective pressure for or against mis-translation. The
strength of this type of selection on synonymous codons would not depend on h and r but, rather, on the
mis-incorporation rate of tRNA (which is far higher
than the mutation rate) and the detriment of mistranslation (which is likely far lower than that ofmost mis-sense mutations). It is diﬃcult at present tomeasure the molecular parameters of tRNA misincorporation and its ﬁtness eﬀects; so it is unclear
how much of a volatility signal arises from mistranslation avoidance versus selection on mis-sense
mutations. Whatever the strength of the mis-translation signal, however, the volatility of a gene would
still reﬂect the degree to which there is selection toconserve, or not to conserve, the (translated) protein
sequence.

Practical Limitations of Volatility
In this paper we have presented a theoretical analysis
of codon usage and selection under the Fisher-Wrightmodel. The Fisher-Wright model accounts for only
the most basic of processes that inﬂuence molecular
evolution. Nevertheless, there is great deal of empirical evidence indicating that the volatility of a gene,
controlling for its ammo acid sequence, is indeed
correlated with the selective constraint it experiences.
Aside from highly signiﬁcant correlations between
volatility p-values and dN/dS in bacterial species andyeast (Plotkin et al. 2004), volatility also reﬂects a
range of other features known to correlate withselection on proteins. In S. cerevisiae, for example,
volatility p-values are strongly correlated with theessentiality of genes, the number of their proteinprotein interactions, and the degree to which they arepreserved throughout the eukaryotic kingdom (Plotkin et al. 2006). Furthermore, volatility is signiﬁcantly
elevated among the known antigens and surface proteins (which experience positive selection) in thepathogens Mycobacterium tuberculosis, Plasmodium
falciparum, and inﬂuenza A virus (Plotkin and Dushoﬀ 2003; Plotkin et al. 2004). In fact, elevated volatility has helped researchers to identify new families ofantigenic proteins expressed on the surface of P. falciparum and infected erythrocytes (Winter et al. 2005).
In addition, volatility is signiﬁcantly depressed in thegenes essential for growth of M. tuberculosis, as wellas in the genes conserved between related Mycobacterium species (Plotkin et al. 2004). Therefore, despite
potential confounding sources of codon bias thatcannot at present be controlled for with appropriate
accuracy, in practice volatility-based methods produce estimates of selection pressures that are consistent with our understanding of protein evolution overa diverse range of micro-organisms.
Several authors have suggested that the empiricalresults above are artifacts caused by the length oramino acid composition of rapidly evolving proteins(Nielsen and Hubisz 2005; Dagan and Graur 2004;
Friedman and Hughes 2005; Sharp 2005; Stoletzki
et al. 2005). These suggestions are patently incorrect.
The empirical results summarized above are all basedon ‘‘volatility p-values,’’ which control exactly foreach geneÕs amino acid sequence (Plotkin et al. 2004).
A geneÕs amino acid content or length cannot possibly bias its volatility p-value toward zero or one.
Related claims that results using volatility are artifacts of horizontal gene transfers (Sharp 2005), although potentially valid, are inconsistent with
empirical evidence (see Supplementary Information).

652

Nevertheless, there remain several factors that will
likely limit the applicability of volatility in practice.
Chief among these concerns is selection on synonymous sites for translational speed or accuracy, as
discussed above. Despite the fact that volatility pvalues are signiﬁcantly correlated with dN/dS evenwhen controlling for expression levels and CAI inyeast, translational selection certainly has the potential to conﬂate synonymous signatures of selection onproteins. This concern is all the more problematicconsidering that translational selection is strongest inorganisms that experience large eﬀective population
sizes, which are required for volatility to detect negative selection.
In addition, our analysis has clariﬁed and quantiﬁed a second, practical limitation of volatility:
namely, tests for selection based on synonymous signals are not very powerful, requiring many sites underselection. As a result, except for rare circumstances,
volatility cannot be expected to reliably detect positive
selection on individual genes. Perhaps it is for thisreason that most of the strong empirical results using
volatility have been observed for groups of genes,
such as families of antigens (Plotkin et al. 2004).
Our analysis has also clariﬁed and quantiﬁed athird, practical limitation of volatility: the need forlarge eﬀective population sizes or mutation rates.
Volatility has virtually no power to detect negative
selection when h per site is much smaller than unity.
Estimates of h for microbes vary widely, and periodsof hyper mutation will increase the relevant value of heven when neutral site heterozygosity is small. But as
we have shown, synonymous codon usage will not
reﬂect selection on proteins in organisms, such as
macroscopic species and probably some microbes as
well, that consistently experience small population
sizes and mutation rates.
Acknowledgments.
We thank Daniel Fisher, Andrew Murray,and Michael Turelli for their input during the preparation of themanuscript. We also thank an anonymous referee for substantial
conceptual input. J.B.P. acknowledges support from the Harvard
Society of Fellows, the Milton Fund, and the Burroughs Wellcome
Fund. M.M.D. acknowledges support from a Merck Award for
Genome-Related Research.

References
Akashi H (2001) Gene expression and molecular evolution. Curr
Opin Genet Dev 11:660–666
Akashi H, Schaeﬀer SW (1997) Natural selection and the frequency
distribution of ‘‘silent’’ DNA polymorphism in Drosophila.
Genetics 146:295–307
Anisimova M, Bielawski JP, Yang Z (2001) The accuracy andpower of likelihood ratio tests to detect positive selection atamino acid sites. Mol Biol Evol 18:1585–1592

Berg O (1996) Selection intensity for codon bias and the eﬀective
population size of Escherichia coli. Genetics 142:1379–1382
Bjedov I, Tenaillon O, Gerard B, et al. (2003) Stress-induced
mutagenesis in bacteria. Science 300:1404–1409
Bulmer M (1991) The selection-mutation-drift theory of synonymous codon usage. Genetics 129:897–907
Bustamante CD, Wakely J, Sawyer S, Hartl DL (2001) Directional
selection and the site-frequency spectrum. Genetics 159:1779–
1788
Chamary JV, Parmley JL, Hurst LD (2006) Hearing silence: nonneutral evolution at synonymous sites in mammals. Nature Rev
Genet 7:98–108
Charlesworth B, Morgan MT, Charlesworth D (1993) The eﬀect ofdeleterious mutations on neutral molecular variation. Genetics
134:1289–1303
Clark A, Glanowski S, Nielsen R, Thomas P, Kejariwal A, MA
MT, Tanenbaum D, Civello D, Lu F, B BM, Ferriera S, Wang
G, Zheng X, White T, Sninsky J, Adams M, Cargill M (2003)
Inferring nonneutral evolution from human-chimp-mouseorthologous gene trios. Science 302:1960–1963
Coghlan A, Wolfe KH (2000) Relationship of codon bias tomRNA concentration and protein length in Saccharomyces
cerevisiae. Yeast 16:1131–1145
Crow JF, Kimura M (1970) An introduction to population genetics
theory. Burgess, Minneapolis
Dagan T, Graur D (2004) The comparative method rules! codon
volatility cannot detect positive Darwinian selection using asingle genome sequence. Mol Biol Evol 22:496–500
Debry R, Marzluﬀ WF (1994) Selection on silent sites in the rodent
H3 historic gene family. Genetics 138:191–202
Denamur E, Lecointre G, Darlu P, OTenaillon CA, Sayada C,
Sunjevaric I (2000) Evolutionary implications of the frequenthorizontal transfer of mismatch repair genes. Cell 103:711–721
Drummond DA, Bloom JD, Adami C, Wilke CO, Arnold FH(2005) Why highly expressed proteins evolve slowly. Proc Natl
Acad Sci USA 102:14338–14343
Drummond DA, Raval A, Wilke CO (2006) A single determinant
dominates the rate of yeast protein evolution. Mol Biol Evol
23:327–337
Ewens W (2004) Mathematical populations genetics I. SpringerVerlag, New York
Fleischmann RD, Alland D, Eisen JA, Carpenter L, White O,
Peterson J, DeBoy R, Dodson R, Gwinn M, Haft D, Hickey E,
Kolonay JF, Nelson WC, Umayam LA, Ermolaeva M, Salzberg SL, Delcher A, Utterback T, Weidman J, Khouri H, Gill J,
Mikula A, Bishai W, Jacobs WR, Venter JC, Fraser CM (2002)
Whole-genome comparison of Mycobacterium tuberculosis
clinical and laboratory strains. J Bacteriol 184:5479–5490
Friedman R, Hughes AL (2005) Codon volatility as an indicator of
positive selection: Data from eukaryotic genome comparisons.
Mol Biol Evol 22:542–543
Gibbons RJ, Kapsimalis B (1967) Estimates of the overall rate of
growth of the intenstinal microﬂora for hamsters, Guinea pigs,
and mice. J Bacteriol 93:510–512
Gillespie J (2001) Is the population size of a species relevant to its
evolution? Evolution 55:2161–2169
Giraud A, Radman M, Matic I, Taddei F (2001) The rise and fallof mutator bacteria. Curr Opin Microbiol 4:582–585
Golding GB, Strobeck C (1982) Expected frequencies of codon use
as a function of mutation rates and codon ﬁtnesses. J Mol Evol
18:379–386
Goldman N, Yang Z (1994) Codon-based model of nucleotide
substitution for protein-coding DNA sequences. Mol Biol Evol
11:725–736
Hahn MW, Mezey JG, Begun DJ, Gillespie JH, Kern AD, LangleyCH, Moyle LC (2005) Codon bias and selection on single
genomes. Nature 433:E1

653
Hartl DL, Sawyer SA (1994) Selection intensity for codon bias.
Genetics 138:227–234
Higgs P (1994) Error thresholds and stationary mutant distributions in multi-locus diploid genetics models. Genet Res Cambr
63:63–78
Hirsh AE, Fraser HB, Wall DP (2005) Adjusting for selection onsynonymous sites in estimates of evolutionary distance. Mol
Biol Evol 22:174–177
Ikemura T (1981) Correlation between the abundance of Escherichia coli transfer-RNAs and the occurrence of the respectivecodons in its protein. J Mol Biol 146:1–21
Kellis M, Patterson N, Endrizzi M, Birren B, Lander E (2003)
Sequencing and comparison of yeast species to identify genes
and regulatory elements. Nature 423:241–254
Kimura M, Crow JF (1964) The number of alleles that can be
maintained in a ﬁnite population. Genetics 49:725–738
King JL, Jukes TH (1969) Non-Darwinian evolution. Science
164:788
Konopka AJ (1985) Theory of degenerate coding and informationalparameters of protein coding genes. Biochimie 67:455–468
Kreitman M (2000) Methods to detect selection in populations withapplications to the human. Annu Rev Genomics. Hum Genet
1:539–559
LeClerc J, Li B, Payne WL, Cebula TA (1996) High mutation
frequencies among Escherichia coli and Salmonella pathogens.
Science 274:1208–1211
Li WH (1993) Unbiased estimation of the rates of synonymous andnonsynonymous substitution. J Mol Evol 36:96–99
Lynch M, Conery JS (2003) The origins of genome complexity.
Science 302:1401–4404
Maynard Smith J, Haigh J (1974) The hitch-hiking eﬀect of afavorable gene. Genet Res Cambr 23:23–25
McDonald JH, Kreitman M (1991) Adaptive protein evolution at
the ADH locus in Drosophila. Nature 351:652–654
Miyata T, Miyazawa S, Yashunaga T (1979) Two types of aminoacid substitutions in protein evolution. J Mol Evol 12:219–236
Myers LA, Ancel FD, Lachmann M (2005) Evolution of genetic
potential. PloS Comput Biol 1:236–243
Nagylaki T (1992) Introduction to theoretical population genetics.
Springer, Berlin
Nielsen R, Hubisz M (2005) Detecting selection needs comparative
data. Nature 433:E6
Notley-McRobb L, Seeto S, Ferenci T (2001) Enrichment andelimination of mutY mutators in Escherichia coli populations.
Genetics 162:1955–1062
Ochman H, Elwyn S, Moran NA (1999) Calibrating bacterial
evolution. Proc Natl Acad Sci USA 96:12638–12643
Oliver A, Canton R, Campo P, Baquero F, Blazquez J (2000) Highfrequency of hypermutable Pseudomonas aeruginosa in cysticﬁbrosis lung infection. Science 288:1251–1254
Pal C, Papp B, Hurst LD (2001) Highly expressed genes in yeast
evolve slowly. Genetics 158:927–931
Plotkin JB, Dushoﬀ J (2003) Codon bias and frequency-dependent
selection on the hemagglutinin epitopes of Inﬂuenza A virus.
Proc Natl Acad Sci USA 100:7152–7157
Plotkin JB, Dushoﬀ J, Fraser HB (2004) Detecting selection using asingle genome sequence of M. tuberculosis and P falciparum.
Nature 248:942–946
Plotkin JB, Dushoﬀ J, Fraser HB (2005) Codon bias and selection
on single genomes: reply. Nature 433:E7–E8
Plotkin JB, Fraser HB, Dushoﬀ J (2006) Natural selection on the
genome of Saccharomyces cerevisiae (in preparation)

Sanjuan R, Moya A, Elena S (2004) The distribution of ﬁtness
eﬀects caused by single-nucleotide substitutions in an RNA
virus. Proc Natl Acad Sci USA 101:8396–8401
Sawyer SA, Hartl DL (1992) Population genetics of polymorphismand divergence. Genetics 132:1161–1176
Sharp PM (2005) Gene ‘‘volatility’’ is most unlikely to reveal
adaptation. Mol Biol Evol 22:807–809
Sharp PM, Li WH (1987) The codon adaptation index: A measureof directional synonymous codon usage bias, and its potential
applications. Nucleic Acids Res 15:1281–1295
Simonsen KL, Churchill GA, Aquadro CF (1995) Poperties ofstatistical tests of neutrality for DNA polymorphism data.
Genetics 141:413–429
Sorensen M, Kurland C, Pedersen S (1989) Codon usage determines translation rate in Escherichia coli. J Mol Biol 207:365–
377
Stoletzki N, Welch J, Hermisson J, Eyre-Walker A (2005) A dissection of volatility in yeast. Mol Biol Evol 22:2022–2026
Tajima F (1996) The amount of DNA polymorphism maintained ina ﬁnite population when the neutral mutation rate varies amongsites. Genetics 143:1457–1465
Tang H, Wyckoﬀ GJ, Lu J, Wu C (2004) Universal evolutionary
index for amino acid changes. Mol Biol Evol 21:1548–1556
Tenaillon O, Denamur E, Matic I (2004) Evolutionary signiﬁcanceof stress induced mutagenesis in bacteria. Trends Microbiol
12:264–270
Thompson CJ, McBride JL (1973) On EigenÕs theory of the selforganization of matter and the evolution of biological macromolecules. Math Biosci 21:127–142
van Nimwegen E, Crutchﬁeld J, Huynen M (1999) Neutral evolution of mutational robustness. Proc Natl Acad Sci USA
96:9716–9820
Wertman K, Drubin D, Botstein D (1992) Systematic mutational
analysis of the yeast ACT1 gene. Genetics 132:337–350
Wilke C (2001) Adaptive evolution on neutral networks. Bull Math
Biol 63:715–130
Winter G, Kawai S, Haeggstrom M, Kaneko O, vonEuler A,
Kawazu S, Palm D, Fernandez V, Walgren M (2005) SURFIN is a polymorphic antigen expression on Plasmodium falciparum merozoites and infected erythrocytes. J Exp Med
20l:1853–1863
Wlocha DM, Szafranieca K, Bortsb RH, Korona R (2001) Direct
estimate of the mutation rate and the distribution of ﬁtness
eﬀects in the yeast Saccharomyces cerevisiae. Genetics 159:441–
452
Wright S (1931) Evolution in Mendelian populations. Genetics
16:97–159
Yampolsky LY, Stoltzfus A (2004) The exchangability of aminoacids in proteins. Genetics 170:1459–1472
Yang Z (2000) Maximum likelihood estimation on large phylogenies and analysis of adaptive evolution in human inﬂuenza virus
A. J Mol Evol 51:423–432
Yang Z, Nielsen R, Goldman N, Pedersen AMK (2000) Codonsubstitution models for heterogeneous selection pressure atamino acid sites. Genetics 155:431–449
Zhang J (2004) On the evolution of codon volatility. Genetics
16S:495–501
Zuckerkandl E, Pauling L (1965) Molecules as documents of evolutionary history. J Theor Biol 8:357–366
Zyle C, DeVisser J (2001) Estimates of the rate and distribution ofﬁtness eﬀects of spontaneous mutation in Saccharomyces cerevisiae. Genetics 157:53–61