Statistical methods for detecting molecular adaptation

2. Measuring selection using the nonsynonymous/synonymous (d_N/d_S) rate ratio

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 ⁵. 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 ⁶), 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 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.

3. Estimation of d_N and d_S between two sequences

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 5., 9., 10., 11., 12., 13., 14., 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 ⁵, and Nei and Gojobori ⁹, 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 ¹⁰. Efforts have been taken to incorporate the transition/transversion rate bias (Box 1) when counting sites and differences 10., 11., 12., 13., 14.. The effect of biased codon usage has largely been ignored ¹⁷; however, extreme codon-usage bias can have devastating effects on the stimation of d _N and d _S (see the next section) 15., 18.. A recent ad hoc method ¹⁵ 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 k 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 k; correcting for multiple hits; and weighting pathways of change between codons.

Statistical tests can be used to test whether d _N is signifi-cantly higher than d _S. For approximate methods, a normal approximation is applied to d _N–d _S. For ML, a likelihood-ratio 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 χ² 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 α²-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 ¹⁵ 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 ¹⁹.

• •

  Ignoring the transition/transversion rate bias leads to underestimation of S, overestimation of d _S and underestimation of the ω ratio ¹⁰.

• •

  Codon-usage bias in these data has the opposite 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 ¹⁵, 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 ²¹.

• •

  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.

4. Detecting lineage-specific episodes of darwinian selection

If, for most of the time, a gene evolves under purifying selection but is occasionally subject to episodes of adaptive change ²², 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 ²³ 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 ²⁶.

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 mode 12., 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 ³⁰, 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 ²⁷.

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

5. Detecting amino acid sites under darwinian selection

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 ³⁵ 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 ²⁴, 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 20., 36., 37..

Nielsen and Yang ³⁶ 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. ³⁸ used this test to identify several sites under strong positive selection in the nef gene of HIV, whereas both pairwise comparison and sliding-window 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 ³⁷. Thus, Yang et al. ³⁷ 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 ³⁷. These models identified positive selection in six out of ten genes the authors analysed. Fig. 1 shows the use of a discrete model (M3) with three classes to identify sites under diversifying selection in abalone sperm lysin ³⁹.

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 ³⁸ and env ⁴⁰ 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 ²⁰.

6. Limitations of current methods and future directions

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 ex-cessive 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 χ² 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.