Multispecies Coalescent Process is a stochastic process model that describes the genealogical relationships for a sample of DNA sequences taken from several species.[1] [2] It represents the application of coalescent theory to the case of multiple species. The multispecies coalescent results in cases where the relationships among species for an individual gene (the gene tree) can differ from the broader history of the species (the species tree). It has important implications for the theory and practice of phylogenetics and for understanding genome evolution.
A gene tree is a binary graph that describes the evolutionary relationships between a sample of sequences for a non-recombining locus. A species tree describes the evolutionary relationships between a set of species, assuming tree-like evolution. However, several processes can lead to discordance between gene trees and species trees. The Multispecies Coalescent model provides a framework for inferring species phylogenies while accounting for ancestral polymorphism and gene tree-species tree conflict. The process is also called the Censored Coalescent.[1]
Besides species tree estimation, the multispecies coalescent model also provides a framework for using genomic data to address a number of biological problems, such as estimation of species divergence times, population sizes of ancestral species, species delimitation, and inference of cross-species gene flow.[3]
If we consider a rooted three-taxon tree, the simplest non-trivial phylogenetic tree, there are three different tree topologies but four possible gene trees.[4] The existence of four distinct gene trees despite the smaller number of topologies reflects the fact that there are topologically identical gene tree that differ in their coalescent times. In the type 1 tree the alleles in species A and B coalesce after the speciation event that separated the A-B lineage from the C lineage. In the type 2 tree the alleles in species A and B coalesce before the speciation event that separated the A-B lineage from the C lineage (in other words, the type 2 tree is a deep coalescence tree). The type 1 and type 2 gene trees are both congruent with the species tree. The other two gene trees differ from the species tree; the two discordant gene trees are also deep coalescence trees.
The distribution of times to coalescence is actually continuous for all of these trees. In other words, the exact coalescent time for any two loci with the same gene tree may differ. However, it is convenient to break up the trees based on whether the coalescence occurred before or after the earliest speciation event.
Given the internal branch length in coalescent units it is straightforward to calculate the probability of each gene tree.[5] For diploid organisms the branch length in coalescent units is the number of generations between the speciation events divided by twice the effective population size. Since all three of the deep coalescence tree are equiprobable and two of those deep coalescence tree are discordant it is easy to see that the probability that a rooted three-taxon gene tree will be congruent with the species tree is:
\begin{align} P(congruence)&=1-
2 | |
3 |
\exp(-T)=1-
2 | \exp(- | |
3 |
t | |
2Ne |
) \end{align}
The phenomenon of hemiplasy is a natural extension of the basic idea underlying gene tree-species tree discordance. If we consider the distribution of some character that disagrees with the species tree it might reflect homoplasy (multiple independent origins of the character or a single origin followed by multiple losses) or it could reflect hemiplasy (a single origin of the trait that is associated with a gene tree that disagrees with the species tree).
The phenomenon called incomplete lineage sorting (often abbreviated ILS in the scientific literatures[7]) is linked to the phenomenon. If we examine the illustration of hemiplasy with using a rooted four-taxon tree (see image to the right) the lineage between the common ancestor of taxa A, B, and C and the common ancestor of taxa A and B must be polymorphic for the allele with the derived trait (e.g., a transposable element insertion[8]) and the allele with the ancestral trait. The concept of incomplete lineage sorting ultimately reflects on persistence of polymorphisms across one or more speciation events.
The probability density of the gene trees under the multispecies coalescent model is discussed along with its use for parameter estimation using multi-locus sequence data.
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.[9] [10]
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,, is assumed known and fixed in the analysis (Figure 1). Let
D=\{Di\}
{Di}
i
i=1,2,\ldots,L
L
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
\tauHC
\tauHCG
\tauHCGO
\thetaH
\thetaC
\thetaHC
\thetaHCG
\thetaHCGO
The divergence times (
\tau
Therefore, the parameters are
\Theta=\{\thetaH,\thetaC,\thetaHC,\thetaHCG,\thetaHCGO,\tauHC,\tauHCG,\tauHCGO\}
The joint distribution of
f(Ti,ti\mid\Theta)
For each population, the genealogy is traced backward in time, until the end of the population at time
\tau
(m)
(n)
m=3,n=2,
\tau=\tauHC
n\geq1
n
With one time unit defined as the time taken to accumulate one mutation per site, any two lineages coalesce at the rate
2 | |
\theta |
tj
j
j-1
f(t | ||||
|
2 | \exp\{- | |
\theta |
j(j-1) | |
2 |
2 | |
\theta |
tj\}, j=m,m-1,\ldots,n+1
If
n\geq1
\tau
\tau-(tm+tm-1+\ldots+tn+1)
\exp\{- | n(n-1) |
\theta |
[\tau-(tm+tm-1+\ldots+tn+1)]
n=1
(Note: One should recall that the probability of no events over time interval
t
λ
e-λ
n
λ=
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
j
1/\binom{j}{2}=2/j(j-1), 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
tm,tm+1,\ldots,tn+1
m | ||
\prod | [ | |
j=n+1 |
2 | \exp\{- | |
\theta |
j(j-1) | |
\theta |
tj\}]\exp\{-
n(n-1) | |
\theta |
(\tau-(tm+tm+1+\ldots+tn+1))\}
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,[11] we have
\begin{align} f(Gi\mid\Theta)&
(H) | |
=[2/\theta | |
3/\theta |
H\}\exp\{-2(\tauHC
(H) | |
-t | |
3)/\theta |
H\}]\\ &{} x [2/\thetaC
(C) | |
\exp\{-2t | |
2/\theta |
C\}]\\ &{} x [2/\thetaHC
HC | |
\exp\{-6t | |
3/\theta |
HC\}] x [2/\thetaHC
HC | |
\exp\{-2t | |
2/\theta |
HC\}]\\ &{} x [\exp\{-2(\tauHCG-\tauHG
HC | |
-(t | |
2))/\theta |
HCG\}]\\ &{} x [2/\thetaHCGO
HCGO | |
\exp\{-6t | |
3/\theta |
HCGO\}] x [2/\thetaHCGO
HCGO | |
\exp\{-2t | |
2/\theta |
HCGO\}] \end{align}
The gene genealogy
Gi
i
Ti
ti
\Theta
Gi=\{Ti,ti\}
f(G\mid\Theta)=\prodif(Gi\mid\Theta)=\prodif(Ti,ti\mid\Theta)
where
f(Gi\mid\Theta)=f(Ti,ti\mid\Theta)
i
The probability of data
Di
f(Di\midGi)
f(D\midG)=\prodif(Di\midGi)
The likelihood function or the probability of the sequence data given the parameters
\Theta
f(D\mid\Theta)=\intf(D\midG)f(G\mid\Theta)dG,
where the integration represents summation over all possible gene tree topologies (
Ti
ti
In Bayesian inference, we assign a prior on the parameters,
f(\Theta)
f(\Theta\midD)=\intf(\Theta,G\midD)dG,
where again the integration represents summation over all possible gene tree topologies (
Ti
ti
f(\Theta,G\midD)\proptof(D\midG)f(G\mid\Theta)f(\Theta).
The above assumes that the species tree is fixed. In species-tree estimation, the species tree (
S
f(S,\Theta,G\midD)\proptof(D\midG)f(G\midS,\Theta)f(\Theta)f(S),
where
f(S)
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.[14] 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, the coordinated NNI, SPR and NodeSlider moves.[15] [16]
Consider for example the case of two species (A and B) and two sequences at each locus, with a sequence divergence time
ti
i
ti<\tau
i
\tau
ti
\tau
ti
\tau
ti
ti
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.[17] [18]
The basic multispecies coalescent model can be extended in a number of ways to accommodate major factors of the biological process of reproduction and drift.[9] [10] For example, incorporating continuous-time migration leads to the MSC+M (for MSC with migration) model, also known as the isolation-with-migration or IM models.[19] [20] Incorporating episodic hybridization/introgression leads to the MSC with introgression (MSci)[21] or multispecies-network-coalescent (MSNC) model.[22] [23]
The multispecies coalescent has profound implications for the theory and practice of molecular phylogenetics.[24] [25] Since individual gene trees can differ from the species tree one cannot estimate the tree for a single locus and assume that the gene tree correspond the species tree. In fact, one can be virtually certain that any individual gene tree will differ from the species tree for at least some relationships when any reasonable number of taxa are considered. However, gene tree-species tree discordance has an impact on the theory and practice of species tree estimation that goes beyond the simple observation that one cannot use a single gene tree to estimate the species tree because there is a part of parameter space where the most frequent gene tree is incongruent with the species tree. This part of parameter space is called the anomaly zone[26] and any discordant gene trees that are more expected to arise more often than the gene tree. that matches the species tree are called anomalous gene trees.
The existence of the anomaly zone implies that one cannot simply estimate a large number of gene trees and assume the gene tree recovered the largest number of times is the species tree. Of course, estimating the species tree by a "democratic vote" of gene trees would only work for a limited number of taxa outside of the anomaly zone given the extremely large number of phylogenetic trees that are possible.[27] However, the existence of the anomalous gene trees also means that simple methods for combining gene trees, like the majority rule extended ("greedy") consensus method or the matrix representation with parsimony (MRP) supertree[28] [29] approach, will not be consistent estimators of the species tree[30] [31] (i.e., they will be misleading). Simply generating the majority-rule consensus tree for the gene trees, where groups that are present in at least 50% of gene trees are retained, will not be misleading as long as a sufficient number of gene trees are used. However, this ability of the majority-rule consensus tree for a set of gene trees to avoid incorrect clades comes at the cost of having unresolved groups.
Simulations have shown that there are parts of species tree parameter space where maximum likelihood estimates of phylogeny are incorrect trees with increasing probability as the amount of data analyzed increases.[32] This is important because the "concatenation approach," where multiple sequence alignments from different loci are concatenated to form a single large supermatrix alignment that is then used for maximum likelihood (or Bayesian MCMC) analysis, is both easy to implement and commonly used in empirical studies. This represents a case of model misspecification because the concatenation approach implicitly assumes that all gene trees have the same topology.[33] Indeed, it has now been proven that analyses of data generated under the multispecies coalescent using maximum likelihood analysis of a concatenated data are not guaranteed to converge on the true species tree as the number of loci used for the analysis increases[34] [35] [36] (i.e., maximum likelihood concatenation is statistically inconsistent).
There are two basic approaches for phylogenetic estimation in the multispecies coalescent framework: 1) full-likelihood or full-data methods which operate on multilocus sequence alignments directly, including both maximum likelihood and Bayesian methods, and 2) summary methods, which use a summary of the original sequence data, including the two-step methods that use estimated gene trees as summary input and SVDQuartets, which use site pattern counts pooled over loci as summary input.
ASTRAL | ASTRAL (Accurate Species TRee ALgorithm) summarizes a set of gene trees using a quartet method generate an estimate of the species tree with coalescent branch lengths and support values (local posterior probabilities[37]) | Summary | Mirarab et al. (2014);[38] Zhang et al. (2018)[39] | |
ASTRID | ASTRID (Accurate Species TRees from Internode Distances) is an extension of the NJst method.[40] ASTRID/NJst is a summary species tree method that calculates the internode distances from a set of input gene trees. A distance method like neighbor joining or minimum evolution is then used to estimate the species tree from those distances. Note that ASTRID/NJst is not consistent under a model of missing data[41] | Summary | Vachaspati and Warnow (2015)[42] | |
BPP | Bayesian MCMC software package for inferring phylogeny and divergence times among populations under the multispecies coalescent process; also includes method for species delimitation | Full likelihood | Yang et al. (2015);[43] Flouri et al. (2018)[44] | |
STACEY | Bayesian MCMC software package for inferring phylogeny and divergence times among populations under the multispecies coalescent process; minimal clusters (samples assumed to belong to the same species according to the model) are sampled during the MCMC without the need to change parameters space | Full likelihood | Jones et al. (2015);[45] Jones GR (2018)[46] | |
| Bayesian MCMC software package for inferring phylogeny and divergence times among populations under the multispecies coalescent process. Implemented as part of the BEAST software package (pronounced Star BEAST) | Full likelihood | Heled and Drummond (2010)[47] | |
MP-EST | Accepts a set of gene trees as input and generates the maximum pseudolikelihood estimate of the species tree | Summary | Liu et al. (2010)[48] | |
SVDquartets (implemented in PAUP*) | PAUP* is a general phylogenetic estimation package that implements many methods. SVDquartets is a method that has shown to be statistically consistent for data generated given the multispecies coalescent | Summary/Site-pattern method | Chifman and Kubatko (2014)[49] |