Effect of Protein Structure on Evolution of Cotranslational Folding

Model connecting protein translation and folding to cellular fitness

To simulate evolution of lattice protein sequences, we built a model relating outcomes of protein translation and folding to cellular fitness. A protein undergoing translation can reach its native state during or after translation (Fig. 1 A). After translation, a free protein not in its native state is vulnerable to degradation or aggregation with other proteins (29), processes assumed here to be irreversible; aggregation is effectively irreversible if disaggregation is slow (30). We consider the survival probability of a protein in the cellular environment over time. Let $S\left(t\right)$ be the probability that a protein survives to time $t$, with $t=0$ being the time the protein leaves the ribosome and $S\left(0\right)=1$. For a protein not in its native state, there is some effective rate $k_d$ for the protein to be degraded or to aggregate. Based on this model, $S\left(t\right)$ evolves according to the following differential equation:

where $\theta \left(t\right)=\left\{0,1\right\}$ is an indicator function whose value is 1 when the protein is in its native state at time $t$. Under this model, $S\left(t\right)$ does not decrease if the protein is in its native state. The half-lives of thermodynamically unstable proteins are as short as a few minutes (31, 32, 33), compared with hours or days for stable proteins (31,33,34), justifying this assumption. A protein thus has a higher survival probability if it folds quickly and remains stably folded.

We next model how $S\left(t\right)$ affects the activity of the protein over time. We assume a protein can perform its biological function only when in its native state and if it has not been degraded. Based on these assumptions, the total cumulative activity of the protein (e.g., enzymatic output), $A_{total}$, is described probabilistically by the following equation:

where $k_a$ is an activity rate constant (which we set to 1), and $T$ is some long timescale corresponding to the period of time that the protein is biologically relevant, such as the length of the cell cycle.

Fig. 1 B illustrates the relationship between protein folding, survival, and activity for a single protein folding trajectory. In this particular trajectory, the protein folds posttranslationally. During translation, the protein is not vulnerable to degradation, and $S\left(t\right)$ is 1. After release from the ribosome, $S\left(t\right)$ decreases for every time unit the protein is not in the native state. $S\left(t\right)$ decreases at a higher rate during the initial passage to the folded state and then decreases more slowly after the protein enters the native state energy basin and fluctuates in and out of the native state. Protein activity only begins to accumulate after the protein reaches the native state. $A_{total}$ corresponds to cumulative activity at a later time T.

To reduce the amount of computation required for evaluating the protein activity function (Eq. 4), we assume proteins fluctuate on fast timescales in and out of the native state, occupying the native state with probability $P_{\mathit{nat}}=⟨\theta \left(t\right)⟩$. After simplifications and applying this fast-fluctuation assumption (see Extended Methods in the Supporting Materials and Methods), Eq. 4 becomes the following equation:

where $t^{\ast }$ is the first passage time to the native state. According to Eq. 5, $A_{total}$ decreases exponentially with $t^{\ast }$. On the other hand, the relationship between $P_{nat}$ and $A_{total}$ is more complex. For realistic values of $k_d$ and with $T-t^{\ast }\approx T$, $k_{d}T\gg 1$. If $P_{nat}$ is close to 1 such that the entire argument of the exponential is small, $A_{total}$ is proportional to $P_{nat}$. Otherwise, $A_{total}$ is proportional to $P_{nat}/\left(1-P_{nat}\right)$.

Finally, we relate $A_{total}$ to cellular fitness, $f$. For convenience, we divide $A_{total}$ by total time $T$ to rescale it to the range $\left[0,1\right]$. We treat the protein as essential to cellular growth, and therefore its activity is related to cellular fitness. We use a metabolic flux-type equation to relate total protein activity to $f$ (35, 36, 37, 38):

where $A_0$ is a constant which sets the value of $A_{total}/T$ where fitness is half maximal. Together, Eqs. 5 and 6 formulate how protein folding kinetics and stability determine fitness in our model.

Evolutionary simulations using a lattice protein model

To explore how proteins evolve under prototypical functional selection, we ran evolutionary simulations that fix or reject mutations in a protein sequence based on the measured fitness. The evolutionary simulation scheme is illustrated in Fig. 1 C. We simulate evolution according to a discrete-generation monoclonal model in which, in each generation, a single arising mutation either fixes (takes over the entire population) or is lost. Our model organism has only a single gene corresponding to the protein under investigation. In every generation, a mutation is made to the current sequence, and the fitness of the trial sequence, $f\text{'}$, is evaluated using protein folding simulations. A selection coefficient is calculated as $s=\left(f\text{'}-f\right)/f$ (27), and the mutation is fixed with probability

where $N$ is the population size. Eq. 7 comes from classical population genetics (39).

In the fitness assessment, we use MC simulations of lattice model proteins undergoing translation. A lattice protein, as illustrated in Fig. 1 A, treats protein residues as a connected set of vertices on a cubic lattice. Our model uses 20 amino acid types whose interaction energy is given by a 20 × 20 interaction matrix (40).

The MC simulations of translation and folding have two phases: translation and posttranslation. During translation, MC dynamics alternates with elongation of the nascent chain at the C-terminus. The ribosome is not explicitly modeled; rather, the nascent chain C-terminus is simulated as if connected to a straight chain of infinite length (as illustrated in Fig. 1 A), representing unextruded residues in a ribosomal channel. There are no energetic interactions between the protein and the untranslated residues, but the channel does exclude a volume that is one lattice unit in width. The protein remains tethered for one additional elongation interval after the final residue is added, representing the ribosome release interval. This treatment of lattice protein translation matches what was used in a previous study (41). During the posttranslation phase, the protein is no longer tethered and has no conformational restrictions. Fig. S2 shows example folding trajectories for sequences studied in this work. Because MC simulations are stochastic, multiple trajectories are used to estimate $t^{\ast }$ and $P_{nat}$, and an average $A_{total}$ is used in Eq. 6. Additional details are described in the Extended Methods in the Supporting Materials and Methods.

The fitness function, as defined by Eqs. 5 and 6, by selecting for protein activity, effectively includes selection on folding kinetics and stability of proteins undergoing translation. As controls, we also ran the same evolutionary simulations under two alternative evolutionary scenarios in which we changed how we assessed fitness. In the first alternative scenario, we skip lattice protein translation. Instead, fitness assessments use MC simulations that begin with full-length proteins in a fully extended conformation, mimicking in vitro refolding and thereby making in vitro folding speed a determinant of fitness. The second alternative scenario ignores first passage time to the native state. In this case, MC simulations begin with full-length proteins already in their native conformations. $t^{\ast }$ in Eq. 5 is set to 0, and MC simulations only measure $P_{nat}$. Fitness in this scenario, therefore, depends on stability and does not depend on folding rate. We refer to sequences evolved without translation as “evolved, no translation,” and we refer to sequences starting in the native state as “evolved, no folding.” The degradation rate $k_d$ and other simulation parameters remained unchanged for these alternative evolutionary scenarios.

Simulation parameters and analysis

The simulation time unit, $t$, is defined as $t=\mathrm{MC}\mathrm{step}/\mathrm{protein}\mathrm{length}$, which accounts for the local nature of the MC move set. The key simulation parameters are the elongation interval and degradation rate $k_d$. Simulation parameters were chosen so that ratios of timescales between translation, protein folding, and degradation are biologically reasonable. An explanation of how simulation parameters were selected is given in the Extended Methods, and a listing of the parameters is shown in Table S1.

A key characteristic of different lattice proteins simulated in this work is the topology of the native structure that the proteins fold to. Different structures differ in the degree to which residues forming native contacts are separated in primary sequence. Contact order is defined as

where $L$ is the length of the protein, $N$ is the number of contacts, and $\Delta S_{i,j}$ is the separation in primary sequence between contacting residues $i$ and $j$ (28).

Nine lattice protein native structures were selected from the representative 10,000 structure subset of 27-mers used in previous works (42,43). Each native structure arranges the chain in a 3 × 3 × 3 cubic native fold. Three each of low, medium, and high contact order structures were chosen. Initial sequences for evolution were designed to be thermodynamically stable in the selected native conformations via Z-score optimization (44, 45, 46); all initial sequence Z-scores were below −50. Table S2 shows unevolved and evolved sequences for the protein structures used in this study.

Key quantities measured in simulations are first passage time to the native state, $t^{\ast }$, folding stability, $P_{nat}$, and native contacts formed. $P_{nat}$ is measured as the proportion of steps that the protein is in its native conformation, from the time that the protein reaches the native state until the end of the simulation. Proteins must be exactly in their native conformations to be considered native. For first passage time to the native state, $t=0$ is defined as the start of the posttranslational phase of simulation in which the protein chain has no conformational restrictions. Native contact counts are either normalized by the maximal possible number of native contacts at a particular chain length or by the number of native contacts in the full-length protein (28 for all lattice proteins in this work) and are typically reported as average values for each nascent chain length. The results presented in this work focus on for nascent chain lengths of 15–27.

Proteins evolve improved stability and kinetics

Beginning with sequences designed for native state thermodynamic stability, we initiated evolutionary simulations. The evolutionary trajectories of the nine proteins are shown across Figs. S3 and S4. In each trajectory, over the course of 1000 mutation attempts, a number of mutations were fixed. Occasionally, deleterious mutations were fixed because fitness evaluations using MC simulations are stochastic. Nonetheless, the fitnesses, folding stabilities, and folding speeds of the evolved sequences are improved over the unevolved sequences. The remaining discussion focuses on the evolved sequences from the end of evolutionary simulation and comparison with unevolved sequences or evolved sequences from alternative evolutionary scenarios.

The overall outcomes of evolutionary simulations are shown in Fig. 2, which compares the fitnesses, folding stabilities, folding times, and native energies of the unevolved, initial sequences to those of the evolved sequences. The total protein activity (Eq. 5), which determines fitness (Eq. 6), is determined by the combination of stability ($P_{nat}$) and folding kinetics (the distribution of $t^{\ast }$). Here, protein folding stability is shown in terms of the two-state free energy $\left(\Delta F/kT\right)\equiv -\ln \left(P_{nat}/1-P_{nat}\right)$ to illustrate differences in folding stability for $P_{nat}$ close to 1 more clearly. Because fitness calculations use the entire ensemble of first passage times, folding time distributions are illustrated using boxplots, with $t=0$ defined as the moment of release from the ribosome.

Initial, unevolved sequences have moderate stabilities that improve with evolution. The unevolved sequences also fold slowly relative to the degradation timescale $\left(1/k_d\right)$ of 200,000 time units. All medium and high contact order unevolved sequences have nonzero median first passage times, meaning that proteins fold posttranslationally in the majority of folding trajectories. The longest first passage times for unevolved sequences are greater than 10⁹ time units, indicating long-lived unfolded or misfolded states. In comparison, evolved sequences have first passage times that mostly fall within the degradation timescale. Evolution of folding times to be within the protein degradation timescale has been predicted by a previous study (47). The only evolved sequence with a nonzero median first passage time is that of structure 4, at 9000 time units. This indicates that cotranslational folding, technically defined as reaching the native state before release from the ribosome, occurs in the majority of folding trajectories for evolved sequences. The section that follows (Two Opposing Strategies for Reaching the Native State) discusses folding during translation in greater detail.

To probe the effect of different selection pressures, evolutionary simulations were also run under two alternative scenarios. The first alternative scenario, “no translation,” simulates in vitro protein folding, which starts full-length proteins in fully extended conformations. The second alternative scenario, “no folding,” starts simulations with proteins already in their native conformations (effectively setting $t^{\ast }$ in Eq. 5 to 0). The properties of sequences obtained from evolution under these two alternative fitness evaluation scenarios are shown in Fig. S5. Note that although evolved sequences were obtained under different evolutionary scenarios, for purpose of comparison, the fitnesses for all sequences shown in Fig. S5 were determined using the same, regular fitness evaluation involving MC simulations of translation and folding. Evolving for in vitro refolding kinetics (“no translation” scenario) results in proteins that fold in the context of translation just as fast as or even faster than the regular sequences which evolved with translation. The evolved “no translation” sequences, however, (except that of structure 5) are less stable than sequences evolved with translation. On the other hand, the results of evolution under the “no folding” scenario shows that long-lived kinetic traps that hamper folding are not eliminated under an evolutionary scenario where folding rates do not affect protein activity and fitness.

The bottom-most plot in Fig. 2 shows the native state energies of unevolved and evolved sequences. Evolved sequences for structures 1, 2, 3, and 6 have native energies that are about 25% lower than those of the evolved sequences for structures 4, 5, 7, 8, and 9. These native energies reflect how these proteins fold during translation—whether these proteins fold early on or late during translation—as we will discuss next.

Two opposing strategies for reaching the native state

Model proteins evolved sufficiently fast kinetics such that the majority of folding trajectories for evolved sequences achieve the native state before release from the ribosome (Fig. 2). To understand the nature of cotranslational folding in our model proteins, we examined folding trajectories and found that proteins have one of two different folding mechanisms. We first characterized trajectories by the fraction of possible native contacts formed at each nascent chain length, ${\langle Q\rangle }_L$. This measure normalizes the number of native contacts formed at a particular chain length by the maximal number of native contacts that can possibly be formed at that chain length. Q is averaged over all trajectory samples at each chain length and over all trajectories, providing an ensemble-average folding trajectory.

Examination of ${\langle Q\rangle }_L$ shows that the behavior of evolved proteins falls into two groups, as illustrated in Fig. 3. Group 1 consists of low contact order structures 1, 2, and 3 as well as medium contact order structure 6. For evolved sequences in Group 1, ${\langle Q\rangle }_L$ is close to 1 for chain lengths beyond 16 (Fig. 3 A). This shows that Group 1 proteins have cotranslational folding mechanism in which native-like conformations establish early on during translation. Group 2 consists of medium contact order structures 4 and 5 and high contact order structures 7, 8, and 9. Evolved sequences in Group 2 do not develop high ${\langle Q\rangle }_L$ values until nascent chains grow to about 25 residues in length (Fig. 3 B). Thus, Group 1 proteins fold early on during translation, and Group 2 proteins fold toward the end of translation.

There are also differences in how sequences evolved, as illustrated by differences in ${\langle Q\rangle }_L$ between evolved and unevolved sequences. For Group 1, ${\langle Q\rangle }_L$ either increased or remained unchanged as a result of evolution (Fig. 3 C). Minimal change in ${\langle Q\rangle }_L$ reflects cases in which ${\langle Q\rangle }_L$ is already close to 1 for unevolved sequences (structures 1 and 3). Interestingly, for Group 2, evolved sequences other than the sequence folding to structure 9 have lower ${\langle Q\rangle }_L$ values at intermediate nascent chain lengths 15–22 compared with those of unevolved sequences (Fig. 3 D) (Mann-Whitney U test between per-trajectory values at each length for each sequence, p < 0.0001). Evolved sequences form fewer native contacts at those lengths.

To further understand these two kinds of proteins, we examined individual folding trajectories for each structure. Here, we switch to quantifying folding by native contact counts, to illustrate development of native structure during translation. Native contacts are normalized by the total number of native contacts for the full-length native conformation, which is 28 for all lattice proteins in this work. We focus on structures 2 and 7 as representatives of Group 1 and Group 2, respectively, because evolved sequences for structures 2 and 7 have the highest folding stability out of all nine evolved sequences. Individual folding trajectories for unevolved sequences folding to structures 2 and 7 are shown in Fig. 4, A and B, respectively, and the corresponding trajectories for evolved sequences are shown in Fig. 4, C and D, respectively. For each trajectory, native contacts are averaged over all samples collected at each nascent chain length. Folding trajectories for unevolved and evolved sequences of all nine native structures are shown Fig. S6.

Sequences folding to structure 2 can stably occupy native-like conformations beginning at a nascent chain length of 16 residues (Fig. 4, A and C). There is an apparent bimodality to the folding trajectories, with the majority of trajectories occupying native-like conformations during translation, and a minority of trajectories in partially-native states. Most trajectories reach native-like states when the nascent chain is between 15 and 21 residues in length. Such trajectories show a steady buildup of native contacts as protein length increases. Trajectories that fail to fold by length 21, however, mostly fold posttranslationally. For instance, for the evolved sequence, 82% of folding trajectories reach the native state by length 21, but of the remaining trajectories, 67.5% remain unfolded at the end of translation. This demonstrates kinetic partitioning (48); the additionally translated residues produce kinetic traps in the folding energy landscape. The difference between the unevolved and evolved sequences is that for the evolved sequence, a greater proportion of trajectories reach the native state before the nascent chain reaches 21 residues in length, resulting in a higher proportion of cotranslational folding. Other proteins in Group 1 are similar: a majority of trajectories (>75%) achieve the native state before release from the ribosome, with folding largely proceeding through a steady increase of native contacts during translation (Fig. S6).

In contrast, for sequences folding to structure 7, folding does not occur until the nascent chain reaches a length of 25 residues (Fig. 4, B and D). Compared with the unevolved sequence, the evolved sequence also forms fewer native contacts at nascent chain lengths below 23 residues. Once the nascent chain passes 25 residues in length, however, evolved sequence trajectories show rapid folding to native-like conformations. Other proteins in Group 2 similarly fold only toward the end of translation (Fig. S6). Upon folding, the full native fold, minus just a few contacts, is achieved. Thus, although many folding trajectories for sequences in this second group technically exhibit cotranslational folding (folding before the end of translation), the folding process is closer to what would occur for full-length proteins.

Contact order is not a perfect predictor of cotranslational folding, as seen by the split of our medium contact order structures among Group 1 and Group 2. In particular, the evolved sequence for medium contact order structure 6 folds early on during translation. To explore this, we constructed two-dimensional contact maps in which average frequencies of residue-residue contacts at particular nascent chain lengths during translation are averaged across all folding trajectories. Figs. S7–S9 show these contact map-based plots of folding trajectories for unevolved and evolved sequences folding to structures 2, 6, and 7, respectively. For sequences folding to structures 2 and 6 (Figs. S7 and S8), evolution strengthened native contacts and weakened nonnative contacts. Both structures 2 and 6 equally support forming 17 native contacts at a nascent chain length of 20. By this point during translation, nascent chains adopt native-like conformations, from which the rest of the native structure forms. These native-like conformations are cotranslational folding intermediates, stable states on the path to the full-length native structure, similar to intermediates predicted in our previous work (20). Contact order roughly predicts whether such cotranslational folding intermediates can form. In contrast, for structure 7 (Fig. S9), only 10 native contacts can be formed when the nascent chain is 20 residues long, and the nascent chain instead forms more nonnative contacts. For sequences folding to structure 7, evolution weakened both native and nonnative contacts at shorter nascent chain lengths. These observations explain how the evolved sequence for structure 7 forms fewer native contacts at intermediate nascent chain lengths than does the unevolved sequence (Fig. 4 B versus Fig. 4 D).

We find further distinctions between Group 1 and Group 2 proteins when examining energetics as a function of nascent chain length. The free energies of native conformations are shown in Fig. S10. Mirroring observations made from analyzing native contacts and contact maps, for Group 1, evolution increased or maintained the stability of native conformations. On the other hand, for Group 2, evolution destabilized native conformations at shorter nascent chain lengths. This destabilization of native conformations is reflected in native state energies as well. Fig. S10 also shows the energies of native conformations for nascent chain lengths 15–27. For proteins in Group 2, the C-terminal residues for the evolved sequences contribute a greater fraction of the stabilization of the native state than is the case for the unevolved sequences. This pattern explains why evolved native energies for structures folding either early or late during translation are different in magnitude (Fig. 2, bottom). Overall, we find that the native structure of a protein determines how many native contacts are available at a nascent chain length, which decides whether the nascent chain can stably fold. This in turn influences whether evolution strengthens or weakens contacts made by residues at particular nascent chain lengths.

Kinetic characterizations contrast cotranslational folding with in vitro folding

For Group 1 proteins, folding trajectories that do not reach the native conformation when the nascent chain is 15–20 residues long mostly fold posttranslationally (Figs. 4, A and C and S6). This suggests that folding kinetics slows with increasing nascent chain length for these proteins. To investigate this, in vitro first passage times to the native state—the full-length protein starts in an extended conformation, and translation is not modeled—were measured. We compare in vitro first passage times of unevolved and evolved sequences with those of sequences obtained from evolution under the “no translation” scenario (Fig. 5); the latter sequences, by design, are optimized for in vitro folding. Compared with the unevolved sequences, both evolved and evolved, “no translation” sequences have faster folding kinetics. More significant, the four evolved sequences in Group 1 (structures 1, 2, 3, and 6) have slower first passage times than those of their “no translation” counterparts and those of evolved sequences in Group 2 (structures 4, 5, 7, 8, and 9). The in vitro first passage times of the evolved sequences that fold early on during translation are only moderately improved compared to those of their initial, unevolved counterparts.

These differences in kinetics between sequences reflect different selection pressures during evolution. When undergoing translation, evolved sequences in Group 1 fold to stable, native-like conformations at lengths of 15–20 residues. Additional translated residues then add to an existing native structure. Any slow-folding intermediates that form when folding from the fully unfolded state are thereby avoided. One consequence demonstrated here is that proteins that have evolved to fold cotranslationally have slow in vitro folding kinetics. Vectorial synthesis reduces the selection pressure for fast folding kinetics when proteins can start folding cotranslationally.

We further characterized the kinetics of our sequences by measuring first passage times to native conformations at chain lengths of 15–27 residues. We fit our data to a simple three-state, three-parameter model (see Extended Methods and Fig. S1). The fitted kinetic parameters are shown in Fig. S11. These results show how proteins in Group 1 have fast folding kinetics at intermediate chain lengths and slower kinetics as chain length increases.

Proteins that fold early on during translation benefit from midsequence slow codons

Thus far, our studies have examined nonsynonymous sequence changes, but another aspect of protein translation is that codon identity influences translation rates, which can affect protein folding efficiency (49). Furthermore, slowly translating rare codons have been associated with cotranslational folding (19,20). We next investigated the effects of changing the elongation intervals for different codon positions. Here, we restrict our investigation only to evolved lattice protein sequences.

Estimated translation rates for different codons in E. coli differ by up to an order of magnitude (50). We performed MC simulations of translation and folding in which the elongation interval for individual residues was increased 10-fold and measured the proportion of trajectories in which proteins folded cotranslationally. Note that increasing the elongation interval for the Nth codon means that the nascent chain spends additional time at a length of N-1 residues. Results from these simulations for evolved sequences folding to structures 2 and 7 are shown in Fig. 6, A and B, respectively, and results for all evolved sequences are available in Fig. S12. For nearly all sequences, increasing the elongation interval of codons at the C-terminus increases the proportion of cotranslational folding (Fig. S12), with proteins in Group 2 showing more substantial increases in cotranslational folding. Increased cotranslational folding due to slowly translated C-terminal codons is a somewhat trivial effect under our model, however, because folding at a nascent chain length of 25 or 26 residues is not very different from folding as a full-length 27-residue protein. Only sequences folding to structures 2 and 6, from Group 1, show increases in cotranslational folding from increased elongation intervals at midsequence positions. These positions reflect nascent chain lengths at which cotranslational folding intermediates become stable and resembles how rare codons are positioned before putative cotranslational folding intermediates in real proteins (20).

To confirm that increases in cotranslational folding proportion could provide a fitness advantage and therefore be selected by evolution, we performed evolutionary simulations in which mutations had the effect of slowing translation at a specific position in the sequence, mimicking the effect of synonymous mutation to a rare codon. In these evolutionary simulations, simulation trajectories were stopped once a single mutation was fixed (see Extended Methods). The distribution of fixed “synonymous mutations” for evolved sequences folding to structures 2 and 7 are shown in Fig. 6, C and D, respectively; the distributions deviate significantly from the neutral expectation of a uniform distribution (χ² test, p < 0.0001). As predicted by the cotranslational folding proportions in Fig. 6, A and B, synonymous substitutions to slower codons are selected for by our evolutionary simulations.

Our model results suggest that proteins that can fold early on during translation benefit from translational slowing at specific midsequence positions. Because our model proteins that fold early on during translation are also lower in contact order, we wondered whether a bioinformatic signature of cotranslational folding, conserved rare codons, would be more likely to be found in genes coding proteins with lower contact orders. A recent study from our group identified conserved rare codons in E. coli (20). The study moreover examined structurally characterized E. coli proteins and found that conserved rare codons are frequently positioned downstream of predicted cotranslational folding intermediates. We measured the contact orders of protein substructures preceding rare codons identified in this previous study. Here, a substructure is defined as the portion of the native structure from the N-terminus to 30 residues before the location of a codon of interest, a length which accounts for the ribosome exit tunnel (51).

For each gene with rare codons, we measured the contact order of the substructure corresponding to the first evolutionarily conserved rare codon (excluding N-terminal-rare codons). We then compared this distribution of contact orders to control distributions generated by measuring the contact orders of substructures preceding random positions in genes without conserved rare codons. This analysis was performed at multiple p-value thresholds for determining evolutionary conservation of rare codons (see Extended Methods). We found that protein substructures preceding rare codons have lower contact orders than those of protein substructures preceding randomly drawn positions (Fig. 6 E). Statistical significance declines with decreasing p-value threshold as the number of genes with qualifying rare codon regions decreases. The most statistically significant difference is found at a p-value threshold for rare codon conservation of 10⁻³ (Mann-Whitney U test (two-sided), p = 0.0019). The distributions of contact orders at this conservation threshold are shown in Fig. 6 F and have medians of 0.2038 and 0.2338 for substructures preceding rare codons and random substructures, respectively. Our simulation results suggest a mechanism for this observation: structures lower in contact order support evolution of a sequential, cotranslational folding pathway, facilitated by slowly translated rare codons.