Functionality of tRNAs encoded in a mobile genetic element from an acidophilic bacterium

ABSTRACT

The genome of the acidophilic, bioleaching bacterium Acidithiobacillus ferrooxidans, strain ATCC 23270, contains 95 predicted tRNA genes. Thirty-six of these genes (all 20 species) are clustered within an actively excising integrative-conjugative element (ICEAfe1). We speculated that these tRNA genes might have a role in adapting the bacterial tRNA pool to the codon usage of ICEAfe1 genes. To answer this question, we performed theoretical calculations of the global tRNA adaptation index to the entire A. ferrooxidans genome with and without the ICEAfe1 encoded tRNA genes. Based on these calculations, we observed that tRNAs encoded in ICEAfe1 negatively contribute to adapt the tRNA pool to the codon use in A. ferrooxidans. Although some of the tRNAs encoded in ICEAfe1 are functional in aminoacylation or protein synthesis, we found that they are expressed at low levels. These findings, along with the identification of a tRNA-like RNA encoded in the same cluster, led us to speculate that tRNA genes encoded in the mobile genetic element ICEAfe1 might have acquired mutations that would result in either inactivation or the acquisition of new functions.

KEYWORDS:

• Transfer RNA
• horizontal gene transfer
• mobile genetic elements
• acidophilic bacterium
• translation
• phylogeny

Introduction

Fidelity of gene expression is crucial for living cells. Organisms spend a major portion of their resources to carry out this process accurately, responsively and efficiently. Transfer RNAs (tRNAs) are molecules that perform the pivotal role of adaptor, translating the genetic information coded in the language of nucleotide triplets into the amino acid alphabet of protein products. The role of tRNAs in protein biosynthesis is well known. tRNAs interact with aminoacyl-tRNA synthetases to induce faithful aminoacylation with the cognate amino acid and deliver the aminoacyl-tRNA to the ribosome by interacting with initiation and elongation factors.^1,2 Based on the existence of 61 coding triplets and the potential wobble between the codon in the mRNA and the anticodon in tRNAs, it is known that at least 31 tRNAs are required for protein biosynthesis to take place. The number of tRNA genes is near the lower limit in a few microorganisms, but in most known organisms the number of encoded tRNAs far exceeds the minimum. In bacteria and archaea numbers range from the minimum to 167 tRNA genes.³ In eukaryotes, encoded tRNAs range from 36 in Plasmodium to more than 500 in mammals and plants, and even thousands in some fishes.⁴ Probably some member of the families of tRNAs are overrepresented and the origin and role of this overabundance is a question to be solved. The apparent redundancy of the genetic code, in which most amino acids can be translated by more than one (synonymous) codon, might account in part for such an overabundance of tRNAs. The existence of iso-species of tRNAs with differential codon-anticodon or ribosomal affinities might influence the translation efficiency and accuracy of synonymous codons.⁵ Thus, evolutionary processes can finely tune translation efficiency and accuracy through variations in the tRNA pool, while maintaining the amino acid sequence of proteins.

Horizontal gene transfer (HGT) is a major force in bacterial evolution,⁶ involving the lateral transfer of genetic information from one organism to another, even across diverse lineages.^7,8 Several factors influence the success of the transfer and fixation of a gene or group of genes in a recipient organism: the toxicity of the encoded proteins, the capacity of genes to integrate into the chromosome, protein interactions of the gene product(s) and the potential advantages that the transferred gene(s) confer to the recipient organism. Recent findings from genomic analyses suggest that an additional and in fact crucial determinant of a gene's chances of being fixed through HGT is the adaptation of the codon usage to the tRNA pool of the recipient organism. It seems that the more similar the tRNA pools between microorganisms, the greater the tendency to share genes through HGT⁹. In this context, the presence of tRNA genes in viral genomes is an intriguing issue. It has been proposed that viral tRNA genes might contribute to adapting the bacterial tRNA pool to the viral genome, facilitating the expression of viral genes in the context of the bacterial tRNA complement. This strategy is particularly relevant when there are marked differences in the GC content between viral and bacterial genomes, as has been proposed for cyanophages.¹⁰

Acidithiobacillus ferrooxidans is a Gram-negative, acidophilic γ-proteobacterium belonging to the consortium of microorganisms that participate in the bioleaching of minerals. Because of its capacity to oxidize either ferrous ions or reduced sulfur compounds, this microorganism provides a model for studying the biochemistry and physiology of bioleaching.¹¹ The genome sequences from 2 A. ferrooxidans strains (ATCC 23270 [GenBank: CP001219] and ATCC 53993 [GenBank: CP001132]) are available in public databases. Comparison of the 2 sequences reveals the presence of DNA segments in each genome not present in the other. In ATCC 23270, a 300-kbp actively excising integrative-conjugative element (ICEAfe1) was identified. Since the genes required to conjugational transfer are encoded and transcribed, we believed that this genetic element might have the capacity to transfer from A. ferrooxidans to other microorganisms.¹² A mobile genetic element (MGE) inserted in a similar genomic context was found in Acidithiobacillus caldus.¹³ Bioinformatic analysis reveals a cluster of 36 tRNA genes in the central region of ICEAfe1.¹⁴ A similar array of tRNA genes was found in Acidithiobacillus ferrivorans, a close relative to A. ferrooxidans.¹⁵ We have previously demonstrated that a tRNA^Glu encoded in ICEAfe1 from A. ferrooxidans is aminoacylated by the glutamyl-tRNA synthetase.¹⁶ These observations raised the question as to what role these tRNA genes, clustered in a MGE, might play in A. ferrooxidans. To shed light on this question, in this work we predicted the effect of these genes on the global translation efficiency of the microorganism. Our data indicate that, from a theoretical point of view, the tRNA gene cluster might negatively contribute to translation efficiency in A. ferrooxidans. Although we observed that some tRNAs encoded in ICEAfe1 may participate in protein synthesis, their predicted negative effect on translation might be compensated for by reduced expression of the tRNA genes.

Results

Effect of tRNAs encoded in ICEAfe1 on the global adaptation of mRNAs to the A. ferrooxidans tRNA pool

Bioinformatic analysis of the 2 sequenced genomes from A. ferrooxidans available in public databases revealed the existence of 2 sets of tRNA genes in strain ATCC 23270. One set contains 59 genes shared with strain ATCC 53993 as well as a tRNA^Ser gene predicted to be interrupted by a class I intron. The other set, containing 36 genes that decode all 20 amino acids (although not all codons can be translated by this set of tRNAs), is confined to a 10-kbp DNA segment encoded within ICEAfe1 (S1 Fig. 1).¹⁴ A compilation of the tRNA families encoded in the chromosome and the ICEAfe1 is shown in S2 Table 1. The ARAGORN bioinformatic platform predicted that 13 out of 95 were intron-containing tRNA genes. Based on RT-PCR analyses and sequencing of the resulting products (data not shown), we concluded that tRNA^Ser _{UGA 75} is the only intron-containing tRNA gene that is transcribed and processed as a functional tRNA (see below). A pseudo-tRNA (tRNA-OTH) and a tmRNA were also predicted from this analysis. As ICEAfe1 is a MGE not present in strain ATCC 53993 or other tested strains,¹² we believe that the extra set of tRNA genes was either acquired by duplication events from the chromosomal genes or by HGT. These hypotheses were approached by 2 procedures. In the first case, it was expected that tRNA genes located in ICEAfe1 should be more closely related to the rest of tRNA genes located in the chromosome from Acidithiobacillus ferroxidans ATCC 23270 than with tRNA genes present in other bacteria. To test this hypothesis we search for homologous of the ICEAfe1 tRNAs in the GtRNAdb³ which contains information of 247,731 tRNAs predicted in 4,039 bacterial and 184 archaeal genomes. This analysis showed that tRNA genes present in ICEAfe1 have their closest homologous in tRNA genes present in bacteria from different genera (Acidithiobacillus, Arthrobacter, Beijerinckia, Brachybacterium, Dehalococcoides, Deinococcus, Desulfarculus, Desulfobulbus, Desulfotomaculum, Dictyoglomus, Escherichia, Fimbriimonas, Geobacter, Kineococcus, Kitasatospora, Leuconostoc, Mesorhizobium, Myxococcus, Nitrospira, Oscillatoria, Pristionchus, Rhodocyclaceae, Rickettsia, Sulfurihydrogenibium, Sulfuritalea, Synechococcus, Terriglobus, Weissella, Xanthobacter, Xylanimonas, and Yersinia). Best scores were achieved by the genes encoding tRNA^Val ₃₉, tRNA^His ₄₇ and tRNA^Ala ₅₈ that matched with tRNA genes present in Sulfuritalea hydrogenivorans, Yersinia ruckeri and Terriglobus saanensis, respectively (S3 Table 2). To perform a phylogenetic reconstruction, the tRNA genes from the GtRNAdb, tRNA genes from draft genomes from members of Acidithiobacillus genus together with sequences from ICEAfe1 and from other microorganisms were used.

Figure 1. Cladogram of genes encoding for bacterial tRNAs from acidophilic and other microorganisms: Phylogenetic reconstruction of bacterial genes encoding for tRNAs that charge lysine, arginine, serine and alanine. The colored symbols represent the tRNAs present in ICEAfe1 from A. ferrooxidans ATCC 23270 (red circle) or chromosomal tRNAs from, Acidithiobacillus sp GGI221, A. ferrooxidans strains ATCC 53993, 23270, YQH1, Hel18, BY0502; A. ferrivorans strains SS3, CF27, DLC5, PQ33, YL15; A. thiooxidans strains 19377, A01, Licanantay, DXSW, BY02, ZBY, A02, DMC, GD13, JYC17; A. caldus SM-1, MTH04, DX, ZBY, ZJ, S1 (green rhombus) and chromosomal tRNAs from Leptospirillum ferriphilum ML-04, L. ferrooxidans C2–3 (blue rhombus). tRNAs from other bacteria were obtained from tRNA database. To improve the visualization, clades including tRNAs not related with acidithiobacillus genera were collapsed and are shown in gray lines. Clades that contain mainly tRNAs from Acidithiobacillus genus were collapsed, named as groups and signaled with green lines when they contain tRNA(s) that belong to ICEAfe1.

Phylogenetic reconstructions showed that tRNAs from Acidithiobacillus cluster together forming groups that contain either tRNAs from different Acidithiobacillus, or from a particular specie (Fig. 1, S4 Table 3). Four tRNA gene families (tRNA^Ser, tRNA^Lys, tRNA^Ala and tRNA^Arg) encoded in the ICEAfe1 from A. ferrooxidans were subjected to phylogenetic analysis. Each one showed a different relationship with the tRNAs encoded in other acidophiles or in other bacteria. tRNA^Ala clustered together with tRNAs present in bacteria that belong to the genus Geobacter, Desulfovibrio, Bordetella and Nitrosomona instead of tRNA genes encoded by other acidophilic bacteria.

tRNAs^Lys clustered in 2 different groups. tRNA^Lys ₃₁ clustered with tRNA^Lys encoded in A. thiooxidans (Group IIa, strains: JYC17, ZBY, GD13, DMC, A02, BY02, DXSW, A01), A. ferroxidans DLC5 and A. ferrivorans CF27. The closest to tRNA^Lys ₃₁ is encoded in A. ferrivorans, sharing 99% identity. tRNA^Lys ₄₉ clustered together with tRNAs encoded in A. thiooxidans (Group IV, strains: JYC17, ZBY, DMC, A02, BY02, A01, and Licanantay). However, despite clustering together, they only share 82% identity with 50% coverage in average.

tRNA^Arg genes seem be more related to other acidophiles. tRNA^Arg ₃₆ gene is close to the group XIII that includes tRNA genes from A. thiooxidans (strains: JYC15, ZBY, GD13, DMC, A02, BY02, DxSW and A01), Acidithiobacillus sp GGI221, A. ferrivorans CF27 and A. ferrooxidans DLC5, sharing on average 86% identity. tRNA^Arg ₂₆ and tRNA^Arg ₂₇ both clustered with tRNA genes encoded in A. ferrooxidans DLC5, however tRNA^Arg ₂₆ shared 93% identity with tRNA genes from strain DLC5 while tRNA^Arg ₂₇ shared only 84% identity.

tRNA^Ser genes also showed a complex pattern, tRNA^Ser ₄₈ gene clustered with tRNA genes present in Agrobacterium tumefaciens, Deinococcus radiodurans and Caulobacter crescentus, although with low bootstrap. tRNA^Ser ₅₅ clustered together with tRNA genes from Leptospira interrogans, and Halobacterium, 2 genus not related with the environment of A. ferroxidans. In contrast, tRNA^Ser ₅₂ was part of a cluster (Group IV) that includes tRNAs from A. ferrivorans CF27, Acidithiobacillus GGI221, A. thiooxidans (Strains: A01, Licanantay, BY02, DMC, JYC17, and ZBY). tRNA^Ser ₅₂ is closer to tRNA genes from A. thiooxidans than tRNA genes of A. ferrivorans CF27, sharing 89% and 80% identity, respectively. tRNA^Ser ₅₃ is close to the group II (A. thiooxidans strains: A01, DXSW, BY01, A02, DMC, GD13, JYC17 and ZBY) and near to tRNAs present in A. ferrivorans CF27. However tRNA^Ser ₅₃ gene has only 79% identity with tRNA genes located in A. thiooxidans and 83% with tRNA gene from A. ferrivorans CF27 but with 51% coverage only.

These results led us to think that tRNA genes encoded in ICEAfe1 are not paralogous of tRNA genes located in the rest of the chromosome from A. ferroxidans ATCC 23270, but could be orthologous of some tRNA genes present in other Acidithiobacillus such A. thiooxidans, or A. ferrivorans, as well as other microorganisms and were acquired by several HGT events.

We then asked whether the additional set of tRNA genes provides properties that contribute to this organism's fitness to the environment. The translation rate of a particular mRNA depends in part on its concentration and the abundance of tRNAs that translate it.¹⁷ A means to theoretically evaluate the adaptation of mRNAs to the tRNA pool of the cell is to measure the tRNA Adaptation Index (tAI).¹⁸ This index is defined as the geometric mean of the availability of tRNAs that decode each codon in certain mRNA. Highly-translated mRNAs generally have codons that match anticodons of highly-represented tRNAs. Note that the availability of tRNAs is defined as the copy number of genes, assuming that all tRNA genes are uniformly expressed.¹⁸ To predict whether the tRNAs encoded in ICEAfe1 might influence the translation efficiency of A. ferrooxidans ATCC 23270, we estimated the global tAI of mRNAs in this strain (Fig. 2). The presence of the tRNA genes clustered in ICEAfe1 negatively contributes to the measured tAI for mRNAs from A. ferrooxidans ATCC 23270. Excluding these genes from the analysis (virtual mutant), the tAI shifted significantly toward higher values (Fig. 2A). The log-ratios of tAI values for the wild-type and the virtual mutant revealed a tendency toward values below 0, indicating that tRNA genes from ICEAfe1 have a negative influence on tAI (Fig. 2B). A smaller effect on tAI was observed (ratio closer to 0) when the analysis was performed exclusively for the ORFs encoded in ICEAfe1. This result might imply that the mRNAs encoded in this genetic element are better adapted to the tRNAs present in ICEAfe1. A comparison of the codon usage between the chromosomal DNA and ICEAfe1 revealed that they are similar one to each other, probably as the result of co-evolution of the MGE and the host genome (S5 Table 4). A similar cluster of tRNA genes was found in A. ferrivorans strain CF27, a close relative to A. ferrooxidans. Although not all tRNA genes present in ICEAfe1 are encoded in this cluster, a similar synteny was found indicating probably a common origin.¹⁵ A global analysis in the genus Acidithiobacilli revealed that codon usage of these 2 species is very similar (S6 Fig. 2), thus we believe that the effect on tAI of the cluster of tRNAs should also be similar in the 2 species.

Figure 2. tRNA adaptation index (tAI) of genes from Acidithiobacillus ferrooxidans strain ATCC 23270 (A)and B. Gene frequency vs. tAI value for wild-type (wt) or virtual mutant (vm, tRNA genes from ICEAfe1 were excluded from calculations). (C) and D,. gene frequency vs. log-ratio of the tAI values for the wild-type and virtual mutant, as in (A) and B, except that the calculations were performed based on the actual tRNA levels obtained from RNA-seq data (Experiment 1). Total genes represent genes from chromosome and ICEAfe1. ICE genes are genes from ICEAfe1 only.

Measurements of tAI values reported in the previous section considered the expression of tRNAs as a function of the number of tRNA genes in the genome (as it is defined).¹⁸ Since the additional set of tRNA genes in ICEAfe1 has a negative influence on the theoretical tAI values, we asked whether the actual tRNA levels might counteract the effect on tAIs. To assess this possible effect, we first used RNA-seq to measure the level of tRNAs from A. ferrooxidans grown in standard laboratory conditions using ferrous ions as the electron source. Surprisingly, tRNAs encoded in ICEAfe1 showed low levels (reads) compared with the rest of the tRNAs encoded in the bacterial chromosome under the growth conditions analyzed (Fig. 3). These results were qualitatively independent of the method used for RNA preparation (S2 Table 1). Although the RNA-seq findings might be attributable to effects other than the real level of tRNAs (such as different chemical modifications of tRNAs), these results were consistent with the data obtained by Northern blot analysis (see below, Fig. 4 and S7 Fig. 3). Thus, we recalculated tAI values considering the level of each tRNA according to RNA-seq data. As the expression of tRNAs encoded in ICEAfe1 is generally low, we expected that the recalculated tAI values would be greater than the theoretical values calculated based on an estimation of tRNA levels as a function of the number of tRNA genes. As shown in Fig. 2C, the tAI values calculated using RNA-seq data (Experiment 1 Materials and Methods) were lower than those obtained based on tRNA gene copy number, probably as a consequence of the marked difference between the expression levels (around 150 times using RNA-seq data, and 10 times using gene copy number). However, the log of tAI ratios between the wild-type and the virtual mutant substantially shifted to positive values, implying that the low expression level of tRNAs encoded in ICEAfe1 may not hamper translation efficiency. Similar results were obtained when only the ORFs encoded in ICEAfe1 were considered in the analysis, supporting the global data obtained (Fig. 2D). An interpretation of these results is that the theoretical negative effect of ICEAfe1 tRNAs, deduced from tAI values, can be compensated for by a reduction in the actual level of these tRNAs. Northern blot analysis of tRNAs obtained from A ferrooxidans ATCC 23270 cultured in ferrous ions or sulfur as the electron source as well as in the presence of mitomycin C, a DNA damage inducer that increases the excision and the level of several mRNAs encoded in ICEAfe1¹² revealed no changes in the level of some tRNAs encoded in the MGE (data not shown).

Figure 3. Expression level of tRNAs from A. ferrooxidans. Number of reads for each tRNA from RNA-seq, ordered as they are organized in the genome. Data of reads for each tRNA are in S2 Table 1. Horizontal bar points out the tRNAs encoded within ICEAfe1.

Figure 4. Functional analysis of tRNA^Ser _UGA. a (Upper panel). Schematic representation of the genes (arrows) serTi, encoding the tRNA^Ser _{UGA 75}-containing intron (light gray), and serT, encoding the canonical tRNA^Ser _{UGA 55} (ICEAfe1, white bar; chromosome, black line. Not at scale), from A. ferrooxidans strain ATCC 23270. b (Lower panel) Northern blot analysis of periodate-treated total tRNA from A. ferrooxidans before or after (OH-/OH⁺) incubation with Tris-acetate, pH 9, followed by electrophoresis under denaturing conditions, revealed with specific probes for tRNA^Ser _{UG A 75}, tRNA ^Ser _{UGA 55} or tRNA^Gln _UUG (chromosomal tRNA used as control). A tRNA preparation from strain ATCC 53993 that lacks ICEAfe1 and therefore contains only the chromosomally-encoded version of tRNA^Ser _UGA was also used as a control.

Functionality of tRNAs encoded in ICEAfe1

Based on the calculated tAI values, it appears that the theoretical negative effect of tRNAs encoded in ICEAfe1 on the global expression of A. ferrooxidans mRNAs is counteracted by a low expression level of these tRNAs, as observed with RNA-seq. This effect might be the result of defects on the transcription promoters of these genes. Bioinformatic analysis using neural network promoter prediction¹⁹ revealed that only a few σ⁷⁰ transcription promoters can be predicted in the DNA region encoding the tRNA genes (data not shown). Whether the tRNA genes encoded in the ICEAfe1 are transcribed as a single unit is to be demonstrated (see RNA-seq data of experiment 2, Table S2 Table 1).

Alternatively, an accumulation of mutations that make the tRNAs non-functional might also contribute to reducing their predicted negative effect on translation. Previous data from our laboratory revealed that the tRNA^Glu _{UUC 34} encoded within ICEAfe1 is aminoacylated in vitro by GluRS1, one of the 2 GluRSs present in A. ferrooxidans, and that the Glu-tRNA^Glu formed is bound to EF-Tu, suggesting that it is functional in protein synthesis.¹⁶ We extended the functional analysis to in vivo aminoacylation of some tRNAs encoded in ICEAfe1. Based on RNA-seq data, 5 tRNAs were selected for further analysis: tRNA^Ser ₅₅, tRNA^Trp _61, tRNA^Val ₃₉, tRNA^Ala ₅₈ and the putative tRNA-OTH₆₀ (see below). These tRNAs gave 73, 111, 138, 398 and 885 reads, respectively, on RNA-seq (Fig. 3 and S2 Table 1, experiment 1). tRNA-OTH gave the highest number of reads among ICEAfe1-encoded tRNAs.

tRNA^Ser _(UGA) is relevant to the analysis since it is the only tRNA that can decode the UCA codon. Bioinformatic analysis revealed the presence of 2 genes encoding different versions of this tRNA in A. ferrooxidans ATCC 23270. One gene (serT) encodes a canonical tRNA (tRNA^Ser ₅₅) within ICEAfe1. The second copy (serTi), encoding tRNA^Ser ₇₅, is localized in the chromosome and encodes a predicted tRNA precursor containing a putative class I intron disrupting the anticodon (Fig. 4A). Both genes are transcribed, although at different levels, as shown by RNA-seq (73 vs. 727 reads on RNA-seq respectively; Fig. 3, S2 Table 1) and Northern blot analysis (Fig. 4B). Both mature products (sizes are as expected for the predicted mature products) are aminoacylated in vivo in A. ferrooxidans (Fig. 4B). Complementation of a thermosensitive mutation in serT from Salmonella enterica serovar Typhimurium strain TH 3238 revealed that both versions of A. ferrooxidans tRNA^Ser _UGA are functional in protein synthesis (data not shown). Taken together, these data led us to conclude that although it is expressed at a lower level compared with the chromosomal version, the tRNA^Ser ₅₅ encoded in ICEAfe1 is probably functional in protein synthesis in A. ferrooxidans. Global analysis of the distribution of the UCA codon in A. ferrooxidans genome revealed that this codon is represented at a higher level in genes encoding DNA-related functions as well as hypothetical and unknown functions.¹⁴ These types of genes represent around 81% of the genes encoded within ICEAfe1.¹² Thus, tRNA^Ser ₅₅ expression might be relevant for the expression of the gene products encoded in ICEAfe1. Additionally, tRNA^Val ₃₉, tRNA^Ala ₅₈ and tRNA^Trp ₆₁, distributed in different locations within the cluster of tRNAs encoded in ICEAfe1, were also functional in in vivo aminoacylation in A. ferrooxidans under the selected growth conditions (S7 Fig. 3).

Based on the tRNAscan bioinformatic platform to predict tRNA genes,²⁰ an RNA from the tRNA gene cluster in ICEAfe1 (encoded between coordinates 1035589−1035660) was predicted as a tRNA-like RNA (tRNA-OTH, S2 Table 1) The predicted secondary structure of this RNA is shown in Fig. 5A. An analysis based on the ARAGORN platform for the prediction of intron-containing tRNAs²¹ revealed that an approximately 150-nucleotide-long intron-containing tRNA-like RNA is encoded in this genomic region. The predicted secondary structure of this RNA is shown in Fig. 5B. Based on circular RT-PCR, as well as RNA-seq, we observed that an approximately 140-nucleotide-long RNA including the predicted intron is the predominant RNA species transcribed from this locus in ICEAfe1. Secondary structure prediction based on the MFold program²² revealed that this RNA might be formed by the fusion of 2 putative tRNA genes to form a new RNA with an as-yet-unknown function (Fig. 5C). As expected for a tRNA-like RNA, no aminoacylation was observed for the RNA product (data not shown). BLAST searches against the entire GenBank nucleotide database revealed, among others, striking similarities of this RNA to plasmid-borne sequences from E. coli and K. pneumoniae and genomic island-borne sequences from Mesorhizobium loti. As with ICEAfe1 in A. ferrooxidans, these RNAs form parts of clusters along with other tRNA genes in these organisms. Whether this tRNA-like RNA plays a role in ICEAfe1 requires to be investigated.

Figure 5. Predicted structure of the RNA encoded by tRNA-OTH. Schematic representation of the predicted RNA based on tRNAscan (A), ARAGORN (B) and MFold (C). Note that RNAs predicted in A and B are contained in C.

Taken together, the data reveal that most tRNAs encoded within ICEAfe1 are expressed at low levels probably as result of genetic defects to compensate potential negative contributions to translation in A. ferrooxidans.

Discussion

The major function of tRNAs is transporting amino acids to the ribosome to decode the mRNA for accurate protein translation after aminoacylation by aminoacyl-tRNA synthetases with the cognate amino acid. The remarkably high number of tRNA genes in many organisms, particularly in eukaryotes, well in excess of the requirements of protein translation, led us to explore the role of this apparent redundancy of tRNA in cells. An interesting discussion is taking place in the literature about how the adaptation of mRNA codons to the tRNA complement affects translation efficiency.^{5, 9, 10} It is believed that genes that are highly expressed tend to have a higher adaptation of their codons to the tRNA pool than genes with lower expression levels. Synonymous codons (codons that decode the same amino acid) likely play different roles in translation. Synonymous codons may be unevenly distributed among the various mRNAs, tRNAs might bind to the codons with differential affinity, and the tRNAs might have different concentrations within the cell. Moreover, the tRNA content might change under different conditions, altering the translation speed of certain mRNAs and possibly the folding of the encoded proteins.⁵ The presence of tRNA genes in many viral genomes has long been an intriguing issue for researchers. It is believed that these genes contribute to adapting the bacterial tRNA pool to the viral genome. This strategy allows for efficient translation of viral genes as an alternative strategy to the massive adaptation of the viral DNA sequence to the bacterial genome. Both strategies have been observed in cyanophages¹⁰ allowing the virus to persist in bacterial cells. The discovery of giant viruses carrying several genes involved in translation of genetic information²³ supports the idea that viruses might contain genetic information relevant for expression of their own genome.

An analysis of tRNA clusters in members of acidithiobacilli revealed that an array similar to that in ICEAfe1 is present in A. ferrivorans, another acidophilic bacterium belonging to the consortium of bioleaching microorganisms. Our phylogenetic data on the cluster of tRNA genes encoded in ICEAfe1 from A. ferrooxidans ATCC 23270 is consistent with the data presented by Tran et al.¹⁵ that includes the tRNA gene cluster from A. ferrivorans. Taken together, these data predict that these 2 arrays were acquired by HGT. These predictions led us to consider that the cluster of tRNA genes from A. ferrooxidans might contribute to adapting the tRNA pool of the bacterium to improve the translation efficiency of genes encoded within ICEAfe1. Our theoretical data, based on the tAI results calculated from the tRNA gene dosage, contradicts such an interpretation. Thus, if lowered global tRNA adaptation indices imply reduced fitness of the bacterium to the environment, one might expect the bacterium to display strategies to reduce the negative effect of the additional tRNA genes. Several alternatives for reducing this negative effect can be predicted, from the inactivation of tRNA genes to the complete excision and loss of the entire ICEAfe1. As implied by the fact that other strains of A. ferrooxidans do not contain ICEAfe1¹² and thus do not carry this cluster of tRNAs, we assume that these genes are not essential to the bacterium. However, although ICEAfe1 is active in excision, it is maintained in A. ferrooxidans without any detectable loss, implying that the genes encoded in this MGE likely contribute to the bacterium's fitness under laboratory conditions. The finding that all tested tRNAs encoded in ICEAfe1 are functional in either aminoacylation or protein synthesis supports the idea that these genes remain active. The fact that the majority of the tRNA genes encoded within ICEAfe1 are poorly expressed might imply that either the transcriptional promoters are poorly active or that these tRNAs are not required under the tested conditions. Their functions might be confined to certain specific conditions yet to be identified. Reduced transcription of the tRNA genes might account for the low expression levels of the majority of these tRNAs without affecting their functionality. Strikingly tRNA-OTH is the most highly expressed among all of the tRNAs encoded in ICEAfe1. Since its secondary structure resembles 2 fused tRNAs, we speculate that this tRNA-like RNA might have evolved from canonical tRNA genes that acquired a novel function in a process yet to be identified. Other tRNA-like RNAs that are part of a longer RNA regulate the activity of genes related to MGE in Bacillus cereus.²⁴ In conclusion, based on the theoretical data reported in this article, we propose that the tRNA genes in ICEAfe1 might produce a negative effect on global translation efficiency in A. ferrooxidans. However, despite the predicted global negative effect on tAI, some of these tRNA might contribute to the fitness of the cell, either through roles other than protein synthesis or by increasing the translation efficiency of genes encoded within ICEAfe1. Fixation of ICEAfe1 in the bacterial genome may be the result of the balance between the negative effects of the tRNAs and compensatory effects exerted by other gene(s) encoded in this element that might improve the fitness of the organism for its ecological niche.

Materials and methods

Bacterial strains and recombinant DNA technology

A. ferrooxidans strain ATCC 23270 was grown at 30°C to the late exponential phase in modified 9K culture medium at pH 1.6,²⁵ supplemented with FeSO₄·7H ₂O (120 mM) or sulfur at 1% as electron sources. All DNA manipulations were performed according to standard protocols.²⁶

RNA preparation for sequencing analysis (RNA-seq)

RNA isolation

Two methods were used for the RNA preparation. Total RNA was isolated from cells using either the RNeasy Mini Kit (Qiagen) or the TrueSeq Small RNA (Illumina) preparation kits (experiments 1 and 2 respectively). Two RNA samples were obtained from exponentially-growing A. ferrooxidans cultures to be used in each method. DNA contamination was eliminated with DNase I (ROCHE) according to manufacturer instructions. The RNA was stored at −80°C.

Library preparation and Illumina Solexa sequencing

When necessary, rRNA was eliminated from the total RNA using the MICROBExpress kit (Ambion). RNA integrity was electrophoretically analyzed. Libraries were generated according to the mRNA Sequencing Sample preparation guide (Illumina Solexa). Briefly, RNA was fragmented and reverse-transcribed using random primers, followed by double-stranded cDNA synthesis. An ‘A’ nucleotide was added at the 3′ end to prepare the DNA fragments for the ligation of adaptors, which have a single ‘T’ nucleotide overhanging at the 3′ end. This process prepared the fragments for hybridization to a single-read flow cell. The ligation products were gel-purified to select a size range of templates for downstream enrichment. The cDNA in the library was PCR-amplified with 2 primers that anneal to the ends of the adaptors. Cluster generation was performed according to manufacturer protocol (Illumina), and 36-nt single-end reads were generated on a Solexa Genome Analyzer at the Molecular and Cellular Imaging Center at Ohio State University, Columbus, United States of America or at AG All Genomics, Chile

Bioinformatic analysis of RNA-seq data

From experiment 1 (RNA extracted with RNeasy Mini Kit was used for RNA preparation, a total of 24,978,516 and 25,194,496 reads were generated from each RNA sample, respectively. We filtered the Solexa RNA sequencing data set by eliminating low-quality reads. Data filtering was performed using the Q-Solexa algorithm (score = −10 log10 (p / ( 1 -p)). Reads were aligned to the A. ferrooxidans ATCC 23270 reference genome obtained from the NCBI GenBank annotation (NC_011761) using the BLASTN program²⁷ with default settings, except for the following modifications: an aligned read required a minimum length of 15 nucleotides (coverage ≥ 70%). If a read could be aligned to multiple targets, it was automatically assigned to each target.

From experiment 2 (TrueSeq Small RNA protocol was used for RNA preparation), a total of 54,574,333 and 60,399,655 reads were generated from each RNA sample, respectively. The Trimmomatic software²⁸ was used to eliminating low-quality reads, short sequences and highly repeated sequences, obtaining a total of 50,208,386 and 54,359,690 reads respectively. Reads were used to construct a reference genome index using Bowtie-build2.²⁹ Next, alignment of the file GFF to the A. ferrooxidans ATCC 23270 reference genome obtained from NCBI using the TopHat software³⁰ with default settings was performed. Finally Cuffmerge software was run³¹ in the assemblies to create the annotated transcriptome.

The A. ferrooxidans ATCC 23270 genome contains a total of 3,007 predicted mRNAs and 86 non-coding RNAs (80 tRNAs and 6 rRNAs). To identify possible genes encoding tRNAs that contain introns, we used the ARAGORN platform²¹

Calculation of tRNA adaptation index (tAI)

We calculated tAI according to the method established by dos Reis¹⁸ (Equations 1 and 2), estimating the relative abundance of the tRNAs as a function of the number of copies of each tRNA gene in the genome of A. ferrooxidans ATCC 23270. The number of copies of each tRNA gene in the genome from A. ferrooxidans ATCC 23270 was determined using the tRNAscan-SE program²⁰ and ARAGORN platform²¹ in Linux.

(1) $W_i={\displaystyle \sum _{j=1}^{n_i}(1-S_{ij})tGN{C}_{ij}}$(1)

W_i is the relative adaptiveness of the i^th codon in the tRNA pool; n_i is the number of tRNA isoacceptors that recognize the i^th codon; tGCN_ij is the number of copies of the j^th tRNA gene that recognize the i^th codon; and S_ij is the selective constraint on the efficiency of codon-anticodon pairing.

The adaptation of a gene to the tRNA pool is calculated according to Equation 2

(2) $tA{I}_g={\left({\displaystyle \prod _{k=1}^{l_g}{w}_{ikg}}\right)}^{1/l_g}$(2)

W_i is defined as the quotient of W_i and W_max (W_i/W_max); ikg is the codon defined by the k^th triplet of gene g; and l_g is the length of gene g in codifying codons.

Calculation of tAI using RNA-seq data

The tRNA adaptation index was calculated as detailed above, with the exception that Equation 1 was modified into Equation 3, using empirical data on tRNA expression obtained from RNA-seq:

(3) $W_i={\displaystyle \sum _{j=1}^{n_i}(1-S_{ij})t{E}_{ij}}$(3)

tE_ij is the tRNA expression in arbitrary units of the j^th tRNA gene that recognizes the i^th codon.

The tRNA adaptation index, determined using either theoretical or empirical tRNA expression, was calculated using a Python script.³²

Detection of in vivo aminoacylated/deacylated tRNA

To determine the level of aminoacylated/deacylated tRNAs in vivo, total tRNA was extracted from A. ferrooxidans as described previously.³³ Briefly, periodate oxidation followed by β-elimination, electrophoresis on denaturing polyacrylamide gel and Northern blot were performed as described previously³⁴

Northern blot

DNA probes for Northern blot complementary to tRNA^Ser _{UGA 55} (serT from ICEAfe1), tRNA^Ser _{UGA 7)} (serTi from the chromosome containing the intron), and for tRNA^Val ₃₉, tRNA^Ala ₅₈ and tRNA^Trp ₆₁ from ICEAfe1 were designed. The following Northern blot protocol was used to analyze serT and serTi. RNA was separated on denaturing polyacrylamide gels (8%, 8 M urea) and transferred to a nylon membrane (Hybond -N⁺, Amersham) by a semi-dry transfer at 1 mA per cm² for 2 hours. RNA was fixed to the membrane at 75°C for 2 h or overnight. To label the DNA probe, polynucleotide kinase (New England Biolabs) and [γ³²P]-ATP(PerkinElmer) were used according to manufacturer instructions. The nylon membranes were incubated in 20 mL of pre-hybridization buffer (6X SSC, 0.1% SDS, 10X Denhardt's solution, 0.1 mg/mL salmon sperm DNA) at 45°C for 2 to 6 hours; this buffer was replaced with 10 mL of hybridization buffer (6X SSC, 0.1% SDS, 10X Denhardt's solution) and 2.5 µmol of the specific radiolabelled probe. The hybridization was incubated at 45°C for 6 to 20 h. Next, 3 incubations with wash buffer 1 (2X SSC and 0.5% SDS) at 45°C for 15 min and 2 washes with wash buffer 2 (2X SSC and 0.1% SDS) were performed. After the last washing step, the membrane was exposed to a screen for 2−72 hours and revealed by Phospho Imager scanning.

Alternatively, for tRNA^Val ₃₉, tRNA^Ala ₅₈ and tRNA^Trp ₆₁, the following procedure was used for Northern blot detection. The DIG Oligonucleotide Tailing Kit, 2nd Generation (Roche) was used according to manufacturer instructions to label the oligonucleotide probes with digoxigenin for generating short tails. Nylon membranes with RNA attached were incubated in 20 mL of pre-hybridization buffer at 50°C for 1 hr; this buffer was replaced with 20 mL of hybridization buffer (6X SSC, 0.1% SDS, 10X Denhardt's solution) and the specific labeled probes at 10 fmol/ml. The hybridizations were incubated at 50°C for 16 hr. Next, 3 washes with 50 mL of wash buffer at 50°C for 30 min were performed as described above. After the last wash, the probe-target hybrids were detected with an enzyme-linked immunoassay. All washing steps were performed at room temperature by shaking. The membranes were washed with 50 mL of wash buffer (0.1 M maleic acid, 0.15 M NaCl, pH 7.5, 0.3% (v/v) TWEEN 20) for 2 min. The blocking step was performed, using 50 mL of blocking solution (0.1 M maleic acid, 0.15 M NaCl, pH 7.5, 1% casein), for 30 min before the addition of the antibody solution (75 mU/ml anti-digoxigenin-AP (Roche) in blocking solution). The membranes were incubated with the antibody solution for 30 min and then washed twice with wash buffer. The membranes were equilibrated for 3 min in detection buffer (0.1 M Tris-HCl, 0.1 M NaCl, pH 9.5), and the BCIP/NBT Kit (Life Technologies) was used for chromogenic detection of probe-target hybrids.

Identification of tRNAs in Draft genomes

The tRNAs were identified using tRNAScan-SE against the draft genomes of Acidithiobacillus sp GGI221 (AEFB01, GCA_000179815.2), A. caldus SM1 (NC_015850.1), A. thiooxidans 19377 (AFOH01, GCA_000227215.2), A. thiooxidans 19377(2011) (AFOH01, GCF_000227215.2), A. thiooxidans A01 (AZMO01, GCA_000559045.1), A. ferrivorans CF27 (CCCS02, GCA_000750615.1), A. thiooxidans Licanantay (JMEB01, GCA_000709715.1), A. ferroxidans DLC5 (JNNH01, GCA_000732185.1), A. ferrooxidans YQH1 (LJBT01, GCA_001418795.1), A. ferrooxidans Hel18 (LQRJ01, GCA_001559335.1), A. ferrooxidans BY0502 (LVXZ01, GCA_001652185.1), A. ferrivorans PQ33 (LVZL01, GCA_001857665.1), A. thiooxidans DXSW (LWRY01, GCA_001705805.1), A. ferrivorans SS3 (NC_015942.1, GCA_000214095.3), A. thiooxidans BY02 (LWRZ01, GCA_001705725.1), A. thiooxidans A02 (LWSA01, GCA_001705645.1), A. thiooxidans DMC (LWSB01, GCA_001705625.1), A. thiooxidans GD13 (LWSC01, GCA_001705695.1), A. thiooxidans JYC17 (LWSD01, GCA_001705755.1), A. caldus MTH04 (LXQG01, GCA_001650235.1), A. caldus DX (LZYE01, GCA_001756675.1), A. caldus ZBY(LZYF01, GCA_001756725.1), A. caldus ZJ (LZYG01, GCA_001756745.1), A. caldus S1 (LZYH01, GCA_001756745.1), A. thiooxidans ZBY (LZYI01, GCA_001756595.1), A. ferrivorans YL15 (MASQ01, GCA_001685225.1).

Analysis of codon usage

For each genome codon usage was determined by JEMBOSS,³⁵ using the cusp algorithms and the coding sequences of each genomes. Coding sequences from non annotated genomes of Acidithiobacilli were determined using RASTserver,³⁶ by mean of automatic annotation. Comparison and visualization was performed using MeV software³⁷ by Hierarchical clustering

Phylogenetic reconstruction

The evolutionary history was inferred using the neighbor-joining method,³⁸ using a bootstrap of 1000 replicates.³⁹ The evolutionary distances were computed with the maximum composite likelihood method⁴⁰ using the number of base substitutions per site as units. The rate variation among sites was modeled with a gamma distribution (shape parameter = 1). All ambiguous positions were removed for each sequence pair. The evolutionary history and further analyses were performed using MEGA6 software.⁴¹

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Acknowledgments

We thank Dr. Michael Ibba and Dr. Assaf Katz for critical discussion of the manuscript. We also thank Dr. T. Hughes at the Department of Biology, University of Utah, Salt Lake City, USA for providing S. enterica serovar Typhimurium strain TH3238.

This work was supported by Fondo Nacional de Ciencias y Tecnología, Chile under grants1110203 and 1150834 to O.O. (http://www.fondecyt.cl); and the Fundación de Innovación Agraria, Chile under grants PYT20120056 and USA 1555 to MT. (http://www.fia.cl). A graduate fellowship from Programa de Mejoramiento de la Educación Superior, Chile was granted to P.A. A graduate fellowship from USACH was granted to J.M.

Additional information

Funding

This work was supported by Fondo Nacional de Ciencias y Tecnología, Chile under grants1110203 and 1150834 to O.O. (http://www.fondecyt.cl); and the Fundación de Innovación Agraria, Chile under grants PYT20120056 and USA 1555 to MT. (http://www.fia.cl). A graduate fellowship from Programa de Mejoramiento de la Educación Superior, Chile was granted to P.A. A graduate fellowship from USACH was granted to J.M.