The Effect of Mutation and Selection on Codon Adaptation in Escherichia coli Bacteriophage

Studying phage codon adaptation is important not only for understanding the process of translation elongation, but also for reengineering phages for medical and industrial purposes. To evaluate the effect of mutation and selection on phage codon usage, we developed an index to measure selection imposed by host translation machinery, based on the difference in codon usage between all host genes and highly expressed host genes. We developed linear and nonlinear models to estimate the C→T mutation bias in different phage lineages and to evaluate the relative effect of mutation and host selection on phage codon usage. C→T-biased mutations occur more frequently in single-stranded DNA (ssDNA) phages than in double-stranded DNA (dsDNA) phages and affect not only synonymous codon usage, but also nonsynonymous substitutions at second codon positions, especially in ssDNA phages. The host translation machinery affects codon adaptation in both dsDNA and ssDNA phages, with a stronger effect on dsDNA phages than on ssDNA phages. Strand asymmetry with the associated local variation in mutation bias can significantly interfere with codon adaptation in both dsDNA and ssDNA phages.

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 Ranjan et al. 2007;Sau 2007;Abedon et al. 2011) or remove infectious biofilms (Azeredo and Sutherland 2008) 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 (Grosjean et al. 1978;Gouy 1987;Kunisawa et al. 1998;Sahu et al. 2005;Carbone 2008;Lucks et al. 2008) 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 (Duncan and Miller 1980;Lindahl 1993;Xia and Yuen 2005). In particular, the spontaneous deamination rate is 100 times higher in ssDNA than in double-stranded DNA (dsDNA) based on experimental evidence (Frederico et al. 1990), 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 (Umemura et al. 2002;Shackelton et al. 2005;Xia and Yuen 2005;Shackelton and Holmes 2006;Holmes 2008, 2009). Oxidative deamination leading to high C/U/T transitional mutation rates has been reported in ssDNA phage M13 (Kreutzer and Essigmann 1998).

The Effect of C/T Mutation Bias
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 proteincoding 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., 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., 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 P U values; e.g., P 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 P 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
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, 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 (Mitra et al. 1977(Mitra et al. , 1979Nasvall et al. 2007). 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 f as a proxy for U friendliness. The Ala codon family is U-friendly and has a positive f-value, whereas the Phe codon family is U-hostile and has a negative f-value ( Figure 1). f takes values between 21 and 1 and is equivalent to the Pearson correlation coefficient for continuous variables. Because f 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 .
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
The development of the f-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., Because P U.i can be readily computed from viral proteincoding genes, and f i can be derived from codon frequencies of host HEGs and non-HEGs (Figure 1), we can use Equation 3 to quantify the relative importance of mutation (B C/T ) and selection (f) on the codon usage bias of Y-ending codon families. If f 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 f), P U.i will become more dependent on bf i . Similarly, if B C/T becomes very large (i.e., very strong mutation bias), then bf 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 f i and approach 0 with decreasing f i . This implies a sigmoidal relationship between P U and f. For this reason, we have also fitted the following sigmoid function, 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 f = 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 f) 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, f 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 f 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 f 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.

Genomic data and processing
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 (Xia 2013b). All phage genomes were searched for encoded tRNAs by using the tRNAscan-SE Search Server (Schattner et al. 2005). The local TC skew plot, with the TC skew computed as (N T 2 N C )/(N T + N C ), Figure 1 Rationale of using the f-coefficient as a proxy for U friendliness, based on the codon frequencies (CF) between highly expressed genes (HEGs) and non-HEGs from E. coli. f can take values within the range between 21 and 1.
where N i is the number of nucleotides i along a moving window, was generated from DAMBE (Xia 2013b). All statistical analyses were done with SAS (SAS Institute 1994), 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 (Rice et al. 2000). 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
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) (Sharp and Li 1987) with the improved implementation (Xia 2007) and the effective number of codons (N c ) (Wright 1990) with the improved implementation (Sun et al. 2013). All these indexes were computed using DAMBE (Xia 2013b). 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
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 (Xu and Hao 2009) 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 k21 , a 2 a 3 . . . a k , and a 2 a 3 . . . a k21 . The CVTree method has been implemented in the most recent version of DAMBE (Xia 2013b). We used a k value of 5, which has been recommended for viral genomes (Xu and Hao 2009). 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
Codon preference by the E. coli translation machinery: f The f-values for E. coli Y-ending codons are generally small, ranging from 20.1512 to 0.1424 (Table 1). Among the 16 Y-ending codon families and subfamilies, 7 are U-friendly (with f . 0, with the mean = 0.0861) and 9 are U-hostile (with f , 0, with the mean = 20.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 Table 1 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 f as a measure of codon preference of the host translation machinery (a large f corresponds to greater preference of U-ending codons over C-ending codons) 7 U-friendly codon families, but will lead to deterioration in the 9 U-hostile codon families (Table 1).

Effect of mutation and selection on codon usage of E. coli ssDNA phages
If selection in favor of U-ending codons by the host translation machinery (f) is efficient, then we expect P U to increase with f. This expectation is consistent with data from E. coli Enterobacteria phage G4 (NC_001420) showing P U increasing roughly linearly with f ( Figure 2). Fitting the linear model in Equation 3 (Table 2). These would be the P U values when selection mediated by the host translation machinery is absent. With the slopes in Table 2, 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 P U values defined in Equation 2, which is true (Figure 3).
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 are all significantly greater than the observed 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 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). 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 f, i.e., large P U in U-friendly codon families (with large f-values) and small P U in U-hostile codon families (with small f-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 f. We therefore expect the correlation between P U and f (R) to decrease with increasing B C/T . This expectation is consistent with the Figure 2 Relationship between P U (the proportion of U-ending codons in Y-ending codon families) and f (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. Table 2 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 f for each phage species, and the statistical significance (two-tailed P) of the relationship  (Figure 4), although there is one outlying point (Enterobacteria phage If1; NC_001954) for which we offer an explanation later. The four ssDNA phage species (phage I2-2, IF1, Ike, and M13) with low correlation between P U and f (red diamonds in Figure 4) have codon usages significantly correlated with each other, which suggests that they might be phylogenetically related. The tree built with the CVTree algorithm (Xu and Hao 2009) implemented in DAMBE (Xia 2013b) does cluster these four species, all belonging to Inoviridae, into a monophyletic taxon ( Figure 5). Other viral proteomic trees (Rohwer and Edwards 2002;Edwards and Rohwer 2005) 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).
Because of the phylogenetic structure ( Figure 5), 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) 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) 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 and performed independent-contrasts analysis (Felsenstein 1985(Felsenstein , 2004 implemented in DAMBE (Xia 2013a, pp. 24-29;Xia 2013b) 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 (Benjamini and Hochberg 1995;Benjamini and Yekutieli 2001) 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 (Xia 2012b) and implemented in DAMBE (Xia 2013b).
Effect of mutation and selection, as well as evolutionary history, on codon usage of E. coli dsDNA phages Some dsDNA phages show strong response to selection by the host translation machinery (f), e.g., NC_010324 phage Phieco32, with P U strongly dependent on f ( Figure 6). The estimated B C/T spans a wide range (Table 3), 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). For example, phage BP-4795 (NC_004813), phage cdtI (NC_009514), and 10 other dsDNA phages have negative slopes between P U and f (Table 3), 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.
We noted that all phages that have negative correlation between P U and f 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), 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). Here we offer two explanations for the lack of codon adaptation in these 12 phages with empirical substantiation.
Phylogenetic inertia: 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 = 20.5638, Table 3). 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 (Bamford et al. 1995;Grahn et al. 2006). 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 Figure 5 Phylogenetic tree of ssDNA phages reconstructed by using the CVTree method (Xu and Hao 2009) implemented in Xia (2013b). 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 f), estimated B C/T , and number of tRNA genes (N_tRNA) in the phage genome. and will not lead to a negative b (Table 3). 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 (Bamford et al. 1995). 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 gramnegative 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 (Xu and Hao 2009) 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) 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 f (r = 0.9436). It is likely that the ancestor of podoviruses evolved good codon adaptation to an E. colilike 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.
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 f 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 Enav et al. 2012). 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, Table 3). Strong heterogeneity in codon adaptation among myoviruses is also visible among species under cluster F (Figure 8), 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 (Marin and Xia 2008;Xia 2012a,c).  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 f for each phage species, and the statistical significance (two-tailed P) of the relationship a Significant at the 0.05 level when experimentwise error rate is controlled by the false discovery rate method. b Significant with the more conservative approach of Benjamini and Yekutieli (2001). c Significant with the approach of Benjamini and Hochberg (1995).
Strand asymmetry, local mutation bias, and phage codon adaptation: 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 (20.6610, 20.6273, and 20.5957, respectively, Table 3) exhibit strong local TC skew, defined as (N T 2 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).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. 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). A similar response of nonsynonymous mutation rate to directional mutation pressure has also been documented in several other studies (Sueoka 1961;Lobry 2004;Urbina et al. 2006). 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). 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.
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 windowspecific skew values as an index of strand asymmetry (I SA ). Take the six dsDNA phage species in Figure 9, 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). 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 relationship between R and I SA in Figure 11 offers an explanation for the outlying point in Figure 4, 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). In short, the mutation effect (on the x-axis) for Enterobacteria phage If1 is underestimated in Figure 4, which does not take local strand asymmetry into consideration. When local strand asymmetry is accounted for (Figure 11), the point is shifted rightward along the x-axis to its proper location.
We may conclude that selection on Y-ending codons, represented by f, detectable in dsDNA phages as the mean Figure 7 Phage species that do not exhibit codon adaptation to E. coli (negative correlation between P U and f) 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.
of the b values (Table 3) 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).
Other factors that may contribute to phage codon usage 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) are scat-tered 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 f 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 f 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 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 (Xia and Yuen 2005).
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 (Xia 1996;Xia and Palidwor 2005;Marin and Xia 2008). 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 (Xia 2003) 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 3 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 Figure 9 TC skew, defined as (N T 2 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 f.  (= 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 (Xia 1998), which is well predicted by tRNA copy number (Percudani et al. 1997). 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). 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.
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 (Xia et al. 2011) 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 and Table 3), 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 (Sharp and Li 1987), which has been known to be a good measure of translation efficiency (Comeron and Aguade 1996;Duret and Mouchiroud 1999;Coghlan and Wolfe 2000), 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 (f) is a proper measure of selection imposed by host translation machinery on phage codon usage, and our linear and nonlinear models

Figure 11
The effect of selection on Y-ending codons, measured by the correlation (R) between P U and f, 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). 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.