Ancestral Genome Organization: An Alignment Approach

1. Introduction

During evolution, genomes continually accumulate mutations. In addition to base mutations and short insertions or deletions, genome-scale changes affect the overall gene content and organization of a genome. Evidence of these latter kinds of changes are observed by comparing the completely sequenced and annotated genomes of related species. Genome-scale changes can be subdivided into two categories: (1) the rearrangement operations that shuffle gene order such as inversion, and (2) the content-modifying operations that affect the number of gene copies, such as gene insertion, loss and duplication. In particular, gene duplication is a fundamental process in the evolution of species (Ohno 1970), especially in eukaryotes (Blomme et al., 2006; Cotton and Page, 2005; Eichler and Sankoff, 2003; Hahn et al., 2007; Lynch and Conery, 2000; Wapinski et al., 2007), where it is believed to play a leading role for the creation of novel gene function. In parallel, gene losses through pseudogenization and segmental deletions, appear generally to maintain a minimum number of functional gene copies (Blomme et al., 2006; Cotton and Page, 2005; Demuth et al., 2006; Eichler and Sankoff, 2003; Hahn et al., 2007; Lynch and Conery, 2000; Ohno, 1970). Transfer RNAs (tRNAs) are typical examples of gene families that are continually duplicated and lost (Bermudez-Santana et al., 2010; Rogers et al., 2010; Tang et al., 2009; Withers et al., 2006). Indeed, tRNA clusters (or operons in microbial genomes) are highly dynamic and unstable genomic regions. In Escherichia coli for example, the rate of tRNA gene duplication/loss events has been estimated to be about one event every 1.5 million years (Bermudez-Santana et al., 2010; Withers et al., 2006).

One of the main goals of comparative genomics is to infer evolutionary histories of gene families, based on the comparison of the genomic organization of extant species. Having an evolutionary perspective of gene families is a key step towards answering many fundamental biological questions. For example, tRNAs are essential to establishing a direct link between codons and their translation into amino-acids. Understanding how the content and organization of tRNAs evolve is essential to the understanding of the translational machinery, and in particular, the variation in codon usage among species (Dong et al., 2006; Kanaya et al., 1999).

In the genome rearrangement approach to comparative genomics, a genome is modeled as one or many (in case of many chromosomes) linear or circular sequences of genes (or other building blocks of a genome). When each gene is present exactly once in a genome, sequences can be represented as permutations. In the most realistic version of the rearrangement problem, a sign is associated with a gene, representing its transcriptional orientation. The pioneering work of Hannenhalli and Pevzner (Hannenhalli and Pevzner, 1995, 1999) has led to efficient algorithms for computing the inversion and/or translocation distance between two signed permutations. Since then, many other algorithms have been developed to compare permutations subject to various rearrangement operations and based on different distance measures. These algorithms have then been used from a phylogenetic perspective to infer ancestral permutations (Bourque and Pevzner, 2002; Chauve and Tannier, 2008; Ma et al., 2007; Moret et al., 2001; Sankoff and Blanchette, 1997) and evolutionary scenarios on a species tree. An extra degree of difficulty is introduced in the case of sequences containing multiple copies of the same gene, as the one-to-one correspondence between copies is not established in advance. A review of the methods used for comparing two ordered gene sequences with duplicates can be found in El-Mabrouk (2005) and Fertin et al. (2009). They can be grouped into two main classes. The “Match-and-Prune” model aims at transforming strings into permutations, so as to minimize a rearrangement distance between the resulting permutations. On the other hand, the “Block Edit” model consists of performing the minimum number of “allowed” rearrangement and content-modifying operations required to transform one string into the other. Most studied distances and ancestral inference problems in this category are NP-complete (Fertin et al., 2009).

In this paper, we focus on comparing two ordered gene sequences with duplicates. We develop our methodology for a model of evolution restricted to content-modifying operations, more specifically to duplication and loss, and then show how to extend it to other content-modifying operations, and to non-overlapping inversions. From a combinatorial point of view, the main consequence of ignoring rearrangement operations is the fact that gene organization is preserved, which allows us to reformulate the problem of comparing two gene orders as an alignment problem. Notice however that only the most recent events that have not been obscured by subsequent events are visible in such an alignment. Another advantage of a model restricted to duplication and loss is that a unique ancestral genome can directly be inferred from a duplication-loss scenario attached to a given alignment, as duplications and losses are asymmetrical operations that are applicable to one of the two aligned sequences.

Although alignments are a priori simpler to handle than rearrangements, there is no direct way of inferring optimal alignments together with a related duplication-loss scenario for two gene orders, as detailed in Section 4. Even our simpler goal of finding an alignment is fraught with difficulty as a naive branch-and-bound approach to compute such an alignment is non-trivial; trying all possible alignments with all possible duplication and loss scenarios for each alignment is hardly practicable. As it is not even clear how, given an alignment, we can assign duplications and losses in a parsimonious manner, we present in Section 4.1 a pseudo-boolean linear programming (PBLP) approach to search for the optimal alignment along with an optimal scenario of duplications and losses. The disadvantage of the approach is that, in the worst case, an exponential number of steps could be used by our algorithm. On the other hand, we show in Section 6.2 that for real data, and larger simulated genomes, the running times are quite reasonable. Further, the PBLP is flexible in that a multitude of weighting schemes for losses and duplications could be employed to, for example, favor certain duplications over others, or allow for gene conversion.

In Section 5, we extend our initial model of evolution to handle additional content-modifying operations, as well as non-overlapping inversions, and we show how to relax some constraints regarding the visible (i.e., non-intersecting) nature of operations. In Section 6.1, we apply our algorithm in a phylogenetic context to infer the evolution of stable RNA (tRNA and rRNA) gene content and organization in various genomes from the genus Bacillus, a so-called “low G+C” gram-positive clade of Firmicutes that includes the model bacterium B. subtilis as well as the agent of anthrax. Stable RNA operon organization in this group is interesting because it has relatively fewer operons that are much larger and contain more segmental duplicates than other bacterial groups. We obtained results leading to various biological insights, such as more accurate quantification of ribosomal RNA operon proliferation, their role in altering tRNA gene content, and evidence of tRNA gene class conversion.

2. Research Context

The evolution of g genomes is often represented by a phylogenetic (or species) tree T, binary or not, with exactly g leaves, each representing a different genome. When such a species tree T is known for a set of species, then we can use the gene order information of the present-day genomes to infer gene order information of ancestral genomes identified with each of the internal nodes of the tree. This problem is known in the literature as the “small” phylogeny problem, in contrast to the “large” phylogeny problem which is one of finding the actual phylogenetic tree T.

Although our methods may be extended to arbitrary genomes, we consider single chromosomal (circular or linear) genomes, represented as gene orders with duplicates. More precisely, given an alphabet Σ where each character represents a specific gene family, a genome or string is a sequence of characters from Σ where each character may appear many times. To simplify our explanation we ignore the orientation (sign) of the genes. However, our methodology is easily extendable to inclusion of this information. For example, given Σ = {a, b, c, d, e}, A = “ababcd” is a genome containing two gene copies from the gene family identified by a, two genes from the gene family b, and a single gene from each family c and d. A gene in a genome A is a singleton if it appears exactly once in A (for example c and d in A), and a duplicate otherwise (a and b in A above).

Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document} be a set of “allowed” evolutionary operations. The set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document} may include organizational operations such as Reversals (R) and Transpositions (T), and content-modifying operations such as Duplications (D), Losses (L) or Insertions (I). For example, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document}  = {R, D, L} is the set of operations in an evolutionary model involving reversals, duplications and losses. In the next section, we will formally define the operations involved in our model of evolution.

Given a genome A, a mutation on A is characterized by an operation O from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document} , the substring of A that is affected by the mutation, as well as possibly other characteristics such as the position of the re-inserted removed (in case of transposition) or duplicated substring. For simplicity, consider a mutation O(k) to be characterized solely by the operation O from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document} , and the size k of the substring affected by the mutation. Consider c(O(k)) to be a cost function defined on mutations. Finally, given two genomes A and X, an evolutionary history O_A→X from A to X is a sequence of mutations (possible of length 0) transforming A into X.

Let A, X be two strings on Σ with A being a potential ancestor of X, meaning that there is at least one evolutionary history \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$O_{A \rightarrow X} = \{O_1 (k_1) , \ldots , O_l (k_l) \} $$ \end{document} from A to X. Then the cost of O_A→X is: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C (O_{A \rightarrow X}) = \sum_{i = 1}^l \; c (O_i (k_i))\end{align*} \end{document}

Now let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}_{A \rightarrow X}$$ \end{document} be the set of possible histories transforming A into X. Then we define: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C (A \rightarrow X) = \min_{O_{A \rightarrow X} \in {\cal O}_{A \rightarrow X}} C (O_{A \rightarrow X})\end{align*} \end{document}

Then, the small phylogeny problem can be formulated as one of finding strings at internal nodes of a given tree T that minimize the total cost: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}C (T) = \sum_{\hbox {all branches} \; b_i \; \hbox{of} \; T} C (X_{i , 1} \rightarrow X_{i , 2})\end{align*} \end{document}

where X_i,1, X_i,2 are the strings labeling the nodes of T adjacent to the branch b_i, with the node labeled X_i,1 being the parent of the node labeled X_i,2.

For most restrictions on genome structure and models of evolution, the simplest version of the small phylogeny problem—the median of three genomes—is NP-hard (Caprara, 1999; Pe'er and Shamir, 1998; Tannier et al., 2009). When duplicate genes are present in the genomes, even finding minimum distances between two genomes is almost always an NP-Hard task (Jiang, 2010). In this paper, we focus on cherries of a species tree (i.e., on subtrees with two leaves). The optimization problem we consider can be formulated as follows:

Two Species Small Phylogeny Problem:

Input: Two genomes X and Y.

Output: A potential common ancestor A of X and Y minimizing C(A → X) + C(A → Y).

Solving the Two Species Small Phylogeny Problem (2-SPP) can be seen as a first step towards solving the problem on a given phylogenetic tree T. The most natural heuristic to the Small Phylogeny Problem, that we will call the SPP-heuristic, is to traverse T depth-first, and to compute successive ancestors of pairs of nodes. Such a heuristic can be used as the initialization step of the steinerization method for SPP (Sankoff and Blanchette, 1997; Blanchette et al., 1997). The sets of all optimal solutions output by an algorithm for the 2-SPP applied to all pairs of nodes of T (in a depth-first traversal) can alternatively be used in an iterative local optimization method, such as the dynamic programming method developed in Kovac et al. (2011).

3. The Duplication And Loss Model Of Evolution

Our evolutionary model accounts for two operations, Duplication (denoted D) and Loss (denoted L). In other words \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O} = \{D , L \} $$ \end{document} , where D and L are defined as follows. Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [i \ldots i + k]$$ \end{document} denote the substring \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_{i} X_{i + 1} \cdots X_{i + k}$$ \end{document} of X.

− D: A Duplication of size k+1 on \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X = X_1 \cdots X_{i} \cdots X_{i + k} \cdots X_j X_{j + 1} \cdots X_n$$ \end{document} is an operation that copies a substring \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [i \ldots i + k]$$ \end{document} to a location j of X outside the interval [i, i + k] (i.e., preceding i or following i + k). In the latter case, D transforms X into \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}X^{\prime} = X_1 \cdots \underline{X_{i} \cdots X_{i + k}} \cdots X_{j - 1} \underline{X_{i} \cdots X_{i + k}} X_{j + 1} \cdots X_n\end{align*} \end{document}

We call the original copy \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [i \ldots i + k]$$ \end{document} the origin, and the copied string the product of the duplication D.

− L: A Loss of size k is an operation that removes a substring of size k from X.

Notice that gene insertions could be considered in our model as well. In particular, our linear programming solution is applicable to an evolutionary model involving insertions, in addition to duplications and losses. We ignore insertions for two main reasons: (1) insertions and losses are two symmetrical operations that can be interchanged in an evolutionary scenario. Distinguishing between insertions and losses may be possible on a phylogeny, but cannot be done by comparing two genomes; (2) gene insertions are usually due to lateral gene transfer, which may be rare events compared to nucleotide-level mutations that eventually transform a gene into a pseudogene.

As duplication and loss are content-modifying operations that do not shuffle gene order, the Two Species Small Phylogeny Problem can be posed as an alignment problem. However, the only operations that are “visible” on an alignment are the events on an evolutionary history that are not obscured by subsequent events. Moreover, as duplications and losses are asymmetrical operations, an alignment of two genomes X and Y does not reflect an evolutionary path from X to Y (as operations going back to a common ancestor are not defined), but rather two paths going from a common ancestor to both X and Y. A precise definition follows.

Definition 1. Let X and Y be two genomes. A visible history of X and Y is a triplet (A, O_A→X, O_A→Y) where A is a potential ancestor of both X and Y, and O_A→X (respectively O_A→Y) are evolutionary histories from A to X (respectively from A to Y) verifying the following property: Let D be a duplication in O_A→X or O_A→Y copying a substring S. Let S₁ be the origin and S₂ be the product of D. Then D is not followed by any other operation inserting (by duplication) genes inside S₁ or S₂, or removing (by loss) genes from S₁ or S₂. We call a visible ancestor of X and Y a genome A belonging to a visible history (A, O_A→X, O_A→Y) of X and Y.

We now define an alignment of two genomes.

Definition 2. Let X be a string on Σ, and let Σ⁻ be the alphabet Σ augmented with an additional character “−”. An extension of A is a string A⁻ on Σ⁻ such that removing all occurrences of the character “−” from A⁻ leads to the string A.

Definition 3. Let X and Y be two strings on Σ. An alignment of size α of X and Y is a pair (X⁻, Y⁻) extending (X, Y) such that |X⁻| = |Y⁻| = α, and for each i, 1 ≤ i ≤ α, the two following properties hold:

− If X_i⁻≠“−” and Y _i ⁻≠“−” then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_i^{-} = Y_i^{-}$$ \end{document} ;

−  \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_i^{-}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y_i^{-}$$ \end{document} cannot be both equal to “−”.

Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A} = (X^{-} , Y^{-})$$ \end{document} be an alignment of X and Y of size α. It can be seen as a 2 × α matrix, where the ith column \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}_i$$ \end{document} of the alignment is just the ith column of the matrix. A column is a match iff it does not contain the character ‘−’, and a gap otherwise. A gap \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left[{X_i \atop -} \right]$$ \end{document} is either part of a loss in Y, or part of a duplication in X (only possible if the character X_i is a duplicate in X). The same holds for the column \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\left[{- \atop Y_j} \right]$$ \end{document} . An interpretation of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}$$ \end{document} as a sequence of duplications and losses is called a labeling of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}$$ \end{document} . The cost of a labeled alignment is the sum of costs of all underlying operations.

As duplications and losses are asymmetric operations that are applied explicitly to one of the two strings, each labeled alignment \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}$$ \end{document} of X and Y leads to a unique common ancestor A for X and Y. The following theorem shows that this, and the converse is true (Fig. 1).

OPEN IN VIEWER

FIG. 1.

Left: a labeled alignment between two strings X = ababacd and Y = abaad. Right: the ancestor A and two histories respectively from A to X and from A to Y obtained from this alignment. The order of operations in the history from A to Y is arbitrary.

Theorem 1. Given two genomes X and Y, there is a one-to-one correspondence between labeled alignments of X and Y and visible ancestors of X and Y.

Proof. Let A be a visible ancestor of X and Y and (A, O_A→X, O_A→Y) be a visible history of X and Y. Then construct a labeled alignment \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}$$ \end{document} of X and Y as follows:    ■

1. Initialization: Define the two strings X⁻ = Y⁻ = A on Σ⁻, and define \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}$$ \end{document} as an alignment with all matches between X⁻ and Y⁻ (i.e., self-alignment of A).

As (A, O_A→X, O_A→Y) is a visible history of X and Y, by definition the origins and products of duplications remain unchanged by subsequent operations on each of O_A→X and O_A→Y. Therefore, all intermediate labellings remain valid in the final alignment. Therefore, the constructive method described above leads to a labeled alignment of X and Y.

On the other hand, a substring \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_i \cdots X_j$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y_i \cdots Y_j$$ \end{document} ) that is labeled as a duplication in X (resp. Y) should not be present in A, as it is duplicated on the branch from A to X (resp. from A to Y). Also, a substring \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X_i \cdots X_j$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y_i \cdots Y_j$$ \end{document} ) that is labeled as a loss in Y (resp. X) should be present in A, as it is lost on the branch from A to Y (resp. from A to X). This implies an obvious algorithm to reconstruct a unique ancestor.     ■

In other words, Theorem 1 states that the Two Species Small Phylogeny Problem reduces in the case of the Duplication-Loss model of evolution to the following optimization problem.

Duplication-Loss Alignment Problem:

Input: Two genomes X and Y on Σ.

Output: A labeled alignment of X and Y of minimum cost.

4. Method

Although alignments are a priori simpler to handle than rearrangements, a straightforward way to solve the Duplication-Loss Alignment Problem is not known. We show in the following paragraphs, that a direct approach based on dynamic programming leads, at best, to an efficient heuristic, with no guarantee of optimality.

Let X be a genome of size n and Y be a genome of size m. Denote by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [1 \ldots i]$$ \end{document} the prefix of size i of X, and by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y [1 \ldots j]$$ \end{document} the prefix of size j of Y. Let C(i, j) be the minimum cost of a labeled alignment of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [1 \ldots i]$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y [1 \ldots j]$$ \end{document} . Then the problem is to compute C(m, n).

DP: A natural idea would be to consider a dynamic programming approach (DP), computing C(i, j), for all 1 ≤ i ≤ n and all 1 ≤ j ≤ m. Consider the variables M(i, j), D_X(i, j), D_Y (i, j), L_X(i, j) and L_Y (i, j) which reflect the minimum cost of an alignment \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}_{i , j}$$ \end{document} of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [1 \ldots i]$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y [1 \ldots j]$$ \end{document} satisfying respectively, the constraint that the last column of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal A}_{i , j}$$ \end{document} is a match, a duplication in X, a duplication in Y, a loss in X, or a loss in Y. Consider the following recursive formulae. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} - M (i , j) = \begin{cases} \begin{matrix} C (i - 1 , j - 1) \quad {\rm if} \ X [i] = Y [j] \\ + \infty \hfill \quad \quad \quad{\rm otherwise} \hfill \end{matrix} \end{cases} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}- L_X (i , j) = \min \nolimits_{0 \leq k \leq i - 1} [C (k , j) + c (L (i - k))]\end{align*} \end{document}

(the corresponding formula holds for L_Y (i, j)) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} - D_X (i , j) = \begin{cases} \begin{matrix} + \infty \hfill\quad\quad\qquad\qquad\qquad\qquad\qquad\qquad {\rm if} \ X [i] {\hbox {is a singleton}} \\ \min \nolimits_{l \leq k \leq i - 1} [C (k , j) + c (D (i - k))] \quad{\rm otherwise ,}\hfill \end{matrix} \end{cases} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [l \ldots i]$$ \end{document} is the longest suffix of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [1 \ldots i]$$ \end{document} that is a duplication

(the corresponding formula holds for D_Y (i, j)).

The recursions for D_X and D_Y imply that duplicated segments are always inserted to the right of the origin. Unfortunately, such an assumption cannot be made while maintaining optimality of the alignment. For example, given the cost c(D(k)) = 1 and c(L(k)) = k, the optimal labeled alignment of S₁ = abxabxab and S₂ = xabx aligns ab of S₂ with the second ab of S₁, leading to an optimal history with two duplications inserting the second ab of S₁ to its left and to its right. Such an optimal scenario cannot be recovered by DP.

DP-2WAY: As a consequence of the last paragraph, consider the two-way dynamic programming approach DP-2WAY that computes D_X(i, j) (resp. D_Y (i, j)) by looking for the longest suffix of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$X [1 \ldots i]$$ \end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$Y [1 \ldots j]$$ \end{document} ) that is a duplication in the whole genome X (resp. Y). Unfortunately, DP-2WAY may lead to invalid cyclic evolutionary scenarios, as the same scenario may involve two duplications: one with origin S₁ and product S₂, and one with origin S₂ and product S₁, where S₁ and S₂ are two duplicated strings. This is described in more detail in Section 4.1.

DP-2WAY-UNLABELED: The problem mentioned above with the output of DP-2WAY is not necessarily the alignment itself, but rather the label of the alignment. As a consequence, one may think about a method, DP-2WAY-UNLABELED, that would consider the unlabeled alignment output by DP-2WAY, and label it in an optimal way (e.g., find an evolutionary scenario of minimum cost that is in agreement with the alignment). Notice first that the problem of finding a most parsimonious labeling of a given alignment, is not a priori an easy problem, and there is no direct and simple way to do it. Moreover, although DP-2WAY-UNLABELED is likely to be a good heuristic algorithm to the Duplication-Loss Alignment Problem, it would not be an exact algorithm, as an optimal cyclic alignment is not guaranteed to have a valid labeling leading to an optimal labeled alignment. Figure 2 shows such an example; an optimal cyclic duplication and loss scenario can be achieved by both alignments (5 operations), while the optimal acyclic scenario can only be achieved by the alignment of Figure 2b.

OPEN IN VIEWER

FIG. 2.

Alignments for strings X = “zxyzxyaxbwaxb” and Y = “zxyxwb”. We consider the following cost: c(D(k)) = 1 and c(L(k)) = k for any integer k. Matches are denoted by a vertical bar, losses denoted by an “L”, and duplications denoted by bars, brackets, and arrows. Alignment (a) yields 6 operations (3 duplications and 3 losses) and implies ancestral sequence “zxyxwaxbz”, while (b) yields 5 operations (3 duplications and 2 losses) and implies ancestral sequence “zxyaxwbz”.

5. Extended Models Of Evolution

This section shows how the PBLP may be adapted to handle additional content-modifying operations, as well as some rearrangement events. We will also consider the possibility of applying general distance methods to specific substrings for which we may want to consider less visible events. For any considered set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal O}$$ \end{document} of allowed operations, the goal will be to add the appropriate equations to the PBLP in order to output an alignment reflecting a most parsimonious visible history.

Before describing our new operations, we need to define the notion of a visible history in a more general context, as Definition 1 is restricted to visible duplications. Definition 4 given bellow is general. Recall that a duplication operation has an origin and a product. In the case of rearrangement operations and losses, we consider the origin to be the empty string, and the product to be the resulting string. So in the case of inversions, the product is the inverted string, and in the case of losses both the origin and product are empty.

Definition 4. Let X and Y be two genomes. A visible history of X and Y is a triplet (A, O_A→X, O_A→Y) where A is a potential ancestor of both X and Y, and O _{A →X} (respectively O_A→Y) are evolutionary histories from A to X (respectively from A to Y) verifying the following property: Let O be an operation in O _{A →X}. Let S₁ be the origin and S₂ be the product of O. Then O is not followed by any other operation modifying (by insertion, deletion, substitution or rearrangement) S₁ or S₂.

6. Applications

6.1. Evolution of stable RNA gene content and organization in Bacillus

The stable RNAs are chiefly transfer RNAs (tRNAs) and ribosomal RNAs (rRNAs), which are essential in the process of translating messenger RNAs (mRNAs) into protein sequences. They are usually grouped in the genome within clusters (or operons in the case of microbial genomes), representing highly repetitive regions, causing genomic instability through illegitimate homologous recombination. In consequence, stable RNA families are rapidly evolving by duplication and loss (Bermudez-Santana et al., 2010; Rogers et al., 2010; Withers et al., 2006).

We applied Algorithm 1 (without the extensions of Section 5) in a phylogenetic context, using the SPP-heuristic described at the end of Section 2, to analyze the stable RNA content and organization of 5 Bacillus lineages: Bacillus cereus ATCC 14579 (NC4722), Bacillus cereus E33L (NC6274), Bacillus anthracis (NC7530), Bacillus licheniformis ATCC 14580 (NC6322) and Bacillus subtilis (NC964). The overall number of represented RNA families in these genomes is around 40, and the total number of RNAs in each genome is around 120. Our PBLP algorithm processes each pair of these genomes in a few seconds. We used the following cost for duplications and losses: c(D(k)) = 1 and c(L(k)) = k, for any integer k representing the size of an operation. The phylogeny in Figure 3 reflects the NCBI taxonomy. Each leaf is labeled by a block representation of the corresponding genome. Details on each colored block is given in Figure 4.

OPEN IN VIEWER

FIG. 3.

An inferred evolutionary history for the five Bacillus lineages identified with each of the five leaves of the tree. Circular bacterial genomes have been linearized according to their origin of replication (e.g., the endpoints of each genome is its origin of replication). Bar length is proportional to the number of genes in the corresponding cluster. A key for the bars is given in Figure 4, except for white bars that represent regions that are perfectly aligned inside the two groups (cereus, anthracis) and (licheniformis, subtilis), but not between the two groups. More specifically, the 27 duplications and losses reported on the top of the tree are obtained from the alignment of these white regions. Finally, each Δ represents a loss and each ∇ is an insertion.

OPEN IN VIEWER

FIG. 4.

The RNA gene clusters represented by each colored bar of Figure 3. Ribosomal RNAs are identified by their families' names: 16S, 23S, 5S. Each tRNA is identified by its anticodon preceded by the letter corresponding to the corresponding amino-acid. The symbol ’,’ (respectively ’;’) separates two genes that are inside (resp. not inside) the same operon.

The costs and evolutionary scenarios given in Figure 3 are those output by our algorithm after interpretation. In particular, the five Bacillus genomes all show a large inverted segment in the region to the left of the origin of replication (the right part of each linearized representation in Figure 3). As the algorithm without extensions has been used, we preprocessed the genomes by inverting this segment. The genome representations given in Figure 3 are, however, the true ones. Signs given below the red bars represent their true orientations. Consequently, the duplication of the right-most red bar to the left-most position should be interpreted as an inverted duplication that occurred around the origin of replication. On the other hand, some duplications of the red bar have been reported by our algorithm as two separate duplications of two segments separated by a single gene. After careful consideration, these pairs of duplications are more likely a single duplication obscured by subsequent substitution, functional shift or loss of a single tRNA gene. Also, when appropriate (e.g., when a lone gene is positioned in a lone genome), we interpreted some of the losses in our alignments as insertions.

The ancestral genomes given in Figure 3 are those output by our algorithm. They reflect two separate inverted duplications that would have occurred independently in each of the two groups (cereus, anthracis) and (licheniformis, subtilis). We could alternatively infer that the ancestral Bacillus genome already contained both the origin and product of the inverted duplication. The consequence would be the simultaneous loss of the leftmost red bar in three of the five considered genomes. Moreover, nucleotide sequence alignment of the red bars in subtilis and cereus ATCC reveal a higher conservation of pairs of paralogous bars versus orthologous ones, which may indicate that inverted duplications are recent. Whatever the situation is, our inference of a duplication around the origin of replication is in agreement with the observation that has been largely reported in the literature that bacterial genomes have a tendency to preserve a symmetry around the replication origin and terminus (Eisen et al., 2000; Tillier and Collins, 2000; Ajana et al., 2002). The results also show independent proliferation of ribosomal RNA gene-containing operons in Bacillus, which has been associated to selection for increased growth rate (Ardell and Kirsebom, 2005). They also show that in Bacillus, growth-selection on ribosomal RNA operon expansions may significantly alter tRNA gene content as well. The results given in Figure 5 also suggest that some tRNA genes may have been affected by substitutions leading to conversions of function. Such tRNA functional shifts have been detected in metazoan mitochondrial genomes (Rawlings et al., 2003) and bacteria (Saks and Conery, 2007).

OPEN IN VIEWER

FIG. 5.

(a) Genomes of Figure 3 restricted to their red bars; (b) An alignment of all red bars, reflecting one gene insertion (or alternatively 5 gene losses) and substitutions, indicated by ‘*’.

7. Conclusion

We have considered the two species small phylogeny problem for an evolutionary model reduced to non-overlapping (visible) content-modifying operations. We provided some extensions to the model by incorporating overlapping operations, and inversions. Although exponential running times are possible in the worst case, our pseudo-boolean linear programming algorithm turns out to be fast on real datasets, such as the RNA gene repertoire of bacterial genomes. We have also explored avenues for developing efficient non-optimal heuristics. As described in Section 4, a dynamic programming approach can be used to infer a reasonable, though not necessarily optimal unlabeled alignment of two genomes, or a labeled but possibly cyclic alignment. A recent investigation on the complexity of these problems has revealed that the minimum labeling alignment problem is APX-hard [Dondi and El-Mabrouk(2012)]. In future work the apparently easy “core” of finding a possibly cyclic alignment may be leveraged to use the Lagrangian relaxation technique—in this case each constraint against a cyclic duplication set would correspond to a modification of the objective function in our dynamic program.

Application to the RNA gene content of five Bacillus lineages has pointed out a number of interesting biological mechanisms that have to be further investigated. Probably the most interesting observations are the substitutions occurring among the tRNAs. Indeed, as illustrated in Figure 5, some positions indicate a shift of the anticodon, and thus potentially a shift of the original function of the tRNA. Alternatively, such divergence may just be an indication that the tRNA is in the process of losing its function and becoming a pseudogene. Taking into consideration not only the anticodon but rather the entire tRNA sequence, and comparing with additional lineages, may allow to conclude to one or the other alternative. Also other investigations on sequence alignment are required to test our hypothesis that two inverted duplications around the origin of replication have occurred independently on each group (cereus, anthracis) and (licheniformis, subtilis), rather than a single duplication preceding the Bacillus ancestor, followed by losses of the same regions in each of the two monophyletic groups.

8. Alignment Details

Ribosomal RNAs are grouped in bacteria into three families: the 16S, 5S, and 23S rRNAs. As the major role of tRNAs is to serve as adapters between codons along the mRNA and the corresponding amino acids, we group them according to their anticodon. More precisely, the four-letter designation starts with one letter indicating functional class (either an IUPAC one-letter code for a charged amino acid, “X” for initiator or “J” for a special class of isoleucine tRNA) followed by an anticodon sequence in a DNA alphabet. The full alignment as given by our program is available as Supplementary Material at www.liebertpub.com/cmb/.