A probabilistic model for gene content evolution with duplication, loss, and horizontal transfer

We introduce a Markov model for the evolution of a gene family along a phylogeny. The model includes parameters for the rates of horizontal gene transfer, gene duplication, and gene loss, in addition to branch lengths in the phylogeny. The likelihood for the changes in the size of a gene family across different organisms can be calculated in O(N+hM^2) time and O(N+M^2) space, where N is the number of organisms, $h$ is the height of the phylogeny, and M is the sum of family sizes. We apply the model to the evolution of gene content in Preoteobacteria using the gene families in the COG (Clusters of Orthologous Groups) database.


Introduction
At this time, 257 microbial genomes are sequenced, about twice as many are soon to be completed, and 20 complete eukaryotic genomes are publicly available (http://www.ncbi.nlm.nih.gov/Genomes/). These numbers continue to grow in an exponential pace with advances of technology and savvy (Green 2001). The wealth of genome sequence data already caused a revolution in molecular evolution methods (Wolfe and Li 2003;Delsuc et al. 2005). A few years ago, scientific studies have necessarily focused on nucleotide-level differences between orthologous genes, mainly because of the technical and financial limitations on DNA sequence collection. With the increasing amount of whole genome information it becomes possible to analyze genome-scale differences between organisms, and to identify the evolutionary forces responsible for these changes. In particular, sizes of gene families can be compared with the aim of better understanding adaptive evolutionary mechanisms and organismal phylogeny. Several studies point to the fact that gene content may carry sufficient phylogenetic signal for the construction of evolutionary trees (Fitz-Gibbon and House 1999;Snel et al. 1999;Tekaia et al. 1999;Lin and Gerstein 2000;Clarke et al. 2002;Korbel et al. 2002;Dutilh et al. 2004;Gu and Zhang 2004;Huson and Steel 2004;Lake and Rivera 2004). Comparative analyses of genome-wide protein domain content (Lin and Gerstein 2000;Yang et al. 2005;Deeds et al. 2005) have also provided important insights into evolution. Gene content and similar features have been used to construct viral (Montague and Hutchison III 2000;Herniou et al. 2001), microbial (Snel et al. 1999;Fitz-Gibbon and House 1999;Gu and Zhang 2004), and universal trees (Tekaia et al. 1999;Simonson et al. 2005;Yang et al. 2005). Comparative gene content analysis is also used to estimate ancestral genome composition Mirkin et al. 2003). The presence-absence pattern of homologs in different organisms, the so-called phyletic pattern (Koonin and Galperin 2002;Tatusov et al. 2003), provides clues about gene function (Pellegrini et al. 1999) and evolution of metabolic pathways (Mirkin et al. 2003).
There are two principal methodological issues in the analysis of gene content. First, one needs to decide how homologous gene are selected, and secondly, an appropriate computational technique must be chosen to analyze the data. A possible choice for compiling a data set is to use pairwise orthologs and compute a score for each pair of genomes (Snel et al. 1999;Tekaia et al. 1999;Korbel et al. 2002;Clarke et al. 2002;Dutilh et al. 2004). The matrix of pairwise scores of shared gene content is then amenable to analysis by distance-based methods of phylogeny construction. An alternative is proposed in (Lake and Rivera 2004): for each pair of genomes, the presence or absence of homologs with respect to a third reference genome is noted, and the pairs of presence-absence sequences are used to calculate a pairwise distance matrix. A third approach is to first compute families of homologous genes, and then to record the absence or presence of each family in the genomes (Fitz-Gibbon and House 1999;Lin and Gerstein 2000;Jordan et al. 2001;Wolf et al. 2001). The resulting absence-presence data can be treated as a set of 0-1 sequences and further analyzed with traditional parsimony or distance-based methods. Some specialized parsimony methods have been developed for the direct purpose of analyzing gene absencepresence data (Mirkin et al. 2003;Kunin and Ouzounis 2003). Instead of simply noting presence-absence, the gene family sizes give an even richer signal for evolutionary analyses Huson and Steel 2004;Hahn et al. 2005). It is that latter type of data that we model here.
A number of processes shape the gene content of an organism. New genes may be created by duplication of an existing gene, horizontal transfer from a different lineage, and fusion/fission . It has been widely debated how the extent of horizontal gene transfer (HGT) compares to vertical inheritance (Jordan et al. 2001;Snel et al. 2002;Gogarten et al. 2002;Kurland et al. 2003;Kunin et al. 2005;Ge et al. 2005;Simonson et al. 2005). It is clear that horizontal gene transfer plays a major role in microbial evolution (Boucher et al. 2003), but there is still need for adequate mathematical models in which that role can be measured.
This paper describes a probabilistic model for gene content evolution. Specifically, we model the evolution of gene families along a phylogeny. The model includes gene duplication, gene loss, and horizontal transfer as mechanisms that determine gene family evolution. We also show how to compute exact likelihoods for gene family sizes in different organisms. A few probabilistic models were proposed for gene content evolution, but they are less general than ours. Usual stochastic models work with two parameters. Gu and Zhang (Gu and Zhang 2004;Gu et al. 2005) rely on a model that includes gene loss and gene duplication but no other modes of gene genesis. They showed how gene family sizes can be used to define additive distances in such a model. Interestingly enough, the data can be reduced to a three-letter alphabet for the purposes of distance calculations: only 0, 1 or "many" homologs per family need to be counted. The distance relies on an estimate of the rate parameters, which is obtained through likeli-hood optimization. Hahn et al. (Hahn et al. 2005) developed an alternative likelihood-based approach for the same two-parameter model. Huson and Steel (Huson and Steel 2004) analyzed a two-parameter model that accounts for gene loss and horizontal transfer but not for gene duplication. They derived a distance measure based on gene family sizes using likelihood maximization arguments. They further showed that other distance measures based on shared gene content (Snel et al. 1999) have inferior accuracy in phylogeny reconstruction than either Dollo parsimony or their own distance measure. Karev et al. developed a rich probabilistic model of gene content evolution in a series of papers (Karev et al. 2002;Karev et al. 2003;Karev et al. 2004). The model explains the distribution of gene family sizes found in different organisms. It is, however, too general for exact detailed calculations, and for likelihood computations in particular.
To our knowledge, no tractable stochastic model was introduced yet that accounts for horizontal transfer, gene loss, and duplication. These processes cannot be modeled by using only two parameters because the intensity of gene loss and duplication depend on the size of a gene family, but the tempo by which genes are acquired through horizontal transfer has a constant component. Among other applications, a model that accounts for duplication and transfer is useful in analyzing the evolution of metabolic networks: do new paths evolve by gene duplication and adaptive selection, or by accommodating genes with new functions via horizontal gene transfer? This paper introduces a probabilistic model for the evolution of homologs on a phylogeny. Specifically, we model the evolution of a single gene family along the phylogeny, where different processes may add new genes to the family or erase members of it, and arrive at the family sizes observed at the terminal taxa. We provide an algorithm that can compute analytically the likelihood of gene family sizes in different organisms, given an evolutionary tree. The algorithm calculates the likelihood of family sizes in O(N + M 2 h) time where M is the total number of genes in the family, N is the number of genomes, and h is the height of the tree. The tree height is at most linear in N, and on average, it is O( √ N ) or O(log N) for uniform or Yule-Harding distribution of random trees. In contrast, the methods of (Gu and Zhang 2004) and (Huson and Steel 2004) compute distances between every pair of organisms, which takes quadratic time in N. The likelihood calculations of (Hahn et al. 2005) take cubic time in M, and involve the evaluation of infinite sums that are truncated heuristically.
The article is organized in the following manner. Section 2 introduces our stochastic model of gene content evolution, and describes formulas for computing various associated probabilities, including likelihood. The formulas are used in an algorithm described in Section 3. Section 4 describes our initial experiments in modeling gene content evolution in 51 Proteobacteria and 3555 gene families from the COG (Clusters of Orthologous Groups) database (Tatusov et al. 2003). Section 5 concludes the paper.

Mathematical model
Let T be a phylogenetic tree over a set of species S. The tree T is a rooted tree with vertex set V (T ) and edge set E(T ), in which leaves are bijectively labeled with elements of S. Every edge e has a length t e > 0. We are interested in modeling the evolution of a gene family. The family size changes along the edges: genes may be duplicated, lost, or gained from an unknown source. We model the evolution of gene counts (family size) at the tree nodes: the gene count at every node u ∈ V (T ) is a random variable χ(u) that can take non-negative integer values. In addition to the tree with its edge lengths, three parameters determine the joint distribution of the gene counts: a duplication rate λ, a loss rate µ, and a transfer rate κ. The loss rate accounts for all possible mechanisms of gene loss, including deletion and pseudogenization. The transfer rate accounts for processes of gene genesis, including HGT from another lineage in the same tree, or HGT from an unknown organism. In our model, the evolution of the gene counts on a branch follows a linear birth-and-death process (Feller 1950) parametrized by λ, κ, and µ. Let {X(t) : t ≥ 0} denote the continuous-time Markov process formed by the gene counts along an edge uv: χ(u) = X(0) and χ(v) = X(t uv ). The transition probabilities of the process are the following: In other words, every existing gene produces an offspring through duplication with an intensity of λ, or disappears with an intensity of µ, and new genes are Figure 1: Galton-Watson forest that shows the evolution of genes in the same family along a tree edge. The top line represents the ancestral genome with three genes; the bottom line represents the descendant genome, in which there are five family members. Symbol o represents the source from which genes might be transferred horizontally, symbols ⋆ represent copies of the gene in the genome at the beginning and the end of the investigated time span t. Each o or ⋆ in the ancestral genome is the root of a Galton-Watson tree. Note that the physical order of genes in the genomes is immaterial: here they are simply drawn next to each other for clarity.
acquired with an intensity of κ, independently from the number of existing genes.
The histories of individual genes on an edge form a Galton-Watson forest, see Figure 1. The figure illustrates a scenario where the gene family increases from three to five genes. The counts at the branch endpoints are the result of many duplication, transfer and loss events. The change involves three horizontally transferred genes, from among which one survives, another one does not, and the third one produces two surviving paralogs.
While it is not too difficult to calculate the probabilities for any particular gene count on a branch (see §2.1), the likelihood L of observed gene counts at the leaves involves an infinite number of possible gene counts at intermediate nodes: where γ(·) is the distribution at the root, and the summation over the m x vectors takes all values in agreement with the gene counts at the leaves in the input data. Our main technique for computing the likelihood is to restrict the computation to genes that have at least one surviving descendant at the leaves. In what follows we develop the formulas to compute the likelihood.

Basic transition probabilities
First we analyze the blocks of homologs at a node that have a common origin. One block is formed by the genes that trace back to a horizontal transfer event on the branch from the parent. Each other block is the set of paralogs with the same ancestor at the parent. The homologs in Figure 1 form four blocks: a block of size three that comprises the descendants of the horizontally transferred genes, a block of size zero for the deceased parental gene, and two blocks of size one. The independent birth-and-death processes associated with the blocks have been thoroughly analyzed in the statistical literature.
Definition 1. Define the following basic transition probabilities for gene count evolution on a branch. Let h t (n) denote the probability that there are n genes of foreign origin [not inherited from the parent] after time t. Let g t (n) denote the probability that a single gene has n copies after time t.
Theorem 1. The basic transition probabilities can be written as follows.
Proof. The size of the HGT block of homologs follows a birth-and-death process with constant rate κ of immigration and no emigration. The transition probabilities of (2) for such a process were analyzed by Karlin and McGregor (Karlin and McGregor 1958). Blocks of paralogs evolve by a simple birth-and-death process: the transition probabilities of (3) are derived in, e.g., (Feller 1950).

Gene extinction and survival
Definition 2. A surviving gene at a node x is such that it has at least one modern descendant at the leaves below x.
Let D x denote the probability that a gene present at node x is not surviving, i.e., that it has no modern descendants.
Lemma 1. The extinction probability D x can be calculated as follows. If x is a leaf, then D x = 0. Otherwise, let x be the parent of x 1 , x 2 , . . . , x k .
where t j is the length of the branch leading from x to x j .
Proof. For leaves, the statement is trivial. When x is not a leaf, condition on the gene count at x j : Plugging in g t (m) from Eq. (3) and replacing the infinite series with a closed form gives (4).

Effective transition probabilities
We introduce two new probabilities, denoted by H x (n) and G x (n). They account for the number n of surviving genes at node x, either acquired through horizontal transfer, or through duplication and loss from a single gene. The effective transition probabilities are related to h t (n), and g t (n), but account for eventual extinction below node x. A formal definition follows.
Definition 3. Let y be a non-root node and x its parent. Define the following effective transition probabilities. Let H y (n) denote the probability that at node y there are n surviving genes of foreign origin, i.e., that have no ancestor at x. Let G y (n) denote the probability that a single gene at x has n surviving copies at node y.
Lemma 2. Let y be a non-root node, let x be its ancestor, and let t be the length of the edge xy. The effective transition probabilities can be written as follows.
Proof. We condition on the number of genes at node y (whether or not they survive).
Using Eq. (2) leads to an infinite series that can be simplified to get (5). Similarly, write Taking the values of g t (n + i) from Eq. (3) and simplifying the resulting infinite series yields (6).

Number of surviving genes on a branch
Definition 4. Let y be a non-root node, and let x be its ancestor. Let p y (m|n) denote the survival probability defined as the probability of the event that there are m surviving genes at node y under the condition that there are n genes at node x (not necessarily surviving).
Lemma 3. The survival probabilities can be computed as follows.
Proof. For p y (m|0) and p y (0|n), the equations are straightforward. Otherwise, we condition on the surviving copies of a single gene at y: Now, using that G y (i + 1) = αG y (i) whenever i > 0, and comparing (9) for p y (m|n) and p y (m − 1|n), we can write p y (m|n) in a recursive form as shown.

Conditional likelihoods
Definition 5. Let x be a node in the tree. Define the conditional likelihood L x (n) for all n as the probability of having the observed gene counts at the leaves in the subtree rooted at x, under the condition that there are n surviving copies at x.
Theorem 2. The conditional likelihoods can be calculated as follows. In the case when x is a leaf, L x (n) = 1 if n is the observed gene count at x, otherwise the likelihood is 0. If x is not a leaf, and has children x 1 , x 2 , . . . , x k , then the following recursions hold.
where M j is the sum of gene counts at the leaves in the subtree rooted at x j . If n > k j=1 M j , then L x (n) = 0.
Proof. For a leaf node, or for n > k j=1 M j , the theorem is trivial. Otherwise, consider the likelihood ℓ x (n) of the observed gene counts at the leaves in the subtree rooted at x, conditioned on the event that there are n genes present at x, which may or may not survive. We write the likelihood in two ways. First, by conditioning on the number of surviving genes at the children, Secondly, by conditioning on the number of surviving genes at x, Now, rearranging the equality of the two right-hand sides gives the desired result.
Remark. Clearly, the gene counts M x of Theorem 2 are easily computed for all x. If m(x) is the gene count for every leaf x then . , x k are the children of x. (13)

Likelihood
It is assumed that the family size at the root is distributed according to the equilibrium probabilities: Theorem 3. Let M be the total number of genes at the leaves. The likelihood of the observed gene counts equals Proof. By summing the likelihoods conditioned on the surviving genes at the root, Now, plugging in the values of γ(·) from Eq. (14) and replacing the infinite series by a closed form gives the theorem's formula.

Algorithm
This section employs the formulas of Section 2 in a dynamic programming algorithm to compute the likelihood exactly. More precisely, the algorithm computes the likelihood of gene counts at the tree leaves, given the duplication rate λ, the transfer rate κ, and the loss rate µ. Algorithm Com-puteLikelihood below proceeds by a depth-first traversal; the necessary variables are calculated from the leaves towards the root. Let m(u) denote the gene count at every leaf u. 9 Compute the likelihood L at the root using Eq. (15).
Theorem 4 below analyzes the algorithm's complexity in terms of the topology of T . In particular, it uses the notions of height of a node x, defined as the number of edges on the path leading from the root to x, levels of nodes, which are sets of nodes with the same height, and height of the tree, which is the maximum of the leaf heights. Clearly, x M x ≤ M if the summation goes over x for which their subtrees do not overlap, such as nodes at the same level. Now, x M 2 x ≤ ( x M x ) 2 ≤ M 2 , and, thus, O(M 2 ) time is spent on each level. Therefore, the total time spent in the loop of Line 4 is ·O(hM 2 ). Line 9 takes O(M) time.
In order to obtain the space complexity result, notice that at the end of the loop in Line 8 the computed variables for the children of x are not needed anymore. Therefore, the nodes for which p x (·|·) is needed are such that their subtrees do not overlap. By the same type of argument as with time spent on a level, the number of variables that need to be kept in memory is O(M 2 ).

Gene content evolution in Proteobacteria
Proteobacteria form one of the most diverse groups of Prokaryotes. Proteobacteria are an excellent model case for studying genome content evolution: they include pathogens, endosymbionts, and free-living organisms. Genome sizes vary tenfold within this group, and horizontal transfer is abundant (Gogarten et al. 2002). Their phylogeny is still not resolved to satisfaction (Lerat et al. 2003;Boussau et al. 2004;Herbeck et al. 2005;Belda et al. 2005). We used 51 Proteobacteria in the first application of our likelihood method. Gene counts were based on the newer version (Tatusov et al. 2003) of the COG database. Each COG is a manually curated protein family of homologs. The COGs are classified into 23 functional categories. For each of our 51 Protebacteria, the number of genes in each COG family was established. There are 3555 COG families that have at least one member in the organisms. (The organisms are listed in the Appendix.) The data set was provided to us by Csaba Pál and Martin Lercher (Pál et al. 2005). The purpose of applying the likelihood method was not to carry out in-depth data analysis, but rather to get a first impression of our method's performance on realistic data.
First we optimized the branch lengths and the λ, κ parameters while keeping µ = 1.0 to fix the scaling. In a second pass, we clustered the COG  17 2 5 18 19 11 7 22 31 27 19 15 23 8 27 9 98 220  583   1 10 8 1 8 13 15 7 14 8 30 28 13 7 2 8 28 18 79 117  405   1 14 22 2 0 18 6 3 13 7 12 22 21 1 2 4 10 13 55 61  274   3 18 23 2 6 13 13 4 1 8 22 30 14 3 3 9 9 18 81  modification, protein turnover and chaperones, C-energy production and conversion, G-carbohydrate transport and metabolism, E-amino acid transport and metabolism, F-nucleotide transport and metabolism, H-coenzyme transport and metabolism, I-lipid transport and metabolism, P-inorganic ion transport and metabolism, Q-secondary metabolites biosynthesis, transport and catabolism, R-general function prediction only, S-function unknown. The "size" columns gives the number of COGs in each rate group. (The numbers in one row do not always add up to the value in the "size" column because some COGs have more than one functional assignment.) families with different rates in different groups. The groups were established in several iterations of Expectation Maximization: in an E-step, each family was assigned to the best group (the one whose rates give the highest likelihood), in an M-step, rates were optimized within each group separately to maximize the likelihood of the COG gene counts within the group's families. Figure 2 shows the rates in different groups (Groups 0-8), as well as the distribution of COG functional classes across clusters. The picture shows that various rate groups are needed to describe the evolution of the families. While the results and the methodology still need a thorough critical assessment, some interesting patterns already emerge. About 19% of the families are very stable (Group 8), this includes the large majority of genes involved in translation (category J) such as tRNA synthetases and ribosomal proteins, and cell cycle control (category D). About one in nine families fall into the groups with large horizontal transfer rates (Groups 1 and 7), while one in three families are in groups with very low transfer rates. There are many categories where duplication plays only a minor role. For instance, the evolution of cell motility (category N), and various metabolic functions (F,H,I) seem to be shaped mainly by horizontal transfer and loss.

Conclusion
We presented the first three-parameter model of gene content evolution, along with a fast algorithm for computing likelihoods. We implemented parameter optimization and a gene family clustering method and carried out a pilot experiment using COG family sizes in 51 Proteobacteria.
We modeled gene family evolution by a birth-and-death process. It was shown that birth-and-death processes (as opposed to concerted evolution) appropriately represent family evolution in at least some gene families (Nei et al. 1997;Michelmore and Meyers 1998;Piontkivska et al. 2002). In addition, birth-and-death processes of various complexity explain the observed power-law behavior of gene family sizes (Karev et al. 2002;Karev et al. 2003;Karev et al. 2004;Reed and Hughes 2004). In order to develop a truly realistic likelihood model, rate variation must be permitted across lineages and families. Our formulas can be readily adapted to branch-dependent rates. The challenge lies rather in the parametrization: introducing four parameters (three rates and branch length) for every tree edge and every family will lead to overfitting. A possible solution is to work with parameters that depend only on the gene family and parameters that depend only on the branch. These two sets of parameters are combined for each branch and family to infer branch-specific rates. We are now working on developing adequate rate-variation models, in a clustering-based approach as we did in Section 4, and by imposing rate distribution functions. In another line of extension, we are investigating the coupling of the model with sequence evolution models to enable a finer modeling of homologies than simple counts. By scoring the similarity of genes within the same family, one can arrive to a finer likelihood model of gene content evolution.
It is interesting to point out that while the mathematical model assigns a non-zero probability to the case when the gene family has no members at any of the leaves, a family with no extant genes is not included usually in the data. Consequently, likelihood methods tend to underestimate the extent of gene losses. The situation is similar to what is encountered in likelihood models of intron evolution and a possible remedy is discussed in (Csűrös 2005).
This paper focuses on the core algorithmic problems of likelihood computations in a biologically realistic model of gene content evolution. The presented likelihood algorithm can be utilized in a number of contexts. The computations can be used in parameter optimization for estimating duplication, loss, and transfer rates in different gene families. By comparing the maximum likelihood values achieved with different evolutionary tree topologies, organismal phylogeny can be derived based on gene content. "Unusual" branches with excess transfer, loss, etc., can be identified by examining the likelihoods, adapting an idea of (Hahn et al. 2005). The conditional likelihoods of §2.5 can be used in likelihood-based computations of ancestral gene content, similarly to standard methods employed in case of molecular sequences (Pupko et al. 2000). The likelihood computation enables also the sampling of different trees in a Bayesian Markov Chain Monte Carlo approach. We believe that our approach to computing exact likelihoods efficiently in the three-parameter model will find many applications in comparative gene content analysis.