Global analysis of protein synthesis in Flavobacterium johnsoniae reveals the use of Kozak-like sequences in diverse bacteria

Cell culture and library preparation

Flavobacterium johnsoniae strain UW101 was grown at 30°C to mid-logarithmic phase in rich CYE medium (28). To halt translation, cultures were rapidly chilled by mixing with crushed ice. Cells were collected by centrifugation and lysed via a freeze-thaw method (29). Ribosome-protected and total RNA fragments were isolated from the lysate, essentially as described previously (30). Briefly, to obtain the ribosome footprints, lysates were treated with 3 U/μl micrococcal nuclease (Roche) in the presence of 5 mM CaCl₂ for 1.5 h at 25°C, and 70S monosomes were purified via sucrose gradient sedimentation. Proteins were extracted with phenol/CHCl₃, RNA was precipitated with ethanol, and then RNA fragments in the 20–35 nt size range were gel purified. In parallel, cell lysates were directly extracted with TRIzol™ (Life Technologies), total RNA was subjected to limited alkaline hydrolysis, and fragments in the 30–50 nt range were gel purified. Corresponding cDNA libraries (ribo-seq and RNA-seq) were prepared as described (31).

High-throughput sequencing, read processing and reproducibility

Sequencing was performed using a HiSeq2500 in the Ohio State Comprehensive Cancer Center Genomics Shared Resource in 50 basepair (bp) single-end mode. The three RNA-seq replicates initially contained 31 621 290, 30 635 429 and 29 764 896 single-end reads, and the three ribo-seq replicates initially contained 28 364 922, 17 598 576 and 27 203 601 single-end reads. 3′ adapter sequences were trimmed from the single-end reads using Skewer (version 0.2.2) (32) and its default settings, with the exception of the minimum allowed read length (‘−’) being set to 15. Sequences were then aligned to the F. johnsoniae genome (NCBI Reference Sequence: NC 009441.1) using the Bowtie 2 aligner (version 2.2.9) (33), using its default settings. Alignments were then converted to the BAM format using SAMtools (version 0.1.18) (34). Reads corresponding to ribosomal RNA (rRNA) and transfer RNA (tRNA) were removed from alignment files using the BEDtools (version 2.26.0) utility ‘intersect’ (35).

For the RNA-seq replicates, alignment efficiencies were 88.4% (65.2% multimapped), 88.1% (62.2% multimapped), and 88.5% (63.3% multimapped), of which 79.6%, 78.0% and 77.9% aligned to rRNA/tRNA genes. For ribo-seq replicates, alignment efficiencies were 81.6% (77.2% multimapped), 72.4% (68.7% multimapped) and 62.6% (60.1% multimapped), with 77.0%, 68.5% and 59.9% aligning to rRNA/tRNA genes. Analysis of gene expression (transcription and translation) entailed counting aligned RNA-seq fragments and ribo-seq footprints that overlapped coding sequence annotations using HTSeq (version 0.6.0) (36). Replicate-to-replicate comparisons of read counts per gene showed that the RNA-seq datasets were reproducible, with Spearman's correlation coefficients (r) of 0.982, 0.981 and 0.980 (Supplementary Figure S1). Analogous replicate-to-replicate comparisons showed that the ribo-seq datasets were also reproducible, with Spearman's correlation coefficients (r) of 0.958, 0.963 and 0.963 (Supplementary Figure S1).

Calculation of average ribosome density (ARD)

For this study, we chose the highly expressed genes to act as a representative subset. To do this, we combined the data of the replicates and summed the per gene counts for RNA-seq fragments and for ribo-seq footprints. Then we rank-ordered (highest to lowest) the entire set of 5138 protein coding genes of F. johnsoniae by the number of RNA-seq fragments per gene length. From this list we arbitrarily selected the top third (1712 genes) as our representative set. Gene counts for both RNA-seq and ribo-seq were then normalized by the total number of counts per million from this representative set. Then, for each gene, we calculated the average ribosome density (ARD), the ratio of normalized ribo-seq footprint counts to normalized RNA-seq fragment counts, a parameter that largely reflects initiation rate (37–39). The genes were rank-ordered based on ARD and divided into octiles of 214 genes.

Analysis of mRNA secondary structure

Regions around the annotated start codon of each analyzed gene were folded in silico and the average pairing probability for each nucleotide (nt) position was calculated. The ViennaRNA Package program “RNAfold" (40), with default settings, was used to calculate the pairing probabilities for each position. These pairing probabilities were then averaged over each set of genes.

Comparative analysis of E. coli and B. subtilis

Published ribosome profiling (ribo-seq and RNA-seq) datasets from wild-type E. coli (GSE88725; (41)) and B. subtilis (GSE50870; (42)) cells grown under optimal conditions were similarly analyzed. We aligned and enumerated gene counts for E. coli and B. subtilis using the reference genome and annotations of NCBI NC 000913.3 and NCBI NC 000964.3, respectively. Analyzed gene sets contained 1440 genes (180 genes per octile) and 1392 genes (174 genes per octile), respectively.

LacZ reporter gene experiments

DNA fragments, each containing a derivative of the P_ant promoter (43) immediately upstream of the 5′ portion of gene 32 of bacteriophage T4 (−21 to 13, where zero corresponds to the first nt of the start codon), without or with a basepair substitution, were generated by cassette mutagenesis (44), using oligonucleotides listed in Supplementary Table S1. Fragments were digested with EcoRI and BamHI, and cloned into the same sites of pRS552 (45). The resulting gene 32-lacZ fusions were transferred to λRS45 by homologous recombination (45) and integrated into the chromosome of E. coli strains BW25113 (46), CSH142 (47) and KLF3027 (48). A PCR-based method was used to confirm that each of the resulting strains contained a single prophage (49). β-Galactosidase activity was measured as described previously (50).

GFP reporter gene experiments

The expression vector pSCH710, containing an IPTG-inducible σ^ABfr-dependent promoter, was constructed in several steps. Overlap extension PCR (OE-PCR; (51)) was used to generate a DNA fragment containing the ompA promoter (P_ompA) of F. johnsoniae (17,19) flanked by the operator sequences lacO3 and lacO1 of E. coli. This product was cloned into pGEM-T easy (Promega), generating pSCH685. Another chimeric DNA fragment containing the promoter for rplC (P_rplC) from F. hibernum (18) upstream of the lacI gene of E. coli was generated by OE-PCR and cloned into pGEM-T easy, creating pSCH696. The Sac II–Sph I fragment of pSCH696 was then moved into the same sites of pSCH685, generating pSCH697. Finally, the Kpn I–Bam HI fragment of pSCH697, containing P_rplC-lacI and lacO3-P_ompA-lacO1 in divergent orientations, was inserted into the same sites of pFj29 (17) to create pSCH710.

For TIR mutational analysis in F. johnsoniae, a ∼200 bp DNA fragment encompassing the TIR region of the EF-Tu gene (Fjoh_1936; chromosome region 2 238 784–2 239 021) was amplified from the F. johnsoniae chromosome and cloned into pDW01, a variant of pSCH710 in which the start codon of gfp is replaced by a Xho I site and six glycine codons. The resulting plasmid, pZM100, contains the 5′ portion of the EF-Tu gene translationally fused to the gfp reporter. Mutations in the TIR were introduced by Gibson Assembly (52), using the primers listed in Supplementary Table S1. Plasmids were moved into F. johnsoniae UW101 via tri-parental mating (28). For each strain, cells from overnight cultures were diluted 100-fold into CYE medium (3 ml) with erythromycin (100 μg/ml), and with or without IPTG (1 mM). Cultures were grown at 30°C for 6 h with shaking. Cells were pelleted, washed twice with PBS [137 mM NaCl, 2.7 mM KCl, 10 mM Na₂HPO₄, 1.8 mM K₂HPO₄ (pH 7.4)] and diluted in the same buffer to an OD₆₀₀ of 0.4. For each cell suspension, fluorescence was measured with a Fluorolog-3 (HORIBA) spectrofluorometer, using an excitation wavelength of 481 nm and emission wavelength of 507 nm. Background, determined by measuring fluorescence of untransformed F. johnsoniae grown in parallel, was subtracted in each case. Reported values, Relative Fluorescence Units (RFU), represent fluorescence intensity per OD₆₀₀.

Ribosome profiling

F. johnsoniae strain UW101 was grown at optimum temperature to mid-logarithmic phase in rich CYE medium. Cells were rapidly chilled to halt translation, sedimented, and lysed. Total RNA and ribosome-protected mRNA fragments were isolated in parallel from the lysate, and corresponding cDNA libraries were prepared and subjected to high-throughput sequencing. Replicate-to-replicate comparisons of read counts per gene showed that the RNA-seq and ribo-seq datasets were reproducible, with Spearman's correlation coefficients (r) in all cases greater than 0.98 and 0.95, respectively (Supplementary Figure S1). For this study, we chose the most highly transcribed genes (top third; 1712 genes) as a representative set for most analyses.

RNA-seq data indicates that most transcripts are leadered

Leaderless mRNAs, characterized by an AUG start codon at (or nearly at) the 5′ end, are prevalent in certain Bacteria and Archaea (53–55). As the the absence of SD sequences could stem from the absence of leaders in F. johnsoniae, we used Rockhopper (version 2.03) (56) to predict transcriptional start sites (TSSs) from the RNA-seq data. Of the 1712 genes in our representative set, TSSs were predicted for 1532 (Supplementary Table S2). A sequence resembling the –7 promoter element characteristic of the Bacteroidetes was obvious just upstream from many of the predicted TSSs. We extracted the 50 nt immediately upstream of the predicted TSSs and used BioProspector (57) to search for promoters. 434 putative –7 elements were identified. For each, the distance to the TSS was calculated, yielding the histogram of Figure 1A. Many of the sequences have 4 or 5 nt between the –7 element and the start site, consistent with bona fide promoters (18,20). These 241 sequences (deemed probable promoters) were aligned with respect to the –7 element, and the nucleotide frequencies at each position were plotted (Figure 1B). Six nucleotides of the –7 element (T-12, A-11, T-8, T-7, T-6, G-5) are nearly invariant in this set of sequences (Figure 1B and C), consistent with previous work that defined the consensus TANNTTTG (17,18,20). As the occurrence of G at positions –10 and –9 is virtually zero, the –7 element consensus should be revised to TAHHTTTG (where H represents any nucleotide but G). A sequence motif corresponding to the previously described –33 element is much less evident. In fact, AT-rich motifs centered at –18 and –28 (consecutive turns of the DNA helix from the –7 element) are more pronounced than any motif near –33 (Figure 1C). Finally, for these genes with a probable promoter, mRNA leader length was calculated (Supplementary Table S2). Most (87%) of the predicted transcripts have a leader of more than 10 nt, and the overall average leader length is 33 nt. Thus, leaderless mRNA translation cannot explain the absence of SD elements in F. johnsoniae.

Translation efficiency is tuned in part by mRNA sequence near the start codon

We calculated average ribosome density (ARD) for our representative set of 1712 genes in F. johnsoniae (Supplementary Table S2). A number of studies indicate that ARD (also termed ‘translation efficiency’ or ‘TE’) largely reflects the rate of initiation (37–39). Genes were rank-ordered by ARD and divided into eight groups, where octile 1 represents highest ARD and octile 8 represents lowest. Features of the TIR were then compared from octile to octile, in an effort to identify those that impact initiation.

We looked at nucleotide frequencies at positions –30 to +30 of the TIR (where 0 corresponds to nt 1 of the start codon). A motif indicated by overrepresentation of adenines from –14 to –11 was evident (Figure 2A), as noted previously (5,58). Enrichment of As in this region correlated with ARD (Supplementary Figure S2A), with significantly higher frequencies of A at positions –12 and –13 in octile-1 genes compared to octile-8 genes (Figure 2B). This region of mRNA is predicted based on structural studies to lie in close proximity to ribosomal protein S1 when the start codon occupies the P site (59–62). Adenine frequencies at positions 0, –3 and –6 were also significantly higher in octile-1 versus octile-8 genes and generally correlated with ARD across all genes analyzed (Figure 2B, Supplementary Figure S2A, Figure 3). The frequency of U at positions 0 and –3 was significantly lower in octile-1 genes compared to octile-8 genes (Supplementary Figure S3A). Frequencies of C and G are generally lower across the TIR (Figure 2A), reflecting the low G+C content of the genome (34%) (22), and showed virtually no significant differences in octile-1 versus octile-eight genes (Supplementary Figure S3B and C). These data suggest that translation initiation is tuned in part by the first nucleotide of the start codon (position 0) and nucleotides upstream (positions –3, –6, –12, –13), with adenines in all cases acting as positive determinants.

For comparison, we similarly analyzed data from E. coli and B. subtilis to assess which nucleotides in the TIR influence initiation rate (Figure 2C and D, Supplementary Figures S4 and S5, Supplementary Tables S3 and S4). All-octile plots showed a prominent peak of guanines centered at -9 and –11 for E. coli and B. subtilis, respectively, indicative of the Shine–Dalgarno element (Figure 2C, Supplementary Figure S5). Yet comparisons of octile-1 versus octile-8 genes showed little or no significant difference in G enrichment in this region (Supplementary Figures S4 and S5). This is consistent with evidence that the strength of the SD is a poor predictor of initiation rate for endogenous genes in both E. coli and B. subtilis (39,63,64). For E. coli, adenines at positions 0, –3 and –6 were significantly more enriched in octile-1 versus octile-8 genes (Figure 2D), reminiscent of the F. johnsoniae case. Moreover, at these positions, A enrichment trends with ARD across all octiles (Figure 3, Supplementary Figure S2B). For B. subtilis, sequence differences between octile-1 and octile-8 genes were less evident, with the most compelling difference being higher enrichment of A at position 0 of octile-1 genes (Supplementary Figure S5). These data suggest that, in all three organisms, the first nucleotide of the start codon helps tune translation initiation. In F. johnsoniae and E. coli, nucleotide identity at positions –3 and –6 also matters, with adenines at both positions stimulating translation.

To test the functional importance of nucleotides at positions -3 and -6, we used reporter gene fusions. In E. coli, we fused the TIR of gene 32 of phage T4 to lacZ, and integrated the fusion in single-copy into the chromosome. Analogous strains with the TIR mutated at position –3 or –6 were also made, and β-galactosidase (β-gal) activity of all strains was quantified (Figure 4A). The control (WT) TIR directed the highest level of β-gal production. Substitution of A-3 with either pyrimidine reduced the level of β-gal by ∼40%, whereas a smaller decrease (∼10%) was seen for A-3G. Mutation of A-6 was more deleterious, with decreases of ∼70%, 70% and 40% observed for A-6C, A-6U and A-6G variants, respectively.

In F. johnsoniae, we analyzed the TIR of gene Fjoh_1936, which encodes EF-Tu. Various derivatives of this TIR were fused to gfp, downstream of an IPTG-inducible promoter, on a plasmid shuttle vector. These plasmids were moved by conjugation into F. johnsoniae, and GFP fluorescence of the resulting strains was quantified (Figure 4B). Cells harboring the control (WT) construct and grown in the presence of IPTG gave the highest level of fluorescence. Substitution of A-3 with C, U, or G reduced GFP by 15, 30 and 70%, respectively. In this case, A-3G was particularly detrimental, suggesting an effect of the mutation beyond the loss of A, such as altered TIR folding. Mutation of A-6 also decreased GFP—by ∼15% for C/U and by 40% for G. Cells grown in the absence of IPTG exhibited much lower levels of fluorescence, yet effects of the TIR mutations were virtually identical (Supplementary Figure S6). These data show that both A-3 and A-6 can act as positive determinants of translation initiation in E. coli and in F. johnsoniae, as predicted from the ribosome profiling data.

The TIR of Fjoh_1936 has a run of adenines from –14 to –10, so we also targeted this region. Single substitutions of A-12 or A-13 reduced GFP production by ∼30%, and various double mutations reduced translation further (by 50–70%). Quadruple mutations replacing nt –14 to –11 with UCCC or UUUU reduced GFP by 90% and 60%, respectively. Why these mutations had differential effects remains unclear. But the fact that the latter mutation is no more detrimental than various double mutations suggests that the central adenines (A-12, A-13) are most critical, in line with the profiling data. Importantly, these data show that A-13 and A-12 can stimulate translation in F. johnsoniae.

Roles of A-3 and A-6 in E. coli do not depend on the β-hairpin of ribosomal protein uS7

Cryo-EM structures of 30S preinitiation complexes (PICs) containing fMet-tRNA (PIC-2C, PIC-3, PIC-4, PIC-III) show a similar conformation of mRNA (65). In these complexes, mRNA nucleotides -3 to -6 contact ribosomal protein uS7, with the turn of its conserved β-hairpin packing against the A-3 base (Supplementary Figure S7). To test whether this interaction is responsible for discrimination of A-3, we moved each of the TIR-lacZ fusions into a strain of E. coli (KLF3027) in which the β-hairpin of uS7 is truncated (i.e. residues R77-Y84 are replaced with two glycines) (48) and into the isogenic control strain (CSH142). Essentially the same pattern of β-galactosidase activity was seen for the fusions, regardless of the genetic background (BW25113, CSH142 or KLF3027) (Supplementary Figure S8, Figure 4A). These data argue against direct recognition of A-3 by the ribosome.

Translation efficiency is tuned in part by mRNA secondary structure

A growing body of evidence indicates that TIR secondary structure plays a key role in initiation, with structured mRNA elements usually inhibiting the process (66–70). To assess the impact of structure, we folded each TIR in silico and then calculated the average probability of pairing for each nucleotide position. In F. johnsoniae, nucleotides in the vicinity of the start codon showed a reduced propensity for pairing, and this region of lowered structure (–20 to 20) surrounds the start codon in a fairly symmetrical fashion (Figure 5). Pairing probabilities for these nucleotides tended to be lower for octile-1 versus octile-8 genes, with P values near or beyond the Bonferroni-corrected significance threshold at several positions just upstream of the start codon. Similar trends were seen for E. coli and B. subtilis (Supplementary Figure S9). These data provide evidence that, in all three organisms, tuning of translation initiation is controlled in part by secondary structure within the TIR. One caveat here is that we considered only local folding (–100 to +100), based on in silico predictions, and hence our results may underestimate the role of secondary structure.

Underrepresentation of AUG trinucleotides other than the start codon in the TIR

SD-ASD pairing is believed to facilitate start codon selection. Given this, it follows that an alternative mechanism(s) must be in place in organisms such as F. johnsoniae to enable start codon recognition. One simple way F. johnsoniae could reduce spurious initiation events would be to eliminate accessible AUG trinucleotides other than the start codon in each TIR. To investigate this idea, we quantified the number of AUG trinucleotides in the vicinity of the start codon relative to trinucleotides of the same base composition—GAU, GUA and AGU. For controls, we computed the analogous parameters at the gene midpoint.

In F. johnsoniae, the trinucleotide AUG is clearly underrepresented immediately upstream (−21 to −1) and downstream (4–24) of the start codon (Table 1). AUG/GAU ratios are ∼0.6 in the TIR windows, compared to ∼0.9 at midgene, suggesting that AUG trinucleotides in the vicinity of the start codon have been selected against in F. johnsoniae. We similarly analyzed E. coli and B. subtilis and found AUG trinucleotides to be more prevalent in the TIR (Table 1). In E. coli, the AUG/GAU ratios in the TIR (0.778, 0.757) are smaller than in the midgene windows (0.905, 0.901), suggesting some degree of underrepresentation. For B. subtilis, similar trends are observed, although AUG trinucleotides are more abundant in all windows. Qualitatively similar reductions in relative AUG frequency between TIR and midgene were found when normalizing to GUA and AGU. However, as these two trinucleotides contain the dinucleotides UA and AG, which are also present in stop codons, they are actually suppressed in the midgene windows, making these results harder to interpret.

We then extended this line of inquiry to include other Bacteroidetes, Proteobacteria and Firmicutes. We chose 10–12 representatives of each phylum (Supplementary Table S1) and conducted the same analysis as described above. For each phylum, we obtained P values using a one-sample t-test to assess any difference from one. To compare inter-phyla differences, we also performed a two-sample t-test for each combination of the three phyla. Furthermore, we adjusted P values with a Bonferroni correction for multiple test comparisons.

The AUG/GAU ratios were significantly different from 1 (and <1) (highlighted in blue) in the upstream window of the TIR for all three phyla (Figure 6). The ratio median was lowest for Bacteroidetes (0.56), intermediate for Proteobacteria (0.65), and largest for Firmicutes (0.75). In the downstream TIR window, the Bacteroidetes and Proteobacteria had AUG/GAU ratios significantly different from 1, and <1, which cannot be said of the Firmicutes. For both TIR windows, there was a significant difference in AUG/GAU between the Bacteroidetes and Firmicutes. In the midgene windows, AUG/GAU ratios for all phyla were statistically equivalent to 1 (Figure 6). Again, the data for the other two normalizing trinucleotides (Supplementary Figure S10) is harder to interpret due to their midgene suppression. Nevertheless, AUG frequency is significantly distinct from 1 (and <1) in the upstream TIR window for both normalizations and in the downstream TIR window for the AUG/AGU ratio in the Bacteroidetes. Collectively, these data show that the trinucleotide AUG tends to be underrepresented in the neighborhood of the start codon across all representatives of the Bacteroidetes. There is some indication of the same trend in the other phyla, particularly Proteobacteria. However, the degree of AUG underrepresentation in the TIR is largest in the case of the Bacteroidetes.

Ribo-seq coverage predicts protein stoichiometries in F. johnsoniae

Weissman and colleagues have shown that relative protein levels in E. coli can be predicted with remarkable accuracy based on ribo-seq data (39). Because, to our knowledge, ours are the first ribosome profiling data for the Bacteroidetes, we calculated the relative rates of protein production in F. johnsoniae, as ribo-seq coverage per gene length, for all genes (Supplementary Table S2). These values estimate the stoichiometries of protein chains within complexes, such as those of the gliding-motility machinery, which includes a type IX secretion system (Supplementary Figure S11).