Evidence of evolutionary selection for co-translational folding

TL;DR

This study provides genome-wide evidence that evolutionary selection favors specific translation rates at certain gene regions to optimize co-translational protein folding, highlighting the importance of kinetics in protein self-assembly.

Contribution

The paper introduces a statistical method to identify conserved, slowly translated codon regions linked to co-translational folding intermediates across many proteins.

Findings

Regions enriched in slow codons are associated with folding intermediates.

Domain boundaries explain only a small part of conserved slow codon regions.

Synonymous substitutions are often selected to optimize translation for folding.

Abstract

Recent experiments and simulations have demonstrated that proteins can fold on the ribosome. However, the extent and generality of fitness effects resulting from co-translational folding remain open questions. Here we report a genome-wide analysis that uncovers evidence of evolutionary selection for co-translational folding. We describe a robust statistical approach to identify loci within genes that are both significantly enriched in slowly translated codons and evolutionarily conserved. Surprisingly, we find that domain boundaries can explain only a small fraction of these conserved loci. Instead, we propose that regions enriched in slowly translated codons are associated with co-translational folding intermediates, which may be smaller than a single domain. We show that the intermediates predicted by a native-centric model of co-translational folding account for the majority of these loci across more than 500 E. coli proteins. By making a direct connection to protein folding, this analysis provides strong evidence that many synonymous substitutions have been selected to optimize translation rates at specific locations within genes. More generally, our results indicate that kinetics, and not just thermodynamics, can significantly alter the efficiency of self-assembly in a biological context.

Tables

Table 1. Table S1: Rare codons in the E. coli genome. A codon is considered to be rare if the abundance-weighted relative usage, as defined in Eq.  S1 , is less than 10%. For comparison, the highest relative-usage codon for each amino acid is shown in the rightmost column.

Amino acid	Codon	Rel. usage	Max. rel. usage
G	GGA	0.0387	0.4474 (GGC)
	GGG	0.0811
I	ATA	0.0161	0.6062 (ATC)
L	CTA	0.0269	0.6785 (CTG)
	CTC	0.0880
	CTT	0.0853
	TTA	0.0495
	TTG	0.0718
P	CCC	0.0842	0.6344 (CCG)
R	AGA	0.0097	0.5283 (CGT)
	AGG	0.0043
	CGA	0.0170
	CGG	0.0260
S	AGT	0.0539	0.3198 (TCT)
	TCA	0.0547
	TCG	0.0657
T	ACA	0.0840	0.5224 (ACC)

Table 2. Table S2: The prokaryotic genomes used in the multiple-sequence alignments. The middle columns show the average codon and amino-acid identities of all genes relative to both the E. coli genome and the other species in the alignment. When computing codon and amino-acid identities, we have excluded genes whose lengths differ from the E. coli homolog by more than 20%, and we have ignored all insertions and deletions.

	Codon identity		Amino-acid identity
Genome ID	E. coli	Average	E. coli	Average	Species
NC_000913.3	1.000	0.390	1.000	0.682	Escherichia coli K-12 MG1655
NC_007606.1	0.775	0.376	0.830	0.654	Shigella dysenteriae
NZ_CP007557.1	0.473	0.382	0.803	0.687	Citrobacter freundii
NC_003198.1	0.448	0.397	0.767	0.681	Salmonella enterica
NC_014121.1	0.420	0.380	0.739	0.675	Enterobacter cloacae
NC_015663.1	0.407	0.384	0.726	0.675	Enterobacter aerogenes
NZ_CP007215.1	0.406	0.375	0.722	0.673	Enterobacter sacchari
NC_016845.1	0.405	0.384	0.724	0.671	Klebsiella pneumoniae
NC_009778.1	0.384	0.372	0.694	0.659	Cronobacter sakazakii
NZ_CP009450.1	0.382	0.369	0.702	0.658	Pluralibacter gergoviae
NZ_CP009451.1	0.373	0.359	0.692	0.655	Cedecea neteri
NC_017910.1	0.349	0.342	0.658	0.631	Shimwellia blattae
NZ_CP009454.1	0.306	0.318	0.591	0.598	Pantoea rwandensis
NC_016818.1	0.304	0.310	0.590	0.592	Rahnella aquatilis
NC_014500.1	0.301	0.313	0.582	0.586	Dickeya dadantii
NC_008800.1	0.284	0.292	0.594	0.604	Yersinia enterocolitica
NC_003143.1	0.273	0.283	0.578	0.592	Yersinia pestis
NC_022546.1	0.268	0.289	0.502	0.520	Plautia stali symbiont (Enterobacter)

Equations

$$p_{\text{use}}^s(c|a)=\frac{{\sum }_g{x}_{sg}{\sum }_{i=1}^{L_{sg}}1(c_{sgi}=c)}{{\sum }_g{x}_{sg}{\sum }_{i=1}^{L_{sg}}1(a_{sgi}=a)},$$

$${\lambda }_g\equiv M^{-1}\sum _{s=1}^M{L}_{sg}^{-1}\sum _{i=31}^{L_{sg}}1[p_{\text{use}}^{sgi}(c|a)\le 0.1].$$

$$n_{\text{enr}}^{sg}\equiv \min \left(l>0|1-p_{\text{Poisson}}^{}({\lambda }_g,l,n_s)\le 0.15\right),$$

$$p_{\text{rare}}^{s,\lambda }(a)\equiv \frac{{\sum }_{\{g{\}}_{\lambda }}{\sum }_{i=31}^{L_{sg}}1(a_{sgi}=a)1[p_{\text{use}}^{sgi}(c|a)\le 0.1]}{{\sum }_{\{g{\}}_{\lambda }}{\sum }_{i=31}^{L_{sg}}1(a_{sgi}=a)},$$

$$p_{\text{rare}}^{sg}(a)=\gamma p_{\text{rare}}^{s,{\lambda }_g}(a),$$

$$M^{-1}\sum _{s=1}^M{L}_{sg}^{-1}\sum _{i=31}^{L_{sg}}\frac{\gamma p_{\text{rare}}^{s,{\lambda }_g}(a_{sgi})}{1+(\gamma -1)p_{\text{rare}}^{s,{\lambda }_g}(a_{sgi})}={\lambda }_g.$$

$$p_{\text{enr}}^{sg}=\sum _{l=n_{\text{enr}}^{sg}}^{n_s}\sum _{\Omega }^{\text{perm.}}\prod _{i=1}^l{p}_{\text{rare}}^{sg}(a_{s{\Omega }_i})\prod _{j=l+1}^{n_s}[1-p_{\text{rare}}^{sg}(a_{s{\Omega }_j})],$$

$$p_{\text{neu}}^g(m_{\text{obs}},M)=\sum _{m=m_{\text{obs}}}^M\sum _{\sigma }^{\text{perm.}}\prod _{i=1}^m{p}_{\text{enr}}^{{\sigma }_{i}g}\prod _{j=m+1}^M(1-p_{\text{enr}}^{{\sigma }_{j}g}),$$

$${ϵ}_{uv}=-{\left({\alpha }_{\text{helix}}\right)}^{1_{uv}^{\text{helix}}}\left[n_{uv}^{\text{nc}}+{\alpha }_{\text{hb}}{1}_{uv}^{\text{hb}}\right],$$

$$\frac{f(g)}{k_{\text{B}}T}=\frac{1}{2}\sum _{u,v\in g}\frac{{ϵ}_{uv}}{k_{\text{B}}T}+\frac{\mu }{k_{\text{B}}T}(N_g-1)-\frac{\Delta S_{\text{l}}(g)}{k_{\text{B}}},$$

$$\frac{\Delta S_{\text{l}}(l)}{k_{\text{B}}}\equiv \{\begin{array}{ll}-\frac{|l|\mu }{k_{\text{B}}T}& \text{if }|l|\le b,\\ -\frac{b\mu }{k_{\text{B}}T}-\frac{3}{2}\left[\ln \frac{|l|}{b}+\frac{r(l{)}^2}{b^2|l|}\right]& \text{if }|l|>b,\end{array}$$

$$\frac{F_{\mathbf{r}}(L)}{k_{\text{B}}T}$$

$$\frac{F_{\text{mf}}(\{{\rho }_u\})}{k_{\text{B}}T}$$

$$=$$

$$q_u=1+\exp \left(-\tilde{\mu }-\sum _{v\in g}{\widetilde{ϵ}}_{uv}{\rho }_v\right).$$

$${\widetilde{ϵ}}_{uv}$$

$${\tilde{\mu }}_u$$

$$⟨\Delta S_{\text{l}}(l){⟩}_{\{{\rho }_u\}}=\frac{{\sum }_{\{(u,v){\}}_l}{\rho }_u{p}_{\text{l}}([u,v]){\rho }_v{e}^{\frac{\Delta S_{\text{l}}\left([u,v]\right)}{k_{\text{B}}}}{S}_{\text{l}}\left([u,v]\right)}{{\sum }_{\{(u,v){\}}_l}{\rho }_u{p}_{\text{l}}([u,v]){\rho }_v{e}^{\frac{\Delta S_{\text{l}}\left([u,v]\right)}{k_{\text{B}}}}},$$

$$\frac{F(L)}{k_{\text{B}}T}\simeq \underset{\mathbf{r}}{\min }\left[\sum _{(a,b)\in \mathbf{r}}\frac{F_{\text{mf}}(\{{\rho }_u{\}}_{ab})}{k_{\text{B}}T}\right],$$

Full text

Introduction

Many proteins can begin folding to their native states before their synthesis is complete Komar (2009); Pechmann et al. (2013). As much as one-third of a bacterial proteome is believed to fold co-translationally Ciryam et al. (2013), with an even higher percentage likely in more slowly translated eukaryotic proteomes. Numerous experiments on both natural and engineered amino-acid sequences have shown that folding during synthesis can have profound effects: compared to denatured and refolded chains, co-translationally folded proteins may be less prone to misfolding Netzer and Hartl (1997); Frydman et al. (1999); Komar et al. (1999); Kim et al. (2015); Siller et al. (2010); Ugrinov and Clark (2010); Agashe et al. (2013); Clark and King (2001), aggregation Evans et al. (2008) and degradation Zhang et al. (2009), or they may preferentially adopt alternate stable structures Sander et al. (2014); Buhr et al. (2016); Zhou et al. (2013). Because the timescales for protein synthesis and folding are often similar O’Brien et al. (2014); Nissley et al. (2016), it is clear that the rate of translation can be used to tune the self-assembly of peptide chains in vivo Xu et al. (2013); Kimchi-Sarfaty et al. (2007). To this point, however, there exists little evidence that evolution has selected specifically for efficient co-translational folding kinetics across any substantial fraction of an organism’s proteome.

In this work, we provide evidence that evolutionary selection has tuned protein-translation rates to optimize co-translational folding pathways. Our approach is motivated by the hypothesis that pauses during protein synthesis may be beneficial for promoting the formation of native structure. By increasing the separation between the timescales for folding and translation, such pauses may promote the assembly of on-pathway intermediates, which, in turn, template the growth of further native structure. Many experimental and computational studies have shown that protein-folding naturally proceeds in a step-wise manner via structurally distinct intermediates Hartl and Hayer-Hartl (2009); Jacobs and Shakhnovich (2016), and that cooperative folding cannot commence until a minimal number of residues have emerged from the ribosome exit tunnel Eichmann et al. (2010); Holtkamp et al. (2015); Elcock (2006); O’Brien et al. (2012). These general findings suggest that any beneficial pauses during synthesis should occur at specific locations within an amino-acid sequence.

Using a coarse-grained model of co-translational folding, we find that translational pauses tend to be associated with stable, native-like co-translational folding intermediates. The relevant folding intermediates are typically not complete structural domains, as has often been assumed Jacobson and Clark (2016), and may be distinct from intermediates that are observed when refolding from a denatured ensemble. By comparing putative translational pause sites with a neutral model that accounts for gene-specific codon usage, we show that evolutionarily optimized co-translational folding is a widespread feature of the E. coli genome. Our results therefore highlight the extent to which evolution has tuned the self-assembly pathways, and not just the native structures, of complex biomolecules.

Results

Unbiased identification of slowly translated regions

Our analysis of beneficial pauses in protein synthesis relies on the identification of regions within mRNA transcripts that are enriched in ‘rare’ codons (SI Appendix, Table S1), i.e. codons that are used substantially less often than alternate synonymous codons in highly expressed genes Sharp and Li (1987). Despite numerous attempts to predict codon-specific translation rates based on physical factors Dana and Tuller (2014); Li et al. (2012); Pelechano et al. (2015); Weinberg et al. (2016), such as tRNA concentrations, translation-speed estimates based on relative-usage metrics Sharp and Li (1987); Clarke IV and Clark (2008) remain among the most accurate Yu et al. (2015); Chaney and Clark (2015); Spencer et al. (2012). Thus, using codon rarity as a proxy for translation speed, we can look for pauses in synthesis by identifying regions in a mRNA transcript that are locally enriched in rare codons.

However, an appropriate neutral model must account for two potential sources of synonymous codon-usage bias at the level of an individual gene. First, we controlled for the overall rare-codon usage in a gene, which is defined as the fraction of rare codons in the entire transcript (SI Appendix, Fig. S1). Multiple factors have been hypothesized to contribute to the overall degree of codon adaptation of each gene, including evolutionary selection for rapid synthesis, accurate translation and the stability of mRNA transcripts Plotkin and Kudla (2011). By taking a gene’s average codon usage into account, we instead pick out regions that are locally enriched in rare codons relative to the gene-specific background. Second, we accounted for synonymous-codon bias due to the amino-acid composition of the protein sequence. Assuming that amino-acid sequences are under stronger selection pressure and can thus be considered immutable, we estimated the average rare-codon frequencies for each amino-acid type among all genes with a similar level of rare-codon usage. Having controlled for the overall rare-codon usage and the amino-acid sequence, we modeled neutral codon usage as a Bernoulli process with sequence-dependent rare-codon probabilities (see SI Appendix, Sec. 1A).

Evaluation of evolutionary conservation

Next, we assessed the functional importance of local rare-codon enrichment by looking for conservation of rare-codon usage across multiple-sequence alignments (Figure 1). We extended the neutral model described above to 18 sufficiently diverged prokaryotic genomes, with rare-codon definitions and gene-specific rare-codon probabilities computed for each genome independently. Here our approach differs from conventional conservation analyses, because we are interested in the enrichment of rare codons within contiguous 15-codon segments of a transcript, as opposed to the codon usage at each aligned site Pechmann and Frydman (2013); Widmann et al. (2008). By examining conservation of rare-codon enrichment, we can identify local regions that do not align precisely but nevertheless result in translational pauses at similar places within the protein sequence. This approach also allows for a meaningful comparison of the local rare-codon enrichment in sequence alignments that contain insertions and deletions. Our choice of a 15-codon enrichment region is comparable to the length of a typical element of protein secondary structure, and we verified that regions with widths of 10 and 20 codons yield similar results. In contrast, larger enrichment regions defined on the basis of complete domains rarely differ significantly from the background rare-codon usage, while analyses of single aligned sites tend not to produce statistically significant results.

To be relevant for co-translational folding, putative slowly translated regions must meet two criteria: a high degree of conservation of slowly translated codons, and a low probability of such an occurrence in the neutral model. For a region to be considered both enriched and conserved, we required that the local concentration of rare codons deviate from the background distribution by approximately one standard deviation in at least 75% of the sequences in the alignment; Fig. S2 of the SI Appendix shows that our results are robust with respect to this conservation threshold. We then computed an associated p-value that reports the probability, within the neutral model, of randomly generating at least the observed number of enriched regions from reverse translations, i.e., by sampling synonymous sequences using the aligned amino-acid sequences and a probabilistic model of the codon usage for each amino-acid type (see SI Appendix, Sec. 1A). This second criterion is central to our findings, as we shall discuss below. We emphasize that these criteria are distinct: depending on the amino-acid identities, it is possible to observe low p-values without significant rare-codon enrichment relative to the background, and vice versa. Consequently, both criteria must be satisfied to constitute evidence for evolutionary selection.

Our analysis reveals numerous rare-codon enrichment loci in the E. coli genome that are inconsistent with the neutral model, and are thus likely to be a result of evolutionary selection (SI Appendix, Fig. S3a). Although these regions occur throughout the mRNA transcripts, their locations are biased towards both the 5’ and 3’-ends (SI Appendix, Fig. S3b). While these trends have been noted previously Clarke IV and Clark (2010), our analysis confirms that the increased probability of rare-codon enrichment at the 3’-end is evolutionarily conserved and is not a consequence of the amino-acid sequences. Furthermore, we find that these biases become more pronounced as we lower the p-value threshold used for comparison with the neutral model (SI Appendix, Fig. S3b), suggesting that any false positives from our analysis are relatively evenly distributed throughout the transcripts. We also analyzed the codon-level similarity among the genomes in our alignments and verified that these results reflect conservation of rarity as opposed to conservation of specific codons (SI Appendix, Fig. S4).

Comparison with predicted co-translational folding pathways

To probe the potential consequences of local rare-codon enrichment for protein folding, we next examined the formation of native-like intermediates during protein synthesis. A large body of simulation evidence Trovato and O’Brien (2016) has shown that intermediates must be stable at equilibrium in order to be sampled with high probability during co-translational folding and are only likely to form when the folding rate is fast relative to the protein elongation rate. Therefore, while an intermediate’s equilibrium free energy does not completely determine whether it will appear on a co-translational folding pathway, we assume that stability at equilibrium is necessary for a pause in translation to promote the development of native structure.

Here we applied a coarse-grained model Jacobs and Shakhnovich (2016) to predict the formation of stable partial structures during nascent-chain elongation. Importantly, this model captures the tertiary structure of nascent chains and does not assume that domains fold cooperatively or independently. To model co-translational folding, a nascent chain of length $L$ is allowed to form native contacts among the first $L$ residues of the full protein. We then computed the minimum free energy of a nascent chain, relative to an unfolded ensemble, using a mean-field theory based on the protein’s native structure (see SI Appendix, Sec. 1B). This approach captures the opposing contributions to the free energy from energetically favorable native contacts and the configurational entropy of an unfolded chain. We used a native-centric energy function that emphasizes hydrogen bonds and contacts between larger residues Jacobs and Shakhnovich (2016), while the thermodynamic stability of the native state is fixed based on the full protein length Ghosh and Dill (2009). We show in Fig. S5 of the SI Appendix that tuning the native-state stability does not significantly affect the results of our analysis.

Our calculations predict that, in general, native structure forms discontinuously during nascent-chain elongation. In the example shown in the lower panel of Figure 2, decreases in the nascent-chain free energy occur at distinct chain lengths. These sudden drops correspond to the appearance of stable intermediates with native-like tertiary structure. In contrast, at chain lengths corresponding to the intervening plateaus, the nascent-chain free energy remains constant because the newly synthesized residues cannot form sufficient stabilizing contacts with any existing tertiary structure. Unsurprisingly, the probability of finding a stable on-pathway intermediate increases as synthesis nears completion (SI Appendix, Fig. S6).

We are now in a position to test the relationship between translational pausing and the formation of native-like intermediates. The ribosome exit tunnel is widely believed to conceal between 30 and 40 amino acids Gloge et al. (2014); Chaney and Clark (2015), although a greater number may be accommodated in partially helical conformations Lu and Deutsch (2005). In addition, some tertiary structure formation may commence within the exit-tunnel vestibule Nilsson et al. (2015). A beneficial pause in synthesis should therefore be separated from a co-translational intermediate by a distance that is roughly equivalent to the exit-tunnel length (see Methods). An example of this correspondence is shown in Figure 2, where a putative translational pause is located approximately 30 residues downstream of the formation of a predicted intermediate. However, we emphasize that, according to the present hypothesis, the formation of an intermediate is necessary but not sufficient to expect that a translational pause would be beneficial. For example, intermediates that fold quickly relative to the average translation rate or appear less than the exit-tunnel distance from the end of the protein are unlikely to be accompanied by a productive pause.

Conserved, enriched regions associate with predicted co-translational folding intermediates

By applying this analysis to a set of approximately 500 E. coli proteins with known native structures, we find widespread support for our co-translational folding hypothesis. In particular, we find that the co-translational folding intermediates predicted by our coarse-grained model account for a significant proportion ($\gtrsim 50\%$) of the putative slowly translated regions (Figure 3). Most importantly, we find that the fraction of rare-codon-enriched regions that can be explained by our model increases consistently as we reduce the p-value threshold for establishing evolutionary conservation. In other words, the predictive power of our model improves as false positives related to the random clustering of rare codons are preferentially eliminated. This trend is also robust with respect to variations in the precise definition of codon rarity (SI Appendix, Fig. S7).

We further tested the sensitivity of our co-translational folding predictions by repeating the above analysis with randomized control sequences, which preserve the total number of pause sites at each p-value threshold but uniformly distribute their locations across the transcripts (Figure 3; see Methods). Although a significant fraction ($\sim\!35$%) of the fictitious pause sites in the randomized sequences can also be explained by our model, likely due to chance overlaps with predicted intermediates, the difference between the genomic and randomized data increases markedly at lower p-value thresholds (one-sided ${p<10^{-7}}$ at neutral-model p-value thresholds below 0.01; see SI Appendix, Fig. S8a). Two alternative controls (SI Appendix, Fig. S9), in which the randomized pause sites are drawn from a non-uniform distribution with a 3’-end bias or obtained directly from reverse translations, verify that our results are not solely a consequence of the 3’-end rare-codon bias in the mRNA transcripts or the amino-acid sequences of the proteins.

Next, we performed inverse tests to assess whether co-translational folding intermediates are preferentially associated with putative translational pauses. However, because the formation of an intermediate is not in itself a sufficient condition for a translational pause to be beneficial, we find that the overall frequency of such associations is small relative to the number of predicted intermediates (SI Appendix, Fig. S10). We therefore computed the odds ratio of finding conserved, rare-codon-enriched regions just downstream of a predicted intermediate, as opposed to elsewhere in a mRNA transcript. The results shown in Figure 4 confirm that the association between folding intermediates and translational pause sites is highly significant (one-sided ${p<10^{-7}}$ at neutral-model p-value thresholds below 0.01; see SI Appendix, Fig. S8b) and, importantly, is not related to the overall frequency of predicted co-translational intermediates. Here again, the predictive power of our model shows a strong dependence on the p-value threshold used for screening putative pause sites. In contrast, tests with randomized control sequences do not deviate from an odds ratio of unity.

We also applied our analysis to structural domain boundaries, which have previously been suggested to play a role in coordinating co-translational folding Purvis et al. (1987). Nevertheless, in agreement with more recent works Jacobson and Clark (2016), we find little evidence of selection for translational pausing at domain boundaries. For these comparisons, we used domain definitions for approximately 800 E. coli proteins from the SCOP database Murzin et al. (1995). Figure 3 shows that domain boundaries explain a much smaller fraction ($\lesssim 12\%$) of the putative pause sites than our folding model. Furthermore, the predictive power of the domain-boundary hypothesis does not vary with the p-value threshold, and the odds ratios are nearly indistinguishable from the randomized controls (Figure 4). These conclusions also hold for various related hypotheses: instead of assuming that a domain must be completely synthesized before folding, we tested models where native structure begins to form either at a fixed number of residues prior to the domain boundary or at a fixed percentage of the domain length (SI Appendix, Fig. S11). In all cases, the correspondence between the domain boundaries and the conserved, rare-codon enriched loci is significantly weaker than the results of our co-translational folding model. While these findings do not imply that domain boundaries are irrelevant for co-translational folding, we can conclude that the domain-boundary hypothesis is insufficient to explain the vast majority of conserved, slowly translated regions.

Discussion

By integrating a multiple-sequence analysis of synonymous codon conservation with protein-folding theory, we have shown that highly conserved rare-codon clusters preferentially associate with predicted co-translational folding intermediates. The putative pause sites in the E. coli genome that are both evolutionarily conserved and unaccounted for by the neutral model systematically appear downstream of predicted co-translational folding intermediates at distances that are similar to the length of the ribosome exit tunnel. Our large-scale study therefore supports the hypothesis that beneficial pauses during protein synthesis follow key steps in the assembly of native structure. Comparisons with randomized control sequences confirm that our observations are highly significant.

This analysis of co-translational folding pathways, as opposed to elements of the static native structure, provides new insights into the interplay between translation and the self-assembly kinetics of nascent proteins. The stabilization of a partial structure often occurs well before a native domain is completely synthesized, especially in cases where the domain comprises more than 200 residues. In particular, sub-domain co-translational folding intermediates typically appear when sufficient tertiary contacts are available to compensate for the loss of chain entropy that is required for folding. Overall, a relatively small fraction ($\lesssim 15\%$) of all predicted intermediates are followed by conserved translational pauses (SI Appendix, Fig. S10), but the association between folding intermediates and conserved pauses is highly significant (Figure 4). This observation is consistent with our hypothesis, since the effect of a pause depends on the relative timescales for translation and folding, as well as potential interference due to non-native interactions. In addition, this observation explains why pause sites are not preferentially associated with domain boundaries: although fully synthesized domains may be stable on the ribosome, the prior formation of a partial-chain intermediate is likely to affect the subsequent folding rates for other parts of the protein. As a result, the entire co-translational folding pathway must be considered when interpreting the effect of a pause in translation. We anticipate that an optimal translation protocol Sharma et al. (2016) could be predicted with knowledge of the substructure-specific folding and translation rates, as well as their propensities for forming non-native interactions, including interactions with the surface of the ribosome Nilsson et al. (2017). In addition, an optimal translation protocol is likely to be affected by the presence of misfolded intermediates, which may be avoided by increasing the local translation rate O’Brien et al. (2014); Trovato and O’Brien (2017).

The approach that we have taken in this work improves upon earlier studies of rare-codon usage, which have addressed alternative hypotheses regarding translational pausing but yielded mixed results Chaney and Clark (2015); Thanaraj and Argos (1996); Pechmann and Frydman (2013); Lee et al. (2010); Chartier et al. (2012); Brunak and Engelbrecht (1996); Zhou et al. (2009). In addition to our distinct focus on co-translational folding pathways, our conclusions are more robust due to our use of a multiple-sequence analysis to detect evolutionary conservation, as well as our formulation of a neutral model that controls for both amino-acid composition and the inherent codon-usage variability across genes. The statistical significance of our results is further increased by the much larger sample size used here.

While this paper was under review, we became aware of a contemporaneous study Chaney et al. (2017) that identifies conserved rare-codon clusters via a complementary statistical analysis. The authors also observe extensive rare-codon conservation across mRNA transcripts and similarly find no evidence of enrichment near domain boundaries.

Synonymous substitutions can also affect protein synthesis through mechanisms that are unrelated to protein folding, most notably via changes to mRNA secondary structure and stability Plotkin and Kudla (2011). However, many experimental studies have shown that these effects originate predominantly from substitutions near the 5’-end of the mRNA transcripts, and typically modulate the total protein production as opposed to the protein quality Kudla et al. (2009); Goodman et al. (2013). Such mRNA-specific effects are thus a likely explanation for the observed 5’-end bias in rare-codon enrichment, where variations in translation speed are unlikely to play a role in co-translational folding. Consequently, we have excluded N-terminal rare codons from our analysis. In addition to rare-codon usage, various studies have proposed that additional factors, such as interactions between the nascent chain and the ribosome exit tunnel Chaney and Clark (2015) or the presence of internal Shine–Dalgarno motifs Li et al. (2012), can affect translation rates. It is likely that a more complete picture of sequence-dependent translation kinetics will enable further refinements to the co-translational folding model presented here.

In conclusion, our study highlights the importance of optimal kinetic pathways for efficient biomolecular self-assembly. Although a protein’s amino-acid sequence entirely determines its thermodynamically stable structure, it is becoming increasingly clear that synonymous mutations are not always silent. Our analysis provides strong evidence that evolutionary selection has tuned local translation rates to improve the efficiency of co-translational protein folding. Further work is needed to understand the relationship between genome-wide codon usage and translation rates and to improve the prediction of co-translational folding intermediates, including those that contain significant amounts of non-native structure. Nevertheless, our results indicate that folding kinetics play a role in evolutionary selection and suggest that similar relationships may exist for other biological self-assembly phenomena, such as the assembly of macromolecular complexes.

Methods

We constructed alignments based on the amino-acid sequences of homologous genes from 18 prokaryotic species with between 50 and 85% average amino-acid sequence identity to E. coli (SI Appendix, Table S2). We then computed p-values associated with rare-codon enriched regions, assuming biased reverse translations and a gene-specific model for the probability of each amino-acid type being encoded by a rare codon. Consensus crystal structures were constructed for 511 non-membrane E. coli proteins with 500 residues or fewer using Protein Databank Berman et al. (2000) entries containing complete structures for sequences with at least 95% amino-acid identity to the E. coli gene. SCOP domain assignments were obtained from Murzin et al. (1995) for all proteins with at most 500 residues. Due to the uncertainty in the number of amino acids that are concealed in the ribosome exit tunnel and the potential for steric interactions between folding intermediates and the ribosome, we consider a rare-codon-enriched region to be associated with a folding intermediate if the enriched region is anywhere between 20 and 60 codons downstream from the position at which an intermediate first becomes stable, ignoring enriched regions within the first 80 codons of a transcript. An intermediate is identified whenever the monotonic co-translational free-energy profile decreases by more than ${1~{}k_{\text{B}}T}$ relative to the previous free-energy plateau; for example, see the pattern of alternating plateaus and precipitous free-energy decreases in the lower panel of Figure 2. To generate the randomized control sequences from which the control distributions in Figures 3 and 4 were calculated, we sampled locations for fictitious rare-codon-enriched regions from a uniform distribution over each mRNA transcript, excluding the first 80 codons. This uniform distribution was normalized such that the expected number of fictitious enriched regions is equal to the total number of observed enriched regions at each p-value threshold. See SI Appendix, Sec. S1 for complete details of all methodologies. Essential data are provided in the SI Materials. All code necessary to reproduce these results is available at https://faculty.chemistry.harvard.edu/shakhnovich/software.

Acknowledgements.

We thank Sanchari Bhattacharyya and Michael Manhart for many insightful discussions. This work was supported by NIH grants R01GM124044 and F32GM116231.

Supplementary Information

S1 Extended methods

S1.1 Neutral model and statistics of rare-codon enrichment

In this section, we describe a neutral model of gene-specific synonymous codon usage. For each genome, we define a set of rare codons that are used significantly less frequently than alternative synonymous codons in the most highly expressed genes Sharp and Li (1987). Using experimentally determined protein abundances to account for the differing expression levels among genes, we compute the relative synonymous codon frequencies for each amino-acid type,

[TABLE]

where $c$ is a codon for an amino acid $a$, $x_{sg}$ is the protein abundance of gene $g$ in genome $s$, $i$ is an index that runs over all coding positions up to the protein length $L_{sg}$, and ${\bm{1}(\cdot)}$ is the indicator function. We then define the set of rare codons as those codons whose protein-abundance-weighted relative usage, ${p^{s}_{\text{use}}(c|a)}$, is less than 10%. The relative usages of E. coli rare codons, as determined by Eq. S1, are shown in Table S1. We use a composite data set for the protein abundances in E. coli Wang et al. (2012) and assume that these abundances are similar for all prokaryotic genomes in our alignment (Table S2). The rare-codon definitions turn out to be comparable, but not identical, for these genomes.

Next, we address the absolute enrichment of rare codons in a local region of a gene, without considering the amino-acid sequence. We calculate the average fraction of rare codons in the $M$ aligned sequences, which all have at least 50% amino-acid identity with respect to the E. coli gene and differ in length by no more than 20%,

[TABLE]

Here and below, we exclude the first 30 codons to avoid the 5’-end bias, and we define ${p^{sgi}_{\text{use}}(c|a)\equiv p^{s}_{\text{use}}(c_{sgi}|a_{sgi})}$ for brevity. We then define a local-region enrichment threshold $n^{sg}_{\text{enr}}$ for each aligned sequence by considering a Poisson process with rare-codon-usage rate $\lambda_{g}$,

[TABLE]

where ${p^{\phantom{\alpha}}_{\text{Poisson}}(\lambda,l,n)}$ is the cumulative distribution function of $l$ events occurring in $n$ trials given a rate $\lambda$. The cumulative probability $0.15$ corresponds to an approximately one-standard-deviation fluctuation in the local rare-codon usage. Every local region has a fixed width within the E. coli sequence, in which case ${n_{s}=n}$. (We have chosen ${n=15}$ in the main text.) However, due to insertions and deletions, the number of codons that align to a given $n$-codon region of the E. coli sequence may be different for the other sequences in the alignment. Therefore, when applying Eq. S3 to these other sequences, we use the number of codons $n_{s}$ that align to each $n$-codon region of the E. coli sequence, including any insertions. With this approach, we can identify structurally equivalent regions that are enriched in rare codons in each aligned sequence.

We can now calculate a p-value that accounts for biases in rare-codon usage due to the local amino-acid composition. To do so, we estimate the probability of using a rare codon for each amino-acid type in a given gene by analyzing the relative codon-usage bias of all genes with a similar overall rare-codon usage $\lambda$,

[TABLE]

where ${\{g\}_{\lambda}}$ is the set of genes with similar $\lambda$. In practice, we implement Eq. S4 by sorting all genes according to their overall rare-codon usages and splitting them into ten groups, each comprising approximately 400 genes. However, because the average codon usage is not exactly the same for all genes in each of these groups, we then adjust the rare-codon probabilities slightly to match $\lambda_{g}$,

[TABLE]

where the scaling factor $\gamma$ is chosen such that

[TABLE]

Eq. S6 is nearly linear in $\gamma$ while ensuring that all probabilities are positive.

Finally, the neutral-model probability that a local region contains at least $n^{sg}_{\text{enr}}$ rare out of $n_{s}$ codons is

[TABLE]

where ${\{\Omega\}}$ is the set of all permutations of the codon positions within the given region and $\Omega_{i}$ is the position at index $i$ in the permutation $\Omega$. The probability of observing at least $m_{\text{obs}}$ enriched sequences is then

[TABLE]

where ${\{\sigma\}}$ is the set of all orderings of the $M$ sequences, and $\sigma_{i}$ is the sequence at the $i$th index of the ordering $\sigma$. We require ${m_{\text{obs}}/M\geq 0.75}$ for an enriched region to be considered conserved (see Figure S2). Eq. S8 is the neutral-model p-value used throughout this work.

S1.2 Quasi-equilibrium co-translational folding model

In this section, we adapt the native-centric coarse-grained model described in Ref. Jacobs and Shakhnovich, 2016 to predict free-energy landscapes for co-translational folding (Supplementary Methods Figure 1a). A native contact between a pair of residues is formed when at least one pair of heavy atoms is separated by at most 4 Å in the consensus crystal structure and the residues are at least three positions apart in the primary sequence. The effective potential energy between residues $u$ and $v$ depends on the number of such heavy-atom contacts, ${n^{\text{nc}}_{uv}\geq 1}$; whether the pair of residues forms a hydrogen bond, $\bm{1}^{\text{hb}}_{uv}$; and whether the contact is part of an $\alpha$-helix, $\bm{1}^{\text{helix}}_{uv}$:

[TABLE]

with ${\alpha_{\text{helix}}=5/8}$ and ${\alpha_{\text{hb}}=16}$.

In this discrete model, a microscopic configuration is defined by the set of native contacts that are formed. The free energy of such a configuration is determined by opposing energetic and entropic contributions: native contacts are stabilizing according to Eq. S9, and configurational entropy is gained when native contacts are broken. Furthermore, sets of native contacts can be decomposed into independent structured regions (Supplementary Methods Figure 1b). Within each structured region, all possible native contacts are formed among the participating residues; however, different structured regions are linked by chains of disordered residues and can thus move independently. We therefore evaluate the free energy of each structured region separately. The free energy of a configuration $g$ with a single structured region is

[TABLE]

where $u$ and $v$ are residues that form native contacts in the configuration $g$, $N_{g}$ is the number of contact-forming residues, $k_{\text{B}}$ is the Boltzmann constant, $T$ is the absolute temperature, and ${\mu=2k_{\text{B}}T}$ is the entropic penalty associated with native-contact formation Jacobs and Shakhnovich (2016). Sequences of disordered residues that begin and end in the same structured region are referred to as loops (Supplementary Methods Figure 1b–c). The total loop entropy, ${\Delta S_{\text{l}}(g)}$, is the sum of ${\Delta S_{\text{l}}(l)}$ for all loops ${\{l\}_{g}}$ that terminate in the structured region,

[TABLE]

where ${|l|}$ is the number of disordered residues in the loop, ${r(l)}$ is the distance between the loop endpoints, and we assume Gaussian statistics for all loops longer than the Kuhn length ${b=2}$ residues.

Our goal is to compute the free energy of a nascent chain of length $L$, ${F(L)}$. First, we note that the free energy of a polypeptide chain in this model is equal to the sum of the free energies of one or more structured regions, which do not interact with one another. To minimize the free energy ${F(L)}$, any combination of non-overlapping structured regions using the first $L$ residues is permitted. For example, various possible sets of structured regions are illustrated in Supplementary Methods Figure 1c: native contacts are formed within each of the structured regions in a configuration, but no contacts are formed between residues in different structured regions. The set of allowed sets of structured regions can be written as ${\bm{r}\!\equiv\!\{(a_{1},b_{1}),\ldots,\!(a_{k},b_{k})|\,a_{1}\!\geq\!0,a_{i}\!\geq\!b_{i-1}\forall i\!\in\![2,k],b_{k}\!\leq\!L\}}$, where the pair ${(a_{i},b_{i})}$ indicates the first and last residues that can participate in the structured region with index $i$, and the number of structured regions is ${k\geq 0}$. Then, to account for the fact that many distinct microstates are typically consistent with a given structured region, we write the free energy of a nascent chain with a specific set of structured regions $\bm{r}$ as

[TABLE]

where ${\{g\}_{ab}}$ indicates the set of all microstates that have a single structured region incorporating residues in the range ${[a,b]}$, and ${F_{ab}}$ is the free energy of this ensemble. Unfortunately, directly calculating ${F_{ab}}$ requires a Monte Carlo approach Jacobs and Shakhnovich (2016), which can be computationally intensive.

We therefore introduce a mean-field approximation to compute the free energy of a single structured region, ${F_{ab}}$. We assume that all configurations in this ensemble contain a single structured region with an essentially equivalent set of loops (i.e., the loops grow or shrink by adding or removing adjacent residues, but otherwise the structured region remains unchanged; see Supplementary Methods Figure 1d). Instead of considering the set of discrete configurations ${\{g\}}$, we now characterize the ensemble by the average ordering of each residue in this set of microscopic configurations, ${\{\rho_{u}\}}$. Ignoring correlations in the residues beyond nearest-neighbor native contacts, we can write the mean-field free energy of this ensemble of microscopic configurations as the sum of the average energetic and entropic contributions,

[TABLE]

where $q_{u}$ is the single-residue mean-field partition function

[TABLE]

Native contacts are correlated on a length scale comparable to the Kuhn length due to the finite size of the residues and the chain connectivity Jacobs and Shakhnovich (2016). We have excluded native contacts between residues that are fewer than three positions apart in the sequence when constructing the energy matrix ${\{\epsilon_{uv}\}}$, because such contacts are likely to be present in both the denatured and folded ensembles. Nevertheless, we need to account for local geometric correlations in native contact formation. To impose correlations between residues that are very close in the sequence, we add nearest and next-nearest-neighbor couplings to ${\{\epsilon_{uv}\}}$ and ${\{\mu_{u}\}}$,

[TABLE]

where $\delta$ is the Kronecker delta, ${\bm{1}_{r}(v)}$ indicates whether the residue $v$ can form native contacts in the structured region, and we assume that ${u<v}$ and ${\bm{1}_{r}(u)=1}$ for notational simplicity. We have chosen coupling constants ${\lambda_{1}=0.8}$ and ${\lambda_{2}=0.4}$ to obtain comparable results to the Monte Carlo simulations reported in Ref. Jacobs and Shakhnovich, 2016. Choosing ${\lambda_{1},\lambda_{2}>0}$ tends to increase the average ordering at the free-energy minimum for all residues that form multiple contacts in the structured region.

To determine ${F_{\text{mf}}}$ and ${\{\rho_{u}\}}$, we first minimize the mean-field free energy without considering the loop term, ${\Delta S_{\text{l}}}$. Eq. S1.2 can be minimized numerically using a standard non-linear, multi-dimensional minimization algorithm. The probability of forming a native contact between residues $u$ and $v$ in the mean-field model is then ${p^{\text{contact}}_{uv}=\rho_{u}\bm{1}(\epsilon_{uv}<0)\rho_{v}}$; example contact maps showing ${\{p^{\text{contact}}_{uv}\}}$ at the free-energy minimum can be found in Supplementary Methods Figure 2. While it is not guaranteed that Eq. S1.2 will have a single minimum, in practice we find that, excluding the trivial solution, ${\{\rho_{u}\}=0}$, this is almost always the case. Such behavior is expected from our previous study Jacobs and Shakhnovich (2016) of the model presented in Eq. S10.

Next, we compute the ensemble-averaged loop entropy, ${\langle\Delta S_{\text{l}}\rangle_{\{\rho_{u}\}}}$, perturbatively by assuming that this contribution to Eq. S1.2 does not alter ${\{\rho_{u}\}}$ at the free-energy minimum. To allow disordered loops to grow or shrink by adding or removing adjacent residues, we calculate the Boltzmann-weighted average loop entropy for all possible pairs of terminal residues ${\{(u,v)\}_{l}}$ for each loop $l$,

[TABLE]

where ${p_{\text{l}}([u,\!v])\equiv\prod_{u>t>v}(1-\rho_{t})}$, and ${\Delta S_{\text{l}}\left([u,\!v]\right)}$ implies Eq. S11 computed with the loop endpoints ${(u,v)}$.

Finally, we can calculate the minimum free energy of a nascent chain in this mean-field approximation. We consider native contacts among the first $L$ residues (Supplementary Methods Figure 2a). We then minimize the mean-field free energy by finding the lowest-free-energy set of structured regions,

[TABLE]

where ${\{\rho_{u}\}_{ab}}$ implies that ${\rho_{u}=0}$ for all ${u<a}$ and ${u>b}$. At each chain length, ${F(L)}$ reports the minimum free energy relative to a random coil of length $L$ (Supplementary Methods Figure 2b). The free energy defined in Eq. S19 is used for identifying co-translational folding intermediates as described in the main text.

S1.3 Analysis of co-translational folding pathways

Due to the uncertainty in the number of amino acids that are concealed in the ribosome exit tunnel and the potential for steric interactions between folding intermediates and the ribosome, we consider a rare-codon-enriched region to be associated with a folding intermediate if the enriched region is anywhere between 20 and 60 codons downstream from the position at which an intermediate first becomes stable. We ignore enriched regions within the first 80 codons of a transcript, since the conservation of rare codons near the 5’-end of the mRNA transcripts is believed to originate from other factors, such as efficient translation initiation. A co-translational folding intermediate is identified whenever the monotonic co-translational free-energy profile decreases by more than ${1~{}k_{\text{B}}T}$ relative to the previous free-energy plateau; for example, see the pattern of alternating plateaus and precipitous free-energy decreases in the lower panel of Figure 2. In Figure 3, the ‘predicted fraction’ is the number of rare-codon-enriched regions that are preceded by a folding intermediate, $n_{\text{RI}}$, divided by the total number of enriched regions, $n_{\text{R}}$. To compute the odds ratios in Figure 4, we first determined the fraction, $f_{\text{I}}$, of all codons in our data set that are 20 to 60 positions downstream of a predicted intermediate. The ‘presence’ odds ratio is then defined as ${(n_{\text{RI}}/n_{\text{R}})/f_{\text{I}}}$, and the ‘absence’ odds ratio is ${(1-n_{\text{RI}}/n_{\text{R}})/(1-f_{\text{I}})}$.

To generate the randomized control sequences from which the control distributions in Figures 3 and 4 were calculated, we sampled locations for fictitious rare-codon-enriched regions from a uniform distribution over each mRNA transcript, excluding the first 80 codons. This uniform distribution was normalized such that the expected number of fictitious enriched regions is equal to the total number of observed enriched regions at each p-value threshold.

S2 Supplementary tables and figures