A slight mismatch between a gene’s codon usage and the cellular tRNA supply is beneficial

Domain: chemistry
Published: 2026-02-17
Editorial outcome: major_revision

Abstract
Each amino acid except two is encoded by multiple synonymous codons, but at unequal frequencies. Such codon usage bias (CUB) is observable in almost all species, and commonly assumed as the result of natural selection towards an optimal CUB that matches the cellular tRNA supply. Here we hypothesize instead that the optimal CUB of a gene should slightly mismatch the tRNA supply to avoid excessive translational costs, while ensuring adequate functional payoff. By modifying the CUB of a resistance gene expressed in bacteria under antibiotic selection, we demonstrate that a small mismatch with the tRNA supply confers faster bacterial growth than those with minimized or large CUB-tRNA mismatches. Intriguingly, the optimal degree of CUB-tRNA mismatch increases as the resistance gene becomes less important in media with lower antibiotic concentrations, which is explainable by our model as a shift in the balance between the gene's functional payoff and translational cost. Furthermore, genomic analyses in model organisms suggest that the optimal degree of CUB-tRNA mismatch is larger for endogenous genes with lower functional importance and higher mRNA abundance, respectively supporting the impact of functional payoff and translational cost. Finally, we find that mutations increasing or decreasing the CUB-tRNA mismatch of native genes are both predominantly deleterious, such that the CUB-tRNA mismatch is likely selectively maintained rather than minimized to that achievable in the presence of genetic drift and mutational bias. These results challenge the commonly assumed unidirectional selection on CUB and highlight the CUB-modulated balance between functional payoff and translational cost.

Introduction
Eighteen out of twenty amino acids are encoded by two or more synonymous codons. The mutation of one codon to its synonymous codon, or synonymous mutation, is traditionally regarded as functionally neutral, since it does not alter the sequence of a protein. However, there is currently a debate regarding whether and why synonymous mutations may be non-neutral1-3. In this context, some “preferred” codons are used more frequently in the genome than their synonymous counterparts, a phenomenon commonly known as codon usage bias (CUB)4. Other than some special cases affecting only limited neighboring regions5-8, the prevailing hypothesis for the widely observed CUB is the selection-mutation-drift model9,10, which proposes that translational selection favors an optimal CUB, whereas deviations from the optimal CUB persist as a result of genetic drift (or in other words, the drift-barrier11) and mutational bias. For translational selection, it is commonly assumed that the optimal CUB maximally matches the cellular tRNA supply so that the translation is more efficient and/or more accurate10,12-15. In this context, we can quantify the CUB-tRNA mismatch using a previously16 proposed metric, D, which is based on the Euclidean distance between synonymous codons’ demand and supply of tRNA (see Methods). Accordingly, the aforementioned assumption can be expressed as unidirectional translational selection towards the optimal D = 0 (Figure 1a). Although this assumption has not been directly tested, it is often considered supported by the widely observed patterns that highly expressed genes, for which translational selection is stronger than lowly expressed genes within the same genome, tend to exhibit a better CUB match for the cellular tRNA supply10,13.

Nevertheless, as we16,17 and others18-20 have recently showed, the heterologous expression of exogenous genes with CUB overly matching the tRNA supply, and therefore excessively small D, may be detrimental to the fitness of the host cell. One plausible molecular mechanism behind such phenomenon is that the expression of exogenous genes trans-regulates the expression of other genes in the cell via tRNA depletion. An exogenous gene with excessively small D will have a very high protein expression level due to an increased abundance of mRNA21,22 and/or more efficient translation20,23-25, but at the same time, will impose a heavy translational cost on the cell, thus hindering the translation of other genes within the same cell16,20. In theory, the same logic applies to endogenous genes native to the cell, since they also trans-regulates the expression of other genes26. However, endogenous genes are critically different from exogenous genes in that they contribute to the native function or fitness of the cell. Consequently, a D-decreasing mutation in an endogenous gene would both exacerbate its translational cost and enhance its functional payoff for the cell. Whether translational selection on CUB is dominated by the functional payoff and thus unidirectional towards maximally matching the tRNA supply (the “conventional model” in Figure 1a), or a balance between the translational cost and functional payoff so that a slight mismatch with the tRNA supply would be optimal (the “novel model” in Figure 1b), remains to be determined.

To answer this question, we evaluated the functional payoff and translational load of 53 synonymous variants of gentamicin resistance gene GmR exhibited on Escherichia coli MG1655 strain grown in Luria-Bertani (LB) media containing 120 μg/ml gentamicin, where GmR mimics an endogenous gene in terms of contributing to cellular fitness. As compared to variants with moderately small D, those with minimal D have greater functional payoffs, higher translational cost, and a lower overall fitness, suggesting that the optimal CUB does not correspond to a perfect match with tRNA supply. Similar experiments conducted under various gentamicin concentrations (ranging from 0 to 180 μg/ml) revealed that the optimal D increases as the GmR gene becomes functionally less important in media with lower gentamicin concentration, indicating that optimal CUB indeed shifts with the balance between GmR’s functional payoff and translational cost. Additional experiments using Ribo-Seq, tRNA supply modification, or the ampicillin resistance gene AmpR provided consistent support. To further test whether our new model applies to evolution in nature, we analyzed endogenous genes of several model organisms under the assumption that the wildtype codon usage has been evolutionarily optimized to the maximum extent possible in the presence of genetic drift. We found that the optimal D decreases for genes with greater functional importance and increases for genes with higher translational costs, respectively supporting their roles in determining the optimal D of actual endogenous genes. Lastly, we found that a substantial proportion of the point mutations occurred during neutral evolution of the mutation accumulation (MA) lines reduces D, and that D-decreasing and D-increasing mutations in representative genes are equally likely to be deleterious. Both phenomena are compatible with the existence of substantial D-increasing selection in the wild-type genome. Collectively, these results challenged the conventional view that translational selection for codon usage is unidirectional towards a maximum match with cellular tRNA supply, and suggested that a slight mismatch to tRNA supply is beneficial, a previously unrecognized complexity in the evolution of codon usage bias.

Results
A codon usage that maximally matches the cellular tRNA supply is not optimal
We used the gentamicin resistance gene GmR as a model to examine how the codon usage of a gene alters its translational cost and functional payoff, as well as the growth rate of the cell. To this end, a collection of synonymous variants of GmR consisting of three distinct sequences for each of 19 unique D values were designed. Fifty-two of them, together with the wildtype GmR, were successfully synthesized and individually constructed into an expression cassette (see below), covering a wide range of CUB from a near-perfect match with (D = 0.078) to a large mismatch with tRNA supply (D = 0.967) (see Figure 2a and Supplementary Data 1). Variants with maximal mismatch (the theoretical maximal D is 1.414) were not considered since their expression is most likely too low for an observable functional payoff. All GmR variants were connected to a strong promoter J23119, fused with mCherry for the direct readout of its protein expression level, and then inserted to an expression vector containing a YFP driven by a weaker promoter J23114 (Figure 2b, see Methods and Supplementary Table 1). Here the YFP is used as a measure of the overall translational efficiency of the cell and is therefore used to inversely estimate the translational cost imposed by GmR. We individually transformed the 53 vectors into Escherichia coli MG1655 strain, and then cultured these strains in Luria-Bertani (LB) media supplemented with 120 μg/ml gentamicin. This concentration is below the minimum inhibitory concentrations (AMICs) (see Methods and Supplementary Table 2) of 48 out of the 53 strains, thereby ensuring the measurability for the majority of the strains.

The 53 synonymous variants of GmR were divided into three groups of minimum, small and large mismatch between the codon usage and cellular tRNA supply (Figure 2a). It was found that the variants with minimum D showed the highest level of GmR expression as indicated by mCherry fluorescence (Figure 2c), which is consistent with the cis-regulatory effect of CUB on the gene’s expression and functional payoff for the cell. As for CUB’s trans-regulatory effect for other genes, the YFP fluorescence suggests a greater translational cost for GmR variants with smaller D (Figure 2d). These observations collectively suggest that a D-decreasing mutation would elevate both the functional payoff and the translational cost of a gene. Most importantly, the net functional payoff combining the two effects above, which is approximated by a direct measurement of cellular growth rate, appears to be the highest for the variants that have a small, but not minimum D (Figure 2e). These results suggest that, due to the balance between functional payoff and translational cost, a gene’s optimal CUB should exhibit a slight mismatch with the cellular tRNA supply.

The optimal degree of mismatch between codon usage and tRNA supply depends on the gene’s functional importance
Our model also predicts that genes of different functional importance will have different optimal CUBs (Figure 3a). More specifically, for a gene with a given transcription level, its translational cost is determined by its codon usage (Figure 3a middle row), but is unaffected by its functional importance. However, an important gene would have a greater functional payoff (Figure 3a top row), which allows the gene to afford a more severe translational cost and, therefore, a codon usage more closely matches the cellular tRNA supply. Consequently, the optimal CUB of a more important gene would likely correspond to a smaller D, whereas the optimal CUB of a less important gene would likely correspond to a larger D (Figure 3a bottom row). We emphasize that this prediction applies to synonymous variants of any particular gene as long as the gene exhibits a diminishing (relative to the increase of translational cost) increase of functional payoff as D decreases.

To test this prediction, we cultured the 53 strains expressing synonymous GmR variants in LB media supplemented with gentamicin at 12 different concentrations ranging from 0 to 180 μg/ml, so that the GmR gene would have different functional importance. The functional payoff of the GmR gene is estimated by the GmR-mCherry expression, which is rescaled so that the average functional payoff of the group with large D is 1 in each gentamicin concentration to facilitate comparison across conditions. The translational cost exerted by each GmR variant under a given gentamicin concentration is estimated by the fold reduction of YFP expression relative to the maximum YFP expression observed for the focal concentration of gentamicin, which is also further rescaled so that the average translational cost of the group with large D is always 1 in each gentamicin concentration. We found that the GmR variants with large D always exerted the smallest functional payoff among the three groups (Figure 3b and Figure S1), which is compatible with CUB’s cis-regulatory effect. Consistent with GmR’s greater functional importance in media with higher gentamicin concentrations, variants with minimum D displayed a greater functional payoff therein, with 68.9% of the rank variations explained by the gentamicin concentration (Figure 3c). The average translational cost of the GmR variants with large D is also always the lowest among the three groups regardless the concentration of gentamicin (Figure 3d and Figure S1). However, the dependency of translational cost on functional importance is much weaker than that of functional payoff, so that only 27% of the rank variation of translational cost imposed by the GmR variants with minimum D is explained by the gentamicin concentration (Figure 3e). Most importantly, direct measurement of the fitness indicates that, as gentamicin concentration increases from 0 to 180 μg/ml and GmR becomes more important, the optimal CUB gradually shifts from a large D at concentrations of 0 to 60 μg/ml, to a small D at concentrations of 80 to 120 μg/ml, and finally to a minimum D at concentration of 180 μg/ml (Figure 3f and Figure S1). This observation demonstrates that the optimal CUB of a gene depends on the gene’s functional importance, and when the gene is less important, the optimal CUB should be one that mismatches the cellular tRNA supply to a greater extent. We note that even in media with high gentamicin concentrations, GmR synthesis cannot exceed a physical limit (e.g., as determined by cell resource uptake rate). Functional payoffs and fitness of the group with the minimum D may have been affected by this limit, thus making the difference with the group of minimum/small D less statistically significant. Nevertheless, the translational cost should have been less affected, resulting in a stronger pattern (right-most group in Figure 3d).

We conducted four additional analyses to further strengthen our conclusion. First, under the assumption that natural selection has optimized the supply-demand balance of tRNAs for the whole transcriptome (but not necessarily for individual genes), we estimated the effective tRNA supply for each codon based on the combined transcriptomic codon demand of the top 10% highly expressed genes. However, it may be argued that direct measurement of tRNA transcript abundance would provide a better estimate of tRNA supply. Although it cannot capture factors such as tRNA charging, binding of EF-Tu and GTP, and codon-anticodon affinity, we still attempted to estimate tRNA supply based on high-throughput tRNA sequencing data27. As a result, we found minor changes in the value of D, which did not affect the grouping (Figure S2) and therefore our results. Second, to demonstrate that the translation of GmR incurs a translational cost that is dependent on its codon usage and also affects endogenous genes, we used Ribo-Seq to measure the codon decoding time of wildtype and four E. coli strains expressing GmRs of different D (see Methods). For the two strains expressing GmRs of minimum D, decoding time of codons frequently used by GmR tend to increase, especially in endogenous genes with high translational activity (Figure S3). This pattern is absent in the two strains that express GmRs of large D. Combined with the YFP-based data, they indicate that the translation of GmRs of minimum D imposes a genome-wide translational cost. Third, the GmR experiment only manipulated codon demand, while our model also predicts a similar effect for manipulation of tRNA supply. The supply of all tRNAs can be difficult to manipulate, but we attempted to increase the supply of a specific tRNA using a tRNA gene expression cassette (Figure S4a), and examined two strains whose D changed from the other two groups to the small mismatch group (see Methods). As is predicted by the conventional model, a decrease in D caused by altered tRNA supply is always beneficial (or at least not detrimental in low or intermediate gentamicin concentrations). On the contrary, our model predicts that at intermediate gentamicin concentration, both strains should be fitter after the tRNA supply change, whereas fitness changes would be opposite at low or high gentamicin concentration (Figure S4b). The experimental results support our model (Figure S4c/d), although the fitness changes observed for the strain with increased D are not statistically significant, possibly due to the relatively small change in D. Fourth, to further exclude the possibility that the above observation is caused by gentamicin’s interference with the translation process9, we tested again the fitness effect of synonymous codon usage change using another antibiotic resistance genes, namely AmpR, the ampicillin resistance gene (Figure S5a). We constructed 47 synonymous variants of AmpR (see Supplementary Data 2 and Methods), which were divided into three equal-sized groups of minimum, small, and large D. tested whether altered functional importance can result in a shift in optimal D as predicted by our model. In the presence of ampicillin at 100, 2000, and 4000 μg/ml concentration (see Supplementary Table 3), the fittest groups were respectively found to be large, small and minimum D (Figure S5b), which again directly supports our model. Also, we note that these patterns observed for GmR and AmpR are also consistent with our previous observations that expression of viral genes with CUB excessively similar to host CUB is highly deleterious16, because viral genes in theory should have zero functional payoff to host cells. Combined, the manipulative experiments above demonstrated that the optimal CUB for a gene contributing to cellular fitness is not one that perfectly matches tRNA supply, but rather one that exhibits some degree of mismatch with tRNA supply.

The optimal CUBs of native endogenous genes are affected by both functional payoff and translational cost
We then ask whether the conclusion drawn from the manipulative experiments on resistance genes is applicable to native endogenous genes as well. Based on the same logic as above (Figure 3a), our model predicts that a more important gene should have an optimal CUB that is more closely aligned with cellular tRNA supply, and therefore correspond to a smaller D. To gauge the optimal CUB for the native genes, we assumed that the native codon usages of endogenous genes have been optimized by translational selection to the maximum extent possible in the presence of genetic drift. Here, the relative strength of translational selection and genetic drift varies for different genes within the genome. By following previous practice28, we presumed that genetic drift (determined by effective population size) was constant across genes, and that translational selection was determined by mRNA expression levels. Therefore, our hypothesis predicts a negative correlation between a gene's functional importance and the D value of its native codon usage, especially for highly expressed genes because they are under stronger translational selection that keeps their native codon usage close to their optimal CUB.

Using the genomic data of E. coli strain MG1655, we found that the functional importance of a gene, as estimated by the fitness reduction upon its replacement by KanR (kanamycin-resistance)29, is indeed anticorrelated with its native D value (Pearson’s R = -0.09, P < 10-6 ; Spearman’s ρ = -0.05, P < 0.01). The anticorrelation is stronger among genes whose mRNA expression levels are higher, which cannot be explained by the reduced number of genes since random downsampling to the same number of genes will not generate the same result (Figure 4a. See Figure S6a for tRNA expression-based D). The same results were observed using genomic data30,31 of Saccharomyces cerevisiae strain BY4742, for which functional importance of 4679 genes are also available from replacement by KanR (Figure 4b and Figure S6b). These results collectively supported the notion that genes with greater functional payoff can withstand a more severe translational cost and, therefore, have optimal CUB more closely matches the cellular tRNA supply.

While the above results supported the role of functional payoff in determining the optimal CUB, demonstrating the role of translational cost is more important. This is because, under the conventional model, where the optimal CUB matches tRNA supply perfectly, the anticorrelation between gene importance and native D would also be predicted because they are both confounded by gene expression level32. That is, highly expressed genes will likely have high functional importance and also be subjected to greater translational selection towards the optimal D = 0 (as posited by the conventional model), therefore exhibiting smaller D. Accordingly, the key feature of our model that differentiates it from the conventional model is the translational cost shown by CUB corresponding to excessively small D. More specifically, a D-decreasing mutation in a gene with high transcript abundance should result in a greater increase of translational cost than that in a gene with low transcript abundance, since a greater number of RNA molecules are mutated in the former case. Consequently, highly expressed genes should be more sensitive to reductions in D, resulting in a larger optimal D in highly expressed genes than lowly expressed genes of similar importance (Figure 5a). Once again, this phenomenon should be more pronounced among genes with stronger translational selection, as their native CUB should more closely resemble their optimal CUB.

Based on this logic, we tested the correlation between mRNA expression levels and the D value of native codon usage among E. coli endogenous genes whose expression exceeds a certain threshold. We found that D and expression level are positively correlated among the most highly expressed genes (Figure 5b, purple bars. See Figure S7a for tRNA expression-based D), supporting the effect of translational cost due to excessively small D. The positive correlation is not an artifact of the limited number of genes above the expression threshold, since randomly downsampling all genes to the same number as those retained at each expression threshold will not result in the same correlation (Figure 5b, gray bars. See also Figure S7a). Moreover, endogenous genes in yeast, fly, and human show similar patterns (Figure 5c-e and Figure S7b-d), and these patterns cannot be attributed to potential base composition biases (Figure S8) or the dominance of an outlier amino acid (Figure S9). Interestingly, the maximum correlation achieved is the weakest in human compared to the other three species, which is consistent with the weak translational selection observed for metazoan species33. It is evident from these observations that translational cost strongly influences the optimal D when a native gene’s CUB excessively matches the cellular tRNA supply. Together with the results that demonstrate the effect of functional payoff, we have found evidence for both facets of our model regarding an optimal CUB as dictated by a balance between translational cost and functional payoff.

The mismatch between endogenous gene’s codon usage and cellular tRNA supply is selectively maintained
To further verify our model, we set out to test the direction of selection on CUB of native endogenous genes of E. coli, which is predicted in conventional model as a unidirectional selection towards a perfect match with the cellular tRNA supply, or optimal D = 0 (Figure 6a). Our novel model, however, predicts converging selections from both directions, which maintained an optimal CUB with a slight mismatch with the cellular tRNA supply (Figure 6b). In this context, it is possible to determine the type of selection acting on a trait from its change upon the removal of the selection34. Specifically, evolution in the absence of selection should tend to cause a rise in D if D has been selectively minimized to the drift barrier, as posited by the conventional model (Figure 6a). On the contrary, if natural selection has maintained D at a moderately small level, the probabilities that D will increase or decrease upon the removal of selection are both substantial (Figure 6b).

Following this logic, we extracted the single nucleotide synonymous mutations observed in E. coli and S. cerevisiae mutation accumulation (MA) lines. These mutations were observed in wildtype or mismatch repair deficient strains (∆mutL for E. coli and ∆MSH2 for S. cerevisiae) evolving with minimized natural selection34-37. Take E. coli an example, among the top 15% highly expressed genes, whose sequence supposedly has been nearly optimized by translational selection before MA, we found that 49.4% of single nucleotide mutations are D-increasing (Figure 6c, red arrow). This fraction does not differ from its expected value of 49.2±3%, which is determined by sampling the same number and types of mutations on the same genes at random (Figure 6c, gray histogram. See Methods). This observation is not caused by a lack of statistical power, since statistical significance can be achieved when a “positive control” set of mutations constructed by random mutations in genes with minimum D values was subjected to the same test (Figure S10a. see Methods). D-increasing and D-decreasing mutations found in E. coli MA lines have average effect sizes of 0.00257±0.000196 and 0.00298±0.000273 on D, respectively, which again show no difference from their respective random expectations (Figure 6d-e). These observations made with E. coli are all consistently recapitulated in S. cerevisiae (Figure 6f-h and Figure S10b). Combined, these results strongly suggest that D is not selectively minimized in the ancestor of MA lines.

Additionally, a recent study on the fitness effects of point mutations based on genome editing has enabled us to examine synonymous mutations directly in 21 representative yeast genes1 and distinguish between our model and the conventional model. In particular, while both our model and the conventional model predict that D-increasing synonymous mutations are mostly deleterious, only our model predicts that D-decreasing synonymous mutations are also mostly deleterious. As a result of analyzing all synonymous mutations in the dataset, we found that an overwhelmingly large percentage of both D-decreasing and D-increasing mutations are deleterious (Figure 6i). A more quantitative analysis of the mutations’ fitness effects shows that both D-decreasing and D-increasing mutations are deleterious on average, with greater deleterious effects in genes with higher expression levels (Figure 6j and Figure S11), which is consistent with stronger translational selection38. All these results are recapitulated by tRNA expression-based D (Figure S12). Collectively, these observations strongly suggest that the native codon usage of endogenous genes is selectively maintained at a slight mismatch with cellular tRNA supply, thus supporting our hypothesis and rejecting the conventional model.

Methods
Measuring codon usage preference
Codon usage of a gene was measured by the disparity between the codon demand of the focal gene and the tRNA supply of the cell as follows:
𝐷 =√∑𝑛𝑖 (𝑌 −𝑋 )2 (1)
𝑖 𝑗=1 𝑖𝑗 𝑖𝑗
Here i represents one of the 18 amino acids that has two or more synonymous codons; n is the number of synonymous codons for amino acid i; Y represents the relative demand of synonymous codons, as approximated by the fraction of codon j among all synonymous codons of amino acid i found in the open reading frame (ORF) of the focal gene; X represents the tRNA supply of the cell, approximated by the fraction of codon j among the synonymous codons of amino acid i for the highly expressed genes in the transcriptome of the cell 52. It is important to note that we are using the transcriptome, not the genome, which means the codon usage has been weighted by the expression of those genes before being summed up. Using all genes would not alter our result as the highly expressed genes dominated the transcriptomic codon usage. Essentially, D is the Euclidean distance between the codon demand and tRNA supply among synonymous codons encoding the amino acid i. Finally, D equals the geometric average of all 18 𝐷. It should be noted that cellular tRNA supply in this context encompasses factors such as tRNA expression levels, the rate of tRNA aminoacylation, and the formation of the EFTu•GTP•aminoacylated-tRNA tertiary complex, as well as codon-anticodon binding affinity considering wobble base pairs, etc. In other words, D measures how well the translational demand of a codon matches the cell’s effective supply of tRNAs, which might not be accurately reflected by tRNA expression profile, and is a better fit for our codon-centered analysis than a tRNA-centered approach. Nevertheless, we extracted tRNA expression levels of E. coli27, S. cerevisiae20, H. Sapiens and D. melanogaster53 from previous publications and re-calculated D using tRNA expression-based X , and obtained qualitatively unchanged results (Figure S2, S6, S7 and S12).

Construction of E. coli strains
A set of 53 synonymous mutants were designed based on the reference sequence of GmR (encoding gentamicin-3-acetyltransferase), with similar guanine and cytosine content (GC%, 46% ~ 55%) but 19 different levels of codon usage bias (D ranges from 0.078 to 0.967. Figure 2a and Supplementary Data 1). All synonymous mutations among these 53 mutants are located in the nucleotides 60 to 531 of the open reading frame (ORF). This design facilitates easier construction of expression cassettes for the mutants (see below), and avoids modifying the first 57 nucleotides of the ORF, which may have an impact on protein expression in a tRNA-independent manner 54,55.

A template reporter cassette (Figure 2b and Supplementary Table 1) consisting of a GmR fused with mCherry in the forward strand and a YFP in the reverse strand was synthesized and ligated to the pUC57 vector (Supplementary Table 1). The other synonymous GmR mutants were directly synthesized from the 37th nucleotide in the ORF to the 22nd nucleotide downstream of its stop codon (the stop codon itself is excluded, giving rise to a total length of 517 nucleotides), and amplified using primers GmRm-F and GmRn-R listed in Supplementary Table 4. Using the ClonExpress II One Step Cloning Kit (Vazyme, no. C112) according to the manufacturer's instructions, these segments were ligated to other elements on the template reporter cassette (Figure 2b and Supplementary Table 1), which were amplified from the template using primers PUC-F and PUC-R in Supplementary Table 4. These elements include (i) a strong constitutive promoter (J23119) that drives the expression of GmR-mCherry fusion protein, ensuring that GRmR has a sufficient impact on cellular translation; (ii) a YFP gene serving as a readout for the overall translational efficiency of the cell, which is driven by a relatively weak promoter (J23114)56 so that its impact on host tRNA supply is negligible. The expression cassette also include two efficient and strong ribosome binding sites (RBS)57,58, as well as the bidirectional E. coli tonB-P14 transcription terminator59. The 53 constructed vectors were individually transformed into E. coli MG1655 cells by the typical calcium chloride method60, and selected on LB with 50 μg/mL ampicillin. The MICs of gentamicin for the 53 constructed strains were determined using the agar dilution method and interpreted using breakpoints defined by the Clinical and Laboratory Standards Institute61.

We aimed to overexpress one tRNA by inserting it into the high copy number PUC57 vector, and we assumed this would increase the supply of this specific tRNA at least twice as much as the supply of the most abundance synonymous tRNAs in the wild-type strain (I.e., set the 𝑋 of the affected codon 𝑗 to twice the maximum 𝑋 of synonymous codons for amino acid 𝑖, then rescale all 𝑋 of amino acid 𝑖). We selected the argWCCT tRNA gene and two specific GmR variants (D = 0.328, in the minimum D group, and D = 0.778, in the large D group. See Supplementary Data 1), as these variants will be moved into the small D group as a result of the tRNA overexpression. Via homologous recombination (ClonExpress Ultra One Step Cloning Kit，Vazyme，C115), we inserted the argWCCT tRNA gene, driven by its native promoter, into the respective GmR plasmids, resulting in two tRNA-overexpressed strains. Primers used are listed in Supplementary Table 4.

We repeated the assessment for optimal D for different functional importance of the gene using AmpR, the ampicillin resistance gene. Specifically, we replaced the AmpR region native to the pUC57-GmR vector by the Mobilized Colistin Resistance-1 gene (MCR-1, see Supplementary Table 1), which confers polymyxin resistance, via homologous recombination. Next, the 47 synonymous variants of AmpR (Supplementary Data 2) were synthesized and similarly used to replace the GmR region in the vector (Supplementary Table 1) through homologous recombination. Using polymyxin as a selection marker, we successfully obtained 47 strains expression synonymous variants of the AmpR. Primers used are listed in Supplementary Table 4.

Fluorescence reporter gene assays
To evaluate the cis-regulatory effect of GmR’s codon usage on GmR expression, as well as the trans-regulatory effects on other endogenous genes within the cell, the 53 strains were first cultured in LB media supplemented with 50 μg/ml ampicillin (Sangon Biotech, A100304-0005) and 0, 10, 20, 30, 40, 60, 70, 80, 100, 120, 140, 160, or 180 μg/ml gentamicin (BioFroxx, 1146GR001). Then the expression levels of GmR-mCherry and YFP proteins were measured in >300,000 cells in log-phase growth for each of the 53 strains using flow cytometry (Attune N×T, Life Technologies). The fluorescence signals of mCherry and YFP were evaluated with a 15 nm bandpass centered on 620 nm and a 30 nm bandpass centered on 530 nm, respectively. The experiments were conducted with three biological and three technological replicates. E. coli cells with mCherry and YFP fluorescence signals ten times greater than those of the MG1655 negative control cells were retained for further analyses. The forward scatter (FSC), as well as the fluorescence signals of mCherry and YFP were measured for all cells. The expression levels of the fluorescent proteins were determined by dividing their fluorescence signals by FSC. The fluorescence levels were not measured for the AmpR strains since ampicillin compromises membrane integrity and causes cell content leakage, making fluorescence levels unreliable for determining the real cost and payoff.

Growth Rate Measurements
A spectrophotometer was used to measure the growth rates of the 53 strains cultured in LB media supplemented with 50 μg/ml ampicillin and 0, 10, 20, 30, 40, 60, 80, 100, 120, 140, 160, or 180 μg/ml gentamicin. Cells were cultured in LB medium at 37℃ and 250 rpm overnight, then diluted to OD = 0.15~0.2 and transferred to 96-well plates. Growth was monitored on Epoch2 microplate reader (BioTek) at 37℃ for 12 h by measuring absorbance (OD ) every 10 min. The doubling time (DT) was calculated by contrasting the starting and ending OD of a given time period via the following equation62
𝑙𝑛(2)
DT= (2)
(
ln(OD𝑡2)−ln(OD𝑡1)
)
𝑡2−𝑡1
Here, OD and OD are, respectively, the starting and ending OD of the focal time period from t to t . Considering the fact that different strains and gentamicin selections might exhibit different exponential growth phases in terms of OD ranges, it is imperative to identify the exponential growth phase of each strain using a unified, subjective method. To this end, for each sample, all combinations of OD measurements in three or more consecutive time points within the range of 0.2 < OD < 0.6 were fitted with linear models using the “lm” function in R. A combination with the most significant P value and a R2 value among the top 40% of all combinations is considered to be best described by exponential growth. The minimum DT in this exponential phase was used as the DT of the sample under consideration. The DT of each strain was measured in three biological and three technological replicates. Finally, the relative Wrightian fitness of each strain is calculated by dividing the mean DT of all strains cultured in gentamicin-free LB by the mean DT of the focal strain. This same procedure was used to measure the growth rate of AmpR strains cultured in LB media supplemented with 3 μg/ml polymyxin and 100, 2000, or 4000 μg/ml ampicillin.

The relative copy numbers of GmR plasmid and the abundance of YFP mRNA
The 53 GmR strains were lysed as previously described 63,64. Specifically, the strains were cultured overnight at 37°C and 250 rpm in 4 ml of LB medium containing 50 μg/ml of ampicillin and 0 or 60 μg/ml gentamicin until saturation. They were diluted to an OD of 0.15~0.2 and then grown to the logarithmic phase. After centrifuging the bacterial cultures, the supernatants were discarded. A sterile PBS solution (pH 7.4) was used to resuspend the cell pellets to an OD of 1. Subsequently, cells were lysed, via 95ºC incubation for 10 minutes followed by immediate freezing in −20°C. The β-lactamase gene (bla) on the plasmid and D-1-deoxyxylulose 5-phosphate synthase gene (dxs) on the E. coli chromosome were amplified by qPCR, using 3 μl suspension of the lysed cells, iTaq Universal SYBR Green Supermix (Bio-Rad, 1725122), and previously described65 primers (Bla-F and Bla-R, Dxs-F and Dxs-R, listed in Supplementary Table 4). The relative copy number of the plasmid was determined by the qPCR-based abundances of bla relative to that of dxs.

To determine the YFP mRNA abundance, total RNA samples were isolated from the log phase of all strains in LB media using the Bacteria RNA Extraction Kit (Vazyme, no. R403) according to the manufacturer’s instructions. Complementary DNAs were reversely transcribed from RNAs using PrimeScript RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, no. RR047A) according to the manufacturer’s instructions. Finally, RT–qPCR was carried out with cDNAs and primers (Supplementary Table 4) using the iTaq Universal SYBR Green Supermix 1000 Rxn (Bio-Rad, no. 1725122) according to the manufacturer’s instructions, measured as abundance relative to that of dxs.

Ribosome profiling (Ribo-Seq)
Single colonies of E. coli strains were inoculated and cultured overnight at 37°C with shaking at 220 rpm. They were then subcultured into fresh medium and incubated under the same conditions (37°C, 220 rpm) for 4–5 hours until the OD reached 0.6–0.7. A portion of the culture was tested for mycoplasma using the Mycolor One-Step Mycoplasma Detector (Vazyme, no. D201). Mycoplasma -free samples were treated with chloramphenicol (Selleck, no. S1677) at a final concentration of 100 μg/ml. The cultures were then incubated again for 10 minutes under the same condition. After cooling on ice, the cells were centrifuged at 4°C and 3200 × g for 10 minutes to collect cell pellets. After being flash frozen in liquid nitrogen for 10 minutes, the pellets were sent to Novogene for Ribo-Seq library preparation and sequencing.

Raw Ribo-Seq data were processed by removing adapters using Cutadapt66, aligning to the reference genome (GCF_000005845.2) and corresponding GmRs using STAR67, and quality-controlled with RiboParser68 to determine the read (ribosome-protected fragment, or RPF) counts per codon for each gene. Translational activity (TA) of a gene was defined as its average number of RPFs per codon after excluding the 20% highest and 20% lowest RPF values. This trimming removes extreme pauses and provides a length-normalized indicator for the translational activity of the gene, which combines both mRNA abundance and translational initiation rate per mRNA molecule. The codon consumption of a gene was calculated by multiplying the number of codon occurrences in the gene by the corresponding gene’s TA. The relative decoding time for each codon was calculated in two steps. First, the normalized ribosome density of a codon in a gene is calculated by dividing its average RPF counts within the gene by the gene’s TA. Second, the relative decoding time was calculated by averaging the normalized ribosome density of the same codon from different genes. Finally, the relative decoding times estimated from the GmR-expressing strains were divided by those estimated from the wildtype strain to reflect the changes in decoding time caused by the expression of GmR (Figure S3).

Genomic analysis for model organisms
Four model organisms were analyzed, namely E. coli, S. cerevisiae, Drosophila melanogaster, and Homo sapiens. Annotations and sequences of protein-coding genes are downloaded from EnsEMBL69 and NCBI70. The gene expression levels were derived from previous studies71-74. Using only the longest transcript with expression levels available for each gene, we obtained 3832, 5520, 13555, and 17212 genes for E. coli, S. cerevisiae, D. melanogaster, and H. sapiens, respectively. For the calculation of D values of genes in each species, only the top 10% highly expressed genes are considered for the X of the host (see above). Using all genes or tRNA expression levels would not change the conclusion.

The growth rate of E. coli gene deletion strains was defined as the average values of strains cultured in 30 environments obtained from a previous study29. The growth rate of S. cerevisiae gene deletion strains was defined as the average values of strains cultured in 11 environments (excluding UV-treatment, which is lethal to most deletion strains) from www.yeastphenome.org31. In both species, the relative growth rate was estimated by dividing the growth rate of a gene deletion strain by the growth rate of the wildtype strain. We then calculated the functional importance of each gene in E. coli and S. cerevisiae by subtracting 1 (the relative growth rate of the wildtype strain) from the corresponding gene deletion strain.

Analysis of mutation accumulation lines
The propensity of D-increasing mutations in the mutation accumulation (MA) lines is compared with its random expectation. Single nucleotide synonymous mutations occurring in MA lines are extracted from previous reports34,36,37. By comparing the D value of the mutant sequence to that of the wildtype, each mutation is classified as either D-increasing or D-decreasing. Considering that selection intensifies as gene expression level increases32,38, we restricted the analysis to m highly expressed genes above a certain threshold, and calculated the observed fraction of D-increasing mutations in MA lines. To approximate the random expectation of the observed fraction, for each highly expressed gene with n single nucleotide mutations, n other single nucleotide mutations were randomly selected from the focal gene. We required that the randomly selected mutations matched the observed mutations in terms of mutation types (one of A/T→C/G, A/T→G/C, A/T→T/A. C/G→A/T, C/G→G/C, C/G→T/A). We combined these random mutations from all m highly expressed genes into an “expectation” set of mutations. Using 10,000 such expectation sets, we estimated the random expectation of the fraction of D-increasing mutations. To demonstrate that this test has sufficient power to detect potential enrichment of D-increasing mutations in MA lines, we extracted the m genes with the smallest D values. Such constructed “positive control set” were subjected to the aforementioned observed (the positive control set) versus expectation (expectation sets reconstructed according to the positive control set) test for the enrichment of D-increasing mutations (Figure S10a/b). Another “negative control set” is constructed (Figure S10c/d) by a group of weakly expressed genes (defined as 20% least expressed genes) that matched the m highly expressed genes in terms of the number of genes, number of mutations of each type in a gene and GC content (±15%) of a gene. Because they have not undergone sufficiently strong translational selection, they are considered "negative controls" in the sense that they should show no enrichment of D-increasing mutations compared to random expectation.

Analysis of the fitness effects of synonymous mutations
A previous study1 measured growth rates in YPD for yeast strains containing one of 1866 synonymous mutations in 21 representative genes. The growth rates were obtained from this previous study. Mutations that appeared to be deleterious in all four experimental replicates were considered deleterious mutations in our study.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an increasingly important line of investigation. They report that their hypothesis challenges the conventional\nassumption that codon usage bias (CUB) is evolutionarily selected to represent the tRNA supply pool. The authors state that\nwhile matching the CUB to the tRNA pool increases the expression levels of the gene, it depletes the tRNA pool and\nimpacts the cell's fitness. In essence, the authors claim that the degree of the CUB-tRNA match depends on the functional\nbenefit vs fitness trade-off arising from the role of the gene in the cell. The more important the role of the gene in cell\nfunction/survival, the better the CUB-tRNA match. To this end, the authors create an Euclidean-distance metric “D”\nrepresenting the difference between a CUB and the representative tRNA pool. The closer “D” is to zero the more identical a\ngiven gene is coded based on the cellular tRNA pool.\nThe authors synthesized 3 groups of small, medium, and large “D” scored sequences of synonymous gentamicin resistance\ngenes, GmR. They fused mCherry to GmR and measured the functional payoff of the gene as the expression levels of GmR\ngene as a function of mCherry expression. They quantified changes in a secondary YFP gene to represent translation cost\nand changes in the exponential phase doubling time of the cells to represent overall fitness cost.\nI have serious concerns about the assertions, assumptions, choice of experimental system, and data interpretation. The\nbiggest flaws are:\n1) The authors do not quantify tRNA pools, which makes the title & premise of the paper completely misleading.\n2) The authors use gentamycin, a translational inhibitor in a study to understand translation elongation. Much of their data is\nnot interpretable since gentamycin alters the translation that they are claiming to evaluate.\n3) The authors keep claiming that the conventionally accepted paradigm is the selective pressure on genes to decrease the\nmismatch between their demand and tRNA pools unidirectionally. No citation is ever provided. Further, this statement is a\nred herring since this view is not widely accepted and therefore there is nothing to disprove. There is no singular selective\npressure for all genes. Different genes have different selective pressures on them, depending on their function, when they\nare expressed, and in what abundance they are needed. A heat shock protein, which is needed under highly specific\nconditions will not have the same selective pressure on it as a glycolytic enzyme. Thus, the basis of their paper in disproving\na conventionally held view is completely unfounded.\nThe specifics are below:\nMajor Comments:\nLine 77-94: Why was Gentamicin antibiotic selected for this study? This is an interesting choice of antibiotic, given that\ngentamicin is known to exacerbate codon-tRNA misreads and lead to the incorporation of the wrong amino acids as well as\nslow down elongation. It would be good to discuss how reducing the stringency of tRNA-codon interactions may change the\nresults when studying codon influence on gene expression. Could the use of gentamicin create non-linear effects in the\ndata, given that an optimal match to tRNA becomes more important in the presence of an antibiotic that exacerbates errors\nwith mismatched codon-tRNA? Do you have data to support or refute this? This is a near-fatal flaw of the paper. How do you\ninterpret the role of translation elongation when the reporter system is based on perturbing translation elongation? A reporter\nsystem like ampicillin or another antibiotic that does not affect transcription or translation should have been used. It’s almost\nimpossible to interpret the results in light of the fact that gentamycin is a translation inhibitor.\nLine 100: The authors did not measure tRNA supply at all. They only used CUB of highly expressed genes! The title and\npremise of this paper are completely misleading and is not what supported by the data presented in the paper. Based on the\ndefinition of D, it seems what the authors are calling tRNA supply is similar to the traditional codon adaptation index (CAI) for\nhighly expressed genes. There are numerous lines of evidence that CAI is not a measure of tRNA pools, nor is there are\nvery strong correlation between the two. tRNA levels have previously been reported for E. coli (10.1006/jmbi.1996.0428,\n10.1080/15476286.2020.1790871). The authors should either measure tRNA levels directly to show it matches CAI or\nremove any mention of the metric they are comparing with D is tRNA pool. The use of CAI as a proxy for tRNA pools is\ncompletely unacceptable.\nLine 111: Placing the YFP and GmR genes on the same plasmid as a bidirectional promoter system (convergent system)\ncan lead to confounding effects during unwinding and transcription. This implies that the transcription of one gene can\nimpact that of another. Did the authors ensure/prove that their GmR constructs only have trans effects on the YFP\ntranslation? Did they prove that the GmR constructs do not cause cis effects on the transcription of YFP at the DNA level? It's\npossible for instance that re-codes can affect plasmid structure and thus gene expression of other ORFs on the plasmid, as\nwell as supercoiling, etc. The authors need to include sufficient controls (e.g., no RBS or moving YFP to the\ngenome/different plasmid) to show that their results are indeed trans effects and not cis effects.\nLine 358-373: A key issue with this paper concerns the calculation of the distance “D” metric.\nThe authors do not use a metric that relates to tRNA supply. Their metric D is the difference between demand by highly\nexpressed genes (Y, a genomic frequency of a specific codon used in a sequence) and supply (X, the frequency of the same\ncodon only in highly expressed genes). Effectively, they create a weight for every codon by subtracting the two frequencies\nabove (Y-X) and then they calculate the geometric mean of all the codons and their weights. The authors claim that X or the\n“tRNA” supply is a representation of “effective” codon supply encompassing tRNA expression levels, tRNA aminoacylation\nrate, etc. But I disagree. They practically created another Codon Adaptation Index in a different way. They did not calibrate\ntheir “D” score meaning that lower D levels may not really be representing the lower tRNA supply. Effectively, for “D” to be\nuseful it has to show that it represents the entire range of tRNA supply somewhat linearly for the authors to be able to claim\nthat lower “D” levels close to zero indeed corresponds to the maximum concordance of demand (codon being read) to\nsupply (tRNA level). I am pessimistic that their metric represents tRNA supply the same way across its entire range.\nLine 49-51: “the prevailing hypothesis …” This sentence is poorly written and hard to understand.\nAuthors need to provide references for the prevailing hypotheses asserting that “translational selection is towards an optimal\nCUB”, “Deviations are due to mutational bias or genetic drift”.\nLine 57-60: It would be easier to understand the author’s claims if they provided a supplementary figure for this assertion\ninstead of citing some papers.\nLine 72-76: The authors’ definition of the “conventional model” is not supported. This statement is a red herring giving the\nauthors the opportunity to disprove a hypothesis that is already quite ridiculous. In no way is the prevailing/conventional\nwisdom that genes drift toward high CAI. Indeed, if that were the case, there would be no codon usage bias among different\ngenes within an organism. The drift toward or away from matching tRNA supply would be highly dependent on 1) how\nimportant the protein is to the cell’s survival, 2) what concentration of the protein is required for cell function and 3) under\nwhat conditions is the protein essential. To state that drift toward a low mismatch between tRNA supply and demand is\nridiculous. The entire premise of this manuscript is based on a red herring.\nLine 73-76: There is unlikely to be a single answer for this since codon bias toward high and low D exists of different\nendogenous genes. Again, this depends on the cellular need for the specific function under the conditions tested. There is\nnot likely to be a single answer. This general premise is highly flawed. The selective pressure is different for different genes.\nLine 109-117: Authors are not controlling for changes in 5' region that are translation initiation effects. This confounds\nresults. It's not clear whether the effects are due to translation initiation (controlled by mRNA structure & ribosome binding) or\nelongation as the authors want to claim. Currently, the authors do not prove that the effects are not from transcription which\ncan misrepresent the translation results. The authors need to show that the results of expression level changes cannot be\nexplained by translation initiation changes.\nSimilarly, the authors do not control for transcription rates. In bacteria, transcription and translation are coupled. The authors\nshould show that their data is not explained simply by changes in mRNA levels of their reporter if they want to claim that it is\nall translation elongation effects.\nLine 192-195: Authors do not show the data for the KanR gene replacement experiments that is necessary for proving their\nclaim. Needs more information on how functional importance was quantified. Furthermore, Spearman’s correlation of -0.08 is\nextremely weak almost zero. Authors need to show the data spread used for their correlations. Just providing a correlation is\nunacceptable.\nLine 207-213: This statement is not novel, is intuitive and expected, and is not a contribution of this work. It is already well-\naccepted that high transcription and translation rates are highly burdensome to cells. See PMID: 27989436.\nLine 217-219: What is the threshold authors set for expression and why? The authors need to show the expression data\ndistribution used for obtaining their correlational data in Figure 5.\nLine 243 and 282-284: E. coli is a gut microbe. This means that evolutionary pressure for optimizing its genes has largely\nbeen driven by the gut microenvironment to match its demands. The authors culture their E. coli in LB media. This means\nthat there will likely be a mismatch in the mutational accumulation experiments relevant to D-increasing or D-decreasing\neffects. It is incorrect to assume that the selective pressure in lab conditions in LB reflects true selective pressure for E. coli.\nLine 303: Authors need to clarify what they mean by “waste of tRNAs for unpreferred codons”.\nRegarding Figure 3.b\n- Why does at high Gm concentrations the small group (theoretically with less GmR protein than the minimum group) exhibit\na higher functional payoff?\n- Why is everything normalized to the large group? This implies you are assuming unidirectionality in codon-match and\nfunctional payoff.\n- The data contradicts the authors’ assumption that as the importance of a gene increases, D (tRNA-CUB supply-demand\ndifferential) decreases. At concentrations above 140 ug/ml we see that the small group is preferred. (In text line 154 you say that\nonly 68.9% of the variations are explained by this assumption).\nFigure 3.c-e\n- It is unclear what the data points here represent. Are they average or individual single recodes? Why is there only 12 points\nhere? Shouldn’t you have the minimum and small groups here (38 points)? Authors need to provide information about their\nstatistical test in the figure legend.\nFigure 3.f\n- Relevant to result interpretations of line 160-163:\n- I disagree that at concentrations above 120 ng/ul minimum D is selected. Mainly because the statistic is weakly probable.\nFrom 140 ng/ul all groups most likely have the same relative fitness. This contradicts your statement that minimum D is best\nat the highest functional importance (highest concentrations).\nFigure 4\n- This correlation does not support any direct like to the hypothesis that genes with smaller D scores have a higher functional\npayoff. Because the function is related to the direct role of the gene and not how it is coded. In order to support your\nhypothesis, you need to have direct fitness data where you measure the fitness cost associated with recoding a specific vital\ngene and replacing it with the native gene to determine the link between functional payoff and translation cost.\n- With regards to line 211-213 An important gene does not have to be highly expressed. You cannot establish the\nimportance of genes by only looking at highly expressed genes and measuring expression indirectly by the D index.\nFigure 5\n- In regards to line 216: Because D and expression level are anticorrelated does not imply anything about translational cost.\nYou are only looking at functional payoff. Translation cost can only be measured experimentally.\nQuestion on the methods:\n- How did the authors control for the sample size reduction of highly expressed genes? Correlation scores tend to shrink in\nabsolute values as sample sizes grow. Conversely, a strong correlation at a small sample size may not be true when the\nsample of the correlation is increased.\nFigure 6\n- While this evidence does support the authors' claim that a slight mismatch with the tRNA supply is selectively maintained,\nthe authors do not quantify the strength of the D-changing mutations. In other words, what is the distribution of D score\nchange in each D-increasing/-decreasing mutation group?\nAnother note regarding the authors' quantification of cellular tRNA:\nIs the index for quantifying the cellular tRNA pool linear in its entire range? Perhaps if this index becomes non-linear in\neither boundary it may falsely inform that D decreasing mutations are deleterious at highly expressed genes.", "reviewer_2": "The manuscript submitted by Chen, Liu, Zhou and co-authors is dedicated to the study of translational selection.\nSpecifically, the authors propose a new model to explain the evolution of codon usage bias (abbreviated CUB; this can be\nroughly defined as the departure from a random distribution among synonymous codons), according to which the optimal\ncodon usage does not (always) perfectly match the abundance of tRNAs in the cell. According to this model, a mismatch\nbetween the codon usage bias of a given gene and the tRNA supply can be beneficial because it reduces the cost imposed\nby the translation of the gene. In other words, a highly expressed and highly translated gene would be costly because its\ntranslation would hamper the translation of other genes.\nTo support this model, the authors mainly rely on genome editing experiments performed in E. coli. Specifically, they insert\ncassettes carrying a gentamicin-resistance gene (GmR) with different codon usage biases. The GmR gene is fused to a\nmCherry coding sequence, which enables the authors to measure the expression level of the GmR-mCherry fusion through\nfluorescence measurements. The cassettes also contain a YFP gene in the opposite orientation. The expression level of the\nYFP gene, again estimated through fluorescence measurements, enables the authors to estimate the translational cost of\nthe CUB of the GmR gene as the reduction in YFP protein levels. The authors perform these experiments at different\nconcentrations of gentamicin and they also measure the growth rate of the different E. coli clones, thus obtaining a more\ndirect estimation of the fitness effect of the different constructs. In addition to these main experiments, the authors also\nanalyze previously published data for mutation accumulation or gene replacement experiments in E. coli and yeast, as well\nas codon usage data for various species.\nI find the model proposed by the authors to be interesting and reasonable. It brings an alternative to the classical model\nwhich posits (even in species where translational selection is efficient, that is, species with high effective population size\nsuch as E. coli, yeast or fruitfly) translational selection cannot optimize the codon usage bias to perfectly match the tRNA\nabundance of the cell because of genetic drift. However, it must be stated that the classical model already incorporates the\nidea of a difference between gene categories (the translational selection coefficient is expectely to vary among genes,\ndepending on the optimum protein quantity produced by each gene, and for some genes the translational selection\ncoefficient is indeed expected to be very low). More importantly, I am concerned about the validity of several assumptions\nmade in this article and of several analyses performed in this article.\nFirst, the authors state that the D metric that they define is a measure of the mismatch between the CUB of a gene and the\ntRNA supply of the cell. It is in fact a measure of the mismatch between the CUB of a given gene and the CUB of the highly\nexpressed genes (Methods, lines 358-373). The authors assume that the CUB of highly expressed genes has been\noptimized to match the tRNA supply, although they then affirm that the optimal CUB should not perfectly match the tRNA\nsupply, even for highly expressed genes. This is circular and illogical. There are better ways to approximate the tRNAs\nsupply without assuming that the CUB for highly expressed genes has been optimize to match it. If the authors are unable to\nmeasure directly tRNA abundance, they could use an approximation based on tRNA gene copy number, taking into account\nwobble rules etc. As this affirmation is the basis of the new model proposed by the authors, it is very important to make sure\nthat D indeed measures the distance compared to the CUB that would perfectly match the tRNA supply of the cell.\nSecond, the authors measure the translational cost of the CUB of the GmR artificial constructs by measuring the\nfluorescence level of the exogenous YFP present in the construct. However, they do not ascertain that variations in YFP\nfluorescence levels are indeed to variations in YFP translation, rather than for example in variations in YFP gene\ntranscription. They also do not show that the variations in YFP fluorescence are causally connected to the translation of the\nGmR-mCherry gene. The authors and others have shown before that variations in CUB can also affect the mRNA level of\ngenes, mainly due to mRNA stability. Additional effects ca be imagined at the transcriptional level, given the overall\ndifferences in base composition of the CUB variants. As the YFP gene is inserted in close proximity to the GmR-mCherry\ngene, on the reverse strand, transcriptional interference can occur between the two genes. Thus, their claim of a translational\ncost is not sufficiently supported. Furthermore, the translational cost should be evaluated on endogenous genes. As this is\nan important claim of the manuscript, it should be properly verified, by evaluating mRNA levels and translation rates in the\ndifferent constructs (for example with RNA sequencing and riboosome profiling assays).\nThird, the fact that the different E. coli strains that they used could display variations in plasmid copy number is quite\nconcerning. I am not convinced that the lack of correlation between the copy number and D is sufficient to ensure that the\nresults are not biased. The authors should verify that the plasmid copy number variation is not a confounding factor for\nfitness analysis and for translational cost analyses.\nFourth, the fact that the authors observe a positive correlation between D and expression levels among the most highly\nexpressed genes (top 10%, top 5% and top 3%) in all 4 analyzed species (Figure 4) is puzzling. In human, translational\nselection is highly inefficient, unlike in the other 3 species. The fact that the same effect is observed in human, E. coli,\nyeast and Drosophila rather points to a more general effect, potentially a confounding effect of base composition. This should\nbe examined carefully.\nFifth, it seems that the mutation accumulation analysis was performed using both synonymous and nonsynonymous\nmutations (only nonsense mutations seem to have been excluded, Methods, line 464 to 466). This does not make sense - if\none wants to measure the effect of the D, only synonymous mutations should be considered. It is also not clear why this\nanalysis was only performed on the 10% most highly expressed genes (especially given that D measures the distance to the\ncodon usage of the most highly expressed genes). Here, a comparison between gene expression classes would be useful,\nto test the hypothesis that the CUB of the highly expressed genes has indeed been optimized. Here again, D should be\ncomputed with respect to the codon usage that perfectly matches the tRNA abundance (not to the codon usage of highly\nexpressed genes).\nFinally, the authors might be interested in this preprint, which shows that translational selection is surprisingly rare in\nmetazoans: https://www.biorxiv.org/content/10.1101/2024.07.22.604600v1 .", "reviewer_3": "The paper titled “A slight mismatch between a gene’s codon usage and the cellular\ntRNA supply is beneficial” describes a very thorough work with a thought-provoking notion. Unlike the current dogma in the\nfield, which assumes that optimal translation efficiency is achieved when codon usage perfectly matches the tRNA pool of\nthe cell, a slight mismatch might be optimal. A mechanism based on feedbacl on the tRNA pool is suggested. A proposed\nmechanism is based on an alleged trans-regulation and tRNA depletion.\nWhile previously there were some indications for this effect in unneeded genes, now they look for genes with fitness effect,\neither a drug resistance gene or natural genes. This is a major step forward as it allows them to balance costs against\nfunctional payoff.\nThis is an excellent work. It is thought-provoking, very well done and deep in it evolutionary approach. It is also very well\nwritten and technically sound. I would ultimately recommend publication in NC. Yet, if I am contribute to the robustness of the\nconclusion and to the discussion that this paper will surely evoke, here are some important points.\nI’m accepting to sign the report with my name: Yitzhak (Tzachi) Pilpel.\n1. A critical experiment could be to manipulate the tRNA pool and then test, what I believe is a prediction of their theory. That\nfor different tRNA pools, alternative members of the collection of 53 synonymous variants would optimal. This is a clear\nprediction of their model since for each tRN pool other variants might be at optimal D. If they showed that I would be\nconvinced better. If the result would show that same variants are good at different tRNA pools then their conclusion is not\nsupported. In other words, they have manipulated so far the codons – the demand, but never the tRNA – the supply.\n2. The proposed mechanism that could mediate their effect is that of tRNA depletion, “In theory, the same logic applies to\nendogenous genes native to the cell, since they also trans-regulates the expression of other genes via tRNA depletion” –\nRef #23. First, I don’t think reference #23 shows tRNA depletion. This is a misunderstanding of what was done in Frumkin et\nal. The increase in tRNA supply in ref #23 was very different, first because it was done artificially by the researchers (not\nspontaneously by the cell), and second, because the tRNA levels themselves were not changed, but rather tRNA anticodons\nwere mutated. So it’s incorrect to ascribe to ref#23 a trans-regulation of the tRNA pool. That said, I think we are left without a\ntrue reference or a proof, that the tRNA pool is indeed trans-regulated in response to change in codon demand. So I can see\nseveral options: 1. The authords measure tRNA pool with some existing tRNA seq protocols and show the effect; 2. Find a\nreference that truly shows this effect, or 3. Acknowledge that we don’t know if such effect indeed occurs at all\n3. Their calculation of the D measure – the extent of deviation from the tRNA pool and codon usage is critical, especially\nsince the whole argument is based on the difference between D that approach zero and D that is slightly above zero. My\nproblem is with their estimation of the tRNA pool which is needed for equation in line 358. Particularly they say “Xij\nrepresents the tRNA supply of the cell, approximated by the fraction of codon j among the synonymous codons of amino acid\ni for the highly expressed genes in the transcriptome of the cell” ref #47. I’m not sure if and how ref #47 provides us with a\ntRNA pool. So a few comments on this: 1. I’d like to see their results reproduced with other, more conventional definitions of\nthe tRNA pool. Again, they could do tRNA seq. There’s also the classical tAI – based tRNA pool deduced from tRNA gene\ncopy number. Is their D optimized to be slightly above zero also when using these more conventional estimations of the\ntRNA pool? The way in which ref #47 appears to be estimating the tRNA pool is affected by codon demand, this is\npotentially problematic.\n4. Their estimates for genes functional importance in case of natural genes is based on fitness upon replacement with\nKan_R. I’m not sure if this measure is adequate because it represents a total deletion of a gene, while here we ask about\nreduction in its translation efficiency. I’m not sure how to solve this problem, but I’d be better convinced if they had additional\nmeasures of genes’ functional importance. Perhaps by employing bioinformatics of evolutionary conservation of genes they\ncould generate an estimation of functional importance but not open total deletion."}, "author_response": "Point-by-point response to reviewers’ comments\nReviewer #1:\nOverall Comment:\nThis is an increasingly important line of investigation. They report that their hypothesis\nchallenges the conventional assumption that codon usage bias (CUB) is evolutionarily\nselected to represent the tRNA supply pool. The authors state that while matching the\nCUB to the tRNA pool increases the expression levels of the gene, it depletes the tRNA\npool and impacts the cell's fitness. In essence, the authors claim that the degree of the\nCUB-tRNA match depends on the functional benefit vs fitness trade-off arising from\nthe role of the gene in the cell. The more important the role of the gene in cell\nfunction/survival, the better the CUB-tRNA match. To this end, the authors create an\nEuclidean-distance metric “D” representing the difference between a CUB and the\nrepresentative tRNA pool. The closer “D” is to zero the more identical a given gene is\ncoded based on the cellular tRNA pool.\nThe authors synthesized 3 groups of small, medium, and large “D” scored sequences of\nsynonymous gentamicin resistance genes, GmR. They fused mCherry to GmR and\nmeasured the functional payoff of the gene as the expression levels of GmR gene as a\nfunction of mCherry expression. They quantified changes in a secondary YFP gene to\nrepresent translation cost and changes in the exponential phase doubling time of the\ncells to represent overall fitness cost.\nI have serious concerns about the assertions, assumptions, choice of experimental\nsystem, and data interpretation. The biggest flaws are:\n1) The authors do not quantify tRNA pools, which makes the title & premise of the\npaper completely misleading.\n2) The authors use gentamycin, a translational inhibitor in a study to understand\ntranslation elongation. Much of their data is not interpretable since gentamycin alters\nthe translation that they are claiming to evaluate.\n3) The authors keep claiming that the conventionally accepted paradigm is the selective\npressure on genes to decrease the mismatch between their demand and tRNA pools\nunidirectionally. No citation is ever provided. Further, this statement is a red herring\nsince this view is not widely accepted and therefore there is nothing to disprove. There\nis no singular selective pressure for all genes. Different genes have different selective\npressures on them, depending on their function, when they are expressed, and in what\nabundance they are needed. A heat shock protein, which is needed under highly specific\nconditions will not have the same selective pressure on it as a glycolytic enzyme. Thus,\nthe basis of their paper in disproving a conventionally held view is completely\nunfounded.\nResponse:\nPlease accept our sincere thanks for your time, your thoughtful reading, and very critical\ncomments. These suggestions are taken very seriously as we believe that addressing\nthem would significantly improve our study. Since these topics are detailed below in\nindividual comment paragraphs, we will address them in a point-by-point manner\nbelow to avoid clumsiness and repetition. More specifically, the concern (1) is\naddressed in our responses to Comments 2 & 4 & 22, concern (2) in Comment 1, and\nconcern (3) in Comment 5 & 7 & 8 below.\nComment 1\nLine 77-94: Why was Gentamicin antibiotic selected for this study? This is an\ninteresting choice of antibiotic, given that gentamicin is known to exacerbate codon-\ntRNA misreads and lead to the incorporation of the wrong amino acids as well as\nslow down elongation. It would be good to discuss how reducing the stringency of tRNA-\ncodon interactions may change the results when studying codon influence on gene\nexpression. Could the use of gentamicin create non-linear effects in the data, given that\nan optimal match to tRNA becomes more important in the presence of an antibiotic that\nexacerbates errors with mismatched codon-tRNA? Do you have data to support or\nrefute this? This is a near-fatal flaw of the paper. How do you interpret the role of\ntranslation elongation when the reporter system is based on perturbing translation\nelongation? A reporter system like ampicillin or another antibiotic that does not affect\ntranscription or translation should have been used. It’s almost impossible to interpret\nthe results in light of the fact that gentamycin is a translation inhibitor.\nResponse:\nThank you very much for your insightful comment. We would like to respond from\nthree angles below.\n1. “Why was Gentamicin antibiotic selected for this study?”\nTo test our model, we require a gene with modifiable functional importance in order\nto evaluate the relationship between functional importance and optimal codon usage.\nWe decided to use antibiotic resistance genes, whose functional importance can be\neasily altered by adjusting the antibiotic concentration in the culture media.\nAdditionally, the gene should be relatively short to reduce the cost of synthesizing its\nsynonymous variants and minimize the complexity of the experimental procedure. The\nGmR gene was selected from a variety of commonly used resistance genes due to its\nshort sequence length of only 531 bases.\n2. “Could the use of gentamicin create non-linear effects in the data, given that an\noptimal match to tRNA becomes more important in the presence of an antibiotic that\nexacerbates errors with mismatched codon-tRNA?”\nWhile we agree that Gentamicin may enhance the importance of matching tRNA\nsupply, the non-linearity cannot be solely attributed Gentamicin. In fact, most models\nof gene expression assume a nonlinear diminishing return for a gene’s functional output\n(e.g. PMID:20884723), and such nonlinearity is indeed observed for the fitness effect\nof increased expression (e.g. PMID: 27545349).\nMore importantly, while GmR function may explain the observed patterns of\nfunctional payoff, it cannot explain the translational costs. In particular, if we consider\na perfect match between GmR codon usage and the cellular tRNA supply, its resultant\nincreased GmR expression and enhanced effect on alleviating Gm's toxicity, then YFP\nexpression should increase, or, in other words, lower translational cost in strains with\nGmRs of smaller D. This is the opposite of what we found for translational cost, which\nwas generally higher in strains with GmRs of minimum D (Fig. 2d). This pattern of\ntranslational cost, however, can be explained by our proposed model of tRNA depletion\ncaused by excessive expression of GmRs of minimum D. Hence, the data refute the\npossibility of a GmR-specific artifact, and they are more consistent with our proposed\nmodel.\n3. “A reporter system like ampicillin or another antibiotic that does not affect\ntranscription or translation should have been used”\nPer your suggestion, we designed and synthesized 47 synonymous variants of\nAmpR, the ampicillin resistance gene, and constructed expression cassettes by\nreplacing GmR in Figure 2b with one of these variants, resulting in 47 synonymous\nAmpR strains. These 47 strains were then tested for fitness at various ampicillin\nconcentrations. As a result, we observed a similar pattern to that observed with GmR.\nThat is, at low/intermediate/high concentrations of ampicillin (100/2000/4000 ug/mL),\noptimal fitness is achieved by cells expressing AmpR with large/small/minimum D. We\nhave now added this result as Supplementary Figure S5, and is pasted below. Please\nnote that the fluorescence levels were not measured for the AmpR strains since\nampicillin compromises membrane integrity and causes cell content leakage, making\nfluorescence levels unreliable for determining the real cost and payoff.\nFigure S5. The optimal D for AmpR\n(a) A total of 47 synonymous AmpR variants with different D values were\ndesigned, which were categorized into three equal-sized groups of minimum,\nsmall and large D. (b) The relative fitness of cells expressing different AmpR\nsynonymous variants in LB media supplemented with 100, 2000, or 4000\nμg/ml ampicillin. The statistical significance of two-tailed Wilcoxon Rank\nSum Tests are indicated as *: P<0.1, **: P < 0.05, or ***: P < 0.01.\nComment 2\nLine 100: The authors did not measure tRNA supply at all. They only used CUB of\nhighly expressed genes! The title and premise of this paper are completely misleading\nand is not what supported by the data presented in the paper. Based on the definition of\nD, it seems what the authors are calling tRNA supply is similar to the traditional codon\nadaptation index (CAI) for highly expressed genes. There are numerous lines of\nevidence that CAI is not a measure of tRNA pools, nor is there are very strong\ncorrelation between the two. tRNA levels have previously been reported for E. coli\n(10.1006/jmbi.1996.0428, 10.1080/15476286.2020.1790871). The authors should\neither measure tRNA levels directly to show it matches CAI or remove any mention of\nthe metric they are comparing with D is tRNA pool. The use of CAI as a proxy for\ntRNA pools is completely unacceptable.\nResponse:\nThank you for your comment. We will explain below the logic behind our estimation\nof tRNA supply, and then demonstrate that results based on tRNA quantification do not\nalter our conclusions. Please also see our response to the related Comment 4 below.\nThe logic behind our estimation of tRNA supply is twofold. First, the actual tRNA\nsupply is neither the copy number of the tRNA gene nor the abundance of tRNA\ntranscripts. Instead, it is the abundance of the ternary complex of aminoacyl-tRNA+EF-\nTu+GTP (or replace EF-Tu with EF1A in eukaryotes). Currently, the abundance of\nternary complexes for different types of tRNAs cannot be systematically measured. In\naddition, the kinetic parameters between different codon-anticodon pairs will also\ninfluence the \"effective\" tRNA supply, further complicating the estimation process.\nSecond, because of the first point, we resorted to the assumption that the cellular tRNA\nsupply and demand (codon usage) of the whole transcriptome have been evolutionarily\noptimized by natural selection to ensure a close match between supply and demand.\nThis assumption allows us to approximate the supply of tRNA by using the codon usage\n(the demand) of the top 10% expressed gene (which dominate the transcriptome) in the\ntranscriptome. It is important to note that we are using the transcriptome, not the\ngenome, which means the codon usage has been weighted by the expression of those\ngenes before being summed up. Furthermore, this assumption of supply-demand match\nhas been empirically supported (PMID: 22479199 and PMCID: PMC12262485), and\nhas been used in our previous publications (PMID: 32123323 and 34977494). It is also\nimportant to note that a supply-demand match on a transcriptome level does not\nnecessarily imply a supply-demand match for individual genes, so this assumption does\nnot contradict our main conclusion.\nNevertheless, per your suggestion, we obtained tRNA expression levels based on\nhigh-throughput tRNA sequencing data from a previous publication (PMID: 31353208,\nGEO accession GSE128812), which were then scaled by codon-anticodon affinity and\nused as the tRNA supply. This metric is similar to the “adaptiveness value” (w) of a\ncodon, as used in the popular tAI metric (PMID:15448185). It was found that the\ngrouping (minimum, small, or large D) in Figure 3 remained unchanged, suggesting\nthat our conclusion is not affected by the change in the method used to quantify tRNA\nsupply. We have now included this result in Supplementary Figure S2, and it is pasted\nbelow.\nFigure S2. Grouping by D of GmR strains remains unchanged if high-\nthroughput tRNA sequencing data are used to estimate tRNA supply\nThe D values of designed GmR synonymous variants were calculated with\ntRNA supply based on either transcriptomic codon usage of the top 10%\nhighly expressed genes (x axis) or tRNA transcript abundance derived from\nhigh-throughput tRNA sequencing (y axis). The black lines indicate the\ngrouping, which is identical for both types of D values. Please note that each\ndot has up to three unique sequences as biological replicates (see Table S2).\nComment 3\nLine 111: Placing the YFP and GmR genes on the same plasmid as a bidirectional\npromoter system (convergent system) can lead to confounding effects during unwinding\nand transcription. This implies that the transcription of one gene can impact that of\nanother. Did the authors ensure/prove that their GmR constructs only have trans effects\non the YFP translation? Did they prove that the GmR constructs do not cause cis effects\non the transcription of YFP at the DNA level? It's possible for instance that re-codes\ncan affect plasmid structure and thus gene expression of other ORFs on the plasmid, as\nwell as supercoiling, etc. The authors need to include sufficient controls (e.g., no RBS\nor moving YFP to the genome/different plasmid) to show that their results are indeed\ntrans effects and not cis effects.\nResponse:\nThank you very much for your critical comment. We have included an insulator\nbetween the terminators of YFP and GmR to prevent transcriptional interference\nbetween the two genes. The Figure 2b pasted below has been updated to explicitly\nindicate the insulator.\nFurthermore, per your comment, we used quantitative RT-PCR to measure YFP mRNA\nexpression levels in the GmR strains growing in LB. We found that YFP mRNA\nexpression did not correlate with GmR’s codon usage (D value). This result has been\nadded to Supplementary Figure S9c and pasted below.\nFigure S9. Correlations between D and functional payoff, translational\ncost or fitness\n(c) The abundance of YFP mRNA in different strains, as determined by\nquantitative RT-PCR (see Methods), is not correlated with codon usage of\nGmR. The Spearman’s Correlation Coefficient ρ of the raw data, as well as\nthe corresponding P value, are indicated.\nComment 4\nLine 358-373: A key issue with this paper concerns the calculation of the distance “D”\nmetric.\nThe authors do not use a metric that relates to tRNA supply. Their metric D is the\ndifference between demand by highly expressed genes (Y, a genomic frequency of a\nspecific codon used in a sequence) and supply (X, the frequency of the same codon\nonly in highly expressed genes). Effectively, they create a weight for every codon by\nsubtracting the two frequencies above (Y-X) and then they calculate the geometric mean\nof all the codons and their weights. The authors claim that X or the “tRNA” supply is a\nrepresentation of “effective” codon supply encompassing tRNA expression levels,\ntRNA aminoacylation rate, etc. But I disagree. They practically created another Codon\nAdaptation Index in a different way. They did not calibrate their “D” score meaning\nthat lower D levels may not really be representing the lower tRNA supply. Effectively,\nfor “D” to be useful it has to show that it represents the entire range of tRNA supply\nsomewhat linearly for the authors to be able to claim that lower “D” levels close to zero\nindeed corresponds to the maximum concordance of demand (codon being read) to\nsupply (tRNA level). I am pessimistic that their metric represents tRNA supply the same\nway across its entire range.\nResponse:\nThis comment raises an imporant point regarding the estimation of tRNA supply in our\nstudy, which is greatly appreciated. Our response to Comment 2 above explains the\nlogic behind the metric D and demonstrates that our results are not affected by using\nhigh-throughut tRNA sequencing estimates of tRNA supply.\nWith regard to the more specific points in this comment, we would like to clarify\nthat we do not claim or assume that D linearly represents the tRNA supply or the tRNA\nsupply-demand mismatch. We do not make conclusions based solely on the linearity\nassumption, such as Pearson's correlation tests with D (we always accompany them\nwith Spearman's rank-based tests). The reasons for such non-linearity can be biological\n(the degree of mismatch is not linear to the phenotypic consequences) and/or technical\n(different amino acids might have different impacts, but we are merely using their\ngeometric average).\nMoreover, the D metric differs from CAI in several ways. For example, increasing\nCAI is selectively advantagous due to improved translational efficiency when more\nprefered codons are used (as assumed by its underlying model), whereas minimizing D\nis advantageous. Also, classic CAI calculation uses the genomic codon usage of the top\n10% highly expressed genes without using their expression levels as a weight to sum\nup codon usage across genes. As for D, we used the transcriptomic, not the genomic,\ncodon usage, which means the codon usage has been weighted by the expression of\nthose genes before being summed. Moreover, CAI is calculated as the geometric mean\nof the relative adaptiveness of each codon site. In other words, its final calculation step\ninvovles taking the L-th root of a number, where L is the length of the protein (or\nnumber of codons in the CDS). But for D, the final calculation step is usually taking\n18-th root, since there are only 18 groups of synonymous codons. This number of 18\nmay be further reduced if some amino acid is completely absent from the protein. There\nare other differences, such as CAI assumes the preferred codon is the best option,\nwhereas D assumes perfect supply-demand balance is the best option. We hope these\nexplanations have clarified D's distinct characteristics, particularly in comparison with\nCAI.\nComment 5\nLine 49-51: “the prevailing hypothesis …” This sentence is poorly written and hard to\nunderstand.\nAuthors need to provide references for the prevailing hypotheses asserting that\n“translational selection is towards an optimal CUB”, “Deviations are due to mutational\nbias or genetic drift”.\nResponse:\nThank you and we apologize for the confusion. We have now included two highly cited\nreferences that describe the prevailing hypothesis. For example, in [Hershberg and\nPetrov. Selection on Codon Bias. Annual Review Genetics. 2008] (1195 citations on\nGoogle Scholar), it is stated that\nExplanations for the existence of codon bias fall into two general classes.\nAccording to the selectionist explanation, codon bias contributes to the\nefficiency and/or the accuracy of protein expression and is thus generated and\nmaintained by selection. The mutational or neutral explanation, by contrast,\nposits that codon bias exists because of nonrandomness in the mutational\npatterns.\n…\nGiven the evidence that both mutational pressures and selection are\ninvolved in the phenomenon of codon bias, the current accepted model is the\nmajor codon preference model, also known as the mutation-selection-drift\nbalance model of codon bias. This model proposes that selection favors the\nmajor (or preferred) codons over minor codons. However, mutational pressure\nand genetic drift allow the minor codons to persist.\nWe have also cited another classic paper that proposes this prevailing model [The\nSelection-Mutation-Drift Theory of Synonymous Codon Usage. Michael Bulmer.\nGenetics. 1991] (1244 citations on Google Scholar), and also revised this sentence,\nwhich is pasted below. It is our hope that the reviewer finds this statement more\nunderstandable in explaining the current prevailing model.\nOther than some special cases affecting only limited neighboring regions5-\n8, the prevailing hypothesis for the widely observed CUB is the slection-\nmutation-drift model9,10, which proposes that translational selection favors an\noptimal CUB, whereas deviations from the optimal CUB persist as a result of\ngenetic drift and mutational bias.\nComment 6\nLine 57-60: It would be easier to understand the author’s claims if they provided a\nsupplementary figure for this assertion instead of citing some papers.\nResponse:\nThank you and we agree with you. Our claim is summarized in Figure 1, which is\npasted below.\nFigure 1. Two competing models for the translational selection on\nsynonymous codon usage bias.\n(a) The conventional model assumes unidirectional selection towards the\nminimal mismatch between a gene’s codon usage and the cellular tRNA\nsupply. (b) We propose a novel model in which the optimal CUB mismatches\nthe cellular tRNA slightly, thereby preventing excessive translational costs. In\nboth panels, the x axis represents D, a measure of the disparity between a\ngene’s codon usage and the cellular tRNA supply (see Methods), and the y\naxis represents cellular fitness. The black curves represent the x-y relationship\nassumed in the two models, while the red dot, the gray arrow, and the\nhorizontal dashed line represent the optimal CUB, the direction of\ntranslational selection, and the genetic drift/drift barrier, respectively.\nComment 7\nLine 72-76: The authors’ definition of the “conventional model” is not supported. This\nstatement is a red herring giving the authors the opportunity to disprove a hypothesis\nthat is already quite ridiculous. In no way is the prevailing/conventional wisdom that\ngenes drift toward high CAI. Indeed, if that were the case, there would be no codon\nusage bias among different genes within an organism. The drift toward or away from\nmatching tRNA supply would be highly dependent on 1) how important the protein is\nto the cell’s survival, 2) what concentration of the protein is required for cell function\nand 3) under what conditions is the protein essential. To state that drift toward a low\nmismatch between tRNA supply and demand is ridiculous. The entire premise of this\nmanuscript is based on a red herring.\nResponse:\nThank you for your comment! We do not say or imply genetic drift toward high CAI,\nor toward a low mismatch between tRNA supply and demand in line 72-76, or anywhere\nelse in our manuscript. Both our proposed model and the conventional model consider\ngenetic drift to be a counteracting force to selection (therefore constituting a “barrier”\nto selection). Our model differs from the conventional model in terms of the direction\nof selection, not direction of genetic drift. Our study are always concerned with\ndetermining what is the optimal D (= 0 or not). Genetic drift (and possibly mutation\nbias too) is then always viewed as a “barrier” that prevents selection from reaching the\noptimal D. Please see our Figure 1 and Figure 6a/b for a depiction of this logic.\nComment 8\nLine 73-76: There is unlikely to be a single answer for this since codon bias toward\nhigh and low D exists of different endogenous genes. Again, this depends on the cellular\nneed for the specific function under the conditions tested. There is not likely to be a\nsingle answer. This general premise is highly flawed. The selective pressure is different\nfor different genes.\nResponse:\nThank you ! As we have explained above, the conventional model posits that selective\nmaximizes translational output, and drift/mutation prevents that from occurring,\nresulting in a different CUB (but not optimal CUB) for every gene. Our argument is\nthat the direction of selection is not to maximize translational output, but instead to\nbalance the functional payoff and cost of translation, resulting in different\nbalancing/optimal CUB depending on the functional importance of the protein.\nComment 9\nLine 109-117: Authors are not controlling for changes in 5' region that are translation\ninitiation effects. This confounds results. It's not clear whether the effects are due to\ntranslation initiation (controlled by mRNA structure & ribosome binding) or elongation\nas the authors want to claim. Currently, the authors do not prove that the effects are not\nfrom transcription which can misrepresent the translation results. The authors need to\nshow that the results of expression level changes cannot be explained by translation\ninitiation changes.\nSimilarly, the authors do not control for transcription rates. In bacteria, transcription\nand translation are coupled. The authors should show that their data is not explained\nsimply by changes in mRNA levels of their reporter if they want to claim that it is all\ntranslation elongation effects.\nResponse:\nThank you for this critical point. As for the controlling of translational initiation effects,\nwe would like to clarify that in designing the synonymous mutants for both GmR and\nAmpR, we deliberately left the first 19 codons of the 5′ region unmodified (details in\nTable S2 and S4). This region of “translational ramp” has previously been found to\nregulate translational initiation (PMID: 20403328). Below is a multiple sequence\nalignment of this region in the designed synonmyous GmR variants. The sequence\nname to the left indicates the D value and the index of biological repeat of the variant\n(three different sequences for each D value, but some of them excluded because they\nfailed construction).\nSimilarly, in our additional experiment with AmpR, the multiple sequence alignment of\nthis region in the designed synonmyous AmpR variants appears below.\nRegarding the possibility that our results might be confounded by transcription\nchanges, we would like to separately discussed the issue for GmR and for YFP. As for\nYFP, we measured YFP mRNA expression using quantitative RT-PCR. We found no\ncorrelation between YFP transcript abundance and the D value of GmR, indicating that\ntranscription rate changes cannot explain the observed translational cost. The related\nresults have been updated in our manuscript, as shown in our response to your Comment\n3.\nFor GmR, however, we did NOT argue for the complete exclusion of transcription\nlevel and initiation differences. This is because, for example, codon usage bias can\ndirectly affect the abundance of transcripts independent of promoter activity\n(PMID:28961875 and 27671647). It is also possible that codon usage downstream of\nthe translational ramp resulted in a ribosome “traffic jam” that ultimately decreased\neffective initiation rate (PMID: 25148683). Consequently, the increased cost associated\nwith GmRs with minimum D can be explained by (i) accelerated elongation and\ntherefore tRNA consumption, (ii) increased effective translational initiation and\ntherefore tRNA consumption per mRNA molecule, (iii) increased mRNA half-lives and\ntherefore total tRNA consumption accumulated across the entire life cycle of the mRNA,\nor a combination of these factors. We are, however, arguing that regardless of the\nspecific molecular mechanisms, the cost of increased translational output of GmR due\nto reduced D is readily observable by the reduction of YFP, and also the fitness.\nWe have cited the above related papers and explicitly mentioned these potential molecular\nmechanisms in line 65-68, as pasted below.\nAn exogenous gene with excessively small D will have a very high protein\nexpression level due to an increased abundance of mRNA21,22 and/or more\nefficient translation20,23-25, but at the same time, will impose a heavy\ntranslational cost on the cell, thus hindering the translation of other genes\nwithin the same cell16,20.\nComment 10\nLine 192-195: Authors do not show the data for the KanR gene replacement\nexperiments that is necessary for proving their claim. Needs more information on how\nfunctional importance was quantified. Furthermore, Spearman’s correlation of -0.08 is\nextremely weak almost zero. Authors need to show the data spread used for their\ncorrelations. Just providing a correlation is unacceptable.\nResponse:\nThank you for this comment. The experimental data used in this section were taken\nfrom previously published studies. Therein, the KanR gene was used to replace various\ntarget genes, and the fitness of cells after each gene’s deletion was measured. The\nimportance of each gene to the cell is calculated as 1 minus the fitness of the cell with\nthe focal gene deleted. The scatter plot of all genes is shown below.\nWe fully acknowledge that the Spearman correlation between gene importance and D\nis weak across all genes. This likely reflects the fact that lowly expressed genes (i.e.\nmajority of the genes) are under weaker selection, so that their observed D cannot\naccurately approximate their optimal D (our model predicts a correlation between\nimportance and optimal D, but NOT between importance and observed D). In contrast,\nhighly expressed genes are subjected to stronger selection, which drives their observed\nD closer to the optimal D. Therefore, a correlation between importance and observed D\nis expected only for highly expressed genes. In other words, in line 192-195, the more\nimportant pattern is the strengthening of the expected correlation when only highly\nexpressed genes are included. We have explicitly explained this logic in line 199-211,\nas pasted below.\n...our model predicts that a more important gene should have an optimal\nCUB that is more closely aligned with cellular tRNA supply, and therefore\ncorrespond to a smaller D. To gauge the optimal CUB for the native genes, we\nassumed that the native codon usages of endogenous genes have been\noptimized by translational selection to the maximum extent possible in the\npresence of genetic drift. Here, the relative strength of translational selection\nand genetic drift varies for different genes within the genome. By following\nprevious practice27, we presumed that genetic drift (determined by effective\npopulation size) was constant across genes, and that translational selection\nwas determined by mRNA expression levels. Therefore, our hypothesis predicts a\nnegative correlation between a gene's functional importance and the D value\nof its native codon usage, especially for highly expressed genes because they\nare under stronger translational selection that keeps their native codon usage\nclose to their optimal CUB.\nComment 11\nLine 207-213: This statement is not novel, is intuitive and expected, and is not a\ncontribution of this work. It is already well-accepted that high transcription and\ntranslation rates are highly burdensome to cells. See PMID: 27989436.\nResponse:\nThank you for this comment! We agree that the existence of translational cost is not\nnovel. The novelty here is that the optimal CUB does not maximize translational\nefficiency (rejects the conventional model), but rather balances the functional payoff\nand translational cost (supports the novel model). In other words, the conventional\nmodel on CUB evolution did not take translational cost into account. Please also note\nthat the senior author of paper you mentioned (PMID: 27989436), Dr. Yitzhak Pilpel,\npraised the novelty of our proposed model in comparison to the conventional model in\nhis signed referee report on our manuscript. His original comments on the novelty are\npasted below.\nThe paper titled “A slight mismatch between a gene’s codon usage and the\ncellular tRNA supply is beneficial” describes a very thorough work with a\nthought-provoking notion. Unlike the current dogma in the field, which\nassumes that optimal translation efficiency is achieved when codon usage\nperfectly matches the tRNA pool of the cell, a slight mismatch might be\noptimal. A mechanism based on feedbacl on the tRNA pool is suggested. A\nproposed mechanism is based on an alleged trans-regulation and tRNA\ndepletion.\nWhile previously there were some indications for this effect in unneeded\ngenes, now they look for genes with fitness effect, either a drug resistance gene\nor natural genes. This is a major step forward as it allows them to balance costs\nagainst functional payoff.\nThis is an excellent work. It is thought-provoking, very well done and deep\nin it evolutionary approach. It is also very well written and technically sound.\nComment 12\nLine 217-219: What is the threshold authors set for expression and why? The authors\nneed to show the expression data distribution used for obtaining their correlational data\nin Figure 5.\nResponse:\nThank you for your comment. The logic behind expression thresholding is that natural\nselection (towards optimal D) should be stronger for highly expressed genes. As a result,\nonly the highly expressed genes should show a D close to optimal D, but not lowly\nexpressed genes. Nevertheless, no empirical information is available to determine the\nprecise expression that will produce an observed D which is close enough (for\ndistinguishing conventional from novel models) to the optimal D. Following this logic,\nthese section and Figure 5b-e test a major prediction of our model, namely a correlation\nbetween optimal D and expression, which is only detectable for highly expressed genes\nsince their observed D is closer to optimal D. So we are simply raising the expression\nthreshold before the number of remaining genes becomes insufficient for statistical\nanalysis.\nPer your suggestion, we have now shown some scatter plots (one from each species)\nunderlying the correlations, and compared the observed correlations for highly\nexpressed genes with random gene sets sampled to the same size as the thresholded\nhighly expressed genes. The updated Figure 5b-e are pasted below\nFigure 5. The effect of translational cost on native codon usage of E. coli\nendogenous genes\n(b-e) For E. coli (b), S. cerevisiae (c), D. melanogaster (d) and H. sapiens (e),\nthe Spearman’s Rank Correlation between a gene’s mRNA expression level\nand D of its native CUB (y axis) is estimated among endogenous genes whose\nexpression levels exceeds a given threshold (x axis). Purple bars indicate\nresults from the real set of genes above the expression threshold, while gray\nbars represent results from a randomly sampled set of genes consisting of the\nsame number of genes as the real set. Error bars indicate standard deviations\nof the correlation assessed from 1,000 bootstraps (purple) or random samples\n(gray). An inset scatter plot represents one correlation from each panel, with\nhistogram to the right indicating the distribution of gene expression. Statistical\nsignificance of Wilcox signed-rank test between real and random gene sets is\nindicated as *: P<0.1, **: P < 0.05, or ***: P < 0.01.\nComment 13\nLine 243 and 282-284: E. coli is a gut microbe. This means that evolutionary pressure\nfor optimizing its genes has largely been driven by the gut microenvironment to match\nits demands. The authors culture their E. coli in LB media. This means that there will\nlikely be a mismatch in the mutational accumulation experiments relevant to D-\nincreasing or D-decreasing effects. It is incorrect to assume that the selective pressure\nin lab conditions in LB reflects true selective pressure for E. coli.\nResponse:\nWe sincerely appreciate your comment. The rationale behind the mutation\naccumulation (MA) experiment is that there is a general absence of natural selection in\nMA lines. If the ancestral wildtype strain has been subjected to natural selection towards\nsome particular phenotype, the absence of selection (in MA) could reveal the forces\n(genetic drift and/or mutation bias) that are counteracting natural selection. By\nanalyzing the direction and effect of such counteracting forces as revealed in the MA\nlines, we can identify the type of selection that has been acting on the ancestral wildtype\nstrain. As an analogy, when a balloon's inlet is released (absence of selection), if air\nrapidly exits from inside the balloon (direction of genetic drift/mutation bias), we can\ninfer that the balloon was previously filled with compressed air (selection applied to\nancestral wildtype strain).\nMore specifically in our study, we are not assuming any selective pressure under\nlab conditions in LB. Instead, we assume a lack of selection in MA lines grown in LB,\nso that we can determine what type of selection the wildtype E. coli experienced in their\nnormal (gut) environment based on the contrasting patterns observed in MA lines\n(without selection) and wildtype E.coli cells (with selection). A similar approach has\nbeen used to determine the direction of selection for mutation rates (PMID:34193872).\nComment 14\nLine 303: Authors need to clarify what they mean by “waste of tRNAs for unpreferred codons”.\nResponse:\nWe are sorry for the confusion. Here the “waste of tRNA for unpreferred codons” means\nthat for a gene with CAI=1, all codons are preferred codons, so that the tRNAs for\nunpreferred codons are never used, therefore a waste. We have modified this sentence\nas below.\nFor example, the optimal CAI (=1), in which preferred codons are\nmaximally used, cannot in theory be considered optimal from a supply-\ndemand perspective, since the supply of tRNA for unpreferred codons is never\nused.\nComment 15\nRegarding Figure 3.b\n- Why does at high Gm concentrations the small group (theoretically with less GmR\nprotein than the minimum group) exhibit a higher functional payoff?\n- Why is everything normalized to the large group? This implies you are assuming\nunidirectionality in codon-match and functional payoff.\n- The data contradicts the authors’ assumption that as the importance of a gene increases,\nD (tRNA-CUB supply-demand differential) decreases. At concentrations above 140\nug/ml we see that the small group is preferred. (In text line 154 you say that only 68.9%\nof the variations are explained by this assumption).\nResponse:\nThank you for these three questions.\nFor 1, the payoff difference between the small group and the minimum group is\nNOT statistically significant at high concentration (>140 ug/mL). The observed minor\ndifferences are likely attributable to experimental noise rather than biologically\nmeaningful effects. Also, even in media with high gentamicin concentrations, GmR\nsynthesis cannot exceed a physical limit (e.g., as determined by cell resource uptake\nrate). Functional payoffs and fitness of the group with the minimum D may have been\naffected by this limit, thus making the difference with the group of minimum/small D\nless statistically significant.\nFor 2, we are confused by your statement that we were “assuming unidirectionality\nin codon-match and functional payoff”. We did not make that assumption. More\nimportantly, since the normalization is only for different groups under the same\ngentamicin concentration, the results of the comparison between the three groups under\neach concentration will remain qualitatively unchanged regardless of which group we\nused as the standard for normalization. For example, the large group will be the fittest\nof the three groups at 0 μg/ml, while the small group will be the fittest at 120 μg/ml,\nregardless of the way these measurements are normalized under the same concentration.\nFor 3, we respectfully disagree. The small group is not preferred at 180 μg/ml\ngentamicin, as it is not significantly fitter than the minimum group (Figure 3f). The\n68.9% variation explained refers to the rank variation in functional payoff of the\nminimum group explainable by gene importance (Figure 3c), not variation of D\nexplained by gene importance. Please also see our response to the next comment 16.\nComment 16\nFigure 3.c-e\n- It is unclear what the data points here represent. Are they average or individual single\nrecodes? Why is there only 12 points here? Shouldn’t you have the minimum and small\ngroups here (38 points)? Authors need to provide information about their statistical test\nin the figure legend.\nResponse:\nThank you and we are sorry for the confusion. Each data point in panel c (d) represents\nthe functional payoff (translational cost) of the minimum group as compared to the large\ngroup, at a specific drug concentration. There are therefore only 12 points, one for each\nconcentration. We have revised this figure to clarify this and the statistical test used,\nwhich is pasted below.\nIt should be noted that if the same analysis were performed using the small group\ninstead of the minimum group, the results remained qualitatively unchanged, as pasted\nbelow. We decided to only present the results of the minimum group.\nComment 17\nFigure 3.f\n- Relevant to result interpretations of line 160-163:\n- I disagree that at concentrations above 120 ng/ul minimum D is selected. Mainly\nbecause the statistic is weakly probable. From 140 ng/ul all groups most likely have the\nsame relative fitness. This contradicts your statement that minimum D is best at the\nhighest functional importance (highest concentrations).\nResponse:\nThank you for this question. The optimality of minimum D is also supported by two\nadditional results. In Supplementary Figure S1, where at concentration 160 and 180\nμg/ml, fitness is significantly negatively correlated with D. Thus, in the presence of\nhigh concentrations of gentamicin, enhancing GmR expression by minimizing D is\nmost likely to be beneficial. Also, in response to your Comment 1, we have conducted\nadditional experiments using AmpR, the ampicillin resistance gene. Consistent with our\nmodel, we observed that at high ampicillin concentration (4000 μg/ml), minimum D is\nthe fittest group. More details on this is given in response to your Comment 1. And this\nresult has been added to the new Supplementary Figure S5.\nWe have chosen to present the grouped result in the main figure and the correlation\nin the supplementary. The rationale for this is twofold. First, it seems obvious that,\nat high gentamicin concentrations, GmRs with a smaller D and a higher expression level\nshould be fitter than those with a larger D and a lower expression level, unless the\nexpression is too high for the physical limit defined by the resource uptake rate of the\ncell. Second, and more importantly, the most critical pattern for our argument is the\nsmall mismatch group being fitter than both the minium mismatch group and the large\nmismatch group at intermediate drug concentrations. The current main figure with three\ngroups illustrates this most clearly.\nComment 18\nFigure 4\n- This correlation does not support any direct like to the hypothesis that genes with\nsmaller D scores have a higher functional payoff. Because the function is related to the\ndirect role of the gene and not how it is coded. In order to support your hypothesis, you\nneed to have direct fitness data where you measure the fitness cost associated with\nrecoding a specific vital gene and replacing it with the native gene to determine the link\nbetween functional payoff and translation cost.\n- With regards to line 211-213 An important gene does not have to be highly expressed.\nYou cannot establish the importance of genes by only looking at highly expressed genes\nand measuring expression indirectly by the D index.\nResponse:\nThank you for giving us the opputunity to clarify this.\nFor 1, we are not hypothesizing that genes with smaller D have a higher functional\npayoff (although we believe this is likely to be the case, since small D gives higher\nexpression and, therefore, higher functional payoff). Instead, our hypothesise is that a\nmore important gene should have a D value that is more closely aligned with the cellular\ntRNA supply. There are two critical differences between your comment and our\nhypothesis. First, a gene’s importance is determined by its functional payoff minus its\ntranslational cost, so a more important gene is not necessarily one with a higher\nfunctional payoff, but rather one with a lower translational cost. Second, our hypothesis\nposits functional importance as the cause and optimal D as the consequence, rather than\nthe other way around.\nWe agree your suggested experiment could be useful, and your reasoning behind\nthis design coincides quite well with our manipulative experiments using recoded\nantibiotic resistance genes. However, because similar experiments have not been\nconducted for every genes in a genome yet, we turn instead to an indirect test, which is\ndescribed in the Result section, as pasted below.\n… our model predicts that a more important gene should have an optimal\nCUB that is more closely aligned with cellular tRNA supply, and therefore\ncorrespond to a smaller D. To gauge the optimal CUB for the native genes, we\nassumed that the native codon usages of endogenous genes have been\noptimized by translational selection to the maximum extent possible in the\npresence of genetic drift. Here, the relative strength of translational selection\nand genetic drift varies for different genes within the genome. By following\nprevious practice27, we presumed that genetic drift (determined by effective\npopulation size) was constant across genes, and that translational selection\nwas determined by mRNA expression levels. Therefore, our hypothesis\npredicts a negative correlation between a gene's functional importance and the\nD value of its native codon usage, especially for highly expressed genes\nbecause they are under stronger translational selection that keeps their native codon usage\nclose to their optimal CUB.\nWe hope this explanation helps clarify the logic behind Figure 4.\nFor 2, as you point out, and we agree, an important gene does not necessarily have\nto be highly expressed. We are not measuring expression by D, nor are we suggesting\nthat highly expressed genes are always important. Instead, we are stating that highly\nexpressed genes are more sensitive to reductions in D. This is because they, by\ndefinition, have more transcripts. A reductions in D in a highly expressed gene will\ntherefore have a greater impact on the tRNA supply-demand than a same reduction in\nD in a lowly expressed gene. This then predicts a larger optimal D in highly expressed\ngenes than lowly expressed genes of similar importance. Please note here, expression\nis the cause and optimal D is the consequence, not the other way around. We have\nclarified this point as pasted below.\nMore specifically, a D-decreasing mutation in a gene with high transcript\nabundance should result in a greater increase of translational cost than that in\na gene with low transcript abundance, since a greater number of RNA\nmolecules are mutated in the former case. Consequently, highly expressed\ngenes should be more sensitive to reductions in D, resulting in a larger optimal\nD in highly expressed genes than lowly expressed genes of similar importance\n(Figure 5a).\nPlease also see our response to the next Comment 19.\nComment 19\nFigure 5\n- In regards to line 216: Because D and expression level are anticorrelated does not\nimply anything about translational cost. You are only looking at functional payoff.\nTranslation cost can only be measured experimentally.\nResponse:\nThank you. To illustrate our logic, consider two genes with identical D values, but\nidentical transcriptional levels: Gene 1 has 1000 mRNA molecules, while Gene 2 has\nonly 1. If a mutation that reduces D to the same level in both genes, the extra\ntranslational resources required by Gene 1 would be 1,000 times greater than those\nrequired by Gene 2. This example illustrates that, genes with higher mRNA abundance\nare more sensitive to reductions in D, therefore they tend to have a larger optimal D.\nIn contrast, genes with lower expression levels should be less sensitive to reduction in D,\nthus they tend to have a smaller optimal D values (although their observed D are usually\nnot sufficiently close to their optimal D sine translation selection is weak on these\ngenes).\nWe would like to emphasize that this prediction of a positive correlation between\ngene expression and optimal D among highly expressed genes is unique to our novel\nhypothesis. If only functional payoff is considered, such as in the conventional model,\nthe optimal D should be negative correlated with gene expression (because selection\nminimize D to enhance translational efficiency/accuracy in highly expressed genes).\nWe have also illustrated this logic in Figure 5a.\nComment 20\nQuestion on the methods:\n- How did the authors control for the sample size reduction of highly expressed genes?\nCorrelation scores tend to shrink in absolute values as sample sizes grow. Conversely,\na strong correlation at a small sample size may not be true when the sample of the\ncorrelation is increased.\nResponse:\nThank you for this important methodological question. To address potential bias\nassociated with varying sample sizes, particularly the reduction in the number of genes\nwhen the analyses is limited to highly expressed genes, we randomly down-sampled\nthe genes to the same size as those retained at each expression threshold. Through this\napproach, we are able toassess whether the observed correlations remain consistent\nacross comparable sample sizes, thus ensuring that differences in group size do not\nsignificantly affect the strength and significance of the correlations. This additional\nanalysis further supports our conclusions in both Figure 4 and Figure 5, where\nsubsetting by expression threshold is conducted. The revised Figure 4 and Figure 5 is\npasted below.\nFigure 4. The effect of functional payoff on native codon usage of E. coli\nendogenous genes\n(a-b) The Spearman’s Rank Correlation between a gene’s functional\nimportance and D of its native CUB (y axis) is estimated for E. coli (a) or S.\ncerevisiae (b) endogenous genes whose expression levels exceed a specific\nthreshold (x axis). Error bars indicate standard deviations of the correlation\nassessed from 1,000 bootstraps (purple) or random samples (gray). An inset\nscatter plot represents one correlation from each panel. Statistical significance\nof Wilcox signed-rank test between real and random gene sets is indicated as\n*: <0.1, **: P < 0.05, or ***: P < 0.01.\nFigure 5. The effect of translational cost on native codon usage of E. coli\nendogenous genes\n(a) A diagram illustrating why the optimal D would be higher for genes with\ngreater translational costs. (b-e) For E. coli (b), S. cerevisiae (c), D.\nmelanogaster (d) and H. sapiens (e), the Spearman’s Rank Correlation\nbetween a gene’s mRNA expression level and D of its native CUB (y axis) is\nestimated among endogenous genes whose expression levels exceeds a given\nthreshold (x axis). Purple bars indicate results from the real set of genes above\nthe expression threshold, while gray bars represent results from a randomly\nsampled set of genes consisting of the same number of genes as the real set.\nError bars indicate standard deviations of the correlation assessed from 1,000\nbootstraps (purple) or random samples (gray). Statistical significance of\nWilcox signed-rank test between real and random gene sets is indicated as *:\nP<0.1, **: P < 0.05, or ***: P < 0.01.\nComment 21\nFigure 6\n- While this evidence does support the authors' claim that a slight mismatch with the\ntRNA supply is selectively maintained, the authors do not quantify the strength of the\nD-changing mutations. In other words, what is the distribution of D score change in\neach D-increasing/-decreasing mutation group?\nResponse:\nThank you for this question. As for Figure 6c-h, the effect size (absolute value of D\nscore changes) is analyzed in d/e/g/h. For results in Figure 6i and j, we have now added\na Supplementary Figure S8 to show the distribution of D changes for each mutation\ngroup, which is pasted below. Since the effect size of D-decreasing and D-increasing\nmutations is similar, the observation in Figure 6i and j is not confounded by the effect\nsize difference between the two types of mutations.\nFigure S8. The distribution of D changes in D-decreasing and D-\nincreasing mutations\nViolin plots showing the changes in D due to the mutations listed in Figure\n6j. The effect size of D-decreasing and D-increasing mutations are similar\n(Wilcox rank sum test P = 0.10).\nComment 22\nAnother note regarding the authors' quantification of cellular tRNA:\nIs the index for quantifying the cellular tRNA pool linear in its entire range? Perhaps if\nthis index becomes non-linear in either boundary it may falsely inform that D\ndecreasing mutations are deleterious at highly expressed genes.\nResponse:\nThank you for this comment. We did not claim or assume that D linearly represents the\ntRNA supply or the tRNA supply-demand mismatch. We did not make conclusions\nbased solely on the linearity assumption, such as Pearson's correlation tests with D (we\nalways accompany them with Spearman's rank-based tests). The reasons for such non-\nlinearity can be biological (the degree of mismatch is not linear to the phenotypic\nconsequences) and/or technical (different amino acids might have different impacts, but\nwe are merely using their geometric average).\nAlso, we have explained (in response to your Comment 2) the rationale behind\nestimating tRNA supply using transcriptomic codon usage of highly expressed genes,\nand our conclusion remained unchanged if high-throughput tRNA sequencing data is\nused to approximate tRNA supply instead.\nMore specifically for this comment (presumably you are referring to Figure 6i),\nwe repeated it with D calculated using tRNA gene copy number in wildtype yeast (from\nGtRNAdb). As shown in the figure pasted below, the observation that an\noverwhelmingly large percentage of both D-decreasing and D-increasing mutations are\ndeleterious remains unchanged. As a result, we exclude the possibility that the non-\nlinearity of D might be responsible for this observation.\nRevewer #2:\nOverall Comment:\nThe manuscript submitted by Chen, Liu, Zhou and co-authors is dedicated to the study\nof translational selection. Specifically, the authors propose a new model to explain the\nevolution of codon usage bias (abbreviated CUB; this can be roughly defined as the\ndeparture from a random distribution among synonymous codons), according to which\nthe optimal codon usage does not (always) perfectly match the abundance of tRNAs in\nthe cell. According to this model, a mismatch between the codon usage bias of a given\ngene and the tRNA supply can be beneficial because it reduces the cost imposed by the\ntranslation of the gene. In other words, a highly expressed and highly translated gene\nwould be costly because its translation would hamper the translation of other genes.\nTo support this model, the authors mainly rely on genome editing experiments\nperformed in E. coli. Specifically, they insert cassettes carrying a gentamicin-resistance\ngene (GmR) with different codon usage biases. The GmR gene is fused to a mCherry\ncoding sequence, which enables the authors to measure the expression level of the\nGmR-mCherry fusion through fluorescence measurements. The cassettes also contain\na YFP gene in the opposite orientation. The expression level of the YFP gene, again\nestimated through fluorescence measurements, enables the authors to estimate the\ntranslational cost of the CUB of the GmR gene as the reduction in YFP protein levels.\nThe authors perform these experiments at different concentrations of gentamicin and\nthey also measure the growth rate of the different E. coli clones, thus obtaining a more\ndirect estimation of the fitness effect of the different constructs. In addition to these\nmain experiments, the authors also analyze previously published data for mutation\naccumulation or gene replacement experiments in E. coli and yeast, as well as codon\nusage data for various species.\nI find the model proposed by the authors to be interesting and reasonable. It brings an\nalternative to the classical model which posits (even in species where translational\nselection is efficient, that is, species with high effective population size such as E. coli,\nyeast or fruitfly) translational selection cannot optimize the codon usage bias to\nperfectly match the tRNA abundance of the cell because of genetic drift. However, it\nmust be stated that the classical model already incorporates the idea of a difference\nbetween gene categories (the translational selection coefficient is expectely to vary\namong genes, depending on the optimum protein quantity produced by each gene, and\nfor some genes the translational selection coefficient is indeed expected to be very low).\nMore importantly, I am concerned about the validity of several assumptions made in\nthis article and of several analyses performed in this article.\nResponse:\nWe are truly grateful for your positive assessment of our work. In accordance with your\nsuggestions, we have conducted additional analyses and incorporated new data to\naddress the points raised. Please find our point-to-point responses below.\nComment 1\nFirst, the authors state that the D metric that they define is a measure of the mismatch\nbetween the CUB of a gene and the tRNA supply of the cell. It is in fact a measure of\nthe mismatch between the CUB of a given gene and the CUB of the highly expressed\ngenes (Methods, lines 358-373). The authors assume that the CUB of highly expressed\ngenes has been optimized to match the tRNA supply, although they then affirm that the\noptimal CUB should not perfectly match the tRNA supply, even for highly expressed\ngenes. This is circular and illogical. There are better ways to approximate the tRNAs\nsupply without assuming that the CUB for highly expressed genes has been optimize to\nmatch it. If the authors are unable to measure directly tRNA abundance, they could use\nan approximation based on tRNA gene copy number, taking into account wobble rules\netc. As this affirmation is the basis of the new model proposed by the authors, it is very\nimportant to make sure that D indeed measures the distance compared to the CUB that\nwould perfectly match the tRNA supply of the cell.\nResponse:\nThank you very much for your suggestion. We would like to address from two angles:\n1. The assumption on CUB of highly expressed genes has been optimized to match the\ntRNA supply\nThe logic behind our estimation of tRNA supply is twofold. First, the actual tRNA\nsupply is neither the copy number of the tRNA gene nor the abundance of tRNA\ntranscripts. Instead, it is the abundance of the ternary complex of aminoacyl-tRNA+EF-\nTu+GTP (or replace EF-Tu with EF1A in eukaryotes). Currently, the abundance of\nternary complexes for different types of tRNAs cannot be systematically measured. In\naddition, the kinetic parameters between different codon-anticodon pairs will also\ninfluence the \"effective\" tRNA supply, further complicating the estimation process.\nSecond, because of the first point, we resorted to the assumption that the cellular tRNA\nsupply and demand (codon usage) of the whole transcriptome have been evolutionarily\noptimized by natural selection to ensure a close match between supply and demand.\nThis assumption allows us to approximate the supply of tRNA by using the codon usage\n(the demand) of the top 10% expressed gene (which dominate the transcriptome) in the\ntranscriptome. It is important to note that we are using the transcriptome, not the\ngenome, which means the codon usage has been weighted by the expression of those\ngenes before being summed up. Furthermore, this assumption of supply-demand match\nhas been empirically supported (PMID: 22479199 and PMCID: PMC12262485), and\nhas been used in our previous publications (PMID: 32123323 and 34977494). It is also\nimportant to note that a supply-demand match on a transcriptome level does not\nnecessarily imply a supply-demand match for individual genes. Therefore, this\nassumption is neither circular nor illogical.\n2. Using measurements of tRNA transcript abundance for tRNA supply does not change\nour conclusion.\nPer your suggestion, we obtained tRNA expression levels based on high-\nthroughput tRNA sequencing data from a previous publication (PMID: 31353208, GEO\naccession GSE128812), which were then scaled by codon-anticodon affinity and used\nas the tRNA supply. This metric is similar to the “adaptiveness value” (w) of a codon,\nas used in the popular tAI metric (PMID:15448185). It was found that the grouping\n(minimum, small, or large D) in Figure 3 remained unchanged, suggesting that our\nconclusion is not affected by the change in the method used to quantify tRNA supply.\nWe have now included this result in Supplementary Figure S2, and it is pasted below.\nFigure S2. Grouping by D of GmR strains remains unchanged if high-\nthroughput tRNA sequencing data are used to estimate tRNA supply\nThe D values of designed GmR synonymous variants were calculated with\ntRNA supply based on either transcriptomic codon usage of the top 10%\nhighly expressed genes (x axis) or tRNA transcript abundance derived from\nhigh-throughput tRNA sequencing (y axis). The black lines indicate the\ngrouping, which is identical for both types of D values. Please note that each\ndot has up to three unique sequences as biological replicates (see Table S2).\nComment 2\nSecond, the authors measure the translational cost of the CUB of the GmR artificial\nconstructs by measuring the fluorescence level of the exogenous YFP present in the\nconstruct. However, they do not ascertain that variations in YFP fluorescence levels are\nindeed to variations in YFP translation, rather than for example in variations in YFP\ngene transcription. They also do not show that the variations in YFP fluorescence are\ncausally connected to the translation of the GmR-mCherry gene. The authors and others\nhave shown before that variations in CUB can also affect the mRNA level of genes,\nmainly due to mRNA stability. Additional effects ca be imagined at the transcriptional\nlevel, given the overall differences in base composition of the CUB variants. As the\nYFP gene is inserted in close proximity to the GmR-mCherry gene, on the reverse\nstrand, transcriptional interference can occur between the two genes. Thus, their claim\nof a translational cost is not sufficiently supported. Furthermore, the translational cost\nshould be evaluated on endogenous genes. As this is an important claim of the\nmanuscript, it should be properly verified, by evaluating mRNA levels and translation\nrates in the different constructs (for example with RNA sequencing and riboosome\nprofiling assays).\nResponse:\nThank you very much for this insightful comment. Please allow us to address these\npoints separately.\nTo prevent transcriptional interference between YFP and GmR, we have included\nan insulator between the terminators of YFP and GmR to prevent transcriptional\ninterference between the two genes. The Figure 2b pasted below has been updated to\nexplicitly indicate the insulator.\nFurthermore, we used quantitative RT-PCR to measure YFP mRNA expression levels\nin the GmR strains growing in LB. We found that YFP mRNA expression did not\ncorrelate with GmR’s codon usage (D value). This result has been added to\nSupplementary Figure S9c and pasted below.\nFigure S9. Correlations between D and functional payoff, translational\ncost or fitness\n(c) The abundance of YFP mRNA in different strains, as determined by\nquantitative RT-PCR (see Methods), is not correlated with codon usage of\nGmR. The Spearman’s Correlation Coefficient ρ of the raw data, as well as\nthe corresponding P value, are indicated.\nTo evaluate the effect of GmR’s translational cost on endogenous genes, we\nperformed ribosome profiling (Ribo-seq) on four strains with varying D values (0.127,\n0.278, 0.778, and 0.878) of the GmR gene. It is observed that when the GmR gene has\na minimum D (D = 0.127, 0.278), the codons frequently used by GmR tend to have a\nlonger decoding time, leading to a significant positive correlation between decoding\ntime and GmR codon consumption among genes with sufficient Ribo-seq coverage.\nThis observation suggests that tRNA depletion by codons frequently used in GmR is\nslowing down the elongation of the same codons in other endogenous genes. In contrast,\nstrains expressing GmRs of large D did not show a similar correlation. This result has\nnow been added to the Result section and as Supplementary Figure S3, which is\npasted below. Please note that RNA-seq is not required here as the relative decoding\ntime is always normalized within the same gene first, and regardless of the RNA-seq\nresults, all codons within the same gene will have the same transcription level.\nFigure S3. Slowed decoding as a result of GmR translation\nThe decoding time of each codon in endogenous E. coli genes was measured\nby Ribo-Seq in wildtype and four E. coli strains expressing GmRs of different\nD (see Methods). The Spearman’s correlation coefficient between GmR’s\nconsumption of a codon and its relative decoding time as compared to the\nsame codon in wildtype E. coli (y axis) is calculated for endogenous genes\nabove a given minimum value of translational activity (x axis) (see Methods).\nTogether, these results suggest that genes with minimum D impose a greater burden on\nthe translation of other endogenous genes, whereas those with large D contribute\nminimally to the overall translational burden.\nComment 3\nThird, the fact that the different E. coli strains that they used could display variations in\nplasmid copy number is quite concerning. I am not convinced that the lack of\ncorrelation between the copy number and D is sufficient to ensure that the results are\nnot biased. The authors should verify that the plasmid copy number variation is not a\nconfounding factor for fitness analysis and for translational cost analyses.\nResponse:\nThank you for your comment. According to your suggestion, we calculated the partial\nSpearman's correlation between the copy number and functional payoff/translational\ncost/relative fitness while controlling for the plasmid copy number. The results shown\nin Figure S9b, and pasted below, are nearly identical to those in Figure S1, which are\nthe same correlations without controlling for the plasmid copy number. This suggests\nthat plasmid copy number variation is not a confounding factor in these analyses.\nFigure S9. Copy number of GmR genes identified by quantitative real-time PCR.\n(a) The dots and error bars represent the average and standard error of the copy number\nof GmR genes relative to dxs as measured by qPCR. The strains growing in LB medium\nsupplemented with 0 μg/ml (left panel) or 60 μg/ml (right panel) gentamicin is analyzed.\nThe Spearman’s Rank Correlation Coefficient (ρ) and corresponding P value are shown\nat the top left corner. (b) The partial Spearman’s correlation between D and functional\npayoff/translational cost/fitness controlling for the measured plasmid copy number.\nPlease also note that, we did NOT argue for the complete exclusion of translational\ncost contributed by increased plasmid copy number or increased GmR mRNA level.\nThis is because, for example, codon usage bias can directly affect the abundance of\ntranscripts independent of promoter activity (PMID:28961875 and 27671647). It is also\npossible that in high gentamicin concentration, strains with a high plasmid copy number\nexhibit growth advantage for its high baseline GmR expression. Essentially, the\nincreased cost of GmR can be mediated by increased DNA abundance, mRNA\nabundance, and translational activity per mRNA. As long as these, or other potential\nmechanisms, simultaneously mediated the increased functional payoff AND\ntranslational cost, they are compatible with our model. We are, however, arguing that\nregardless of the specific molecular mechanisms, the cost of increased translational\noutput of GmR due to reduced D is readily observable by the reduction of YFP,\nincreased decoding time of codons frequently used by GmR, and also the fitness.\nComment 4\nFourth, the fact that the authors observe a positive correlation between D and expression\nlevels among the most highly expressed genes (top 10%, top 5% and top 3%) in all 4\nanalyzed species (Figure 4) is puzzling. In human, translational selection is highly\ninefficient, unlike in the other 3 species. The fact that the same effect is observed in\nhuman, E. coli, yeast and Drosophila rather points to a more general effect, potentially\na confounding effect of base composition. This should be examined carefully.\nResponse:\nThank you for this insightful comment. We believe you mean Figure 5, as Figure 4 is\nabout gene importance and D, and only includes two species. We agree that human\nshould have a less efficient translational selection, for it has a relatively small effective\npopulation size. And this is consistent with the observation that for most highly\nexpressed genes, the positive correlation between D and expression levels are weaker\nin human (ρ < 0.2), but stronger in other species (ρ > 0.2).\nTo further exclude the potential confounding effect of base composition, we\nrepeated the observed correlations by partial correlations that controlled A%, C%, G%\nT% or GC% among those genes. As shown in Figure S6, which is pasted below,\ncontrolling for base compositions does not affect our observation of a positive\ncorrelation between gene expression and D among highly expressed genes.\nFigure S6. The correlation between D and expression for endogenous\ngenes are not confounded by base composition bias\nEach column of panels is similar to Figure 5b-e, except that the Spearman’s\ncorrelation in Figure 5b-e is replaced with partial Spearman’s correlation\ncontrolling for the percentage of Adenine (A%, in panels a/f/k/p), Thymine\n(T%, in panels b/g/l/q), Guanine (G%, in panels c/h/m/r), Cytosine (C%, in\npanels d/i/n/s) or G+C% (in panels e/j/o/t) in the sequences.\nComment 5\nFifth, it seems that the mutation accumulation analysis was performed using both\nsynonymous and nonsynonymous mutations (only nonsense mutations seem to have\nbeen excluded, Methods, line 464 to 466). This does not make sense - if one wants to\nmeasure the effect of the D, only synonymous mutations should be considered. It is also\nnot clear why this analysis was only performed on the 10% most highly expressed genes\n(especially given that D measures the distance to the codon usage of the most highly\nexpressed genes). Here, a comparison between gene expression classes would be useful,\nto test the hypothesis that the CUB of the highly expressed genes has indeed been\noptimized. Here again, D should be computed with respect to the codon usage that\nperfectly matches the tRNA abundance (not to the codon usage of highly expressed\ngenes).\nResponse:\nThank you very much for this comment. We now only consider synonymous mutations\nfrom MA lines. The result supports our model, and is pasted below. We have also\nupdated the Method section to reflect this.\nAlso, we would like to explain why this analysis is limited to the 10% most highly\nexpressed genes. The rationale behind the mutation accumulation (MA) experiment is\nthat there is a general absence of natural selection in MA lines. If a gene in the ancestral\nwildtype strain has been subjected to natural selection towards some particular\nphenotype, the absence of selection (in MA) could reveal the forces (genetic drift and/or\nmutation bias) that are counteracting natural selection. By analyzing the direction and\neffect of such counteracting forces as revealed in the MA lines, we can identify the type\nof selection that has been acting on the ancestral wildtype strain. As an analogy, when\na balloon's inlet is released (absence of selection), if air rapidly exits from inside the\nballoon (direction of genetic drift/mutation bias), we can infer that the balloon was\npreviously filled with compressed air (selection applied to ancestral wildtype strain).\nBased on this logic, studying lowly expressed genes would be useless, since their\nancestral sequences have not been subjected to strong selection on D, as suggested by\nresults shown in Figures 4 and 5.\nAs for the question regarding the definition of D, please refer to our response to\nyour Comment 1. We hope you found these explanations satisfactory.\nComment 6\nFinally, the authors might be interested in this preprint, which shows that translational\nselection is surprisingly rare in metazoans:\nhttps://www.biorxiv.org/content/10.1101/2024.07.22.604600v1 .\nResponse:\nThank you very much for letting us know about this interesting paper. We have cited\nthis paper in the related context, which is pasted below.\nWe found that D and expression level are positively correlated among the\nmost highly expressed genes (Figure 5b, purple bars), supporting the effect\nof translational cost due to excessively small D. The positive correlation is not\nan artifact of the limited number of genes above the expression threshold,\nsince randomly downsampling all genes to the same number as those retained\nat each expression threshold will not result in the same correlation (Figure 5b,\ngray bars). Moreover, endogenous genes in yeast, fly, and human show similar\npatterns (Figure 5c-e), and these patterns cannot be attributed to potential\nbase composition biases (Figure S6). Interestingly, the maximum correlation\nachieved is the weakest in human compared to the other three species, which\nis consistent with the weak translational selection observed for metazoan\nspecies33.\nComment 7 (Remarks on code availability\nI haven't yet reviewed the code but I plan to do so for a revised version of the manuscript.\nResponse:\nThank you and we have updated our GitHub repo with the code for the additional\nanalyses during this revision.\nReviewer #3:\nOverall Comment:\nThe paper titled “A slight mismatch between a gene’s codon usage and the cellular\ntRNA supply is beneficial” describes a very thorough work with a thought-provoking\nnotion. Unlike the current dogma in the field, which assumes that optimal translation\nefficiency is achieved when codon usage perfectly matches the tRNA pool of the cell,\na slight mismatch might be optimal. A mechanism based on feedbacl on the tRNA pool\nis suggested. A proposed mechanism is based on an alleged trans-regulation and tRNA\ndepletion.\nWhile previously there were some indications for this effect in unneeded genes, now\nthey look for genes with fitness effect, either a drug resistance gene or natural genes.\nThis is a major step forward as it allows them to balance costs against functional payoff.\nThis is an excellent work. It is thought-provoking, very well done and deep in it\nevolutionary approach. It is also very well written and technically sound. I would\nultimately recommend publication in NC. Yet, if I am contribute to the robustness of\nthe conclusion and to the discussion that this paper will surely evoke, here are some\nimportant points.\nI’m accepting to sign the report with my name: Yitzhak (Tzachi) Pilpel.\nResponse:\nWe are profoundly grateful for your generous and encouraging evaluation of our work,\nespecially from a respected scholar like yourself.\nComment 1\n1. A critical experiment could be to manipulate the tRNA pool and then test, what I\nbelieve is a prediction of their theory. That for different tRNA pools, alternative\nmembers of the collection of 53 synonymous variants would optimal. This is a clear\nprediction of their model since for each tRN pool other variants might be at optimal D.\nIf they showed that I would be convinced better. If the result would show that same\nvariants are good at different tRNA pools then their conclusion is not supported. In\nother words, they have manipulated so far the codons – the demand, but never the tRNA\n– the supply.\nResponse:\nThank you very much for the critical and very constructive comment! While\nmanipulating all tRNA genes is difficult, following your suggestion, we aimed to\noverexpress one tRNA by inserting it into the high copy number PUC57 vector. We\nassumed this would increase the supply of this specific tRNA at least twice as much as\nthe supply of the most abundance synonymous tRNAs in the wild-type strain (I.e., set\nthe 𝑋 of the affected codon 𝑗 to twice the maximum 𝑋 of synonymous codons\n𝑖𝑗 𝑖𝑗\nfor amino acid 𝑖, then rescale all 𝑋 of amino acid 𝑖).\n𝑖𝑗\nWe selected the argWCCT tRNA\ngene and two specific GmR variants (D = 0.328, in the minimum D group, and D =\n0.778, in the large D group), as these variants will be moved into the small D group as\na result of the tRNA overexpression. Via homologous recombination, we inserted the\nargWCCT tRNA gene, driven by its native promoter, into the respective GmR plasmids\n(Figure S4a, as pasted below), resulting in two tRNA-overexpressed strains.\nAs is predicted by the conventional model, a decrease in D caused by altered tRNA\nsupply is always beneficial (or at least not detrimental in low or intermediate gentamicin\nconcentrations). On the contrary, our model predicts that at intermediate gentamicin\nconcentration, both strains should be fitter after the tRNA supply change, whereas\nfitness changes would be opposite at low or high gentamicin concentration (Figure\nS4b). The experimental results support our model (Figure S4c/d), although the fitness\nchanges observed for the strain with increased D are not statistically significant,\npossibly due to the relatively small change in D.\nFigure S4. Fitness change upon manipulating tRNA supply\n(a) A schematic illustration of the expression vector structure with (bottom)\nor without (top) the additional tRNA gene argW. (b) Our model’s prediction\nfor fitness changes following an increase in argW supply. This is essentially\nsame as the bottom panel of Figure 3a. (c and d) The relative fitness of two\nGmR variants in cultures with 10, 120, and 180 μg/ml gentamicin before and\nafter argW tRNA overexpression, which resulted in a D values change from\n0.328 to 0.338 (c) and from 0.778 to 0.632 (d). The statistical significance of\ntwo-tailed Wilcoxon Rank Sum Tests are indicated. The rightmost bar in d\nindicates no observable growth for the native strain under 180 μg/ml\ngentamicin.\nCritically, these results clearly demonstrate a case where changes in tRNA supply will\nalso lead to a change in relative fitness in the direction predict by our model, which\nfurther strengthen our model.\nComment 2\n2. The proposed mechanism that could mediate their effect is that of tRNA depletion,\n“In theory, the same logic applies to endogenous genes native to the cell, since they also\ntrans-regulates the expression of other genes via tRNA depletion” – Ref #23. First, I\ndon’t think reference #23 shows tRNA depletion. This is a misunderstanding of what\nwas done in Frumkin et al. The increase in tRNA supply in ref #23 was very different,\nfirst because it was done artificially by the researchers (not spontaneously by the cell),\nand second, because the tRNA levels themselves were not changed, but rather tRNA\nanticodons were mutated. So it’s incorrect to ascribe to ref#23 a trans-regulation of the\ntRNA pool. That said, I think we are left without a true reference or a proof, that the\ntRNA pool is indeed trans-regulated in response to change in codon demand. So I can\nsee several options: 1. The authors measure tRNA pool with some existing tRNA seq\nprotocols and show the effect; 2. Find a reference that truly shows this effect, or 3.\nAcknowledge that we don’t know if such effect indeed occurs at all\nResponse:\nThank you for this insightful comment. We agree that neither our study nor Frumkin et\nal. directly observed tRNA depletion. Here, we believe that measuring tRNA transcript\nabundance is not useful, since several factors other than the tRNA transcript abundance\ncan affect the “effective” tRNA supply (see our response to the next Comment 3 for a\nmore detailed explanation on this point). Instead, we shifted the focus a little bit\ndownstream, to investigate whether the decoding time of the codons frequently used by\nGmR is indeed slowed down when GmR has a minimum D. This directly tests the\ntranslational burden exerted by the GmR on the level of individual codons, which is\nmost parsimoniously explained by a shortage of the translation-ready tRNA. To this\nend, we performed ribosome profiling (Ribo-seq) on four strains with varying D values\n(0.127, 0.278, 0.778, and 0.878) of the GmR gene. It is observed that when the GmR\ngene has a minimum D (minimum mismatch with the tRNA pool, D = 0.127 or 0.278),\nthe codons frequently used by GmR tend to have a longer decoding time, which gives\nrise to a significant positive correlation between decoding time and GmR codon\nconsumption. On the contrary, strains expressing GmR genes of large D did not exhibit\nthe same correlation. This result has now been added as Supplementary Figure S3, and\npasted below. Note that RNA-seq is not necessary here since the relative decoding time\nis always normalized within the same gene first, and regardless of the RNA-seq results,\nall codons within the same gene will have the same transcription level.\nWe have removed “tRNA depletion” when citing Frumkin et al., and also explicitly\nmentioned that we have not proved the existance of tRNA depletion in Discussion,\nwhich is pasted below. We hope these revision have better elucidate this point.\n… in terms of the molecular mechanism behind our observation, our working\nmodel is tRNA depletion triggered by excessive translation of genes with\nminimum D. Although the results from our experiments with GmR and AmpR,\nas well as those obtained from tRNA overexpression and Ribo-Seq, can\nparsimoniously be explained by this model, it has not been proven\nconclusively. Moreover, we did not rule out other molecular mechnisms that\ncould contribute to GmR’s translational cost, such as increased plasmid copy\nnumber, transcription level, or translational initiation rate. Future research\nwill be necessary to fully elucidate the molecular mechanisms involved.\nComment 3\n3. Their calculation of the D measure – the extent of deviation from the tRNA pool and\ncodon usage is critical, especially since the whole argument is based on the difference\nbetween D that approach zero and D that is slightly above zero. My problem is with\ntheir estimation of the tRNA pool which is needed for equation in line 358. Particularly\nthey say “Xij represents the tRNA supply of the cell, approximated by the fraction of\ncodon j among the synonymous codons of amino acid i for the highly expressed genes\nin the transcriptome of the cell” ref #47. I’m not sure if and how ref #47 provides us with a\ntRNA pool. So a few comments on this: 1. I’d like to see their results reproduced with other, more conventional definitions of\nthe tRNA pool. Again, they could do tRNA seq. There’s also the classical tAI – based tRNA pool deduced from tRNA gene\ncopy number. Is their D optimized to be slightly above zero also when using these more conventional estimations of the\ntRNA pool? The way in which ref #47 appears to be estimating the tRNA pool is affected by codon demand, this is\npotentially problematic.\nResponse:\nThank you very much for this important comment. We would like to address from two\nangles.\n1. The assumption on CUB of highly expressed genes has been optimized to match the\ntRNA supply\nThe logic behind our estimation of tRNA supply is twofold. First, the actual tRNA\nsupply is neither the copy number of the tRNA gene nor the abundance of tRNA\ntranscripts. Instead, it is the abundance of the ternary complex of aminoacyl-tRNA+EF-\nTu+GTP (or replace EF-Tu with EF1A in eukaryotes). Currently, the abundance of\nternary complexes for different types of tRNAs cannot be systematically measured. In\naddition, the kinetic parameters between different codon-anticodon pairs will also\ninfluence the \"effective\" tRNA supply, further complicating the estimation process.\nSecond, because of the first point, we resorted to the assumption that the cellular tRNA\nsupply and demand (codon usage) of the whole transcriptome have been evolutionarily\noptimized by natural selection to ensure a close match between supply and demand.\nThis assumption allows us to approximate the supply of tRNA by using the codon usage\n(the demand) of the top 10% expressed gene (which dominate the transcriptome) in the\ntranscriptome. It is important to note that we are using the transcriptome, not the\ngenome, which means the codon usage has been weighted by the expression of those\ngenes before being summed up. Furthermore, this assumption of supply-demand match\nhas been empirically supported (PMID: 22479199 and PMCID: PMC12262485), and\nhas been used in our previous publications (PMID: 32123323 and 34977494). It is also\nimportant to note that a supply-demand match on a transcriptome level does not\nnecessarily imply a supply-demand match for individual genes. Therefore, this\nassumption is neither circular nor illogical.\n2. Using measurements of tRNA transcript abundance for tRNA supply does not change\nour conclusion.\nPer your suggestion, we obtained tRNA expression levels based on high-\nthroughput tRNA sequencing data from a previous publication (PMID: 31353208, GEO\naccession GSE128812), which were then scaled by codon-anticodon affinity and used\nas the tRNA supply. This metric is similar to the “adaptiveness value” (w) of a codon,\nas used in the popular tAI metric (PMID:15448185). It was found that the grouping\n(minimum, small, or large D) in Figure 3 remained unchanged, suggesting that our\nconclusion is not affected by the change in the method used to quantify tRNA supply.\nWe have now included this result in Supplementary Figure S2, and it is pasted below.\nFigure S2. Grouping by D of GmR strains remains unchanged if high-\nthroughput tRNA sequencing data are used to estimate tRNA supply\nThe D values of designed GmR synonymous variants were calculated with\ntRNA supply based on either transcriptomic codon usage of the top 10%\nhighly expressed genes (x axis) or tRNA transcript abundance derived from\nhigh-throughput tRNA sequencing (y axis). The black lines indicate the\ngrouping, which is identical for both types of D values. Please note that each\ndot has up to three unique sequences as biological replicates (see Table S2).\nComment 4\n4. Their estimates for genes functional importance in case of natural genes is based on\nfitness upon replacement with Kan_R. I’m not sure if this measure is adequate because\nit represents a total deletion of a gene, while here we ask about reduction in its\ntranslation efficiency. I’m not sure how to solve this problem, but I’d be better\nconvinced if they had additional measures of genes’ functional importance. Perhaps by\nemploying bioinformatics of evolutionary conservation of genes they could generate an\nestimation of functional importance but not open total deletion.\nResponse:\nThank you very much for your comment. We agree that using fitness after gene deletion\nis not ideal. Ideally, we would like to see fitness data for the same gene with slightly\nreduced expression. However, this type of data is not currently available at the genomic\nlevel. Therefore, we have resorted to the crude indirect approximation currently used in\nour study. We did not try evolutionary conservation here because it is determined by\nfunctional constraints instead of functional importance (PMID:26055156).\nPoint-by-point response to reviewers’ comments\nReviewer #1:\nOverall Comment:\nThe authors are commended on thoroughly addressing my comments.\nResponse:\nWe thank you for your thorough evaluation of our work and your very constructive suggestions.\nReviewer #2:\nOverall Comment:\nFollowing the first round of reviews, Chen and collaborators have considerably improved their\nmanuscript entitled \" A slight mismatch between a gene’s codon usage and the cellular tRNA supply\nis beneficial\", by addressing most of reviewers' comments. I will go through the main comments\nbelow.\nResponse:\nPlease accept our sincere thanks for acknowledging our improvement. Below, we will respond to\nyour further comments one-by-one.\nComment 1\nIn my view, they have satisfactorily addressed Reviewer 1's main objection related to the use of\ngentamycin, which has potential confounding effects as it affects translation. They have performed\nadditional experiments using ampicillin and an ampicillin-resistant gene in their genetic editing\nwork, which confirm their original observations.\nResponse:\nThank you for your support.\nComment 2\nThe second objection, brought forward by all three reviewers, is that their measure of distance with\nrespect to the tRNA supply (noted D in the manuscript) is an imperfect approximation. All three\nreviewers highlighted the fact that D is rather a measure of distance with respect to the synonymous\ncodon frequencies of the 10% most highly expressed genes, and not a direct measure of the distance\nwith respect to the actual tRNA supply. The authors argue that the two should match, under the\nhypothesis that translational selection has optimized the codon usage bias of the 10% most highly\nexpressed genes to match the tRNA supply. That is indeed what is expected, under the\n\"conventional\" translational selection model (which is disputed in this article), for species where\ntranslational selection is strong - that is, species with high effective population size. However, in\nspecies with low Ne, codon usage bias does not match tRNA supply even for highly expressed genes\n- again, see Bénitière et al. Nevertheless, the authors have followed the reviewers' suggestion and\nincluded in the revised manuscript a measure of tRNA supply based on tRNA gene expression levels.\nThey show that their classification of GmR versions into categories (low D, medium D and high D)\ndoes not change with this definition of tRNA supply. This confirms the validity of most of their\nresults, but not all - indeed, the analyses presented in Figure 4, 5 and 6 (which go beyond the analysis\nof E coli strains with different GmR versions) have not been redone with this measure of D. As these\nanalyses are important for the authors' model, they need to be validated with this new measure of D.\nResponse:\nThank you for your question. We have now redone all relevant results in Figure 4, 5 and 6\nusing D calculated based on the tRNA expression, whose results remained qualitatively unchanged.\nPlease note that we have also slightly updated the results derived from the mutation accumulation\nlines in order to better reflect mutation biases. More specifically, when generating random mutations\nto assess the null distribution, not only is the number of mutations matched, but also the type of\nmutations matched. This renders the original positive control for MA data in Figure S7 (currently\nFigure S10) insignificant, indicating that the sample size is insufficient for a proper test. As a result,\nwe were forced to include 15% most highly expressed genes. The tRNA expression-based results\ncorresponding to Figure 4/5/6 are pasted below.\nFigure S6. The effect of functional payoff on native codon usage of E. coli endogenous\ngenes replicated by tRNA expression-based D\nThese panels are similar to Figure 4a-b, except that the tRNA supply (𝑋 in D) is\n𝑖𝑗\nestimated using tRNA expression level instead of the codon usage of the 10% most highly\nexpressed genes.\nFigure S7. The correlation between D and expression for endogenous genes is\nrecapitulated by tRNA expression-based D\nThese panels are similar to Figure 5b-e, except that the tRNA supply (𝑋 in D) is\n𝑖𝑗\nestimated using tRNA expression level instead of the codon usage of the 10% most highly\nexpressed genes.\nFigure S11. Selective maintenance of the mismatch between endogenous gene’s codon\nusage and cellular tRNA supply is recapitulated by tRNA expression-based D\nThese panels are similar to Figure 6c-h, except that the tRNA supply (𝑋 in D) is\n𝑖𝑗\nestimated using tRNA expression level instead of the codon usage of the 10% most highly\nexpressed genes.\nComment 3\nOn the same point, I still have trouble with the circularity of the argument: the authors rely on a\nmeasure of D that only makes sense under the \"conventional\" translational selection theory to refute\nthe \"conventional\" translational selection theory. Under their model, codon usage bias should match\nthe tRNA supply only for highly expressed *and essential/important* genes, for which the\nfunctional payoff is more important than the translational cost. The authors should analyze whether\ntheir analyses are still valid when defining D based on the codon usage of the most highly expressed\nand essential/important genes. I am not suggesting that this new D measure should replace the one\noriginally proposed by the authors throughout the manuscript, as that would be circular as well - but\nrather I am asking for additional analyses, to be presented in the supplementary material.\nResponse:\nThank you for your insightful question. Please allow us to clarify that our argument is not\ncircular, and that there is a subtle difference between the conventional translational selection model\nand our assumption underlying the definition of D. Specifically, we assume that the 10% most highly\nexpressed genes COLLECTIVELY match the tRNA supply, but individually they may not\nnecessarily match it. By contrast, the conventional translational selection model assumes that\nselection optimizes the codon usage bias of INDIVIDUAL genes so as to match the tRNA supply.\nTo illustrate this, we plotted the synonymous codon frequency of each amino acid (i.e., fraction of\nsites of an amino acid being encoded by one of its synonymous codon) for each of the 10% most\nhighly expressed genes. In addition, we plotted the expression-weighted average of the synonymous\ncodon frequency as a red triangle. As shown in Figure RR1 below, the codon frequencies of\nindividual genes all deviate from the collective codon frequencies in one way or another,\ndemonstrating that the collective codon usage of top 10% highly expressed genes may differ from\nthose of the individual genes. In other words, the assumption underlying D (about the collective\ncodon usage) is different from the prediction of the conventional translational selection model\n(about codon usage of individual genes).\nFigure RR1. Synonymous codon frequencies of 10% most highly expressed genes\nThe frequency of synonymous codons in individual genes (gray points) and their\nexpression-weighted averages (red triangles) are plotted. The size of the gray point reflects\nthe expression level of the focal gene, as indicated by the legend at the bottom left. For\namino acids with two synonymous codons (bottom two rows of panels), the two codon\nfrequencies are directly shown on both axes. For amino acids with three or more\nsynonymous codons (top row of panels), the first two principal components of the codon\nfrequencies are shown, while the loadings of frequencies of the individual synonymous\ncodons are indicated by blue arrows.\nIn addition, we would like to emphasize that we are not advocating that the most important\ngenes should perfectly match the supply of tRNA. Rather, we believe that even highly expressed\nimportant genes should exhibit a slight mismatch with the tRNA supply. This has now been clarified\nin the Discussion.\n… In this work, we further demonstrated that this concept can be extended to endogenous\ngenes native to the genome. In other words, to ensure optimal cellular growth, any gene\nexpressed in a cell, including those that are functionally very important, should have a\nCUB that slightly deviates from the cellular tRNA supply, though the specific direction of\ndeviation may differ for each gene. …\nComment 4\nThe authors have addressed my comment regarding the translational effects on endogenous genes\nby performing Ribo-seq and analyzing the decoding time of codons that are frequent in the GmR.\nThey show that in E coli strains where the GmR has low D, there is a positive correlation between\ndecoding time and the frequency of codons in the GmR sequence, in particularly for genes that are\nhighly translated. I appreciate the effort that the authors have put into performing the Ribo-seq\nexperiments and I am satisfied with the result.\nResponse:\nThanks again for your constructive suggestions that pushed us to improve our analysis.\nComment 5\nRegarding the mutation accumulation experiments, the authors have followed my suggestion to only\ninclude synonymous mutations. However, I am puzzled as to why they argue that it is not useful to\nlook at weakly expressed genes. In fact, weakly expressed genes should provide a \"control\" of sorts,\nprecisely because their CUB is not expected to have been optimized by purifying selection. I would\nlike to see a comparison between highly expressed and weakly expressed genes for the analyses\npresented in Figure 6.\nResponse:\nThank you for this clarification. We agree with you that weakly expressed genes can serve as a type\nof control. We therefore randomly sampled a group of weakly expressed genes (defined as 20% least\nexpressed genes) that matched the highly expressed genes in terms of number of genes, number of\nmutations of each type (to control mutation bias) and GC content. Using the same analyses as for\nhighly expressed genes, we found no evidence for bias for D-increasing synonymous mutations.\nNote that this result should be explained by a lack of strong translational selection in weakly\nexpressed genes rather than by the selective maintenance of the slight mismatch between CUB and\ntRNA supply. The related results are now included as Supplemental Figure S9c-d, which is pasted\nbelow.\nFigure S9. The fraction of D-increasing mutations is significantly higher than random\nexpectation among genes with minimal D, but not among weakly expressed genes.\n(a) This is similar to Figure 6c and f, except that mutations among genes with 10% lowest\nD values were analyzed. This is a positive control for Figure 6c and f, demonstrating that\nthe observation therein is not due to insufficient statistical power. See Methods for details.\n(b) This is similar to panel a, except that tRNA supply (𝑋 in D) is estimated using tRNA\n𝑖𝑗\nexpression level instead of the codon usage of the 10% most highly expressed genes. (c-\nd) These are similar to panel a-b, except that a group of lowly expressed genes (defined as\n20% least expressed genes) that matched the highly expressed genes in Figure 6c and f\nwere analyzed as a control. The number of genes, number of mutations of each type in a\ngene and GC content of a gene (±15%) were matched.\nComment 6\nRegarding the analysis presented in Figure 5 (correlation between D and expression level, among\nhighly expressed genes), I am still not entirely convinced by the authors' interpretation, mainly\nbecause of the results observed for human. The corerlation between D and expression level is 0.15\nfor human and 0.13 for yeast for the top 10% most highly expressed genes, although in yeast\ntranslational selection (be it \"conventional\" or according to the authors' new model) should be much\nmore efficient than in humans. Because D is defined as the distance with respect to the CUB of the\n10% most highly expressed genes, and that the CUB is computed based on the transcriptome not on\nthe genome, I wonder if part of this correlation may be artefactual (though I confess that I have\ntrouble imagining what the confounding factor might be!). Again, here it is important to redo the\nanalysis with the D measure that is based on tRNA abundance. Some protein-coding gene families\n(e.g. ribosomal proteins) are greatly over-represented among the most highly expressed genes.\nCould that induce a biased *amino-acid* usage, that could explain some of this correlation? Perhaps\nif one could compute a D value by amino-acid, and check if all amino-acid D values show the same\npattern?\nResponse:\nThanks for your comment. We are in agreement with you that species with stronger translational\nselection should exhibit a stronger correlation between D and expression level. But the issue here is\nthat gene expression distributions may vary between species, making it difficult to directly compare\nquantile thresholds (x axis in Figure 5). Accordingly, we believe that the maximum correlation\nachieved should be a better reflection of the strength of translational selection, which is >0.2 in E.\ncoli, S. cerevisiae and D. melanogaster, but never reaches 0.2 in human.\nFor the re-analyses using D derived from tRNA expression, please refer to the above Comment\n2. We also followed your suggestion to separately analyze D for each amino acids, i.e. D.\nThe results\ni\nhave been added as Supplemental Figure S9 (too big to paste inline here). Essentially, all of these\nresults show the same pattern as our main results based on D.\nComment 7\nI am otherwise satisfied by authors' reply to my suggestions (including my questions related to YFP\ntranscription level, plasmid copy number variations and GC content effects). Overall I believe that\nthis is a valuable manuscript and that it should ultimately be published in Nature Communications,\nfollowing an additional round of review that would address the remaining points listed above.\nResponse:\nThanks again for your very thorough evaluation of our work !\nPoint-by-point response to reviewers’ comments\nReviewer #2:\nOverall Comment:\nI thank again the authors for their thorough reply. I remain puzzled by some of the results (in\nparticular the apparent high translational selection levels in humans, according to the authors'\nobservations), but overall I do not have any further objections to publication. This manuscript\ndeserves to be published and further examined by the scientific community.\nResponse:\nWe thank you for your very constructive suggestions and your support of our work."}]

Decision letter