An improved estimation of tRNA expression to better elucidate the coevolution between tRNA abundance and codon usage in bacteria

The degree to which codon usage can be explained by tRNA abundance in bacterial species is often inadequate, partly because differential tRNA abundance is often approximated by tRNA copy numbers. To better understand the coevolution between tRNA abundance and codon usage, we provide a better estimate of tRNA abundance by profiling tRNA mapped reads (tRNA tpm) using publicly available RNA Sequencing data. To emphasize the feasibility of our approach, we demonstrate that tRNA tpm is consistent with tRNA abundances derived from RNA fingerprinting experiments in Escherichia coli, Bacillus subtilis, and Salmonella enterica. Furthermore, we do not observe an appreciable reduction in tRNA sequencing efficiency due to post-transcriptional methylations in the seven bacteria studied. To determine optimal codons, we calculate codon usage in highly and lowly expressed genes determined by protein per transcript. We found that tRNA tpm is sensitive to identify more translationally optimal codons than gene copy number and early tRNA fingerprinting abundances. Additionally, tRNA tpm improves the predictive power of tRNA adaptation index over codon preference. Our results suggest that dependence of codon usage on tRNA availability is not always associated with species growth-rate. Conversely, tRNA availability is better optimized to codon usage in fast-growing than slow-growing species.


Results
RNA-Seq mappings are consistent with tRNA abundance estimates. To accurately profile tRNA transcripts in tpm, we processed RNA-Seq data to remove adapters and low-quality sequences (See Materials and Methods for more detail) and quantified reads mapped to all unique tRNA sequences (Supplementary File S1) in tRNA tpm. To demonstrate the fidelity of tRNA tpm in bacteria, we compared these values with RNA fingerprinting abundances (Fig. 1, Supplementary File S2) previously reported in E. coli 10,42 , S. enterica 10 , and B. subtilis 11 . In all cases, tRNA tpm correlates with RNA fingerprinting abundance ( Fig. 1: R 2 > 0.4, P < 0.05).
Documented tRNA methylation does not appreciably affect tRNA sequencing in bacteria studied. To investigate the potential effects of site-specific tRNA methylation on standard RNA-Seq experiments in bacteria, we visualized the RNA-Seq read depths for all seven species studied (Fig. 2, Supplementary File S1). In E. coli, read depths before and after documented tRNA methylation sites (18, 32, 34, 37, 46, and 54) in GenBank annotation (NC_000913, Supplementary File S1) do not vary substantially, and we observe no partial tRNA mappings (Fig. 2a-d) contrary to the "hard-stops" previously described in both yeast and human tRNAs in the absence of demethylation treatment 31 . Additionally, tpm values associated with the set of tRNAs that can be potentially modified at five or all six documented methylation sites do not differ substantially from the set of tRNAs that can be potentially modified at four or less sites (Fig. 2d,e; two-tailed Student's t-test with unequal variance: P = 0.477). We define tRNAs that can potentially be methylated at >4 sites as heavily methylated with respect to other tRNAs. This cut-off was chosen because it divides the 50 unique tRNA sequences into two subsets of roughly equal size. Similarly, hard-stops were not observed at documented methylated sites in all six other bacteria studied (Supplementary File S1).
An improved estimation of codon preference using tRNA tpm. To investigate how well tRNA tpm explains codon preference across bacterial lineages, we first designed an integrative approach to determine www.nature.com/scientificreports www.nature.com/scientificreports/ translationally optimal codons. We use two criteria to define what constitutes a translationally optimal codon: 1) it has the highest relative synonymous codon usage (RSCU) 45 in HEGs, and 2) its RSCU in HEGs is greater than that in lowly expressed genes (LEGs). These two criteria combine the specifications of optimal codons defined by two previous studies: the first criteria is used by Novoa,et al. 19 and the second is adapted from Rocha 24 (HEGs vs. others). These criteria are also consistent with those used in the calculation of the Codon Adaptation Index 4 and Index of Translation Elongation 5 . By our definition, there can be only one codon that is most translationally optimal within each synonymous group, consistent with the original definition of optimal codon determined by F op 10 . Thus, the characterized translationally optimal codons will increase elongation efficiency relative to others. The availability of protein abundance data in parts per million (ppm) for species studied herein and mRNA transcript abundance (in tpm) determined using kallisto 26 enables us to calculate protein per transcript (ppm/tpm) in the identification of HEGs and LEGs (see Materials and Methods for more detail).
To establish readable isoacceptor tRNA content, we consider all cognate and near-cognate interactions, as well as those enabled by anticodon modifications. Most studies 10,11,18,42 consider anticodon-codon cognate (e.g., tRNA Arg UGC reading GCA) and near-cognate (e.g., tRNA Arg UGC reading GCG) pairings; while some 19,24 also consider pairings allowed due to anticodon modifications (e.g., tRNA Arg UmGC reading GCC and GCU) which increases the predictive influence of tRNA content on codon usage 19 . In order to explain RSCU, we adapt the idea of Relative tRNA Gene frequency from Novoa, et al. 19 to derive Relative tRNA Usage (RTU) (see Materials and Methods for more detail). By considering tRNA abundance for synonymous codons, RTU improves tRNA tpm as an estimator of tRNA abundance (Fig. 3). In particular, the correlation between tRNA tpm and average tRNA abundance 42 in E. coli ( Fig. 1a; R 2 = 0.498, P < 0.0001) improves if we consider their RTU values ( Fig. 3a: R 2 = 0.646, P < 0.0001). Furthermore, both RTU values correlate with codon usage (Fig. 3b).
Using RTU, we estimate how well translationally optimal codons match codons with the highest tRNA availability by adapting the four rules of codon-anticodon constraint 10 . Rule one states that codon usage is constrained by tRNA availability 46 . Rules two to four focus on specific base pairing efficiencies and they describe that cognate codon-anticodon pairs are generally more efficient and preferable relative to near-cognate wobble pairs 22,[47][48][49][50] . Hence, we first rank synonymous codons by highest RTU, and then rank by cognate tRNA abundances (Supplementary File S2). Supplementary Fig. S1 provides a flowchart explaining the approach to predict translationally optimal codons by tRNA availability. For example, applying this two-step identification approach for Threonine codons in E. coli, we first select ACC and ACU since they both have the highest RTU and are both readable by the same tRNAs (tRNA Thr GGU and tRNA Thr UGU ) due to anticodon modification by ADATs 19 . Next, we rank codon preference by cognate tRNA abundance: tRNA tpm of the cognate tRNA Thr GGU for ACC is 46537.7, and ACU is not decoded by a cognate tRNA. Thus, the predicted optimal codon based on tRNA availability and pairing constraints is ACC in Threonine. Indeed, ACC is the translationally optimal codon for Threonine www.nature.com/scientificreports www.nature.com/scientificreports/ based on our definition: it has the highest RSCU in HEGs (2.111), and (2) its RSCU HEG is higher than RSCU LEG (1.468). We extended this two-step approach to predict translationally optimal codons in all seven species (Table 1).
Many studies have demonstrated that bacterial tRNA abundance fluctuates due to growth phase and culture conditions such as temperature and media 42,[51][52][53] . Hence, for each species, we have acquired tRNA tpm values from two independent, but experimentally consistent (i.e., all strains are wildtype, and all cultures are taken during log-phase growth) RNA-Seq datasets ( Table 2) that were prepared for sequencing on the same platform (Illumina) to verify consistency in predicting translationally optimal codons using tRNA tpm (Table 1). Implementing tRNA tpm in tAI calculation improves the non-parametric S correlation. We calculated tAI using gene copy number and tRNA tpm (See Materials and Methods for more detail, Supplementary  Fig. S2), and Table 3 shows the non-parametric regression S which reflects the correlation between tAI values and effective number of codons 54 (corrected for silent substitutions 18 ). For six out of seven species studied, S' calculated with tRNA tpm is higher than S calculated using tRNA gene copy number (Table 3). Hence, tAI' values calculated using tRNA tpm better correlate with codon usage than tAI values obtained from tRNA gene copy number (two-tailed Student's t-test with unequal variance: P < 0.0001; Supplementary Fig. S2) in all species except L. interrogans. Note that S values calculated herein use non-hypothetical and non-pseudo genes, whereas those originally calculated (S 0 ) 18 use all coding DNA sequences, although value ranks stay consistent ( Table 3).
Abundance of tRNA depends on codon usage in fast-growing species. Codon usage preferences can also drive up tRNA content 6,24 . For each tRNA species (distinguished by anticodon), we define its readable codon usage as the sum usage of its cognate and near-cognate codons (e.g., the readable codon usage of tRNA Ala GGC = GCC usage + GCU usage). Readable codon usage in HEGs better explains tRNA tpm in fast-growing species than in slow-growing species. More specifically, readable codon usage better correlates with tRNA tpm in E. coli, S. enterica, and B. subtilis, than in B. thetaiotaomicron, Synechocystis sp., M. tuberculosis, and L. interrogans ( Supplementary Fig. S2).

Discussion
Two approaches have been taken to characterize the effect of tRNA on codon usage. The first tests whether gain or loss of tRNA genes will lead to predicted changes in codon usage [55][56][57] . In tunicates and bivalves, an additional tRNA Met/UAU gene is present in the mitochondrial genome. One would expect that the additional tRNA Met/UAU would favor increased usage of AUA codons, and this expectation is empirically substantiated 55,56 . One may also reason that, if a bacteriophage encodes many tRNA genes in its own genome, especially when these tRNAs are www.nature.com/scientificreports www.nature.com/scientificreports/ rare in the host, then the phage codon usage will be less dependent on the host tRNA pool. This expectation is also consistent with empirical evidence 16,58 . The second approach quantifies within-species association of codon usage with tRNA abundance which is often approximated by tRNA gene copy number, but this proxy has two key shortcomings. First, we do not know if it is generally true that tRNA gene copy number serves as a good proxy for tRNA abundance. Second, some bacterial species, such as M. tuberculosis, have only a single tRNA gene for all anticodons. In such cases, we cannot use tRNA gene copy number as a proxy of abundance because of its lack of variability. Thus, accurate quantification of tRNA abundance is crucial for detecting codon adaptation.
Our RNA-Seq-based analysis shows that snapshots of bacterial tRNA pools can be captured using RNA-Seq profiling in spite of claims that standard Illumina RNA Sequencing protocols inefficiently quantify tRNAs in eukaryotes 32,43,44 due to a number of modifications that increase the stability of tRNA secondary structure. While some tRNA modifications are shared among the three kingdoms of life, others are kingdom specific 59 . Nonetheless, tRNA post-transcriptional methylation is extensive in Eubacteria 59 , and we acknowledge that tRNA tpm values may underestimate tRNA abundances since tRNA is highly structured and difficult to denature. Despite this, we do not observe any drastic read count drops or hard-stops in mapped reads at or flanking documented methylated sites in the seven species studied here, nor do we see partially mapped tRNA regions (Fig. 2a-c; Supplementary File S1) that were commonly observed in untreated tRNA sequencing data in Eukaryotes 31 . A caveat is the presence of a hard-stop in most tRNAs at site 50 for B. subtilis (SRX2804667), Synechocystis sp. (SRX4145044) and L. interrogans (SRX2448246) (Supplementary File S1). However, there is no documented methyltransferase activity acting at this site for these species, and the drop in mapped reads is only observed in one of the two SRX datasets retrieved for each species. Nonetheless, we acknowledge that tRNA sequencing may be potentially enhanced with demethylase treatments 31,32 that should be employed in future tRNA-Seq studies in bacteria.
Because eukaryotic tRNA read mapping abundances are considerably higher in demethylated samples than untreated samples 31 , we expected that our read mapping abundances may be similarly impacted by tRNA methylation. In light of this, we considered all annotated methylation sites in bacterial tRNAs, even though they may not be methylated at all times due to structural constraints. Surprisingly, our observations reveal that read mapping abundances of E. coli tRNAs that are potentially heavily methylated do not differ from those that are susceptible to methylation at fewer than five methylation sites (Fig. 2e). This suggests that the requirement for demethylase treatment prior to tRNA sequencing in bacterial species with functional AlkB demethylase homologs may be more relaxed, especially since demethylation treatments in current eukaryotic tRNA sequencing approaches are bacterial AlkB-facilitated 31,32 . www.nature.com/scientificreports www.nature.com/scientificreports/ In all bacteria, we combined our RNA quantification approach with available protein abundances to determine the most translationally optimal codon in 21 codon groups (Table 1) based on codon usage differences between HEG and LEG subsets. These subsets of genes are established based on protein per transcript (Supplementary  Table S1, File S3), with the aim of providing a more accurate estimate of translation efficiency from protein  Table 1. Translationally optimal codons estimated for synonymous groups in seven bacterial species. Bold are translationally optimal codons that also have the highest tRNA availability estimated by tpm from two independent RNA-Seq datasets. 1 Two amino acids (Met and Trp) are omitted because they are each encoded by a single codon. 2 The 6-fold degenerate codon families (Leu, Arg, and Ser) are broken into 2 and 4-fold families because of differences in the first codon base. 3 An optimal codon cannot be determined by our definition (e.g., RSCU HEG < RSCU LEG violates the second criterion). * Predicted preferred codons match optimal codons in only one of the two RNA-Seq data analyzed. a Translationally optimal codons match optimal codons determined by F op in Ikemura 10 and Kanaya, et al. 11 . b Translationally optimal codons identified using average tRNA abundance from RNA fingerprinting approach, retrieved from Dong, et al. 42 .  Table 2. The seven bacterial species studied herein due to their availability of protein abundance, growth rate and RNA-Seq data. a Information on species growth-rate are retrieved from Rocha 24 . * All selected datasets have matching species strains between protein abundance, GtRNAdb, NCBI and RNA-Seq data, with the exception of L. interrogans (SRX2448246) due to the lack of a second RNA-Seq data in GEO Datasets for strain Fiocruz L1-130. All cultures in RNA-Seq experiments were isolated during log phase of growth. The first listed SRX dataset was selected for all analyses, and both were used for Table 1 and Fig. S4. www.nature.com/scientificreports www.nature.com/scientificreports/ abundance by decoupling rates of transcription. However, species-specific translationally optimal codons cannot always be established because the two described criteria are violated in some groups (e.g., the codon with the highest RSCU is more over-represented in LEGs than HEGs). In these cases, there is no evidence that the most abundantly used synonymous codon would contribute to increase translation efficiency. In particular, L. interrogans represents a slow-growing species wherein codon optimization is very poor (only seven translationally optimal codons can be characterized in 21 codon groups).

Species
Our tRNA quantification approach better predicts translationally optimal codons over F op 10 . In E. coli, synonymous codons with the highest tRNA tpm (ranked by highest TPU followed by highest cognate tRNA abundance) match 15 out of 17 translationally optimal codons, but 13 out of 17 when we replace tRNA tpm with averaged RNA fingerprinting abundance retrieved from Dong, et al. 42 (Table 1). In contrast, F op determines 12 such translationally optimal codons 10 of which AGA (Arg) is the only codon that our method does not predict to be optimal (Table 1). Additionally, all translationally optimal codons determined herein are consistent with optimal codons determined by F op , except in the Serine 4-fold family where F op predicts UCC whereas we predict UCU to be optimal. In the case of B. subtilis, both tRNA tpm and F op determine six translationally optimal codons (Table 1). It is worth mentioning that F op determines 16 optimal codons in B. subtilis 11 , but most may not be translationally optimal. For example, CCA (Pro) was determined to be an optimal codon, but CCG (Pro) is substantially more preferred than CCA (Pro) (RSCU in HEGs are 1.735 and 0.782, respectively).
Nonetheless, bacterial tRNA abundance may not fully explain the variation in usage of all 61 sense codons ( Supplementary Fig. S4). First, codon preference cannot always be inferred reliably from tRNA gene redundancy or experimentally measured tRNA abundance. For example, inosine is expected to pair best with C and U, but less with A (presumably because of the bulky I/A pairing involving two purines) 60 . Second, what matters in translation elongation is the availability of charged tRNAs. It is difficult to determine the level of charged tRNAs, and researchers typically would use transcriptionally determined tRNAs or even the number of tRNA genes in the genome as a proxy of charged tRNAs. Unfortunately, the abundance of tRNAs does not always reflect the abundance of charged tRNA 61 . Lastly, other factors such as mutation bias 21,62-65 may exert more pressure on codon usage in certain species.
Conversely, the variation in tRNA tpm is better explained by codon usage (Supplementary Fig. S3) in fast-growing (E. coli, B. subtilis and S. enterica) than in slow-growing species (B. thetaiotaomicron, L. interrogans, M. tuberculosis and Synechocystis sp.). This result supports the theory that tRNA translation machinery is better optimized to codon usage in fast-growing than slow-growing species 9,24 . Indeed, duplicating tRNA genes is an effective way to elevate transcript abundance in species that grow and replicate rapidly [10][11][12]42 , but not in slow-growing species (Supplementary Fig. S5).
One potentially important implementation of tRNA tpm is in the calculation of tAI. Our results (Table 3, Supplementary Fig. S2) show that tAI' (calculated using tRNA tpm) better explains effective number of codons than tAI (calculated using tRNA gene copy number) for all species studied except L. interrogans. Considering S and S 0 calculated using tRNA gene copy numbers, their differences are likely due to our usage of the subset of non-hypothetical and non-pseudo genes that have protein abundance values 30 (S) whereas all DNA coding sequences (including hypothetical and pseudo genes) were used in the original calculation 18 (S 0 ). Additionally, both the GtRNAdb 25,66 and DNA coding sequences (GenBank annotations) have been continuously curated since 2004. Lastly, only wobble base pairings were considered in the original introduction of tAI 18,67 ; whereas we have also considered possible anticodon modifications 19,24 that further relax codon pairing. These differences improve the calculation of the S correlation using tRNA copy numbers, notably in B. subtilis and M. tuberculosis. In contrast, the originally calculated negative S 0 correlation for B. subtilis was a major shortcoming of the tAI method 18 and was criticized 21 for suggesting a lack of selective pressure exerted by tRNA abundance on codon preference in this species.
In the case of M. tuberculosis, our tRNA quantification approach is much more sensitive to determining tRNA-mediated codon bias than tAI. The S correlations are consistently the lowest for M. tuberculosis (Table 3), yet we identified 17 out of 19 translationally optimal codons using tRNA tpm (Table 1). Our method recaptures the "weak but significant codon usage preference" previous reported 21,68 in this slow-growing species, and show that the degree to which tRNA availability explains optimal codon usage is species-specific and does not always depend on growth-rate.  Table 3. Non-parametric regression S correlations between tAI values and effective number of codons. a S values retrieved from dos Reis, et al. 18 , calculated using all coding DNA sequences. b S values calculated using tRNA gene copy number, using genes having non-zero protein abundances. c S values calculated using tRNA tpm, using genes having non-zero protein abundances.
www.nature.com/scientificreports www.nature.com/scientificreports/ We studied the coevolution between codon usage and tRNA abundance in three fast-growing species (E. coli, S. enterica, and B. subtilis) and four slow-growing species (B. thetaiotaomicron, L. interrogans, M. tuberculosis, and Synechocystis sp.). Our findings indicate that tRNA quantification by tpm offers better predictions of translationally optimal codons over F op in E. coli, and improves the calculation of tAI to better reflect codon preference in all species studied except L. interrogans. The usage of translationally optimal codons can be well explained by relative tRNA tpm in E. coli and S. enterica; however, both tRNA tpm and RNA fingerprinting abundances 11 offer weaker explanations for codon preference in B. subtilis. The influence of tRNA availability on codon bias is not always stronger in fast-growing species, and optimal codons can be well explained by tRNA content in certain slow-growing species such as M. tuberculosis. Conversely, the tRNA translation machinery is better optimized to codon usage in HEGs of fast-growing than slow-growing species.

Materials and Methods
Processing genomic, proteomic and RNA-seq data. We retrieved the annotated genomes (Table 2) for three fast growing species (E. coli, B. subtilis, and S. enterica) and four slow-growing species (B. thetaiotaomicron, L. interrogans, M. tuberculosis, and Synechocystis sp.) in GenBank format from the National Center for Biotechnology Information (NCBI) database (http://www.ncbi.nlm.nih.gov). For each species, all documented tRNA methyltransferase genes were retrieved from GenBank annotations (Supplementary File S1).
Protein abundance data corresponding with each of these species were retrieved from PaxDb 4.0 30 https:// pax-db.org/ and abundance values were associated with GeneIDs retrieved from the genomes using DAMBE 7 69 . The integrated protein abundance dataset was taken when available.
RNA-Seq runs of wildtype species were fetched from GEO DataSets (https://www.ncbi.nlm.nih.gov/gds/) in FASTQ format. FASTQ files were converted to FASTQ+ format using ARSDA 1.1 27 in order to reduce file sizes by grouping identical reads under a single ID while retaining the copy number for each read (in the format S<read#>_<copy#>). The FASTQ+ data was then processed using both CutAdapt 1.17 70 and Trimmomatic 0.38 71 to remove flanking adapter sequences and purge low quality reads. For experiments that use the oligo(dT)-adapter primer for cDNA synthesis, RNA fragments are first poly-adenylated at the 3′ end. In these cases, we set CutAdapt to recognize "AAAAA". We also used CutAdapt to recognize and remove all possible adapters in experiments that used custom adapters, with 10% mismatch error rate. When the adapters conformed to standard Illumina protocols, we simply used the "ILLUMINACLIP" function built into Trimmomatic to trim the relevant library of adapters from reads (ILLUMINACLIP: <adapters.fa> :2:20:6). After adapters were trimmed, we retained reads that were a minimum of 25 nt long (-m 25 in CutAdapt or MINLEN:25 in Trimmomatic) to mitigate bias in expression levels 72 . We filtered trimmed reads to remove poor quality sequences with average Phred scores lower than 30 (0.1% probability of a base calling error) 73 .

RNA-Seq read mapping for tRNAs and mRNAs.
We retrieved the sequences of all genomically encoded tRNAs for each organism from the Genomic tRNA Database (GtRNAdb 2.0 66 ) in FASTA format and removed predicted pseudo-tRNAs and those with unspecified anticodons. The FASTA files containing tRNA sequences were read into DAMBE to represent identical sequences with one ID indicating the number of identical copies. Since mature tRNAs are modified to have 5′-CCA-3′ appended to their 3′ end, we manually added CCA to sequences lacking this motif 32 . The modified tRNA FASTA files for all species were indexed and the associated RNA-Seq reads from processed FASTQ+ files were pseudo-aligned to each tRNA index and tRNA tpm was quantified using kallisto v0.44.0 26 . The tRNA pseudo-alignments for each species were subsequently sorted and site-specific depth values were generated for each tRNA using the sorted pseudo-alignments via the 'sort' and 'depth' commands from SAMtools 74 , respectively.
Similarly, all non-pseudo and non-hypothetical DNA coding sequences with non-zero protein abundance values were retrieved using DAMBE in FASTA format and indexed. The associated RNA-Seq reads from processed FASTQ+ files were pseudo-aligned to each mRNA index and mRNA tpm values quantified using kallisto.

Determination of translationally highly and lowly expressed genes by protein per transcript.
Protein per transcript (ppm/tpm) was estimated by taking gene protein abundance (ppm) divided by its mRNA tpm, for both RNA-Seq datasets in each species except B. subtilis. For B. subtilis, protein per transcript was obtained with only SRX515181 dataset, because SRX2804667 MiSeq experimental protocol effectively removes large transcripts to study tRNAs 75 . From each dataset, genes with top and bottom 30% ppm/tpm values were selected (Supplementary File S3). A gene is considered to be highly expressed if the gene ID is found in both gene sets for each species; the same was done to identify lowly expressed genes from the bottom 30% ppm/tpm gene sets (Supplementary File S3, Table S1). To verify the validity of this approach, we determined the number of ribosomal protein (30S and 50S subunit) genes that are present in each gene sets. We observed a great number of ribosomal protein genes in the top 30% ppm/tpm gene sets and nearly none in the bottom 30% ppm/tpm gene sets (Supplementary Table S1). This is expected because ribosomal protein genes are commonly accepted and used as highly expressed genes 4,24 .
Computation of relative synonymous codon usage and tRNA usage metrics. Relative synonymous codon usage (RSCU) 45 values were computed for each species by loading HEGs and LEGs (Supplementary File S3) into DAMBE and selecting "Seq. Analysis" > "Codon Usage" > "Relative synonymous codon usage". DAMBE's implementation of the RSCU computation automatically splits 6-fold degenerate codon families into a 2-fold and 4-fold degenerate family based on difference at the first codon position.
To acquire relative tRNA usage for each synonymous codon (RTU), we adapt the RSCU formula (1) in the same way that Relative tRNA Gene frequency was employed by Novoa, et al. 19 using tRNA gene copy number: www.nature.com/scientificreports www.nature.com/scientificreports/ where i is any codon within a 2 or 4-fold degenerate codon family, and n is the total number of codons in the synonymous group. RSCU and RTU are both calculated by breaking 6-fold codon families (Arginine, Leucine, and Serine) into 2 and 4-fold groups (e.g., 4-fold CGN (arg) and 2-fold AGN (arg)). Methionine and Tryptophan have been omitted from RSCU and RTU calculation because they are encoded by a single codon (RSCU and RTU = 1). Similarly, codon groups with RTU = 1 have been omitted from plots when applicable ( Supplementary Fig. S4), because RTU will not estimate RSCU for these codons. All correlation coefficients (R 2 ) are calculated by taking the square of Pearson's correlation r (Figs 1 and 3, Supplementary Figs S3-5).
Computation of tAI and correlation S. We first calculated tAI using the original formulation of the model 18,67 which considers the copy number of all isoacceptor tRNAs for each codon via the author's tAI R package version 0.2 (https://github.com/mariodosreis/tai) for all species in this study. We additionally computed a modified version of tAI (tAI') which uses the summed tpm values associated with codon-specific isoacceptor tRNAs in lieu of tRNA copy number. Rather than using all annotated DNA coding sequences in these calculations, as was done originally 18 , we only considered genes with non-zero protein abundance values from the integrated datasets stored in PaxDb. The tAI and tAI′ values for each species were plotted ( Supplementary Fig. S2) against the effective number of codons corrected for silent substitutions at the third codon position (f[GC3s] -Nc) to determine the S and S' correlation coefficients, respectively.

Data Availability
Supplementary File S1 contains RNA-Seq read depths and tRNA methylation profile. Supplementary File S2 contains identifications of translationally optimal codons and tRNA abundances (gene copy, tpm and fingerprinting data) and tRNA quantification approaches. Supplementary File S3 contains protein per transcript data, protein abundance information, identified translationally HEGs and LEGs. Supplementary File S4 contains Supplementary Figs S1-S5 and Table S1.