DAMBE6: New Tools for Microbial Genomics, Phylogenetics, and Molecular Evolution

1DAMBE (Data analysis for molecular biology and evolution) is a comprehensive software package for sequence manipulation and analysis featuring a user-friendly interface and a variety of analytical functions in bioinformatics, phylogenetics, and descriptive and comparative genomics. It is often listed as one of the most widely used software packages in molecular phylogenetics. Version 6 of DAMBE (DAMBE6) added several new functions in genomic evolution and phylogenetics since DAMBE5 and updated and improved a number of existing functions.

2Translation initiation is often rate-limiting in bacteria and in bacteriophage. Translation initiation signals on mRNA in prokaryotes include the start codon decoded by fMet-tRNAfMet and Shine–Dalgarno sequence (SD) binding to the anti-SD (aSD) sequences at the 3′ end of small subunit ribosomal RNA (ssu rRNA).

3I will first clarify what constitute a good SD/aSD pairing. Structural determination showed that fMet-tRNAfMet and translation initiation factors can bind to 30S ribosome synergistically (binding of one facilitates the binding of others). The function of SD/aSD binding is to juxtapose the start codon against anticodon of fMet-tRNAfMet (Figure 1a). While many genes have their SD being AGGAGGU or part of it, many have different SDs (Figure 1b). Each SD has its specific optimal distance (D) between the SD and start codon, for example, the optimal D is D1 for SD1 and D2 for SD2 in Figure 1a. One real case involves Escherichia coli rpsQ gene (Figure 1c) which has 2 putative SDs, AAGG and GGUG (Figure 1c). However, we note that the 2 SDs in Figure 1a have the same DtoStart defined as the distance between the 3′ end of ssu rRNA and the start codon (Figure 1a). DtoStart is strongly constrained within a narrow range in a variety of bacterial species from the gram-negative E. coli to the gram-positive Bacillus subtilis, suggesting that DtoStart is a better and more general index for measuring optimal positioning of SD/aSD pairing than D1 or D2 in Figure 1a (because D1 or D2 are SD-specific). The 2 putative SDs in the rpsQ gene (Figure 1c) bind to different aSDs but both have similar DtoStart values (15 and 14, respectively, Figure 1c). It is meaningless to state that the optimal distance between an SD and start codon is 5 or 10 nucleotides (nt) without specifying what SD is. An SD can be very close to the start codon, as in the case of pflB with only 4 nt in between (Figure 1) or far apart as in the case of adk with 11 nt in between (Figure 1b). What is common among all of them is that they all have similar DtoStart (Figure 1b, c).

4Given an annotated prokaryotic genome, DAMBE can 1) extract the sequence upstream of each coding sequences (CDSs), for example, 20 nt immediately upstream of the initiation codon, 2) identify the putative SDs in all protein-coding sequences, and 3) output a variety of summary statistics to show which gene has a strong and optimally positioned SD/aSD. This is illustrated with data from E. coli (Figure 2). Most E. coli SD/aSD matches have DtoStart = 13 (Figure 2a). The frequency increases sharply from DtoStart = 11 to DtoStart = 12, but decreases more gradually on the right side (Figure 2a). This feature is common among diverse bacterial species. Most E. coli SDs are confined within a narrow range within 20 nt upstream of the start codon (Figure 2b). Three most frequent E. coli SDs are AGGA, GGAG, and GAGG which overlap to form the longer and better known motif of AGGAGG (Table 1). An overwhelming majority of SDs are 4-nt long (Figure 2d), although many studies suggest that longer SDs are more efficient in localizing the start codon.

5The 13 nt at the 3′ end of E. coil ssu RNA are differentially involved in SD/aSD, with some sites (e.g., UCCUC at sites 3–7) involved in SD/aSD pairing more frequently than others. A nucleotide substitution at one of these sites will affect SD/aSD pairing (and consequently translation initiation) of thousands of genes. We therefore expect these sites to be extremely conserved due to the constraints of so many genes. A corollary from this framework of reasoning is that other sites constrained by fewer SD/aSD pairing would be more tolerated.

6Translation initiation of most prokaryotic genes requires well-positioned SD/aSD base-pairing, although SD/aSD base-pairing is not always essential for translation in E. coli and Chlamydomonas reinhardtii chloroplasts, and for translating leaderless genes that have no SD sequence. The strength and position of SD/aSD base-pairing do strongly affect translation initiation in many genes. The tools offered in DAMBE facilitate large-scale study of SD/aSD coevolution as different species do have different 3′ end of small subunit rRNA (3′ TAIL) demanding different SD/aSD pairing dynamics.

7SD and the start codon constitutes key translation initiation signals on mRNA to be recognized by ribosomes and initiation tRNA, respectively. Having these signals embedded in secondary structure decreases translation initiation efficiency, especially in highly expressed genes.

8DAMBE includes functions to extract CDSs and their upstream and downstream sequences, and can use a sliding window to compute minimum folding energy (MFE) which measures stability of local secondary structure. The MFE profile shows a dramatic decrease in secondary structure around the start codon and the surrounding region, which is particularly pronounced in highly expressed genes (Figure 3a). A similar trend is observed in stop codons, but it is overshadowed by a much stronger increase in secondary structure stability about 30 nt downstream of the stop codon (Figure 3b), which is likely due to the hairpin involved in the rho-independent termination.

9In prokaryotes, some genes are closely spaced with sequence configurations such as -AUGA- (where UGA is the stop codon of the upstream gene and AUG is the start codon of the downstream gene) and -UAAUG- (where UAA is the stop codon of the upstream gene and AUG is the start codon of the downstream gene). The patterns in Figure 3 are from genes with an intergenic sequence of at least 100 nt to avoid confounding the MFE pattern near the start codon and that at the stop codon. DAMBE can optionally include gene location information in the output sequence file.

10Many gene-specific codon usage indices have already been formulated and improved, including CAI, tAI and several indices that are based on CDSs only, such as Nc and its improved versions. The first 2 have been used frequently as proxies for translation elongation efficiency, but they both have major problems.

11The problem with tAI is that we often cannot infer which tRNA favors which synonymous codon. For example, inosine is expected to pair best with C and U, less with A (partly because of the bulky I/A pairing involving 2 purines), and not with G, but this is not true with tRNAVal/IAC from rabbit liver which pairs better with GUG codon than with other synonymous codons. Similarly, the B. subtilis genome codes a tRNAAla/GGC for decoding GCY codons, but the GCC codon which forms Watson–Crick base pair with the anticodon is not used as frequently as the GCU codon which wobble-pairs with the anticodon. Furthermore, codon–anticodon base pairing is known to be context-dependent, for example, a wobble cmo5U in the anticodon of tRNAPro, tRNAAla, and tRNAVal can read all 4 synonymous codons in the respective codon family, but the same cmo5U in tRNAThr cannot read C-ending codons. For this reason, the optimal codon usage is likely better approximated by the codon usage of highly expressed genes than what we can infer based on codon–anticodon pairing.

12CAI also has problems. In particular, it ignores background mutation bias which can result in misinterpretation of tRNA-mediated selection. Take for example the Ala codon subfamily GCR (where R stands for either A or G). The frequencies of GCA and GCG in E. coli HEGs, as compiled and distributed with EMBOSS, are 1973 and 2654, respectively, which may lead one to think that E. coli translation machinery prefer GCG over GCA. However, GCA is relatively more frequent in E. coli HEGs than in E. coli non-HEGs. This suggests that mutation bias favors GCG, but tRNA-mediated selection favors GCA. This interpretation is corroborated by the E. coli genome encoding three tRNAArg genes for GCR codons, all with a UGC anticodon forming perfect Watson–Crick base pair with codon GCA.

13DAMBE implements a new index of translation elongation (ITE) which incorporates both tRNA-mediated selection and background mutation bias and fits protein production better than CAI or tAI. CAI is a special case of ITE when background mutation bias is absent. There are 4 variations of ITE with different treatment of synonymous codon families (Figure 4). The first is to treat R-ending and Y-ending codon groups as if they are separate codon families, with reasons for such a treatment outlined before. The second (the default) is to separate compound 8-fold codon families into 2 separate 4-fold codon families, and 6-fold codon families into 2 codon families with 4 and 2 synonymous codons each. Such separation is reasonable because the 4-codon and 2-codon synonymous families are translated by different tRNAs. The third is to lump all synonymous codons into one codon family. The fourth is to use only R-ending codons because, in some species such as E. coli, codon bias is strong in R-ending synonymous codons but weak in Y-ending synonymous codons. ITE has been used to facilitate studies on translation initiation and elongation in bacteriophages and coevolution between stop codons and release factors in bacteria.

14Pairwise sequence alignment (PSA) by dynamic programming is guaranteed to generate one of the optimal alignments, but multiple sequence alignment (MSA) by dynamic programming is often not practical. The commonly used progressive alignment along a guide tree often results in poor alignment for highly diverged sequences in spite of many iterations to update the guide tree and the alignment, plaguing all subsequent phylogenetic analysis. One way to avoid this problem is to use only PSA to reconstruct phylogenetic trees, which can only be done with distance-based methods. DAMBE implements such a phylogenetic method (PhyPA) based only on pairwise alignment. I compared the accuracy of PhyPA against the combination of maximum likelihood method and MSA (the ML+MSA approach), using nucleotide, amino acid, and codon sequences simulated with different topologies and tree lengths. Surprisingly, the fast PhyPA method consistently outperforms the slow ML+MSA approach for highly diverged sequences even when all optimization options were turned on for the ML+MSA approach. Only when sequences are not highly diverged (i.e., when a reliable MSA can be obtained) does the ML+MSA approach outperforms PhyPA. The PhyPA method implemented in DAMBE also includes 2 approaches making use of multi-gene data sets to derive phylogenetic support for subtrees equivalent to resampling techniques such as bootstrapping and jackknifing.

15PhyPA can also be used to characterize phylogenetic structure of gene families or identify pseudogenes. For example, after using PSI-BLAST to obtain hundreds or even thousands of sequences with remote homology, one can use PhyPA to reconstruct a phylogenetic tree based on pairwise alignment and the resulting tree structure will reflect the number of, and relationship among, gene families. PhyPA can also be used to quickly identify candidate pseudogenes. For example, zebrafish (Danio rerio) has 12292 tRNA genes according to GtRNAdb (http://gtrnadb.ucsc.edu/). Phylogenetic analysis with PhyPA revealed many tRNAs to have extraordinarily long branches and most likely are tRNA pseudogenes. This is exemplified by partial phylogenetic trees from 1478 tRNALys genes (Figure 5, top) and from 1162 tRNAGly genes (Figure 5, bottom).

16Many other functions in DAMBE have been updated and improved. A variety of statistical tests have been added to position weight matrix for motif characterization which has been applied in characterizing the splicing signal strength in yeast and vertebrates. The function for handling multiple files in a number of phylogenetic analyses is particularly useful with sequence simulation that includes indels. Such simulations often generate a large number of unaligned sequence files. With only a few clicks, DAMBE will be able to align all these files, reconstruct phylogenetic trees and compare the differences between the resulting trees and the true tree used in sequence simulation. Other functions that have been improved include hidden Markov models for protein secondary structure prediction, Gibbs sampler for de novo motif discovery, and self-organizing map for nonlinear clustering of transcriptomic data (;, p. 231–250).

17In short, DAMBE is a comprehensive software workbench in molecular biology, phylogenetics, and evolution, with new functions continuously added to empower researchers to perform leading-edge data analysis in prokaryotic genomic data to solve practical research problems. DAMBE is user-friendly with a variety of graphic functions, which makes it ideal not only for research, but also for teaching. DAMBE is available free of charge from http://dambe.bio.uottawa.ca, where a set of laboratory tutorials designed for teaching can be found. DAMBE is a Windows program, but may run on Linux and Macintosh computers.

18This work was supported by the Discovery Grant of Natural Science and Engineering Research Council of Canada (NSERC, RGPIN/261252–2013).