Proc. R. Soc. B (2008) 275, 2455–2464
doi:10.1098/rspb.2008.0521
Published online 22 July 2008

Directionality in the evolution of influenza
A haemagglutinin
Sergey Kryazhimskiy1,*, Georgii A. Bazykin2, Joshua Plotkin1
and Jonathan Dushoff 3
1
Department of Biology, University of Pennsylvania, Philadelphia, PA 19104, USA
Institute for Information Transmission Problems (Kharkevich Institute) of the Russian Academy of Sciences,
Moscow 127994, Russia
3
Department of Biology, McMaster University, Hamilton, Ontario, Canada L8S 4K1

2

The evolution of haemagglutinin (HA), an important influenza virus antigen, has been the subject ofintensive research for more than two decades. Many characteristics of HA’s sequence evolution are
captured by standard Markov chain substitution models. Such models assign equal fitness to all accessibleamino acids at a site. We show, however, that such models strongly underestimate the number ofhomoplastic amino acid substitutions during the course of HA’s evolution, i.e. substitutions that
repeatedly give rise to the same amino acid at a site. We develop statistics to detect individual homoplastic
events and find that they preferentially occur at positively selected epitopic sites. Our results suggest thatthe evolution of the influenza A HA, including evolution by positive selection, is strongly affected by thelong-term site-specific preferences for individual amino acids.
Keywords: directional selection; dN/dS; haemagglutinin; homoplasy; influenza A

1. INTRODUCTION
Influenza viruses offer an extraordinary opportunity for
improving our understanding of molecular evolution.
Several hundred complete genomes of the influenza A
virus have been sequenced, as well as several thousand
variants of its primary surface antigen, haemagglutinin(HA). Over the past four decades, roughly 20% of sites,
concentrated within the epitopic regions of the HA1
domain of HA, have undergone amino acid substitutions,representing the equivalent of millions of years of evolutionary change in a typical mammalian protein (Carroll 2003).
These substitutions were driven primarily by selection toevade the antibody response in the host population(Fitch et al. 1991; Nelson & Holmes 2007).
Aside from recent studies of recombination (Lindstrom
et al. 2004; Holmes et al. 2005), most research on theevolution of HA has focused on identifying the sites thatexperience positive selection for amino acid substitutions(Fitch et al. 1997; Bush et al. 1999a,b; Yang 2000; Plotkin &Dushoff 2003; Suzuki 2006; Wolf et al. 2006; Shih et al.
2007) among the majority of sites that evolve undernegative selection. One of the key tools used for detecting
genes or sites under positive selection is the concept of thedN/dS ratio, the ratio of the rates of non-synonymous(amino acid altering) and synonymous (amino acidpreserving) substitutions along a phylogenetic tree. Assuming that the rate of synonymous substitutions is a good
approximation of the neutral standard, a dN/dS ratio* Author and address for correspondence: 204K Lynch Laboratory,
433 S University Avenue, Philadelphia, PA 19104, USA(skr@sas.upenn.edu).
Electronic supplementary material is available at http://dx.doi.org/10.
1098/rspb.2008.0521 or via http://journals.royalsociety.org.
Received 16 April 2008
Accepted 23 June 2008

exceeding unity is an indication of positive Darwinian
selection (Ina & Gojobori 1994; Yang & Bielawski 2000).
This idea stems from a prediction of the neutral theory thatthe rate of amino acid substitutions in a gene must be
smaller than or equal to the rate of synonymous
substitutions (Kimura 1977, 1983). Later, this idea wascarried over to subregions of genes and even individualsites. In several recent studies, a dN/dS analysis was applied
to HA sequences, and sites with dN/dS ratios significantly
greater than unity were identified (Ina & Gojobori 1994;
Bush et al. 1999a; Yang 2000; Suzuki 2006; Wolf et al.
2006), suggesting that positive selection plays an importantrole in the evolution of influenza. Although this researchsheds light on HA evolution, and although it may help tocalibrate appropriate vaccines (Bush et al. 1999b; Plotkin
et al. 2002), a simple catalogue of sites experiencing positive
selection fails to address the full spectrum of possibilities forthe action of natural selection on HA. In particular, an
important question—which amino acids were fixed in HA
due to natural selection as opposed to random drift—hasremained largely untouched in the literature (but see Wolfet al. 2006; Shih et al. 2007). In this work, we attempt totake an initial step in this direction.
It is instructive to re-examine the substitutions models
in which the dN/dS concept is developed. Such models,
introduced by Goldman & Yang (1994) and by Muse &Gaut (1994), are based on the theory of finite-state,continuous-time Markov chains. This general approachhas proved extremely fruitful, and various modifications of
the original models have been suggested ( Nielsen & Yang
1998; Yang & Nielsen 2000, 2002, 2008; Yang et al. 2000;
Forsberg & Christiansen 2003; Guindon et al. 2004;
Mayrose et al. 2007). We refer to this whole family ofmodels as ‘codon-based Markov chain models’ or simply‘Markov chain models’. The key feature of these models is

2455This journal is q 2008 The Royal Society

2456 S. Kryazhimskiy et al.

Directionality in evolution of HA

that codons rather than nucleotides are treated as theevolving unit in the coding sequences. Each position in the
evolving sequence can assume one of 61 states (stopcodons are usually excluded), and transitions between
states occur with certain rates. The synonymous substitution rates equal the corresponding mutation rates,
and the non-synonymous substitution rates equal themutation rates scaled by the dN/dS ratio (usually labelledas u), so that sites with dN/dS ratios smaller (greater) thanunity evolve, at the amino acid level, slower (faster) than
expected from neutrality.
Sites under positive selection are typically identified byreconstructing the phylogenetic history among sampled
sequences, and inferring the evolutionary parameters of a
Markov chain model, either within the maximumlikelihood or Bayesian setting. A statistical test is thenapplied to determine whether the dN/dS ratios of any sites
significantly exceed unity.
However, what exactly does it mean that a site evolves
according to a Markov chain model with us1?
Standard Markov chain substitution models assume
independence of sites and assign an u value to each site. In
such models, substitution rates towards different aminoacids are entirely determined by the mutation rates, andno preference is given to one substitution over another.
This simplification ignores an important aspect of real
protein evolution, the fact that different amino acids are
preferred at different sites. More generally, one can picture
a site as evolving on a fitness landscape whereby eachamino acid at the site provides the organism with certain
fitness, given a genetic background. This marginal fitness
is what we call the ‘fitness of amino acid X at site n’. The
fitness landscape of a site is not static: the fitness of anamino acid can change depending on the environment,
genetic background, etc.
At a site, substitutions towards amino acids with higher
fitness values will occur faster than lower fitness values.
One straightforward consequence of this is the well-known
fact that substitutions between biochemically similaramino acids occur more readily than dissimilar aminoacids (Miyata & Yasunaga 1980). Although many widelyused Markov chain models ignore biochemical properties,
models have been developed that account for them(Goldman & Yang 1994; Yang et al. 1998; Sainudiin
et al. 2005; Wong et al. 2006).
Another, more serious, consequence of the aforementioned simplification is the fact that Markov chain models
with a single u value per site describe a site evolving on afitness landscape that changes with each substitution event.
In particular, in the ‘negative selection mode’ (u!1), theamino acid currently fixed in the population provides a local
fitness maximum, whereas in the ‘continual positive
selection mode’ (uO1), it provides a local fitness
minimum. If a site is characterized by u!1 (uO1), then
any new arising amino acid at the site is assigned a negative(positive) selection coefficient (Nielsen & Yang 2003).
Even though it is conceivable that fitness landscapes ofsome proteins are dynamic, one rarely expects this to be thecase at the level of individual sites. Indeed, consider theevolution of the influenza A HA, a canonical example ofevolution on a highly dynamic fitness landscape shaped byconstantly changing immune status of the human population. Epitopic sites in HA are described by a Markov
chain model with uO1 (Bush et al. 1999b; Yang 2000).
Proc. R. Soc. B (2008)

Suppose that amino acid A was recently substituted byamino acid B at such site, for example, because amino acid
B provided a protein conformation that was less recognizable by human antibodies. According to the Markov chain
model with uO1, immediately after the substitution
occurs, amino acid B becomes deleterious and any aminoacid other than B is preferred at the site. This implies that,as soon as A is replaced by B, there is immediately pressure
to change again, to any other amino acid (including A),
even though the immune status of the population has
hardly changed. Analogously, consider a site that evolves
under negative selection and is described by a Markov
chain model with u!1. Suppose that a certain amino acidA is currently fixed in the population and provides a local
fitness maximum. Any other amino acid at this site istherefore deleterious. Nevertheless, in a finite population, A
can be substituted by some other amino acid B. In thenegative selection regime, one would expect to observe a
rapid reversion to the locally optimal A. However, in a
Markov chain model with u!1 such substitution isunlikely, and B becomes the local fitness maximum. Thisimplies that the fitness landscape has changed, and thesubstitution A/B was, in fact, adaptive.
These examples illustrate that standard Markov chain
models may describe unrealistic evolutionary trajectories
for individual sites and/or lead to wrong interpretations ofevolution under an elevated/decreased dN/dS ratio. In
reality, the fitness landscape of a site is expected to be lessdynamic than assumed by such Markov chain models. In
order to treat this problem, considerable effort has been
made in recent years to create substitution models that
attempt to use more fully the information about the fitness
landscape on which a protein evolves. This can be done atvarious levels of detail and complexity. In the mostcomplex models, substitution rates are specified basedon the three-dimensional protein structure ( Fornasari
et al. 2002; Rodrigue et al. 2005; Robinson et al. 2007). In
another type of model, sites are treated independently but
a distinct fitness value is ascribed to each amino acid ateach site (Halpern & Bruno 1998). In yet another type ofmodel, sites are grouped into classes and different fitness
landscapes are modelled in different site classes ( Thorne
et al. 1996; Koshi & Goldstein 1998; Koshi et al. 1999;
Dimmic et al. 2000; Lartillot & Philippe 2004; Sainudiin
et al. 2005; Wong et al. 2006). These models have beenused to more accurately estimate evolutionary distances
between sequences (Halpern & Bruno 1998), reconstruct
phylogenies ( Thorne et al. 1996; Koshi et al. 1999;
Lartillot & Philippe 2004) and infer the dN/dS ratios(Sainudiin et al. 2005; Wong et al. 2006; Robinson et al.
2007). Although some of these models, notably those by
Halpern & Bruno (1998) and Koshi & Goldstein (1998),
are capable of capturing evolution at a site on a static
fitness landscape, no conclusions about the fitness landscapes on which real proteins evolve have been drawn fromthese models.
Very recently, Kosakovsky Pond et al. (in press)
developed a maximum-likelihood (ML) approach fordetecting directional selection for amino acid substitutions
and applied it to the influenza A sequence data. Theirmethod, although potentially very useful, has two important limitations. First, it uses an amino acid- rather thancodon-based Markov chain model and, as such, may notproperly account for the underlying mutational biases,Directionality in evolution of HA
which may potentially lead to spuriously significant
results. Second, the method lacks power because, inorder to detect directional selection towards a certain
amino acid, it requires an increase in the overall frequencyof that amino acid in the entire gene.
Here we take an alternative approach to studying
directionality in the evolution of an amino acid site in a
protein. We use an independent site codon-based Markov
chain substitution model with site-specific dN/dS ratios as
our null model and quantify the departures from it thatstem from the violation of the assumption of a fitness
landscape that changes with every substitution event. In
reality, the changes in the fitness landscape of a site areunlikely to be precisely synchronized with the substitution
events, so that a specific amino acid or a set of specificamino acids are selectively advantageous at a site for a
prolonged period of time. If this set is small, convergent
and/or parallel evolution may be observed ( Yeager et al.
1997; Zhang & Kumar 1997). This fact can be used todevelop an analogue of the dN/dS test to detect not only
sites but also amino acids at those sites that are substituted
more frequently than expected under neutrality. A similar
approach has been suggested earlier by Chen et al.
(2004) who have detected the majority of known as wellas many previously unknown amino acid substitutions
that are likely to have arisen in HIV as a response to drug
therapy. Their approach was later formalized within the
ML framework by Seoighe et al. (2007). In this paper,
we provide a characterization of the fitness landscapeof influenza A HA. Our approach differs from that byChen et al. and Seoighe et al., in which we (i) take into
account the phylogenetic relationships between sequences
to infer the amino acid substitution opportunities and(ii) quantify the departures of the HA evolution from a
Markov chain model with site-specific dN/dS values rather
than from a purely neutral model.
Unlike the dN/dS measure that discriminates rapidly
evolving sites from slowly evolving sites, our method isaimed at detecting individual amino acids that have
evolved under directional selection at specific sites.
Directed substitutions can occur at either rapidly or
slowly evolving sites, and we do not a priori expect aconsistent relationship between a site’s dN/dS ratio andthe presence of directionally selected amino acids.

2. MATERIAL AND METHODS(a) Data
We downloaded all HA nucleotide sequences from human
influenza A virus subtype H3N2, which were available in the
NCBI database of April 2006. We excluded isolates that were
processed in chicken eggs (Bush et al. 2000). The remaining
sequences were aligned using CLUSTAL W v. 1.83 ( Thompson
et al. 1994) and coding regions for the HA1 part wereextracted. The alignment length was 987 nucleotides.
Occasional gaps were filled if more than 90% of sequences
agreed on the symbol at the gap position, otherwise thesequence with a gap was excluded from further analysis. The
resulting alignment contained 1249 sequences. The list of
sequences is available upon request.
(b) Phylogeny reconstruction
We reconstructed parsimonious phylogenetic trees for the HA
protein of subtype H3N2 using the parsimony ratchet
Proc. R. Soc. B (2008)

S. Kryazhimskiy et al.

2457

algorithm ( Nixon 1999) built on top of PAUP (Swofford
2002). We choose the maximum parsimony (MP) algorithm
in favour of ML because (i) ML algorithms, especially those
based on realistic evolutionary models, are unreasonably time
consuming for such a large number of sequences and (ii)theoretically, the MP algorithm should perform well on thisdataset because the sequences are sufficiently similar thatmultiple substitutions in a single position on a single branch
are rare. However, our results are robust with respect to themethod of phylogeny reconstruction (see the electronic
supplementary material).
The following parameters were used for the parsimony
ratchet: fraction of sites whose weight was set to zero in the‘jump step’, 0.1; number of trees kept in memory, from 1 to 2;
and number of iterations, 30 ( Nixon 1999). The minimum
tree length obtained was 3349. We found 34 islands
containing trees with this score. We randomly selected 20
trees from these islands and performed our analysis on each ofthem. The selected 20 trees were rooted by designating theisolate V01085/Aichi/1968 as an out-group. For each of thesetrees, we inferred internal nodes by the parsimony algorithm
using PAUP with the ACCTRAN option. It turned out that theroot sequence was the same for all 20 trees and was equal to
the sequence of the out-group V01085/Aichi/1968.
(c) Simulation of sequence evolution
We simulate the evolution of the HA along the phylogeny using amodified version of the codon-based Goldman–Yang (GY)model with site-specific dN/dS ratios (Goldman & Yang 1994;
Yang 2000). Our algorithm takes the following objects as input:(i) a rooted phylogenetic tree with branch lengths that correspond to the number of nucleotide differences, (ii) the sequenceat the root node, (iii) the 4!3 mutation matrix, (iv) the listof site-specific dN/dS ratios, uk, kZ1, 2, ., L, and (v) thecodon-position correction coefficients, li , iZ1, 2, 3 (see below).
Starting from the root of the tree, the algorithm recursively
generates sequences for all other nodes. A sequence at a nodeis generated by the following rule, given that the sequence of
its ancestor is already known. Consider a node A connected toits ancestral node B by a branch of length n; in other words, A
differs from B in n nucleotides. To generate the sequence in B,
we change exactly n different nucleotides of sequence in A. As
a result, the simulated tree differs from the original tree onlyin the sequences at the nodes, while the original topology andbranch lengths are preserved.
We distribute n nucleotide changes along the sequence
as follows. Suppose m (0%m!n) changes have alreadybeen distributed along the sequence and fell on nucleotide
positions from the set Km Z fk 1 ; k 2 ; .; km g. To simulate
the mC1th change, we first assign the following weights
wmC1 ðxk / yk Þ to each possible mutation of nucleotide xk
at position k (kZ1, 2, ., 3L) in the parental sequence( yk is any nucleotide other than xk):
wmC1 ðxk /yk Þ
80
is a stop codon;0
if k 2 Km or cpðkÞ
>
>
<
 0 


if A cpðkÞ Z A cpðkÞ ;Z rðxk /yk Þ
>
>
 0 


:
sA cpðkÞ ;upðkÞ lqðkÞ rðxk /yk Þ if A cpðkÞ
ð2:1Þ
where, ci is the codon at the position i in the ancestral
sequence; p(k)Zb(kK1)/3cC1 is the amino acid position
corresponding to nucleotide position k; q(k)Z(kK1) mod0
3C1 is the position of nucleotide k in the codon cp(k); cpðkÞ
is2458 S. Kryazhimskiy et al.

Directionality in evolution of HA

a codon that equals cp(k) at all positions except for q(k) where
it has nucleotide yk instead of xk; and A(c) is the amino acidencoded by codon c.
Then, the mC1th change becomes the change of
nucleotide xk with nucleotide yk at position k with probability
proportional to the weight of the corresponding change,pðxk / yk Þ Z 3L
P

wmC1 ðxk / yk Þ:P
wmC1 ðxk / yk Þ

kZ1 yk : yksxk

After the change falls at position k mC1, we set KmC1Z
Km g fkmC1 g and repeat the procedure until all n substitutions
have been distributed.
This procedure is similar to the original GY model, withthe differences as follows.
(i) GY implements a continuous-time Markov process,drawing the number of changes that occurred on abranch from a distribution whose mean equals thebranch length. As a result, the branch lengths of thesimulated tree generally differ from those of the real tree.
This difference could lead to biases in some of the kindsof comparisons between observed and simulated data
presented below. Instead, we distribute the same exact
number of changes along a branch in the simulation as
observed in the data. Therefore, the simulated branch
lengths exactly match those in the data.
(ii) We have included an additional parameter, the codonposition correction coefficient l that weights amino acidsubstitutions according to the codon position in which asubstitution occurs (see equation (2.1)). The modelwith this parameter captures the evolution of HA better
than the model in which the probability of a nonsynonymous substitution is independent of the position
in a codon (see the electronic supplementary material).
We refer to the evolutionary model with codon-position
correction coefficient as mGYCl for ‘modified GY withcodon-position corrections’, and to the model wherel1 Z l2 Z l3 Z 1 as mGY.
The code implementing the mGYCl model written in
OBJECTIVE CAML is available upon request.
(d) Statistics(i) Excess-of-substitution statistics
The excess-of-substitutions (ES) statistic is designed todetect the individual amino acid substitutions that occurredwith unusually high frequency. Unlike the homoplasy statistic(see the electronic supplementary material), it controls forvariation in mutational opportunities of different substitutions. The ES statistic is elevated by either homoplastic
or individual substitutions that occur at a rate higher than
expected on the basis of their mutational opportunities.
First, we define the opportunity r(c/Y ) for a sense codonc to change by a single nucleotide substitution to any codon
that encodes amino acid Y other than A(c),Xrðc/ YÞ Z
lfðc;c0 Þ rcod ðc/ c 0 Þ;
c 0 2MðcÞ:Aðc 0 ÞZYwhere M(c)is the set of one-mutational neighbours of codon c
excluding the stop codons; f(c,c 0 ) is the position at whichcodons c and c 0 differ; and rcod ðc/ c 0 Þ is the rate of thenucleotide mutation that transforms codon c into codon c 0 ;Proc. R. Soc. B (2008)

for example, if r(A/G) is the mutation rate between adenine
and guanine, then rcod ðAAT/ AGTÞZ rðA/ GÞ. We putr(c/Y )Z0 if A(c)ZY. Next, the opportunity for amino acidX to be substituted by amino acid Y at site k on the whole tree
is given by
X ðnÞ
Ok ðX/ YÞ Z ukl i rðck ðiÞ/ YÞ;i2Bk ðXÞ

where, Bk ðXÞ is the set of branches of the tree whose parental
node sequence encodes amino acid X at position k; l ðnÞ
i is thetotal number of non-synonymous substitutions that occurredat branch i. Then,X
Ok ð/YÞ Z
Ok ðX/ YÞX

is the opportunity of gaining amino acid Y at site k.
Finally, let nk ð/XÞ be the number of observed (orinferred) substitutions towards amino acid X at site k, andP
let nk Z X nk ð/XÞ be the total number of amino acidsubstitutions at site k. Then,O ð/XÞ
Ek ð/XÞ Z P kn
Ok ð/XÞ kX

is the number of substitutions leading to amino acid X at site
k expected under the assumptions of the mGYCl model.
We define the ES statistic for site k and amino acid X assES ðk; XÞ Z nk ð/XÞK Ek ð/XÞ:A positive ES statistic implies that the number of amino acidsubstitutions resulting in amino acid X at site k is larger than
expected from the mGYCl model.
We define the grand excess-of-substitutions (GES)
statistic as
"
#1=2
L X
X
2
SGES ZsES ðk; XÞ
:
kZ1 X

P

BecauseX sES ðk; XÞ h 0 for each site k, a value of GES
significantly different from 0 implies an excess of substitutions
into one or several amino acids and a deficit of substitutions
into other amino acids, compared to what is expected from
opportunities Ok.
To separately compute the GES statistic for epitopes(denoted by SGESjepi ) and non-epitopes (denoted by
SGESjnepi ), the outer summation is taken over the corresponding sites.
(ii) Finding selectively advantageous amino acids
We use the ES statistic described above to identify aminoacids that were selectively advantageous or disadvantageous
at individual sites. For that, we group together all potential
amino acid substitutions at a site that have the same target
amino acid and for which the opportunity on the phylogenetic
tree is non-zero (such substitutions may or may not haveactually occurred on the phylogenetic tree). For example, allpotential amino acid substitutions that resulted in asparagine
at site 23, no matter what the substituted amino acid was,would fall in the same group. We only include a group in
further analysis if it was observed for each of the 20 trees; thisis true for 2490 groups (more than 98% of all observed
groups). Each such group of potential substitutions is
associated with a corresponding value of the ES statistic. By
comparison with the distribution of the ES statistic expectedfrom the mGYCl model, we identify the groups of potentialDirectionality in evolution of HA
substitutions resulting in the same amino acid whose ES
statistic (averaged over 20 trees) is significantly elevated or
depressed. Since we perform multiple hypothesis testing, we
expect to find N$x false positives if we are looking at the x quantile of the null distribution, where N is the number oftests. Under the stringent Bonferroni correction, we treat as
significant only those target amino acids whose ES statistic
falls into the x/N-quantile of the null distribution.

3. RESULTS
Our analysis of HA sequence data is designed to quantify
aspects of the HA fitness landscape that is not captured bya simple Markov chain model of substitutions. Our basic
methodology is straightforward: we design statistics thatquantify aspects of protein evolution, such as directionality, and we assess the significance of these statistics bycomparing the observed data with a null distribution
generated by simulating a Markov chain substitution
model. The substitution model used to generate the null
distribution is a modified version of the GY codon-basedmodel, which we call the mGYCl model. Under thismodel, the probabilities of synonymous substitutions at acodon site are determined by a neutral mutation matrix(constant across the gene), and the probabilities of nonsynonymous substitutions are determined by the mutation
matrix, the site-specific u values and codon-position
correction coefficients to u (see §2). Thus, the probabilities of non-synonymous substitutions differ between
different codon sites and the three positions in thecodon, but no consideration is given to the identitiesof the target amino acids. For the sake of comparison,
we also use the mGY model in which the probabilities of
non-synonymous substitutions are equal between threecodon positions.
In order to fit the parameters of our null model, we
reconstruct 20 equally parsimonious phylogenies for 1249
sequences of the HA1 part of the HA protein ( Wilson &Cox 1990) sampled from human patients between 1968
and 2005. From each of the reconstructed trees, we infer
the sequences at internal nodes, the mutation rates, thesite-specific dN/dS ratios and the codon-position correction coefficients (see the electronic supplementarymaterial). The inferred parameter values can be found in
table S1 in the electronic supplementary material. Then,
using these parameters and the inferred root sequence, we
produce null distributions of various statistics bysimulating the mGYCl substitution model 50 times
along each tree, obtaining a total of 1000 replicate
sequences of all leaf and internal nodes.
We compared the features of real HA protein evolution
as inferred from the 20 parsimonious trees with thefeatures of the null distribution generated by the mGYClmodel using three types of statistics. First, we assessed thedegree of HA protein conservation by measuring the
average Miyata distance for an amino acid substitution, as
well as the Simpson diversity index for the amino acid
distribution at a typical site in HA (see the electronic
supplementary material). Second, we used the GES
statistic to assess the overall overabundance of substitutions
in HA that gave rise to the same amino acid at a site. Third,
we used the ES statistic to find particular amino acids that
arose at certain sites more frequently than expected underthe mGYCl model.
Proc. R. Soc. B (2008)

S. Kryazhimskiy et al.

2459(a) Model consistency
The mGYCl model captures both the synonymous
evolution and the coarse aspects of non-synonymous
evolution of the HA gene. It provides a good fit to theobserved numbers of synonymous and non-synonymous
substitutions on the tree, the numbers of synonymous andnon-synonymous substitutions that fall at different codon
positions, and the time trajectory of the nucleotide
composition at fourfold degenerate sites and second
position sites (see the electronic supplementary material).
The mGY model, which does not control for differences inrates of non-synonymous substitutions between the threecodon positions, fails to capture the distribution ofsubstitutions between the positions; in particular, it
overestimates the number of non-synonymous substitutions at the second codon position (see the electronic
supplementary material).
(b) Protein conservation
We computed the average Miyata distance of an aminoacid substitution for the mGYCl simulations and the real
HA protein. The average Miyata distance of an amino acidsubstitution in the HA protein as a whole is significantly
smaller than in the null distribution, which ignores
physicochemical constraints ( p!0.05; see the electronic
supplementary material). This result holds for nonepitopic regions of the protein ( p!0.01), but is notstatistically significant for the epitopes ( pZ0.16). Hereand below, p-values are computed for the average (over 20
parsimonious trees) value of each statistic, and all tests are
two tailed unless noted otherwise.
We also measured protein conservation using the
average Simpson index of the amino acid frequency
distribution at each site (see the electronic supplementarymaterial). The inverse of the Simpson index quantifies theeffective number of different amino acids that have beenpresent at a typical site during the course of its evolution(Simpson 1949). In other words, if the average Simpson
index is equal to 1, then most sites have seen many aminoacids in the course of evolution; while if the Simpson index
is close to 1, then most sites have been dominated by asingle amino acid. The latter situation would indicate
strong protein conservation. The mGYCl model captures
the degree of conservation of the HA protein as a whole
relatively well ( pZ0.2; see the electronic supplementarymaterial). However, separating sites into epitopes andnon-epitopes reveals an important pattern: non-epitopic
sites are strongly conserved (the Simpson index issignificantly greater than expected, p!10K3 ) while theepitopic sites are not (the Simpson index does not
significantly differ from expectation, pZ0.14). The
effective number of amino acids at a typical epitopic site
is 1.20 while at a non-epitopic site it is 1.02.
(c) Directional selection
At each amino acid site, an average of 2.0 different aminoacids arose by non-synonymous substitutions at least once
somewhere on the phylogenetic tree, compared to 2.3 insimulations of the mGYCl model (standard errors are lessthan 0.1). 334.5G0.3 (49%) of these amino acids arose on
average more than once, compared with 364.1G0.4 (48%)
in simulations. In other words, the actual pattern of
substitutions during HA evolution exhibits characteristics
of directional selection: substitutions tend to lead to a2460 S. Kryazhimskiy et al.

Directionality in evolution of HA

specific set of amino acids, producing a smaller number ofamino acids observed at each site.
In order to estimate whether the substitution rates
towards certain amino acids are significantly higher in thedata than expected under the mGYCl model, we use
the ES and GES statistics. The ES statistic is defined as thedifference between the observed number of substitutions
towards a specific amino acid at a site and the corresponding number expected if the substitution probabilities were
determined according to the mGYCl model (see §2 fordetails). This statistic is elevated if substitutions leading to
the amino acid under consideration at a given site occurred
more frequently than expected. The GES statistic is thesum of squares of the ES values for each amino acid at eachsite; it is elevated if many amino acid substitutions in thewhole gene occurred at unexpectedly high rates.
Using the GES statistic, we find that significantly more
substitutions gave rise to the same amino acids at a site in
the data than in simulations. The value of the GES statistic
is significantly larger in the data than in simulations, for
the HA protein as a whole and for both the epitopic andthe non-epitopic parts of the HA protein ( p!10K3;
figure 1). We also explored other statistics that measure
the same quantity, and the results are qualitatively thesame (see the electronic supplementary material).
Using the ES statistic, we identified the individual aminoacids that were significantly more and less likely to be thetargets for non-synonymous substitutions at a given sitethan expected from the corresponding mutational opportunities, site-specific dN/dS ratios and codon-position
corrections. The ES statistic is significantly elevated at the
0.1 per cent level (one-tailed test without the Bonferroni
correction) for 28 amino acids at 25 sites (table 1); theES statistic is significantly depressed at the 0.1 per centlevel for 19 amino acids at 13 sites (table 1). The expected
number of false positives in each case is 2.5. Therefore,
most (although probably not all) of the target aminoacids presented in table 1 have, in fact, been favoured
or disfavoured by natural selection. In each of the cases
presented in table 1 with the elevated ES statistic, the target
amino acid arose more than once on the phylogenetic tree.
Of the 25 reported sites (84%), 21 are located in theepitopic regions of the protein; 20 sites (80%) have dN/dS
values exceeding 1, indicative of positive selection as the
predominant evolutionary force (table 1).
The presented list of sites with at least one amino acidwith an elevated ES statistic partially overlaps withpreviously published lists of sites with an elevated dN/dS
ratio (Bush et al. 1999a; Yang 2000; Suzuki 2006). We do
not compare these lists explicitly because the dN/dS andthe ES statistics measure different quantities, and we haveno a priori expectation for their relationship. Even though
we have excluded from our analysis all variants that were
denoted as being cultured in chicken eggs, some of theamino acids listed in table 1 (e.g. lysine at sites 145 and156, serine at site 186) were previously identified as being
adaptations of the influenza virus to growth in eggs
(Meyer et al. 1993). This is not surprising, as it is known
that some mutations first observed in egg cultured variants
also exist in circulating human strains (Rocha et al. 1993).
Figure 2 shows the phylogenetic distribution of lysine atsite 145, which has the highest value of the ES statistic andprobably evolved under directional selection: it was notpresent in the majority of the virus population until theProc. R. Soc. B (2008)

0.2

(a)

0.1
73.61
0
20

0.2

30

40

50
SGSE

60

70

80

(b)

0.1
68.50
0
20

0.2

30

40

50
SGSE|epi

60

70

80

40

50

60

(c)

0.1
26.96
0

10

20

30
SGSE|nepi

Figure 1. GES statistic. Grey bars represent the histogram for
the GES statistic distribution obtained from the simulation
based on the mGYCl model. The corresponding best-fit
Gaussian curves are indicated in black. Dark grey triangles
show the values of the GES statistic for the real HA protein as
inferred from each of the 20 MP trees, with the mean value
shown above the triangles. Data for all sites pulled together
are shown in (a), for epitopic sites in (b) and for the nonepitopic sites in (c).

second half of the 1990s, when it became fixed after
emerging repeatedly in different subpopulations.

4. DISCUSSION
Our simulation null model of HA sequence evolution isbased on a detailed codon-based substitution modelthat is similar to the widely used Markov chain-basedmodels. As in most such models, the direction ofsubstitutions is independent of the identity of the target
amino acid. The simulated data faithfully reproduce mostof the coarse characteristics of evolution observed in
the true empirical data (see the electronic supplementarymaterial). However, the null model fails to adequately
describe finer aspects of non-synonymous evolution. Our
analysis has focused on the discrepancies between the data
and the model in the probabilities of different amino acidsubstitutions at a site.
(a) Protein conservation
Owing to functional constraints, the protein shape isgenerally more conserved in the course of evolution thanin the underlying amino acid sequence (Lesk & Chothia1980). For the same reason, we might expect that
substitutions between biochemically dissimilar amino

Directionality in evolution of HA
Table 1. Amino acids whose average ES statistic is significantly
elevated or depressed (grey shade) at the 0.1% level, without
the Bonferroni correction (one-tailed tests). (The expected
number of false positives is 2.5. Notation k denotes the site;
dN/dS denotes the average dN/dS ratio for the site; X denotes
the amino acid at the site; nk ð/XÞ denotes the average numberof observed substitutions leading to amino acid X at site k; andsES ðk; XÞ denotes the average value of the ES statistic
corresponding to amino acid X at site k. Averages are takenover 20 parsimonious trees. Standard errors are not shown
because they are typically small relative to the averages (meanrelative standard error is 2.4%, maximum 13.8%). In italics are
the sites and amino acids whose average ES statistic issignificantly elevated or depressed at the 5% level, with the
Bonferroni correction (one-tailed tests).)site (k)epitope

3

dN/dS
1.8

92
106
121
124
133
137

E
D
A
A
A

1.5
0.7
1.4
1.6
2.0
2.3

138

A

2.2

142
144
145

A
A
A

1.1
1.5
5.0

156

B

2.8

159

B

2.2

167

D

1.5

186

B

2.2

190
193
194

B
B
B

1.0
1.5
2.0

213
220

D

0.6
1.7

226

D

6.7

229

D

3.1

233
248

D

0.9
1.6

275

C

0.9

amino
acid (X) nk ð/XÞ
Ile
Pro
Thr
Val
Thr
Gly
Asn
Ser
His
Ser
Val
Gly
Asp
Lys
Asn
Arg
Asp
Glu
Lys
Arg
Asn
His
Ala
Ile
Ile
Arg
Asn
Val
Lys
Ile
Pro
Val
Gly
Lys
Gln
Ile
Leu
Met
Phe
Pro
Ala
Ile
Lys
His
Ile
Ala
Asp

10.2
2.0
6.0
7.0
11.0
9.5
12.1
7.6
1.0
13.0
0.0
8.0
9.2
32.4
11.6
2.1
0.0
0.0
16.7
0.5
8.3
2.0
11.0
1.0
19.9
0.0
2.0
7.0
6.0
13.0
0.0
4.9
13.2
1.0
16.9
32.4
8.2
2.5
1.0
1.0
0.0
15.7
9.4
9.0
13.0
0.0
10.2

2461

50

sES ðk; XÞ
6.6
K4.9
5.2
4.7
8.3
6.1
6.2
6.2
K5.1
9.2
K7.5
6.4
7.4
25.7
4.8
K6.1
K6.5
K7.7
11.8
K5.8
5.0
K5.4
5.9
K4.3
17.4
K5.6
K7.7
6.4
5.1
8.9
K8.2
4.6
7.6
K8.8
15.2
14.1
6.7
K4.0
K4.3
K10.3
K10.8
11.2
K8.2
4.7
8.0
K4.7
5.3

acids would occur less frequently than between similaramino acids (Miyata & Yasunaga 1980).
Proc. R. Soc. B (2008)

S. Kryazhimskiy et al.

*
Figure 2. Homoplastic substitutions resulting in lysine at site145. The phylogenetic tree of HA is shown, with the root
marked by an asterisk. Branches where the sequence at thedescendent node has lysine at site 145 are in black. On thistree, lysine appeared 33 times at site 145. Branch lengths aremeasured in nucleotide substitutions.

Our results indicate that non-synonymous substitutions
at the second codon positions of HA occurred lessfrequently than expected under the mGY substitution
model, which does not account for differences in the rates of
non-synonymous substitutions between the three codon
positions (see the electronic supplementary material). Thismismatch indicates that different codon positions evolve, on
average, under different selective pressures, probably due tothe structure of the genetic code. Indeed, non-synonymouschanges at the second amino acid position tend to change
the properties of the encoded amino acid more radicallythan non-synonymous changes at the first or third positions(e.g. Haig & Hurst 1991; Urbina et al. 2006).
The mGYCl model provides a better overall fit to
the data (see the electronic supplementary material),
as it controls for the differences in probabilities of nonsynonymous substitutions between the three codon
positions. However, even under this model, amino acidsubstitutions are significantly more conservative with regard
to the physicochemical properties than expected from theirmutational opportunities (see the electronic supplementarymaterial). The conservation is lower in epitopic than in nonepitopic regions of HA (see the electronic supplementarymaterial). This conforms with previous observations (Ina &Gojobori 1994; Bush et al. 1999b) that epitopes are lessconstrained than the rest of the HA protein, but alsosuggests that the mGYCl model is not sufficient, even after
fitting site-specific dN/dS values and codon-position
corrections, to capture the amino acid level patterns of
HA protein evolution.
(b) Directional selection
The direction of amino acid substitutions has site-specificbiases. Specifically, in both epitopic and non-epitopic

2462 S. Kryazhimskiy et al.

Directionality in evolution of HA

sites, homoplasies are more frequent than expected: aminoacid substitutions at a site tend to repeatedly give rise
to identical amino acids, whereas other specific aminoacids are avoided (figure 1; table 1; see the electronic
supplementary material). These deviations are not causedby the difference in the number of non-synonymous
substitutions between the data and the simulations, since
the simulations approximately preserve this number (seethe electronic supplementary material).
Although we can pinpoint the amino acids that arose
more frequently than expected (table 1), it is possible thatsome of these amino acids were driven to high frequencies
by hitch-hiking along with other selected amino acids
rather than due to direct selection. It is difficult to assessthe degree to which hitch-hiking affects our site-by-sitestatistics. However, since the hitch-hiking hypothesis
implies that selective sweeps occurred at least at somesites, it is unlikely that hitch-hiking accounts for all
significant results presented in table 1.
Therefore, the differences between model and data
have to do with natural selection affecting the variability ofamino acids at a site. These differences could be due toeither negative (selective constraint) or positive selection.
Under a conceivable negative selection scenario, if just a
few amino acids are permitted at an amino acid site,
substitutions between these amino acids predominate(Bazykin et al. 2007), leading to more observed cases ofhomoplastic evolution in the data than in the simulation.
On the other hand, positive selection could also lead torepeated substitutions towards identical novel aminoacids, either due to repeated changes of the selective
landscape or the independent response of different
evolving lineages to the similar selective pressure.
Several arguments suggest that the observed statisticalresults are primarily due to the action of positive, rather
than negative, selection. First, homoplastic substitutions at
rates higher than expected on the basis of their mutational
opportunities, as shown by the significant elevation of the
ES statistic, are unlikely to be due to negative selection.
Second, most of the sites with high ES statistics listed intable 1, and all of the sites for which the results are mostsignificant (indicated in italic in table 1), have dN/dS values
greater than 1, suggesting that positive selection is the
predominant regime at those sites. Finally, most of thesesites are located in the epitopes (table 1), for which theimportance of positive selection is well known (Bush et al.
1999a; Yang 2000). Therefore, it is reasonable to assume
that most amino acids with the positive ES statistic listed
in table 1 were positively selected, at the sites specified.
Our results, therefore, suggest that the shape of thefitness landscape at each site has a major effect on
the evolution of HA by positive selection, in contrast
to the unrealistic assumption of the identical fitness oftarget amino acids made by most Markov chain models.
Other authors recently came to similar conclusions(Kosakovsky Pond et al. in press). Positive selection is
directional, in the sense that the substitution probabilities
into different amino acids are non-uniform. Frequently,
directional selection leads to repeated (homoplastic)
substitutions into the same amino acid at a site and/or
substitutions into amino acids that are unlikely to be
based on mutation probabilities alone. Ascribing theselection coefficients to individual amino acids at a site isa subject for future research.
Proc. R. Soc. B (2008)

S.K. gratefully acknowledges support by the Burroughs
Wellcome Fund Training programme in Biological Dynamics(no 1001782) and the DARPA grant HR0011-05-1-0057.
G.A.B. was partially supported by the Molecular and Cellular
Biology programme of the Russian Academy of Sciences.
J.B.P. and S.K. were funded by grants from the Burroughs
Wellcome Fund and the James S. McDonnell Foundation.
J.D. gratefully acknowledges support by the Natural Sciencesand Engineering Research Council of Canada and the
DARPA grant HR0011-05-1-0057.

REFERENCES
Bazykin, G. A., Kondrashov, F. A., Brudno, M., Poliakov, A.,
Dubchak, I. & Kondrashov, A. S. 2007 Extensive
parallelism in protein evolution. Biol Direct 2, 20.
(doi:10.1186/1745-6150-2-20)
Bush, R. M., Bender, C. A., Subbarao, K., Cox, N. J. & Fitch,
W. M. 1999a Predicting the evolution of human influenza
A. Science 286, 1921–1925. (doi:10.1126/science.286.
5446.1921)
Bush, R. M., Fitch, W. M., Bender, C. A. & Cox, N. J. 1999b
Positive selection on the H3 hemagglutinin gene of human
influenza virus A. Mol. Biol. Evol. 16, 1457–1465.
Bush, R. M., Smith, C. B., Cox, N. J. & Fitch, W. M. 2000
Effects of passage history and sampling bias on phylogenetic reconstruction of human influenza A evolution. Proc.
Natl Acad. Sci. USA 97, 6974–6980. (doi:10.1073/pnas.
97.13.6974)
Carroll, S. B. 2003 Genetics and the making of Homo sapiens.
Nature 422, 849–857. (doi:10.1038/nature01495)
Chen, L., Perlina, A. & Lee, C. J. 2004 Positive selection
detection in 40,000 human immunodeficiency virus(HIV ) type 1 sequences automatically identifies drug
resistance and positive fitness mutations in HIV protease
and reverse transcriptase. J. Virol. 78, 3722–3732. (doi:10.
1128/JVI.78.7.3722-3732.2004)
Dimmic, M. W., Mindell, D. P. & Goldstein, R. A. 2000
Modeling evolution at the protein level using an adjustableamino acid fitness model. Pac. Symp. Biocomput. 5, 18–29.
Fitch, W. M., Leiter, J. M. E., Li, X. & Palese, P. 1991
Positive Darwinian evolution in human influenza A
viruses. Proc. Natl Acad. Sci. USA 88, 4270–4274.
(doi:10.1073/pnas.88.10.4270)
Fitch, W. M., Bush, R. M., Bender, C. A. & Cox, N. J. 1997
Long term trends in the evolution of H(3)NA1 human
influenza type A. Proc. Natl Acad. Sci. USA 94,
7712–7718. (doi:10.1073/pnas.94.15.7712)
Fornasari, M. S., Parisi, G. & Echave, J. 2002 Site-specificamino acid replacement matrices from structurally
constrained protein evolution simulations. Mol. Biol.
Evol. 19, 352–356.
Forsberg, R. & Christiansen, F. B. 2003 A codon-basedmodel of host-specific selection in parasites, with anapplication to the influenza A virus, Mol. Biol. Evol. 20,
1252–1259. (doi:10.1093/molbev/msg149)
Goldman, N. & Yang, Z. 1994 A codon-based model of
nucleotide substitution for protein-coding DNA sequences.
Mol. Biol. Evol. 11, 725–736.
Guindon, S., Rodrigo, A. G., Dyer, K. A. & Huelsenbeck,J. P. 2004 Modeling the site-specific variation of selection
patterns along lineages. Proc. Natl Acad. Sci. USA 101,
12 957–12 962. (doi:10.1073/pnas.0402177101)
Haig, D. & Hurst, L. D. 1991 A quantitative measure of error
minimization in the genetic code. J. Mol. Evol. 33,
412–417. (doi:10.1007/BF02103132)
Halpern, A. L. & Bruno, W. J. 1998 Evolutionary distances
for protein-coding sequences: modeling site-specificresidue frequencies. Mol. Biol. Evol. 15, 910–917.

Directionality in evolution of HA
Holmes, E. C. et al. 2005 Whole-genome analysis of human
influenza A virus reveals multiple persistent lineages andreassortment among recent H3N2 viruses. PLoS Biol. 3,
e300. (doi:10.1371/journal.pbio.0030300)
Ina, Y. & Gojobori, T. 1994 Statistical analysis of nucleotide
sequences of the hemagglutinin gene of human influenza
A viruses. Proc. Natl Acad. Sci. USA 91, 8388–8392.
(doi:10.1073/pnas.91.18.8388)
Kimura, M. 1977 Preponderance of synonymous changes asevidence for the neutral theory of molecular evolution.
Nature 267, 275–276. (doi:10.1038/267275a0)
Kimura, M. 1983 The neutral theory of molecular evolution.
Cambridge, UK: Cambridge University Press.
Kosakovsky Pond, S. L., Poon, A. F. Y., Brown, A. J. L. &Frost, S. D. W. In press. A maximum likelihood methodfor detecting directional evolution in protein sequences
and its application to influenza A virus, Mol. Biol. Evol.
(doi:10.1093/molbev/msn123)
Koshi, J. M. & Goldstein, R. A. 1998 Models of natural
mutations including site heterogeneity. Proteins 32,
289–295. (doi:10.1002/(SICI )1097-0134(19980815)32:
3!289::AID-PROT4O3.0.CO;2-D)
Koshi, J. M., Mindell, D. P. & Goldstein, R. A. 1999 Using
physical-chemistry-based substitution models in phylogenetic analysis of HIV-1 subtypes. Mol. Biol. Evol. 16,
173–179.
Lartillot, N. & Philippe, H. 2004 A Bayesian mixture model
for across-sites heterogeneities in the amino-acid replacement process. Mol. Biol. Evol. 21, 1095–1109. (doi:10.
1093/molbev/msh112)
Lesk, A. M. & Chothia, C. 1980 How different amino acidsequences determine similar protein structures: thestructure and evolutionary dynamics of the globins.
J. Mol. Biol. 136, 225–270. (doi:10.1016/00222836(80)90373-3)
Lindstrom, S. E., Cox, N. J. & Klimov, A. 2004 Genetic
analysis of human H2N2 and early H3N2 influenza
viruses, 1957–1972: evidence for genetic divergence andmultiple reassortment events. Virology 328, 101–119.
(doi:10.1016/j.virol.2004.06.009)
Mayrose, I., Doron-Faigenboim, A., Bacharach, E. & Pupko,
T. 2007 Towards realistic codon models: among sitevariability and dependency of synonymous and nonsynononymous rates. Bioinformatics 23, i319–i327. (doi:10.
1093/bioinformatics/btm176)
Meyer, W. J., Wood, J. M., Major, D., Robertson, J. S.,
Webster, R. G. & Katz, J. M. 1993 Influence of host cellmediated variation on the international surveillance of
influenza A (H3N2) viruses. Virology 196, 130–137.
(doi:10.1006/viro.1993.1461)
Miyata, T. & Yasunaga, T. 1980 Molecular evolution ofmRNA: a method for estimating evolutionary rates ofsynonymous and amino acid substitutions from homologous nucleotide sequences and its application. J. Mol.
Evol. 16, 23–36. (doi:10.1007/BF01732067)
Muse, S. V. & Gaut, B. S. 1994 A likelihood approach forcomparing synonymous and nonsynonymous nucleotide
substitution rates, with application to the chloroplast
genome, Mol. Biol. Evol. 11, 715–724.
Nelson, M. I. & Holmes, E. C. 2007 The evolution of
epidemic influenza. Nat. Rev. Genet. 8, 196–205. (doi:10.
1038/nrg2053)
Nielsen, R. & Yang, Z. 1998 Likelihood models for detecting
positively selected amino acid sites and applications to theHIV-1 envelope gene. Genetics 148, 929–936.
Nielsen, R. & Yang, Z. 2003 Estimating the distribution of
selection coefficients from phylogenetic data with applications to mitochondrial and viral DNA. Mol. Biol. Evol.
20, 1231–1239. (doi:10.1093/molbev/msg147)
Proc. R. Soc. B (2008)

S. Kryazhimskiy et al.

2463

Nixon, K. C. 1999 The parsimony ratchet, a new method for
rapid parsimony analysis. Cladistics 15, 407–414. (doi:10.
1111/j.1096-0031.1999.tb00277.x)
Plotkin, J. & Dushoff, J. 2003 Codon bias and frequencydependent selection on the hemagglutinin epitopes ofinfluenza A virus. Proc. Natl Acad. Sci. USA 100,
7152–7157. (doi:10.1073/pnas.1132114100)
Plotkin, J., Dushoff, J. & Levin, S. A. 2002 Hemagglutinin
sequence clusters and the antigenic evolution of influenza
A virus. Proc. Natl Acad. Sci. USA 99, 6263–6268. (doi:10.
1073/pnas.082110799)
Robinson, D. M., Jones, D. T., Kishino, H., Goldman, N. &Thorne, J. 2007 Protein evolution with dependence
among codons due to tertiary structure. Mol. Biol. Evol.
20, 1692–1704. (doi:10.1093/molbev/msg184)
Rocha, E. P., Xu, X., Hall, H. E., Allen, J. R., Regnery, H. L.
& Cox, N. J. 1993 Comparison of 10 influenza A (H1N1
and H3N2) haemagglutinin sequences obtained directlyfrom clinical specimens to those of MDCK cell- and egggrown viruses. J. Gen. Virol. 74, 2513–2518. (doi:10.1099/
0022-1317-74-11-2513)
Rodrigue, N., Lartillot, N., Bryant, D. & Philippe, H. 2005
Site interdependence attributed to tertiary structure inamino acid sequence evolution. Gene 347, 207–217.
(doi:10.1016/j.gene.2004.12.011)
Sainudiin, R., Wong, W. S. W., Yogeeswaran, K., Nasrallah,
J. B., Yang, Z. & Nielsen, R. 2005 Detecting site-specificphysicochemical selective pressures: applications to theclass I HLA of the human major histocompatibility
complex and the SRK of the plant sporophytic selfincompatibility system. J. Mol. Evol. 60, 315–326. (doi:10.
1007/s00239-004-0153-1)
Seoighe, C. et al. 2007 A model of directional selection applied
to the evolution of drug resistance in HIV-1. Mol. Biol.
Evol. 24, 1025–1031. (doi:10.1093/molbev/msm021)
Shih, A. C.-C., Hsiao, T.-C., Ho, M.-S. & Li, W.-H. 2007
Simultaneous amino acid substitutions at antigenic sites
drive influenza A hemagglutinin evolution. Proc. Natl
Acad. Sci. USA 104, 6283–6288. (doi:10.1073/pnas.
0701396104)
Simpson, E. H. 1949 Measurement of diversity. Nature 163,
688. (doi:10.1038/163688a0)
Suzuki, Y. 2006 Natural selection on the influenza A virus
genome. Mol. Biol. Evol. 23, 1902–1911. (doi:10.1093/
molbev/msl050)
Swofford, D. L. 2002 PAUP . Phylogenetic analysis usingparsimony (and other methods), v. 4. Sunderland, MA:
Sinauer Associates.
Thompson, J. D., Higgins, D. G. & Gibson, T. J. 1994
CLUSTAL W: improving the sensitivity of progressive
multiple sequence alignment through sequence weighting,position-specific gap penalties and weight matrix choice.
Nucleic Acids Res. 22, 4673–4680. (doi:10.1093/nar/22.22.
4673)
Thorne, J. L., Goldman, N. & Jones, D. T. 1996 Combining
protein evolution and secondary structure. Mol. Biol. Evol.
13, 666–673.
Urbina, D., Tang, B. & Higgs, P. G. 2006 The response ofamino acid frequencies to directional mutation pressure inmitochondrial genome sequences is related to physicalproperties of the amino acids and to the structure of thegenetic code. J. Mol. Evol. 62, 340–361. (doi:10.1007/
s00239-005-0051-1)
Wilson, I. A. & Cox, N. J. 1990 Structural basis of immune
recognition of influenza virus hemagglutinin. Annu. Rev.
Immunol. 8, 737–771. (doi:10.1146/annurev.iy.08.040
190.003513)
Wolf, Y. I., Viboud, C., Holmes, E. C., Koonin, E. V. &Lipman, D. J. 2006 Long intervals of stasis punctuated bybursts of positive selection in the seasonal evolution of

2464 S. Kryazhimskiy et al.

Directionality in evolution of HA

influenza A virus. Biol. Direct 1, 34–53. (doi:10.1186/
1745-6150-1-34)
Wong, W. S. W., Sainudiin, R. & Nielsen, R. 2006
Identification of physicochemical selective pressure on
protein encoding nucleotide sequences. BMC Bioinform.
7, 148. (doi:10.1186/1471-2105-7-148)
Yang, Z. 2000 Maximum likelihood estimation on large
phylogenies and analysis of adaptive evolution in human
influenza A. J. Mol. Evol. 51, 423 –432.
Yang, Z. & Bielawski, J. P. 2000 Statistical methods for
detecting molecular adaptation. Trends Ecol. Evol. 15,
496–503. (doi:10.1016/S0169-5347(00)01994-7)
Yang, Z. & Nielsen, R. 2000 Estimating synonymous andnonsynonymous substitution rates under realistic
evolutionary models. Mol. Biol. Evol. 17, 32–43.
Yang, Z. & Nielsen, R. 2002 Codon-substitution models for
detecting molecular adaptation at individual sites along
specific lineages. Mol. Biol. Evol. 19, 908–917.

Proc. R. Soc. B (2008)

Yang, Z. & Nielsen, R. 2008 Mutation–selection models ofcodon substitution and their use to estimate selective
strengths on codon usage. Mol. Biol. Evol. 25, 568–579.
(doi:10.1093/molbev/msm284)
Yang, Z., Nielsen, R. & Hasegawa, M. 1998 Models of aminoacid substitution and applications to mitochondrial
protein evolution. Mol. Biol. Evol. 15, 1600–1611.
Yang, Z., Nielsen, R., Goldman, N. & Pedersen, A.-M. K.
2000 Codon-substitution models for heterogeneous
selection pressure at amino acid sites. Genetics 155,
431–449.
Yeager, M., Kumar, S. & Hughes, A. L. 1997 Sequence
convergence in the peptide-binding region of primate androdent MHC class Ib molecules. Mol. Biol. Evol. 14,
1035–1041.
Zhang, J. & Kumar, S. 1997 Detection of convergent andparallel evolution at the amino acid sequence level. Mol.
Biol. Evol. 14, 527–536.