THEORY

Model of codon substitution: A simplified version of the codon substitution model proposed by Goldman and Yang (1994) will be used in this paper. Stop codons are not allowed in the gene sequence, and substitutions between sense codons are described by a continuous time Markov process. The instantaneous substitution rate from codon i to j (i ≠ j) is given by Qij={0,if the two codons differ at more than one position,πj,for synonymous transversion,κπj,for synonymous transition,ωπj,for nonsynonymous transversion,ωκπj,for nonsynonymous transition,} (1) where parameter κ is the transition/transversion rate ratio, ω is the nonsynonymous/synonymous rate ratio, and π_j is the equilibrium frequency of codon j, calculated using the nucleotide frequencies at the three codon positions. Note that the relationship holds that ω = d_N/d_S. The unit of evolution under such a model is the codon and, in this paper, the term “site” refers to a codon or amino acid instead of a nucleotide.

Neutral model: This model assumes two categories of sites in the gene. The first category includes neutral sites where nonsynonymous mutations are neutral, with the d_N/d_S ratio ω₁ = 1. The substitution rate from codon i to j for a site in this category is Qij(1)={0,if the two codons differ at more than one position,πj,for transversion,κπj,for transition.} (2) The second category includes the conserved sites, where nonsynonymous mutations are deleterious and eliminated by selection, and only synonymous substitutions are possible, so that ω₂ = 0. The substitution rate from codon i to j (i ≠ j) at such a site is thus Qij(2)={πj,for synonymous transversion,κπj,for synonymoums transition,0,otherwise.} (3) Under this model, the synonymous substitution rate is constant among sites, while the nonsynonymous rate is variable. At the conserved sites, the nonsynonymous rate is zero, and at the neutral sites, the nonsynonymous rate is equal to the synonymous rate.

Let the proportions of codon sites in the two categories be p₁ and p₂ = 1 − p₁. Let n be the number of sites (codons) in the sequence and the data at site h be x_h (h = 1, 2, …, n). Because we do not know a priori which category a site belongs to, the probability of observing data x_h is an average over the two possibilities: f(xh)=∑k=12pkf(xh∣ωk)=p1f(xh∣ω1)+p2f(xh∣ω2), (4) where f(x_h|ω_k) is the probability of observing data x_h given that site h is from category k (k = 1, 2), with nonsynonymous/synonymous rate ratio ω_k. This conditional probability can be calculated for a given phylogenetic tree and branch lengths according to the method of Goldman and Yang (1994). The log likelihood is a sum over all n sites in the sequence 1=∑k=1nlog{f(xh)}. (5)

It may be noted that the structure of the model is similar to models of variable evolutionary rates among nucleotide or amino acid sites developed previously (Yang 1993, 1994, 1995). Evolutionary rates under these variable rates models, and the d_N/d_S ratios (ω) reflecting selection intensity under the present model, are both formulated as random variables and integrated out in the likelihood function. Calculation of the likelihood function under the codon-based model can be easily adapted from the algorithm for variable evolutionary rates among sites (Yang 1994). The matrix of transition probabilities over time t, P(t) = e^Qt, is calculated through diagonalization of the rate matrix Q = {Q_ij}; a standard numerical algorithm is used for this purpose (Goldman and Yang 1994). P(t) needs to be calculated separately for each branch (with length t) and site category (with parameter ω). Parameters in the model include the branch lengths (t), the transition/transversion rate ratio (κ), and the proportions of the two categories of sites (p₁ and p₂, with p₁ + p₂ = 1). These parameters are estimated by maximum likelihood using numerical optimization algorithms. The codon frequency parameters (π_j) are calculated using the observed nucleotide frequencies at the three codon positions.

Positive selection model: The neutral model can be extended by adding an extra category of positively selected sites (with ω₃ > 1). Nonsynonymous mutations at such sites thus have higher probabilities of fixation than synonymous mutations. Let the proportions of sites in the three categories be p₁, p₂, and p₃ (with p₁ + p₂ + p₃ = 1), and let the corresponding nonsynonymous/synonymous rate ratios be ω₁ = 1, ω₂ = 0, and ω₃ > 1. The probability of observing data (x_h) at site h is then f(xh)=∑k=13pkf(xh∣ωk) (6) The log likelihood function can be calculated similarly to that under the neutral model with two categories of sites.

The positive-selection model is an extension of the neutral model, with two more parameters. Twice the log likelihood difference between the two models can be compared with a χ² distribution with d.f. = 2. This constitutes a likelihood ratio test of neutrality against an alternative model of positive selection.

In practice, ω₃ is optimized in the entire region from zero to infinity. In this case, it is appropriate to call the model a positive-selection model only if ω₃ > 1, as an estimate of ω₃ smaller than one provides no evidence for positive selection.

Because one might expect a continuum of the d_N/d_S rate ratio among sites even when no positive selection operates, we consider a variation of the neutral model in which the class ω₀ = 1 is replaced by a truncated gamma distribution on (0, 1). A proportion p₁ of sites have rates from the truncated gamma distribution, while a proportion p₂ (= 1 − p₁) of sites are conserved and have ω₂ = 0. The truncated gamma distribution with shape parameter α and scale parameter β has the following density f(w)=e−βwwα−1∫01e−βwwα−1dw,0<w≤1. (7)

To avoid the use of too many parameters, we fix β to be equal to α. Use of other values for β is found to give roughly the same likelihood. For efficient computation, the truncated gamma distribution was approximated by five rate categories with equal probabilities (see Yang 1994). We will refer to this model as the “continuous neutral model.” In the corresponding positive-selection model, we add an extra category of sites with proportion p₃, and d_N/d_S ratio ω₃. This will be referred to as the “continuous positive-selection model.” Notice that when α = ∞, the two new models with continuous substitution rates are identical to the neutral and positive-selection models described above.

Identification of positively selected amino acid sites: When parameters of the positive-selection model are estimated, an empirical Bayes' approach can be used to infer which category the site most likely belongs to. This method is similar to the approach of Yang (1994) for estimating rate categories for nucleotide sites under a model of variable evolutionary rates among sites. The posterior probability that a site with data x_h belongs to category k (with rate ratio ω_k) is given by prob(ωk∣xh)=pkf(xh∣ωk)f(xh)=pkf(xh∣ωk)Σi=13p1f(xh∣ωI). (8) The category k that maximizes the posterior probability is the most likely category for the site. Positively selected sites (i.e., sites belonging to the third category with ω₃ > 1) may be identified this way. The posterior probabilities provide a measure of accuracy of that inference.

APPLICATIONS TO THE HIV-1 ENVELOPE GENE

To examine the utilities of the models developed in this paper, we analyze the DNA sequence data of the HIV-1 envelope genes published by Holmes et al. (1992). The sequences are from viral variants in years 3 through 7 after infection of one single patient. No sequence variation was detected in year 0, and no data were available for years 1 and 2 (Holmeset al. 1992). Each sequence contains 77 codon sites from the third hypervariable region (V3) and flanking regions of the external glycoprotein (gp120). There has been considerable interest in possible selective factors acting on this gene because it encodes known targets for cytotoxic T lymphocytes and neutralizing antibodies. For example, the V3 region is suggested to be under positive diversifying selection during the course of infection (Boenhoefferet al. 1995). Boenhoeffer et al. found that the average d_N/d_S rate ratios for each year, estimated using the method of Nei and Gojobori (1986) for pairwise sequence comparison, were greater than one and decreased over the years, suggesting that positive selection may be important in the diversification of the HIV-1 envelope genes and that the selection pressure may have been relaxed during the cause of the viral infection.

There are 15, 11, 23, 15, and 13 sequences for years 3, 4, 5, 6, and 7, respectively, and the number of distinct sequences is 13, 11, 17, 15, and 12, respectively. As our codon-based models involve heavy computation, which makes it unfeasible to perform an analysis on a phylogeny of all the 77 sequences simultaneously, we analyze sequences from different years as if they were separate data sets. Only distinct sequences are used. For data of each year, the likelihood method was used under both simple and sophisticated nucleotide substitution models to perform heuristic tree searches to identify candidate trees. Candidate topologies found in this way often share long interior branches, for which the statistical support is strong, while the details of the topology may be different in different analyses. Results (not shown) suggest that our codon-based analysis, to be presented below, is not sensitive to the assumed topology of the phylogenetic tree, as long as the long interior branches are preserved in the topology. Two candidate topologies for each data set are used in later analysis. They give essentially identical parameter estimates, and results obtained from only one of them are reported in this paper. The tree topologies used are not presented, but are available from the authors upon request. For each data set, several starting values were used in the iteration. This was done as a protection against the existence of multiple local optima in the likelihood function. In this study, all starting values for a given data set resulted in the same optimum.

Besides the neutral and positive-selection models of this paper, two additional codon-based models developed previously (Goldman and Yang 1994) are also used in the analysis. Both of these models assume a constant d_N/d_S ratio among sites. The first, referred to as the “Goldman and Yang” model in Table 1, applies the same rate matrix (Q specified by Equation 1) to all sites. The second, referred to as the “gamma-rates” model in Table 1, assumes gamma-distributed rates among codon sites (Yang 1994; Goldman and Yang 1994). In this model, the substitution rate from codon i to j at a codon site with rate r is rQ_ij, with Q_ij given by Equation 1. A codon with a high rate tends to have both high synonymous and high nonsynonymous rates. The shape parameter α of the gamma distribution is inversely related to the extent of rate variation.

Comparison of models and tests of positive selection: The log likelihood values under different codon-based models are listed in Table 1 for data of each year. The gamma model reduces to the Goldman and Yang model when the gamma shape parameter α → ∞. Comparison of twice the log likelihood difference between the two models with a χ² distribution with 1 d.f. suggests that the gamma model provides a significantly better fit to data of each year. This result is in accordance with previous analyses (Goldman and Yang 1994), demonstrating that rate variation is an important factor in the evolution of coding sequences.

The neutral model and the Goldman and Yang model are not nested, and so cannot be tested using a χ² approximation. Nevertheless, their likelihood values are comparable. The two models have the same number of parameters but the log likelihood value under the neutral model is higher than that under the Goldman and Yang model by 5–12 units (Table 1). The neutral model is thus a much more realistic representation of the evolutionary dynamics of the HIV-1 envelope gene. This is also true for year 3, even though the average rate of nonsynonymous substitution is much higher than the average rate of synonymous substitution, suggesting that the neutral model should provide a poor fit to the data. Except for data of year 3, the fit of the neutral model is almost as good as, or even slightly better than, the gamma-rates model.

The positive-selection model includes two more parameters than the neutral model. These two models are nested, and twice the log likelihood difference can be compared with a χ² distribution with d.f. = 2 to test whether the positive-selection model provides a better fit to data than the neutral model. The difference is significant (P < 5%) for data of every year, except year 7 (for which P ≈ 9%). Clearly, the neutral model is inadequate to describe these data, and positive selection has been operating in the evolution of this viral gene. The Goldman and Yang model is also a special case of the positive-selection model with the constraint that p₁ = p₂ = 0, and p₃ = 1. Twice the log likelihood difference between the two models ranges from 24 to 37 for different years. These differences are all significant (P < 1% with d.f. = 2), and the Goldman and Yang model with a constant d_N/d_S ratio among sites is rejected in favor of the positive-selection model for data of each year.

The same conclusion is reached when rates in the range between 0 and 1 are allowed in the neutral model (continuous neutral model). Indeed, the neutral model with, and without, continuously distributed mutation rates, modeled by the truncated gamma distribution, have the same likelihood values. Similarly, the positive-selection model, with and without continuous mutation rates, have very similar likelihood values. We have also fitted another set of neutral and positive-selection models by removing the class with ω₂ = 0 in the continuous mutation models. The neutral model thus constructed, which assumes that all sites have the d_N/d_S rate ratio ω from the truncated gamma distribution, fits the data much more poorly than the neutral model specified by Equations 2 and 3. The corresponding positive-selection model produces likelihood values very similar to those under the simple three-class positive-selection model. Therefore, likelihood ratio tests using this set of models all suggest positive selection acting on the gene (results not shown). To summarize, the statistical support for positive selection on the envelope gene appears to be rather insensitive to the assumed distribution of selection intensity among sites.

The positive-selection model and the gamma-rates model are not nested and cannot be compared using a χ² approximation. However, the positive selection model has higher likelihood value than the gamma-rates model for data of each year. Although the significance of the differences is uncertain, the results suggest that the positive-selection model provides the most realistic description of the evolution of the analyzed sequences. The reason appears to be that in the gamma-rates model the same distribution of rates is applied to both the synonymous and nonsynonymous substitutions. However, most of the rate variation appears to be caused by variation in the selection intensity among nonsynonymous substitutions.

Variation of nonsynonymous substitution rates among years: Maximum likelihood estimates of parameters obtained under different codon-based models are listed in Table 1 for data of each year. Estimates of the transition/transversion rate ratio (κ) are more variable among years than among methods, and range from two to six, indicating that transitions occur much more frequently than transversions.

Notice that the neutral model, with and without a truncated gamma distribution, provides identical parameter estimates. The reason for this may be that the true rates are highly bimodal with a mode at ω = 0 and a mode at ω > 1. Very little of the probability mass appears to be located in the region between zero and one. Likewise, for the selective models, almost no improvement in the likelihood is obtained by allowing nonsynonymous rates in the region between zero and one. In the following we will therefore concentrate on the results of the neutral and selection models that do not allow nonsynonymous rates between ω = 0 and ω = 1.

Considerable differences exist in estimates of ω between year 3 and the remaining years. For years 4 through 7 under the Goldman and Yang model, the estimate of ω is <1. Likewise, estimates of p₂ (proportion of conserved sites with ω₂ = 0) under both the neutral and the positive-selection models range from 0.55 to 0.67 among years 4 through 7, suggesting that a majority of amino acid sites in the protein are conserved in years 4 through 7. Estimates of p₃ (proportion of positively selected sites with ω₃ > 1) under the positive-selection model range from 0.04 to 0.22, indicating that only a few sites in the sequence are under positive diversifying selection at any particular time.

Boenhoeffer et al. (1995) calculated the average d_N/d_S ratio among pairwise comparisons of sequences in the same year, and noted a decrease in the average d_N/d_S ratio over the years. They suggested that relaxed selection pressure during the course of infection may be responsible. Nielsen (1997) suggested that the change in the inferred nonsynonymous/synonymous rate ratios using a model of constant selection intensity may partly be due to the presence of hypervariable nonsynonymous sites. Ignoring nonsynonymous rate variation leads to underestimation of the nonsynonymous/synonymous rate ratio at high levels of sequence divergence. The results of the present study show that hyper-variable nonsynonymous sites are in fact present in the genes. However, even after the nonsynonymous rate variation has been accounted for, a change in the selection intensity is still observed between year 3 and the subsequent years. Estimates for year 3 under the positive-selection model suggest a high proportion (p₃ = 21.6%) of sites under strong positive selection (ω₃ = 29.5). Estimates of both p₃ and ω₃ for years 4 through 7 are much smaller, indicating that the effect of positive selection is stronger in year 3 than in the later years. These results are in general agreement with the conclusion of Boenhoeffer et al. (1995), although the standard errors of our parameter estimates (not shown) are large.

Identification of positively selected sites: Because the positive-selection model provides a better fit to data of every year than the neutral model, Equation 7 is used to infer the most likely site category (with the associated d_N/d_S ratio) at each codon (amino acid) site. The posterior probabilities are also calculated. The results are shown in Figure 1, A–E, for years 3 through 7. Consistent with the parameter estimates (Table 1), more sites are inferred to be under positive selection in year 3 than in years 4 through 7. Only site 77 was identified with very high posterior probability to be under positive selection in all five years. Site 32 was identified to be under positive selection in years 3, 5, 6, and 7, while site 44 may be under positive selection in years 3, 4, 5, and 7. Apart from these sites, there seems to be considerable variation in the inferred sites for positive selection from year to year. The variation over the years may in part be due to rapid fixations of old mutations, and arrivals of new ones, as Holmes et al. (1992) observed drastic changes of the most abundant sequence over the years.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

(A–E) Estimates of nonsynonymous/synonymous rate ratios (ω in Equation 1) for sites in the V3 region of the HIV-1 envelope gene. Positively selected sites are labeled black if the posterior probability is larger than 95%, and gray if otherwise. Note that the scales are different among years. (F) Estimates of amino acid replacement rates under the empirical model of Jones et al. (1992) with gamma distributed rates among sites (Yang 1994). These are relative rates with the mean to be one. The location of the V3 region is indicated.

Our analysis suggests that there may be as many positively selected sites (sites potentially under positive selection) in the flanking regions as in the V3 region. The single site that is identified to be under positive selection in all 5 years is site 77, outside the V3 region. The results suggest that selection may be acting on a broader region of the gene sequence than previously suggested (e.g., Willeyet al. 1986). Furthermore, a considerable amount of the variability in the V3 region appears to be neutral.

It should be noted that positively selected sites identified in this way are not equivalent to highly variable sites. A positively selected site is characterized by both a high variability and an excess of nonsynonymous to synonymous substitutions. To illustrate this point, evolutionary rates at amino acid sites were estimated using the protein sequences translated from the DNA sequences of years 6 and 7 (Figure 1F). The empirical model of Jones et al. (1992) with gamma-distributed rates among sites was assumed. The method of Yang (1994, 1995) was then used to estimate the amino acid replacement rates at sites. Notice that there are many more sites identified as having high rates than there are sites undergoing positive selection. Our method, which assumes variable selection intensity instead of variable amino acid replacement rates, should therefore have a higher precision in identifying target sites for selection.

Analysis of Cytochrome Oxidase II: a case of no positive selection: To illustrate the performance of the method in a case where positive selection is not expected to have strongly influenced the evolution of the DNA sequence, we also applied the method to 10 vertebrate sequences of the mitochondrial Cytochrome Oxidase II (COII) gene. Sequences from 10 vertebrate species are used (the data set is described in Cummingset al. 1995). This COII gene has been extensively analyzed in the literature, and to the knowledge of the authors, there have been no suggestions of codons with elevated nonsynonymous rates. The results of the analysis are in Table 2. Note that an estimate of ω₃ < 1 was obtained for this data set. This suggests that there may not be a large category of amino acid sites in this gene, where nonsynonymous substitutions occur at a higher rate than synonymous substitutions. This result is in contrast to those obtained from the HIV-1 envelope gene.