Evolutionary conservation of codon optimality reveals hidden signatures of co-translational folding

A normalized translational efficiency scale

The classical view of translational efficiency (cTE), which is the known tRNA adaptation index¹⁹ (tAI), does not incorporate the cellular dynamics of tRNAs, driven by a trade-off between tRNA supply and demand^{1, 21} (Fig. 1b). We hypothesized that a translational efficiency scale that reflects this competition for the cellular pool of tRNAs may better capture the biological forces shaping the codon bias. We thus normalized the cellular tRNA abundances and selective constraints on codon-tRNA interactions as defined in tAI¹⁹ by the codon usage (see Methods, Supplementary Fig. 1a). How often a codon is translated in the cell depends on the codon frequencies in the mRNAs, the abundances of the mRNAs that are attached to ribosomes, and the densities of ribosomes on the specific mRNAs²⁴. We verified that mRNA abundances alone serve as a sufficient and readily available proxy for the calculation of the codon usage (Supplementary Fig. 1b–d). In this normalized translational efficiency (nTE) scale, codons are considered optimal if the relative availability of cognate tRNAs exceeds their relative usage.

Comparing the nTE to the cTE scale (Fig. 1c–e, Supplementary Fig. 2a, b), the nTE scale has a more shallow plateau-like middle region with two distinct tails of high and low efficiency codons. This stands in contrast to the cTE that scale increases almost linearly (Fig. 1e). The tail of low efficiency codons and plateau-like middle region are unique to the nTE scale, i.e. the ratio of corresponding tRNA availability and codon usage (Supplementary Fig. 2c). Importantly, this analysis suggests that tRNA supply and demand are closely matched for most codons at steady-state expression, likely supporting cost-effective proteome maintenance. Also interesting, the nTE scale contains more optimal codons that encode hydrophobic amino acids, which resonates with observations that optimal codons are associated with structurally sensitive, buried sites²⁰. The most efficient new optimal codons in the normalized scale encode glycine, the smallest thus most neutral amino acid, and arginine, cysteine, and proline, which are important gatekeeper residues that can prevent aberrant protein aggregation²⁵. The higher fraction of optimal codons encoding gatekeeper and hydrophobic residues revealed by the nTE scale likely reflects the need for increased fidelity during translation of sensitive polypeptide segments critical to correct folding and avoidance of aggregation^{20, 26}. In addition, the codons encoding the most abundant amino acid glutamine are the most balanced between supply and demand in the nTE scale. The only non-degenerate amino acids methionine and tryptophan, which are nonoptimal in the cTE scale, are optimal in the nTE scale. Importantly, the nTE scale does not, unlike the cTE scale, correlate with the codon bias; the nTE scale is also not directly correlated with the cTE scale (Fig. 1f). Thus, while the sets of optimal codons before and after normalization largely overlap, the nTE scale provides a new and independent metric of translational efficiency that reflects the cellular competition for tRNAs^{1, 21}.

A conserved translational efficiency “dip” at the start of mRNAs

Evolutionary conservation is a strong indication of functional importance. Thus, any uniform link between translational efficiency patterns in the mRNA sequences and the folding of the encoded polypeptides would expectedly be evolutionarily conserved. We undertook a systematic analysis of the conservation of codon optimality across ten closely related yeasts (see Methods). For this, we computed the nTE scale for all ten yeasts; in all cases we observed the characteristic shape of the nTE scales as observed for S. cerevisiae. Previous work using the cTE scale revealed an evolutionarily conserved region of low translational efficiency at the beginning of coding sequences, termed the “ramp”, spanning the first ca. 35–50 codons¹⁸. We thus next tested if the new nTE scale can also detect the “ramp”. The average nTE profile of the S. cerevisiae genome indeed validated the initial region of low translational efficiency (Fig. 2a), but showed that it is only ca. 10 codons long, i.e. much shorter than found using the cTE scale (Fig. 2b). Importantly, this short “dip” was evolutionarily conserved in all other analyzed yeasts (Supplementary Fig. 3). While the longer “ramp” is not observed in translational efficiency profiles of individual genes using either the cTE or nTE scales (Fig. 2c–e), this characteristic “dip” can be observed in almost all individual translational efficiency profiles, both in highly expressed genes, e.g. hydroxylase LIA1 or AAA ATPase CDC48 (Fig. 2c, d), and in lowly expressed RNA polymerase II transcription factor gene TFB3 (Fig. 2e). Interestingly, the length of the evolutionarily conserved “dip” corresponds to roughly the distance from the peptidyltransferase center to a constriction site within the ribosome^{10, 11} (Fig. 2f). The constriction site consists of ribosomal proteins L4 and L17, which are thought to sense nascent chain conformations to signal to the outside of the ribosome^{27, 28}. The finding of only a short dip of low translational efficiency also stands in good agreement with the experimental determination of ribosome densities along mRNAs in yeast, which is dominated by systematic ribosome pausing at the very start of the coding sequences²⁴. Importantly, this analysis indicates that the combination of the nTE scale and evolutionary analysis can reveal sequence signatures important for translation.

Site-specific evolutionary conservation of codon optimality

A functional link between the positioning of optimal and nonoptimal codons along the mRNA and co-translational folding would predict that codon optimality is conserved in a site-specific manner, that relates to how nascent chains fold. To test for site-specific conservation of codon optimality we constructed sequence alignments across orthologs of ten closely related yeasts, and for each gene computed a conservation score of codon optimality for each position (see Methods, Fig. 3a, Supplementary Fig 4a, b). Randomizing the alignments yielded the distributions of random conservation scores, which allowed us to determine the alignment-specific minimal conservation scores found by less than 5% chance (i.e. significance thresholds with p < 0.05). Sites in the biological alignment that have higher conservation scores than these significance thresholds are thus considered significantly conserved. If the number of significantly conserved sites in an observed alignment exceeds the corresponding number in the randomized alignment, than these positions in the sequence alignment must be under direct selective pressure for site-specific conservation of codon optimality (Fig. 3b). Because expression alone can explain the fraction of optimal codons in mRNA sequences²⁹, we devised a randomization procedure that maintains the overall codon composition in each sequence. To avoid a persistent expression bias, all analyses were performed for two curated sets of 404 alignments of high and 302 alignments of low mRNA expression levels respectively (see Methods, Supplementary Fig. 4c). Importantly, we also employed independent randomization procedures that take into account the distribution of optimal and nonoptimal codons for each amino acid in the genetic code. These latter additional analyses demonstrated that the conservation of codon optimality is completely independent of amino acid biases (Supplementary Fig. 5a–g). Furthermore, these results are independent of the 5′ coding regions (Supplementary Fig. 5h). Together, these independent analyses clearly indicate that there is site-specific evolutionary conservation of codon optimality regardless of amino acid bias or expression level.

We found that codon optimality is under selection in almost 80% of the low expression, and over 90% of the high expression genes (Fig. 3c). This is consistent with the observation that highly expressed proteins are generally more conserved³⁰. The significance thresholds for optimal and nonoptimal codons in high and low expression genes distribute homogeneously for the nTE scale, and reflect the higher content of optimal codons in highly expressed genes as expected (Fig. 3d). Strikingly, both optimal and nonoptimal codons are found conserved in equal measure in both highly and lowly expressed genes. This challenges the view that mostly optimal codons are selected for in the context of translation. Instead, our findings indicate that codon optimality is not only tuned to expression but rather fulfills an evolutionarily selected function in protein biogenesis^6–8. Of note, the same analysis using cTE to define codon optimality produces similar results for conservation of optimal codons (Fig. 3e, f). However, the cTE scale reduces the significance of conserved nonoptimal codons (Fig. 3e, f). While optimal codons are important for translational fidelity, nonoptimal codons play a key role in modulating the speed of elongation. For this reason, the nTE scale is better suited to assess the possible link between codon optimality and co-translational folding.

Conserved hidden signatures of co-translational folding

Having determined that optimal and nonoptimal codons are generally conserved raises the question regarding their role in co-translational folding events. We thus mapped the codon conservation profiles from the alignments described above onto the corresponding S. cerevisiae protein sequences and structures (e.g. Fig. 4a). Both conserved optimal and nonoptimal codons are generally distributed throughout the mRNA sequences, and often appear in clusters (Fig. 4a). We next tested for statistical associations between conserved optimal and nonoptimal codons and secondary structure elements of the encoded nascent chains (Fig. 4b, c). A control analysis using randomized synonymous codons confirmed the independence of our results from amino acid biases (Supplementary Fig. 6a, b).

We found distinct patterns of codon optimality conservation depending on the secondary structure of the encoded polypeptides, for both high and low expression proteins. Predicted α-helices are enriched in both conserved optimal and conserved nonoptimal codons, independent of expression (Fig 4b). In contrast, β-sheets are enriched in conserved optimal but depleted in conserved nonoptimal codons, in both highly and lowly expressed genes (Fig 4b). Of note, coil regions are always depleted of conserved optimal codons. Conserved nonoptimal codons are weakly enriched in highly expressed genes, and depleted in lowly expressed genes (Fig. 4b). Interestingly, α-helices can already form co-translationally even within the ribosomal tunnel^{10, 11}, while coil regions, comprising loops that fold near the exit of the ribosomal tunnel, have been shown to play key roles in co-translational protein folding⁷. In contrast, β-sheet containing domains are topologically discontinuous and must await synthesis to begin folding. Furthermore, β-sheets are characterized by their high content of hydrophobic residues, the presence of gatekeepers, and a high aggregation propensity, thus, the general strong enrichment of conserved optimal as well as depletion of conserved nonoptimal codons could primarily serve to reduce the risk of phenotypic missense mutations leading to aggregation.

Since hydrophobicity is linked to both protein folding and aggregation, we also tested for associations between conserved codon optimality and hydrophobicity. As expected, conserved optimal codons are enriched in hydrophobic regions, and conserved nonoptimal codons are depleted. This observation is stronger for highly expressed genes (Supplementary Fig. 6c), likely owing to a greater need for translational fidelity in these abundant proteins^{3, 20, 26}.

We next considered only sites that appear in clusters. This analysis provided a more stringent test for our results that validated and increased the significance of all the above associations between codon optimality and secondary structure propensity (Fig. 4c). In particular, the enrichment of conserved nonoptimal codons appearing in clusters is much stronger in α-helices at both high and low expression levels, and in coil regions of highly expressed genes.

To validate these observations with confirmed secondary structures, we analyzed 357 experimentally determined protein structures from the protein data bank (PDB). Importantly, PDB structures allow to distinguish between more defined structural elements such as hydrogen-bonded turns that often connect the more regular elements of secondary structure, α-helix and β-strand, and are important for folding into the final structure. This analysis indicated that conserved nonoptimal codons are strongly enriched in turns (Fig. 5a, Supplementary Fig. 6d, e). The enrichment of conserved nonoptimal codons only in α-helices and in turns is remarkable, since these secondary structure elements can form co-translationally^{10, 11}. These findings are not detected when using the classical definition of codon optimality (Fig. 5b). All the above findings for conserved optimal codons are weaker using cTE and all significant associations to conserved nonoptimal codons are lost (Fig. 5b). This may explain why the link between codon optimality and co-translational folding has been hard to detect. Individual examples support our findings. For instance, the all-helical Myosin light chain 1 shows alternating patterns of conserved optimal and nonoptimal codons in all its helices (Fig. 5c). The 20S proteasome subunit G is characterized by a clear conservation of optimal codons in its structural core of β-sheets, complemented by conservation of both optimal and nonoptimal codons in its helices and turns (Fig. 5c).

What stands out from this analysis is the enrichment of both conserved optimal and nonoptimal codons in α-helices. The formation of α-helices is one of the elementary steps in protein folding, characterized by complex kinetics due low cooperativity that depend on an initial nucleation step³¹. Furthermore, this is the main folding event that has been found to occur inside the ribosome exit tunnel^{11, 32, 33}; indeed helix formation has been experimentally observed to occur very early in the ribosomal tunnel, between the peptidyltransferase center and the constriction site²⁸. Strikingly, for the α-helices in the experimental protein structures of S. cerevisiae, we found a distinct alternating pattern of codon optimality. Specifically, we observed a preference for a nonoptimal codon at the transition into the helix, followed by an enrichment of optimal codons at position 1 and 4, interspersed with a strong preference for nonoptimal codons at positions 2 and 3 (Fig. 5d). Moreover, the significance of this profile only extends across the first full helix turn, independent of helix length (Fig. 5d). This suggests that codon optimality may be evolutionarily selected to tune the translation and folding rates of helices early in their entry into the ribosomal tunnel. The ribosomal tunnel has been found to possess distinct folding zones that may facilitate helix formation deep inside the ribosome³³, with the strongest compaction into secondary structure observed proximal to the peptidyltransferase center³³. Helix formation inside the ribosomal tunnel appears strongly sequence dependent, but cannot be explained by sequence hydrophobicity or helical propensity alone³³. The specialized environment within the tunnel is proposed to create a rugged solvation landscape that may slow down folding³⁴. We speculate that the evolutionarily conserved patterns of codon optimality in helices could facilitate specific interactions with the exit tunnel wall, and may even assist helix nucleation inside the exit tunnel.

Data sources

Genomic sequences and ortholog assignments for S. cerevisiae, C. glabrata, D. hansenii, K. lactis, S. bayanus, S. kluyveri, S. mikatae, S. paradoxus, S. pombe, and Y. lipolytica were retrieved from the Broad Institute (http://www.broadinstitute.org/regev/orthogroups)⁴⁰, and alignments of the genetic sequences between orthologs were computed with ClustalW⁴¹ via the corresponding amino acid sequences. tRNA counts were retrieved from the tRNA database (http://www.gtrnadb.ucsc.edu) if available, otherwise predicted with tRNAscan-SE from the genomic sequences⁴². mRNA expression levels for S. cerevisiae were obtained from⁴³, for S. pombe from⁴⁴, and for all other yeasts from⁴⁵. 1574 alignments contained more than 7 sequences, our requirement to obtain meaningful conservation scores. To remove any intrinsic expression bias, we normalized the expression levels across the yeasts by quantile normalization, and selected 500 alignments each of high and low expression and with the lowest inner expression divergence measured by the standard deviation across orthologs. Removing alignments with more than 30% gaps yielded a set of 404 alignments of highly, and 302 alignments of lowly expressed genes. The mapping of PDB structures to S. cerevisiae genes was obtained from⁴⁶, yielding a curated set of 357 alignments.

Protein secondary structures were predicted with PSIPRED⁴⁷. RNA secondary structures were predicted with the Vienna package⁴⁸. For alignments with assigned PDB structure, secondary structure and relative accessible surface area were extracted with the DSSP program⁴⁹. Sequence hydrophobicity was computed using the Kyte&Doolittle scale.

Translational efficiency and codon optimality

The classical translational effiency cTE_i for each codon is the published tRNA adaptation index (tAI) as computed with the codonR program¹⁹. It estimates the tRNA availability for each codon i from a weighted sum of the gene copy numbers tGCN_ij of the matching tRNA isoacceptors j under a selective constraint s_ij on the efficiency of the codon–anticodon coupling, incorporating Crick’s wobble rules¹⁹: n_i

The selective constraint on codon–anticodon interactions s_ij is 0 for cognate tRNAs, and small for wobble interactions¹⁹. The overall efficiency of a codon is thus given by the sum of the contributions of the recognizing individual tRNAs under consideration of specific selective constraints based on the codon–anticodon interaction¹⁹. The division of the individual translational efficiencies W_i by the maximum translational efficiency W_max linearly rescales all translational efficiencies so that the maximum value is 1.

The codon usage cu_i was defined as relative estimate of how often each codon is translated. It is derived from the number of occurrences of each codon in an ORF, weighted by the corresponding transcript abundance, and summed up over all ORFs (Supplementary Fig. 1a). For codon i this is the sum of the counts c_ij of the codon i in gene j, weighted by the transcript abundance a_j of gene j, considering all genes in the genome g. For comparability, the codon usage is also rescale to have a maximum value of 1.

In this work, the normalized translational efficiency nTE_i is subsequently defined as the ratio of tRNA availability cTE_i (supply), which is based on cellular tRNA abundance and selective constraints for wobble interactions, and codon usage cu_i (demand), linearly rescaled to have a maximum value of 1.

Codons i with cTE_i ≥ cu_i are considered optimal, and nonoptimal otherwise. For comparison, we used the set of classical optimal codons reported in²⁰. They are those found significantly enriched in the highest expressed genes by a Chi-square test²⁰. Average translational efficiency profiles were computed as described in reference 18.

Randomization procedure & significant sites

We tested for evolutionary conservation of codon optimality in 10 closely related yeasts. The conservation score S of optimal or nonoptimal codons at any given position i is defined as S_i = n_i/ N_i, where n_i is the number of optimal or nonoptimal codons at position i, and N_i the total number of aligned codons at position i respectively. We only considered alignments with at least 7 orthologous sequences, and a minimum of 5 codons has to be aligned for a conservation score to be computed at that position. Each alignment of orthologs was randomized 1000 times by individually shuffling each sequence to keep its original composition and maintain the individual codon bias. We employed additional randomization schemes to verify independence of amino acid biases (Supplementary Fig. 5,6). From the distribution of the conservation scores in the randomized alignments, we extracted alignment-specific the minimal conservation scores that are observed at less than 5% chance as significance thresholds. We considered sites with higher conservation scores than the significance thresholds as significantly conserved. If more optimal and nonoptimal codons respectively are significantly conserved in the biological alignment than in the randomized alignment, the site-specific evolutionary conservation of optimal and nonoptimal codons respectively is the result of selective pressure.

Individual and average translational efficiency profiles

Individual translational efficiency profiles were computed with the nTE and cTE scales, and smoothed with a sliding window of size 15, the size of the immediate ribosome footprint¹⁸. Average translational efficiency profiles were computed by aligning all genes at the start codon, and subsequently calculating the average translational efficiency for each position as described in¹⁸. Thus, the average profiles only reflect those fluctuations of the individual profiles that are present in all sequences. We randomly reshuffled all sequences to calculate the mean and standard deviation for each position¹⁸.

Statistical testing

All statistical testing of associations was performed using the Cochran–Mantel–Haenszel test in the statistical computing environment R (http://www.r-project.org), and as described in reference 20. The enrichment of optimal or nonoptimal codons at specific helix positions in S. cerevisiae sequences was tested with Fisher’s exact test and corrected for multiple-testing with the Benjamini & Hochberg (1995) method. All test statistics and definitions of optimal and nonoptimal codons are listed in the Supplementary Information (Supplementary Tables 1–8).

Data availability

All datasets of this study are available at http://www.stanford.edu/group/frydman/codons.