Translation termination is a crucial step in protein synthesis that, in most organisms, is triggered by three stop codons; TAA, TGA and TAG. These three stop codons are thought to be functionally equivalent in the broad sense of effective translation termination. Additional functions, such as coding for extra amino acids, effects only a tiny fraction of all codons [1], and these stop codons can be interchanged [2, 3] or even lost [4–8] without obvious functional consequences. Indeed, one of the motivations in a recent experimental study of genome-wide codon replacement in selecting to substitute all TAG stop codons in Escherichia coli, rather than making synonymous substitutions, was the rationale that synonymous “codon utilization bias has been shown to affect translation efficiency” [3] suggesting that in the author’s opinion stop codon substitution may have fewer functional consequences than synonymous substitution. Thus, at present there is broad consensus that three stop codons are functionally equivalent and interchanging stop codons is not expected to have functional or selective consequences. In that case substitutions between different stop codons should be neutral, such that the rate of evolution between stop codons should be broadly equivalent to the synonymous rate of evolution and the stop codon frequency should be governed by similar selective and mutational forces that govern nucleotide usage in synonymous sites.

The hypothesis of selective equivalence of stop codons has not been rigorously tested and, contrary to the general expectation, there are data that suggest that stop codon may not be entirely synonymous. Firstly, translation termination efficiency may be nucleotide context dependent [9–13]. Second, TAG and TGA stop codon frequencies in bacterial genomes with different GC-contents are strikingly different (see Figure 1 in [14]), such that TGA frequency increases with genomic GC-content while TAG is GC-content independent. Here, we study stop codon frequency and evolution in bacterial genomes to gain an understanding of whether or not stop codons are used indiscriminately without any fitness costs. We compare rates of stop codon evolution to the rate of synonymous evolution and apply a simple population genetics model formulated by Bulmer [15] to stop codon frequency and nucleotide content in bacterial genomes.

The relationship between stop codon frequency and G-content, or GC-content as reported previously [14], is one of the most striking and unambiguous patterns in bacterial genome composition (Figure 1). Here, we have developed a simple model that captures all of the major observations of stop codon distribution across bacterial genomes. However, as with many theoretical treatments our model necessarily makes several simplifying assumptions. First, we ignore the effects of neighboring nucleotides while, in principle, co-evolution of the nucleotide immediately after the stop codon may be effecting TAG frequency more than that of TGA and TAA frequencies through either contextually-dependent mutational effects [18, 19] or by affecting the efficiency of stop-codon recognition [9–12]. However, the lack of a relationship of TAG frequency and G content does not depend on the nucleotide immediately after the stop codon (Figure 6) indicating that the nucleotide context is not an important factor in explaining the pattern of stop codon usage (Figure 1). Second, the same stop codons can be subject to different selection pressures in different genes due to difference in the levels of expression [20] or other factors. To alleviate the fears that a more general model, one that takes into account the distribution of selection coefficients in a genome, would substantially alter our conclusions we have analyzed the effect of assuming a distribution of selection coefficients. We have shown that if a given frequency of the TAG stop codon is explained by a distribution of selection coefficients then the expected value of the average selection of such a distribution would have to be greater or equal to a selection coefficient that is uniform across all TAG codons and leads to the same TAG frequency in the genome (see Methods). In other words, differences in the strength of selection across different TAG codons make our argument stronger that, on average, the TAG stop codon is unpreferred.

Using a population genetics model modified to describe stop codon and guanine frequencies we demonstrated that stop codon usage can be explained when selection is acting specifically on TAG. The predicted selection regime on TAG, S ₂ , has three properties: it is relatively weak, with N _e s between −0.5 and 1.5, nucleotide content dependent and is positive when G-content <16% and negative when G-content is >16%. The predicted selection strength is weak, on the order of 1/N _e , which is not strong enough to severely restrict the rate of evolution of stop codons. Indeed, such weak selection on individual alleles can be overpowered by genetic drift, which may result in the large variability of stop codon frequencies in our data (Figure 1B). Alternatively, the observed variability of stop codon frequencies relative to the average expectation (compare Figure 1A and 1B) may be due to slight changes in selection pressure on TAG and the rates of A < − > G mutation between different species (Figure 3A).

The G-content dependence of the selection follows from the roughly constant TAG frequency relative to G-content. Yet at this point, there are no known molecular mechanisms that may explain why TAG stop codon has different selective consequences depending on nucleotide content. One possibility, however, is the dependence of translation termination efficiency on the nucleotide context in the vicinity of the TAG stop codon. Bacteria generally code for two release factors (RF), RF1 that recognizes TAA and TAG stop codons and RF2 that recognizes TAA and TGA [21]. Thus, the prediction of the context-dependence hypothesis is that the efficiency of RF1 is GC-context dependent while RF2 functions independent of nucleotide context. Empirical evidence may be necessary to confirm or refute this hypothesis, however, given the relatively weak nature of the selection the differences in translation termination efficiency may be too small to be easily detected in the laboratory. The possibility of the molecular mechanism involving elongation termination factors, however, is left necessarily uncertain by conflicting data from other species. Eukaryota, that have only one release factor for all three stop codons [22, 23], and chloroplast genomes that have retained orthologs of both release factors [24], show a clear increase of TAG frequency with higher GC-content (Additional file 3: Figures S3 and Additional file 4: Figure S4). Clearly, further experimental work is likely necessary to elucidate the molecular mechanisms behind selection on TAG stop codon in bacteria.

Within the framework of our model it is possible to compare the fitness impacts of different stop codons depending on genomic nucleotide content. Regardless of selection on G-content itself (S ₁ ) the difference in fitness between TAG and TGA stop codon is defined by S ₂ (Figure 2). Thus, regardless of the value of S ₁ our data signify that TAG stop codon is always less fit that the TGA stop codon for G-content >16%. Comparing the relative fitness of TAG and TAA, however, involves both S₁ and S₂, with their sum being the difference in relative fitness of these two stop codons (Figure 2). Within the model, G-content depends on relative rates of mutation A < − > G and S ₁ , and we cannot disentangle the contribution of mutation (${\mu }_1$and ${\mu }_2$) versus selection (S ₁ ) so we cannot analytically estimate the value of S ₁ . However, we can define the range of values of these parameters for specific G-content.

To identify the plausible range of S ₁ values in bacterial genomes we investigate the values of ${\mu }_2/{\mu }_1$ and S ₁ for G-content of 16%. We find that in order for the G < − > A twofold sites to maintain 16% G-content either the rate of G > A mutation must be at least five times larger than the rate of A > G mutations or, S ₁ must be positive (Figure 7). Thus, for genomes with ~16% G-content TAG can be the stop codon with the highest fitness only if G > A rate of mutation is five times higher than the rate of A > G mutation. For G-content of 5% TAA would confer higher fitness than TAG if S₁ < 0.5. When S₁ = 0 the ratio of G < − > A mutations (${\mu }_2/{\mu }_1$) must be ~19, while a smaller ratio implies a positive S₁ (Figure 7). When S₁ = 0.5 the ratio ${\mu }_2/{\mu }_1$ must be at least 11, implying that if ${\mu }_2/{\mu }_1$ < 11 then S₁ > 0.5 and TAA is more fit than TAG (Figure 7).

Is there any evidence that G > A can be five or ten times faster than A > G mutations? Mutational biases against GC-content that have been measured were shown to be always less than tenfold in favor of AT-content and less than fivefold for 151 out of a total of 154 species considered in two separate studies (25, 26). Similarly, weak selection acting on GC-content has been postulated by several researchers (25,26). Given this evidence it is unlikely that the observed GC-content can be explained solely by G < − > A mutational biases and, therefore, S ₁ is positive and > > 0.5 for G-content ≈ 5% and >0 for G-content ≈ 16%. Thus, for G-content < 16% the TAG stop codon is expected to be less fit than the TAA stop codon.

The relative fitness of TAG to TAA and TGA stop codons can thus be described as follows. When G-content is >16% TAG has lower fitness than TGA. As long as S ₁  > −S ₂ for G-content <16% then TAG has lower fitness than TAA in bacterial genomes with G-content < 16%. Because S ₁  > −S ₂ is expected to hold for G-content < 16% given the mutation parameters observed in nature [25, 26] it follows that TAG is a striking example of a global suboptimal codon, such that the substitution of TAG into either TAA or TGA for any bacterial species would lead to an increase of fitness. The use of suboptimal synonymous codons in bacteria is a well-documented phenomena, however, the exact codons that are suboptimal differ substantially between different species (see [27] for review). To our knowledge, the observation that one codon with synonymous function to other codons is always worse in such a large group of organisms, bacteria, is the first example of a global sub-optimality of the genetic table. The sub-optimal organization of the genetic table revealed here provides a striking counterexample to the remarkable optimization of the genetic code with respect to error minimization [28–30].

All available complete bacterial genomes were downloaded from NCBI website and 736 of those that utilize the standard genetic code were used for the analysis (See Additional file 5). Plasmid sequences were excluded. All available pairs of closely related genomes from the ATGC database [31], of which there were 11 pairs with 0.03 < K _S  < 0.22 were used to measure A < − > G synonymous transition rates (K _AG ) and rates of stop codon evolution (K _stop ). Orthologues were constructed using two-directional best BLAST [32] hit approach and aligned using MUSCLE [33]. To obtain K _AG we looked at the number of synonymous differences between three pairs of codons: CAA and CAG, AAA and AAG and GAA and GAG. The expected number of substitutions occurred was estimated using Jukes-Cantor model [34]. The same method was applied to estimate the number of substitutions between stop codons with the only difference that the number of synonymous sites for TAA codon is twice as high as the number of synonymous sites for TAG and TGA codons. In order to obtain rates of TAG codon evolution the substitutions have to be polarized and for that the third organism was added to the 11 pairs of the genomes such that the synonymous distance between sister species 0.02 < K _S  < 0.15 and between sister species and outgroup 0.04 < K _S  < 0.62. Substitutions were polarized using simple parsimony approach.

To show that a distribution of selection coefficients for the same stop codon across different genes can only increase the differences between average selection coefficients of stop codons we proved the following conjecture. A given frequency of TAG codon in the genome can be explained by an equal strength of selection acting on all TAG codons in the genome ($S_0$) or a distribution of selection coefficients across different codons with an expected value of the distribution ($\overline{S}$). For any given observed frequency of the TAG codon in the genome $S_0\le \overline{S}$, such that the average strength of selection in a distribution is larger when different codons are under different selection pressures. We consider the case where selection on each TAG stop codon is a discrete random variable which assumes the value S _i with the probability p _i . In this case we use S _i as discrete values of a distribution of selection coefficients on TAG stop codons in different genes in the same genome, while S ₁ and S ₂ were used as fixed values of the selection coefficients for all genes across a single genome. In this case for any selection S _i the expected number of the sites under this selection is N _i  = p _i *N _stop , the frequency of TAG is $f_i^{\mathit{\text{TAG}}}=\frac{1}{1+\frac{\text{exp}\left(S_i\right)}{f_G}}$ and the number of TAG stop codon is $N_i^{\mathit{\text{TAG}}}=f_i^{\mathit{\text{TAG}}}*p_i*N_{\mathit{\text{stop}}}$. The observed frequency of TAG in the genome is $f_{\mathit{\text{TAG}}}=\frac{{\displaystyle \sum _i{N}_i^{\mathit{\text{TAG}}}}}{N_{\mathit{\text{stop}}}}={\displaystyle \sum _i{p}_i*f_i^{\mathit{\text{TAG}}}}$ and the value of selection S ₀ acting on TAG sites is estimated from the formula $f_{\mathit{\text{TAG}}}={\displaystyle \sum _i{p}_i*f_i^{\mathit{\text{TAG}}}}=\frac{1}{1+\frac{\text{exp}{S}_0}{f_G}}$. Taking into account that the second derivative of f, $f^{″}=\frac{\frac{\text{exp}S}{f_G}\left(\frac{\text{exp}S}{f_G}-1\right)}{{\left(1+\frac{\text{exp}S}{f_G}\right)}^3}\ge 0$, if S ≥ 1nf_G, the Jensen’s inequality ${\displaystyle \sum _i{p}_i*f_i^{\mathit{\text{TAG}}}}\ge f_{\mathit{\text{TAG}}}\left(\overline{S}\right)$ holds, or $\frac{1}{1+\frac{\text{exp}\overline{S}}{f_G}}\le \frac{1}{1+\frac{\text{exp}{S}_0}{f_G}}$ and $S_0\le \overline{S}$. The only condition for this inequality to hold is S ≥ 1nf_G, which is a reasonable assumption taking into account the fact that out of 736 genomes analyzed S₀ ≥ 1nf_G for 734 (Additional file 6: Figure S5).

Reviewer 1: Dr Mikhail Gelfand, Institute for Information Transmission Problems, RAS, Bolshoi Karetny per. 19, Moscow 127994, Russia and Faculty of Bioengineering and Bioinformatics, Moscow State University, Vorobievy Gory 1-73, Moscow 119992, Russia.gelfand@iitp.ru

The authors present a model explaining the following observation: while the use of the UGA stop codon depends on G-content, the UAG frequency is almost constant in genomes with highly diverse G-content. While I see no problems with the observations and the model, I have some editorial comments and questions.

The authors state several times – starting with the very first sentence of the abstract – that the usage of stop codons has not been rigorously studied. This is not correct. In the 90’s, several papers considered the usage of stop codons and its dependence on the local context, including tandem stops and tetranucleotides involving stop-codons. I think these papers should be mentioned.

Author response: Indeed, the term “usage” in this context is not very precise. We acknowledge that there have been studies of stop codon usage in the local context, that is to say that some stop codons have a preferred local context, however, in this manuscript we discuss only the evolution and genomic frequencies of the three different stop codons, which to our knowledge has not been rigorously considered previously. We cite some of the relevant literature and use the word “frequency” which we believe is not as ambiguous as “usage” in this context.

How the 11 studied genome pairs were selected?

Author response: We selected all genome triplets with 0.03 < K_S < 0.22 that were available in the ATGC database. We now report this in the Methods section.

Is the G/A content the same in the 3rd codon position in all codon pairs? If not, why this is a good parameter?

Author response: There are three pairs of two-fold degenerated codon families: AAG/A, GAG/A, CAG/A. G-content at the third position of every pair is indeed highly correlated with overall G-content (see the figure below).

Dependency between G content in the third position of two-fold degenerated codon families and overall G content for AAG/A (blue), GAG/A (red), CAG/A (green).

And in any case, what are the reasons to suspect that the selection regime in the amino-acid-encoding codons is the same as in the stops (the former may depend on concentrations of tRNAs and the codon-anticodon interactions; the latter, on interactions with the release factors). What about the A/G choice in the four-fold codon families?

Author response: Indeed, we have created the model based on this assumption because it allowed us to reduce the number of parameters and make the system of equations solvable. However, we can also show that this assumption does not affect our main result that the TAG codon is selectively disadvantageous. Specifically, from system of equations (4) it follows that exp(S ₂ ) = f _TGA /f _TAG . Thus, we can solve for the selective impact of TAG (S ₂ ) solely based on the frequencies of TAG and TGA without making the assumption that the selective regime is the same in stop and amino acid codons. Since S2 is positive for almost the entire range of G content it follows that the TAG codon provides a selective disadvantage relative to the TGA codon. Unfortunately, we cannot estimate S ₂ by comparing the frequencies of TAG and TAA codons because we cannot independently estimate the $\frac{{\mu }_2}{{\mu }_1}\text{exp}\left(S_1\right)$ component of f_TAA from (4). We now present the new estimate of S ₂ in Figure 5 and the main text.

The reasoning in page 6 is not clearly presented, and misprints add to the confusion. How is formula S2 = ln ((fG(1-fTAG))/fTAG) used? Do I understand it correctly that the next formula S2 = ln (3.6fG + 0.4) results from a fit to observations (comparison of genome pairs)? – I think, this should be explained more explicitly.

Author response: Yes, this is what we mean, and we rewrote this section to hopefully make this clearer.

By the way, the two formulas for S2, theoretical and observed ones, yield a dependence between fG and fTAG – does it hold?

Author response: Yes, there is a slight dependence as can be seen from Figure 1.

Finally, reference to equation ( 5 ) in the preceding paragraph should be about equation ( 4 ), and the sentence “S2 has a clear G-content dependence is well approximated…” probably should be “S2 has a clear G-content dependence that is well approximated…” .

Author response: If the referee means this sentence “Thus, selection on G-content,, affects only G-content itself and does not change the form of the relationship between G frequency and stop codon usage as is evident from expressions (4).” then we mean that in the system of equations (4) G-content (f(taa,tga,tag) does not depend on S1. The other typo is corrected.

Polarization of substitutions using parsimony may be dangerous if there is selection towards a specific, preferred nucleotide: in some cases two parallel nonpreferred-to-preferred substitutions may occur, and they will be interpreted as a single preferred-to-nonpreferred substitution, hence skewing the substitution statistics.

Author response: This is true, however, these data has been obtained for a number of species with different GC-content and low sequence divergence. Therefore, we believe that it is unlikely that the use of parsimony have produced a systematic error of substantial effect that jeopardizes our conclusions.

Reviewer 2: Dr. Arcady Mushegian, Stowers Institute for Medical Research, Kansas City, Missouri, United States of America and Department of Microbiology, Kansas University Medical Center, Kansas City, Kansas, United States of America.arm@stowers.org

The manuscript by Povolotskaya et al. puts forward a simple model of nucleotide substitutions in the stop codons in bacteria, and tests it against the genome-wide data. One of the main conclusions is that TAG may be globally suboptimal, with each of the remaining two codons turning out more fit under different values of GC content.

One biological explanation of these data may be in the phenomenon of overlapping ORFs in bacterial operons. TAG is the only codon that does not accommodate a minimal overlap, whereas TAA can give one kind of stop-start codon overlap (TAATG) and TGA even two kinds (ATGA and TGATG). Perhaps if the authors restricted their sample to the termination codons in the last (or only) genes in operons, they would see much less difference between fitness of those two and TAG?

Author response: The idea that the observed pattern of stop codon frequency in bacterial genomes is explained by gene overlap has occurred to us as well. However, we observe the same relationship between G-content and stop codon frequency in overlapping and non-overlapping genes. We now report these data in a new figure that is Additional file 2 Figure S2 in the new version of the manuscript. We have considered only tail-to-tail overlaps due to a much higher certainty of stop codon annotation compared to the uncertainty in the annotation of many start codons.

Reviewer 3: Dr. Shamil Sunyaev, Dr. Shamil Sunyaev, Division of Genetics, Brigham and Women's Hospital, Harvard Medical School, 77 Ave. Louis Pasteur, Boston MA 02115, USA. ssunyaev@rics.bwh.harvard.edu

This manuscript presents an analysis of stop codon usage in bacterial species.

The authors report that TAG codon is un-preferred in most bacterial species and that its frequency does not depend on GC content. They suggest presence of weak selection against TAG codon due to unknown mechanism. One potential mechanism may involve dependency of efficiency of one of the release factors on GC content. I find the results of great interest. I only have two minor technical comments.

1) The analysis is based on Bulmer equations, which hold only if evolution is mutation limited. It would be great to briefly discuss applicability of this model to a wide variety of bacterial species.

Author response: Bulmer’s model assumes that the fate of a new mutation is decided independently of other mutations, that is to say that generally only one mutation is segregating in the population at the same time. This is certainly true if we consider only mutations in stop codons. In most bacterial genomes there are 2–5 thousand protein coding genes making it rather unlikely that more than one stop codon polymorphism is segregating at the same time.

2) Approximation of selection coefficient against TAG codon as a sum of contributions due to selection against GC content (S1) and selection against this specific codon (S2) ignores the S1*S2 term. It is OK if both selective forces are assumed to be small. It would be great if this assumption would be spelled out.

Author response: The referee is absolutely correct, we assume that both of the selective forces are small. We have added an explicit statement to this effect in the text.