RESULTS

Nucleotide (codon) frequencies in Drosophila nuclear genes: GC3 varied substantially among the 83 genes, ranging from 28 to 93%. This variation was characteristic of each dataset (DmDpDsub, 44-90%; DmDsimDy, 52-93%; DmDp, 36-90%; DmDv, 45-88%; and DmDsim, 28-91%). Chi-square tests suggested significant heterogeneity in nucleotide composition among lineages for 4 of 20 genes in DmDpDsub (Gapdh2, Gpdh, RpII215, and ry), 8 of 35 genes in DmDp [Gapdh2, Gpdh, l(2)gl, Rh2, RpII215, ry, Tl, and Ubx], and 11 of 39 genes in DmDv [Adh, Amy-p, Cdc37, fu, gbb, Gpdh, Kr, l(2)tid, lama, nos, and Sry-beta].

Consistent with patterns of nucleotide bias, codon usage varied greatly among genes; the ENC values ranged from extreme bias at 26.7 (GstD1) to no bias at 61 (sc). Each dataset exhibited substantial variation in codon bias among genes (DmDpDsub, 27.3-59.8; DmDsimDy, 26.7-61; DmDp, 28.6-56.6; DmDv, 29.2-61; and DmDsim, 26.7-61). Codon bias was not related to heterogeneity in nucleotide composition among lineages. For example, genes with homogeneous nucleotide composition among lineages exhibited both highly biased (e.g., gene Eno, ENC = 28.6) and unbiased (e.g., gene ac, ENC = 56.1) codon usage.

Lineage-specific patterns of synonymous and nonsynonymous substitution: ML estimation of d_N and d_S was carried out for each gene in the DmDpDsub and DmDsimDy datasets (Table 1). Constancy of nonsynonymous/synonymous rate ratios (d_N/d_S) was tested, and d_N/d_S ratios were largely homogeneous over lineages (Table 1). Homogeneity was rejected in only 3 of the 20 genes in DmDpDsub and 3 of the 11 genes in DmDsimDy. This trend differs from that observed in mammals, where over half of genes surveyed exhibited significant variation in selective constraints among lineages (Yang and Nielsen 1998; Bielawskiet al. 2000).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

—The relationship between d_N and d_S (a) and between d_S and ENC (b) for Drosophila pseudoobscura. Estimates of d_S and d_N were obtained from lineage-specific analyses.

Genes exhibiting high synonymous or nonsynonymous rates in one lineage tended to exhibit high rates in other lineages as well. For example, the coefficient of determination (r²) of d_S estimates between Dm and Dp is r² = 0.5098 (P = 0.0004), and the coefficient of determination of d_N is r² = 0.5005 (P = 0.0005). Similar patterns were reported for mammalian genes (Bulmeret al. 1991; Mouchiroudet al. 1995). Substantial variation in d_N/d_S was found among genes, as they are under very different selective constraints (Table 1).

The relationship between synonymous and nonsynonymous substitution rates was evaluated by linear correlation using rate estimates of Table 1 for five species. Estimates of d_S were positively correlated with d_N in every lineage (Dp, r² = 0.4643, P = 0.0009; Dsub, r² = 0.7312, P ⪡ 0.0001; Dsim, r² = 0.5722, P = 0.0071; Dy, r² = 0.3945, P = 0.0385; Dm, r² = 0.4587, P ⪡ 0.0001). The plot for D. pseudoobscura is presented as an example (Figure 1a). These results are consistent with a previous analysis of substitution rates in Drosophila (Comeron and Kreitman 1998). However, this study represents a wider sample of Drosophila species and the correlations reported here are stronger than previously reported.

Synonymous substitution rates were independent of codon bias, as evaluated by linear regression of d_S and ENC, in every lineage (Dp, r² = 0.0004, P = 0.7911; Dsub, r² = 0.0138, P = 0.6218; Dsim, r² = 0.151, P = 0.2376; Dy, r² = 0.0045, P = 0.8446; Dm, r² = 0.0079, P = 0.6405). The plot for D. pseudoobscura is presented as an example (Figure 1b). All previous analyses differ from ours in suggesting synonymous substitution rates are negatively correlated with codon bias (Shieldset al. 1988; Sharp and Li 1989; Moriyama and Hartl 1993; Powell and Moriyama 1997). Previous studies, however, also differ from ours in using approximate estimation methods, which ignore codon usage bias. Furthermore, in this study we estimated d_S per lineage rather than between pairs of sequences. In the next section, we explore what factors may be responsible for this conflict.

Reconciling differences between ML and approximate methods: Models employed in ML estimation of substitution rates explicitly assumed that nucleotide/codon frequencies were at equilibrium (Goldman and Yang 1994). To see whether violation of this assumption influenced our correlation analysis, we reevaluated the relationship between d_S and ENC using only stationary genes. To keep more genes in the dataset and to facilitate direct comparison with approximate methods, we apply the ML method to two-taxa datasets. This pairwise ML analysis led to the same conclusion as the phylogeny-based ML analyses; i.e., synonymous rates were independent of codon bias (Figure 2). However, when either NG or Comeron’s (1995) method was used, d_S was negatively correlated with ENC for both the DmDv (NG, r² = 0.5036, P ⪡ 0.0001; Comeron, r² = 0.4615, P = 0.0001) and DmDp (NG, r² = 0.6754, P ⪡ 0.0001; Comeron, r² = 0.6127, P ⪡ 0.0001) datasets (Figure 2), a result consistent with previous studies (Sharp and Li 1989; Comeron and Aguadé 1996; Powell and Moriyama 1997). However, in the DmDsim dataset, the approximate methods lead to the same conclusion as the ML method; d_S was independent of ENC (NG, r² = 0.0995, P = 0.1332; Comeron, r² = 0.0583, P = 0.2557; Figure 2). Because results were similar for the two approximate methods, only results for NG are presented in Figure 2. On the basis of these results we conclude that the conflict between approximate methods and ML cannot be attributed to use of a pairwise approach, the sample of genes, or nonstationarity of nucleotide frequencies.

To examine the effect of transition/transversion bias and codon bias on estimation of d_S, we changed parameters of the codon model to reanalyze the two-taxa datasets by ML. The effect of ignoring transition/transversion bias was evaluated by setting κ = 1, but allowing for biased codon usage. Results obtained using this codon model did not differ from those obtained using the full codon model; i.e., there was no significant relationship between d_S and ENC for DmDv (r² = 0.140, P = 0.060) and DmDp (r² = 0.051, P = 0.256). The effect of ignoring codon bias was evaluated by assuming equal codon frequencies, but allowing for transition/transversion bias. This codon model produced results different from those obtained previously and matched the approximate methods, with a significant correlation between d_S and ENC in DmDv (r² = 0.435, P = 0.0005) and DmDp (r² = 0.6281, P ⪡ 0.0001). These findings suggest that the different treatment of codon bias was responsible for the conflict between ML and approximate methods in the relationship between d_S and ENC.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

—The relationship between pairwise estimates of d_S and mean ENC (a-f). Pairwise estimates of d_S were computed using both ML (Goldman and Yang 1994) and NG (Nei and Gojobori 1986). DmDv indicates pairwise comparisons between D. melanogaster and D. virilis; DmDp, pairwise comparisons between D. melanogaster and D. pseudoobscura; and DmDsim, pairwise comparisons between D. melanogaster and D. simulans.

However, in one pairwise comparison (DmDsim), approximate and ML methods were in agreement, yet codon usage in this dataset (mean ENC = 43) was just as biased as in the other two sets of genes (DmDv, mean ENC = 44; DmDp, mean ENC = 41) where approximate and ML methods differed. Hence, there must have been an additional factor. In both the DmDv and DmDp datasets, where approximate methods produced a positive relationship between d_S and ENC, uncorrected percent sequence divergences (p) were very large (DmDv, mean p = 21.6 ± 6.0%; DmDp, mean p = 16.3 ± 5.8%). In the DmDsim dataset, where d_S was independent of ENC, divergences were very low (mean p = 3.2 ± 1.4%). This pattern suggested a possible “saturation” effect. Consequently, we hypothesized that differences between approximate and ML methods might be related to the combined effects of sequence divergence and codon usage bias.

Simulation studies: Simulations were performed to evaluate the hypothesis that improper treatment of both divergent sequences and codon usage bias by approximate methods could have led to a significant positive correlation between d_S and ENC when, in reality, they were independent. Codon sequences were simulated under values of t chosen to reflect the range exhibited among Drosophila genes, and codon frequencies were modeled after observed frequencies of Drosophila genes (see materials and methods). Simulated codon sequences were analyzed using both NG and ML (F3 × 4), and these results were compared with the true values employed in the simulation (Figure 3).

Both sequence divergence and codon bias had an effect on estimates of d_S (Figure 3). Only at low sequence divergences (t < 1.6) and relatively unbiased codon usage (ENC = 53) did both NG and ML (F3 × 4) produce unbiased estimation of d_S. With one exception (Figure 3a), d_S was underestimated by both methods with increasing levels of sequence divergence (Figure 3, b-h), and the NG method involved a much more serious estimation bias than ML. Because equilibrium codon frequencies obtained using the F3 × 4 model (9 frequency parameters) do not perfectly match the empirical Drosophila frequencies, some degree of error was expected for ML. Use of the more parameter-rich model (60 frequency parameters) produced estimates of d_S essentially identical to the true values (data not shown). Differences between the analysis of real data using the NG and ML (F3 × 4) methods are consistent with differences observed in our simulation study; i.e., ML estimates of d_S were larger than estimates obtained using NG. These simulations further suggested that the ML (F3 × 4) estimates of d_S for Drosophila genes are, themselves, likely to be underestimates of the true values. Although the F3 × 4 model is biased in cases of extreme codon bias, this model was acceptable over a wide range of codon biases and consistently outperformed the NG method.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

—Estimates of d_S in simulated data by NG (○) and ML (F3 × 4; ▵) plotted against sequence divergence t (the expected number of nucleotide substitutions per codon). Data for a pair of sequences, each of 1 million codons, were simulated for different codon bias measured by ENC.

Plots of d_S (Figure 3) are reminiscent of the “saturation effect” on plots of uncorrected sequence divergence. For d_S, however, the “ceiling” appears to be related to levels of codon usage bias; the effect is relatively minor when codon usage is unbiased and extreme when codon usage is highly biased (Figure 3). For example, in the most extreme case of codon bias (ENC = 28), the estimate of d_S obtained using the NG method peaked at 0.5 (t > 0.8), although true values of d_S range from 0.95 when t = 0.4 to 6.6 when t = 2.8. The only exception to the general pattern occurred when codon usage was relatively unbiased (Figure 3a). In this case, when 1.2 < t < 2.0, NG overestimated rather than underestimated d_S, and when t > 2.0 the method yielded invalid estimates of d_S. Interestingly, Muse (1996) suggested that even when there is no codon usage bias, ad hoc corrections for multiple substitutions applied in the approximate methods could introduce substantial bias.

The dramatic effect of codon usage bias on estimation of d_S is caused mainly by its effect on counting of synonymous (S) and nonsynonymous (N) sites. When codon usage is biased but ignored (Figure 3, b-h), the number of synonymous sites (S) was overestimated, leading to an underestimation of d_S. For example, in the most extreme case (ENC = 28, t = 2.8), most codons end with C or G and most changes at the third position are transversions between C and G, which are more likely to be nonsynonymous than random changes between nucleotides. As a result, the proportion of synonymous sites is as low as S = 8.53%. The NG method assumes unbiased codon usage or equal rates of change between nucleotides and expects more frequent transitional or synonymous changes, giving S = 23.6%, with almost a fourfold difference. When there is little codon bias (ENC = 53), NG reliably estimated S (NG, S = 23.4%; true value, S = 23.3%). Thus the underestimation of d_S by NG is more serious for more extreme codon bias (Figure 3). In Figure 4a, the NG estimates of d_S corresponding to different sequence divergences (t) from the simulation (Figure 3) were plotted against ENC. Although d_S and ENC are independent, NG incorrectly indicated a positive correlation, due to the underestimation of d_S at high divergences and extreme codon bias. The ML method (not plotted) indicated no correlation between d_S and ENC (r² = 0.00006, P = 0.95). The pattern seen in the simulated data is consistent with that in the real data. This is most clear when NG estimates of d_S from the real data (DmDp, Figure 2d) were superimposed onto the approximate estimates from the simulations (Figure 4b). These data suggest that the significantly positive correlation between d_S and ENC indicated by the approximate methods is an artifact of these methods’ failure to properly account for the combined effects of codon bias and sequence divergence.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

—(a) Approximate estimates of d_S by NG obtained in the simulation (Figure 3) plotted against ENC. Multiple estimates of d_S for each ENC correspond to the seven values of t. Note the true d_S is independent of ENC, but a positive correlation is suggested by NG. (b) The same data as in a but with NG estimates from the pairwise comparisons between Drosophila melanogaster and D. pseudoobscura superimposed (○).