Relation between Codon Bias and Protein Abundance

To test the idea that ribosome load is the main driving factor for biased codon usage [14], we developed an evolution model (Eq. (6) in Methods) that describes the competition between mutation among synonymous codons and selection for the fast codons. The selection is defined [by Eq. (5)] as a dependence of the growth rate on the translation rates of individual codons (or the time a ribosome spends at each codon) [30]. The fitness cost of a slow codon is incurred every time this codon is read by a ribosome, and the overall fitness cost is proportional to the abundance of the protein encoded by the genetic sequence to which this codon belongs. The model thus describes a fitness landscape that decreases linearly from a single peak (sometimes called a Mount-Fuji landscape).

Before we discuss the relation between codon usage and protein abundance that results from our model, we note that, in principle, codon usage and the choice of the preferred codons (or the tRNA concentrations, which are a major factor determining that choice) co-evolve [30]. Our model does not account for such co-evolution, based on the following rationale: the shift from one preferred codon to another would affect many proteins and thus carries a much bigger fitness cost than individual synonymous mutations. We therefore expect codon preferences (and tRNA concentrations) to evolve on much longer timescales than the codon usage, so that the evolution of codon usage occurs essentially for fixed codon preferences. A similar argument has recently been supported by the observation that tRNA gene copy numbers (used as proxies for tRNA concentrations) are strongly correlated amongst different yeast species [33]. Here we analyzed the copy numbers of tRNA genes of 12 enteric and closely related bacteria. We found them to be highly correlated as well (Pearson correlation coefficients >0.7, Figure 1 A and B). Nevertheless, the tRNA gene copy numbers display a surprising amount of variability between closely related strains. We therefore considered a second measure of codon preferences, the ‘preferred codons’ for each amino acid as defined by Hershberg and Petrov [34], which are identified as those that get most strongly enriched with increasing codon bias among all genes of a species. For most of the bacteria considered here, the preferred codons are the same for all 18 amino acids despite a large fraction of synonymous mutations (Figure 1C). Differences in the preferred codons are only seen for the most distantly related species (Figure 1 C and D), which also exhibit the largest frequency of synonymous and non-synonymous mutations. These results support the general picture that the usage of synonymous codons evolves on a shorter timescale than the assignment of preferred codons.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Evolution of tRNA concentrations and assignment of preferred codons.

A: Correlation of tRNA gene copy numbers between E. coli MG1655 and S. typhimurium. For tRNAs with different gene copy numbers in the two organisms, the anticodon is indicated. B: Correlation coefficients of tRNA gene copy numbers of 11 enteric and related bacteria relative to E. coli MG1655. C: Number of unchanged ‘preferred codons’ as defined by Hershberg and Petrov [34] as a function of the frequency of synonymous mutations for the same 11 bacteria relative to E. coli MG1655. D: Phylogram obtained from an alignment of the rpoB sequences with ClustalW [60].

https://doi.org/10.1371/journal.pone.0048542.g001

As the translation rates are known only for a subset of codons, we classify all codons as either ‘fast’ or ‘slow’ according to their use in ribosomal protein genes (Table S1). This classification is clearly a simplification, but has the advantage that only genomic sequence information is required and no additional information such as cellular tRNA concentrations or measured translation rates for individual codons. This classification is consistent with the measured translation rates, where these are known, and it agrees with other classification schemes used in the literature (see Table S1 and Methods). Based on this classification, the fraction of slow codons in the nucleotide sequence of a gene arises as a natural measure of codon bias. displays a simple dependence on the abundance of the protein (N_p) encoded by that gene, that can be described by the interpolation formula.(1)

This dependence exhibits two regimes: For weak selection (low protein abundance), synonymous codons are used randomly, and slow codons occur with frequency r₀, while for strong selection (protein abundance higher than a threshold abundance N₀), the frequency of slow codons depends inversely on the protein abundance.

Even though it is known, qualitatively, for a long time that rare codons are particularly rare in abundant proteins [4], [19], we can now quantitatively compare the predicted dependence to genomic [35] and proteomic data [36], [37], [38], [39], [40] using E. coli as a model system. We analyze the relation between and both at (i) the full protein sequence level using individual protein abundance and (ii) individual amino acid level (the next section) by examining the usage of synonymous codons for each amino acid across all measureable proteins in the organism.

We compared the prediction from our model with five different protein expression data sets (Figures 2 and S1, Table 1). Figure 2 shows the comparison using the largest protein abundance data of Ishihama et al. [36]. Each dot in the plot represents one protein (through N_p and r_p), the red points show binned data, i.e. averages of r_p over groups of proteins with similar abundance N_p (within an up to 3-fold range of protein abundance), and the solid green line shows the prediction from the model for the best fit of the parameters and . There is good agreement of the binned data with the model especially for low to medium protein abundance. At the same time individual proteins scatter widely (most within two-fold as marked by the dashed lines) around the predicted behavior. The deviation between the binned data and the model for very abundant proteins (N_p>10000) could be specific to the abundance data of Ishihama et al., as it is not obvious in the other abundance data sets, which otherwise exhibit rather similar behavior (Figure S1). Abundance values for high-abundance proteins are generally somewhat uncertain, as indicted by the wide range of abundance numbers for different ribosomal proteins (which should have the same abundance in the cell), but this scatter is particularly pronounced in the Ishihama et al. data set. The protein abundance at which strong selection sets in, characterized by the parameter in Eq. (1) is found to be a few thousand molecules per cell (N₀≈5030 with the data of Ishihama et al.; the precise value depends on which protein abundance data set we use, see Table 1). This value of N₀ falls within the range obtained from microscopic parameter estimates (see Methods). For low protein abundance, our fit leads to . This is lower than the expected value based on fully random usage of synonymous codons (, see Methods), providing a first indication that not all codon bias is captured by the ribosome load hypothesis.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Data sets for protein abundance and parameters of fit with the abundance-codon bias relation from Eq. (6).

https://doi.org/10.1371/journal.pone.0048542.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Relation between codon bias and protein abundance.

Black dots show individual proteins with abundance data (N_p, protein copy number per cell) taken from the data set of Ishihama et al. [36]. The red points depict averages over proteins of similar abundance (see the Supporting Text for a discussion of the error bars). The solid green line is fit of the model (Eq. 7) for individual proteins in the data set. The dashed green lines enclose a two-fold deviation from the fit. Codon bias is measured by the fraction of slow codons () in the protein coding sequence accounting for all amino acids except Met, Trp, Glu, and Lys.

https://doi.org/10.1371/journal.pone.0048542.g002

For comparison, we also used a gene expression data set from an RNA microarray experiment (from ref. [41], data for rich medium) and analyzed it in the same way (Figure S2). Using mRNA levels as a proxy for protein abundance has the disadvantage that it implicitly assumes that all genes are translated at the same rate, while in reality there is a broad distribution of the ratio of proteins per mRNA [37]. Indeed the main difference to the results with proteomic data (Figs. 2 and S1) is the even greater scatter of the fraction of slow codons for highly expressed genes.

Analysis for Individual Amino Acids

We also performed the same analysis for each individual amino acid encoded by more than one codon (i.e. excluding Trp and Met). Figs. 3A and B show two examples, Asn and Lys (discussed below); results for all amino acids are shown in Figure S3. For every amino acid, we determined and (Fig. 3C–F) by fitting with Eq. (1) using the two abundance data sets of Ishihama et al. [36] and Lu et al. [37]. We compared the values of (which measures the bias of codon usage for low abundance proteins, i.e. in the weak selection limit of our model) as obtained from the fit with the expectations from random usage of synonymous codons. The expected and observed values are strongly correlated, but one can see that the values from the fit tend to be smaller than expected (Figure 3 C and E). This indicates that codon usage is not fully random even for low-abundance proteins.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Analysis of individual amino acids for codon usage bias and protein abundance.

A and B: Fraction of slow codons (r_p) for Asn and Lys plotted against the protein abundance [36]. As in Figure 2, black dots show data for individual proteins; red dots show averages over proteins of similar abundance; the solid green line is a fit based on Eq. 7. Dashed green lines mark the two-fold range from the model to the measured protein abundance. C: Values of r₀, the parameter characterizing codon usage bias in the weak-selection limit plotted versus the expected values based on the assumption of complete random codon usage. Red dots are for amino acids encoded by two synonymous codons, black dots for those with more than two synonymous codons. D: Histograms of values obtained for the threshold parameter N₀ that characterizes the protein abundance above which strong selection sets in. Abundance data in (C) and (D) are from Ishihama et al. [36]. (E) and (F) contain the same measures as (C) and (D) using abundance data from Lu et al. [37].

https://doi.org/10.1371/journal.pone.0048542.g003

The parameter indicates the threshold in protein abundance where selection becomes dominant. For most amino acids, is found to be in the range of a few 1000 proteins/cell (Figure 3 D and F). The fact the values for are similar for all amino acids provides a rationale for pooling different amino acids together. An exception are the amino acids glutamate and lysine, for which we find essentially no or only very weak abundance-dependence of codon bias (Figures S3A and 3B), corresponding to much larger N₀ than for other amino acids (Figure 3 D and F). Both however exhibit abundance-independent codon bias, with values of 0.28 (Glu) and 0.22 (Lys). For this reason, Glu and Lys were excluded in the calculations of the fraction of slow codons () of Figure 2 and Figure S1.

Effects of Sequence Length

In the analysis shown so far, we have excluded the first 50 amino acids of each protein sequence because several studies have found that the initial part of coding sequences is enriched in rare codons [42], [43], a phenomenon attributed to either a requirement for weak secondary structures in the mRNA [23] or the need for a slow on-ramp to limit the maximal translation rate and prevent ribosome traffic jams [43], [44]. Our analysis confirms this observation: We find that for most genes, the fraction of slow codons () is slightly (<2-fold) higher in the initial sequence than in the rest of the sequence (Figure 4A). Overall, we do not see a dependence of the fraction of slow codons (r_p) on the length of the protein sequence, both for full sequences and for sequences without the first 50 amino acids (correlation coefficients R = −0.01 and 0.02, respectively, see Figure 4B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Analysis of sequence-length dependent effects.

A: Fraction of slow codons () in the initial part of a gene sequence (first 50 amino acid) vs. fraction of slow codons in the rest of the sequence. Inset: Histogram of the ratio of these two fractions of slow codons. B: Dependence of the fraction of slow codons on the sequence length with (black) and without (red) the initial 50 amino acids. C and D: Length-dependence of the mean square deviation of the codon usage in a protein sequence from expected values based on the fitted average, using the abundance data of Ishihama et al. [36] (C) and Lu et al. [37] (D).

https://doi.org/10.1371/journal.pone.0048542.g004

The observation that binning the proteins with similar abundance reduces the scatter of the data and improves the agreement of the data with the model suggests that such scatter in individual proteins may simply be due to insufficient sampling as the sequences of individual proteins may be too short. Accordingly, the deviation of the observed fraction of slow codons () from the prediction should be lower for larger proteins (longer amino acid sequences). We proceeded to calculate the mean square deviation between the data and model prediction for each sequence and tested for correlations with the sequence length (Figure 4C and D). However, we found only very weak negative correlation (correlation coefficients R = −0.25 and −0.03 for the data sets of Ishihama et al. [36] and Lu et al. [37], respectively). We therefore conclude that the deviations of the fraction of slow codons () observed in individual proteins from the predictions of our model are not random variations due to limited sampling, but rather reflect some other properties of these proteins or their mRNA not captured by our model (see also the calculation of error estimates in Methods). Nevertheless, these deviations are averaged out if proteins of similar abundance are pooled together.

Predicting Protein Abundance

Finally, the relation between the frequency of slow codons and protein abundance that results from the evolution model can be used to predict protein abundance from its gene sequence. Our analysis suggests the following procedure: First, codons have to be classified as either fast or slow according to their usage in ribosomal protein genes as described in Methods. No information beyond sequence information (such as tRNA concentrations) is needed. Second, the fraction of slow codons, , in each protein is determined. The predicted abundance of the protein can then be obtained from rewriting Eq. (1) as.(2)provided that estimates for and are known. However , for which only a relatively broad range can be estimated (see Methods) is only required for absolute protein copy numbers, as gives protein abundance relative to the abundance threshold for strong selection. r₀ can be estimated from random codon usage, weighted with the amino acid frequencies (for E. coli, this leads to , see Methods). We predicted the abundance for all E. coli proteins and compare them with the abundance data from Ishihama et al. [36] and Lu et al. [37] in Figure 5A and B. The two quantities are clearly correlated with correlation coefficient R = 0.66 and R = 0.54, respectively (calculated excluding all proteins for which the predicted abundance values are negative or infinite). The correlation is comparable to the correlations among different published abundance data sets (e.g., Ishihama et al. – Lu et al.: R = 0.62, Ishihama et al. – Pedersen et al.: R = 0.69, see Table S2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Predicting protein abundance from sequence data.

Predicted protein abundance (relative to the selection threshold N₀) as obtained from Eq. 2, vs. the measured abundance data of Ishihama et al. [36] (A) and Lu et al. [37] (B).

https://doi.org/10.1371/journal.pone.0048542.g005

Fitness Landscape

In a fast growing unicellular organism, a characteristic fitness measure is the growth rate μ. The effect of translation speed on μ, can be calculated by the following argument given first by Bulmer [30], [31]: If ribosome load is growth limiting in the sense that essentially all ribosomes are translating all the time, then the doubling time τ of the cells is given by the time these ribosomes need to double all protein components of the cell,(3)

In this equation, the sum runs over all proteins in the cell. T_p is the time it takes to translate one single protein p, N_p is the copy number of this protein, and N_Rb is the number of ribosomes per cell. (If there is extensive turnover of a protein, the number of copies of that protein that are synthesized during one cell doubling can be considerably larger than the abundance of the protein in the cell. In that case, should be interpreted as the number of copies synthesized during a cell doubling or as the synthesis rate. However, under conditions of rapid cell growth that are relevant here, proteins are typically stable in E. coli [50], [51].) We assume that every change from a fast codon to a slow codon increases T_p by a small time Δt, which in principle can be codon-dependent (as, e.g., in Ref. [32]). As will be justified below, we will make the simplification that all codons can be classified as either fast or slow and that there is only one such translation time difference.

We now consider a single protein and write the doubling time τ as a function of the number n of slow codons in this specific protein as(4)

The growth rate is(5)where τ₀ and μ₀ = ln2/τ₀ are the doubling time and growth rate in the absence of slow codons for the protein under consideration and is defined as the selection coefficient. The approximation in Eq. (5) requires (shown to be the case below) and leads to a linear, Mount Fuji-type fitness landscape with a selection coefficient, s, that is proportional to the abundance of the protein under consideration.

Evolution Model

We now use the expression for fitness given in Eq. (5) in an evolution equation [52] for a population of cells with different codon usage for a particular protein sequence. Assuming all codons are either fast or slow, we write the evolution equation for the number of cells in that population that have n slow codons in the sequence of a particular protein p, N(n), as(6)

The first term on the right hand side describes growth, the second and third term describe beneficial and harmful mutations by decreasing and increasing n, the number of slow codons, respectively. ν₀ is the synonymous mutation rate per codon. C is the total number of codons encoding the amino acid under consideration, so C-1 is the number of possible mutations. C_fast and C_slow = C- C_fast are the number of fast and slow codons. L_p is the total number of codons (length of the protein p).

Several comments are in order: (i) In Eq. (6), N(n) is not normalized to a fixed population size. The equation describes an overall exponentially growing population and is therefore linear. It can thus be solved as an eigenvalue problem that has one positive eigenvalue corresponding to exponential growth of the total population while the distribution of the number of fast and slow codons in the population is stationary. (ii) Here we consider a deterministic evolution equation, thus neglecting genetic drift. For the case C_fast = C_slow = 1, the stochastic effects resulting from a finite population size N_pop have been studied previously. The finite population size was found to effectively modify the mutation rate ν₀ to ν₀+1/(4N_pop) [53]. (iii) The main parameter of the model is , where We provide an estimate of for E. coli below.

Solution of the Evolution Equation

Within our model of codon usage evolution, the fraction of slow codons (or the probability that a given codon is a slow one) in a protein emerges as the natural measure for the usage of synonymous codons. This quantity is calculated from the average number of slow codons via . As codons are independent in our model, can be calculated by considering a single codon, which leads to(7)

We note that for C_fast = C_slow = 1, the exact solution has been given previously [30], [53]. In the limit of weak selection (small or low protein abundance), , which corresponds to random usage of synonymous codons. For strong selection, . Throughout our analysis we use the simpler expression given by Eq. (1) that interpolates between the two limits (plotted in Fig. S4 together with the exact result). It depends on two parameters, and , that describe the codon usage under weak selection and the threshold for protein abundance-dependent selection, respectively. They can be related to the microscopic parameters of our model via the two limiting cases of weak and strong selection, from which we obtain and (where we have used the relation ).

These relations are used to estimate expected values for both N₀ and r₀: In principle, there is one value of r₀ for every amino acid, as r₀ depends on C_slow and C, the numbers of slow and all codons encoding that amino acid (Fig. 3C and E and Table S1). However, we can get an overall estimate for r₀ by averaging these values, weighted with the frequency of the amino acids in all protein-coding sequences in the E. coli genome. This leads to r₀≈0.53. Likewise we can estimate an overall threshold value for selection by ignoring the prefactor C/(C-1), This leads to values of N₀ that are about 1/10 of the abundance of ribosomes (using ). With a ribosome abundance of several 10 0000 molecules per cell [26], this corresponds to an N₀ value of several thousand molecules per cell. This estimate is, however, strongly dependent on the estimate of γ for which we can only estimate a rather wide range (, see below).

Estimate of the Model Parameters

The main parameter of our evolution model is , the ratio of the selection coefficient and the mutation rate, which is proportional to the abundance of the protein under consideration. The dependence of this parameter on the protein abundance is characterized by ν₀ can be estimated by the mutation rate per base which is about 10⁻⁹ per generation [54], under the approximation that all third-position mutations are synonymous. A naïve estimate of the selection parameter, known to overestimate the selection strength [30], is obtained using the doubling time τ₀∼1 hour and taking the translation time difference between synonymous codons to be 10% of the average translation time per codon (∼0.1 seconds [12]). This leads to γ ∼ 10³. This value strongly overestimates the strength of selection mainly because in natural environments, fast exponential growth is not typical. However, selection for fast translation is likely present only during rapid cell doubling. We can therefore improve the estimate by multiplying it with the fraction of total time that the cells actually grow fast. A rough estimate is a few hours per a few days or more [55], [56], so γ would be reduced by a factor of 0.001–0.1. We thus expect to find values of γ in the range 1–100. These estimates of γ also imply that our assumption of is fulfilled.

Alternatively γ can be estimated from codon usage data using the results of our model: One option is to use the sequences of ribosomal proteins () and to determine γ from their average fraction of slow codons. This procedure requires no proteome data and leads to γ≈10. A second option is to obtain γ from our fit of the model to the dependence of codon usage on protein abundance (see Results) via . From that analysis we obtained γ≈1.4 using the abundance data of Lu et al. [37] and γ≈100 using the abundance data of Ishihama et al. [36], both in the range obtained from the microscopic estimate given above. The difference between the two values is mostly due to the different abundance of ribosomal proteins in the two data sets.

Co-evolution of Codon Usage and Codon Preferences (or tRNA Concentration)

To test the hypothesis that evolution of codon preferences or tRNA concentration occurs on larger time scales than the evolution of the usage of individual codons, we obtained tRNA gene copy numbers from the database tRNADB-CE [57] for 12 enteric bacteria and determined the correlations between. Even though tRNA gene copy numbers are not a direct measure of tRNA concentrations or codon preferences, they are correlated [17], [58]. Nevertheless, there are some difficulties associated with the use of tRNA gene copy numbers. In particular, we do not find a systematic variation of the correlation coefficients plotted in Figure 1B with the phylogenetic distance or the frequency of mutations. The latter was estimated by aligning rpoB sequences [a commonly used protein-coding alternative to ribosomal RNA sequences for phylogeny and ecological analysis [59]] using ClustalW [60], [61] and counting synonymous and non-synonymous mutations. Therefore, we also counted changes in the ‘preferred codons’ for each amino acid as defined by Hershberg and Petrov [34], a measure of codon preference that requires only genomic sequence information.

Identification of the Fast and Slow Codons

As mentioned above, the time difference Δt could in principle be measured for each synonymous substitution. However, either absolute or relative speeds of translation have only been measured for a small set of codons [12], [13]. For that reason, we adopt a simplified description in which all codons are classified as either fast or slow according to their usage in the gene of ribosomal proteins, which are known to exhibit strong codon bias [4], [18], [19]. Codons that are preferentially used in ribosomal protein genes are considered as fast codons (Table S1). Compared to other criteria [34], [44], [62] such as classifications according to tRNA concentrations, this criterion has the advantage that it depends only on sequence information and does not require the knowledge of tRNA concentrations or translation speeds. It is therefore readily applicable to a species once its genome is sequenced.

Alternative classifications are based on tRNA concentrations (or tRNA gene copy numbers). There are, however, several difficulties in using tRNA concentrations as a measure of the translation speed of a codon: First, two codons recognized by the same tRNA may be recognized with different affinities. This possibility is indicated by the observation that for 7 out of the 9 amino acids encoded by 2 codons (Asn, Asp, Cys, Glu, His, Lys, Phe), there is only one tRNA that reads both codons, yet there is still a clear “preferred” one between the two. In all of these cases, the preferred codon is the one with perfect codon-anticodon matching (Table S1). Second, it has been shown in vitro that kinetic parameters of the ribosome are codon-dependent even when corrected for tRNA concentrations, although tRNA concentration is likely to account for the biggest effect [15]. Furthermore, if a codon is read by several tRNAs, their concentrations would have to be weighted with the corresponding affinities. Nevertheless, tRNA concentrations are known to correlate with codon preferences [3], [17] and the choice of fast codons according to our criterion of use in ribosomal genes mostly agrees with the choices obtained based on other criteria (Table S1). In particular, we emphasize that whenever we identify a single preferred codon, it agrees with the ‘optimal codon’ in Ikemura’s codon hierarchy [62], which is based on tRNA concentrations together with structural arguments, and with the bioinformatics-based ‘preferred codon’ of Hershberg and Petrov [34]. When we classify several codons encoding the same amino acids as ‘fast’, the ‘preferred’ or ‘optimal’ codon according to these two approaches is always found among our ‘fast codons’, although it is not always the one that is used most frequently in ribosomal protein genes (Table S1).

Average Fraction of Slow Codons and Error Estimates

The relation between protein abundance (N_p) and the fraction of slow codons (r_p) provided by our evolution model applies also to individual codons if one interprets r_p as the probability that a certain codon is a slow codon. Therefore the average fraction of slow codons (e.g., the red points in Figure 2) in a group of proteins of similar abundance can be determined in two ways: We can take the point of view of the proteins and calculate the (weighted) average r_p value over these proteins, with the weights given by the sequence length L_p. Alternatively, we can calculate the fraction of slow codons in the set of all codons of these sequences as an estimate of the probability that such a codon is a slow codon. In that case, we treat the sequences of these proteins as if they were one long sequence and discard all information about the identity of proteins. Because we weight sequences by their length, the two calculations lead to the same result for the average fraction of slow codons. However, they result in different estimates of the accuracy of that average (Table S3). In the first case, we calculate a (weighted) mean of a set of r_p values that exhibit a certain distribution. The standard error of that mean () reflects the standard deviation of that distribution and is given by(8)

This error estimate is indicated by the red error bars in Figures 2, 3 and S1. From the second point of view, we determine a fraction of codons that are slow among a set of codons. The standard error of that fraction is obtained as(9)where is the total number of codons from all sequences. This error estimate is smaller than the first one (Table S2) and is comparable to the size of the red circles in Figure 2. We interpret the difference between these two error estimates as reflecting the fact that the deviations of the fractions of slow codons in the sequences of individual proteins from the average is not due to random sampling, but due to other roles of codon usage, as also suggested by the lack of a length dependence (Figure 4C and D).