Operon mRNAs are organized into ORF-centric structures that predict translation efficiency

1Protein synthesis is the most energetically costly process in bacteria, consuming up to 50% of cellular energy. Thus, to optimize cellular efficiency, the rate of synthesis of each protein must be carefully controlled. In bacteria, operons are central to this process. Open reading frames (ORFs) with related functions are organized into operons that are transcribed as a single mRNA ensuring that operonic genes are transcriptionally co-regulated in response to various conditions. Additionally, the translation of each ORF in the operon is precisely tuned to cellular need. Indeed, the rate of protein production (i.e. the translation efficiency) of adjacent ORFs on a single mRNA can vary by as much as 100-fold, and members of protein complexes encoded on a single mRNA are generally translated in proportion to their stoichiometry. Understanding how the cell achieves optimal energy utilization critically depends on understanding how mRNA sequence features reliably drive ORF-specific translation.

2A number of mRNA features have been identified as contributing to the rate of translation of an ORF. Both the strength and accessibility of a Shine-Dalgarno (SD) sequence upstream from the initiation codon have been implicated in translatability. In support of the importance of SD accessibility, highly stable mRNA structures in direct proximity to the initiation codon diminish translatability and rare codons that disfavor RNA structure are enriched in positions immediately following the translation start site. Moreover, several studies examining either synthetic ORFs with a few bases difference (e.g. alterations to GFP), or fluorescent reporter assays studying the effect of multiple codon changes in the 5’ UTR and N-terminal coding sequences, find that models based on their predicted RNA structure at the translation start site are relatively successful at predicting their differences in translatability. Most recently, codon usage has emerged as an important variable for translation. A large study examining thousands of foreign ORFs concluded that except for the very initial nucleotides of the ORF, codon usage rather than mRNA folding propensity was the critical determinant for translatability.

3While these mRNA features are of clear value for predicting the translatability of exogenously expressed ORFs, several considerations suggest that they may not capture the key features that have evolved to set the translation efficiency of endogenous genes. First, all the high-throughput studies overexpressed the mRNAs they studied, which is known to perturb the charged tRNA pool and introduce biases in codon usage. Second, in the Boel et al. manuscript and in other studies, mRNAs were transcribed by T7 RNA polymerase, which not only elongates significantly faster than E. coli RNA polymerase but also removes the influence of the many endogenous E. coli RNA polymerase binding proteins that modulate its elongation rate. Where examined, such RNAs exhibit altered folding patterns. Thus, these transcripts likely have non-native structure. Third, these studies all used foreign mRNAs, which had not been subjected to evolution for precise tuning in E. coli. Finally, these studies primarily measured protein abundance, a quantity that is dependent on mRNA and protein stability as well as on the efficiency at which each ORF is translated.

4The goal of this work is to understand how E. coli establishes the relative expression of adjacent ORFs on the same mRNA. To accomplish this, we systematically assessed the translational efficiency (TE) of every ORF mRNA and then examined which of its features (e.g. secondary structure, codon usage, and the strength of ribosome binding site) correlated with its translatability. The translation efficiency of endogenous messages in E. coli could be probed with existing global technologies and the effects of codon usage with two metrics, tAI (tRNA adaptation index) and codon influence. However, in vivo mRNA structure has not previously been empirically evaluated at the global level in E. coli. We therefore adapted the dimethyl sulfate (DMS)-seq technique, which uses next-generation sequencing to determine chemical accessibility of RNA to DMS, to E. coli. Our studies point to a critical role of intrinsic ORF-wide differences in mRNA structure in allowing differential translation of ORFs sharing the same operonic mRNAs.

5New genomic technologies enable the determination of RNA structure on a global scale. DMS-seq uses next-generation sequencing to determine chemical accessibility of RNA to DMS (dimethyl sulfate), a reagent that reacts with unpaired adenosine and cytosine nucleotides. We adapted DMS-seq to E. coli to monitor global in vivo RNA structure (Figure 1A). By exploring the effect of coverage on reproducibility, we find that a read coverage of ~15 reads/nucleotide is sufficient for reproducible structure determination (Figure 1B), and used this cutoff in all subsequent analyses. Structures determined from E. coli-adapted DMS-seq are in excellent agreement with both the 16S rRNA crystal structure (Figure 1C), and a mutationally verified E. coli mRNA structure (Figure 1D).

6We quantified the degree of secondary structure for each ORF using the Gini index metric, which measures the variability in reactivity of A and C residues to DMS in the region being examined. A low Gini index indicates a relatively even distribution of DMS-seq reads and occurs in unstructured regions of the mRNA. A high Gini index occurs when a subset of residues is strongly protected from DMS reactivity and indicates a high degree of structure (Figure 1E). We found that the degree of RNA secondary structure varied greatly between ORFs: a small number are nearly as structured as rRNA, whereas some are close to the denatured state (Figure 1F).

7In contrast to the large variation in the degree of secondary structure among ORFs (Figure 1F), individual ORFs generally have fairly consistent Gini scores across their bodies – for example, the first and second halves of each ORF exhibit highly correlated Gini scores (Figure 2A). We tested whether the Gini scores across ORFs remained correlated in the absence of translation in two different conditions. First, we examined Gini scores of in vivo mRNA when translation initiation was inhibited with kasugamycin. We achieved rapid inhibition by using a ΔgcvB mutant, which has enhanced kasugamycin uptake rates (Figure 2—figure supplement 1A; ). Using ∆gcvB mutant is critical for this experiment because kasugamycin uptake by wild-type (WT) cells is slow enough to allow massive degradation of mRNA before ribosomes are cleared (see extended methods for protocol and Figure 2—figure supplement 1B–C for method validation, including demonstrating that ∆gcvB does not alter global mRNA structure). Second, we examined Gini scores of purified mRNA refolded in vitro at 37°C. In both cases, the translation-independent mRNA structures obtained from DMS-seq indicated that Gini scores across the ORF mRNAs remain correlated (Figure 2B–C). This correlation also holds true for computationally predicted mRNA structure of ORFs (Figure 2D). Moreover, the degree/extent of mRNA structures (henceforth referred to as structure) determined in these various ways are highly correlated with each other (Figure 2E–F). We conclude that mRNA is organized in ORF-wide structures that depend on the intrinsic sequence of the mRNA.

8We next examined whether structural correlation extends to adjacent ORFs on the same polycistronic (operonic) mRNA. We considered only those operons in which each ORF has an approximately equivalent mRNA levels, thus excluding those with significant internal promoters or terminators (see Materials and methods). Within operonic (polycistronic) mRNAs, the mRNA structure of adjacent ORFs can differ significantly (Figure 2G and Figure 2—figure supplement 1D), even when the start and stop codons of the adjacent ORFs overlap (Figure 2—figure supplement 1E). Thus, characteristic mRNA structures are a property of individual ORF mRNAs rather than of the entire polycistronic transcript.

9We next explored the relationship between the level of ORF-wide mRNA structure identified above with the TE of that ORF. We previously demonstrated that the overall rate of protein production can be accurately measured by an ORF’s average ribosome footprint density (number of footprints per unit length of the ORF), showing that protein copy number per cell determined from average ribosome footprint density was in superb agreement with that obtained by individually quantifying stable proteins in E. coli. Here, we build on that validated parameter, defining TE as the rate of protein production per mRNA, measured by normalizing average ribosome footprint density of an ORF with its mRNA abundance (i.e. RPKM of mRNA sequencing), with both measurements obtained from the same biological samples (see Materials and methods). Importantly, this metric is not affected by differences in either mRNA or protein abundance or stability.

10We found that the TE’s of E. coli endogenous ORFs in operonic mRNAs were highly negatively correlated with their level of ORF-wide mRNA structure (ρ = −0.75, Figure 3A: well-translated ORFs have less mRNA structure, while poorly translated ORFs have more structure). Consistent with the fact that the ORF-wide mRNA structures of adjacent ORFs in an operon can differ significantly (Figure 2G), the TE’s of adjacent ORFs can also differ significantly (Figure 3B, Figure 3—figure supplement 1A). Notably, ORF pairs with overlapping start and stop codons, believed to be translationally coupled, show essentially as much variability in their relative translation as non-overlapping ORF pairs (p=0.06, K-S test, Figure 3B), suggesting that the extent of coupling is variable. We then expanded this analysis beyond operons to all ORFs and found that the level of mRNA structure and TE are highly anti-correlated on all endogenous open reading frames (ρ = −0.76, Figure 3C). Importantly, the Gini scores of ORFs calculated from control RNA samples without DMS modification were not correlated to TE (ρ = 0.05, Figure 3—figure supplement 1B), indicating that Gini scores calculated from DMS-seq indeed reflect the level of mRNA structure and the potential sequencing bias/noise does not contribute to the correlation between TE and mRNA structure.

11Translation itself can influence mRNA structure as the helicase activity of translating ribosomes is likely to decrease the mRNA structure of highly translated ORFs more than that of poorly translated ORFs. We asked whether TE is correlated solely to the mRNA structure that results from ribosome unwinding or whether it is also correlated to the intrinsic mRNA structure that exists in the absence of translation. We find that when translation is inhibited in vivo (e.g. following kasugamycin treatment), the absolute correlation of TE to structure remains high but decreases somewhat (ρ = −0.58, Figure 3D), and that there is a small further decrease in correlation when mRNAs are refolded in vitro (ρ = −0.48, Figure 3E). Additionally, computationally predicted structures of entire ORFs also show robust correlation to their TE’s (ρ = −0.52, Figure 3F). The results are very similar when we confine ourselves to the 421 ORFs with ≥15 DMS reads/nucleotide in all conditions (Figure 4—figure supplement 2).

12We further dissected the influence of translation on ORF mRNA structure by determining how the difference in Gini score of in vivo mRNA with and without translation is related to its TE. We found that there is a tendency for mRNA to be more structured (higher Gini index) in the absence of translation (Figure 2E) and that mRNAs with the highest TEs had the greatest difference in their Gini’s (ρ = 0.52, Figure 3G). These data are consistent with the idea that unwinding by ribosomes contributes to the in vivo structure of highly translated genes. The decreased correlation of untranslated ORF mRNA structure to TE may result from removing the contribution of unwinding by translating ribosomes.

13Previous work on ORF translatability has pointed to the important role of sequences around the ORF start site. Using the ORFs that are common in all datasets and separated by ≥20 nt from the upstream ORF, we examined the correlation of TE with the level of mRNA structure only around the start site or extending further into the ORF (Figure 3—figure supplement 1C). We find increasing correlation with TE as successively larger regions of the ORF are considered in the structural analysis (−20nt to 40nt, 0 to 60nt, and 0 to 100nt relative to the gene start). Notably, the correlation of TE with extent of structure in either the first or second halves of the ORF are very similar, and the highest correlation is with the Gini of the ORF-wide mRNA structure.

14In toto, these analyses indicate that the linear sequences of bacterial mRNAs encode not only ORFs, but also ORF-wide secondary structures. These structures provide a rough blueprint for the TE of that ORF. Instructions from this blueprint are augmented by ribosomes and additional factors (see Discussion).

15We next examined the ability of the Shine-Dalgarno sequence and codon usage to predict TE. Data for all ORFs are presented in Figure 4, and that for the 421 ORFs in common between conditions are presented in Figure 4—figure supplement 2.

16Consistent with earlier studies, we found that the strength of the Shine-Dalgarno sequence does not have predictive power for TE, even after controlling for structure as measured by Gini index (Figure 4—figure supplement 1A).

17Codon usage, quantified by tAI (tRNA adaptation index) modestly correlates with TE (ρ = 0.34, Figure 4A). Interestingly, codon usage correlates more strongly with the overall rate of translation (i.e. average ribosome footprint density, ρ = 0.61) and ORF mRNA abundance (RPKM of mRNA sequencing, ρ = 0.48) than with their TE’s (Figure 4B–C). In contrast, the Gini score exhibits its highest correlation with TE (ρ = −0.76, Figure 3C) and is poorly correlated with mRNA abundance (ρ = −0.05) (Figure 4—figure supplement 1C). This suggests that codon bias may be evolutionally selected to correspond to the ORF expression level rather than to its translation efficiency. Additionally, there is evidence that codon usage correlates with mRNA half-life in both eukaryotes and prokaryotes. ORF-wide codon usage (tAI) and intrinsic mRNA structure appear to be largely independent variables, as they show little correlation with each other (Figure 4—figure supplement 1D). Although a novel metric quantifying codon influence was highly successful at predicting protein production from overexpressed exogenous genes transcribed by T7 RNA polymerase, it is relatively weakly correlated with TE under physiological conditions for endogenous genes (ρ = 0.29) (Figure 4D). This suggests that the codons providing efficient translation of an over-expressed transgene may differ from the efficient codons for an endogenous gene, as overexpression causes amino acid starvation and concomitant alteration of charged tRNA pools.

18Overall, ORF-wide mRNA secondary structure is by far the strongest and most significant predictor of endogenous TE compared to the other factors discussed above (Figure 4—figure supplement 2). A linear regression model that includes the addition of the Boel metric, tAI, and Shine-Dalgarno sequence strength showed marginal improvement in the predictive power compared to the ORF-wide structure alone (Figure 4—source data 1). Therefore, rather than being a driver for TE, codon optimization may be critical for highly expressed genes due to higher demand for these tRNAs and may play a role in setting the appropriate mRNA half-life.

19Our results thus far indicate that the rules for endogenous translation differ from those for overexpressed genes, particularly in the role of codon choice. Considering the fact that the expression of each tRNA species is tuned to the endogenous usage of its cognate codon(s), overexpressed transgenes are likely to perturb the balance between codon usage and tRNA abundance, creating a global translation defect. To directly test this hypothesis, we evaluated the effects of transgene overexpression containing one codon at a time. We constructed a synthetic gene with only one sense codon after the initiating codon and expressed this minimal ORF to directly assess the influence of a single tRNA and amino acid without additional complications from the protein product. When the minimal ORF contains the rare leucine codon CUA, which has only one cognate tRNA, we observed elevated ribosome occupancy at CUA codons in endogenous genes (Figure 4G). In particular, slow translation at CUA codons in the leuL leader sequence triggers the expression of leucine biosynthetic genes (Figure 4H), whereas overexpressing the minimal ORF with the common leucine codon CUG or without any coding sequence does not change the expression level or ribosome occupancy at leucine codons of endogenous genes (Figure 4G and I). These results suggest that overexpression of a rare codon and not a common codon can deplete the pool of free cognate aa-tRNA molecules, leading to global perturbation of translation. Cells expressing a transgene that contain more rare codons are thus under a different physiological state compared to WT cells solely expressing endogenous genes.

20Bacterial operons are densely packed with ORFs, as the majority of adjacent ORFs (62%) are separated by only 25nt or less (Figure 5A). Our finding that ORF mRNAs have a roughly similar degree of structure (Gini index) throughout their entire length (Figure 2A–D), but that the degree of structure of adjacent ORF mRNAs on polycistronic transcripts can differ significantly (Figure 2G) suggests that mRNA structure undergoes a sharp transition at ORF boundaries.

21We examined the structural organization of mRNA at ORF boundaries in polycistronic mRNAs. We find that the local degree of mRNA folding immediately downstream of the start site correlates with the TE of the downstream gene, but that this correlation rapidly diminishes upstream of the start site. Conversely, local mRNA structure upstream of the start site is only correlated with the TE of the upstream ORF (Figure 5B). This is true not only for mRNAs that are being translated (WT cells; Figure 5B) but also for untranslated mRNAs (kasugamycin-treated cells; Figure 5C), in vitro refolded mRNAs (Figure 5D) and computationally predicted mRNA structures (Figure 5—figure supplement 1A). Thus, mRNA structure undergoes a sharp transition at ORF boundaries, and polycistronic mRNAs consist of distinct ORF-length structural domains.

22The close packing of ORF mRNAs raises the issue of how they maintain distinct structural domains, and suggests that bacterial ORFs may be marked not only by start and stop codons, but also by features that assist within-ORF mRNA folding. To investigate this, we computationally predicted the structure of mRNA extending −250 to +250 nt from the translation start at the boundary of adjacent ORF pairs within the same operon. Because folding algorithms often predict a large ensemble of possible folds for a long stretch of RNA, we constrained the predictions by forcing positions that were highly DMS-modified to be unpaired in the predicted structures.

23Consistent with previous studies, we found a lack of structure in the immediate vicinity of the start sites for most ORFs (Figure 5E). Downstream from this structure-free zone (25–50 nt), endogenous mRNA has a high propensity to base pair with regions further downstream, that is pairing within the same ORF (Figure 5F). Conversely, nucleotides located 25-50nt upstream of the start site have a strong preference for interacting with regions further upstream in the preceding ORF (Figure 5F). Importantly, in vivo mRNA without translating ribosomes and in vitro probed mRNA (Figure 5E–F) also showed such preferences. Thus, the sharp transition in the directionality of base-pairing around start sites is driven by the mRNA sequence itself, promoting ORF-centric units of secondary structure.

24We experimentally investigated the effects of disrupting a region that promotes independent mRNA folding within adjacent ORFs. The dusB-fis operon is composed of a highly structured upstream gene (dusB) and a poorly structured downstream gene (fis) separated by 25 nucleotides. The two ORFs have an ~100 fold difference in TE (Figure 6A). Previous work indicated that the upstream dusB gene has a stem-loop structure near the 3’ end of the gene; that mutationally disrupting the stem-loop (Mutation M3; Figure 6A) decreased translation of fis; and that restoring base pairing by a second mutation (M2) restored fis translation for unknown reasons. After confirming these results (Figure 6B), we performed DMS-seq on WT and mutant cells to determine whether destroying the stem-loop decreased fis translation by reducing the structural isolation of dusB and fis.

25A model of the structure of the dusB-fis interface constrained by DMS-seq data (Figure 6D) indicates that the dusB and fis ORFs are structurally distinct in WT and double mutant (M3/M2 or M3:2 cells) (Figure 6E), but that M3 increases the structure of fis mRNA (Figure 6C and Figure 6—figure supplement 1). In the M3 mutant, the -58 ~ −53nt region (blue) pairs with the +9 ~ +14 nt region of fis (red), rather than forming a stem-loop structure within dusB as it does in WT and M3:2 cells (Figure 6E–F). The increased structure of fis mRNA in the M3 mutant starts at the ORF boundary and spreads across the entire downstream coding region of fis (Figure 6—figure supplement 1). Thus, mutation M3 induces long-range interactions between mRNA of the dusB and fis ORFs, which are normally structurally insulated from each other.

26In toto, our results suggest that specific sequences isolate mRNA folding within adjacent ORFs thereby minimizing structural crosstalk between adjacent ORFs. Disruption of structural boundaries affects both local and long-range mRNA folding, which is likely to be critical for programming the degree of translational isolation between ORFs on the same mRNA.

27Translation is a highly controlled process in bacteria, making it critical to understand the mRNA features contributing to differential translatability. Numerous studies have investigated the important question of which features control protein production from overexpressed, foreign ORF mRNAs, identifying codon usage and local structure around the translation start site as key variables. However, these studies have left open the question of which mRNA features regulate endogenous translation. The importance of this question is highlighted by the observation that the rate of protein production from each ORF in a polycistronic mRNA can vary as much as 100-fold. Our global study now examines this issue. Our principal finding is that ORF mRNAs have modular structures within polycistronic mRNAs and that ORF-wide mRNA structure rather than codon usage correlates most strongly with the translation efficiency of endogenous ORFs.

28Our analysis of mRNA structure revealed the unanticipated finding that operonic mRNAs have modular structures. Each ORF mRNA in the operon has a characteristic degree of structure, with highly correlated Gini scores between their first and second halves. This correlation persists in the absence of translation, when mRNAs are refolded in vitro and when structure is determined computationally (Figure 2A–D). In stark contrast, there is little correlation between the extent of structure in adjacent ORFs (Figure 2G). Additionally, and consistent with earlier computational findings, we observe a small ~25 nucleotide region beginning at the translation start that is more unstructured than the remainder of the ORF (Figure 5E). Thus, polycistronic mRNAs consist of a series of ORF-wide modules each with characteristic but different extents of structure, punctuated by regions of low basepairing at the translation start site (Figure 7). Maintenance of a common degree of structure throughout an ORF suggests that this parameter, like reduced structure at the start of ORFs, is a selective force in the evolution of ORF sequence, providing yet another constraint on mRNA sequence beyond codon adaptation.

29We find that the TE of each ORF correlates very highly, and most strongly with the ORF-wide extent of mRNA structure. We have begun to deconvolute the ‘chicken and egg’ problem of whether mRNA structure is a cause or a consequence of translation by examining the correlation of TE to ORF-wide structure when translation is inhibited. This removes the ribosome contribution but retains vectorial folding, RNA binding proteins and in vivo concentrations of salts and macromolecules. Untranslated mRNA structure is highly correlated with TE but less so than translated mRNA (Figure 3C–D). Moreover, the difference in the mRNA structure of an ORF with and without translation is highly correlated to its TE (Figure 3G). Thus, poorly translated mRNAs, have virtually identical extents of structure with and without translation, but more highly translated RNAs become increasingly more unstructured. Finally, computationally predicted structures or those obtained from in vitro refolded mRNAs correlate somewhat more poorly with TE (ρ = −0.52 or −0.48 respectively) than the structure of untranslated in vivo mRNA (ρ = −0.58) (Figure 3D–F). This suggests that features of the in vivo cell, besides translation by the ribosome, may also affect mRNA structure.

30Taken together, these results suggest that the intrinsic mRNA sequence itself encodes a rough blueprint for the ORF-centric mRNA structures that are predictive of TE. These structures are then amplified by translation and other features of the living cell, in a self-reinforcing loop, to provide the structure that ultimately specifies the translation of each ORF.

31Interestingly, E. coli has chosen to insert a predominance of its low TE ORFs into operons where adjacent genes have moderately high TE’s (for example, the rpsF-priB-rpsR-rplI operon shown in Figure 3—figure supplement 1A). All of its 10 lowest TE ORFs, and 86% of its ORFs in the bottom 10% of TEs are located in operons, compared with 58% of all ORFs. The evolutionary advantage of this arrangement is not known, but may relate either to decreasing transcriptional noise or to mRNA stabilization.

32The necessity for achieving widely different TEs for adjacent ORFs in operons may have driven the evolution of the ORF-centric mRNA folding strategy. As the translation termination codon of most ORFs is separated by less than 25nt of untranslated mRNA from the start site of the downstream ORF, the abundant ribosomes of the highly translated ORF could transiently open the structure of the poorly translated ORF and increase the accessibility of its start site. The propensity for in-ORF mRNA folding at both the beginning and ends of ORFs may prevent the upstream ORF from influencing the structure and hence TE of the downstream ORF, effectively insulating each ORF from its neighbors. We have identified small regions, located about 25–50 nucleotides both downstream and upstream of ORF translation start sites that preferentially pair within their ORFs. These regions may reinforce the folding barriers between adjacent ORF mRNAs, as we demonstrated for the dusB-fis operon. Interestingly, RNA polymerase pausing is enriched at translation start sites and this may reinforce ORF-centric structural insulation by allowing ORFs to fold independently during the pioneer round of translation.

33It is likely that the extent to which adjacent ORFs are insulated has been tuned. Approximately 15% of ORF pairs have overlapping stop and start codons and translational coupling has been demonstrated in some cases. This overlap may enable upstream ribosomes to influence downstream ORF translation by unwinding mRNA structure, thereby promoting translational coupling. Indeed, it is likely that the propensity for in-ORF basepairing is slightly weaker for overlapping ORF pairs than for non-overlapping ORF pairs (Figure 5—figure supplement 1B–C).

34The precise role of modular ORF structures that provide a rough blueprint for TE has not yet been established. It is certainly possible that the mRNA structure of an entire or a significant fraction of the ORF is required to define translation initiation, as has been demonstrated experimentally for rimM. Alternatively, a constant degree of ORF-wide mRNA structure may be the most robust way to ensure the appropriate amount of mRNA structure around the translation initiation site. In support of this idea, a recent study using in situ codon mutation of the E. coli essential gene infA showed that mutations of codons even far downstream from the start of the gene can be deleterious if they disrupt the native 5’ RNA conformation via long-range structural interactions predicted computationally. ORF-wide structures may also play additional roles. For high TE (poorly structured) ORFs, extended lack of structure may provide the landing pad necessary to capture a large pool of non-specifically bound 30S subunits to wait for opening of the SD and start codon, the so-called ‘standby model’ of translation initiation. Additionally, the ORF-centric mRNA folding strategy may have been driven by the necessity for adjacent ORFs to have discrete, often significantly different TEs. Finally, ORF-wide mRNA structures may help set the rate of endonucleolytic cleavage. The function of these modules is an important area for future inquiry.

35Although the TE of endogenous ORFs is primarily predicted by the extent of its mRNA structure, translatability of overexpressed foreign ORFs appears to be strongly driven by codon usage and tRNA limitation. This difference may arise from the fact that codon usage and tRNA abundance are largely balanced under physiological conditions, but become imbalanced when foreign ORFs are overexpressed, and we have directly demonstrated that this is the case (Figure 4G–I). This suggests that synthetic biologists and the cell tune translation in different ways. However, synthetic biologists struggle to robustly program differential translation of ORFs on the same mRNA. Our finding that polycistronic mRNAs consist of ORF-wide modules with set amounts of structure that are insulated from their neighbors may be key to this issue. Design approaches that incorporate appropriate mRNA structures may have the potential to produce the finely tuned synthesis rates observed in natural operons.

36E. coli K-12 MG1655 (RRID:SCR_002433) was used as the WT strain. All culture experiments were performed in MOPS medium supplemented with 0.2% glucose, all amino acids except methionine, vitamins, bases and micronutrients (Teknova, Hollister CA). Cells were grown in an overnight liquid culture at 37°C, diluted to an OD420 = 0.001 in fresh medium and grown until OD420 reached 0.4 where samples were collected. Multiple deletion strains were generated by transduction of FRT-flanked deletion alleles from the Keio collection followed by marker excision by Flp recombinase. All major experiments were biologically repeated for at least twice (see raw data files for sequencing data).

37In the experiment testing the effects of overexpressing the rare CUA leucine codon and the common CUG leucine codon, plasmids with pBR322 origin of replication was constructed to have a mini ORF ATGCTATAA or ATGCTGTAA driven by an IPTG-inducible promoter. The plasmid also contains lacIq to increase the expression of lac repressor. MG1655 containing the control plasmid (without mini ORFs) and MG1655 containing the plasmid with CUA or CUG mini ORF were grown overnight in MOPS rich glucose medium with 100 µg/ml ampicillin, diluted 1:1000 into 250 ml pre-warmed fresh medium containing 1 mM IPTG next morning. Cells were grown at 200 rpm at 37°C and harvested when OD600 reached 0.3 by vacuum filtration.

38The protocol for bacterial ribosome profiling with flash freezing was described. Briefly, 200 mL of cell culture were filtered rapidly and the resulting cell pellet was flash-frozen in liquid nitrogen and combined with 650 µL of frozen lysis buffer (10 mM MgCl2, 100 mM NH4Cl, 20 mM Tris-HCl pH 8.0, 0.1% Nonidet P40, 0.4% Triton X-100, 100 U/mL DNase I (Roche, St. Louis MO), 1 mM chloramphenicol). Cells were pulverized in 10 mL canisters pre-chilled in liquid nitrogen. Lysate containing 0.5 mg of RNA was digested for 1 hr with 750 U of micrococcal nuclease (Roche) at 25°C. The ribosome-protected RNA fragments were isolated using a sucrose gradient followed by hot acid phenol extraction. Library generation was performed using the previously described strategy detailed below.

39For experiments performed in parallel with ribosome profiling, total RNA was phenol extracted from the same lysate that was used for ribosome footprinting. For experiments performed independently of ribosome profilng experiments, and for total mRNA used for in vitro DMS-seq experiments, 4 mL of OD420 = 0.4 culture was added to 500 µL of ice-cold stop solution (475 µL of 100% EtOH and 25 µL acid phenol), vortexed, spun for 2 min at 8000 rpm, and the cell pellet was flash frozen in liquid nitrogen. Total RNA was then hot acid phenol extracted. For mRNA-seq experiments, ribosomal RNA and small RNA were removed from the total RNA with MICROBExpress (Ambion, Grand Island NY) or Ribozero (Epicenter, Madison WI) and MEGAclear (Ambion), respectively. mRNA was randomly fragmented as described. The fragmented mRNA sample was converted to a complementary DNA library with the same strategy as for ribosome footprints.

40The footprints and mRNA fragments were ligated to miRNA cloning linker-1 (IDT) 5rApp/CTGTAGGCACCATCAAT/3ddC/, using a recombinantly expressed truncated T4 RNA ligase 2 K227Q produced in our laboratory. The ligated RNA fragments were reverse transcribed using the primer 5'/5Phos/GATCGTCGGACTGTAGAACTCTGAACCTGTCGGTGGTCGCC GTATCATT/iSp18/CACTCA/iSp18/CAAGCAGAAGACGGCATACGAATTGATG GTGCCTACAG 3'. The resulting cDNA was circularized with CircLigase (Epicentre), and PCR amplification was done as described previously.

41For in vivo DMS modification, 15 mL of exponentially growing E. coli were incubated with 750 µL DMS. Incubation was performed for 2 min at 37°C. For kasugamycin (ksg) experiments, ksg was added to a final concentration of 10 mg/mL to ΔgcvB cells for 2 min at 37°C prior to DMS modification. Untreated ΔgcvB cells were also modified by DMS and collected in parallel as control. DMS was quenched by adding 30 mL 0°C stop solution (30% β-mercaptoethanol, 25% isoamyl alcohol), after which cells were quickly put on ice, collected by centrifugation at 8000 g and 4°C for 2 min, and washed with 8 mL 30% BME solution. Cells were then resuspended in 450 µL total RNA lysis buffer (10 mM EDTA, 50 mM sodium acetate pH 5.5), and total RNA was purified with hot acid phenol (Ambion). For in vitro DMS modifications, mRNA was collected from cells that were not treated with DMS. Two micrograms of mRNA was denatured at 95°C for 2 min, cooled on ice and refolded in 90 µL RNA folding buffer (10 mM Tris pH 8.0, 100 mM NaCl, 6 mM MgCl2) at 37°C for 30 min then incubated in either. 2% DMS for 1 min (95°C) or 4% DMS for 5 min (37°C). The DMS reaction was quenched using 30% BME, 0.3 M sodium acetate pH 5.5 and precipitated with 2 µL GlycoBlue and 1X volume of isopropanol.

42Sequencing libraries were prepared as described. Specifically, DMS treated mRNA samples were denatured for 2 min at 95°C and fragmented at 95°C for 2 min in 1x RNA fragmentation buffer (Zn2+ based, Ambion). The reaction was stopped by adding 1/10 vol of 10X Stop solution (Ambion) and quickly placed on ice. The fragmented RNA was run on a 10% TBU (Tris borate urea) gel for 60 min. Fragments of 60–70 nucleotides in size were visualized by blue light (Invitrogen, Carlsbad CA) and excised. Reverse transcription was performed in a 20 µL volume at 52°C using Superscript III (Invitrogen), and truncated reverse transcription products of 25–45 nucleotides were extracted by gel purification.

43Sequencing was performed on an Illumina HiSeq 2000 or 4000. Sequence alignment with Bowtie v. 0.12.0 mapped the footprint data to the reference genomes NC_000913.fna obtained from the NCBI Reference Sequence Bank. Sequencing data from mutated strains were aligned to appropriately modified genome. For ribosome footprint and mRNA-seq samples, the center residues that were at least 12 nucleotides from either end were given a score of 1/N in which N equals the number of positions leftover after discarding the 5' and 3' ends. For DMS-seq samples, read counts were assigned to the base immediately 5' of the 5' end of each read, which is the base that was DMS-modified.

44Data analysis was performed with custom scripts written for R 2.15.2 and Python 2.6.6. To calculate mRNA abundance, the number of mRNA sequencing reads mapped to a gene, following a Winsorization applied to trim the top and bottom 5% of reads, was divided by the length of the gene to yield the number of reads corresponding to the message per thousand bases of message per million sequencing reads (RPKM). The protein synthesis rate of individual ORFs was measured by average ribosome footprint density of the ORF calculated as described in. First, genes with less than 128 reads mapped and genes with unconventional translation events were excluded from the analysis, which include (1) genes encoding selenoproteins (e.g. fdhF, fdoG, fdnG); (2) proteins with nearly identical coding sequences (e.g. gadA and gadB, ynaE and ydfK, ldrA and ldrC, ybfD and yhhI, tfaR and tfaQ, rzoD and rzoR, pinR and pinQ). Second, sequencing reads from ribosome profiling mapped to the first and last five codons of the gene were excluded to remove effects of translation and termination. Third, correction for the variations in translation elongation rate was done in three steps as described in ): (1) correcting for the elevated ribosome footprint density observed for the first 50–100 codons; (2) correcting for the elevated density at the ribosomal anti-Shine-Dalgarno (aSD) site; (3) correcting for other possible ribosome pausing using 90% Winsorization, by removing the top and bottom 5% of the ribosome profiling signal for each gene. Finally, the average ribosome footprint density of a gene was calculated by dividing the corrected number of mapped ribosome footprint reads by the corrected length of the gene. Translation efficiency of a gene was calculated by normalizing the average ribosome footprint density by the mRNA abundance of the gene (defined above). The average ribosome footprint density (i.e. protein synthesis rate), mRNA abundance, and translation efficiency of genes from different samples are listed in Supplementary file 1–4.

45For identification of unpaired bases, raw DMS-seq data was normalized to the most highly reactive residue after removing outliers by 95% Winsorisation (all data above the 95th percentile is set to the 95th percentile). Bases with DMS-seq signal greater that 20% of the signal on the most highly reactive residue (after normalization) were called ‘unpaired’. For determination of rimM mRNA structures, a Viennafold (http://rna.tbi.univie.ac.at/) minimum free-energy model of the indicated region was generated, constrained by bases experimentally determined to be unpaired in the indicated dataset. Color-coding by DMS signal was done using VARNA (http://varna.lri.fr/).

46The secondary structure models for E. coli ribosomal RNAs were downloaded from Comparative RNA Website and Project database (http://www.rna.icmb.utexas.edu/DAT/3C/Structure/index.php). The crystal structure model was downloaded from Protein Data Bank (http://www.pdb.org, PDB entries 3I1M, 3I1N, 3I1O, and 3I1P). The solvent-accessible surface area was calculated in PyMOL, and DMS was modeled as a sphere with 2.5 Å radius (representing a conservative estimate for accessibility because DMS is a flat molecule). Accessible residues were defined as residues with solvent accessibility area of greater than 2 Å2. Unpaired residues in DMS-seq data were identified as described above. True positive bases were defined as bases that are both unpaired in the secondary structure model and solvent-accessible in the crystal structure model. True negative bases were defined as bases than are paired (A-U or C-G specifically) in the secondary structure model. Accuracy was calculated as the number of true positive bases plus the number of true negative bases divided by all tested bases.

47The R package ‘ineq’ (https://github.com/cran/ineq) was used to calculate Gini indices over As and Cs in the region specified for each experiment. For each DMS-seq sample, Gini indices were calculated only for genes that had greater than an average of 15 reads per nucleotide across the ORF. Genes with discontinuous mRNA-seq reads (due to an early termination event or an internal promoter, 1% of genes) were excluded from the analysis. Specifically, Gini indices were calculated on mRNA-seq data, and a cutoff was created based on two standard deviations from the mean. The Gini indices of genes from different samples were listed in Supplementary file 1–4.

48Adjacent ORFs in annotated operons often have differing levels of mRNA-seq reads, suggesting that they are not always on the same mRNA molecule. To identify adjacent ORFs expressed as a single operon, we assessed mRNA-seq data for equivalent mean message level and for signal continuity, as described below. Equivalent mean message level was assessed by first determining the variability in mean mRNA-seq read density within individual ORFs. There is a single transcript that extends over the entire body of the large majority of ORFs, and so the variability in mean read density level in the first half of each ORF was compared to mean read density in the second half of each ORF, and the variability in this distribution was used to define a cut-off for ORFs on a single message. Adjacent ORFs that fell within a 2σ cutoff in mean level (calculated to be a 1.5-fold difference in mRNA level) were determined to have equivalent mRNA level and were then assessed for signal continuity. Signal continuity was assessed by first determining the distribution of read density in windows within messages. Gini index of mRNA-seq signal was calculated for all 80nt windows within ORF bodies, and the variability in this distribution was again used to define a cutoff for continuous mRNA regions. Gini index were then calculated for 80nt windows tiling the region between adjacent ORFs. Gene pairs that fell within a 2σ cutoff defined by the intra-ORF distribution were considered to be a pair of adjacent ORFs on a single message.

49To determine the directionality of mRNA base pairing at ORF boundaries, sequence from −250 to +250 nt relative to the translation start site of the downstream gene was extracted for each adjacent pair of ORFs. A Viennafold (http://rna.tbi.univie.ac.at/) minimum free-energy model of each 500nt sequence was generated (constrained by DMS-seq data). The predicted probability of each base interacting with every other base in each mRNA structure model was then extracted from the Viennafold output. For each position, the probability of that position base pairing with any position within the upstream or downstream 60nt was then calculated. The ratio between summed upstream over downstream interaction probability across all mRNAs was then calculated for each position.

501 µC of Perkin Elmer EasyTag 35S labeled methionine (Product # NEG709A) was mixed with 5 µL 288 µmol unlabeled methionine and 24 µL media. At the time of capture, 900 µL of culture was added to methionine mix, and was labeled on a shaker for 1 min at 37°C. After labeling, 100 µL of ice-cold 50% trichloracetic acid (TCA) was added to the sample, which was vortexed and placed on ice for at least 20 min to allow precipitation. Samples were then counted by running 100 µL of sample through a 25 mm APFC glass fiber filter (Millipore APFC02500, Hainesport NJ) pre-wetted with 750 µL of 5% TCA on a vacuum stand, and washing three times with 750 µL 5% TCA and three times with 750 µL 100% ethanol. Filters were then placed in MP Ecolume scintillation fluid and counted.

51We used the RBS Calculator established by Salis et al downloaded from http://www.github.com/hsalis/Ribosome-Binding-Site-Calculator-v1.0 to predict the strength of Shine-Dalgarno sequence.

52The measurement of tAI (tRNA adaptation index) was adapted from the previous works, which gauges the availability of tRNAs for each codon within a gene. tAI incorporates different efficiency weights of the wobble interactions between codons and anticodons, with wi is defined as the relative adaptiveness value of codon i of a gene. The final tAI of a gene is the geometric mean of all its codons as shown below.

53ik is the kth codon of the gene and L is the length of the gene (excluding start and stop codons).

54Wild-type, dusB-M3, and dusB-M3:2 cells were grown in MOPS rich medium at 37°C till log phase (OD420 ~0.3). 1 mL cells were collected, resuspended in 30 µL SDS loading buffer, and boiled for 5 min. 10 µL of cell lysate was subject to Blot 12% Bis-Tris plus gel (ThermoFisher scientific, Grand Island NY). Proteins were transferred to a nitrocellulose membrane (BIO-RAD, Hercules CA). The membrane was first incubated with rabbit polyclonal anti-Fis antibody (a kind gift from Dr. Reid C. Johnson at UCLA) and mouse monoclonal anti-RNAP β subunit antibody (abcam, Cambridge MA), and then incubated with goat anti-rabbit IgG IRDye 800CW and anti-mouse IgG IRDye 680RD secondary antibodies (LI-COR, Lincoln NE). The blots were visualized and quantified by an Odyssey imaging system (LI-COR). The amount of Fis protein in each sample was normalized against the amount of RNAP β subunit in the same sample.

55All the processed and raw datasets of sequencing experiments were uploaded to NCBI GEO database with accession number GSE77617.