Composability of regulatory sequences controlling transcription and translation in Escherichia coli

To explore the effects of regulatory element composition systematically, we designed and synthesized all combinations of 114 promoters with 111 RBSs (12,653 constructs in total; one combination resulted in an incompatible restriction site). We used 90 promoters from an existing library from BIOFAB: International Open Facility Advancing Biotechnology, 17 promoters from the Anderson promoter library on the BioBricks registry, 6 promoters from common cloning vectors, and a spacer sequence chosen as a negative control. From RBSs, we used 55 RBSs from the BIOFAB library, 31 from the Anderson BioBrick library, 13 from the Salis RBS Calculator (14) expected to give a range of expression, 12 commonly used RBSs from cloning vectors and the BioBrick Registry, and one sequence chosen as a negative control (reverse complement of canonical RBS sequence).

We synthesized the construct library using Agilent’s oligo library synthesis (OLS) technology (30) and cloned at ∼50-fold coverage into a custom medium-copy vector (pGERC), where the constructs drive expression of superfolder GFP (31) (Fig. S1). pGERC also contains an mCherry (32) reporter under constant expression by P_LTetO-1 (33) to act as a control for extrinsic noise (Fig. 1A). We grew the library to early exponential phase and characterized expression levels by flow cytometry. As expected, cells in the library expressed constant levels of mCherry, whereas expression levels of GFP varied over four orders of magnitude (Fig. 1B). We sequence-verified 282 colonies and found that 154 (55%) were error-free. We measured fluorescence levels of 144 of the unique error-free colonies individually to act as a defined set of controls (Fig. 1C).

We grew the entire pooled library to early exponential phase and performed multiplexed measurements of the steady-state DNA, RNA, and protein levels. We used sequencing, DNASeq and RNASeq, to obtain steady-state DNA and RNA levels, respectively, across the library (12). For obtaining protein levels, we used FlowSeq, which combines fluorescence-activated cell sorting and high-throughput DNA sequencing and is similar in design to recently published work (34, 35). Briefly, we sorted cells into 12 log-spaced bins of varying GFP/mCherry ratios; isolated, amplified, and barcoded DNA from each of the bins; and then used high-throughput sequencing to count the number of constructs that fell into each bin (Fig. 1 A and D). Using the read counts from each of the bins, we reconstructed the average expression level for each construct. Because our library contains a mixture of perfect and imperfect constructs, we only use reads that match the fully designed sequences perfectly, and thus filter out the effects of synthesis error.

Using DNASeq, we detected 98.5% of constructs and there was high concordance between technical replicates (R² = 0.997; Figs. S2 and S3). Most of the missing constructs and constructs with few DNA reads (which prevented accurate RNA level measurements) were expected to have very high expression levels, indicating either growth defects or cloning issues (Figs. S4 and S5). RNA level calculations also showed high concordance between technical replicates (R² = 0.995; Fig. S5). Overall, RNA levels varied by three orders of magnitude, but within a single promoter, the coefficient of variation was only 0.63 (Fig. 2, Left and Fig. S6). RNASeq data also allowed us to identify dominant transcriptional start sites for most promoters (Fig. S7). Eighty-seven percent of all promoters had one dominant start position (>60% of all mapped reads). Two promoters (marked with asterisks in Fig. S7) had very few uniquely mapping reads, did not show a strong start site, and showed unrealistic translation efficiency calculations. These observations indicated that we were missing most of the RNA (but not protein) reads from these promoters, possibly because of transcription starting after the end of the barcode sequence preventing unique identification. The 222 constructs (1.7%) containing these promoters were removed from all analyses.

Using FlowSeq, we were able to reconstruct expected protein levels for 94% of the constructs (Fig. 2, Right). As expected, individual constructs mostly fell into one to three contiguous flow-sorted bins (Fig. S8). The average protein expression levels displayed a large range and were highly correlated with the independently characterized constructs (R² = 0.94; Fig. 3A and Fig. S9). Due to the boundaries of our sorted bins, we determined that accurate quantitation was limited within a maximum and minimum range; 6.5% of the constructs were above and 14% were below this range (Fig. S9). Again, most constructs with missing measurements (insufficient or zero reads) contained combinations of strong promoters and RBSs. We calculated average promoter and RBS strengths by averaging transcription levels and translation efficiency (protein/RNA), respectively (Datasets S1 and S2). Promoters and RBSs were ordered and named based on their relative deviation from the average element (SI Materials and Methods).

Finally, we spiked 42 of the individual clones into a separate library (not analyzed here) and performed DNASeq, RNASeq, and FlowSeq to test reproducibility in biological replicates. Once again, protein levels were highly correlated with the individual measurements (R² = 0.91; Fig. 3C). Reconstructed values for RNA and protein levels also matched well between independent runs (R² = 0.89 and R² = 0.90, respectively; Fig. 3 B and D).

Our large dataset allows us to measure the extent to which combining regulatory elements led to predictable outcomes. Using a simple model for gene expression, where promoter strengths determine RNA levels and RBS strengths determine translation efficiencies, we reconstructed expected expression across all constructs and compared them with measurements (Fig. 4). We find that 80% of RNA levels and 64% of protein levels fall within twofold of the model predictions, and had R² = 0.92 and R² = 0.76 for RNA and protein, respectively (Fig. S10 A and B).

When unexpected levels of expression do occur, they can be quite large; the largest 5% of protein model deviations are off by an average of 13-fold. Such unpredictability makes precise engineering of large systems intractable. The ease and scale of these measurements indicate that rather than using prediction or standardization to construct a single design, we can construct a library to screen for desired expression levels when optimizing large genetic systems. Desired RNA and protein levels for an entire pathway of genes could be chosen from measurements across subsets of promoters and RBSs for each gene. For example, given a desired protein level, we can choose from many sequence-divergent promoter and RBS combinations that achieve desired transcription and translation strengths of GFP (Table 1).

We conducted a more detailed ANOVA (27), where both RNA and protein levels are independently determined by both promoter and RBS identity. This model is able to take into account effects such as the dependency of RNA levels on the translation rate. We found that the model resulted in a modestly better fit (RNA R² = 0.96, protein R² = 0.82; Fig. S10 C and D). Analysis of explained variance showed that 92% of the RNA levels can be explained by the promoter choice, whereas only 4% can be explained by the RBS choice and the remaining 4% are unexplained (Fig. 5A). For protein levels, both promoter choice (54% explained variation) and RBS choice (30%) are important, but a larger portion remains unexplained (16.7%).

To understand better how factors such as RBS choice can affect RNA levels, we examined interactions between RNA and protein levels. For example, several previous studies in Escherichia coli and Bacillus subtilis have shown that for particular model transcripts, increased ribosome binding or occupancy may enhance mRNA stability (36–42). Such studies have been hard to interpret due to the complex interactions between the ribosome, RNA degradation machinery, and transcript. We indeed find a significant and prevalent correlation between mRNA stability and RBS strength across all promoters. Given the size and sequence diversity of our library, it is likely that RBS strength is responsible for increased mRNA levels. Overall, we find that an ∼10-fold increase in translation efficiency correlates to an approximately threefold increase in RNA abundance (Fig. 5B). However, the effect is limited at the extremes; the difference between the weakest and strongest RBSs (an 87-fold increase in translation efficiency) corresponds to only an ∼4.3-fold increase in mRNA. As another example, many groups have found that secondary structure across the 5′ UTR and initial coding sequence can hinder effective translation (14, 43–46). In our data, we find that the correlation between secondary structure free energy across the UTR/GFP interface is significant (Fig. 5C). However, this metric of secondary structure is neither necessary nor sufficient, because many sequences with high secondary structure do not display reductions in expected expression, and vice versa. Improved models for how secondary structure interacts with ribosome binding could increase this correlation (14).

We used E. coli MG1655 (Yale Coli Genetic Stock Center no. 6300) for all experiments. The oligo library was constructed by Agilent Technologies using the OLS process (30). The design of pGERC is based on the synthetic plasmid pZS-123 (33), which allows independent expression from three promoters, and it was synthesized by DNA2.0, Inc. The amplified OLS pool was subcloned into 5α-electrocompetent E. coli (New England Biolabs) (giving an initial library size of ∼600,000 colonies), purified, and retransformed into MG1655, and several aliquots were frozen. Overnight cultures from both pooled experiments and individual clones were first diluted 1,000-fold grown at 30 °C in LB–Miller media shaking at 250 rpm (Infors HT Multitron) for 2–3 h until reaching an OD₆₀₀ of 0.15–0.25. Detailed information can be found in SI Materials and Methods.

From a single 300-mL culture of the library, pellets from four 50-mL aliquots of culture were frozen in liquid nitrogen, with the remaining culture saved for FlowSeq. Two technical replicates of DNA and RNA were isolated using Qiagen DNA and RNA Midiprep Kits. Ribosomal RNA was removed by means of a Ribo-Zero rRNA removal kit for metabacteria (Epicentre). The 5′ triphosphates were monophosphorylated by 5′ polyphosphatase (Epicentre) and then ligated to an RNA adaptor using T4 RNA Ligase (Epicentre). First-strand cDNA was made from a specific primer in superfolder GFP. Both DNA and cDNA were amplified and monitored by real-time PCR to prevent overamplification. Illumina adaptors and barcodes were then added, and sequencing was performed on a HiSeq 2000 in two separate PE100 lanes. A separate library that contained spike-ins from the 42 colonies underwent the same procedure. Detailed information can be found in SI Materials and Methods.

We used 50 mL of the library culture as prepared above for analysis by FlowSeq. We flow-sorted the cells into 12 log-spaced bins in three sequential runs sorting 4 bins each. Cells were then grown overnight to saturation and plasmid-prepped using a Qiagen Miniprep kit. A small aliquot was diluted, regrown, and subjected to flow cytometry to verify proper sorting. All data from library measurements are reported in GFP/RFP ratio units, which range from 1 to 255,000. The 12 minipreps were amplified again by real-time PCR, barcoded, and sequenced on a single-lane PE100 using a HiSeq 2000. Detailed information can be found in SI Materials and Methods.

Reads from all experiments were first aligned using SeqPrep (47) to form paired-end contigs for improved accuracy. Custom software was written to identify unique contigs and map them to library members using Bowtie (48) and grep (global search with the regular expression and printing all matching lines). DNASeq and RNASeq contigs were counted, where reads mapped uniquely and contained less than three mismatches. In addition, DNA contamination from RNASeq reads was identified and removed. Statistics, graphs, and tables were all generated using custom software written in Python, R, and the ggplot2 package (49). Detailed information can be found in SI Materials and Methods. All values used in intermediate and final calculations are enumerated in Dataset S3.