Differential Selection for Translation Efficiency Shapes Translation Machineries in Bacterial Species

1. Introduction

“The dream of a bacterium is to become two bacteria” [1], but the realization of this dream is often delayed by the rate of biosynthesis, especially the rate of translation [2,3,4,5,6,7,8,9], partly because most of the cellular dry weight is contributed by proteins. In well-studied enterobacteria, proteins account for more than half of the cell dry mass [10,11]. In Escherichia coli, the ratio of dry weight to wet weight is 0.2294 [12], and the ratio of protein weight to wet weight is about 0.2 [13]. Thus, proteins contribute about 87% (=0.2/0.2294) to the cellular dry biomass. It is not surprising that the growth rate of E. coli increases with the rate of protein production [11]. One therefore would expect a strong selection on the efficiency of the translation machinery and mRNA features to increase translation efficiency.

1.1. Empirical Studies on Selection Optimizing Translation Efficiency

Efficient translation requires efficient initiation, elongation, and termination, as well as mRNA stability, and limits the rate of biosynthesis in both bacteria [3,6,7,14] and phage [5,15,16,17]. Translation initiation is often the limiting step [2,4,18,19,20]. Efficient initiation in bacteria generally requires (1) AUG as a start codon, (2) a well-positioned base-pairing between Shine-Dalgarno (SD) sequence and the anti-SD (aSD) of the free 3’ end of the small ribosomal rRNA [21,22,23,24,25,26], and (3) no strong secondary structure that would embed the start codon or SD and consequently obscure them from being decoded by initiation tRNA and aSD, respectively [27].

With efficient translation initiation, translation elongation becomes rate-limiting [7,14]. Codon-anticodon adaptation is invariably observed in rapidly replicating organisms [2,6,28,29], especially in highly expressed genes [6,7,28,29,30,31,32,33,34]. The same pattern was also found in phages [5,9,16,17,35]. Experimental replacement of minor codons by major codons or vice versa typically leads to an increased or decreased translation rate in bacteria [7,36,37,38,39,40] and viruses [41,42]. Eukaryotic viruses such as HIV-1 tend to have the codon usage of their early and late genes adapted to their respective tRNA pools [43]. These results have led to the explicit formulation of codon-anticodon coevolution and adaptation theory e.g., [6,44,45,46,47] and codon adaptation indices [7,48,49,50]. Although the relationship between codon usage optimization and protein production was briefly challenged [3], the relationship was fully reestablished by more detailed data analysis [7,14].

Translation termination in bacterial species is mediated by one or two release factors RF1 (decoding UAA and UAG) and RF2 (decoding UAA and UGA). All three stop codons can be misread by tRNAs in bacterial species [51,52]. The readthrough frequency is at least 10^-3 to 10^-2 for UGA in Salmonella typhimurium [53] and E. coli [54,55], and 1.1×10^-4 to 7×10^-3 for UAG, depending on the nature of the downstream nucleotides [52,56,57,58]. The readthrough rate is the lowest for UAA, at frequencies from 9 ×10-4 to less than 1×10-5 [52]. The readthrough error rate in the order of UGA > UAG > UAA has been consistently observed in multiple studies [55,59,60,61,62,63,64]. Empirical evidence supports the hypothesis that highly expressed genes prefer UAA stop codons [65].

Termination efficiency also depends on the nucleotide immediately downstream of the stop codon [66,67,68], leading to the proposal of the tetranucleotide stop signal including the +4 site [66,68,69,70,71,72]. The best-documented case involves the translation of prfB mRNA (encoding RF2) in E. coli in which an inframe UGA stop codon is followed by nucleotide C [73,74,75]. When RF2 is abundant, the inframe UGA is decoded correctly to terminate translation, generating a short non-functional peptide. When RF2 is rare, the inframe UGA is not decoded. A +1 frameshift leads to the translation of GAC at a different coding frame, generating a functional RF2.

Almost all studies reviewed above focus on the optimization of mRNA but not on the translation machinery which consists mainly of ribosomes, translation initiation factors, tRNAs, tRNA-charging enzymes, and release factors, as well as the energy that drives the translation machinery [9,76], pp. 522-635. Given that ribosomes represent parallel translation machines, there are two lines of empirical evidence suggesting that short-generation bacteria (SGB) should produce more ribosomal proteins, rRNA [77,78], and tRNA [79] than LGB. First, bacterial growth rate increases with the number of ribosomes in the cell [80,81,82,83]. Second, the number of ribosomes increases with rRNA abundance [84]. More ribosomes imply a greater need for tRNAs [85]. There should be more rRNA operons and tRNA genes in SGB than in LGB.

1.2. Two Approaches in Studying Selection on Translation Optimization

The first approach to identify the effect of selection on translation optimization is by within-species comparisons, especially by contracting between highly expressed and lowly expressed genes (HEGs and LEGs) in rapid replicating bacteria [5,8,9,31,47,65]. Because AUG and UAA represent stronger start and stop signals, respectively, than their alternatives, one would predict that HEGs should exhibit a stronger preference for AUG and UAA as start and stop codons than LEGs [5,8,9,65]. Because a strong secondary structure embedding translation signals such as start codon, SD sequence, and stop codon, one would predict that HEGs should have a weaker secondary structure in sequence flanking these translation signals than LEGs [15,86]. Because early genes and late genes are translated with dramatically different tRNA pools, the two groups of genes should exhibit different codon optimizations [43]. The studies reviewed above constitute empirical tests of such predictions involving different groups of genes within individual species.

This paper focuses on differential selection for translation optimization among species. Relatively few studies are available, partly because of the difficulty in identifying differential selection. Occasionally, codon usage itself was taken as a proxy of selection [77], which is not a satisfactory approach. However, in at least four cases, such differential selection can be identified with little controversy. The first case involves mitochondrial tRNA and protein-coding genes. Mitochondrial genomes of most multicellular eukaryotes encode only a single tRNA^Met/CAU (where Met is amino acid methionine and CAU is the tRNA anticodon for Watson-Crick base-pairing with Met codon AUG) to decode both Met codons AUG and AUA (through wobble-pairing). In contrast, mitochondrial genomes in some bivalve species and tunicate species encode an additional tRNA^Met/UAU gene which can decode AUA codons without wobble-pairing. One can predict that AUA codons should be used more frequently in mitochondrial protein-coding genes in those bivalve and tunicate species, which is true [29,87,88].

The second case of identifiable differential selection involves stop codons and release factors (RF1 decoding UAA and UAG, and RF2 decoding UAA and UGA) in bacteria [65]. If bacterial species X has more RF2 but less RF1 than species Y, then we would predict more UGA usage in species X than in species Y. This predicted pattern has also been documented in multiple bacterial species [65].

The third case involves phages differing in the presence of a lysogenic phase. Phage genomes in a lysogenic phase generally are not under selection for translation efficiency because the phage genome replicates by hitchhiking on the host genome. In contrast, phages without a lysogenic phase will be more selected for translation optimization. This leads to the prediction that the former should exhibit stronger translation optimization than the latter, which is consistent with the empirical evidence [5,16].

The fourth case involves a key component of the translation machinery, i.e., small subunit rRNA whose functions depend on stable secondary structures. Thermophiles are expected to have longer and more GC-rich stems than mesophiles, and this prediction is empirically supported [89]. The difference in growth temperature between thermophiles and mesophiles also affects the proportion of SD-led genes [90].

1.3. Hypothesized Differential Selection on Translation among Bacterial Species

Bacterial species differ dramatically in generation time. Under optimal culture conditions, the generation time is ~10 min in Vibrio natriegens [91,92], 16-20 min in Vibrio cholerae [93], 20-30 minutes in Escherichia coli [94], 30-70 minutes Bacillus subtilis [95], 103-107 minutes in Haemophilus influenzae [96], about two hours in Mycobacterium smegmatis [97], 4-5 hours in M. abscessus [98], 20-30 hours in M. tuberculosis [99,100,101] and more than seven days in Mycobacterium leprae [100,102]. Such dramatic differences in generation time have been dichotomized into r- and K-selection [103,104]. However, the evolutionary consequence of such differential selection has rarely been explored at the molecular level.

Translation machinery that needs to complete the task of protein synthesis in 10 minutes, as in the case of V. natriegens, should be under stronger selection for optimizing translation than that in seven days, as in the case of M. leprae. The fitness of some parasitic bacteria is more related to survival against host attacks than to efficient translation. For example, M. tuberculosis cells build their private niches essentially isolated from the surrounding environment. This protective environment, while beneficial for the survival of the pathogen, also makes it difficult for the pathogen to gain nutrients and oxygen. If protein production is not the limiting factor in growth and reproduction, then one would not expect strong selection optimizing the translation machinery. Thus, we predict that selection for optimizing translation should leave much stronger signatures in rapidly replicating species such as V. natriegens than in slowly replicating species such as M. tuberculosis and M. leprae.

2. Materials and Methods

We include nine bacterial species with completely sequenced genomes and a well-documented generation time from the shortest to the longest among mesophiles under optimal growth conditions (Table 1). Generation time has typically been quantified only in model organisms (e.g., E. coli as a model species for gram-negative bacteria and Bacillus subtilis as a model species for gram-positive bacteria) and very harmful pathogens, so the sample of nine species is not representative of all bacteria. Also, the optimal growth condition from an experimenter’s perspective may not be the true optimal, so the generation time in Table 1 might be overestimated. However, the rank of the relative length of generation time (Table 1) should be reasonable, and subsequent analysis will be mainly based on the rank of generation time. Our general prediction is that SGB should have more efficient protein-production systems than long-generation bacteria (LGB).

The RefSeq genomic sequences for the nine species were downloaded from GenBank by using the accession numbers in Table 1 Software DAMBE [105] was used to extract coding sequences (CDSs), rRNA, and tRNA genes from the GenBank files. We classified the bacterial host genes into two expression groups: 1) known highly expressed protein-coding genes (HEGs) including small and large ribosomal protein genes, RNA polymerases, and some other genes known to be highly expressed (Supplemental Table S1), and 2) the rest of protein-coding genes not in the HEG group (REST). Pseudogenes were excluded from all analyses.

RNA secondary structure could embed important translation signals such as start codon, SD sequence, and stop codon and prevent them from being decoded by the translation machinery. We measured secondary structure stability by the minimum folding energy (MFE) implemented in DAMBE [105] which uses the Vienna RNA fold library [106] for secondary structure characterization. A sliding window of 40 nt along the CDSs was used to characterize the change of MFE along the sequence. Codon adaptation was measured by the index of translation efficiency (I_TE) which generalizes the conventional codon adaptation index (CAI) [48,49] to accommodate background mutation rate [7].

The base-pairing between the SD sequence on mRNA and the aSD sequence on the small subunit (ssu) rRNA was hypothesized to position the start codon at the P-site to pair with the initiation tRNA [5,9,15,107]. The distance between the end of the ssu rRNA and the start codon (D_toStart) is highly constrained in bacterial genes, especially highly expressed ones. A narrow distribution of D_toStart values suggests a stronger selection on the SD/aSD pairing than a wide distribution. DAMBE [86] implements the calculation of D_toStart

3. Results

3.1. Differential Investment in Translation Machinery

Ribosomes represent parallel protein-production factories and their abundance in E. coli and their abundance increases linearly with growth rate [108,109]. Short-generation species may have not only more efficient factories but also more factories than LGB. Each ribosome features a set of 16S, 23S, and 5S rRNAs that are transcribed from the same operon and processed into individual rRNAs [110]. E. coli has seven rRNA operons (rrnA to rrnE, rrnG, rrnH) with promotors that are almost identical to the -10 and -35 consensus [111,112], suggesting a high demand for rRNA molecules met by both efficient and parallel transcription of multiple rRNA operons. The production of ribosomes in E. coli is limited by rRNA production [113], which explains why E. coli maintains multiple rrn operons in its genome for parallel transcription. A generalization of this would lead to the prediction that SGB should have more rrn operons than LGB. This prediction should also apply to tRNA genes because, with more ribosomes, more tRNA molecules are needed [114].

The two predictions are both supported by empirical evidence (Table 2), with the number of rRNA and tRNA genes decreasing highly significantly with increasing generation time (p < 0.0001 for any rank-based nonparametric tests). The observation that SGB maintains more rRNA and tRNA genes in their genomes than LGB is consistent with the interpretation of stronger selection for translation efficiency in SGB than in LGB.

Note that an organism needs at least one rrn operon for translation, so the number of rrn operons cannot be smaller than 1. Also, there should be a minimum set of tRNA genes to decode all 61 sense codons. One might use the following two equations to model the number of rrn operons and tRNA genes, respectively:

$$N_{rrn}=1+a{e}^{-bx}$$

(1)

$$N_{tRNA}=c+a{e}^{-bx}$$

(2)

where x is RankGT in Table 2 Eq. (1) ensures a minimum N_rrn of 1, and Eq. (2) ensures a minimum N_tRNA of c. The parameters were estimated by the least-squares approach, and the fitted equations were plotted in Figure 1.

3.2. Differential Preference of Start Codon AUG

Among the three canonical start codons (AUG, GUG, and UUG), translation initiation efficiency is consistently in the order of AUG > GUG > UUG [115,116]. HEGs in bacteria and bacteriophages tend to prefer AUG as a start codon [5,9,18,115,116,117]. This is true for all nine bacterial species in which the percentage of AUG (AUG%, Table 3) is consistently higher in HEGs than those in REST (Table 3) which includes all genes not in the HEGs group and consequently contains both highly expressed and lowly expressed genes. This suggests that AUG is the most efficient start codon and is consistent with previous empirical studies based on within-species comparisons [5,8,9].

Given that AUG is the most efficient start codon, one would predict that SGB should use more AUG as start codons than LGB. AUG% in Table 3 is indeed strongly associated with ranked generation time in Table 1 as predicted. However, AUG% is also affected by genomic GC% because GC-rich genomes tend to have more GC-rich codons [29]. In terms of start codon usage, GC-rich genomes may use more GUG as start codons than AUG. Because the four LGB all have higher GC% than the five species with relatively short generations, they may use fewer AUG and more GUG as start codons simply because of their genomic GC-richness. For this reason, it is necessary to include genomic GC as a control variable. Also, HEGs and REST genes may differ in the relationship. A conceptually more comprehensive and explicit model is therefore needed.

Here the dependent variable is AUG%, and independent variables include ranked generation time (RankGT), genomic GC content (GC%), and gene expression (GE with two categories, HEGs and REST encoded as 0 and 1, respectively). The input data (Table 4) is used to fit the model.

A regression analysis of the input data in Table 4 shows that all three independent variables are statistically significant, but interaction terms are not (Table 5). The regression model was then fitted without interaction terms. This reduced model accounts for 82.0% of the variation in the dependent variable AUG%. The two regression equations derived from the regression coefficients in Table 5, one for HEGs and the other for REST, are

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{H}\mathrm{E}\mathrm{G}\mathrm{s}:AUG\%=125.86-2.54RankGT-0.52GC\%$$

(3)

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{R}\mathrm{E}\mathrm{S}\mathrm{T}:AUG\%=116.68-2.54RankGT-0.52GC\%$$

(4)

Eqs. (3) and (4) and Table 5 show that AUG% decreases highly significantly with generation time, and that AUG% is significantly higher in HEGs than REST genes (Table 5). This is consistent with the interpretation of stronger selection operating on SGB than on LGB. That is, a non-AUG start codon mutating to AUG is more strongly favored by natural selection in SGB than in LGB.

We should emphasize that Eqs. (3) and (4) are descriptive models. They do not explicitly prevent AUG% from taking values smaller than 0 or larger than 1. A sigmoidal function would have been more appropriate if there were enough data for parameter estimation.

3.3. Differential Preference of Stop Codon UAA

As reviewed previously, stop codon UAA exhibits the smallest readthrough error rate among the three nonsense codons [55,59,60,61,62,63,64]. Consequently, HEGs favor stop codon UAA over other stop codons in multiple bacterial species [9,65,72]. This is true for all nine bacterial species (data not shown). Given that UAA is the best stop codon, we predict that SGB should use more UAA as stop codons than LGB. Similar to our analysis of the start codon AUG usage, the dependent variable now is UAA%, and independent variables include ranked generation time (RankGT), genomic GC content (GC%), and gene expression (GE with two categories, HEGs and REST encoded as 0 and 1, respectively). The genomic GC% is particularly relevant in studying UAA usage because protein-coding genes in a GC-rich genome tend to use more UGA and UAG stop codons than those in an AT-rich genome [65].

A regression analysis shows that all three independent variables are statistically significant (Table 6). The model accounts for 94.7% of the total variation in UAA%. The two regression equations derived from the regression coefficients in Table 6, one for HEGs and the other for REST, are

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{H}\mathrm{E}\mathrm{G}\mathrm{s}:UAA\%=197.95-3.78RankGT-2.22GC\%$$

(5)

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{R}\mathrm{E}\mathrm{S}\mathrm{T}:UAA\%=180.60-3.78RankGT-2.22GC\%$$

(6)

Eqs. (5) and (6) and Table 6 show that UAA% decreases highly significantly with generation time, and that HEGs use UAA significantly more frequently than REST genes (Table 6, a 17.35326% difference between the two). This is consistent with the interpretation of stronger selection operating on SGB than on LGB. A non-UAA stop codon mutating to UAA is more strongly favored by natural selection in SGB than in LGB.

UAA usage also decreases significantly with genomic GC% (Table 6), which is understandable. As genomic GC% increases, GC-biased mutations would favor UAG and UGA codons against UAA codons [65]. It is for this reason that the genomic GC content needs to be taken into consideration when assessing codon usage bias.

3.4. Differential Selection on Sense Codons

It is difficult to evaluate tRNA-mediated selection on codon usage across species because the codon adaptation index (CAI) [48,49] and the index of translation efficiency (I_TE) [7] are both for comparing genes within-species. The effective number of codons (ENC) [50,118] could potentially be used for among-species comparisons, but all these indices are strongly affected by genomic mutation bias [119,120]. For example, protein-coding genes in strongly GC-biased genomes will have mostly G-ending and C-ending codons, leading to a reduced ENC that has little to do with selection.

It is reasonable to assume that mutation bias affects both HEGs and REST genes. Thus, the difference in mean I_TE between HEGs and REST, i.e.,

$$D_{I_{TE}}=({\overline{I}}_{TE.HEG}-{\overline{I}}_{TE.REST})$$

(7)

, should be relatively independent of mutation bias. $D_{I_{TE}}$ should increase with selection intensity for translation elongation efficiency. If there is no selection for translation elongation efficiency, then $D_{I_{TE}}$ is expected to be 0. With increasing strong selection for translation elongation efficiency, $D_{I_{TE}}$ should also increase because the selection is expected to be stronger for HEGs than for REST genes. In other words, strong selection for translation efficiency in SGB should drive HEGs towards better codon adaptation than the REST genes, thereby increasing $D_{I_{TE}}$.

The $D_{I_{TE}}$ values and its ranks (Rank $D_{I_{TE}}$, Table 7 ) are depends strongly on generation time (RankGT, Table 7). The relationship is best illustrated with two ranked variables, i.e., ranked generation time (RandGT) and ranked $D_{I_{TE}}$ (Figure 2). The fitted regression line accounts for 87.84% of the total variation in ranked $D_{I_{TE}}$ (Figure 2).

3.5. Differential Selection Drives tRNA Adaptation

Codon and tRNA are expected to coevolve and adapt to each other [2,6,28,29], especially in highly expressed genes [6,7,28,29,30,31,32,33,34,121]. Given the better codon optimization in SGB than in LGB (Figure 2 and Table 7), one would predict more tRNA genes for highly used codons than rarely used codons. If we focus on the anticodon of tRNA genes, then the prediction above implies a smaller effective number of anticodons (N_AC), equivalent to the concept of the effective number of codons [50,118], in SGB than in LGB. Specifically, N_AC should increase with RankGT (ranked generation time).

We calculated N_AC in the same way the effective number of codons is calculated [50] except that codons in coding sequences are replaced by anticodons in tRNA genes. The ranked N_AC increases highly significantly with ranked RankGT (Table 8 and Figure 3A), consistent with our prediction that SGB should have a smaller N_AC than LGB. One should note the difference between a codon replacement and an anticodon replacement. A codon replacement may have only a minor effect on the translation of a single gene, but an anticodon replacement would affect the translation of numerous codons. For this reason, anticodons are strongly constrained and much less affected by genomic GC%.

N_AC also exhibits a statistically highly significant (p =0.00015, Figure 3B) positive correlation with Rank $D_{I_{TE}}$ (Table 7) which measures codon adaptation. This is consistent with our prediction. The major codons are strongly favored, more tRNAs with the corresponding decoding anticodon will be needed to translate such overused codons, leading to decreased N_AC.

The analysis above assumes that tRNA gene copy numbers in the bacterial genomes are proportional to the abundance of tRNA molecules in the cell. With the availability of transcriptomic data, it has been found that the assumption is generally true, i.e., the copy number of a tRNA is highly correlated with the transcriptomic representation of the tRNA [122].

3.6. Secondary Structure Stability Near the Start and Stop Codons

Sequences immediately flanking translation signals in bacteria (e.g., Shine-Dalgarno sequence, and start and stop codons) tend to have reduced secondary structure to avoid embedding translation signals in a stable secondary structure [18,27,107,116,123,124]. This pattern has also been observed in bacteriophage genes [5,9,15]. The weakening of the secondary structure near, or immediately upstream of, the start codon has also been observed in highly expressed eukaryotic mRNAs [8,125], especially in mRNAs requiring internal ribosome entry for translation [126].

Secondary structure stability in RNA is typically measured by the minimum folding energy (MFE). MFE equal to 0 means no secondary structure, and a stronger secondary structure corresponds to a more negative MFE value. Experiments involving engineered E. coli genes have shown that translation initiation efficiency depends heavily on MFE of the sequence upstream and including the start codon [3,7,14]. We follow the convention of previous studies [5,9,15] measuring MFE with a sliding window of 40 nt along mRNA sequences to quantify the change of secondary structure stability

Secondary structure stability, as measured by MFE, decreases near the start codon, but the weakest secondary structure is observed slightly upstream of the start codon (Figure 4), corresponding to the SD sequence. This pattern has been observed before in bacteriophage genes and their host genes [5,9,15], and is consistent with the interpretation that a strong secondary structure embedding the SD sequence or the start codon is selected against because it prevents the translation initiation signal (SD and start codon) from being decoded by the aSD sequence and the initiation tRNA, respectively.

One might argue against the interpretation that the reduced secondary structure near the translation initiation is to avoid embedding crucial translation initiation signals such as SD sequences and start codons in a stable secondary structure. The SD sequences are purine-rich and cannot form a secondary structure within them. Thus, the reduced secondary structure near the SD sequence (Figure 4) is a direct consequence of the purine-richness in the SD sequences, with nothing to do with the hypothesized avoidance of secondary structure embedding important translation initiation signals. To test this possibility, we simulated equivalent sequences with a purine-rich core about 10 nt upstream of the start codon, with the rest of the upstream sequences randomly drawn from a nucleotide pool derived from intergenic sequences. The result (not shown) does not support this possibility. The purine-rich core decreases the secondary structure stability only negligibly.

All nine bacterial species exhibited decreased secondary structure (Figure 4), and the minor differences in MFE among species are not straightforward to interpret, partly because MFE is strongly affected by GC content, i.e., stronger secondary structure with higher GC content in mRNA. As we can see from Table 9, H. influenzae has the lowest genomic GC% (38.2%), and its mean MFE is the closest to 0 (which means no secondary structure). In contrast, M. smegmatis and M. tuberculosis have the highest genomic GC%, and their mRNAs tend to have more negative MFE values (Figure 4). For this reason, without controlling for GC%, the observation that the LGB have a stronger secondary structure than the SGB cannot be attributed to a reduced selection in these LGB against stable secondary structure.

Given that a weak secondary structure near the SD sequence and the start codon is favorable (Figure 4), one would predict that SGB should have weaker secondary structures (larger MFE values) than LGB. In order to test this prediction, we characterized the MFE plot of sliding windows for each species by a single value. That is, for each plot in Figure 4, we calculated the mean value for mid-window-sites 46-65 (with the start codon occupying sites 58-60). These sites should include both the SD sequence and the start codon. The mean MEF (Table 9) is now the dependent variable. It is expected to 1) decrease with increasing generation time (RankGT), 2) decrease with increasing GC%, and 3) change with gene expression (GE, i.e., be greater for HEGs than for REST genes). These three independent variables are also included in Table 9.

The best model, which accounts for 98.2% of the total variation in MeanMFE in Table 9, includes the three dependent variables and an interaction term (Table 10). The two-tailed p for RankGT is 0.064 (Table 10). However, because we have an explicit one-tailed prediction of a negative slope (i.e., MeanMFE should decrease with increasing RankGT), p should be half of 0.064, i.e., 0.032. Other regression terms for GC%, GE and their interaction are highly significant (Table 10).

As before, we write the two regression equations separately for HEGs and REST from regression coefficients in Table 10:

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{H}\mathrm{E}\mathrm{G}\mathrm{s}:MeanMFE=0.57408-0.08382RankGT-0.10038GC\%$$

(8)

$$\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{R}\mathrm{E}\mathrm{S}\mathrm{T}:MeanMFE=5.03448-0.08382RankGT-0.21009GC\%$$

(9)

Eqs. (8) and (9) show that, for both HEGs and REST genes, secondary structure stability increases with generation time (MFE becomes more negative with increasing generation time). This is consistent with our prediction that selection against secondary structure near the translation start signals (SD sequence and start codon) is stronger in SGB than LGB.

MeanMFE decreases more sharply with GC% in Eq. (9) than in Eq. (8), i.e., secondary structure stability increases more rapidly with GC% for the REST genes than for HEGs. This is easy to understand if the selection against secondary structure stability is on average stronger in HEGs than in REST genes. MeanMFE decreases by only 0.10038 (Eq. 8) for an increase of GC% by 1% with the strong selection in HEGs, but decreases by 0.21009 (Eq. 9) for the same change in GC% with the relatively weak selection in REST genes.

The selectionist interpretation above does not consider the effect of mutations which offers an alternative interpretation. In general, spontaneous mutations in AT-rich genomes tend to be AT-biased, based on comparison between pseudogenes and their functional counterparts [127,128], on mutation pattern of pathogenic bacteria with relaxed selection [129,130], and on nucleotide bias at the three codon sites across multiple bacterial species [131]. H. influenzae has an AT-rich genome, suggesting AT-biased mutation, in contrast to M. smegmatis which has a GC-rich genome. However, protein-coding genes in both species have SD sequences that are purine-rich (especially G-rich) (Figure 5). The G-rich SD will form base pairs with nearby C nucleotides so MFE will not go to 0. In pseudogenes where selection for maintaining the G-rich SD sequence is absent or in lowly expressed genes where the selection is weak, the AT-rich H. influenzae will lose these G nucleotides in the SD sequence, leading to an MFE closer to 0. Indeed, the MeanMFE value for the 19 pseudogenes in H. influenzae -2.2887, closer to 0 than all MeanMFE values in Table 9. Similarly, the MeanMFE is -3.2660 for HEGs and -2.5466 for REST genes (Table 9). This is consistent with the interpretation that the G-rich SD is more likely to be hit by G→A and G→T mutations and lose G/C base pairs in REST genes than in HEGs. For example, H. influenzae has a GC% of 37.9329% for HEGs but only 33.7986% for REST genes, leading to a MeanMFE value closer to 0 in REST genes than in HEGs.

In short, although the purine rich SD sequences (Figure 5) can hardly form secondary structures within itself, the dramatically increased G nucleotides within the SD sequence could base-pair with the neighboring C nucleotides and contribute to secondary structure stability. If there is no selection maintaining the G-richness in the SD sequences, then these G nucleotides may be replaced by A and T, leading to a further decrease in secondary structure stability. Thus, both selection and mutation could contribute to secondary structure stability in sequences near the translation initiation signals (the SD sequence and the start codon). The models in Eqs. (8) and (9) are therefore oversimplified and should be interpreted with caution. The secondary structure in sequences near the stop codon exhibits a similar pattern as those near the start codon (Figure 6).

The decrease in secondary structure stability may not necessarily be related to the avoidance of embedding SD sequences and start codons. Efficiently translated yeast mRNAs (i.e., mRNAs in polysomes with high ribosome density) often have a short poly(A) tract before the start codon [125], with the poly(A) interpreted as long enough to recruit translation initiation factors but short enough to avoid binding by the poly(A)-binding proteins. However, the presence of poly(A) also weakens the secondary structure stability as a secondary consequence.

4. Discussion

Our assumption that SGB are under stronger selection for translation efficiency than LGB appears to be valid because multiple predictions based on the assumption are consistent with the empirical evidence. First, the number of ribosome RNA operons, as well as the number of tRNA genes, increases as decreasing generation time (Table 2 and Figure 1). Second, AUG was known to be the most efficient start codon, and SGB genes exhibit a stronger preference for AUG as a start codon than LGB, especially in highly expressed genes (HEGs) (Table 3-Table 5). This is also true in the usage of stop codon UAA (Table 6) which is known to be the most efficient termination signal with the smallest readthrough error rate as a stop codon. Third, SGB, especially their HEGs, exhibits better codon and anticodon adaptation than LGB (Table 7 and Figure 2 for codons, and Table 8 and Figure 3). Finally, SGB genes have weaker secondary structure near translation initiation signals than LGB (Figure 4 and Tables Table 9-Table 10). A similar pattern was observed with sequence secondary structure near the stop codon (Figure 6). However, as we discussed previously, the selectionist interpretation is sometimes confounded by biased mutations.

Modern biological research aims to formulate and validate quantitative and mechanistic models. In this context, this study has several shortcomings. The first shortcoming is the inherent inaccuracy of generation time. The experimentally measured generation time sometimes varies widely among strains and among different studies. For example, the generation time is ~21 hours in drug-sensitive strains of M. tuberculosis, but ~35 hours in the multi-drug-resistant strains [132]. This suggests a cost to the pathogen in developing multi-drug resistance, i.e., the resistance is at the cost of longer generations. However, the observation also highlights the inherent variation in measured generation time. Particularly controversial is what generation time to use. For example, the generation time in E. coli is only about 20-30 minutes under favorable culture conditions but could be 15 hours in the natural habitat of the mammalian intestine estimated by the rate of mutation accumulation [133]. For modeling the joint effect of mutation and selection, the generation time under natural conditions would seem more appropriate than that under optimal experimental conditions. Unfortunately, the generation time in bacterial populations in nature cannot be measured accurately. For example, the estimated generation of 15 hours for E. coli [133] is associated with a 95% confidence interval of 0-30 hours. Another uncertainty with generation time is that we do not know if long-generation species such as M. leprae really cannot replicate fast or if they could grow fast but microbiologists have not been able to shift them into the fact-growing mode. Because of the uncertainty generation time among the bacterial species, we ranked generation time in the hope that the relative generation time is in proper order. It is for this reason that we have chosen species widely different generation times characterized experimentally, so that our ranking of the generation time would not be controversial.

The second shortcoming of the study is the small number of bacterial species, partly because of our conscious effort to avoid species with uncertain generation time. There are various tables of bacterial generation times for more than nine bacterial species, but they often do not include original references or do not have fully sequenced genomes. Some bacterial species in Mollicutes have known generation times, e.g., ~6 hours in Mycoplasma pneumoniae [134],[135], as well as fully sequenced genomes. However, their genetic code (i.e., translation table 4) differs from the rest of the bacteria (translation table 11). This complicates comparisons. Consequently, they are not included in this study. The small number of bacterial species included here results in two limitations. First, it does not permit the validation of parameter-rich models. Second, it would not allow a phylogeny-based inference [136,137,138] to alleviate the issue of data dependence. We have previously used such phylogeny-based inference to quantify relationships between body temperature and genome size [139] and between optimal growth temperature and ribosomal RNA secondary structure stability (i.e., the stem length and GC% in the stem-loop structure in rRNA) [89]. In this context, we may highlight two points. First, our results are highly consistent with the prediction that increasing generation time decreases selection on translation efficiency. Second, the bacterial species we use are highly divergent. Even the two most closely related species, M. tuberculosis and M. leprae, have an evolutionary distance of more than 0.3 even for conserved ribosomal protein genes [140], so phylogenetic inertia might not affect the quantified relationships, i.e., the data points might be considered to be roughly independent.

The third shortcoming is the lack of a general conceptual framework of mutation and selection on translation optimization that can be used across species. Protein production from an mRNA depends on the ribosomal recruitment rate, the efficiency in forming the 70S initiation complex, elongation efficiency and accuracy, termination efficiency and accuracy, the stability of the mRNA, differential amino acid and tRNA availability, and the energy level of the cell [4,141,142,143]. How do these variables interact with each other affect protein production? For example, when translation is not efficient, codon usage optimization has little effect on protein production. However, protein production increases significantly with codon optimization in mRNAs with efficient initiation [3,7,14]. Similarly, there are conflicts between maximizing transcription efficiency and translation efficiency. An RNA would be transcribed efficiently if it maximize the usage of the abundant nucleotide A and minimize the rare nucleotide C [144,145]. However, this would drive up the usage of A-ending codons that may not be the optimal codon for translation. There are also factors affecting codon usage that are not related to translation. For example, mammalian species have zinc-finger antiviral proteins (ZAP) against RNA viruses by targeting CpG dinucleotides in the viral RNA [146,147,148,149,150]. Many human RNA viruses, including SARS-CoV-2, exhibit much reduced CpG dinucleotides [151,152]. Most of the CpG reduction occurs at the di-codon configuration, from NNC GNN to NNT GNN (i.e., a synonymous replacement). This NNC to NNU change has nothing to do with codon optimization. What would be the functional relationship among all these variables? The model would be even more complicated if we consider not only HEGs, but also genes whose optimal protein level is not the maximal, such as the autoregulated level of release factor 2 in E. coli [73,74] or many others [153]. We highlight such questions in the hope that they will motivate researchers to address them.

5. Conclusions

Detection of selection on translation efficiency is mostly done by contrasting highly expressed genes and lowly expressed genes within rapidly replicating species. We generalized such studies to understand how differential selection on translation optimization among different species would leave its footprints on the species-specific translation machinery. Our results suggest that selection for translation optimization is stronger in short-generation species than long-generation species, and that this differential selection strongly shapes the evolutionary trajectories of translation machinery in these species, affecting translation initiation, elongation and termination.