Multispecies Coalescent and Quartet Species Trees

The Multispecies Coalescent (MSC) is a probabilistic population-genetic framework that models the genealogical histories of non-recombining genomic loci (gene trees) evolving within the demographic branches of an overarching species tree. Formulated through the integration of J. F. C. Kingman's (1982) n-coalescent with speciation models by Fumio Tajima (1983), Naoyuki Takahata (1989), and Bruce Rannala and Ziheng Yang (2003), the Multispecies Coalescent explains why individual orthologous gene trees frequently disagree with the true species phylogenetics even in the absence of horizontal gene transfer or paralogous gene duplication. Looking backward in time from the present, two homologous gene copies sampled from Sister Species A and B cannot coalesce until they enter their shared ancestral species population AB; however, once inside AB, they do not necessarily coalesce before reaching an even deeper ancestral speciation node ABC.

When ancestral polymorphism persists across consecutive speciation events so that gene lineages fail to coalesce within their immediate ancestral population, the resulting phenomenon is termed Incomplete Lineage Sorting (ILS) or deep coalescence. The probability of ILS depends on the ratio between the internode divergence time Δ*t* (in generations) and the diploid effective population size N_e of the ancestral population, parameterized in coalescent time units τ = Δ*t* / (2*N*_e). When rapid evolutionary radiations combine short internode times Δ*t* with large ancestral population sizes N_e (small τ), discordant gene topologies can outnumber the concordant topology in a region of parameter space known as the anomaly zone. Because unrooted four-taxon quartets have no anomaly zone under the MSC, quartet-based summary algorithms such as ASTRAL provide statistically consistent species-tree inference across genome-scale datasets.

The Three-Taxon Coalescent Process and Deep Coalescence

Consider a rooted three-taxon species tree ((A, B), C) in which sister species A and B diverged at time t₁ in the past, and their common ancestor AB diverged from outgroup species C at time t₂ > t₁. Suppose one haploid gene copy (a, b, c) is sampled from each species, and the ancestral branch connecting the (A, B) speciation event to the (AB, C) root has duration Δ*t* = t₂ − t₁ generations and constant diploid effective population size N_e (containing 2*N*_e gene copies). Under Kingman's continuous-time coalescent approximation, the rate at which two lineages coalesce in a panmictic diploid population of size N_e is 1 / (2*N*_e) per generation, meaning the waiting time to coalescence follows an exponential distribution with scale parameter τ = Δ*t* / (2*N*_e).

Two distinct genealogical outcomes can occur along the internal ancestral branch of length τ:

1. Coalescence within the ancestral branch (AB): Lineages a and b coalesce within the time interval τ before reaching the deeper root node t₂. This event occurs with probability 1 − exp(−τ) = 1 − e^{−Δ*t*/(2*N*_e)} and guarantees the concordant gene tree ((a, b), c).
2. Deep coalescence (Incomplete Lineage Sorting): Lineages a and b fail to coalesce within interval τ, occurring with probability exp(−τ) = e^{−Δ*t*/(2*N*_e)}. Consequently, all three lineages (a, b, and c) enter the deeper ancestral population ABC simultaneously as exchangeable lineages. Within ABC, each of the three possible pairs—(a, b), (a, c), or (b, c)—has an equal 1/3 probability of coalescing first.

Summing the probabilities of these two mutually exclusive pathways yields the fundamental three-taxon gene-tree probability distribution under the Multispecies Coalescent:

*P*[ ((a, b), c) ] = [ 1 − e^−τ ] + (1 / 3) · e^−τ = 1 − (2 / 3) · e^−τ

*P*[ ((a, c), b) ] = *P*[ ((b, c), a) ] = (1 / 3) · e^−τ

Because e^−τ < 1 for every positive branch length τ > 0, the concordant rooted triplet ((a, b), c) is always strictly more probable than either of the two discordant triplets ((a, c), b) or ((b, c), a):

1 − (2 / 3) · e^−τ > (1 / 3) · e^−τ for all τ > 0

Furthermore, the exact symmetry between the two discordant triplets (*P*[((a, c), b)] = *P*[((b, c), a)]) forms the basis of Patterson's D-statistic (the ABBA-BABA test): a statistically significant excess of one discordant topology over the other indicates post-speciation introgression or horizontal gene transfer rather than purely stochastic Incomplete Lineage Sorting.

Anomalous Gene Trees (AGTs) on Four or More Taxa

Although the most probable rooted triplet on three taxa always matches the species tree, James Degnan and Noah Rosenberg (2006) proved that for rooted trees with n ≥ 4 taxa (and unrooted trees with n ≥ 5 taxa), there exist combinations of short internal branch lengths where the single most probable gene tree topology does not match the species tree topology. Such discordant topologies are termed Anomalous Gene Trees (AGTs), and the set of internal branch lengths (x, y) producing them is called the anomaly zone.

Consider an asymmetric (caterpillar) four-taxon rooted species tree (((A, B), C), D) with two internal branches of coalescent lengths x (ancestral to A and B) and y (ancestral to AB and C). If all four lineages a, b, c, and d fail to coalesce in both internal branches and enter the root population ABCD together (which occurs with high probability when x and y are very small, e.g., x < 0.1 and y < 0.1 coalescent units), a symmetric balanced gene tree ((a, b), (c, d)) can be generated by two distinct temporal orderings of coalescent events in the root: either (a, b) coalesces first followed by (c, d), or (c, d) coalesces first followed by (a, b). By contrast, the concordant asymmetric gene tree (((a, b), c), d) requires a single specific ordering: (a, b) must coalesce first, and then the (ab) ancestor must coalesce with c before c coalesces with d. Consequently, in the limit as x, y → 0, each balanced topology has probability 2/18 = 1/9, whereas each asymmetric topology has probability 1/18. When x sits below the upper boundary of the anomaly zone (x ≈ 0.156 coalescent units as y → 0), the discordant balanced gene tree ((a, b), (c, d)) has a higher probability than the true species tree (((a, b), C), D).

Inconsistency of Supermatrix Concatenation

In traditional phylogenomics, alignments from hundreds or thousands of orthologous genes are concatenated end-to-end into a single supermatrix and analyzed with maximum-likelihood phylogenetics. Concatenation implicitly assumes that all sites across the genome share a single underlying topology. Laura Kubatko and L. Lacey Knowles (2007) and Sebastien Roch and Mike Steel (2015) proved that when internal branches fall inside the anomaly zone, maximum-likelihood analysis of concatenated alignments is statistically inconsistent: as the number of sampled loci m → ∞, concatenation converges with 100% bootstrap support onto the incorrect anomalous balanced topology rather than the true species tree.

Absence of Unrooted Quartet Anomaly Zones

To overcome the anomaly zone without the heavy computational burden of full Bayesian co-estimation of gene trees and species trees (such as *BEAST or BPP), summary coalescent methods decompose input n-taxon gene trees into their induced four-taxon unrooted quartets. For any four taxa {*A*, B, C, *D*}, there are exactly three possible unrooted quartet topologies: q₁ = AB | CD, q₂ = AC | BD, and q₃ = AD | BC. Remarkably, both the asymmetric rooted species tree (((A, B), C), D) and the symmetric rooted species tree ((A, B), (C, D)) induce the exact same unrooted quartet topology q₁ = AB | CD, and the probability of observing q₁ in an unrooted gene tree depends only on the single internal branch length τ separating (A, B) from (C, D):

P(AB | CD) = 1 − (2 / 3) · e^−τ > (1 / 3) · e^−τ = P(AC | BD) = P(AD | BC)

Because 1 − (2 / 3)*e*^−τ strictly exceeds (1 / 3)*e*^−τ for every τ > 0, there are no anomalous unrooted four-taxon gene trees under the Multispecies Coalescent.

Dynamic Programming Over Bipartitions in ASTRAL

Siavash Mirarab et al. (2014) leveraged this mathematical property in ASTRAL (Accurate Species TRee ALgorithm), which solves the Maximum Quartet Support Species Tree (MQSST) problem. Given a set of m unrooted gene trees G = {*g*₁, ..., g_m}, ASTRAL seeks a species tree S that maximizes the total number of induced four-taxon quartets shared across all m gene trees:

Score(S) = ∑_{g ∈ G} ∑_{Q ∈ Quartets(S)} I(Q ∈ Quartets(g))

Although the MQSST problem is NP-hard over the unconstrained space of all possible trees, ASTRAL restricts the candidate bipartitions (tripartitions around internal nodes) of S to a set X drawn from the bipartitions observed in the input gene trees (augmented by heuristic quartet additions in ASTRAL-II and ASTRAL-III). Using a dynamic programming recurrence over sub-clades in X without explicitly enumerating all (n choose 4) quartets, ASTRAL-III finds the optimal constrained species tree in O(n · m · |*X*|^1.726) time and computes local posterior probabilities (localPP) and internal branch lengths in coalescent units directly from empirical quartet frequencies f₁ via τ̂ = −ln[ 1.5 · (1 − f₁) ].

How are internal branch lengths estimated in coalescent units from gene-tree quartet frequencies?

Because the expected frequency of the dominant quartet topology under the Multispecies Coalescent is *E*[f₁] = 1 − (2 / 3)*e*^−τ, inverting this relationship yields the method-of-moments estimator τ̂ = −ln[ (3 / 2) · (1 − f₁) ], where f₁ ∈ (1/3, 1] is the empirical fraction of gene trees displaying the concordant quartet around that internal branch. If f₁ = 0.755 across 1,000 loci, the estimated internal branch length is τ̂ = −ln[1.5 × 0.245] = 1.00 coalescent unit (equivalent to 2*N*_e generations).

Why can single-gene barcoding trees disagree with the true species tree even when sequences are long and error-free?

Even an infinitely long, completely error-free alignment of a single non-recombining locus (such as mitochondrial DNA or a single nuclear exon) only reconstructs the genealogical history of that single locus (the gene tree), not the species tree. If the ancestral internode duration was short relative to effective population size (for example, τ = 0.25 coalescent units), there is a 2 × (1/3)*e*^−0.25 = 51.9% probability that the true genealogical history of that locus is discordant with the species tree due to Incomplete Lineage Sorting.

What is the difference between Incomplete Lineage Sorting (ILS) and Horizontal Gene Transfer (HGT)?

Incomplete Lineage Sorting is the vertical retention of ancestral polymorphism across speciation nodes; for a species tree ((A, B), C), ILS produces the two discordant topologies ((a, c), b) and ((b, c), a) in equal expected frequencies (1/3 · e^−τ each) and yields divergence times older than the speciation node t₂. By contrast, horizontal gene transfer or introgression between B and C after speciation creates a significant asymmetry favoring ((b, c), a) and produces sequence divergence times younger than the speciation split.

1. Kingman, J. F. C. (1982). "The coalescent." Stochastic Processes and their Applications, 13(3), 235–248. https://doi.org/10.1016/0304-4149(82)90011-490011-4)
2. Rannala, B., & Yang, Z. (2003). "Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci." Genetics, 164(4), 1645–1656. https://doi.org/10.1093/genetics/164.4.1645
3. Degnan, J. H., & Rosenberg, N. A. (2006). "Discordant of species trees with their most likely gene trees." PLoS Genetics, 2(5), e68. https://doi.org/10.1371/journal.pgen.0020068
4. Mirarab, S., Reaz, R., Bayzid, M. S., Zimmermann, T., Swenson, M. S., & Warnow, T. (2014). "ASTRAL: genome-scale coalescent-based species tree estimation." Bioinformatics, 30(17), i541–i548. https://doi.org/10.1093/bioinformatics/btu462