CODON adaptation has been well documented in bacterial and fungal genomes ([@bib22]; [@bib28], [@bib29]; [@bib64]) as well as in mitochondrial genomes in vertebrates ([@bib66]; [@bib76]) and fungi ([@bib9]; [@bib68]). Optimizing codon usage according to the codon usage of highly expressed host genes has been shown to increase the production of viral proteins ([@bib25]; [@bib41]) or transgenic genes ([@bib26]; [@bib31]; [@bib32]). Studies on codon--anticodon adaptation have progressed in theoretical elaboration ([@bib6], [@bib7]; [@bib64], [@bib68]; [@bib27]; [@bib30]; [@bib42]) as well as in critical tests of alternative theoretical predictions ([@bib63], [@bib66]; [@bib9]; [@bib61]).

Codon--anticodon adaptation has been documented in bacteriophage (referred to as phage hereafter), partly because several phage species have been used to treat human infections ([@bib51]; [@bib45]; [@bib50]; [@bib1]) or remove infectious biofilms ([@bib2]) and need to be reengineered to improve translation efficiency. While phage codon adaptation is shaped mainly by mutation and transfer (t)RNA-mediated selection, previous studies ([@bib24]; [@bib21]; [@bib34]; [@bib48]; [@bib8]; [@bib37]) have focused mainly on the tRNA-mediated selection on codon usage and have not assessed quantitatively the joint effect of mutation and selection on codon usage bias of phages.

Here we aim to elucidate how biased mutation and selection mediated by host translation machinery will alter the trajectory of codon adaptation in phage protein-coding genes. Many DNA phages, especially single-stranded DNA (ssDNA) phages, experience strong C→T mutation bias mediated by spontaneous or enzymatic deamination ([@bib14]; [@bib35]; [@bib75]). In particular, the spontaneous deamination rate is ∼100 times higher in ssDNA than in double-stranded DNA (dsDNA) based on experimental evidence ([@bib20]), which may explain why some ssDNA viruses, including ssDNA phages, evolve much faster than dsDNA viruses, with their evolutionary rate comparable to that of RNA viruses ([@bib59]; [@bib55]; [@bib75]; [@bib54]; [@bib12], [@bib13]). Oxidative deamination leading to high C→U/T transitional mutation rates has been reported in ssDNA phage M13 ([@bib33]).

The Effect of C→T Mutation Bias {#s1}
===============================

When selection is absent, if the C→T mutation bias experienced by a phage genome is strong enough to overcome stochastic fluctuation of codon frequencies of viral protein-coding genes, then all Y-ending codon families and subfamilies (where Y stands for pyrimidine) in viral protein-coding genes will tend to have the same proportion of U-ending codons; *i.e.*,$$P_{\text{U}.i} = \frac{N_{\text{U}.i}}{N_{\text{U}.i} + N_{\text{C}.i}} = B_{\text{C}\rightarrow\text{T}},$$where *N*~U.~*~i~* and *N*~C.~*~i~* are the numbers of codons ending with U or C, respectively, in codon family *i*, and *B*~C→T~ is a constant representing C→T mutation bias (being 0.5 when there is no C→T mutation bias). When *B*~C→T~ increases, *P*~U~ for all codon families will tend to increase synchronously if not checked by selection. Equation 1 represents a purely mutation-only model of codon usage bias in the viral Y-ending codon families.

When the effect of selection on viral codon usage is negligible, *B*~C→T~ can be approximated simply as the average of *P*~U.~*~i~* values for all Y-ending codon families in the viral protein-coding genes; *i.e.*,$$B_{\text{C}\rightarrow\text{T}} = {\overline{P}}_{\text{U}} = \frac{\sum_{i = 1}^{N_{\text{Y}}}P_{\text{U}.i}}{N_{\text{Y}}},$$where *N*~Y~ is the number of codon families with Y-ending codons. For simplicity, we refer to both Y-ending codon families (*e.g.*, Asn codon family AAY) and Y-ending codon subfamilies (*e.g.*, GGY codons in the Gly codon family) as Y-ending codon families.

Some ssDNA phages have high ${\overline{P}}_{\text{U}}$ values; *e.g.*, ${\overline{P}}_{\text{U}}$ at the third codon position of *Chlamydia* phage Chp1 (Microviridae, NC_001741) is 0.9518, with U-ending codons being invariably the most frequent in all Y-ending or N-ending codon families. Some dsDNA phages can also have high ${\overline{P}}_{\text{U}}$, *e.g.*, being 0.9014 at the third codon position of *Clostridium* phage phi3626 (Siphoviridae, NC_003524). However, codon usage bias is almost always the result of both mutation bias and selection.

The Effect of tRNA-Mediated Selection and Its Characterization {#s2}
==============================================================

A bacterial host may have many tRNAs to read the U-ending codons and few to read the C-ending codons in certain codon families. In such a codon family, a U-ending codon is expected to be decoded efficiently (U-friendly); *i.e.*, tRNA-mediated selection will favor U-ending codons. Similarly, we refer to a codon in which C-ending codons can be decoded more efficiently than U-ending codons as U-hostile. A strong C→T mutation bias would accelerate/enhance codon adaptation in a U-friendly codon family, but would go against codon adaptation in U-hostile codon families. Thus, the degree of U friendliness in the host is expected to be a major determinant of phage codon evolution.

How do we measure U friendliness (*i.e.*, selection in favor of U-ending codons)? We develop a simple index, numerically illustrated in [Figure 1](#fig1){ref-type="fig"}, based on the comparison of codon frequency (CF) between highly expressed host genes (HEGs) and all other genes (non-HEGs), designated by CF~HEG~ and CF~non-HEG~, respectively. Take, for example, the Ala (A) and Phe (F) codon families, where the Y-ending codons are translated by tRNA with a wobble G. In the Ala codon family, GCC is more frequent than GCU when all coding sequences (CDSs) are included. This alone may suggest that the host translation machinery favors C-ending codons. However, with only HEGs from *Escherichia coli*, GCC is much less frequent than GCU, suggesting that U-ending codons are more efficiently translated than C-ending codons in the Ala codon family. The mechanistic explanation for this is that GCC can be decoded only by tRNA^Ala/GGC^, whereas GCU can be decoded by both tRNA^Ala/GGC^ and tRNA^Ala/UGC^, where the wobble U in the anticodon is modified to cmo^5^U, which pairs efficiently with U at the third codon position ([@bib39], [@bib40]; [@bib79]). Similarly, the Phe codon family has more UUU codons than UUC codons when all CDSs are included, but fewer UUU than UUC codons when only *E. coli* HEGs are included. This is also expected because these codons are decoded by tRNA^Phe/GAA^, which prefers UUC over UUU. The observation that UUC is preferred by HEGs suggests that the Phe codon family is not U-friendly (it is C friendly). These illustrations lead us to adopt the association coefficient *φ* as a proxy for U friendliness. The Ala codon family is U-friendly and has a positive *φ*-value, whereas the Phe codon family is U-hostile and has a negative *φ*-value ([Figure 1](#fig1){ref-type="fig"}). *φ* takes values between −1 and 1 and is equivalent to the Pearson correlation coefficient for continuous variables. Because *φ~i~* measures the selection (preference of the host machinery) in favor of the U-ending codons, it is expected to be positively correlated with *P*~U.~*~i~*.

![Rationale of using the *φ*-coefficient as a proxy for U friendliness, based on the codon frequencies (CF) between highly expressed genes (HEGs) and non-HEGs from *E. coli*. *φ* can take values within the range between −1 and 1.](301fig1){#fig1}

Should we develop an index of selection based only on the highly expressed genes? The following scenario suggests that we should not. Suppose the codon frequencies of NNC and NNU from HEGs are 80 and 90, respectively, but those for all CDSs are 200 and 600, respectively. A proper interpretation of this scenario is that extremely high T-biased mutation leads to the dominance of NNU codons. However, the host translation machinery prefers C-ending codons and this selection acts against the T-biased mutation so that codon usage in HEGs is not as U-biased as that in all CDSs. If we have codon usage of only HEGs, we may conclude that the host translation machinery prefers U-ending codons.

A Simple Model of the Joint Effect of Mutation and Selection {#s3}
============================================================

The development of the *φ*-coefficient as a measure of selection for each Y-ending codon family allows us to extend the mutation-only model in Equation 1 to include the effect of selection on *P*~U.~*~i~*; *i.e.*,$$P_{\text{U}.i} = B_{\text{C}\rightarrow\text{T}} + b\phi_{i}.$$Because *P*~U.~*~i~* can be readily computed from viral protein-coding genes, and *φ~i~* can be derived from codon frequencies of host HEGs and non-HEGs ([Figure 1](#fig1){ref-type="fig"}), we can use Equation 3 to quantify the relative importance of mutation (*B*~C→T~) and selection (*φ*) on the codon usage bias of Y-ending codon families. If *φ~i~*-values differ little among Y-ending codon families in a host, then *P*~U.~*~i~* will largely depend on *B*~C→T~, and we will observe little variation in *P*~U.~*~i~* values among different codon families. In contrast, with increasing intensity of selection (a large *b*) or increasing variation in preference of U-ending codons by *E. coli* (*i.e.*, large variation in *φ*), *P*~U.~*~i~* will become more dependent on *bφ~i~*. Similarly, if *B*~C→T~ becomes very large (*i.e.*, very strong mutation bias), then *bφ~i~* naturally would become relatively small and we would conclude that the mutation bias is the dominant factor in shaping codon usage in Y-ending codon families.

One may also argue that *P*~U.~*~i~* cannot be \>1 or \<0, so it will asymptotically approach 1 with increasing *φ~i~* and approach 0 with decreasing *φ~i~*. This implies a sigmoidal relationship between *P*~U~ and *φ*. For this reason, we have also fitted the following sigmoid function,$$P_{\text{U}.i} = \frac{1}{1 + Ce^{- D\phi_{i}}},$$where parameters *C* and *D* are constants. The maximum and minimum values for *P*~U~, according to Equation 4, are 1 and 0, respectively. When *φ* = 0 or *D* = 0, the expected *P*~U~ is 1/(1 + *C*), which is equivalent to *B*~C→T~ in Equation 3. In most cases, *B*~C→T~ and 1/(1 + *C*) are nearly identical and we will use *B*~C→T~ to refer to both as an index of C→T mutation bias. *D* measures the benefit of codon adaptation for the phage. If *D* = 0 (*i.e.*, codon adaptation is not important for the phage), then which codon is favored by the host machinery (measured by *φ*) is irrelevant to phage codon usage. If *D* is very large, then even a codon that is weakly favored by the host will be strongly favored by the phage.

Note that for a given viral species, *B*~C→T~ is constant and affects uniformly the codon usage bias of all Y-ending codon families. In contrast, *φ~i~* is specific to individual codon families. *B*~C→T~ is estimated by the intercept of the linear regression model and selection intensity *b* is the slope. Also note that the correlation coefficient between *P*~U~ and *φ* is also a measure of the effect of selection on codon usage bias (a measure of adaptation) in the Y-ending codon families in phages. We interpret adaptation broadly. For example, suppose that a phage species has evolved good codon adaptation to host species A. If the phage subsequently invaded host species B, and if the codon preference in host species B is exactly the same as that in host species A, then we will state that the phage exhibits good codon adaptation to host species B, although it is preadaptation that is applicable here.

We use Equation 3 to characterize mutation bias and selection intensity based on existing genomic data from dsDNA and ssDNA phages and their hosts. We detected the effect of *φ* in most dsDNA and ssDNA phage species. However, increasing C→T mutation bias significantly reduced the effect of selection in ssDNA phages and shifted the phage codon usage away from the optimum. Some *E. coli* phages such as phage PRD1 whose close relatives are all parasitizing gram-positive bacteria may have recently invaded *E. coli* and have codon usage highly different from that of *E. coli* HEGs. Strand asymmetry with the associated local variation in mutation bias (U-biased in one-half of the genome and C-biased in the other half) can significantly interfere with codon adaptation in both dsDNA and ssDNA phages. Much of the variation in codon adaptation among dsDNA phages can be attributed to lineage effects, with some phage lineages having uniformly strong codon adaptation and some other lineages having uniformly weak codon adaptation.

Materials and Methods {#s4}
=====================

Genomic data and processing {#s5}
---------------------------

The genome sequences of 469 dsDNA phages, 41 ssDNA phages, and their corresponding bacterial hosts were downloaded from GenBank, of which 71 have *E. coli* specified as their host in the "/HOST" tag in the "FEATURES" table, including 60 dsDNA phages and 11 ssDNA phages. The CDSs and codon usage data, as well as three codon positions, were extracted using DAMBE ([@bib73]). All phage genomes were searched for encoded tRNAs by using the tRNAscan-SE Search Server ([@bib53]). The local TC skew plot, with the TC skew computed as (*N*~T~ − *N*~C~)/(*N*~T~ + *N*~C~), where *N~i~* is the number of nucleotides *i* along a moving window, was generated from DAMBE ([@bib73]). All statistical analyses were done with SAS ([@bib49]), with the linear regression fitted by the GLM procedure and the sigmoid function by the NLIN procedure.

*E. coli* has 29 strains with RefSeq genomic sequences, but the /HOST tag in a viral genome gives only species name (*i.e.*, *E. coli*), with no strain-specific information. For this reason, all 29 RefSeq genomic sequences were downloaded and *E. coli* codon usage is computed as the average of all CDSs from these 29 genomes. The codon usage of highly expressed *E. coli* genes was compiled in the Eeco_h.cut file distributed with EMBOSS ([@bib46]). It is almost perfectly correlated with our own compilation of codon usage from all *E. coli* ribosomal proteins (which are necessarily highly expressed because of the high density of ribosomes in the cell). There is little variation in codon usage in highly expressed genes among different *E. coli* genomes.

Indexes of codon usage bias {#s6}
---------------------------

While we mainly focus on modeling mutation and selection on *P*~U~ in Equations 3 and 4, two indexes of codon usage bias were used to aid in the interpretation of the results: the codon adaptation index (CAI) ([@bib56]) with the improved implementation ([@bib67]) and the effective number of codons (*N*~c~) ([@bib62]) with the improved implementation ([@bib58]). All these indexes were computed using DAMBE ([@bib73]). For computing the phage CAI, the host highly expressed genes are used as the reference set of genes. Only CDSs with at least 33 codons (99 nt) are included in computing the indexes of codon usage bias to alleviate stochastic noise in computing these indexes with few codons.

Phylogenetic analysis {#s7}
---------------------

Coancestry of phage species is difficult to establish. Although some dsDNA phage genomes are annotated to contain a DNA polymerase gene, the gene sequences from different phage lineages are often not homologous and cannot be aligned. We build phage "phylogenetic" trees by using a composition vector approach called CVTree ([@bib78]) that does not require aligned sequences but implicitly assumes the sharing of ancestral peptides as phylogenetic signals. The method uses amino acid sequences and is conceptually based on the sharing of ancient peptides that give individual evolutionary lineages their uniqueness. Computationally, the method is built upon the similarities in the sharing of words of length *k* (*a*~1~*a*~2~ ... *a~k~*) after subtracting its random expectation based on the frequencies of *a*~1~*a*~2~ ... *a~k~*~−1~, *a*~2~*a*~3~ ... *a~k~*, and *a*~2~*a*~3~ ... *a~k~*~−1~. The CVTree method has been implemented in the most recent version of DAMBE ([@bib73]). We used a *k* value of 5, which has been recommended for viral genomes ([@bib78]). The data for reconstructing phylogenetic trees with the CVTree method are .faa files downloaded from GenBank, with each .faa file containing all annotated amino acid sequences for each phage species.

Results and Discussion {#s8}
======================

Codon preference by the *E. coli* translation machinery: *φ* {#s9}
------------------------------------------------------------

The *φ*-values for *E. coli* Y-ending codons are generally small, ranging from −0.1512 to 0.1424 ([Table 1](#t1){ref-type="table"}). Among the 16 Y-ending codon families and subfamilies, 7 are U-friendly (with *φ* \> 0, with the mean = 0.0861) and 9 are U-hostile (with *φ* \< 0, with the mean = −0.1025). Thus, C-ending codons overall should be slightly favored over U-ending codons to achieve the codon usage pattern of highly expressed host genes, which is consistent with the proportion of U-ending codons in Y-ending codon families in *E. coli* highly expressed genes (*P*~U.~*~E.coli~* = 0.4421). Increased C→T mutation bias will improve codon adaptation for the 7 U-friendly codon families, but will lead to deterioration in the 9 U-hostile codon families ([Table 1](#t1){ref-type="table"}).

###### Codon frequencies (CF) for Y-ending codons in *E. coli*, compiled for highly expressed genes (HEG) and all other genes (non-HEG), together with the gene copy number of tRNA in the genome (strain K12) whose anticodon matches the codon, and *φ* as a measure of codon preference of the host translation machinery (a large *φ* corresponds to greater preference of U-ending codons over C-ending codons)

  AA   Codon   CF~non-HEG~   CF~HEG~   tRNA                                   *φ*
  ---- ------- ------------- --------- -------------------------------------- ---------
  A    GCC     33,463        1,306     2                                      0.1424
  A    GCT     18,526        2,288                                            
  C    TGC     8,397         475       1                                      −0.0362
  C    TGT     6,802         270                                              
  D    GAC     23,226        2,786     3                                      −0.0993
  D    GAT     41,472        2,345                                            
  F    TTC     20,332        2,229     2                                      −0.1478
  F    TTT     29,556        872                                              
  G    GGC     37,418        2,987     4                                      0.0566
  G    GGT     30,154        3,583                                            
  H    CAC     12,144        1,160     1                                      −0.1331
  H    CAT     17,170        477                                              
  I    ATC     30,787        3,488     3                                      −0.1232
  I    ATT     39,788        1,640                                            
  L    CTC     14,591        541       1                                      −0.0353
  L    CTT     14,679        357                                              
  N    AAC     26,674        2,832     4                                      −0.1512
  N    AAT     23,652        539                                              
  P    CCC     7,443         38        1                                      0.1032
  P    CCT     9,235         343                                              
  R    CGC     28,473        1,530                                            0.1011
  R    CGT     25,528        2,995     3[*^a^*](#t1n1){ref-type="table-fn"}   
  S    AGC     20,868        1,015     1                                      −0.0842
  S    AGT     11,802        168                                              
  S    TCC     10,649        1,110     2                                      0.0327
  S    TCT     10,217        1,320                                            
  T    ACC     29,335        2,533     2                                      0.0408
  T    ACT     10,950        1,286                                            
  V    GTC     19,972        824       2                                      0.1262
  V    GTT     22,297        2,669                                            
  Y    TAC     15,094        1,569     3                                      −0.1122
  Y    TAT     21,207        865                                              

*P*-values from a chi-square test of 2 × 2 contingency tables, with the null hypothesis that *φ* = 0, are all \<0.0001.

The anticodon has a wobble A modified to inosine.

Effect of mutation and selection on codon usage of *E. coli* ssDNA phages {#s10}
-------------------------------------------------------------------------

If selection in favor of U-ending codons by the host translation machinery (*φ*) is efficient, then we expect *P*~U~ to increase with *φ*. This expectation is consistent with data from *E. coli* Enterobacteria phage G4 (NC_001420) showing *P*~U~ increasing roughly linearly with *φ* ([Figure 2](#fig2){ref-type="fig"}). Fitting the linear model in Equation 3 results in *B*~C→T~ = 0.5707, and *b* = 0.6609, with the relationship being statistically significant (*P* = 0.0159). Fitting the sigmoid function in Equation 4 yields *C* = 0.7475, *D* = 2.7241, and 1/(1 + *C*) = 0.5722, which is equivalent to *B*~C→T~ in the linear model, *i.e.*, both being the expected *P*~U~ value when *Dφ* in Equation 4 is zero (*i.e.*, no selection on phage codon usage in Y-ending codon families). The predicted values from the linear model in Equation 3 and the nonlinear model in Equation 4 are identical to the first two digits after the decimal point, indicating sufficiency of the linear model.

![Relationship between *P*~U~ (the proportion of U-ending codons in Y-ending codon families) and *φ* (selection in favor of U-ending codons), based on codon usage data from *E. coli* Enterobacteria phage G4 (NC_001420). Also shown is the linear fit to the data. Applying the sigmoid function in Equation 4 generated effectively the same predicted values.](301fig2){#fig2}

The estimated *B*~C→T~ from applying the regression model in Equation 3 to all 11 ssDNA Enterobacteria phages parasitizing *E. coli* varies from 0.5443 to 0.7419 ([Table 2](#t2){ref-type="table"}). These would be the *P*~U~ values when selection mediated by the host translation machinery is absent. With the slopes in [Table 2](#t2){ref-type="table"}, the effect of selection on viral *P*~U~ values is small relative to *B*~C→T~. We thus expect the estimated *B*~C→T~ values to be close to the empirical ${\overline{P}}_{U}$ values defined in Equation 2, which is true ([Figure 3](#fig3){ref-type="fig"}).

###### Results of fitting the linear regression model in Equation 3 to codon usage in ssDNA Enterobacteria phages parasitizing *E. coli*, with viral genome accession number (ACCN), viral genome length (*L*), number of viral genes (*N*~g~), the estimated intercept (*B*~C→T~) and slope (*b*), the Pearson correlation between *P*~U~ and *φ* for each phage species, and the statistical significance (two-tailed *P*) of the relationship

  Phage           ACCN        *L*    *N*~g~   *B*~C→T~   *b*       *R*       *P*
  --------------- ----------- ------ -------- ---------- --------- --------- --------
  Phage alpha3    NC_001330   6087   10       0.7022     0.3790    0.3936    0.1315
  Phage G4        NC_001420   5577   11       0.5707     0.6609    0.5912    0.0159
  Phage ID18      NC_007856   5486   11       0.5876     0.6551    0.5105    0.0433
  Phage ID2       NC_007817   5486   11       0.5443     0.5820    0.5055    0.0458
  Phage phiX174   NC_001422   5386   11       0.6965     0.2840    0.2969    0.2641
  Phage St-1      NC_012868   6094   11       0.6881     0.2895    0.3048    0.2511
  Phage WA13      NC_007821   6068   10       0.7122     0.2190    0.2354    0.3801
  Phage I2-2      NC_001332   6744   9        0.6987     0.3301    0.2649    0.3214
  Phage If1       NC_001954   8454   10       0.6551     −0.0763   −0.0819   0.7629
  Phage Ike       NC_002014   6883   10       0.7419     0.1614    0.1097    0.6858
  Phage M13       NC_003287   6407   10       0.7390     0.1717    0.1899    0.4812

![The average ${\overline{P}}_{\text{U}}$ defined in Equation 2 is similar to *B*~C→T~ estimated from fitting the linear model in Equation 3, based on 11 ssDNA Enterobacteria phages parasitizing *E. coli* ([Table 2](#t2){ref-type="table"}).](301fig3){#fig3}

The standard error associated with the *B*~C→T~ values is on the order of 0.02 (not shown), so that *B*~C→T~ values in [Table 2](#t2){ref-type="table"} are all significantly greater than the observed ${\overline{P}}_{U}$ (= 0.4421) in *E. coli* HEGs. If we assume that the codon usage of *E. coli* HEGs represents the optimum achievable given the counterbalance between mutation and selection, then the large *B*~C→T~ values in [Table 2](#t2){ref-type="table"} suggest that C→T-biased mutation in ssDNA has shifted the codon usage of ssDNA phages away from the optimum.

*B*~C→T~ has a strong effect on the effective number of codons (*N*~c~), as expected. *N*~c~ is at its maximum when *B*~C→T~ ∼ 0.5, but decreases sharply as *B*~C→T~ increases, leading to U-ending codons dominating over C-ending codons. However, *B*~C→T~ has little effect on CAI, partly because *E. coli* translation machinery favors U-ending codons in about half of the Y-ending codon families and C-ending codons in the other half ([Table 1](#t1){ref-type="table"}). A large *B*~C→T~ will increase the frequency of U-ending codons in both the U-friendly and U-hostile codon families. The positive effect in the U-friendly codon families is offset by the negative effect on U-hostile codon families.

A well-adapted codon usage in a phage species in *E. coli* should have *P*~U~ positively and highly correlated with *φ*, *i.e.*, large *P*~U~ in U-friendly codon families (with large *φ*-values) and small *P*~U~ in U-hostile codon families (with small *φ*-values). However, a strong C→T mutation bias (a large *B*~C→T~) will lead to high *P*~U~ in all Y-ending codon families, resulting in reduced correlation between *P*~U~ and *φ*. We therefore expect the correlation between *P*~U~ and *φ* (*R*) to decrease with increasing *B*~C→T~. This expectation is consistent with the empirical data ([Figure 4](#fig4){ref-type="fig"}), although there is one outlying point (Enterobacteria phage If1; NC_001954) for which we offer an explanation later.

![The correlation (*R*) between *P*~U~ and *φ* decreases with increasing *B*~C→T~. The outline point is *E. coli* Enterobacteria phage If1 (NC_001954). The negative association is statistically significant (*P* = 0.0292 with the outlying point included). The four red circles form a monophyletic taxon and the rest form another monophyletic taxon ([Figure 5](#fig5){ref-type="fig"}).](301fig4){#fig4}

The four ssDNA phage species (phage I2-2, IF1, Ike, and M13) with low correlation between *P*~U~ and *φ* (red diamonds in [Figure 4](#fig4){ref-type="fig"}) have codon usages significantly correlated with each other, which suggests that they might be phylogenetically related. The tree built with the CVTree algorithm ([@bib78]) implemented in DAMBE ([@bib73]) does cluster these four species, all belonging to Inoviridae, into a monophyletic taxon ([Figure 5](#fig5){ref-type="fig"}). Other viral proteomic trees ([@bib47]; [@bib16]) also group these four *E. coli* ssDNA phages into the same clade. The other seven ssDNA phages, all belonging to Microviridae, are also clustered into a monophyletic taxon ([Figure 5](#fig5){ref-type="fig"}).

![Phylogenetic tree of ssDNA phages reconstructed by using the CVTree method ([@bib78]) implemented in [@bib73]. The four phage species colored in blue belong to Inoviridae whereas the other seven phages colored in red belong to Microviridae. The Operational Taxonomic Units (OTUs) are formed by a combination of host (the first letter of the genus name and the first four letters of the host species name), phage species name, GenBank accession number, *R* (correlation between *P*~U~ and *φ*), estimated *B*~C→T~, and number of tRNA genes (*N*\_tRNA) in the phage genome.](301fig5){#fig5}

Because of the phylogenetic structure ([Figure 5](#fig5){ref-type="fig"}), one may argue that the 11 points are not statistically independent and question the validity of using conventional regression to test the significance of the negative association between *R* and *B*~C→T~. For example, the ancestor of the seven phages in Microviridae (colored red in [Figure 5](#fig5){ref-type="fig"}) may have codon usage similar to that in *E. coli*, which was then inherited by all seven descendant phage lineages. Similarly, the ancestor of the four phage species in Inoviridae (colored blue in [Figure 5](#fig5){ref-type="fig"}) may have codon usage different from that of *E. coli*, which was then inherited by its four descendant phage lineages. Thus, we would, in an extreme case, have only two data points. To overcome this problem, we used the subtree for the 11 species in [Figure 5](#fig5){ref-type="fig"} and performed independent-contrasts analysis ([@bib18], [@bib19], pp. 435--443) implemented in DAMBE ([@bib72], pp. 24--29; [@bib73]) and found the negative association still significant (*P* \< 0.05).

We may conclude from the results above that the overall effect of selection in ssDNA phage is statistically significant as the mean of the 11 *b* values is significantly \>0 (mean *b* = 0.3324, SE = 0.0685, *t* = 4.8546, degrees of freedom (d.f.) = 10, *P* = 0.0007). However, when the false discovery rate method ([@bib4]; [@bib5]) is used to control for type I error rate involving multiple comparisons, none of the 11 individual *P*-values is statistically significant at the 0.05 significance level. The false discovery rate method has been numerically illustrated ([@bib69]) and implemented in DAMBE ([@bib73]).

Effect of mutation and selection, as well as evolutionary history, on codon usage of *E. coli* dsDNA phages {#s11}
-----------------------------------------------------------------------------------------------------------

Some dsDNA phages show strong response to selection by the host translation machinery (*φ*), *e.g.*, NC_010324 phage Phieco32, with *P*~U~ strongly dependent on *φ* ([Figure 6](#fig6){ref-type="fig"}). The estimated *B*~C→T~ spans a wide range ([Table 3](#t3){ref-type="table"}), but on average is significantly smaller than that of ssDNA phages (*t*-test, *t* = 2.1379, d.f. = 69, two-tailed *P* = 0.0361), suggesting a weaker effect of C→T mutation on codon usage in dsDNA phages than in ssDNA phages. The *b* values also vary substantially ([Table 3](#t3){ref-type="table"}). For example, phage BP-4795 (NC_004813), phage cdtI (NC_009514), and 10 other dsDNA phages have negative slopes between *P*~U~ and *φ* ([Table 3](#t3){ref-type="table"}), contrary to what we would have predicted based on codon adaptation. While it is easy to understand why phages should exhibit codon adaptation (with a positive slope *b*) to the host, it is puzzling why some dsDNA phages do not evolve codon usage similar to that of the host.

![Relationship between *P*~U~ (the proportion of U-ending codons in Y-ending codon families) and *φ* (selection in favor of U-ending codons), based on codon usage data from *E. coli* Enterobacteria phage Phieco32 (NC_010324).](301fig6){#fig6}

###### Results of fitting the linear regression model in Equation 3 to codon usage in dsDNA *E. coli* phages, with viral genome accession number (ACCN), the estimated intercept (*B*~C→T~) and slope (*b*), the Pearson correlation between *P*~U~ and *φ* for each phage species, and the statistical significance (two-tailed *P*) of the relationship

  Phage                        ACCN        *B*~C→T~   *b*       *R*       *P*[*a*](#t3n1){ref-type="table-fn"}
  ---------------------------- ----------- ---------- --------- --------- --------------------------------------------
  Phage 13a                    NC_011045   0.5690     1.3782    0.8381    0.00005[*^b^*](#t3n2){ref-type="table-fn"}
  Phage 285P                   NC_015249   0.5769     1.4737    0.8751    0.00001[*^b^*](#t3n2){ref-type="table-fn"}
  Phage 933W                   NC_000924   0.5543     0.0089    0.0129    0.96205
  Phage BP-4795                NC_004813   0.5210     −0.2133   −0.2768   0.29931
  Phage CC31                   NC_014662   0.7277     1.0780    0.8437    0.00004[*^b^*](#t3n2){ref-type="table-fn"}
  Phage cdtI                   NC_009514   0.5632     −0.2310   −0.2948   0.26777
  Phage EcoDS1                 NC_011042   0.5477     1.4925    0.8077    0.00015[*^b^*](#t3n2){ref-type="table-fn"}
  Phage EPS7                   NC_010583   0.7368     0.8949    0.9011    0.00000[*^b^*](#t3n2){ref-type="table-fn"}
  Phage HK022                  NC_002166   0.4858     0.1603    0.1659    0.53920
  Phage HK97                   NC_002167   0.4951     0.2864    0.2699    0.31203
  Phage IME08                  NC_014260   0.7156     0.6011    0.6583    0.00557[*c*](#t3n3){ref-type="table-fn"}
  Phage JK06                   NC_007291   0.6157     0.7818    0.5302    0.03462
  Phage JS10                   NC_012741   0.7206     0.6104    0.6538    0.00601[*^c^*](#t3n3){ref-type="table-fn"}
  Phage JS98                   NC_010105   0.7200     0.6300    0.6684    0.00464[*^c^*](#t3n3){ref-type="table-fn"}
  Phage JSE                    NC_012740   0.6931     0.6136    0.5988    0.01425[*^c^*](#t3n3){ref-type="table-fn"}
  Phage K1-5                   NC_008152   0.6612     0.9844    0.8118    0.00013[*^b^*](#t3n2){ref-type="table-fn"}
  Phage K1E                    NC_007637   0.6605     0.9377    0.8209    0.00010[*^b^*](#t3n2){ref-type="table-fn"}
  Phage K1F                    NC_007456   0.6605     0.9377    0.8209    0.00010[*^b^*](#t3n2){ref-type="table-fn"}
  Phage lambda                 NC_001416   0.4971     −0.1962   −0.1908   0.47894
  Phage Min27                  NC_010237   0.5504     0.0858    0.1302    0.63088
  Phage Mu                     NC_000929   0.4980     −0.6273   −0.5245   0.03700
  Phage N15                    NC_001901   0.4834     0.0208    0.0231    0.93224
  Phage N4                     NC_008720   0.7556     0.9110    0.8554    0.00002[*^b^*](#t3n2){ref-type="table-fn"}
  Phage P1                     NC_005856   0.5475     0.0698    0.0885    0.74454
  Phage P2                     NC_001895   0.4768     −0.4703   −0.4765   0.06202
  Phage P4                     NC_001609   0.5065     −0.6610   −0.5315   0.03412
  Phage Phi1                   NC_009821   0.6905     0.5815    0.5728    0.02039[*^c^*](#t3n3){ref-type="table-fn"}
  Phage Phieco32               NC_010324   0.7277     1.1889    0.9436    0.00000[*^b^*](#t3n2){ref-type="table-fn"}
  Phage phiEcoM-GJ1            NC_010106   0.6421     1.0372    0.8690    0.00001[*^b^*](#t3n2){ref-type="table-fn"}
  Phage phiP27                 NC_003356   0.5451     −0.3341   −0.3569   0.17473
  Phage PRD1                   NC_001421   0.5324     −0.5638   −0.2643   0.32262
  Phage RB16                   NC_014467   0.6257     1.0301    0.7718    0.00046[*^b^*](#t3n2){ref-type="table-fn"}
  Phage RB49                   NC_005066   0.6919     0.5836    0.5790    0.01878[*^c^*](#t3n3){ref-type="table-fn"}
  Phage RB69                   NC_004928   0.7549     0.6396    0.7414    0.00101[*^b^*](#t3n2){ref-type="table-fn"}
  Phage RTP                    NC_007603   0.6092     0.8781    0.5433    0.02964
  Phage SfV                    NC_003444   0.5316     −0.1209   −0.1431   0.59705
  Phage SPC35                  NC_015269   0.7544     0.6200    0.8677    0.00001[*^b^*](#t3n2){ref-type="table-fn"}
  Phage SSL-2009a              NC_012223   0.3735     0.2555    0.3003    0.25843
  Phage T1                     NC_005833   0.5786     0.9154    0.5723    0.02053[*^c^*](#t3n3){ref-type="table-fn"}
  Phage T3                     NC_003298   0.5336     1.3661    0.8663    0.00001[*^b^*](#t3n2){ref-type="table-fn"}
  Phage T4                     NC_000866   0.7967     0.3770    0.5538    0.02603
  Phage T5                     NC_005859   0.7386     0.5950    0.8399    0.00005[*^b^*](#t3n2){ref-type="table-fn"}
  Phage T7                     NC_001604   0.5641     1.3971    0.8461    0.00004[*^b^*](#t3n2){ref-type="table-fn"}
  Phage TLS                    NC_009540   0.6377     0.5840    0.3261    0.21766
  Phage vB_EcoM-VR7            NC_014792   0.6945     0.4866    0.6418    0.00736[*^c^*](#t3n3){ref-type="table-fn"}
  Phage VT2-Sakai              NC_000902   0.5434     0.0377    0.0517    0.84906
  Phage WV8                    NC_012749   0.7372     1.0841    0.9077    0.00000[*^b^*](#t3n2){ref-type="table-fn"}
  Phage bV_EcoS_AKFV3          NC_017969   0.7444     0.5894    0.8094    0.00015[*^b^*](#t3n2){ref-type="table-fn"}
  Phage D108                   NC_013594   0.4994     −0.5957   −0.4791   0.06045
  Phage HK639                  NC_016158   0.4491     0.2538    0.2437    0.36316
  Phage HK75                   NC_016160   0.4839     0.3029    0.3018    0.25602
  Phage phiV10                 NC_007804   0.5705     0.4618    0.4478    0.08199
  Phage rv5                    NC_011041   0.6296     0.8548    0.6931    0.00291[*^b^*](#t3n2){ref-type="table-fn"}
  Phage vB_EcoM_CBA12          NC_016570   0.5988     0.3685    0.4449    0.08425
  Stx1-converting bac          NC_004913   0.5512     0.0121    0.0165    0.95162
  Phage BA14                   NC_011040   0.5691     1.4057    0.8765    0.00001[*^b^*](#t3n2){ref-type="table-fn"}
  Stx2-converting phage II     NC_004914   0.5479     0.0495    0.0706    0.79503
  Stx2-converting phage 1717   NC_011357   0.5011     −0.1374   −0.1724   0.52309
  Stx2-converting phage 86     NC_008464   0.5681     0.0285    0.0417    0.87813
  Stx2-converting phage I      NC_003525   0.5132     −0.0566   −0.0986   0.71646

Significant at the 0.05 level when experimentwise error rate is controlled by the false discovery rate method.

Significant with the more conservative approach of [@bib5].

Significant with the approach of [@bib4].

We noted that all phages that have negative correlation between *P*~U~ and *φ* share similar codon usages. For example, *P*~U~ values from phage BP-4795 (NC_004813) and those from phage cdtI (NC_009514) have a correlation coefficient of 0.9349. The observation that they have similar codon usages that are different from that of their host increases the plausibility that they may have adapted to a common host that has codon usage different from that of *E. coli* and that they have invaded *E. coli* recently and have not yet had enough time to evolve codon adaptation. The shared correlation among the 12 phages with negative slopes ([Table 3](#t3){ref-type="table"}), summarized in the first principal component (PC1), accounts for 81% of the total variation. All 12 phages with a negative slope have *P*~U~ positively correlated with PC1 ([Figure 7](#fig7){ref-type="fig"}). Here we offer two explanations for the lack of codon adaptation in these 12 phages with empirical substantiation.

![Phage species that do not exhibit codon adaptation to *E. coli* (negative correlation between *P*~U~ and *φ*) nevertheless exhibit codon usage similar to each other in their Y-ending codon families, with the shared correlation summarized in the first principal component (PC1) that accounts for 81% of total variation. The phage genomes are identified by their GenBank accession number and species name.](301fig7){#fig7}

### Phylogenetic inertia: {#s12}

If a phage has only recently invaded *E. coli*, it would have little time to evolve codon adaptation to *E. coli* translation machinery. This explanation may be applicable to phage PRD1, which has the fourth most negative slope (*b* = −0.5638, [Table 3](#t3){ref-type="table"}). Phage PRD1 belongs to the peculiar Tectiviridae family with members parasitizing both gram-negative and gram-positive bacteria. Phage PRD1 is the only species in the family known to parasitize a variety of gram-negative bacteria, including *Salmonella*, *Pseudomonas*, *Escherichia*, *Proteus*, *Vibrio*, *Acinetobacter*, and *Serratia* species ([@bib3]; [@bib23]). This wide host range might lead one to think that the poor codon adaptation of phage PRD1 to *E. coli* is because the phage is not *E. coli* specific. However, other lines of evidence suggest that this is not true. First, the gram-negative hosts that phage PRD1 parasitizes have similar codon usage and adaptation to any one of them or to the average of all of them and will not lead to a negative *b* ([Table 3](#t3){ref-type="table"}). Second, other members of the phage family, *i.e.*, phages PR3, PR4, PR5, L17, and PR772, parasitize gram-positive bacteria. Phage PRD1 is extremely similar to its sister lineages parasitizing gram-positive bacteria. For example, there is only one amino acid difference in the coat protein between PRDl and PR4 ([@bib3]). It is thus quite likely that the ancestor of phage PRD1 parasitizes gram-positive bacteria. The lineage leading to phage PRD1 may have switched to gram-negative bacterial hosts only recently and still has its codon usage similar to that of its ancestral gram-positive bacterial host. In support of this, among 87 bacterial genomes covering major groups of bacterial species, the host species with codon usage most similar to that of phage PRD1 are strains in the gram-positive *Geobacillus* (NC_014206, NC_012793, NC_014650, NC_014915, and NC_013411).

Phylogenetic reconstruction with the CVTree method ([@bib78]) suggests that evolutionary history may have contributed to the differences in codon adaptation in *E. coli* dsDNA phages because codon usage captured by the intercept (*B*~C→T~) and the correlation (*R*) is similar among related phage species. Descendants of node A ([Figure 8](#fig8){ref-type="fig"}) belong to phage family Podoviridae, and all have high *R* and a narrow range of *B*~C→T~ values. None of them encode tRNA genes in their genomes, whereas tRNA genes are present in many other *E. coli* dsDNA phage lineages. Another podovirus, phage Phi eco32 (NC_010324), has the highest correlation between *P*~U~ and *φ* (*r* = 0.9436). It is likely that the ancestor of podoviruses evolved good codon adaptation to an *E. coli*-like bacterial species, which was then inherited by its descendants. Because of the good codon adaptation, there is no need for the phages to carry their own tRNA genes, and all *E. coli* podoviruses studied do not encode tRNA genes in their genomes, except for phage Phi eco32 (NC_010324), which carried one putative tRNA^Arg^ of uncertain function. The sequence is incorrectly annotated because the tRNA^Arg^ sequence cannot be folded into a proper 7-nt anticodon loop for Arg. It also has an extraordinarily long branch when aligned and clustered with any of the *E. coli* tRNA genes, suggesting that it is unlikely to be used by *E. coli* translation machinery even if it is transcribed.

![Phylogenetic tree of dsDNA phages reconstructed by using the CVTree method with *k* (peptide length) = 5. The OTUs are formed by a combination of host (the first letter of the genus name and the first four letters of the host species name), phage species name, GenBank accession number, *R* (correlation between *P*~U~ and *φ*), estimated *B*~C→T~, and number of tRNA genes (tRNA) in the phage genome.](301fig8){#fig8}

In contrast to podoviruses, a number of myoviruses have phage-encoded tRNA genes. For example, Enterobacteria phage WV8 (NC_012749) and *Erwinia* phage phiEa21-4 (NC_011811) have 19 and 23 tRNA genes, respectively (excluding one tRNA pseudogene in phage WV8). Enterobacteria phage WV8 has excellent codon adaptation, with *R* between *P*~U~ and *φ* being 0.9077. One may wonder why Enterobacteria phage WV8, with excellent codon adaptation at least for the Y-ending codon families and subfamilies, should still keep its set of tRNA genes. One possibility is that it is a generalist with more than one host. Previous studies have already suggested an association of host diversity and the number of tRNA genes carried on phage genomes ([@bib52]; [@bib17]). Another possibility is that the Enterobacteria phage WV8 is already in the process of losing its tRNA genes because it has fewer tRNA genes than its sister lineage *Erwinia* phage phiEa21-4, which has 23 tRNA genes. Furthermore, similar to the annotated tRNA^Arg^ "gene" in phage Phi eco32, annotated tRNA genes in phage WV8 are also quite different from their *E. coli* counterparts and may not be used by *E. coli* translation machinery. They may be either nonfunctional or functional in non-*E. coli* hosts.

While phage WV8 exhibits high codon adaptation to *E. coli*, other myoviruses such as those under node C have poor codon adaptation with the correlation either close to zero or negative ([Figure 8](#fig8){ref-type="fig"}, [Table 3](#t3){ref-type="table"}). Strong heterogeneity in codon adaptation among myoviruses is also visible among species under cluster F ([Figure 8](#fig8){ref-type="fig"}), with one myovirus having negative *R* and the rest have positive *R*. Another cluster of species with poor codon adaptation to *E. coli* (having *R* close to 0 or negative) are those under node B, made mainly of lambdoid phages. While the term "lambdoid" is never intended as a taxonomic term, the clustering of these species together suggests phylogenetic affiliations.

The lambdoid phages must have evolved in *E. coli* and *E. coli*-like hosts for a long time, and it would be weak to invoke phylogenetic inertia as an explanation for the lack of concordance of their codon usage with that of *E. coli*. This forces us to reexamine the assumptions of our model in Equations 3 and 4 in search of an alternative explanation for phages with poor codon adaptation to *E. coli*. In the two models specified in Equations 3 and 4, we have assumed a uniform mutation bias that will affect all genes and all Y-ending codon families. However, because of the asymmetric replication of the two DNA strands with associated asymmetric mutation bias, local mutation bias is often not accounted for by *P*~U.~*~i~*~.~ Differences in mutation bias between the two DNA strands has been documented in organisms ranging from viruses and organelles to prokaryotic and eukaryotic genomes ([@bib38]; [@bib69],[@bib71]).

### Strand asymmetry, local mutation bias, and phage codon adaptation: {#s13}

We found that *E. coli* dsDNA phages with negative slopes exhibit stronger local strand asymmetry than those phages with positive slopes. For example, among the 60 dsDNA *E. coli* phages, the three phages (P4, NC_001609; Mu, NC_000929; and D108, NC_013594) with the most negative *b* values (−0.6610, −0.6273, and −0.5957, respectively, [Table 3](#t3){ref-type="table"}) exhibit strong local TC skew, defined as (*N*~T~ − *N*~C~)/(*N*~T~ + *N*~C~), relative to the three phage species (phage BA14, NC_011040; phage 285P, NC_015249; and phage EcoDS1, NC_011042) with the largest *b* values ([Figure 9](#fig9){ref-type="fig"}).These results suggest that codon adaptation may be difficult to achieve when one part of the genome experiences strong T-biased mutation and the other part strong C-biased mutation. Such local mutation bias is obscured if we consider only global codon frequencies over all phage genes.

![TC skew, defined as (*N*~T~ − *N*~C~)/(*N*~T~ + *N*~C~), for three phages with the most negative slopes and three phages with the most positive slopes, plotted over sliding windows (window width = 2000 nt and step size = 200). Phage P4 is much shorter than the others and requires coinfection by phage P2 to complete its lytic cycle. The key shows the phage name with the GenBank accession number and the slope from regression of *P*~U~ on *φ*.](301fig9){#fig9}

How strong local strand asymmetry affects codon adaptation is not immediately obvious, so we offer an illustration here. Take, for example, the phage P4 genome (NC_001609) with 14 genes. We first need to recognize that C→T mutations often lead to not only increased U at the third codon position, but also increased U at the first and second codon positions ([Figure 10](#fig10){ref-type="fig"}). A similar response of nonsynonymous mutation rate to directional mutation pressure has also been documented in several other studies ([@bib57]; [@bib36]; [@bib60]). The sites before genomic position 4500 (five genes) are relatively T rich and those after (nine genes) are relatively C rich (which is obvious from [Figure 9](#fig9){ref-type="fig"}). T-biased mutation reduces codons such as CCY (which is U-friendly as highly expressed *E. coli* genes strongly prefer CCU over CCC), so that CCY are found mostly at T-poor (and C-rich) regions and present mainly as CCC, which are not favored by *E. coli* translation machinery. In contrast, codons such as UUY (which is U-hostile because highly expressed *E. coli* genes strongly prefer UUC over UUU) are found mostly in T-rich regions and present mainly in the unfavored UUU form. Thus, the T-rich region features many unfavored UUU codons and C-rich regions feature many unfavored CCC codons, leading to poor codon adaptation.

![TC skew at the first and second codon positions (SKEW~TC12~) in the 14 genes in phage P4 increases with TC skew at the third codon position (SKEW~TC3~).](301fig10){#fig10}

The effect of strand asymmetry on codon adaptation observed in dsDNA phages is also visible in ssDNA phages. To show this quantitatively, we used the same window size and step size and computed the variance of the window-specific skew values as an index of strand asymmetry (*I*~SA~). Take the six dsDNA phage species in [Figure 9](#fig9){ref-type="fig"}, for example. The *I*~SA~ value would be much greater for the three species with negative slopes than for the three species with positive slopes. The *R* value is highly significantly and negatively correlated with *I*~SA~ for ssDNA phages (*P* = 0.0008, [Figure 11](#fig11){ref-type="fig"}). The same negative relationship between *R* and *I*~SA~ holds for dsDNA phages with *P* \< 0.0001. Thus, mutation bias along different parts of the phage genome in opposite directions can significantly reduce the efficiency of selection on codon usage of both ssDNA and dsDNA phages. To properly assess the effect of mutation bias and selection by the host translation machinery, it is important to apply Equations 3 and 4 to phage genomic segments with relatively homogenous mutation bias.

![The effect of selection on Y-ending codons, measured by the correlation (*R*) between *P*~U~ and *φ*, decreases with the degree of strand asymmetry, measured by the index of strand asymmetry (*I*~SA~, which is the variance of the window-specific TC skew values).](301fig11){#fig11}

The relationship between *R* and *I*~SA~ in [Figure 11](#fig11){ref-type="fig"} offers an explanation for the outlying point in [Figure 4](#fig4){ref-type="fig"}, where Enterobacteria phage If1 (NC_001954) has an *R* value much smaller than expected from the general trend. This phage has the strongest strand asymmetry (*i.e.*, the largest *I*~SA~ value) among ssDNA phages and in this new light is expected to be associated with a low *R* value ([Figure 11](#fig11){ref-type="fig"}). In short, the mutation effect (on the *x*-axis) for Enterobacteria phage If1 is underestimated in [Figure 4](#fig4){ref-type="fig"}, which does not take local strand asymmetry into consideration. When local strand asymmetry is accounted for ([Figure 11](#fig11){ref-type="fig"}), the point is shifted rightward along the *x*-axis to its proper location.

We may conclude that selection on Y-ending codons, represented by *φ*, detectable in dsDNA phages as the mean of the *b* values ([Table 3](#t3){ref-type="table"}) is significantly \>0 (mean *b* = 0.4622, SE = 0.07436, *t* = 6.2156, d.f. = 59, *P* \< 0.0001). When the false discovery method is used to control for type I error rate, about half of the *b* values are statistically significant at the 0.05 level ([Table 3](#t3){ref-type="table"}).

Other factors that may contribute to phage codon usage {#s14}
------------------------------------------------------

One factor that may contribute to phage codon usage bias is phage-encoded tRNA genes. Note that *E. coli* dsDNA phages other than those in clusters A, B, and C ([Figure 8](#fig8){ref-type="fig"}) are scattered in clusters that frequently have phage lineages with phage-encoded tRNA genes. The presence of phage-encoded tRNA genes can alter the host tRNA pool, so that *φ* may no longer reflect the selection on codon usage. For example, if the host tRNA for an NNY codon favors U-ending codons, but phage-encoded tRNA favors C-ending codons, then *φ* would not be a good predictor of phage codon usage. It is noteworthy that phage species in cluster A that do not have phage-encoded tRNA genes all have uniformly high *R* values, suggesting that selection by the host translation machinery may be more effective on phage codon usage in phages with no tRNA encoded in the phage genome. Alteration of the host tRNA pool through selective local tRNA enrichment in favor of viral gene translation has been documented in several viral species, including HIV-1 ([@bib61]) and vaccinia and influenza A ([@bib43]).

While almost all Y-ending codons are translated by tRNAs with a wobble G (except for Ile and Arg codon families where Y-ending codons are decoded by tRNAs with a wobble A chemically modified to inosine), different tRNAs with a wobble G appear to have different codon preferences, with some favoring C-ending codons, some favoring U-ending codons, and some with no detectable preference. At present, such a preference is not well understood and cannot be properly measured. This is in contrast to R-ending codon families that are typically translated by two types of tRNAs, one with a wobble U consistently preferring A-ending codons and the other with a wobble C consistently preferring G-ending codons. To properly assess the effect of phage-encoded tRNAs on phage codon usage, one needs minimally to assess whether the tRNA is actually functional and measure the synonymous codon preference of phage tRNAs. Currently we have no means of doing this bioinformatically.

Although our focus is on the joint effect of mutation and host tRNA-mediated selection on phage codon usage, we are aware of other factors that have been suggested to affect codon usage. Some bacterial hosts live in a high-temperature environment and have relatively high genomic GC. Their phages are also expected to have high GC to maintain genome stability at high temperature ([@bib75]). Different host species may have the 4 nt in quite different concentrations (*e.g.*, nucleotide C is typically rare and A typically abundant) and cytoplasmic parasites or symbionts such as virus or organelles should avoid using rare nucleotides in building their genome and RNA molecules ([@bib63]; [@bib74]; [@bib38]). Such avoidance would also be reflected in codon usage bias. However, these effects are likely additive and not expected to confound the relationships we aim to study here.

Dinucleotide frequencies are often used to explore the presence of site-dependent mutation patterns. Some bacterial species such as *Mycoplasma pulmonis* carry CpG-specific methyltransferase and exhibit strong CpG deficiency ([@bib65]) that could lead to context-dependent codon usage; *e.g.*, C-ending codons are particularly rare if the next codon is GNN (where N stands for any nucleotide). However, neither *E. coli* nor any of its phages carries CpG-specific methyltransferase genes. The ratio *P*~CpG~/(*P*~C~*P*~G~) is close to 1 for both *E. coli* and its phages.

Given nucleotide frequencies *P~i~* (where *i* = A, C, G, or T), the observed dinucleotide frequencies *P~ij~*, assuming random association, are expected to be *P~i~P~j~* (where *i*, *j* = A, C, G, or T). The deviation of the observed frequency from the expected frequency may be expressed as *D~ij~* = (*P~ij~* -- *P~i~P~j~*)/*P~i~P~j~*. A dinucleotide is in surplus if *D~ij~* \> 0 and in deficiency if *D~ij~* \< 0. The 11 ssDNA phages share similar dinucleotide frequencies, all with *D*~AA~ being the highest and *D*~AG~ the lowest. However, the surplus of AA dinucleotides is largely explainable by the usage of AAN codons, especially AAA codons, which is far more than expected. If we designate total number of codons in the 11 ssDNA phages as *N*~T~, then the expected number of AAA codons (designated *E*~AAA~) is 404.7 (= *N*~T~ × *P*~A~^3^) whereas the observed AAA codon (designated *O*~AAA~) is 930. *O*~AAG~, *O*~AAC~, and *O*~AAU~ codons are also greater than *E*~AAG~, *E*~AAC~ and *E*~AAU~, although not as dramatic as AAA codons. In contrast, the deficiency of AG dinucleotides can be attributed to *O*~AGA~, *O*~AGG~, *O*~AGC~, and *O*~AGT~ (= 104, 44, 163, and 179, respectively) all being smaller than *E*~AGA~, *E*~AGG~, *E*~AGC~, and *E*~AGT~ (= 351.5, 305.3, 306.9, and 430.6, respectively).

One may ask why *O*~AAN~ is far greater than *E*~AAN~ whereas *O*~AGN~ is far smaller than *E*~AGN~. One simple explanation invokes tRNA abundance ([@bib64]), which is well predicted by tRNA copy number ([@bib44]). In unicellular organisms such as *E. coli*, *Salmonella typhimurium*, and *Saccharomyces cerevisiae*, the frequency of an amino acid increases with the abundance of tRNA carrying the amino acid. As ssDNA phages do not carry their own tRNA genes and therefore depend entirely on the host tRNA pool to translate phage genes, we expect the frequency of amino acids to increase with the number of their cognate tRNA genes in the *E. coli* genome. This expectation is supported by a strong positive correlation between the frequency of an amino acid and the number of gene copies of the tRNA carrying the amino acid (*r* = 0.8426, *P* \< 0.0001, [Figure 12](#fig12){ref-type="fig"}). This finding explains why there is a surplus of AA dinucleotides because there are six tRNA genes decoding AAA and AAG codons, leading to relative overuse of AAR codons. In contrast, there are only two tRNA genes for AGR codons (coding Arg1) and only one tRNA for the AGY codons (coding for Ser1), which explains the relatively rarity of AGR and AGY codons as well as the deficiency of AG dinucleotides mentioned above.

![Frequency of amino acids increases with the number of tRNA genes in *E. coli* K12 MG1655. Met tRNAs are not included because Met is translated by both initiator and elongator tRNAs while all other amino acids are by elongator tRNAs only.](301fig12){#fig12}

We did not incorporate translation initiation into our model in this article. However, selection on codon adaptation is present only for mRNAs with efficient initiation; *i.e.*, one should expect little selection for codon adaptation in genes whose mRNAs have poor translation initiation efficiency. The availability of ribosomal loading data and their analysis ([@bib77]) may eventually lead to an index of translation initiation to facilitate more vigorous studies of codon adaptation.

We should mention that, although selection on codon usage appears more detectable in dsDNA phages than in ssDNA phages ([Table 2](#t2){ref-type="table"} and [Table 3](#t3){ref-type="table"}), this does not mean that dsDNA phage mRNA will necessarily have better translation efficiency than ssDNA phage mRNA. A previous study suggested that dsDNA phage coat protein genes have higher CAI ([@bib56]), which has been known to be a good measure of translation efficiency ([@bib11]; [@bib15]; [@bib10]), than ssDNA phage genes. However, there is no coat gene that is homologous between dsDNA and ssDNA phages, so the comparison may be between apples and oranges. If all genes are included, then the difference is minimal (mean CAI = 0.4768 for dsDNA phages and 0.4743 for ssDNA phages, excluding the 22 dsDNA phages with phage-encoded tRNA genes) and not statistically significant. The mean CAI value for the 22 dsDNA phages with phage-encoded tRNA genes is even lower, but that is because the phage-encoded tRNAs may allow the phage to use codons frequently used in the phage but rare in the host (*i.e.*, the phage tRNAs reduce the need for the phage to evolve codon usage similar to that of the host). Some ssDNA phages with increasing C→T mutation bias appear to increase the usage of codons in the "U-friendly" codon families, thereby achieving CAI values almost as large as those of dsDNA phages.

In summary, our results show that previous studies on phage codon adaptation are insufficient in at least two ways. First, codon frequencies from either all host CDSs or all highly expressed host genes are insufficient to capture the selection by host translation machinery. Second, it is crucially important to have explicit models to dissect the effect of mutation and selection. Our index (*φ*) is a proper measure of selection imposed by host translation machinery on phage codon usage, and our linear and nonlinear models allow us to estimate the C→T mutation bias (*B*~C→T~) and to evaluate the relative effect of the mutation bias and host translation machinery on phage codon usage. C→T mutations occur more frequently in ssDNA phages than in dsDNA phages and affect not only synonymous codon usage, but also nonsynonymous substitutions, especially in ssDNA phages. dsDNA phages exhibit better codon adaptation to host translation machinery than ssDNA phages, but much of the variation in codon usage may be attributed to phylogenetic inertia. Strand asymmetry strongly influences the efficiency of selection on codon adaptation and needs to be taken into account when studying codon adaptation.

We thank S. Aris-Brosou, N. Corradi, M. Gilchrist, A. Golshani, and J. Wang for discussion and comments. Two anonymous reviews and comments from J. Lawrence resulted in significant improvement of the manuscript. This study is supported by a Natural Science and Engineering Research Council of Canada discovery grant (to X.X.).

Communicating editor: J. Lawrence
