Codon usage patterns in Nematoda: analysis based on over 25 million codons in thirty-two species

A codon usage table for 32 nematode species is presented and suggests that total genomic GC content drives codon usage.


Background
Utilization of the degenerate triplet code for amino acid (AA) translation is neither uniform nor random. In particular, there are distinct patterns among different species and genes. Such patterns can readily be characterized by codon usage, namely the observed percentage occurrence with which each codon is used to encode a given AA. This measure has direct utility in molecular characterization of a species in that it enables efficient degenerate and nondegenerate primer design for cross-species gene cloning, open reading frame determination, and optimal protein expression [1]. Such tools are particularly important with respect to species for which limited molecular information exists. Codon usage also serves as an indicator of molecular evolution [2]. Codon usage bias, namely the degree to which usage departs from uniform use of all available codons for an AA, can be influenced by a number of evolutionary processes. The guanine and cytosine (GC) versus adenine and thymine (AT) composition of the species' genome, probably the product of directional mutation pressure [3,4], is a key driver of both codon usage and AA composition [5,6]. Other factors that influence codon usage may include the relative abundance of isoaccepting tRNAs [7][8][9], especially for highly expressed mRNAs that require translational efficiency [10,11], presence of mRNA secondary structure [12,13], and facilitation of correct co-translational protein folding [14]. Codon usage appears not to be optimized to minimize the impact of errors in translation and replication [15].
Nematodes are a highly abundant and diverse group of organisms that exploit niches from free-living microbivory to plant and animal parasitism. Molecular phylogenies divide nematodes into five major named and numbered clades within which parasitism has arisen multiple times [16]: Dorylaimia (clade I), Enoplia (clade II), Spirurina (clade III), Tylenchina (clade IV), and Rhabditina (clade V). Following the sequencing of the complete genome of the model nematode Caenorhabditis elegans [17], we have begun to catalog the molecular diversity of nematode genomes through the generation of over 250,000 expressed sequence tags (ESTs) from more than 30 nematode species (including 28 parasites) in four clades. Gene expression analyses for several medically and economically important parasites such as filarial, hookworm, and root knot nematode species have been completed [18][19][20][21][22][23] (for reviews [24,25]). Moreover, we recently conducted a meta-analysis of partial genomes across the whole phylum with a focus on the conservation and diversification of encoded protein families [26]. Project information is maintained on several online resources [27][28][29][30]. Now, in the most extensive such study yet performed for any phylum, we extend the above analyses with a comprehensive survey of observed codon usage and bias based on nearly 26 million codons in 32 species of the Nematoda. Because of its completed genome, C. elegans has been the primary species utilized in nematode codon usage studies [31][32][33][34]. Our find-ings provide more complete information for Caenorhabditis based on all 41,782 currently predicted proteins in C. elegans and C. briggsae [35]. Studies for other nematode species have been more limited. Codon usage has been tabulated for a number of parasitic nematodes including filarial species Brugia malayi, Onchocerca volvulus, Wucheria bancrofti, Acanthocheilonema viteae, Dirofilaria immitis [36][37][38][39], Strongyloides stercoralis [40], Ascaris suum [41], Ancylostoma caninum, and Necator americanus [42]. Although Fadiel and coworkers [39] used up to 60 genes per species, sample sizes in the other studies were quite small, typically fewer than 10 representative genes and 5,000 codons per species. In the present study we used an average of 2,350 genes and 270,000 codons per species for the 30 non-Caenorhabditis species. Our results provide the first codon usage tables for 24 of these organisms. Web available automated codon usage databases compiled from GenBank [43] lack almost all of this information because they rely only on full-length protein coding gene sequence submissions rather than the EST data used here.
In analyzing codon distribution in Nematoda, we describe how average usage varies between species and across the phylum. For instance, it has been shown that there is a level of conservation in codon distribution between 'closely' related nematodes such as Brugia malayi and B. pahangi [37] and Brugia and Onchocerca [38]. These relationships do not appear to extend over greater evolutionary distances, for instance between Onchocerca and Caenorhabditis [36]. The evolutionary distance at which conservation of codon usage diminishes has not previously been established [32]. Here we show that codon usage similarity in Nematoda is a relatively short-range phenomenon, generally persisting over the breadth of a genus but then rapidly diminishing within each clade. We also show that the major factor affecting differences in mean codon usage between distantly related species is the coding sequence GC as compared with AT content. GC content also explains much of the observed variation in the effective number of codons, a measure of codon bias, and even differences in AA frequency.
predictions. To reduce noise derived from poor translations, our analysis considered only the longest open reading frame (ORF) translations with strong supporting evidence in the form of similarity to known or predicted proteins (BLASTX cutoff 1 × e -8 ) and retained only the polypeptide aligned portion of the nucleotide sequence. About 75% of the clusters met these criteria, yielding 8,080,057 codons originating from species other than Caenorhabditis, and 25,871,325 total codons from all 32 species including available predictions from C. elegans and C. briggsae. The 18 AA residues with redundant codons gave a total of (18) × C 32,2 = 496 comparisons of codon usage between species. Comprehensive tables of AA composition (Tables 2 and 3) and codon usage (Table 4) for all 32 Nematoda species studied are provided. Below we use these tables to examine, first, variation in AA composition and its relationship to GC content and, second, codon usage and its relationship to GC content.
To examine these variables independent of species relatedness, correlations were calculated using phylogenetically independent contrasts (see Materials and methods, below). The variances of the contrasts were computed for each character as a measure of the variance accumulating per unit branch length. The branch lengths were estimated from the maximum likelihood phylogeny assuming a molecular clock ( Figure 1); by this criterion, the tips of the tree are all equidistant in branch length from its root. Computed contrasts were plotted in all figures representing pair-wise comparisons, and the correlation coefficients were calculated from the paired contrasts. This method is robust to changes in molecular clock assumptions. (Trees calculated without the assumption of a molecular clock are similar in topology but differ in rooting, and branch lengths vary according to amount of base substitution in the 18S rRNA; the clock-based tree provides branch lengths that should estimate most closely the relative durations of branches in evolutionary time. Because independent contrasts are influenced mainly by relative branch lengths, our results should be robust to alternative placements of the root.)

Amino acid composition of nematode proteins and relationship to GC content
AA composition of predicted proteins in nematodes varies among species within a narrow window and is similar to that observed in other organisms (Tables 2 and 3). (Standard deviations in AA usage among nematodes range from 5% to 15% of mean usage, and mean nematode AA usage differs from the mean of four representative organisms by an average of 8%.) Across nematodes, Leu is the most common AA (8.8% of all codons) and Trp the least common (1.1%). Eight AAs contribute an average of more than 6% each to AA content (Ile, Gly, Val, Glu, Ala, Lys, Ser, and Leu); these AAs are also among the most common in the proteomes of other representative species, including humans ( Table 3). As in other taxa [46], nematodes show a correlation between AA usage and the degree of codon degeneracy (R = 0.72).
In nematodes, coding sequence GC content, derived from our EST clusters, varies from 32% to 51% (Table 1) among species, with a mean of 43.6 ± 5.9%. The distribution is biphasic, with a peak at 36% GC and a second peak at 48%. Strongyloides (SS and SR), Meloidogyne (MI, MJ, and so on), and filarial Table 2 Amino acid composition (%) of translations by nematode species

Clade Species
Amino acid Definitions of species two letter codes are provided in Table 1.
parasites (BM, DI, and OV) are the most AT rich (low GC); and NB, PP, and cyst nematodes (GP, GR, and HG) are the most GC rich (approximately 50%). The variation observed in AA composition among species shows a clear relationship to the species' coding sequence GC content. The frequency of AAs encoded by WWN codons (AA, AT, TA, or TT in the first and second nucleotide positions; Asn, Ile, Lys, Try, Phe, and Met) decreases with increasing coding sequence GC content (Figure 2a), whereas the proportion of AAs encoded by SSN codons (GG, GC, CG, and CC; Ala, Arg, Pro, and Gly) increases with higher coding sequence GC content (Figure 2b), and these relationships remain even after removing the effect of evolutionary relationships using phylogenetically independent contrasts. Among AAs, the most uniform and precipitous decrease with increasing GC content was seen with Ile and Tyr whereas the most uniform and rapid increase with higher GC content was seen with Ala and Arg. The trend is less pronounced for other AAs (flatter slope, lower R value). Thr, encoded by four GC/AT 'balanced' codons (ACN), exhibits no change in its frequency with changing GC content (data not shown).

Base composition by codon position in nematode transcripts and relationship to GC content
Codon usage in nematode species was examined by several methods, including comparison of base usage by position (1-3) over all AAs and comparison of codon usage within each AA. Over all AAs, purine (AG) and pyrimidine (TC) usage in positions 1, 2, and 3 is remarkably uniform between species, favoring purines in position 1 (AG 59.6 ± 1.5%), near equal usage in position 2 (AG 50.0 ± 0.8%), and pyrimidines in position 3 (AG 47.9 ± 1.5%; Additional data file 1). Similar values were observed in Schistosoma mansoni (AG 61%, 53%, and 48% in positions 1, 2, and 3, respectively) [1]. GC versus AT usage also varies by position but with much greater variance, with near equal usage in position 1 (50.3% GC) and lower GC usage in positions 2 and 3 (39.1 and 41.4%, respectively), mainly due to greater use of G in position 1 and T in positions 2 and 3 [4].
Additional file 1 Click here for file The variation observed in GC usage by codon position among species exhibits a clear relationship to the species' overall coding sequence GC content. Not surprisingly, both GC1 and GC2 composition increase with higher coding sequence GC3 content ( Figure 3). Specifically, species with high AT content like root-knot Meloidogyne species (MI, MJ, and so on) and filarial worms (BM, DI, and OV) [38,39] are biased toward codons terminating in A or T, whereas species with higher GC content such as NB, PP, cyst nematodes, and whipworms (TM and TV) prefer codons ending with G or C. Differences in calculated GC composition by codon position (1-3) between species are determined both by the species' AA usage (as described above) and the codons used for each AA. For example, Cys was encoded by TGT as much as 85% of the time for the AT-rich Strongyloides genomes, whereas TGC was used up to 60% of the time in GC-rich genomes such as NB, PP, and HG. To compare codon usage more systematically for individual AAs between species, we employed a statistical approach (described in Materials and Methods and in the following section).

Codon usage patterns and relationships to sampling method, nematode phylogeny, and GC content
Similarity in codon usage was quantified and reported as D 100 values for each species and AA compared [47,48] (matrix of D 100 values for each species and AA compared is available in Additional data file 2).
Additional file 2 Click here for file Because analyses of all but two of the nematode species were based on EST-derived partial genomes [26], comparisons were performed to estimate the differences in codon usage pattern that could be expected using EST collections versus gene predictions derived from a fully assembled and annotated genome. Using C. elegans, parallel analyses were performed using either all 22,254 predicted gene products or two EST datasets (CE-A and CE-B) each comprising 10,000 ESTs. Clustering and peptide predictions were performed using the same algorithms as for the other 30 species. The average D 100    Values are given as % per AA, or as numbers for Codons per AA. Definitions of species two letter codes are provided in Table 1. AA, amino acid. Phylogenetic analysis of changes in codon usage using (1antilog [-D]) × 100, interpretable as percentage divergence in overall codon usage (Figure 1), identifies five branches that have accumulated more than 5% change in codon usage. These branches are as follows: the most recent common ancestor of clades III, IVa, and IVb (5.2%); the most recent common ancestor of clade IVa (11.2%); the most recent common ancestor of genus Meloidogyne (6.7%); the most recent common ancestor of genus Globodera (7.3%); and the lineage represented by PP (8.3%). Genera Globodera, Meloidogyne, Pristionchus, and Strongyloides therefore represent the most highly derived patterns of codon usage in nematodes, with the remaining species exhibiting less relatively divergence from an ancestral nematode pattern.

Codon bias in nematode transcripts and relationship to GC content
We used the effective number of codons (ENC) to measure the degree of codon bias for a gene [49]. To ensure that differences in our available data for each species (for instance, cluster number and cluster length) were not creating artifacts in ENC values, quality checks were performed. Unlike measures such as codon bias index, scaled ×2, and intrinsic codon bias index, ENC values should be independent of translated polypeptide length and sample size [49,51], and our analysis confirmed this. No correlation with ENC was observed with either average translated polypeptide length or number of clusters for a species. In fact, SS and SR with the lowest ENC values had above average cluster length and number. As additional confirmation, we randomly selected 2,400 C. elegans genes (the average number of clusters for species other than CE and CB) and calculated ENC based on either full-length genes or genes trimmed to 121 AAs (the average length cluster translation for species other than CE and CB). Differences in the average ENC numbers for these datasets were not statistically significantly different from zero (P > 0.05).
In addition to codon bias, neighboring nucleotides influence the codon observed at a position relative to synonymous codons. The most important nucleotide determining such context dependent codon bias [52][53][54] is the first one following the codon (N1 context) [55,56]. An analysis using the complete genesets of Homo sapiens, Drosophila melanogaster, C. elegans, and Arabidopsis thaliana revealed that 90% of codons have a statistically significant N1 contextdependent codon bias [57]. Using the same method we calculated that, for the 30 nematode species represented by ESTderived codon data, an average of 63% of codons with N1 context have a statistically significant bias (because the R values differed from 1 by more than 3 standard deviations). Fedorov and colleagues [57] showed that their results were not considerably affected by gene sampling. However, for our dataset the calculated CE-A and CE-B N1 context with statistically significant bias was 75% and 83% of the codons, respectively, as compared with 96% when the complete C. elegans gene set was used. Therefore, the extent of significant N1 contextdependent codon bias determined from EST-based codon usage data may change as more complete nematode genomes become available.  [82]. Approximate percentage change in overall codon usage is indicated for five branches inferred to have undergone 5% or more divergence from an ancestral nematode pattern. This analysis identified genera Globodera, Meloidogyne, Pristionchus, and Strongyloides as having the most highly derived patterns of codon usage, and the remaining species as having relatively little net divergence from an ancestral nematode pattern. Definitions of species two letter codes are provided in Table 1

11.2%
percentage noncoding estimate. Although GC content varies across the genome for some organisms (for example, isochors in vertebrates [58]), GC content is fairly uniform across the C. elegans genome [17]; furthermore, as yet there is no evidence of non-uniformity in other nematode genomes. A positive correlation was observed between coding GC3 content and both total GC content and extrapolated noncoding GC content (R = 0.92; Figure 5). Noncoding genomic sequences varied across a wider span of GC values than did coding sequences. In all six nematodes, coding sequences were somewhat more GC rich than were noncoding sequences (2-10%).
Comparison of coding sequence GC versus 3'-untranslated region (UTR) GC also supports this conclusion. Calculated 3'-UTR GC for the 30 species in our EST dataset ranges from 28.6% to 46.1%. Correlation between phylogenetically independent contrasts of coding GC content (Table 1) and 3'-UTR has an R value of 0.81 (data not shown).

Codon usage patterns in abundantly expressed genes and candidate optimal codons
Representation in cDNA library generally correlates with abundance in the original biologic sample [59] although artifacts occur [60,61]. To investigate the difference in the codon usage patterns in highly abundant transcripts as compared with less abundantly expressed genes, as determined by Correlation between phylogenetically independent contrasts of coding sequence GC3 content and AA usage for 25 nematode species Figure 2 Correlation between phylogenetically independent contrasts of coding sequence GC3 content and AA usage for 25 nematode species. (a) AAs lysine (Lys), isoleucine (Ile), asparagine (Asn), and tyrosine (Tyr) are used less frequently as the species' coding sequence GC3 content increases. (b) AAs alanine (Ala), glycine (Gly), arginine (Arg), and proline (Pro) are used more frequently as the coding sequence GC3 content increases. AA, amino acid. Correlation between phylogenetically independent contrasts of the third position GC content (GC3) and that of the first (GC1) and second (GC2) codon positions for 25 nematode species Figure 3 Correlation between phylogenetically independent contrasts of the third position GC content (GC3) and that of the first (GC1) and second (GC2) codon positions for 25 nematode species.
Correlation between phylogenetically independent contrasts of each species' %GC3 and its mean ENC for 25 nematode species Figure 4 Correlation between phylogenetically independent contrasts of each species' %GC3 and its mean ENC for 25 nematode species. ENC, effective number of codons. Contrasts of %GC3 Contrasts of ENC

R=0.70
ESTs, we selected five species, each of which is a member of a different clade. The selected species (AY, MI, OV, SR, and TS) were represented by approximately 3,000 clusters each (range 2,693-3,214), and codon usage tables were generated for subsets of genes from each species: the 20 most abundant clusters versus all remaining clusters, and the 50 most abundant clusters versus all remaining clusters. Results of both comparisons were similar, and for simplicity we discuss only the results based on the comparison of the 50 most abundant versus all remaining clusters. Clusters 51 to about can be described as containing mainly genes with low to moderate expression because transcripts of extremely low abundance are less likely to be represented in EST collections (for instance, neuronal 7-transmembrane receptors). Codon usage tables, AA frequencies, and relative differences between AA usage of the most abundant and less abundant genes are available in Additional data file 5.
Additional file 5 Click here for file D values were calculated across all AAs and the codon usage in each species was generally similar for genes represented by abundant EST clusters and genes represented by low to moderate expression EST clusters. SR exhibited the greatest difference between the two usage patterns (D 100 = 6.15). Additionally, for all the species at least seven AAs were used significantly more frequently in the abundant genes than in the remainder of the genes. For example, although the abundant OV clusters had a Pro composition of 10.5% of all AAs, the rest of the clusters were only 4.4% Pro.
Examining the codon usage frequencies within an AA, an increase in usage has been noted with higher gene expression for specific so-called 'optimal' codons [62,63]. Using the codon usage tables for the top 50 and remaining clusters, we have defined a list of potentially optimal codons with usage that is higher in abundant transcripts by a statistically significant measure. Out of the 59 synonymous codons there were 24, 28, 25, 27, and 23 candidate optimal codons (Table 5) in AY, MI, OV, SR, and TS, respectively. For example, Tyr is encoded by two codons (TAC and TAT); in AY TAC is used 75% of the time in the abundant clusters and 59% of the time in the less abundant clusters. Similar analysis documented about 21 candidate optimal codons in C. elegans for which usage differed significantly when comparing high and low expressed genes [31,33,64]. Confirmation of these candidate codons as truly 'optimal' will require additional investigations, including other means of verifying relative expression levels (for example, microarrays and reverse transcription [RT]-polymerase chain reaction [PCR]).

Discussion
A comprehensive and well supported codon usage table for 32 nematode species across most of the phylum's major clades and based on nearly 26 million codons is now available. Use of large EST datasets provide an excellent resource for determining a species mean codon usage with results that differ only modestly from those obtained from full genomes. In nematodes, codon usage varies widely, as does coding and noncoding GC content of nematode genomes. GC content correlates with AA usage, similarity of codon usage, and codon bias. Codon usage similarity in Nematoda usually persists within a genus but then rapidly diminishes, even within each major clade (clades I-V). Based on EST sampling, differences in codon usage between highly abundant genes and moderately expressed genes are recognizable, and candidate optimal codons can be identified.

GC content, causality, and directional mutation pressure
Correlations between GC content and mean codon usage and mean AA usage similar to those we describe across the phylum Nematoda have been observed in many other species [4,[65][66][67][68][69][70]. Directional mutation pressure is a theory proposed to quantify differences in GC content observed in species [3]. Important variables include the relative values of the mutation rates u (GC/CG → AT/TA change) and v (AT/TA → GC/ CG change). The preponderance of the evidence supports causality of genome GC content, as determined by directional mutation pressure or nucleotide level selective pressure, driving both codon usage and AA composition rather than the reverse. First, in an examination of sequence data from a large number of a bacteria, archaea, and eukaryotes, a model assuming directional mutation and selection at the nucleotide level with different rates of change for each of the three codon positions can explain 71-87% of the variance in codon usage Correlation between coding sequence (transcriptome) %GC3 and genome %GC for six nematode species with extensive available genomic sequence Figure 5 Correlation between coding sequence (transcriptome) %GC3 and genome %GC for six nematode species with extensive available genomic sequence. The green line indicates the coding sequence %GC versus the full genomic %GC. In this case, coding sequence %GC3 is a contributor to the full genomic %GC such that X and Y are not independent variables. The red line indicates the coding sequence %GC3 versus noncoding genomic %GC. In this case, the coding sequence contribution has been removed from genomic totals such that X and Y are independent variables. For BM, TS, HC, and AC, the calculation of noncoding genomic %GC relies on the assumption that the species will have a similar breakdown of coding and noncoding sequence as CE. Assembly and gene calling for the BM, HC, TS, and AC sequences are needed to test this assumption. Definitions of species two letter codes are provided in Table 1. Coding %GC3 Genomic %GC

R=0.92
and 71-79% of the variance in AA composition [5]. Knight and coworkers [5] found that between species an AA's change in frequency in response to GC content is determined by the mean GC content of its codons, whereas a codon's change in frequency is determined by the difference between its GC content and the mean GC content of its synonyms. We observe this result to be generally true across nematodes as well.
Second, an analysis comparing codon usage from eubacterial and archaeal species with complete genomes [6] found that codon usage can be predicted with some accuracy if one knows only the sequence of the species' intergenic sequences from which genome GC content, and context dependent nucleotide bias parameters can be calculated. Using data from six nematode species for which substantial genome sequence data are available, we observed that coding sequence GC3 content correlates with noncoding sequence GC content. This perhaps indicates that, for nematodes too, it should be possible to predict mean codon usage using only knowledge of the intergenic sequences of the species. Our findings are consistent with the model that genome GC content drives both mean codon usage and AA composition.
Little is known about why directional mutation pressure or selective pressure leads to differences in genomic GC content among species [5,6]. Within nematodes we see no pattern based on ecologic niche or other factors that corresponds to GC content. For instance, cyst nematodes (GP, GR, and HG) . Whatever the driving forces, it is important for nematologists to note that they are sufficiently strong not only to change base composition in wobble sites (third position) but also to alter first and second codon positions and even AA sequences -features that are sometimes assumed to be under selective pressure for conservation.

Species' mean codon usage versus optimal codons in abundantly expressed genes
Our use of thousands of genes per species without weighting for abundance of expression has produced a codon usage dataset that probably reflects codon usage for genes with low to moderate abundance of mRNAs. In the case of C. elegans and C. briggsae, our codon usage table reflects the mean of all predicted genes, although this is similar to that observed based on sampling of 10,000 ESTs. At this 'genome-wide' level, genome GC content is a dominant factor. However, codon usage within a species does vary from gene to gene.
Prior studies of C. elegans codon usage have examined codon usage and the role of 'optimal codons' in putatively abundantly expressed genes [31,33,64]. Stenico and coworkers [31] observed differences between usage of specific codons based on 168 known genes, including many highly expressed transcripts (for instance, actin, myosin, collagen, and vitellogenin), and 90 unidentified reading frames (URFs) emerging from sequencing efforts presumed to represent a more 'random sampling' of the genome. Overall, our codon usage results based on the full C. elegans genome are similar to both the results from Stenico and coworkers' 168 known genes (D 100 = 0.97) and the 90 URFs (D 100 = 1.1). Duret and coworkers [33] weighted 15,425 C. elegans genes for expression levels based on their EST abundance and identified 21 favored codons used most frequently in highly expressed genes. In all cases, these optimal codons could be decoded by isoaccepting tRNAs that had the highest gene copy number in the genome, indicating that optimal codons are probably selected for translational efficiency. Likewise, Kanaya and coworkers [64] showed that, in C. elegans, ribosomal proteins and histones, selected as representatives of highly expressed genes, also use optimal codons different from those used by average genes and that these optimal codons correspond to tRNA gene copy number. AA frequencies in abundant C. elegans genes also correspond to isoaccepting tRNA gene copy number (R 2 = 0.67) [33].
Therefore, in C. elegans different pictures emerge of evolutionary forces acting on codons and AAs in low to moderately expressed genes (directional mutation pressure, genome GC content) compared with abundantly expressed genes (optimal codons, tRNA copy number). In other nematodes, it is possible that a similar dichotomy exists, although we cur-rently lack knowledge of tRNA gene copy number, and information on gene expression levels is largely limited to estimations based on EST abundance. Here, we have provided candidate optimal codons in AY, MI, OV, SR, and TS. A more detailed examination of codon usage as it relates to gene expression level in other nematodes will be possible by taking advantage of microarray and RF-PCR confirmation of transcript abundance.

Implications for phylogenetic studies and molecular biology
The extent to which average nematode genes sequences are susceptible to GC or AT shifts should sound a cautionary note for phylogenetic studies of nematode species, genes, and proteins based solely on coding sequences because convergent evolution may create confusing results. Knight and coworkers [5] noted that, 'Pairs of species with convergent GC contents might also evolve convergent protein sequences, especially at functionally unconstrained positions. For example, the frequencies of both lysine and arginine are highly (but oppositely) correlated with GC content, and lysine and arginine can easily substitute for one another in proteins.' In nematodes as well, one can envision exchanges of Lys and Arg ( Figure 2).
For cloning genes of interest from various nematode species, we found that codon usage is a rapidly evolving feature such that codon usage patterns beyond within a genus comparisons are often divergent. Therefore, extrapolating assumed codon usage patterns to unsampled species in nematodes beyond the genus level is unlikely to be successful. At a practical level of species choice, cloning of orthologs and homologs of interest from species that are AT rich and have low ENC values, such as SS and MI with low ENC values, will require fewer degenerate primers than may be needed for more GC rich species such as TC and MI. Transcript abundance is also an important factor because genes suspected of high level expression are likely to exhibit a shift in their codon usage from the species average toward optimal codons selected for translational efficiency.

Conclusion
Extensive sequence datasets from one complete, one draft, and 30 partial genomes across the phylum Nematoda have been used to analyze the conservation and diversification of encoded protein families [26] and the factors effecting codon usage and bias (the present report). The undertaken comprehensive survey of observed codon usage and bias is based on 26 million codons in 32 species, making it the most extensive study for any phylum. Our data indicate that similarity between species in average codon usage is a short range phenomenon, generally rapidly diminishing beyond the genus level. Mapping codon usage changes to the phyla indicates the genera Globodera, Meloidogyne, Pristionchus, and Strongyloides have the most highly derived patterns of codon usage in nematodes, with the remaining species exhibiting less relatively divergence from an ancestral nematode pattern. There was a strong correlation between the exonic GC content and similarity in codon usage. GC content also explains much of the observed variation in the effective number of codons, a measure of codon bias, and even differences in AA frequency. Results from partial genomes assembled from ESTs and complete genomes provide generally good agreement on codon usage, although refinement will be necessary as more sequences become available. EST collections from five species have also been used as a starting point to identify potentially abundant genes and predict optimal codons. These predictions will also be refined using more accurate measures of gene expression, including microarrays and quantitative RT-PCR.

Sequence acquisition and organization
To perform the first meta-analysis of the genomic biology of the phylum Nematoda [26], ESTs from 30 nematode species generated by our laboratories and others were downloaded from the dbEST division of GenBank in May 2003. For consistency, in this accompanying analysis of codon usage we used this dataset for all analyses. Sequences were collated and processed into partial genomes using the PartiGene pipeline [71,72]. Polypeptide translations were predicted using prot4EST [45,72]. Wormpep_dna121 (March 2004; Welcome Trust Sanger Institute, unpublished data) was used for C. elegans analysis, and the hybrid gene set [35] was used for C. briggsae analysis. Mitochondria can have codon usage differing from that of the nuclear genome, and therefore protein coding genes from mitochondrial genomes were eliminated from consideration. Codon usage tables for human, Saccharomyces cerevisiae, and Escherichia coli were derived from the Codon Usage Table Database [73] derived from GenBank Release 140.0 (22 March 2004 [74]).

Phylogenetic correction
Analyses of the relationships among GC content, AA, and codon usage values require statistical correction for the phylogenetic relatedness of the species being studied using phylogenetically independent contrasts [75]. To generate these contrasts, we performed the following procedures. First, to construct a phylogenetic tree independent of the transcriptomic data analyzed in this paper, we aligned 18S ribosomal RNA sequences using Clustalw [76] for all nematodes for which more than 15 kilobases of sequence was available. The 18S sequence GenBank accession numbers are available in Figure 1; the sequences from a priapulid worm and a nematomorph were used as outgroups [16] but excluded from our analysis. Alignments were trimmed to reflect only the overlapping portion of the sequences from all species analyzed. Second, this alignment, containing 1,841 base pairs/species (including gaps) and an alternative alignment excluding any region involving an insertion or deletion (1,423 base pairs/ species remained), was used to estimate phylogenies from the nucleotide sequences by parsimony and maximum likelihood (with and without assumption of a molecular clock) using Phylip [77]. Third, the trees with branch lengths derived from molecular clock-based analysis were used to calculate phylogenetically independent contrasts for our parameters of interest [75]. The Phylip program 'contrasts' was used to compute the phylogenetically independent contrasts using a Brownian-motion model [78,79] of genomic evolution.

Bioinformatics
The Emboss program 'cusp' was used to calculate codon usage in the predicted translations [80]. The ENC [49] was calculated using the Emboss program 'Codon Heterozygosity (Inverse of) in Protein-coding Sequences'. A genetic distance statistic was used to quantify divergence of synonymous codon usage between species [47] follows. Let t j be the number of codons that code for the j th amino acid. We omit analysis of the nondegenerate codons Met and Trp, as well as the 'stop' codon, so that j = 1, 2 ... r, where r = 18. Furthermore, let a ij and b ij be the frequencies of the i th synonymous codon in the j th AA of two organisms A and B, respectively. Then, Nei's difference statistic D is defined as the following: Investigators have used D as an empirical measure of difference, averaged over all r residues, of the codon usage between organisms [48]. There are a total of C 32,2 = 496 meaningful comparisons for the entire collection of 32 species. These results are presented as an N × N square matrix and the values are presented as D × 100. For simplicity in the remainder of the text we will refer to D × 100 as D 100 .
Phylogenetic changes in codon usage were analyzed using the species tree derived from 18S rRNA sequences estimated by maximum likelihood with a molecular clock imposed. Partitioning a matrix of distance values on a phylogenetic tree can estimate amounts of change occurring on each branch, provided that the distance metrics obey the triangle inequality (see the discussion on page 25 of the report by Page and Holmes [81]). Because of its logarithmic operation, Nei's difference statistic D violates the triangle inequality at high values.
For the phylogenetic analysis of codon usage, we substituted for D a distance measure equal to 1 -antilog(-D), which obeys the triangle inequality. Distances were partitioned on the tree topology using the cyclic neighbor-joining algorithm illus- trated by Avise [82], except that the topology was specified by the prior analysis of 18S rRNA sequences.
The ENC was used to measure the degree of codon bias for a gene [49]. Because the ENC statistic is not reliable when analyzing very short sequences (20 AAs or less), 54 short translations out of a total of 70,358 were discarded from these analyses. The relative abundance of nematode codons (per species) having a statistically significant N1 context-dependent codon bias was calculated by computing the R values and the standard deviations, as described by Fedorov and coworkers [57].
Predicted expression level of a transcript (abundant, moderate, rare, and so on) was determined by counting the number of ESTs comprising the cluster. Five species from different clades that have been sampled with at least 10,000 ESTs from several life-stage libraries were selected for these analyses.
Most of the cDNA libraries were constructed using the same protocols [61,83], and although the libraries generally correlate with abundance in the original biologic sample, artifacts can occur. The increase in use of a given codon for an AA in highly expressed genes (optimal codons) was considered significant when the difference of the codon distributions within that AA was statistically significant between datasets (P ≤ 0.01).
To assess the differences in calculated codon usage distributions when using partial (EST-based) as compared with whole genome data, we generated two datasets using C. elegans ESTs and compared them with the curated gene set of C. elegans (Wormpep version 121). Each EST dataset was composed of 10,000 ESTs (approximately the average number of ESTs used for the other 30 species); clustering and peptide predictions were performed using the same algorithms as for the other species.

Additional data files
The following additional data are included with the online version of this article: An Excel file containing a table that shows the nucleotide usage (%) by codon position and nematode species (Additional data file 1); an Excel file containing an N × N square matrix that shows codon usage across degenerate AAs for 25 nematode species reported as D values (Additional data file 2); a PowerPoint file containing a figure that shows the distribution of genes with various degrees of codon usage bias as measured by ENC for three species with approximately the same number of clusters but with different coding GC content (S.