Diversity of Translation Initiation Mechanisms across Bacterial Species Is Driven by Environmental Conditions and Growth Demands.

The Shine-Dalgarno (SD) sequence motif is frequently found upstream of protein coding genes and is thought to be the dominant mechanism of translation initiation used by bacteria. Experimental studies have shown that the SD sequence facilitates start codon recognition and enhances translation initiation by directly interacting with the highly conserved anti-SD sequence on the 30S ribosomal subunit. However, the proportion of SD-led genes within a genome varies across species and the factors governing this variation in translation initiation mechanisms remain largely unknown. Here, we conduct a phylogenetically informed analysis and find that species capable of rapid growth contain a higher proportion of SD-led genes throughout their genomes. We show that SD sequence utilization covaries with a suite of genomic features that are important for efficient translation initiation and elongation. In addition to these endogenous genomic factors, we further show that exogenous environmental factors may influence the evolution of translation initiation mechanisms by finding that thermophilic species contain significantly more SD-led genes than mesophiles. Our results demonstrate that variation in translation initiation mechanisms across bacterial species is predictable and is a consequence of differential life-history strategies related to maximum growth rate and environmental-specific constraints.

Introduction

Translation of a given messenger-RNA (mRNA) into functional protein relies on the ability of the translational apparatus to recognize the proper start codon. Bacteria have evolved several distinct mechanisms to discriminate between potential start codons, with the Shine–Dalgarno (SD) mechanism being the most well-studied (Shine and Dalgarno 1974; Nakagawa et al. 2010). Variants of the SD sequence are frequently found upstream of bacterial start codons where they facilitate translation initiation by hybridizing with the complementary anti-SD (aSD) sequence on the 16S rRNA of the small ribosomal subunit (fig. 1).

. 1.

Sequence entropy quantifies genome-wide SD sequence utilization. (A) Illustration of the anti-Shine–Dalgarno(aSD)::Shine–Dalgarno(SD) sequence mechanism of translation initiation. (B) Representative sequence logos of the 5′ upstream region of all annotated coding sequences for individual genomes displays heterogeneity in sequence entropy within and between species. (C) Illustration of the ΔI metric for Caulobacter crescentus as an example.

For a given gene within an organism, the structural accessibility of the SD sequence, the thermodynamic binding potential between the SD sequence and the aSD sequence, and the exact positioning of the SD sequence relative to the start codon, are all features that can be predictably tuned in order to modulate the translation initiation rate of downstream genes (Barrick et al. 1994; de Smit and van Duin 1994; Vimberg et al. 2007; Salis et al. 2009; Devaraj and Fredrick 2010; Na et al. 2010; Kosuri et al. 2013; Espah Borujeni et al. 2014; Bonde et al. 2016; Espah Borujeni and Salis 2016; Hockenberry et al. 2017). In particular, researchers have shown that transcripts with a strong SD sequence are translated at higher rates resulting in more protein being produced per mRNA, a fact which is particularly important for the design of recombinant protein expression systems (Salis et al. 2009; Na et al. 2010; Bonde et al. 2016). Nevertheless, there are several SD sequence-independent mechanisms that operate in bacteria including leaderless translation and RPS1-mediated translation of unstructured mRNA sequences (Komarova et al. 2005; Chang et al. 2006; Gu et al. 2010; Scharff et al. 2011; Zheng et al. 2011; Keller et al. 2012; Barendt et al. 2013; Cortes et al. 2013; Duval et al. 2013; Kramer et al. 2014; Shell et al. 2015). Recent research also suggests that mechanisms traditionally associated with eukaryotic species such as translational scanning and internal ribosome entry sites may operate in bacterial systems (Colussi et al. 2015; Yamamoto et al. 2016).

The aSD sequence is highly conserved throughout the bacterial domain (though notable exceptions exist) and stronger SD sequence interactions are associated with increased translation efficiency across a wide array of species (Chen et al. 1994; Osada et al. 1999; Sakai et al. 2001; Ma et al. 2002; Starmer et al. 2006; Lim et al. 2012; Omotajo et al. 2015; Hockenberry et al. 2017). We therefore might expect that most species would utilize the SD sequence mechanism to a similar degree. However, there is a wide-diversity of SD sequence utilization between different species. For instance, roughly 90% of Bacillus subtilis genes are preceded by a SD sequence, whereas for Caulobacter crescentus, the comparable number is closer to 50% (Chang et al. 2006; Starmer et al. 2006; Nakagawa et al. 2010; Schrader et al. 2014). SD-like sequence motifs are also underrepresented within the coding sequences of most bacteria—possibly reflecting their role in translational pausing and/or erroneous initiation—and like the diversity of SD sequence utilization, the degree of this underrepresentation varies from species-to-species (Li et al. 2012; Diwan and Agashe 2016; Mohammad et al. 2016; Yang et al. 2016).

Cross-species variation in translation initiation mechanisms may impact genetic isolation and transfer of genetic material, and quantifying the source and extent of variation may prove useful in identifying important genes in a genome or microbial community(Krisko et al. 2014; Omotajo et al. 2015). Further, the synthetic biology community is increasingly targeting both translation-system engineering and biotechnology applications involving less well-studied microbial species (Tauer et al. 2014; Markley et al. 2015; Orelle et al. 2015; Guiziou et al. 2016; Weinstock et al. 2016; Yi et al. 2016). A better understanding of the factors shaping the utilization of different translation initiation mechanisms may ultimately aid in the design of synthetic gene constructs.

Here, we conduct a phylogenetic comparative analysis in order to isolate independent evolutionary events and determine whether any endogenous or exogenous factors are predictive of genome-wide SD sequence utilization. We develop a metric that captures position-dependent sequence preferences within the translation initiation region of a given genome, and demonstrate a strong link between SD sequence utilization and minimum doubling times for 187 species. Furthermore, in a database of 613 phylogenetically diverse bacterial species, we show that genome-wide variation in SD sequence utilization covaries along-side a number of genomic features that are indicative of rapid-growth and efficient translation. Finally, we investigate exogenous environmental constraints and show that SD sequence utilization varies according to optimal growth temperatures.

Results

Quantifying Genome-Wide SD Sequence Utilization

Several techniques have been previously developed to quantify the overall utilization of the aSD::SD mechanism within a given species. In motif-based methods, researchers predefine subsequences closely related to the canonical SD sequence and search a sequence window upstream of each protein coding gene within a given genome to determine the fraction of genes that are preceded by a SD motif (Chang et al. 2006; Omotajo et al. 2015). Similarly, in aSD sequence complementarity based methods, researchers predefine a window upstream of the start codon to consider for each gene, a putative aSD sequence, and a hybridization energy threshold for determining whether a gene is SD-led or not (Osada et al. 1999; Sakai et al. 2001; Ma et al. 2002; Starmer et al. 2006; Nakagawa et al. 2010).

Both of these metrics rely on critical assumptions that may not hold when applied across large sets of phylogenetically diverse organisms. First, they carry an implicit assumption that a SD sequence, regardless of its location relative to the start codon, has the same impact on translation initiation. However, experimental approaches have shown that spacing between the SD sequence and start codon can have dramatic effects on translation initiation rates (Vimberg et al. 2007; Salis et al. 2009; Devaraj and Fredrick 2010; Hockenberry et al. 2017). Second, both methods rest on a dichotomy between SD-led and non-SD-led genes. Although this simplification is useful for describing the phenomenon, an abundance of research has shown that a spectrum of sequence complementarity affects translation initiation in a continuous manner (Vimberg et al. 2007; Salis et al. 2009). Third, bacterial genomes span a range of GC contents, and previous research has shown that it is possible to account for this bias by comparing the proportion of SD-led genes in a genome to an appropriate null model expectation (Nakagawa et al. 2010). We define the following term to summarize SD sequence utilization using the SD-motif based method by: where fSD is the fraction of genes within a genome classified as SD-led and is the expected fraction of SD-led genes derived from repeating this calculation for 500 nucleotide shuffled “genomes” (where 30nts of the 5′-UTR for each gene are randomly permuted, see Materials and Methods). We similarly define , where SD-led genes are defined via hybridization of the putative aSD sequence using a threshold binding energy value of −4.5 kcal/mol (see Materials and Methods).

Sequence Entropy and Its Relationship with SD Sequence Utilization

We sought a complementary approach that would allow us to investigate hundreds of diverse genomes without having to a priori define either an aSD sequence or SD sequence motifs. To measure arbitrary position-specific sequence signals, we extract the 5′ upstream sequences from all annotated protein coding sequences (see Materials and Methods) and sum the “information content” at each position within the region where SD motifs are expected to occur (−20 to −4 relative to the start codon): where p, is the empirical frequency of base k at position i. The right side of this equation is the entropy (calculated from empirical base frequencies at each position) while the component scales the data in line with what is typically represented in sequence logo plots such that individual nucleotide positions have a maximum value of 2 when fully conserved and a minimum of 0 when fully random. We again repeat this process for 500 shuffled “genomes” (as described above) and compare the sequence information from the actual genome to this null expectation:

For a more intuitive interpretation, figure 1 shows sequence logo plots of the 5′-UTRs for three genomes to highlight the variation in sequence preferences between species. Iobs is simply the sum across each position of these nucleotide “stacks.” To account for uneven nucleotide usage, we calculate ΔI by comparing this value to nucleotide shuffled controls (fig. 1). In supplementary figure S1, Supplementary Material online, we simulate sequences to show that ΔI takes larger values when (all else being equal): 1) a larger fraction of simulated genes contain a particular sequence motif, 2) the enriched sequence motif is longer, and 3) the enriched sequence motif is more strictly defined in its position.

We compiled a data set of 613 bacterial species, unique at the genus level, with annotated genome-sequences as well as a previously constructed high-quality phylogenetic tree describing their relatedness (Hug et al. 2016) (see Materials and Methods). In figure 2, we show that while summary methods based on SD-motif and aSD sequence complementarity (ΔfSD and , respectively) are strongly related for a large set of diverse species, this association is not perfect and there is a distinct change in the slope that occurs for the Firmicutes phylum.

. 2.

Relationship between ΔI and existing metrics of SD sequence utilization. (A) Comparison between different ways of summarizing SD sequence utilization; each data point represents a single genome. On the left, we show the relationship between SD motif and aSD sequence complementarity based methods (ΔfSD and ). On the right, we compare ΔI and . The four largest phyla are color-coded according to the legend. Arrows highlight phyla with “anomalous” patterns. (B) Phylogenetic tree illustrating variation in SD sequence utilization across species according to the ΔI metric (indicated by bar plots on concentric rings).

In principle, ΔI may quantify a variety of position-specific sequence signals that may or may not be related to the SD sequences or translation initiation. In practice, we observe that this metric correlates strongly with both methods traditionally used to describe SD sequence utilization (fig. 2 and supplementary fig. S2, Supplementary Material online). However, in the Bacteroidetes phylum, we detect significant variation in ΔI without any apparent variation in either of the other two metrics. These findings are consistent with prior research that identified changes in the aSD sequence region of the 16 S rRNA sequence within this phylum (Lim et al. 2012), and indicate that Bacteroidetes may nevertheless contain position-specific translation initiation sequences. The sequence logos for representative Bacteroidetes species showed A/T rich UTRs (supplementary fig. S3, Supplementary Material online), possibly indicative of translation initiation mediated by ribosomal protein S1. However, we stress that general nucleotide bias is removed in our calculation of ΔI via shuffling and the A/T richness in these species instead has a fairly strong position-specific signal. The fact that the ΔI metric quantifies position-specific translation initiation signals for Bacteroidetes allows us to incorporate them into future analyses (fig. 2, red data points) though we note that much remains to be understood about possible changes to the identity of the SD and aSD sequences in this phyla.

Due to the strong relationship between ΔI and explicit SD sequence methods, we refer to ΔI as quantifying SD sequence utilization throughout the remainder of this manuscript and make explicit note of differences between SD sequence quantification methods when necessary. Consistent with prior research (Nakagawa et al. 2010), we show that SD sequence utilization according to the ΔI metric varies considerably across species while showing broad phyla-specific patterns that should be accounted for when performing statistical tests (fig. 2, bar heights). We additionally found that the UTRs from ribosomal protein coding genes contained more position-dependent sequence information—using either Iobs or ΔI, for nearly all species in our data set—compared with genome-wide UTRs (supplementary fig. S4, Supplementary Material online). This observation is consistent with theory emphasizing the overall importance of efficient translation for the most highly expressed proteins. We also tested whether ΔI varies according to genome size and found no significant association with either genome length or the number of annotated protein coding genes (P = 0.54 and 0.84, respectively).

Translation Initiation and Organismal Growth Demands

In prior research, Vieira-Silva et al. (2010) curated a list of minimum doubling times from the literature for a large number of bacterial species (Roller et al. 2016). Organisms that are capable of rapid growth have high protein production demands during these periods and there are a number of regulatory points that can be bottlenecks for this process. Meeting high translational demands associated with rapid growth requires coordination of a number of processes, and Vieira-Silva et al. (2010) showed that increasing numbers of rRNA and tRNA genes, and increasing codon usage biases among ribosomal mRNAs in individual genomes were all partially predictive of the minimum doubling times of individual species. Subsequently, Yang et al. (2016) showed that the extent of anti-SD sequence binding across coding sequences was also predictive of minimum doubling times reflecting the possible influence of internal SD-like sequences on translational pausing and/or spurious translation initiation events (Yang et al. 2016).

At the individual gene level, translation initiation is an important control point and likely to be rate-limiting for the production of most proteins. Experimental evidence largely supports this conjecture, with alterations to SD sequence utilization, mRNA structural availability, and start codon identities predictably tuning protein production (Kudla et al. 2009; Salis et al. 2009; Goodman et al. 2013; Kosuri et al. 2013; Espah Borujeni et al. 2014, 2016; Espah Borujeni and Salis 2016; Hecht et al. 2017). We thus reasoned that efficient translation initiation across all or subsets of genes may similarly play an important role in meeting protein production demands imposed by rapid growth rates. In addition to genome-wide SD sequence utilization, we also quantified the percentage “ATG” start codons in a genome as well as the average difference in mRNA folding energy surrounding the start codons relative to a control region internal to the gene (see Materials and Methods).

We first replicated several of the findings of Vieira-Silva et al. (2010) and Yang et al. (2016) using Phylogenetically Generalized Least Squares regression (Revell 2010) to account for the lack of independence in species (see Materials and Methods). We verified that rRNA gene counts, tRNA gene counts, a metric of relative codon usage bias (based off the “effective number of codons” [ΔENC], see Materials and Methods), and a metric of internal SD-like sequence binding (see Materials and Methods) are all significantly predictive of minimum doubling times after controlling for phylogenetic effects (F-test, P < 0.002 for all cases, table 1).

Table 1.

Contribution of Several Factors for Predicting Minimum Doubling Times.

Model for Min. Doubling Time	R2	Pagel’s λ [95% CI]	|ΔR^2|
Full model	0.35***	0.91 [0.80, 0.96]	–
ΔENC	0.17***	0.96 [0.92, 0.99]	0.15
ΔI	0.11***	0.97 [0.93, 0.99]	0.12
mRNA folding	0.08***	0.98 [0.95, 0.99]	0.04
Internal SD-like	0.08***	0.98 [0.95, 0.99]	0.03
16S gene counts	0.06**	0.98 [0.95, 0.99]	0.01
tRNA gene counts	0.06***	0.98 [0.95, 0.99]	0.01
ATG start %	0.02	0.98 [0.95, 0.99]	<0.01

Note.—The left column indicates individual variables that we considered for predicting minimum doubling times with the full multivariate model listed at the top. R2 column illustrates the overall goodness-of-fit for individual factors (*** indicates P < 0.001, ** indicates P < 0.01). Pagel’s λ is the fitted phylogenetic signal parameter, which we show with 95% confidence intervals in brackets. The right column illustrates the change in goodness-of-fit from a model that includes all predictors to one that excludes only the variable in the given row.

Next, we turned to metrics related to translation initiation. We found that ΔI calculated over all genes within a genome significantly correlates with minimum doubling times in this set of species (P < 10−5), showing the second strongest correlation of any individual trait that we considered (table 1 and fig. 3). Our metric of mRNA structural stability surrounding the start codons—“mRNA folding” in table 1—similarly showed a significant relationship such that genomes that, on an average, have the weakest mRNA structure surrounding the start codon are more capable of rapid growth. By contrast, the proportion of protein coding genes containing an ATG start codon is not significantly correlated with minimum doubling times (P = 0.056), but showed the expected slope indicative of faster growth for genomes with more ATG start codons.

. 3.

Relationship between SD sequence utilization and organismal growth. (A) ΔI is significantly correlated with minimum observed doubling times for 187 bacterial species. (B) Visualization of the full model listed in table 1 depicting a strong relationship between observed and predicted minimum doubling times. In both plots, individual species data points are colored according to phyla as in figure 2.

In order to test the robustness of these findings and to assess overall predictability of minimum doubling times from these features, we constructed a multivariable Phylogenetic Generalized Least Squares regression model that combines all of the above factors, and found that only relative codon usage biases (ΔENC), SD sequence utilization (ΔI), and mRNA folding energy had statistically significant coefficients in this full model (P < 0.01). Overall, a model containing all factors resulted in R2 = 0.35 (P < 10−13, fig. 3), whereas a more parsimonious model containing only the three factors with statistically significant coefficients resulted in R2 = 0.31 (P < 10−14). Removing either of these features from the full model reduces its predictive power as illustrated in the right column of table 1, which further emphasizes the large contribution of ΔENC and ΔI in particular.

In order to compare our work with prior research, we also conducted a phylogenetically agnostic linear regression model using all of these factors, which yielded R2 = 0.6 (P < 10−15)—though we caution that ignoring the effects of shared ancestry will substantially bias statistical analyses. We further generated the same data as in table 1 using as a metric of SD sequence utilization and found largely similar results with slightly less predictive power overall: R2 = 0.29 for the full model and R2 = 0.09 as an individual predictor compared with 0.35 and 0.11 using ΔI (supplementary table S1, Supplementary Material online). We repeated these analyses excluding species from Bacteroidetes phylum, but this restriction removed only five species and correlations for both ΔI and independently and in the full model remained essentially unchanged.

Interestingly, we looked at the effect of using only using ribosomal protein coding genes to calculate ΔI and found that this metric performed worse at predicting observed minimum doubling times individually (R2 = 0.04) and as part of the full model (R2 = 0.27) (supplementary table S2, Supplementary Material online). The metrics of codon usage bias and internal SD-like sequence usage that we used here are both calculated from the difference between values for ribosomal protein coding genes and all genes within a genome, thus representing relative differences in selection strength acting specifically on ribosomal protein coding genes. However, we observed little predictive power when using a relative difference metric between ribosomal protein coding genes and all genes within a genome—using either Iobs or ΔI to predict minimum doubling times (R2 ≈ 0). Instead, genome-wide calculations of ΔI consistently provided the highest predictive power.

As a final caveat for this data, we do note that the list of species compiled by Vieira-Silva et al. (2010) is phylogenetically biased toward Proteobacteria, contains only culturable species, and is distinct in its composition from the data set analyzed in figure 2.

Relationship between SD Sequence Utilization and Other Translation Efficiency-Associated Traits

Since a coordinated effort between multiple processes is required to maximize protein production in experimental systems, we reasoned that the various genome-wide traits associated with efficient translation are likely to covary with one another across species. In order to test this hypothesis, we assessed the correlation between different definitions of SD sequence utilization and all of the alternative traits listed in table 1 via Phylogenetic Generalized Least Squares regression. In figure 4, we show the results of this analysis, finding that in all cases where a pair of traits is significantly correlated, the correlation is positive (note that the color bar corresponds to R and not R2). Increasing SD sequence utilization is thus significantly associated with an increasing fraction of ATG start codons, less structured mRNA folding around the start codons, increased 16S rRNA and tRNA gene counts, and increasing codon usage bias/avoidance of internal SD-like sequences within ribosomal protein coding genes.

. 4.

SD sequence utilization covaries alongside a suite of translation-related traits and according to optimal growth temperatures. (A) Correlation matrix between listed variables used in table 1 for a set of 613 diverse bacterial species. In all instances of significant correlation, the features covary with one another in the positive direction. (B) SD sequence utilization, quantified using ΔI is significantly higher in thermophiles than in mesophiles. Box limits show 25th and 75th percentiles of the data, whiskers extend to 5th and 95th percentiles, triangles depict the means of each category and red lines highlight the median.

We next tested the overall robustness and universality of these results by independently analyzing these relationships within individual phyla. We specifically looked at the four largest phyla in our data set—Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes—and repeated the analysis from figure 4 using independent model fits. Again, we observed that nearly every significant correlation is in the positive direction (supplementary fig. S5, Supplementary Material online) with the exception of mRNA structure and SD sequence utilization in the Firmicutes phylum. When looking at relationships between variable SD sequence utilization in the Bacteroidetes phylum, ΔI has a significantly positive relationship with several other variables whereas ΔfSD and show no significant relationships apart from with one-another.

Relationship between Translation Initiation Mechanisms and Ecological Factors

Having established that genome-scale SD sequence utilization is part of a suite of traits related to differential organismal growth strategies, we last wanted to assess whether ecological factors pertaining to an organisms habitat may constrain the evolution of SD sequence utilization. The aSD:: SD sequence interaction is thought to increase translation rates by stabilizing the assembly of the translation initiation complex. Since this interaction is based on RNA base-pairing, pairing between longer sequences may be necessary in order to get an equivalent level of stabilization at higher temperatures. As we show in supplementary figure S1, Supplementary Material online, ΔI takes larger values when longer sequences are preferred so we reasoned that there might be an association between this metric and optimal growth temperatures.

Nakagawa et al. (2010) investigated this possibility, but found no association (Nakagawa et al. 2010). By contrast, our phylogenetically informed modeling approach applied to this larger data set (481 of the 613 species in our data set have high-confidence growth temperature annotations) finds that temperature constrains genome-wide SD sequence utilization. Specifically, the genomes of thermophilic species display significantly larger values of ΔI than mesophilic species (fig. 4test P = 0.002 using temperature as a fixed-effect in Phylogenetically Generalized Least Squares modeling) in line with our hypothesis. Although optimal growth temperature is a continuous variable, in practice, measurements for the vast majority of species fall into discrete categories (30 °C and 37 °C). This fact, coupled with the availability of a large database of “mesophile” and “thermophile” annotations motivated our decision to study this effect using discrete temperature categories.

Our finding here illustrates the role that ecological factors relating to growth conditions places on the evolution of genome architectures. We further tested whether there are any systematic differences in SD sequence utilization by looking at free-living/host-associated species, or pathogenic/nonpathogenic species, and found no significant effects for either (P = 0.91 and 0.44, respectively). Nevertheless, it is possible that other ecological and life-history traits that we did not consider here may play an important role in constraining the evolution and utilization of different translation initiation mechanisms.

Discussion

Our study shows a relationship between bacterial translation initiation mechanisms, life-history strategies, and environmental demands faced by individual species. Although we do not specifically address causality, we found that minimum observed doubling times and SD sequence utilization at the genome-scale are significantly correlated (fig. 3). In a larger and more diverse data set of 613 species, we further showed that SD sequence utilization predictably covaries with several other genomic and environmental features, including the number of rRNA genes and optimal growth temperatures. Taken together, our findings demonstrate that organisms with greater translational demands are likely to coevolve a common suite of genomic features that help to maximize translation during periods of rapid growth, and that SD sequence utilization is an important component of this shared genome architecture.

Our analysis throughout is performed in a manner that corrects for the confounding effects of shared ancestry between species, and our phyla specific results illustrate several critical points. First, the sign on the relationships between individual features that we observe is extremely robust, regardless of the phylum or SD sequence utilization summary statistic under consideration (fig. 4 and supplementary fig. S5, Supplementary Material online). Increasing 16S/tRNA gene counts, increasing codon usage biases/avoidance of internal SD-like sequences in ribosomal protein genes, higher ATG start codon usage, and weak mRNA structure surrounding the start codon are associated with increasing SD sequence utilization. Second, we developed a novel metric (ΔI) to measure position-specific translation initiation region sequence preferences in a manner that is independent of the anti-SD or SD sequence for a given species. We found that ΔI is highly correlated with prior estimates of SD sequence utilization, but that this relationship breaks down in the Bacteroidetes phylum, which may reflect novel position-specific sequence preferences in this lineage. Future research applying our methodology to larger data sets may allow researchers to uncover branches within phylogenetic trees where mechanistic differences in the translational apparatus—resulting in differences in the slope and/or sign on the relationships between different features—have evolved.

Despite the overall robustness of our results, there are several caveats to our study that we wish to explicitly highlight for readers. First, the lack of available growth rate data for most phyla prohibits us from looking at the relationships between genomic traits and minimum doubling times on a phyla specific level. The species for which we do have minimum doubling time estimates are biased phylogenetically and represent only a fraction of the known bacterial universe. Second, our findings throughout rely on genome annotations, which frequently rely on the presence of SD sequences to identify the 5′ ends of genes. Such annotations are known to be imperfect, and as many as ≈ 10% of genes in a given genome may have misannotated start codons (Schrader et al. 2014). However, these biases are likely to be uniform—or at least randomly distributed—throughout our data set making it unlikely that systematic differences in start codon annotation between genomes would produce spurious results of the order that we report here. Nevertheless, continued improvement in start codon annotations is an important issue and potential limitation of our study.

Our results add to the body of knowledge showing that a small number of genomic traits—that includes utilization of the SD sequence mechanism—can be used to predict variation in minimum doubling times. However, a critical question that remains unanswered is: what is the causal relationship between different genome-scale phenomena and growth rate control? We emphasize that our study is correlative, and that truly establishing causality on genome-scale patterns will require large-scale experiments that, for instance, systematically alter the SD sequences for hundreds of different genes at the same time and measure the resulting growth rates. Prior research has shown detrimental growth effects from deleting rRNA operons and tRNA genes in different microbial species (Stevenson and Schmidt 2004; Yano et al. 2013; Bloom-Ackermann et al. 2014; Samhita et al. 2014). Additionally, genome-scale engineering efforts to alter codon usage bias patterns have likewise showed increasing doubling times in these strains (Napolitano et al. 2016; Ostrov et al. 2016). Rather than suggest a single limiting feature for bacterial growth control, these results highlight that a number of features are likely necessary but insufficient on their own to increase bacterial growth rates.

We do not wish to suggest that increasing the translation initiation rates for all genes within a genome would result in a faster growing species, but we do hypothesize that decreasing the translation initiation rates for a large number of randomly chosen genes is likely to produce substantially longer doubling times. Why our results showed the strongest correlations with minimum doubling times when assessing genome-wide SD sequence utilization rather than only looking at ribosomal protein coding genes is unknown. However, we speculate that under periods of rapid growth, hundreds of genes (that include ribosomal proteins) are likely required and capable of bottle-necking bacterial growth. Calculating SD sequence utilization on a subset of such genes may produce an even stronger relationship with minimum doubling times, but the identities of these genes are likely to be organism and condition-specific.

We found that measurements of SD sequence utilization outperform more commonly known associations such as the number of rRNA genes at predicting minimum doubling times (R2 = 0.11 vs. 0.06, and increasing feature importance in the full model depicted in table 1) (Roller et al. 2016). We believe that this finding may, in part, be a consequence of the relevant evolutionary time-scales that it takes to alter different traits. Specifically, substantially changing genome-wide SD sequence utilization or codon usage biases would require hundreds or thousands of fixed mutations, respectively. We found that these two traits were the most predictive of minimum doubling times out of all traits that we considered, which is in contrast to traits such as rRNA and tRNA gene counts that can be altered rapidly by deletions/insertions. Thus, the degree of codon usage bias or SD sequence utilization within a genome may better reflect long-term historical forces acting on a species, whereas copy number variations have the potential to better reflect more rapid evolutionary changes.

Finally, we note that—like codon usage biases and in contrast to rRNA and tRNA gene counts—summary statistics based on SD sequence utilization do not require complete genome sequences and therefore may be estimated with partial genome fragments. The results and methods that we present here may thus have important applications in our understanding of novel, uncultivated genomes, environmental meta-genomic sequencing efforts, and the relationship between higher order genome traits and growth strategies (Brown et al. 2016).

Materials and Methods

Data Assembly

We first assembled a database of prokaryotic genomes from NCBI using the GBProks software (https://github.com/hyattpd/gbproks), including only “complete” genomes in our download and subsequent analysis (accessed on: March 10, 2016). From the annotated GenBank files, we excluded pseudogenes and plasmid based sequences from all subsequent analyses and proceeded to compile a data table with several traits for each genome. In addition to SD sequence utilization summary statistics described below, we applied RNAmmer to each genome in order to compile a list of ribosomal-RNA genes, and tRNAscan-SE to assemble a list of the tRNA genes (Lagesen et al. 2007; Lowe and Chan 2016).

We wrote custom scripts to calculate the fraction of annotated coding sequences that begin with “ATG,” as well as the metric of codon usage bias (ΔENC as described in Vieira-Silva et al. 2010). For this latter metric, we first parsed the gene annotations to find ribosomal protein coding genes. We next computed the relative differences in codon usage bias between ribosomal protein coding genes and the rest of the genome, whereby: where “all” and “ribo” refer to all protein coding genes and ribosomal protein coding genes, respectively. We altered the method used to calculate the “effective number of codons” or “ENC” from the one originally used by Vieira-Silva et al. (2010) to better control for GC content differences according to recent metric developed in our lab (manuscript submitted). The interpretation is the same, with larger positive values indicating more codon usage bias in ribosomal protein coding genes relative to the rest of the genome.

We calculated the metric of internal SD-like sequence occurrences as reported in in Yang et al. (2016). Briefly, we first calculate the average aSD sequence binding strength for each hexamer within a given coding sequence. As above, the final metric is then calculated as the difference in average aSD binding energies between all protein coding genes and ribosomal protein genes divided by the value for all genes. Larger value indicate fewer SD-like sequences are present in ribosomal protein coding genes relative to all genes.

Our measurement of genome-wide “mRNA folding,” is calculated as follows. For each gene greater than 150nts in length, we extract a 60-nt long region centered on the start codon (−30 to +30) and calculate the minimum free energy of this segment under ViennaRNA defaults. We next do the same thing for an internal 60-nt segment (which we chose to standardize as nucleotides +90 to +150 for all genes). For each gene, we calculate the difference of these two folding energies and for each genome, we calculate the average of these differences across all genes. Negative values indicate that regions surrounding the start codon are more structured than internal regions, and as values get more positive the region surrounding the start codon is increasingly less structured compared with internal sites.

For data on minimum doubling times, we downloaded the data table from Vieira-Silva et al. (2010), and paired each bacterial species with a complete genome from our database resulting in 187 matched species. To control for shared ancestry in subsequent analyses, we constructed a phylogenetic tree based off the rRNA sequences for this set of species. We first used RNAmmer to extract a randomly chosen 16S and 23S rRNA sequence from each genome, followed by MUSCLE (v3.8.31) on each individual rRNA to produce a multiple-sequence alignments (Edgar 2004). These were concatenated together and we conducted a partitioned analysis using RAxML to construct a final tree. We performed 100 rapid Bootstrap searches, 20 ML searches and selected the best ML tree for subsequent analysis (Stamatakis 2014).

For the larger data set, we instead relied on a previously computed high-quality phylogenetic tree published by Hug et al. (2016) (Hug et al. 2016). We used custom scripts to match entries in this tree with genomes from our complete-genome database, and pruned away all species without a high-quality match resulting in 613 bacterial species in our final data set that were used for subsequent analyses. For temperature annotations, we matched this set of 613 species to the ProTraits database using custom scripts, and restricted our analysis to species with temperature annotations exceeding a precision of 0.9 (equivalent to a FDR < 0.1) (Brbić et al. 2016).

Calculating Summary Statistics of SD Sequence Utilization

The calculation of ΔI is illustrated mathematically in the main text. Here, we only add that the calculation of the randomized sequences for all SD summary statistics is performed by first shuffling the upstream region of each gene between the region −30 to 0 (where +1 is the first base of the start codon). Having shuffled each gene in this manner, we then performed the analysis as discussed in the main text for this shuffled “genome” and repeat this calculation 500 times in order to derive null expectation for fSD, faSD<−4.5 and Iobs.

Next, we elaborate on our calculation of the other two methods for calculating SD sequence utilization. For each genome, we extract the −20 to −4 region upstream of the start codon for each gene. For fSD, we consider a gene as being SD-led if, in this defined region, any of the following motifs appear: “GGAA,” “GGAG,” “GAGG,” “AGGA,” or “AAGG.” We repeat this same process for 500 randomized “genomes” where a randomized genome is defined as noted above (with the nucleotide region from −30 to 0 for each gene shuffled on a per-gene basis) prior to motif search.

For , we perform a nearly identical procedure to the one listed above with the major difference being that instead of searching the upstream region of genes for particular motifs, we evaluate the hybridization energy between each eight nucleotide segment contained within the −20 to −4 region and the putative aSD sequence defined as 5′-ACCUCCUU-3′ using the “cofold” method of the ViennaRNA software package with default parameters. If any sequence binds at a threshold of −4.5 kcal/mol or stronger (i.e., more negative ΔG values), we consider this gene to be SD-led.

Phylogenetically Generalized Least Squares

Throughout this manuscript, we utilize Phylogenetically Generalized Least Squares regression in order to mitigate the effects that arise from shared ancestry in statistical analyses. Our Phylogenetically Generalized Least Squares analysis relies on the most common null model, which assumes a Brownian motion model of trait evolution. For all statistical analyses presented in the paper, we use the R package “caper” and perform a simultaneous maximum-likelihood estimate of Pagel’s λ, a branch length transformation, alongside the coefficients for independent variables of interest. All P values that we report come from the F-test according to these results. For temperature analysis, we assigned “mesophiles” and “thermophiles” a value of 0 and 1, respectively, and performed the equivalent fixed-effect analysis with ΔI as the dependent variable.

Data Availability and Computer Code

Data are provided as a supplementary file, Supplementary Material online and all custom scripts and code that is sufficient to perform the analysis can be found at https://github.com/adamhockenberry/SD-evolution-publication.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.