key: cord-328681-jf2mj16z authors: Yang, Ziheng; Bielawski, Joseph P. title: Statistical methods for detecting molecular adaptation date: 2000-12-01 journal: Trends Ecol Evol DOI: 10.1016/s0169-5347(00)01994-7 sha: doc_id: 328681 cord_uid: jf2mj16z The past few years have seen the development of powerful statistical methods for detecting adaptive molecular evolution. These methods compare synonymous and nonsynonymous substitution rates in protein-coding genes, and regard a nonsynonymous rate elevated above the synonymous rate as evidence for darwinian selection. Numerous cases of molecular adaptation are being identified in various systems from viruses to humans. Although previous analyses averaging rates over sites and time have little power, recent methods designed to detect positive selection at individual sites and lineages have been successful. Here, we summarize recent statistical methods for detecting molecular adaptation, and discuss their limitations and possible improvements. 'I t has been proved remarkably difficult to get compelling evidence for changes in enzymes brought about by selection, not to speak of adaptive changes' 1 . Although Darwin's theory of evolution by natural selection is generally accepted by biologists for morphological traits (including behavioural and physiological), the importance of natural selection in molecular evolution has long been a matter of debate. The neutral theory 2 maintains that most observed molecular variation -both polymorphism within species and divergence between species -is due to random fixation of selectively neutral mutations. Well established cases of molecular adaptation have been rare 3 . Several tests of neutrality have been developed and applied to real data, and although they are powerful enough to reject strict neutrality in many genes, they rarely provide unequivocal evidence for positive darwinian selection. Most convincing cases of adaptive molecular evolution have been identified through comparison of synonymous (silent; d S ) and nonsynonymous (amino acid-changing; d N ) substitution rates in protein-coding DNA sequences, thus providing fascinating case studies of natural selection in action on the protein molecule. Selected examples are listed in Table 1 ; see Hughes 4 for detailed descriptions of many case studies. Here, we summarize recent methodological developments that improve the power to detect adaptive molecular evolution, and examine their strengths and weaknesses , so that they can be used to detect more cases of molecular adaptation. Traditionally, synonymous and nonsynonymous substitution rates (Box 1) are defined in the context of comparing two DNA sequences, with d S and d N as the numbers of synonymous and nonsynonymous substitutions per site, respectively 5 . Thus, the ratio ϭ d N /d S measures the difference between the two rates and is most easily understood from a mathematical description of a codon substitution model (Box 2). If an amino acid change is neutral, it will be fixed at the same rate as a synonymous mutation, with ϭ 1. If the amino acid change is deleterious, purifying selection (Box 1) will reduce its fixation rate, thus Ͻ 1. Only when the amino acid change offers a selective advantage is it fixed at a higher rate than a synonymous mutation, with Ͼ 1. Therefore, an ratio significantly higher than one is convincing evidence for diversifying selection. The codon-based analysis (Box 2) cannot infer whether synonymous substitutions are driven by mutation or selection, but it does not assume that synonymous substitutions are neutral. For example, highly biased codon usage can be caused by both mutational bias and selection (e.g. for translational efficiency 6 ), and can greatly affect synonymous substitution rates. However, by employing parameters j for the frequency of codon j in the model (Box 2), estimation of Ziheng Yang and Joseph P. Bielawski The past few years have seen the development of powerful statistical methods for detecting adaptive molecular evolution. These methods compare synonymous and nonsynonymous substitution rates in protein-coding genes, and regard a nonsynonymous rate elevated above the synonymous rate as evidence for darwinian selection. Numerous cases of molecular adaptation are being identified in various systems from viruses to humans. Although previous analyses averaging rates over sites and time have little power, recent methods designed to detect positive selection at individual sites and lineages have been successful. Here, we summarize recent statistical methods for detecting molecular adaptation, and discuss their limitations and possible improvements. substitution rates will fully account for codon-usage bias (Box 1), irrespective of its source. Because parameter is a measure of selective pressure on a protein, it differentiates codon-based analyses from the more general tests of neutrality proposed in population genetics 7,8 . These general tests often lack the power to determine the sources of the departure from the strict neutral model, such as changes in population size, fluctuating environment or different forms of selection. Two classes of methods have been suggested to estimate d N and d S between two protein-coding DNA sequences. The first class includes over a dozen intuitive methods developed since the early 1980s (Refs 5,9-15). These methods involve the following steps: counting synonymous (S) and nonsynonymous (N) sites in the two sequences, counting synonymous and nonsynonymous differences between the two sequences, and correcting for multiple substitutions at the same site. S and N are defined as the sequence length multiplied by the proportions of synonymous and nonsynonymous changes before selection on the protein 14,16 . Most of these methods make simplistic assumptions about the nucleotide substitution process and also involve ad hoc treatment of the data that cannot be justified 14,15 ; therefore, we refer to these methods of estimating d N and d S as approximate methods. The methods of Miyata and Yasunaga 5 , and Nei and Gojobori 9 , assume an equal rate for transitions (T ↔ C and A ↔ G) and transversions (T,C ↔ A,G), as well as a uniform codon usage. Because transitions at the third 'wobble' position are more likely to be synonymous than transversions, ignoring the transition/ transversion rate ratio leads to underestimation of S and overestimation of N (Ref. 10). Efforts have been taken to incorporate the transition/transversion rate bias (Box 1) when counting sites and differences 10-14 . The effect of Prior probability: the probability of an event (such as a site belonging to a site class) before the collection of data. Positive selection: darwinian selection fixing advantageous mutations with positive selective coefficients. The term is used interchangeably with molecular adaptation and adaptive molecular evolution. Posterior probability: the probability of an event conditional on the observed data, which reflects both the prior assumption and information in the data. Purifying selection: natural selection against deleterious mutations with negative selective coefficients. The term is used interchangeably with negative selection or selective constraints. Synonymous substitution: a nucleotide substitution that does not change the encoded amino acid. Transition/transversion rate bias: unequal substitution rates between nucleotides, with a higher rate for transitions (changes between T and C and between A and G) than transversions (all other changes). biased codon usage has largely been ignored 17 ; however, extreme codon-usage bias can have devastating effects on the estimation of d N and d S (see the next section) 15,18 . A recent ad hoc method 15 incorporates both transition and codon-usage biases. The second class is the maximum likelihood (ML) method based on explicit models of codon substitution (Box 2) 16, 19 . Parameters in the model (i.e. sequence divergence t, transition/transversion rate ratio and the d N /d S ratio ) are estimated from the data by ML, and are used to calculate d N and d S according to their definitions 15,16,20 . A major feature of the method is that the model is formulated at the level of instantaneous rates (where there is no possibility for multiple changes) and that probability theory accomplishes all difficult tasks in one step: estimating mutational parameters, such as ; correcting for multiple hits; and weighting pathways of change between codons. Statistical tests can be used to test whether d N is significantly higher than d S . For approximate methods, a normal approximation is applied to d N Ϫ d S . For ML, a likelihoodratio test can be used. In this case, the null model has fixed at 1, whereas the alternative model estimates as a free parameter. Twice the log-likelihood difference between the two models is compared with a 2 distribution with one degree of freedom to test whether is different from one. Computer simulation has been used to examine the performance of different estimation methods; the findings are consistent with observations made in real data analyses 14,15,19 . We demonstrate the effects of different estimation procedures using human and orangutan ␣ 2 -globin genes ( Table 2) . For comparison, different assumptions are made in ML concerning the transition/transversion rate bias and the codon-usage bias. The simpler models are each rejected when compared with more complex models by likelihood-ratio tests, confirming biased transition rates and codon usage. Thus, estimates from ML accounting for both biases (Model 8, Table 2 ) are expected to be the most reliable. We make the following observations: • Assumptions appear to matter more than methods. The approximate methods and ML produce similar results under similar assumptions. The method of Nei and Gojobori is similar to ML under a model that ignores both transition/transversion bias and codon-usage bias (Model 1, Table 2 ), whereas the methods of Ina and Li are similar to ML under a model accounting for the transition/transversion bias but ignoring codon-usage bias (Model 2, Table 2 ). The method of Yang and Nielsen 15 is similar to ML under a model accounting for both biases (Model 6, Table 2 ). However, for distantly related sequences, ad hoc treatment in approximate methods can lead to serious biases even under the correct assumptions 19 . • Ignoring the transition/transversion rate bias leads to underestimation of S, overestimation of d S and underestimation of the ratio 10 . • Codon-usage bias in these data has the opposite The codon is considered the unit of evolution. The substitution rate from codons i to j (i j) is given as: Parameter is the transition/transversion rate ratio, j is the equilibrium frequency of codon j and (ϭ d N /d S ) measures the selective pressure on the protein. The q ij are relative rates because time and rate are confounded in such an analysis. Given the rate matrix Q ϭ {q ij }, the transition probability matrix over time t is calculated as: where p ij (t) is the probability that codon i becomes codon j after time t. Likelihood calculation on a phylogeny involves summing over all possible codons in extinct ancestors (internal nodes of the tree effect to the transition/transversion bias; ignoring codon-usage bias leads to overestimation of S, underestimation of d S and overestimation of . This gene is extremely GC-rich at the third codon position, with base frequencies at 9% (T), 52% (C), 1% (A) and 37% (G). Most changes at the third position (before selection at the amino acid level) are transversions between C and G. Thus, the number of synonymous sites is less than half that expected under equal base and codon frequencies. Although, in theory, the bias caused by unequal codon frequencies can be in the opposite direction 15 , we have not encountered a real gene showing that pattern. Such codon-usage bias appears to have misled previous analyses examining the relationship between the GC content at silent sites and d S , because those studies ignored the codon-usage bias when estimating d S (Ref. 21). • Different methods can produce different estimates, even when the sequences are highly similar. The sequences used in Table 2 are only about 10% different at silent sites and Ͻ1% different at nonsynonymous sites; however, estimates of are three times different. Because all estimation procedures partition the total numbers of sites and differences into synonymous and nonsynonymous categories, underestimation of one means overestimation of the other, thus resulting in large errors in the ratio. If, for most of the time, a gene evolves under purifying selection but is occasionally subject to episodes of adaptive change 22 , a comparison between two distantly related sequences is unlikely to yield a d N /d S ratio significantly greater than one. Methods have been developed to detect positive selection (Box 1) along specific lineages on a phylogeny. If the gene sequences of the extinct ancestors were known, it would be straightforward to use the pairwise methods discussed above. Thus, Messier and Stewart 23 inferred ancestral lysozyme gene sequences through phylogenetic analysis 24,25 , and used them to calculate d N and d S for each branch in the phylogeny. Their analysis identified two lineages in a primate phylogeny with highly elevated nonsynonymous substitution rates. The same approach was taken in a test of relaxed selective constraint in the rhodopsin gene of cave-dwelling crayfishes 26 . There are also likelihood models that allow different ratios for branches in a phylogeny 18,27 . Using such models, likelihood-ratio tests can be constructed to test hypotheses. For example, the ratio for a predefined lineage can be either fixed at one or estimated as a free parameter. The likelihood values under those two models can be compared, to test whether Ͼ 1 in that lineage. Similarly, a model assuming a single for all lineages (the one-ratio model) can be compared with another model assuming an independent for each lineage (the free-ratio model), to test the neutral prediction that the ratio is identical among lineages 18,27 . It should be noted that variation in the ratio among lineages is a violation of the strictly neutral model 2,18,28,29 , but it is not sufficient evidence for adaptive evolution. In particular, if nonsynonymous mutations are slightly deleterious, they will have a higher probability of fixation in a small population than in a large one 30 , and thus lineages of different population sizes will have different ratios. Besides positive selection, relaxed selective constraint can also elevate the ratio -it might be difficult to distinguish the two if the estimated is not larger than one. Furthermore, it is incorrect to use the free-ratio model to identify lineages of interest and then to perform further tests on the ratios for those lineages using the same data without any correction 27 . The statistical-estimation theory used in the methods discussed in this review can be explained with the following simple hypothetical example. Suppose that a population is an admixture of two groups of people in the proportions 60% and 40%, and a certain disease occurs at a rate of 1% in Group I and of 0.1% in Group II. Suppose a random sample of 100 individuals is taken from the population, what is the probability that three of them carry the disease? The probability that a random individual carries the disease (D) is an average over the two groups (G 1 and G 2 ): p = P(D) = P(G 1 ) × P(DG 1 ) + P(G 2 ) × P(DG 2 ) = 0.6 × 0.01 + 0.4 × 0.001 = 0.0064 (1) Similarly, the probability that an individual does not carry the disease is: The probability that three out of 100 individuals carry the disease is given by the binomial probability: ( If Eqn 3 involves an unknown parameter [such as the rate P(D|G 1 ) in Group I], that parameter can be estimated by maximizing Eqn 3. In that case, Eqn 3 gives the probability of observing the data (sample) and is called the likelihood function. The second question is to calculate the probability that an individual in the sample who carries the disease is from Group I. The Bayes theorem gives this probability as: P(G 1 D) = P(G 1 ) × P(DG 1 )/P(D) = 0.6 × 0.01/0.0064 = 0.94 (4) Note that this is just the proportion of the contribution from Group I to P(D) in Eqn 1. Thus, this individual is most likely to be from Group I. Similarly, a healthy individual in the sample is more likely to be from Group I than from Group II because P(G 1 D -) = P(G 1 ) × P(D -G 1 )/P(D -) = 0.6 × 0.99/0.9936 = 0.5978 and P(G 2 D -) = 1 -P(G 1 D -) = 0.4022 (5) In methods for inferring sites under positive selection 36,37 , we let D in the example be the data at a site and G i be the ith site class with the d N /d S ratio i . The probability of observing data at a site is then an average over the site classes (Eqn 1). The product of such probabilities over sites constitutes the likelihood (Eqn 3), from which we estimate any unknown parameters, such as the branch lengths and parameters in the distribution over sites. After the parameters are estimated, we use the Bayes theorem to calculate the probability that any site, given data at that site, is from each site class (Eqns 4 and 5). Another straightforward application of the theory is ancestral sequence reconstruction; in this case, we replace G i with a reconstruction (characters at interior nodes of the phylogeny) at a site. When we calculate the likelihood function, the probability of data at a site P(D) is a sum over all possible ancestral reconstructions (G i s) (Eqns 1 and 2). After parameters are estimated, the reconstruction that makes the greatest contribution to P(D) is the most likely (Eqns 4 and 5) 24 . The Bayes method discussed here is known as the empirical Bayes, because it uses estimates of parameters and does not account for their sam-pling errors. This might be a concern if parameters are estimated from small samples or if the posterior probabilities are sensitive to parameter estimates. An alternative approach is the hierarchical Bayes method, which accounts for the uncertainty in unknown parameters by averaging over their prior distribution. Note that the reconstructed ancestral sequences 24 , as well as the inferred site classes in the site-class models 36,37 , are pseudo data and involve systematic biases. To appreciate such biases, note that in the previous example, the Bayes calculations (Eqns 4 and 5) predict that each of the 100 individuals in the sample, healthy or sick, are from Group I. Although this is the best prediction, the accuracy is low. If such inferred group identities are used for further statistical analysis, misleading results might follow. Methods based on ancestral reconstruction might not provide reliable statistical tests because they ignore errors and biases in reconstructed ancestral sequences (Box 3). The ML method has the advantage of not relying on reconstructed ancestral sequences. It can also easily incorporate features of DNA sequence evolution, such as the transition/transversion rate bias and codon-usage bias, and is thus based on a more realistic evolutionary model. When likelihood-ratio tests suggest adaptive evolution along certain lineages, ancestral reconstruction might be useful to pinpoint the involved amino acids and to infer ancestral proteins, which can be synthesized and examined in the laboratory 31,32 . The methods discussed so far assume that all amino acid sites are under the same selective pressure, with the same ratio. The analysis effectively averages the ratio across all sites and positive selection is detected only if that average is Ͼ1. This appears to be a conservative test of positive selection because many sites might be under strong purifying selection owing to functional constraint, with the ratio close to zero. A few recent studies addressed this problem. Fitch and colleagues 33,34 used parsimony to reconstruct ancestral DNA sequences, and counted changes at each codon site along branches of the tree. They tested whether the proportion of nonsynonymous substitutions at each site is greater than the average over all sites in the sequence. Suzuki and Gojobori 35 took a more systematic approach. For each site in the sequence, they estimated the numbers of synonymous and nonsynonymous sites and differences along the tree using reconstructed ancestral sequences, and then tested whether the proportion of nonsynonymous substitutions differed from the neutral expectation ( ϭ 1). Suzuki and Gojobori's criterion is more stringent than Fitch et al.'s, because the ratio averaged over sites is almost always Ͻ 1. These methods are expected to require many sequences in the data set so that there are enough changes at individual sites. Furthermore, the reliability of significance values produced by these methods might be affected by the use of ancestral reconstruction, which is most unreliable at the positively selected or variable sites 24 , and by codon composition bias, which is most extreme at a single site. In a likelihood model, it is impractical to use one parameter for each site. The standard approach is to use a statistical distribution to describe the variation of among sites; for example, we might assume several classes of sites in the protein with different ratios 36,37 . The test of positive selection then involves two major steps: first, to test whether sites exist where Ͼ 1, which is achieved by a likelihood-ratio test comparing a model that does not allow for such sites with a more general model that does; and second, to use the Bayes theorem to identify positively selected sites when they exist. Sites having high posterior probabilities (Box 1) for site classes with Ͼ 1 are potential targets of diversifying selection. The theory is explained in Box 3 (Refs 20,36,37) . Nielsen and Yang 36 implemented a likelihood-ratio test based on two simple models. The null model, M1 (neutral), assumes a class of conserved sites with ϭ 0 and another class of neutral sites with ϭ 1. The alternative model, M2 (selection), adds a third class of sites with estimated from the data. (The model codes are those used in the codeml program in the PAML package.) If M2 fits the data significantly better than M1 and the estimated ratio for the third class in M2 is Ͼ1, then some sites are under diversifying selection. Zanotto et al. 38 used this test to identify several sites under strong positive selection in the nef gene of HIV, whereas both pairwise comparison and slidingwindow analysis failed. This comparison was later found to lack power in some genes because M1 does not account for sites with 0 Ͻ Ͻ 1 and the third class in M2 is forced to account for such sites 37 . Thus, Yang et al. 37 implemented several new models. For example, the beta distribution (M7 beta) is a flexible null model with 0 Ͻ Ͻ 1, and can be compared with an alternative that adds an additional site class with estimated (M8 beta&). A general discrete model (M3) was also implemented 37 . These models identified positive selection in six out of ten genes the authors analysed. Figure 1 shows the use of a discrete model (M3) with three classes to identify sites under diversifying selection in abalone sperm lysin 39 . The methods discussed above assume that there are heterogeneous classes of amino acid sites but that we do not know a priori which class each site is from. Such 'fishing-expedition' studies might be useful in generating hypotheses for laboratory investigation because they could identify crucial amino acids whose changes have offered a selective advantage in Nature's evolutionary experiment. For example, amino acid residues under diversifying selection were inferred in analyses of HIV-1 nef (Ref. 38) and env (Ref. 40) genes, which might constitute unidentified viral epitopes. Alternatively, we might wish to test an a priori hypothesis that certain structural and functional domains of the protein are under positive selection. In such cases, likelihood models can be constructed that assign and estimate different parameters for sites from different structural and functional domains 20 . All the methods for detecting positive selection reviewed here appear to be conservative. They detect selection only if d N is higher than d S -selection that does not cause excessive nonsynonymous substitutions, such as balancing selection, might not be detected. The pairwise comparison has little power because it averages the ratio over sites and over time. Methods for detecting selection along lineages work only if the ratio averaged over all sites is Ͼ1. Similarly, the test of positive selection at sites works only if the ratio averaged over all branches is Ͼ1. If adaptive evolution occurs only in a short time interval and affects only a few crucial amino acids, none of the methods is likely to succeed. Constancy of selective pressure at sites appears to be a much more serious assumption than constancy among lineages, especially for genes likely to be under continuous selective pressure, such as the HIV env gene. Indeed, models of variable selective pressures among sites 36,37 have been successful in detecting positive selection, even in a background of overwhelming purifying selection indicated by an average ratio much smaller than one 37, 38, 41, 42 . Models that allow to vary among both lineages and sites should have increased power. The methods discussed here also assume the same ratio for all possible amino acid changes; for example, at a positively selected site, all amino acid changes are assumed to be advantageous, which is unrealistic. Although amino acid substitution rates are known to correlate with their chemical properties, the relationship is poorly understood 43, 44 . It is also not entirely clear how to define positive selection in a model accounting for chemical properties. It will be interesting to perform computer simulations to examine the power of various detection methods and to investigate how this is affected by important factors, such as the size of the gene, sampling of species (sequences) and the level of sequence divergence. Including more sequences in the data should improve the power of site-based analyses. Sequence divergence is also important because neither very similar nor very divergent sequences contain much information. Very divergent sequences might also be associated with problems with alignment and unequal nucleotide compositions in different species. Analyses discussed here, which require information from both synonymous and nonsynonymous substitutions, are expected to have a narrower window of suitable sequence divergences than phylogeny reconstruction. The large-sample 2 approximation to the likelihood-ratio test statistic might also be examined, but limited simulations suggest that typical sequence data (with Ͼ100 codons) are large enough for it to be reliable. For very short genes or gene regions and especially at low sequence divergences, Monte Carlo simulation might be needed to derive the null distribution. The likelihood analysis assumes no recombination within a gene. If recombination occurs, different regions will have different phylogenies. Empirical data analysis suggests that the phylogeny does not have much impact on tests of positive selection and identification of sites, and one might suspect that recombination will not cause false positives by the likelihood-ratio test. However, simulation studies are necessary to understand whether this is the case. T he evolutionary origin of snakes (or Serpentes) has been discussed for over 130 years and their phylogenetic position within squamates is still debated. Around 2700 snake species are alive today and these are divided into three main groups 1-3 (Box 1): tiny fossorial (burrowing) scolecophidians (blindsnakes); anilioids (pipesnakes), which are mostly semi-fossorial; and macrostomatans, which include more familiar taxa, such as boas, pythons, vipers and cobras. In addition to the more obvious diagnostic characters of body elongation, limblessness and jaws that can engulf surprisingly large prey, other key features of snakes include absence of eyelids and external ears, and the presence of deeply forked tongues (linked to their highly attuned and sophisticated chemosensory systems 3 ). Hypotheses concerning snake interrelationships fall into two main groups. For some researchers, snakes descend from terrestrial squamates that developed fossorial (burrowing) habits. Two groups of lizards exhibiting such habitats, amphisbaenians and dibamids, have often been regarded as snakes' closest living relatives 4 . Amphisbaenians, in particular, resemble scaly, loose-skinned earthworms, whose shovel-shaped or wedge-like heads function as soil-shunting devices. Specializations shared by snakes (Fig. 1a) , amphisbaenians (Fig. 1b) and dibamids include loss, reduction and consolidation of skull bones; braincase enclosure; dorsal displacement of jaw-closing muscles; loss or reduction of limbs and girdles; and increased uniformity along the vertebral column. Furthermore, differences between the eyes of lizards and snakes are consistent with a model in which structures that were barely useful in a burrower underwent progressive reduction. Thus, whereas lizards, like humans, distort eye lens curvature to focus on objects, snakes lack ciliary muscles and are compelled to move the entire lens back and forth relative to the retina. Moreover, unlike lizards, snakes lack both a fovea and coloured oil droplets in retinal cells 1 . Alternative hypotheses 5 postulate that snakes are related to mosasauroids (Fig. 1c) : spectacular marine reptiles from the upper half of the Cretaceous period, some 65-100 Mya 6 . Mosasauroids and snakes share reduced ossification of the pelvis and hindlimbs as well as specialized features of the jaw suspension and intramandibular joint kinetics (presence of a hinge allowing a degree of lateral movement within the lower jaw; Fig. 1a,c,d; Fig. 2, red circle) . Phylogenetically, mosasauroids would be the nearest monophyletic sister group of snakes, with varanoid lizards (monitors) as the immediate sister group to this pair. Given this theory of relationships, the latest common ancestor of mosasaurs and snakes has been argued to have been a limbed, aquatic or semiaquatic squamate 5,7-11 . Note that the implied ecological shift from an aquatic to a terrestrial environment in snake ancestry suggests that mosasaurs' (implied) aquatic habits were also primitive for Serpentes. Subsequently, snakes reduced and lost their limbs, although rudiments of the posterior pair remain in some forms, such as pythons. Renewed interest in the origin of snakes has been triggered by the recognition and discovery of three remarkable fossil forms with hind legs. Each of these ancient snakes is around 97 My old and originates from lowermost Upper Cretaceous sediments in the Middle East. Pachyrhachis problematicus, from Israel ( Fig. 1d-f ), rapidly assumed a central position in debates about snake phylogeny 6,7,12 . It has miniature hindlimbs articulated with a rudimentary pelvic girdle (Fig. 1e,f) , but sadly, its feet are missing. Currently described from only two specimens, it Nice snake, shame about the legs Snakes are one of the most extraordinary groups of terrestrial vertebrates, with numerous specializations distinguishing them from other squamates (lizards and their allies). Their musculoskeletal system allows creeping, burrowing, swimming and even gliding, and their predatory habits are aided by chemo-and thermoreceptors, an extraordinary degree of cranial kinesis and, sometimes, powerful venoms. Recent discoveries of indisputable early fossil snakes with posterior legs are generating intense debate about the evolutionary origin of these reptiles. New cladistic analyses dispute the precise significance and phylogenetic placement of these fossils. These conflicting hypotheses imply radically different scenarios of snake origins and relationships with wide biological implications. Rates of conservative and radical nonsynonymous nucleotide substitutions in mammalian nuclear genes Positive selection for colicin diversity in bacteria Coordinated amino acid changes in the evolution of mammalian defensins Molecular phylogeney of Fv1 Positive darwinian selection observed at the variable-region genes of immunoglobulins Pattern of nucleotide substitution at major histocompatibility complex class I loci reveals overdominant selection Identification of target amino acids that affect interactions of fungal polygalacturonases and their plant inhibitors Conserved evolution of the Rh50 gene compared to its homologous Rh blood group gene Positive darwinian selection after gene duplication in primate ribonuclease genes Natural selection promotes divergence of transferrin among salmonid species The evolution of the type I interferon family in mammals Patterns of divergence during evolution of a 1 -proteinase inhibitors in mammals Natural selection on Plasmodium surface proteins Recombination of hepatitis D virus RNA sequences and its implications Positive selection and interallelic recombination at the merozoite surface antigen-1 (MSA-1) locus of Plasmodium falciparum Sequence evolution of the porB gene of Neisseria gonorrhoeae and Neisseria meningitidis: evidence for positive darwinian selection Episodic evolution mediates interspecific transfer of a murine coronavirus Positive darwinian selection on two homologous fertilization proteins: what is the selective pressure driving their divergence? Positive selection and the molecular evolution of a gene of male reproduction, Acp26Aa of Drosophila Reduced nucleotide variability at an androgen-binding protein locus (Abpa) in house mice: evidence for positive natural selection Positive selection and sequence arrangements generate extensive polymorphism in the gamete recognition protein bindin A rapidly evolving homeobox at the site of a hybrid sterility gene Rapid evolution of a homeodomain: evidence for positive selection Rapid evolution of a primate sperm protein: relaxation of functional constraint or positive darwinian selection? Identification of regions in which positive selection may operate in S-RNase of Rosaceae: implications for Sallele-specific recognition sites in S-RNase Evolution of Sry genes Nucleotide sequence evolution at the k-casein locus: evidence for positive selection within the family Bovidae Molecular genetics of ecological diversification: duplication and rapid evolution of toxin genes of the venomous gastropod Conus Accelerated evolution in the protein-coding regions is universal in crotalinae snake venom gland phospholipase A 2 isozyme genes Molecular evolution of the COX7A gene family in primates Molecular evolution of cytochrome c oxidase subunit IV: evidence for positive selection in simian primates Evolution of hemopoietic ligands and their receptors: influence of positive selection on correlated replacements throughout ligand and receptor proteins The molecular evolution of vertebrate growth hormones: a pattern of near-stasis interrupted by sustained bursts of rapid changes Antarctic fish hemoglobins: evidence for adaptive evolution at subzero temperatures Natural selection and the origin of jingwei, a chimeric processed functional gene in Drosophila A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to the chloroplast genome We thank D. Haydon, J. Mallet, T. Ohta, A. Pomiankowski, V. Vacquier, W. Swanson and three anonymous referees for comments. We also thank several users of the PAML package (http://abacus.gene.ucl.ac.uk/software/paml.html), in particular C. Woelk, for comments and suggestions concerning the implementation. This work is supported by grant #31/G10434 from the Biotechnology and Biological Sciences Research Council (UK).