|Home | About | Journals | Submit | Contact Us | Français|
Toxin–antitoxin (TA) systems are key regulators of bacterial persistence, a multidrug-tolerant state found in bacterial species that is a major contributing factor to the growing human health crisis of antibiotic resistance. Type II TA systems consist of two proteins, a toxin and an antitoxin; the toxin is neutralized when they form a complex. The ratio of antitoxin to toxin is significantly greater than 1.0 in the susceptible population (non-persister state), but this ratio is expected to become smaller during persistence. Analysis of multiple datasets (RNA-seq, ribosome profiling) and results from translation initiation rate calculators reveal multiple mechanisms that ensure a high antitoxin-to-toxin ratio in the non-persister state. The regulation mechanisms include both translational and transcriptional regulation. We classified E. coli type II TA systems into four distinct classes based on the mechanism of differential protein production between toxin and antitoxin. We find that the most common regulation mechanism is translational regulation. This classification scheme further refines our understanding of one of the fundamental mechanisms underlying bacterial persistence, especially regarding maintenance of the antitoxin-to-toxin ratio.
Persistence is a metabolically inactive state that enables many bacterial species to maintain a subpopulation of cells that can survive harsh changes in the environment [1,2]. From a human health standpoint, persister cells are a growing problem since the metabolic dormancy that characterizes the persister state results in the persister population being multidrug-tolerant and a major contributing factor to ineffective antibiotic treatments. In addition, it has been suggested that persisters indirectly lead to antibiotic resistance; persister cells survive antibiotic treatment and are then able to acquire antibiotic resistant genes from their neighbors [3,4]. Investigations have revealed that a central regulator of bacterial persistence is a network of multiple toxin–antitoxin (TA) systems [5,6,7]. Evidence suggests that TA systems trigger persistence when rare events allow active toxins to accumulate and affect metabolic dormancy by slowing processes such as translation and transcription. A variety of bacterial species have TA systems, which are classified into types based on the mechanism the antitoxin uses to neutralize its cognate toxin [8,9,10,11]. Types I and III have RNA antitoxins that inhibit toxin synthesis or sequester their toxin, while types II, IV and V have protein antitoxins. Type II TA systems are the only type where the antitoxin protein directly binds to the toxin to form a protein–protein complex which sequesters the toxin and effectively neutralizes it . In Escherichia coli alone, there at least 36 known TA systems, most of which are type two [6,12]. In this work, we were able to classify 10 out of 16 type II TA systems (there was insufficient data to classify the remaining six) found in E. coli K12 str. MG1655.
Type II TA systems are operons that encode genes for two proteins, a stable toxin and an unstable antitoxin. Since the antitoxin protein is rapidly degraded by proteases, it must be continually produced to prevent a buildup of free toxin protein and to maintain the susceptible population (non-persister state) (Figure 1) [11,13,14]. Thus, antitoxin is expected to be produced at a sufficiently higher rate than toxin in non-persister cells [8,15]. However, there is disagreement in the literature as to how this ratio is ensured across type II TA systems. It is often cited that transcriptional regulation is responsible for the higher production rate of antitoxin [8,16,17,18,19,20,21]. Research shows that the RnlAB system contains an internal promoter that is independently regulated and allows for the possibility of differential transcriptional regulation , but differential transcriptional regulation has not been confirmed to regulate the production ratio of antitoxin to toxin for most type II TA systems.
A model of TA systems that includes transcriptional regulation, conditional cooperativity, proposes that the antitoxin when complexed with its cognate toxin autoregulates the whole operon and could lead to some control of the antitoxin-to-toxin ratio [14,21], but this theory depends on antitoxin translation rates being higher than toxin . Another model of TA systems also depends on the translation rate of antitoxin being higher than that of toxin, but does not require conditional cooperativity . Many other operons are known to use translational regulation (different translation rates) to maintain differential protein production within an operon, including the ATPase operon [24,25]. As we will support, many type II TA systems use translational regulation, and transcriptional regulation is not the only mechanism that can explain the higher production rate of the antitoxin protein.
This study examines the possibilities of both transcriptional and translational regulation using bioinformatics approaches that combine diverse datasets and analyses. We analyzed type II TA systems in E. coli K12 str. MG1655 using RNA sequencing (RNA-seq) data from multiple studies and growth conditions, data from ribosome profiling analysis (Ribo-Seq), identified promoters and terminators with experimental data and predictions annotated on EcoCyc , and calculated translation initiation rates with prediction algorithms (TIR calculators). By combining these analyses, we classified E. coli type II TA systems into four classes based on the mechanisms of differential protein production. This classification scheme further refines our understanding of how TA systems maintain differential antitoxin-to-toxin expression and one of the fundamental mechanisms underlying bacterial persistence.
A major result of our investigation is that differential transcriptional regulation of antitoxin and toxin expression is unlikely in many type II TA systems. The key pieces of evidence are a less than two-fold difference of antitoxin-to-toxin mRNA for several systems, and a lack of obvious internal promoters or terminators within the operon sequence that could explain excess antitoxin mRNA. We predict that TA systems with less than a two-fold difference in mRNA expression leverage translational regulation to maintain higher production rates of antitoxin. Our results extend the pattern found for many operonic genes that use translational regulation to maintain differential protein production from a single transcript [24,25,27].
We analyzed 10 different type II TA systems in E. coli (Table 1) using data informing their operon organization, mRNA concentration (RNA-seq), protein synthesis rates (Ribo-Seq), and predicted translation initiation rates (TIR’s). Six type II TA systems lacking substantial RNA-seq data were not included. Understanding of the operon organization and mRNA products for each TA system led to a compelling classification scheme that includes at least four major classes. We find that these classes of TA systems all have mechanisms in place to ensure sufficient production of antitoxin protein relative to toxin protein, though the details vary from class to class. In this section, we present the major findings of our study. We first present the details concerning our TA system classification scheme; then we analyze trends in associated quantitative data. We finally discuss our interpretation of these results in terms of the four identified classes of TA systems.
A survey of 10 TA systems (Table 1) using the online database EcoCyc , revealed consistent patterns based on operon organization, which refers to the order of the genes (whether the toxin is upstream at the 5′ end or downstream at the 3′ end of the operon), and whether there is an additional promoter or transcriptional termination mechanism that can produce multiple distinct mRNA products.
Further analysis of our representative set of TA systems led to one major class and three other classes (Figure 2). The major class (Class 1: FicAT, MazEF, MqsAR, PrlF-YhaV, RelBE, YefM-YoeB) found in our study includes TA systems that produce a single transcript, which results in the condition A ≈ T for mRNA (the concentration of antitoxin coding region is approximately equal to the concentration of toxin coding region). Because these systems are apparently not differentially regulated at the transcriptional level, we predict that different translation rates of the two proteins are responsible for higher antitoxin protein production than toxin. These predictions are confirmed later in this article. The three other classes (Class 2: HicAB; Class 3: DinJ-YafQ, YafNO; and Class 4: RnlAB) are inspired by a few systems that have similar gene organization to Class 1, but other features deviate from Class 1 that allow for transcriptional or post-transcriptional regulation of the antitoxin to toxin protein production ratio. Interestingly, we find the ribosome-binding site (RBS) of the downstream gene is embedded in the upstream gene for nine out of 10 TA systems (only HicAB in Class 2 does not).
Many of these classes were identified by examining mRNA expression from RNA-seq data in combination with promoter identification (see Methods). Class 2 is similar to Class 1, but an additional external promoter (located before the coding regions) is near the toxin RBS. Transcription from this promoter results in a truncated mRNA that has a weakened RBS for toxin translation. Thus, antitoxin is produced at a greater rate than toxin . Class 3 likely has a truncated mRNA according to our analysis of RNA-seq data, but due to the gene organization of this class, it probably has either an early transcriptional terminator or post-transcriptional mRNA degradation to truncate the toxin-coding region. Class 4 has an internal promoter that produces excess antitoxin mRNA. Classes 3 and 4 rely on increased abundance of functional antitoxin-coding mRNA relative to toxin. Thus, a higher TIR for antitoxin relative to toxin is perhaps then not as necessary, in contrast to Class 1.
Our classification scheme (Figure 2) was supported using a diverse-data analysis approach. As part of this analysis, we quantitatively estimated the mRNA abundance for antitoxin and toxin coding regions using publicly available RNA-seq data. A total of 13 different whole transcriptome E. coli K12 str. MG1655 RNA-seq datasets were derived from two different studies with a total of six different experimental conditions that include different growth densities and media (Table S1).
RNA-seq analysis shows that the sequencing coverage (number of reads aligned to the gene normalized by dividing the length of the gene) for seven TA systems had less than a two-fold difference in expression at the mRNA level for a variety of conditions (Table 2, Figure 3), as anticipated for Class 1 and Class 2 TA systems. The two TA systems in Class 3 (DinJ-YafQ, YafNO) consistently had more antitoxin mRNA than toxin, which is also as anticipated, though the two-fold difference is modest. Our one example of a Class 4 system (RnlAB) also had less than a two-fold difference between toxin and antitoxin mRNA expression, but the coverage of functional antitoxin mRNA is higher than functional toxin mRNA (Figure S1).
Sequencing coverage of mRNA was suggestive, in that most TA systems had less than a two-fold difference in coverage of antitoxin and toxin. To test our classification of type II TA systems into four classes based on gene organization and RNA-seq data, we analyzed protein synthesis rates from publicly available Ribo-Seq data. Indeed, the Ribo-Seq data supports our classification. Protein synthesis rates were determined using the Li et al. (2014) open source database that gives protein synthesis rates for E. coli based on Ribo-Seq data  (See Methods). Protein synthesis rates were calculated from Ribo-Seq based on the coverage (counts normalized by gene length) of the ribosomally protected mRNA present after the mRNA has been extracted and treated with RNase . From this analysis, seven out of eight of the TA systems with sufficient coverage (HicAB and FicAT had low confidence, less than 128 mapped reads) had a higher protein synthesis rate for the antitoxin than the toxin (Figure 4), as anticipated. We hypothesize that the systems with less than a two-fold difference in toxin and antitoxin mRNA expression likely use translational regulation to maintain the higher antitoxin production rate. The major exception we found was the synthesis ratio for the RnlAB TA system (Class 4), which we hypothesize may be due to third protein that interacts with this system (see comments in Section 2.5).
An independent theoretical investigation of TIRs for the 10 TA systems was done to assess the robustness of our findings based on Ribo-Seq data. We predicted TIRs using three experimentally-conditioned TIR calculators, of which two are based on detailed thermodynamic models (the RBS Calculator [41,42] and the UTR Designer ), and a third is based on a log-linear regression of E. coli gene expression (Barrick Calculator ). The results from the TIR calculators vary greatly due to the inherent differences between the calculation methods, but for six of the seven systems that have less than a two-fold difference in mRNA expression, at least two out of three calculation methods qualitatively agree that for each system the antitoxin TIR is higher than toxin (Table 3 and Table S2). In all cases, the calculators support our classification of the TA systems and provide a method independent of Ribo-Seq.
Our picture from a diverse-data analysis of 10 TA systems has thus led to the following key interpretations of their behavior:
Class 1 is the most abundant class with six TA systems that use translational regulation to maintain a high antitoxin-to-toxin protein synthesis rate, as seen in the Ribo-Seq analysis (Table 2, Figure 4). The TA systems in this class have less than a two-fold difference in mRNA expression (Table 2, Figure 3) and likely use translation regulation to maintain the antitoxin-to-toxin ratio. The RBS site of the downstream gene in these systems overlaps with the upstream gene’s coding region, and they are not dependent upon the order of those genes to maintain the antitoxin-to-toxin production ratio (Figure 2A). Included in this class is RelBE, which previous researchers have shown to use translational regulation to maintain a high antitoxin-to-toxin ratio . The only TA system in this study to clearly use translational regulation that is not in Class 1 is HicAB, which we have placed in Class 2.
The protein synthesis rates of HicAB could not be determined by the Li et al. (2014) data due to low expression of this system, but our RNA-seq analysis shows that the HicAB system has less than a two-fold difference in mRNA expression, like Class 1. Unlike Class 1, HicAB does not have overlapping genes and the antitoxin is located downstream of the toxin (Figure 2B). Our first analysis showed that translational regulation is unlikely because two out of three TIR calculators predict that the antitoxin TIR is lower than the toxin. Interestingly, a recent study showed that HicAB combines transcriptional and translational regulation by using two different promoters at the beginning of the operon to produce two different transcripts with different toxin translation rates . When using the transcription start site of the second promoter, the results of the RBS Calculator change to predict the antitoxin TIR to be higher than the toxin. This entails in two of the three TIR calculators predicting antitoxin TIR is higher than toxin (Table S2), supporting that this system uses translational regulation and supporting the effectiveness of the TIR calculators when using a consensus of two out of three.
Class 3 contains the remaining two systems that had higher antitoxin protein synthesis rates (Figure 4), but these systems also have two-fold higher antitoxin mRNA expression than toxin (Table 2, Figure 3). The organization of these systems has the antitoxin upstream of the toxin with the RBS of the toxin overlapping the end of the antitoxin-coding region (Figure 2C). While the difference in mRNA would indicate that these systems are transcriptionally regulated, the regulation mechanism apparently does not use a promoter; the antitoxin location upstream of the toxin would result in an additional promoter for the antitoxin reading through the toxin as well. We hypothesize that the lower toxin mRNA level is a result of a truncated mRNA due to an unidentified transcriptional terminator or RNA degradation specific to the toxin mRNA sequence. Several type II toxins are endoribonucleases (including those in Class 3) and they could possibly degrade their own message with a bias toward toxin mRNA degradation (Table 1). Regardless of the mechanism, Class 3 TA systems likely use differential mRNA expression to regulate antitoxin-to-toxin ratios either transcriptionally with terminators or post-transcriptionally through RNA degradation.
Class 4 contains RnlAB, the only TA system that does not have a higher antitoxin protein synthesis rate than toxin (Figure 4). The gene organization of this system has the antitoxin downstream of the toxin with the antitoxin RBS overlapping the end of the toxin-coding region (Figure 2D). This organization means that the RnlAB system can use transcriptional regulation of an internal promoter to express antitoxin mRNA higher than toxin, and analysis of the RNA-seq data supports the presence of an internal promoter at approximately 280 nucleotides upstream from the antitoxin (Figure S1). The RNA-seq analysis is supported by a previous study on the regulation of RnlAB, which experimentally determined an internal promoter to be near this location. The same study establishes that the two promoters are independent and transcriptionally regulated, which would allow for upregulation of antitoxin expression . The internal promoter in this system provides a possible mechanism to upregulate antitoxin production, but the protein synthesis rates indicate that toxin production is higher than antitoxin production.
The toxin homologs RnlA and LsoA (not found in the E. coli strain used in this study) are unlike other type II TA systems because the protein structure of these toxins is different from established structures of others. Additionally, a third protein interacts with the RnlAB system, RNase HI. Recent studies have shown that RNase HI acts as a corepressor of the toxin RnlA in the presence of its antitoxin RnlB . However, when antitoxin is not present in the system (in the absence of RnlB), RNase HI stimulates RnlA activity . This system also acts as an anti-phage response, because during infection RnlA degrades the infecting phage’s mRNA and RnlB concentration decreases allowing RNase HI to activate RnlA and strengthen its activity [45,46]. The third component of the RnlAB system changes the dynamics and is likely the reason that the ratio of antitoxin protein synthesis rate to the toxin is lower than expected.
Currently, many studies on regulation of TA systems focus on conditional cooperativity, which depends on transcriptional regulation by TA autoregulation complexes. Numerous studies on TA autoregulation complexes focus on RelBE, a Class 1 system. The RelBE system is known to produce the antitoxin (RelB) approximately 10-fold faster than the toxin (RelE) , which is further supported by the protein synthesis rate data, which predicts RelB to be produced 7.5 times greater than RelE . Studies on RelBE show that DNA binding of the RelB2E complex is stronger than the RelB2 alone [47,48], and studies suggest that the stronger autoregulation activity of the complex is important to the control of the antitoxin-to-toxin ratio [20,47]. However, in these models, the translation rate of antitoxin being higher than that of the toxin is critical to prevent a constant overabundance of free toxin, and therefore the system is still dependent on translational regulation [16,20]. In contrast to the TA autoregulation complex of the RelBE system, the DinJ-YafQ complex does not have increased DNA binding affinity when compared to DNA binding affinity of the antitoxin (DinJ) alone . Interestingly, the DinJ-YafQ system was placed into Class 3, which our results suggest relies on increased abundance of antitoxin mRNA rather than translation rates. Further studies are needed to fully understand the variety of regulation mechanisms that play a role in the maintenance and control of the antitoxin-to-toxin ratio.
Both translational and transcriptional regulation play important roles in the maintenance of the production ratio of antitoxin-to-toxin in type II TA systems in the non-persister state. Analysis of RNA-seq data reveals that most TA systems have less than a two-fold difference between antitoxin mRNA and toxin under the conditions studied in these datasets (see Table S2). The antitoxin-to-toxin production ratio in Class 1 TA systems, which represent the most abundant class in our study, is not differentially regulated transcriptionally or post-transcriptionally, but it is regulated at the translational level. Class 2 uses a combination of both transcriptional and translational regulation to ensure a higher rate of antitoxin production. Class 3 is likely not regulated translationally but by some mechanism that results in a truncated mRNA and overall greater expression of antitoxin mRNA than toxin. Class 4 contains RnlAB, which is regulated by internal promoters, but the interaction of RNase HI with RnlAB lessens its dependency on a higher production rate of antitoxin.
Our classification of TA systems emphasizes their diversity with respect to gene organization and regulation mechanisms, while other classification systems are based on protein structures and functions [10,13]. The diversity of regulation mechanisms explains the disagreement in the literature as to whether TA systems use transcriptional regulation or some other form of regulation to maintain an antitoxin production rate higher than the cognate toxin. Various studies have focused on specific TA systems and their conclusions have been applied to type II TA systems as a group, but our results emphasize the diversity found within these systems. Further applications of the methods used in this study would be to classify and expand the knowledge of TA systems in other organisms and other differentially expressed operons in bacteria. Using bioinformatics methods alone, the methods in this study can be applied to classify other type II TA systems. Classification of different type II TA systems into these classes should allow future researchers to predict regulation (transcriptional, post-transcriptional, or translational) without the expense and time of RNA-seq and Ribo-Seq experiments.
The annotated E. coli K12 str. MG1655 reference sequence NC_000913.3 was used for all gene sequences in this study.
Selected RNA-seq datasets (Table S1) from GEO NCBI (accession numbers GSE48829 and GSE74809) [50,51] were aligned to the E. coli K12 str. MG1655 reference sequence NC_000913.3 using Geneious v. 10.0.9 . Gene expression was analyzed in Geneious using the calculate expression tool to generate a list of reads mapped to each gene. Coverage for any specific antitoxin or toxin was calculated in Reads Per Kilobase Million (RPKM). Defining as the raw number of reads mapped to the coding sequence (CDS) (with partial reads, reads mapping to more than one region of the genome, counted as 0.5 reads), as the total number of mapped reads for a particular RNA-seq dataset, and as the number of bases in the CDS, then RPKM = . Results were exported to a CSV file, which was then processed by custom Python scripts using the NumPy library .
The ratio of antitoxin to toxin was calculated for each TA system for all 13 datasets, and the median of those 13 ratios was used to determine the A/T ratio. The standard deviations were calculated using the 13 ratios from the datasets (Table 1).
Error estimates for gene coverage were computed as follows. Each RNA-seq dataset contained biological replicates, and we analyzed error by comparing gene expression between replicates (see Figure S2 for an illustration of replicates). Dataset GSE48829 contained triplicate data for one experimental condition, while dataset GSE74809 contained duplicate data for five experimental conditions. For each dataset separately, we ran an error analysis of the log-error (standard error of logarithmic quantities) in two major directions, ratio and magnitude. We choose log-error since this is the most natural representation of the error for features plotted in log-space, as with Figure 3. If and are gene coverages for antitoxin and toxin, respectively, then the error of the log-ratio is the standard error for the quantity . Since when and are close in value, this error also approximates the error of the ratio (the constant -1 does not contribute to the standard error). This error is plotted as error bars along the ratio direction in Figure 3 after scaling by the factor to plot in logarithm base 10 space. A complementary error estimate is that of the log-magnitude, which we define as the standard error of the quantity , and which we plot as error bars along the magnitude direction in Figure 3 after scaling by . Notice that the log-ratio error is zero for quantities with zero uncertainty in the ratio , even though the log-magnitude error may be significant.
We used protein synthesis rates for E. coli calculated by Li et al. 2014  from their Ribo-Seq data according to their open source database. Genes with less than 128 mapped reads were given values of low confidence in this database.
TIRs were predicted using three different translation rate initiation (TIR) calculators with the reference sequence, NC_000913.3. The reverse engineering feature of the Ribosome Binding Site Calculator v2.0 was used [41,42]. Any difference in sequence length using the RBS Calculator v2.0 changed the translation initiation rates because this software accounts for secondary structure. Therefore, the input sequence we used included the entire 5′ UTR region as determined by the transcription start site annotated in the EcoCyc database  to the end of the TA system. The reverse engineering feature of the UTR Designer was used with the required input of −25 to +35 from the start codon for each gene . The third method analyzed the sequence −11 to +1 from the start codon using an equation based on empirical data developed by Barrick et al. in their study of E. coli ribosome binding sites (Barrick Calculator) . We used the agreement of two out of three calculators to determine whether antitoxin TIR is predicted to be higher than toxin.
This work was supported by funds from the National Science Foundation Division of Molecular and Cellular Biosciences, MCB-1330180.
The following are available online at www.mdpi.com/2072-6651/9/7/211/s1, Figure S1: RNA-seq reads from experiment SRX1424838 (GSE74809) mapped to mazEF, and rnlAB; Figure S2: Illustration of RNA-seq coverage for all replicates and conditions used; Figure S3: Global error analysis of RNA-seq coverage; Table S1: Selected RNA-seq experiment numbers and conditions for GSE48829 (grey background) and GSE74809; Table S2: Ratios of the calculated antitoxin to toxin translation initiation rates.
H.S.D. performed RNA-seq data and translation initiation calculator analysis to classify the TA systems and wrote the manuscript. N.C.B. initiated and directed the project. R.V.J. was essential to understanding the methods of RNA-seq data analysis. W.H.M. directed the project and performed statistical analyses for the RNA-seq data. All the authors took an active part in writing and editing the manuscript.
The authors declare no conflict of interest.