Multispecies coalescent process

In the basic multispecies coalescent model, the species phylogeny is assumed to be known. Complete isolation after species divergence, with no migration, hybridization, or introgression is also assumed. We assume no recombination so that all the sites within the locus share the same gene tree (topology and coalescent times). However, the basic model can be extended in different ways to accommodate migration or introgression, population size changes, recombination.[13][14]

The model and implementation of this method can be applied to any species tree. As an example, the species tree of the great apes: humans (H), chimpanzees (C), gorillas (G) and orangutans (O) is considered. The topology of the species tree, (((HC)G)O)), is assumed known and fixed in the analysis (Figure 1).[1] Let ${\displaystyle D=\{D_{i}\}}$ be the entire data set, where ${\displaystyle {D_{i}}}$ represent the sequence alignment at locus ${\displaystyle i}$, with ${\displaystyle i=1,2,\ldots ,L}$ for a total of ${\displaystyle L}$ loci.

The population size of a current species is considered only if more than one individual is sampled from that species at some loci.

The parameters in the model for the example of Figure 1 include the three divergence times ${\displaystyle \tau _{HC}}$, ${\displaystyle \tau _{HCG}}$ and ${\displaystyle \tau _{HCGO}}$ and population size parameters ${\displaystyle \theta _{H}}$ for humans; ${\displaystyle \theta _{C}}$ for chimpanzees; and ${\displaystyle \theta _{HC}}$, ${\displaystyle \theta _{HCG}}$ and ${\displaystyle \theta _{HCGO}}$ for the three ancestral species.

The divergence times (${\displaystyle \tau }$'s) are measured by the expected number of mutations per site from the ancestral node in the species tree to the present time (Figure 1 of Rannala and Yang, 2003).

Therefore, the parameters are ${\displaystyle \Theta =\{\theta _{H},\theta _{C},\theta _{HC},\theta _{HCG},\theta _{HCGO},\tau _{HC},\tau _{HCG},\tau _{HCGO}\}}$.

The joint distribution of ${\displaystyle f(T_{i},t_{i}\mid \Theta )}$ is derived directly in this section.[1] Two sequences from different species can coalesce only in one populations that are ancestral to the two species. For example, sequences H and G can coalesce in populations HCG or HCGO, but not in populations H or HC. The coalescent processes in different populations are different.

For each population, the genealogy is traced backward in time, until the end of the population at time ${\displaystyle \tau }$, and the number of lineages ${\displaystyle (m)}$ entering the population and the number of lineages leaving it ${\displaystyle (n)}$ are recorded. For example, ${\displaystyle m=3,n=2,}$ and ${\displaystyle \tau =\tau _{HC}}$, for population H (Table 1).[1] This process is called a censored coalescent process because the coalescent process for one population may be terminated before all lineages that entered the population have coalesced. If ${\displaystyle n\geq 1}$ the population consists of ${\displaystyle n}$ disconnected subtrees or lineages.

With one time unit defined as the time taken to accumulate one mutation per site, any two lineages coalesce at the rate ${\displaystyle {\frac {2}{\theta }}}$. The waiting time ${\displaystyle t_{j}}$ until the next coalescent event, which reduces the number of lineages from ${\displaystyle j}$ to ${\displaystyle j-1}$ has exponential density

${\displaystyle f(t_{j})={\frac {j(j-1)}{2}}{\frac {2}{\theta }}\exp\{-{\frac {j(j-1)}{2}}{\frac {2}{\theta }}t_{j}\},\quad j=m,m-1,\ldots ,n+1}$

If ${\displaystyle n\geq 1}$, the probability that no coalescent event occurs between the last one and the end of the population at time ${\displaystyle \tau }$; i.e. during the time interval ${\displaystyle \tau -(t_{m}+t_{m-1}+\ldots +t_{n+1})}$. This probability is ${\displaystyle \exp\{-{\frac {n(n-1)}{\theta }}[\tau -(t_{m}+t_{m-1}+\ldots +t_{n+1})]}$ and is 1 if ${\displaystyle n=1}$.

(Note: One should recall that the probability of no events over time interval ${\displaystyle t}$ for a Poisson process with rate ${\displaystyle \lambda }$ is ${\displaystyle e^{-\lambda t}}$. Here the coalescent rate when there are ${\displaystyle n}$ lineages is ${\displaystyle \lambda ={\frac {n(n-1)}{\theta }}}$.)

In addition, to derive the probability of a particular gene tree topology in the population, if a coalescent event occurs in a sample of ${\displaystyle j}$ lineages, the probability that a particular pair of lineages coalesce is ${\displaystyle 1/{\binom {j}{2}}=2/j(j-1),\quad j=m,m-1,\ldots ,n+1}$.

Multiplying these probabilities together, the joint probability distribution of the gene tree topology in the population and its coalescent times ${\displaystyle t_{m},t_{m+1},\ldots ,t_{n+1}}$ as

${\displaystyle \prod _{j=n+1}^{m}{\Big [}{\frac {2}{\theta }}\exp {\Big \{}-{\frac {j(j-1)}{\theta }}t_{j}{\Big \}}{\Big ]}\exp {\Big \{}-{\frac {n(n-1)}{\theta }}(\tau -(t_{m}+t_{m+1}+\ldots +t_{n+1})){\Big \}}}$.

The probability of the gene tree and coalescent times for the locus is the product of such probabilities across all the populations. Therefore, the gene genealogy of Figure 1,[1][15] we have

${\displaystyle {\begin{aligned}f(G_{i}\mid \Theta )&=[2/\theta _{H}\exp\{-6t_{3}^{(H)}/\theta _{H}\}\exp\{-2(\tau _{HC}-t_{3}^{(H)})/\theta _{H}\}]\\&{}\times [2/\theta _{C}\exp\{-2t_{2}^{(C)}/\theta _{C}\}]\\&{}\times [2/\theta _{HC}\exp\{-6t_{3}^{HC}/\theta _{HC}\}]\times [2/\theta _{HC}\exp\{-2t_{2}^{HC}/\theta _{HC}\}]\\&{}\times [\exp\{-2(\tau _{HCG}-\tau _{HG}-(t_{3}^{HC}+t_{2}^{HC}))/\theta _{HCG}\}]\\&{}\times [2/\theta _{HCGO}\exp\{-6t_{3}^{HCGO}/\theta _{HCGO}\}]\times [2/\theta _{HCGO}\exp\{-2t_{2}^{HCGO}/\theta _{HCGO}\}]\end{aligned}}}$

The gene genealogy ${\displaystyle G_{i}}$ at each locus ${\displaystyle i}$ is represented by the tree topology ${\displaystyle T_{i}}$ and the coalescent times ${\displaystyle t_{i}}$. Given the species tree and the parameters ${\displaystyle \Theta }$ on it, the probability distribution of ${\displaystyle G_{i}=\{T_{i},t_{i}\}}$ is specified by the coalescent process as

${\displaystyle f(G\mid \Theta )=\prod _{i}f(G_{i}\mid \Theta )=\prod _{i}f(T_{i},t_{i}\mid \Theta )}$,

where ${\displaystyle f(G_{i}\mid \Theta )=f(T_{i},t_{i}\mid \Theta )}$ is the probability density for the gene tree at locus locus ${\displaystyle i}$,[1] and the product is because we assume that the gene trees are independent given the parameters.

The probability of data ${\displaystyle D_{i}}$ given the gene tree and coalescent times (and thus branch lengths) at the locus, ${\displaystyle f(D_{i}\mid G_{i})}$, is Felsenstein's phylogenetic likelihood.[16] Due to the assumption of independent evolution across the loci,

${\displaystyle f(D\mid G)=\prod _{i}f(D_{i}\mid G_{i})}$

The likelihood function or the probability of the sequence data given the parameters ${\displaystyle \Theta }$ is then an average over the unobserved gene trees

${\displaystyle f(D\mid \Theta )=\int f(D\mid G)f(G\mid \Theta )dG,}$

where the integration represents summation over all possible gene tree topologies (${\displaystyle T_{i}}$) and, for each possible topology at each locus, integration over the coalescent times ${\displaystyle t_{i}}$.[17] This is in general intractable except for very small species trees.

In Bayesian inference, we assign a prior on the parameters, ${\displaystyle f(\Theta )}$, and then the posterior is given as

${\displaystyle f(\Theta \mid D)=\int f(\Theta ,G\mid D)dG,}$

where again the integration represents summation over all possible gene tree topologies (${\displaystyle T_{i}}$) and integration over the coalescent times ${\displaystyle t_{i}}$. In practice this integration over the gene trees is achieved through a Markov chain Monte Carlo algorithm, which samples from the joint conditional distribution of the parameters and the gene trees

${\displaystyle f(\Theta ,G\mid D)\propto f(D\mid G)f(G\mid \Theta )f(\Theta ).}$

The above assumes that the species tree is fixed. In species-tree estimation, the species tree (${\displaystyle S}$) changes as well, so that the joint conditional distribution (from which the MCMC samples) is

${\displaystyle f(S,\Theta ,G\mid D)\propto f(D\mid G)f(G\mid S,\Theta )f(\Theta )f(S),}$

where ${\displaystyle f(S)}$ is the prior on species trees.

As a major departure from two-step summary methods, full-likelihood methods average over the gene trees. This means that they make use of information in the branch lengths (coalescent times) on the gene trees and accommodate their uncertainties (due to limited sequence length in the alignments) at the same time. It also explains why full-likelihood methods are computationally much more demanding than two-step summary methods.

The integration or summation over the gene trees in the definition of the likelihood function above is virtually impossible to compute except for very small species trees with only two or three species.[18] Full-likelihood or full-data methods, based on calculation of the likelihood function on sequence alignments, have thus mostly relied on Markov chain Monte Carlo algorithms. MCMC algorithms under the multispecies coalescent model are similar to those used in Bayesian phylogenetics but are distinctly more complex, mainly due to the fact that the gene trees at multiple loci and the species tree have to be compatible: sequence divergence has to be older than species divergence. As a result, changing the species tree while the gene trees are fixed (or changing a gene tree while the species tree is fixed) leads to inefficient algorithms with poor mixing properties. Considerable efforts have been taken to design smart algorithms that change the species tree and gene trees in a coordinated manner, as in the rubber-band algorithm for changing species divergence times,[1] the coordinated NNI, SPR and NodeSlider moves.[19][20]

Consider for example the case of two species (A and B) and two sequences at each locus, with a sequence divergence time ${\displaystyle t_{i}}$ at locus ${\displaystyle i}$. We have ${\displaystyle t_{i}<\tau }$ for all ${\displaystyle i}$. When we want to change the species divergence time ${\displaystyle \tau }$ within the constraint of the current ${\displaystyle t_{i}}$, we may have very little room for change, as ${\displaystyle \tau }$ may be virtually identical to the smallest of the ${\displaystyle t_{i}}$. The rubber-band algorithm[1] changes ${\displaystyle \tau }$ without consideration of the ${\displaystyle t_{i}}$, and then modifies the ${\displaystyle t_{i}}$ deterministically in the same way that marks on a rubber band move when the rubber band is held from a fixed point pulled towards one end. In general, the rubber-band move guarantees that the ages of nodes in the gene trees are modified so that they remain compatible with the modified species divergence time.

Full likelihood methods tend to reach their limit when the data consist of a few hundred loci, even though more than 10,000 loci have been analyzed in a few published studies.[21][22]