Maximum Likelihood & Bayesian Phylogenetics
Why should a biologist trust one evolutionary tree over another? Both trees might explain the same DNA sequences equally well by eye. This is the problem that model-based phylogenetics was built to solve. Where Cladistics & Maximum Parsimony asks only "which tree requires the fewest changes," likelihood-based methods ask a sharper question. They ask: "which tree, under an explicit model of how mutations accumulate, makes the observed sequences most probable?" That shift — from counting changes to modeling probabilities — is one of the most consequential methodological transitions in twentieth-century biology. The problem matters because DNA evolves unevenly. Some sites mutate fast, some slow. Transitions (A↔G, C↔T) happen more often than transversions. Branch lengths vary enormously across a tree. Parsimony silently assumes none of this complexity. When it is wrong, parsimony can be positively misled. This phenomenon is called long-branch attraction. In it, parsimony converges on an incorrect answer more confidently as more data are added. Maximum likelihood (ML) and Bayesian inference instead build an explicit, testable model of the substitution process. They then ask which tree topology and parameters make the data most probable (ML). Alternatively, they ask which trees are most probable given the data and prior beliefs (Bayesian). Both approaches also give you something parsimony structurally cannot: a native, interpretable measure of statistical confidence in each branch. The big picture is that phylogenetics graduated from a descriptive, algorithmic exercise into a full-blown statistical science. It borrowed the same log-likelihood machinery that underlies population genetics, epidemiology, and modern genomics. This brought costs — heavier computation, model-choice decisions, philosophical debates over priors. But it bought enormous power: the ability to test alternative hypotheses of relationship, to estimate divergence times, to detect selection. It also enabled building the vast, well-supported trees that now scaffold comparative biology. These range from The Tree of Life down to outbreak investigations of a single evolving virus. This article works through the substitution models that make likelihood calculations possible. It walks step by step through Felsenstein's pruning algorithm and Markov chain Monte Carlo (MCMC) tree sampling. It examines the classic simulation and empirical studies that exposed parsimony's failure modes and validated likelihood. It closes with worked numerical examples, common pitfalls, and the modern toolkit these ideas built.
- Maximum likelihood asks: given a substitution model, which tree and parameters make the observed sequence data most probable?
- Bayesian inference goes further, computing a posterior probability distribution over trees using Bayes' theorem and MCMC sampling.
- Felsenstein's pruning algorithm computes tree likelihood efficiently by summing over all unobserved ancestral states at each site.
- Substitution models (JC69, K80, GTR, +Γ, +I) correct for unequal rates and multiple hits that parsimony ignores.
- Long-branch attraction and simulation studies (Felsenstein 1978, Huelsenbeck & Hillis 1993) showed parsimony can be "positively misleading" — model-based methods are statistically consistent.
- Bootstrap support (ML) and posterior probability (Bayesian) quantify confidence in clades, but they are not directly comparable numbers.
Likelihood-based phylogenetics replaces "count the fewest changes" with "compute the probability of the data under an explicit model of evolution," making trees testable statistical hypotheses rather than parsimony optimizations.
The core idea
A phylogenetic likelihood is the probability of observing a given multiple sequence alignment, D, conditional on a tree topology τ, a set of branch lengths v, and a model of nucleotide (or amino acid, or codon) substitution M with parameters θ. Formally, L = Pr(D | τ, v, θ). Maximum likelihood phylogenetics searches the space of topologies and branch lengths to find the combination that maximizes this quantity. Because the sites in an alignment are treated as independent and identically distributed draws under the model (a simplification, since real sites can be correlated by structure or linkage), the total likelihood is the product of per-site likelihoods, and in practice biologists work with the log-likelihood, ln L, which turns that product into a more tractable sum. Bayesian phylogenetics reframes the same problem using Bayes' theorem: Pr(τ, v, θ | D) ∝ Pr(D | τ, v, θ) × Pr(τ, v, θ). The left side, the posterior probability, is the probability of the tree and parameters given the data; the likelihood is exactly the same quantity ML maximizes; and Pr(τ, v, θ) is the prior, the analyst's stated belief about trees and parameters before seeing the data. Because the posterior cannot be computed analytically for realistic models — the space of possible trees is astronomically large — Bayesian phylogenetics uses Markov chain Monte Carlo (MCMC) to sample trees in proportion to their posterior probability, yielding not a single best tree but a distribution of credible trees and a posterior probability for every clade in that distribution. Both frameworks depend on a substitution model: a continuous-time Markov chain describing the instantaneous rate of change from one nucleotide to another. The simplest, Jukes-Cantor (JC69), assumes all four bases are equally frequent and all substitution types equally likely. More realistic models — Kimura's two-parameter model (K80), the general time-reversible model (GTR), and rate-heterogeneity extensions using a gamma distribution (+Γ) or a proportion of invariable sites (+I) — relax these assumptions to match real sequence evolution, where transitions outpace transversions and some sites are functionally constrained. Choosing a model is itself a statistical decision, typically made by likelihood-ratio tests or information criteria (AIC, BIC) before or during tree search.
Mechanism
Step 1 — Build a substitution model. Nucleotide substitution is modeled as a continuous-time Markov process with an instantaneous rate matrix Q, whose off-diagonal entries q_ij give the rate of change from base i to base j. Exponentiating this matrix over a branch length t (measured in expected substitutions per site) gives the transition probability matrix P(t) = e^(Qt), whose entries P_ij(t) are the probability that a site starting as base i ends as base j after evolutionary "time" t. Under JC69, all four bases have equal frequency (0.25 each) and all substitution rates are equal; under GTR, the model allows six distinct exchangeability rates (A↔C, A↔G, A↔T, C↔G, C↔T, G↔T) and unequal base frequencies, at the cost of more parameters to estimate.
Step 2 — Compute the likelihood of a single site on a tree (Felsenstein's pruning algorithm). For a fixed tree with branch lengths, the likelihood of one alignment column is obtained by summing over all possible, unobserved ancestral nucleotide states at every internal node — a computation that would be prohibitively slow done naively (4 possible states at each of many internal nodes) but that Joseph Felsenstein showed in 1981 can be computed efficiently by a bottom-up recursion. At each tip, the "conditional likelihood" vector is 1 for the observed base and 0 for the others. Moving toward the root, at each internal node the conditional likelihood of state i is: L_i(node) = [Σ_j P_ij(t_left) L_j(left child)] × [Σ_k P_ik(t_right) L_k(right child)] This multiplies, over each descendant branch, the probability of transitioning to every possible descendant state weighted by that descendant's own conditional likelihood, then combines the two (or more) children multiplicatively because lineages evolve independently after they split. At the root, the site likelihood is the weighted sum over root states, weighted by the equilibrium base frequencies: L(site) = Σ_i π_i L_i(root). This recursive "pruning" reduces the computational cost from exponential to linear in the number of tips, which is why likelihood phylogenetics became computationally feasible at all.
Step 3 — Multiply across sites, search across trees. The full alignment's log-likelihood is the sum of the per-site log-likelihoods (assuming independence across sites, a simplifying assumption). ML tree search then explores combinations of topology and branch lengths — using hill-climbing heuristics like nearest-neighbor interchange (NNI), subtree pruning and regrafting (SPR), or tree bisection and reconnection (TBR) — to find the tree that maximizes this sum, because exhaustively evaluating all possible trees is impossible beyond a handful of taxa (there are already over 2 million possible unrooted trees for just 10 taxa).
Step 4 — Quantify confidence via bootstrapping (ML) or MCMC (Bayesian). For ML, Felsenstein's 1985 nonparametric bootstrap resamples alignment columns with replacement to generate pseudo-replicate datasets, re-estimates the ML tree for each, and reports the percentage of bootstrap replicates in which each clade reappears — the bootstrap support value. For Bayesian inference, MCMC (typically Metropolis-coupled MCMC, or MC³) proposes small changes to the tree (topology rearrangements, branch length or parameter tweaks), accepts or rejects each proposed state according to the Metropolis-Hastings ratio (which depends on the likelihood ratio and the prior ratio of the proposed versus current state), and — after discarding an initial "burn-in" period before the chain has converged — the sampled trees, sampled in exact proportion to their posterior probability, are summarized into a majority-rule consensus tree, where each clade's frequency across the sampled trees is its posterior probability.
Step 5 — Model selection matters. Before or during inference, biologists compare candidate substitution models (JC69, K80, HKY85, GTR, with or without +Γ rate heterogeneity or +I invariable sites) using likelihood-ratio tests (for nested models) or information criteria like AIC/BIC (for non-nested comparisons). Using an over-simplified model (e.g., JC69 on data with strong transition bias) can bias branch length estimates and, in extreme cases, topology itself; using an over-parameterized model wastes statistical power and increases variance. Rate heterogeneity across sites is nearly universal in real data — some codon positions or structural regions are far more conserved than others — and the discrete-gamma approximation (splitting a continuous gamma distribution of rates into typically 4 rate categories) is now close to a default assumption in serious analyses.
The evidence
Felsenstein (1978): parsimony can be positively misleading. Joseph Felsenstein published a landmark theoretical paper showing that under certain combinations of branch lengths — specifically, when two long branches are separated by short internal branches, a scenario now called long-branch attraction — maximum parsimony converges, with high statistical confidence, on the wrong tree as more sequence data are added. This is because parsimony ignores the possibility of convergent, parallel substitutions on unrelated long branches; when two lineages evolve fast, they independently accumulate the same changes by chance often enough to look like shared derived similarity (synapomorphy), and parsimony has no mechanism to discount this. Crucially, Felsenstein showed this failure is not a small-sample problem that more data fixes — it gets worse with more data, meaning parsimony is statistically inconsistent under these conditions. This result was the single strongest theoretical motivation for adopting explicit substitution models.
Huelsenbeck & Hillis (1993) and related simulation studies. Using computer simulations under known "true" trees with specified branch lengths and substitution models, researchers generated synthetic sequence data and then asked how often parsimony, distance methods, and likelihood recovered the correct topology. In regions of branch-length space prone to long-branch attraction, maximum likelihood substantially outperformed parsimony in recovering the true tree, precisely because ML's substitution model explicitly accounts for multiple hits (superimposed substitutions) and unequal rates, correcting for the very convergence that fools parsimony. These simulations, replicated and extended by many labs through the 1990s and 2000s, became the standard evidentiary basis for preferring likelihood-based methods whenever rate heterogeneity across lineages is suspected.
A real empirical case: the "Long-Branch" rooting of the animal phylogeny. Early parsimony-based analyses of 18S ribosomal RNA sequences in the 1990s at times placed fast-evolving taxa (like certain nematodes or microsporidian fungi) in artifactual positions at the base of animal or eukaryotic trees, apparently due to long-branch attraction between the fast-evolving taxon and the (also long-branched) outgroup used to root the tree. Subsequent re-analyses using explicit substitution models, denser taxon sampling (which breaks up long branches by adding intermediate-rate lineages, a strategy predicted by Felsenstein's theory to reduce long-branch attraction), and multi-gene likelihood or Bayesian analyses repositioned these taxa in agreement with morphological and later genomic evidence — a direct, applied demonstration that the theoretical failure mode Felsenstein described in 1978 was not merely hypothetical.
Felsenstein (1985): the bootstrap. Recognizing that ML trees, like parsimony trees, are point estimates without a native measure of uncertainty, Felsenstein adapted Bradley Efron's statistical bootstrap (1979) to phylogenetics: resample alignment columns with replacement to create pseudo-replicate datasets of the same length, reconstruct a tree for each, and tally how often each clade in the original tree reappears. This bootstrap support percentage became (and largely remains) the standard confidence measure reported on ML trees, though it is now understood to be a conservative measure of repeatability rather than a literal probability that a clade is "true."
Bayesian MCMC adoption (Yang & Rannala 1997; Huelsenbeck & Ronquist's MrBayes, 2001). Ziheng Yang and Bruce Rannala first applied MCMC to sample tree topologies from the posterior distribution in a phylogenetic context, and John Huelsenbeck and Fredrik Ronquist's MrBayes software made Bayesian phylogenetics broadly accessible to working biologists. Simulation and empirical benchmarking studies through the early 2000s showed Bayesian posterior probabilities often (though not always) agreed closely with ML topologies while requiring less computation time than bootstrapping thousands of ML replicates, driving rapid adoption — alongside continued debate (see below) about whether posterior probabilities are systematically too confident.
Worked examples
Example 1 — Computing a simple site likelihood by hand under JC69. Consider a 3-taxon unrooted tree: ((A,B),C), with branch lengths t₁ = 0.1 (A's tip branch), t₂ = 0.1 (B's tip branch), t₃ = 0.05 (internal branch to the ancestor of A+B), and t₄ = 0.2 (C's tip branch). Under JC69, the probability of no change on a branch of length t is P_same(t) = 1/4 + 3/4·e^(−4t/3), and the probability of change to any specific one of the other three bases is P_diff(t) = 1/4 − 1/4·e^(−4t/3). Suppose the observed column is A = "A", B = "A", C = "G". At the internal node (call it N, ancestor of A and B), we sum over N's possible states. For N = "A": contribution = P_same(0.1) [A→A] × P_same(0.1) [B→A, i.e., A→A] = P_same(0.1)². Plugging t = 0.1: e^(−4(0.1)/3) = e^(−0.133) ≈ 0.875, so P_same(0.1) ≈ 0.25 + 0.75(0.875) ≈ 0.906. So the N="A" term ≈ 0.906 × 0.906 ≈ 0.821. For N = any other base (say "G", "C", or "T"), both A and B would each require a change from that state to "A", contributing P_diff(0.1)² ≈ (0.031)² ≈ 0.00097 each — three such states contribute a total of ≈ 0.0029. So the summed conditional likelihood vector at N is approximately (0.821 for N=A; ≈0.0010 for each of N=C, G, T). Now combine with C's branch (t₄ = 0.2, giving P_same(0.2) ≈ 0.25 + 0.75·e^(−0.267) ≈ 0.25 + 0.75(0.766) ≈ 0.824, and P_diff(0.2) ≈ 0.0586) and the internal branch t₃ = 0.05 connecting N to the root. This full recursive combination (root frequencies π_i = 0.25 each under JC69) ultimately multiplies these conditional likelihoods together and sums over root states — the arithmetic is mechanical but illustrates the core point: most of the likelihood mass for this column comes from the "boring" hypothesis that A and B's shared ancestor was simply "A" and C independently changed to "G", exactly the kind of probabilistic reasoning parsimony cannot perform, since parsimony would simply count this column as requiring one change (on C's branch) regardless of branch lengths.
Example 2 — Interpreting discordant support values in a real dataset. A phylogeneticist reconstructs a 40-taxon ML tree for a group of freshwater fish using a GTR+Γ model, obtaining bootstrap support of 55% for a clade uniting genera Barbus and Puntius, while a Bayesian analysis on the same alignment gives that same clade a posterior probability of 0.98. Which number should the student trust, and why the discrepancy? The resolution requires understanding what each number measures. Bootstrap support estimates how often the data itself, when resampled, would recover this clade under repeated sampling — it is a frequentist measure of the clade's robustness to sampling variance in the alignment, and simulation studies (Hillis & Bull, 1993) have shown bootstrap values tend to underestimate the true accuracy of a clade at moderate values (a 70% bootstrap value can correspond to a clade that is correct >95% of the time). Posterior probability, by contrast, is the literal probability the clade is correct given the model and priors are correct — but multiple studies (Suzuki, Glazko & Nei, 2002; Simmons et al., 2004) have shown posterior probabilities computed under simplified models or with model misspecification can be systematically overconfident, especially for short internal branches, sometimes reporting >0.95 support for clades that are not robustly supported by independent data. A responsible interpretation here treats the 55%/0.98 discrepancy as a flag: check for model adequacy (was rate heterogeneity properly modeled?), check whether the two methods used the same substitution model, and seek corroborating evidence (additional genes, morphology) rather than mechanically trusting the larger number.
Common misconceptions & traps
| Misconception | Why it's wrong |
|---|---|
| "The ML tree is the true tree." | It is the tree that maximizes likelihood under a chosen model — model misspecification, insufficient data, or long-branch effects can still mislead ML, though less often than parsimony. |
| "Bayesian posterior probability = probability the clade is real." | Only true if the model and priors are correctly specified; in practice, model misspecification often inflates posterior probabilities, especially on short branches. |
| "Bootstrap support and posterior probability are on the same scale." | They estimate different things (repeatability under resampling vs. posterior probability under a model); a 70% bootstrap and a 0.95 posterior are not directly comparable numbers. |
| "More parameters in a model is always better." | Overparameterized models increase variance in parameter estimates and can reduce power; model selection (AIC/BIC, LRTs) should justify added complexity. |
| "Parsimony and likelihood always agree, so it doesn't matter which you use." | They provably diverge under long-branch attraction and other rate-heterogeneous scenarios — cases the whole model-based framework was built to fix. |
| "A well-supported clade in one gene tree is a well-supported clade in the species tree." | Gene trees can differ from the species tree due to incomplete lineage sorting, introgression, or horizontal transfer — a separate problem addressed by species-tree methods. |
The scientist behind it
Joseph Felsenstein, born in 1942, is the single most influential figure in transforming phylogenetics from a largely algorithmic, parsimony-driven discipline into a statistical one grounded in explicit models of molecular evolution. Trained originally in population genetics — a field steeped in probabilistic reasoning about allele frequencies since the era of Molecular Phylogenetics's founders — Felsenstein brought that same statistical instinct to the reconstruction of evolutionary trees at a time (the 1970s and early 1980s) when parsimony dominated systematics almost unchallenged. His 1978 paper proving that parsimony can be positively misleading under long-branch conditions was, at the time, a deeply controversial claim within a field that had only recently converted many practitioners to cladistic parsimony; it took over a decade of simulation studies and empirical vindication before likelihood methods became mainstream. His 1981 paper introducing the pruning algorithm for computing tree likelihoods, and his free software package PHYLIP (Phylogeny Inference Package, first released in 1980 and continuously updated since), gave the field both the mathematical machinery and the practical tools to actually use likelihood in daily research — PHYLIP remains one of the most widely cited software packages in the history of evolutionary biology.
Origins · 1 minJoseph Felsenstein and how this was discovered
Where this matters today
Public health and outbreak tracking. Real-time Bayesian phylogenetics (using tools built on the same MCMC principles, such as BEAST) reconstructs the transmission history and estimates the time to most recent common ancestor of pathogen outbreaks — HIV, Ebola, and SARS-CoV-2 lineages have all been dated and traced using molecular-clock-calibrated Bayesian trees, directly informing public health interventions and even legal investigations of transmission events.
Drug and vaccine design. ML and Bayesian trees of rapidly evolving viruses (influenza, HIV) identify which circulating strains are diverging fastest from vaccine targets, informing annual influenza vaccine strain selection by tracking antigenic drift along a statistically supported phylogeny.
Conservation genetics. Model-based phylogenies distinguish genuinely distinct evolutionary lineages (candidates for separate conservation management, sometimes called Evolutionarily Significant Units) from populations that merely look different, informing endangered species law and captive breeding programs.
Genome-scale systematics and the assembly of the tree of life. Modern phylogenomic studies analyze thousands of genes simultaneously using ML or Bayesian software (RAxML, IQ-TREE, MrBayes, BEAST, revBayes), reconciling individual gene trees that can disagree due to incomplete lineage sorting — a statistical challenge that itself required extending single-gene likelihood models into multispecies coalescent frameworks, directly building toward the ambitious project of resolving The Tree of Life.
Detecting natural selection. Codon-based likelihood models (comparing the rate of nonsynonymous to synonymous substitution, dN/dS) use the same likelihood machinery described here to test, gene by gene and even site by site, whether natural selection has been positive, purifying, or neutral — connecting phylogenetic methodology directly to functional and adaptive questions in molecular evolution.
Ask the tutor: How does incomplete lineage sorting complicate building a species tree from many individual Bayesian gene trees?
Why can't we just try every possible tree and find the true maximum-likelihood one?
The number of possible unrooted binary trees grows explosively with the number of taxa: for n taxa there are (2n−5)!! possible unrooted trees, which is already over 2 million for just 10 taxa and vastly larger for realistic datasets of 50–500 taxa. Exhaustive search is computationally impossible beyond a handful of species. Instead, ML methods use heuristic search algorithms — hill-climbing procedures like nearest-neighbor interchange (NNI), subtree pruning and regrafting (SPR), and tree bisection and reconnection (TBR) — that start from a reasonable initial tree (often built by a fast distance method) and iteratively rearrange it, keeping any change that improves the likelihood, until no local rearrangement improves the score further. This finds a good, often optimal or near-optimal tree without guaranteeing the global maximum, which is why multiple independent search replicates are standard practice.
Are posterior probabilities and bootstrap percentages interchangeable numbers?
No, and treating them as equivalent is a common student error. Bootstrap support answers a frequentist question: if we resampled columns from this exact alignment with replacement many times and reanalyzed, how often would this clade reappear? It measures repeatability given the specific dataset's information content. Posterior probability answers a Bayesian question: given the model, priors, and data, what is the probability this clade is correct? These are conceptually different quantities with different statistical properties — simulation studies show bootstrap values are often conservative (a 70% value can correspond to a clade correct more than 95% of the time), while posterior probabilities can be overconfident under model misspecification. Never directly compare a 70% bootstrap value to a 0.70 posterior probability as if they meant the same thing.
What exactly is "long-branch attraction" and why does more data make it worse under parsimony?
Long-branch attraction occurs when two lineages evolving quickly (long branches) are incorrectly grouped together in a reconstructed tree, even though they are not each other's closest relatives, because they have independently accumulated similar-looking changes by chance (convergent or parallel substitutions) rather than by shared ancestry. Maximum parsimony counts these coincidental shared states as evidence of relatedness because it cannot distinguish "inherited from a common ancestor" from "arrived at independently." Critically, as more sequence data accumulate, more of these chance similarities pile up between the two long branches, so parsimony becomes more confident in the wrong answer rather than converging on the truth — this is what makes parsimony statistically inconsistent in this scenario, a mathematical property Felsenstein proved in 1978.
Why do we need a substitution model at all — can't we just count sequence differences?
Simply counting the proportion of differing sites (the "p-distance") systematically underestimates true evolutionary divergence because it cannot see multiple substitutions that happened at the same site (a base that mutated twice and reverted, or mutated to a different base than the current comparison, are both invisible to raw counting) — a phenomenon called multiple hits or superimposed substitutions. As two sequences diverge for longer times, the true number of substitutions increasingly exceeds the observed number of differences, and this underestimation is worse on longer branches. Explicit substitution models like JC69 or GTR use the mathematics of Markov chains to correct for this hidden multiple-hit problem, recovering a much more accurate estimate of true evolutionary distance and, more importantly, using the same probabilistic logic to correctly weigh evidence for competing tree topologies.
What is model selection and why does it matter for the final tree?
Model selection is the process of choosing which substitution model (JC69, HKY85, GTR, with or without gamma-distributed rate heterogeneity or a proportion of invariable sites) best fits a given dataset before or during phylogenetic analysis. This matters because using an overly simple model can bias branch length estimates (typically underestimating them) and, in cases with strong rate heterogeneity or transition/transversion bias, can even mislead the tree topology recovered. Model selection is typically done via likelihood-ratio tests for nested models (e.g., testing whether adding gamma rate heterogeneity significantly improves the likelihood) or information criteria like AIC or BIC for comparing non-nested models, balancing model fit against complexity to avoid overfitting.
How does Bayesian MCMC actually decide which trees to keep exploring?
The Markov chain proposes a small change to the current tree — a topology rearrangement, a branch length adjustment, or a substitution parameter tweak — and then computes the Metropolis-Hastings acceptance ratio, which compares the posterior probability (likelihood times prior) of the proposed state to that of the current state. If the proposed state is more probable, it is always accepted; if less probable, it is accepted only with a probability equal to that ratio, allowing the chain to occasionally move "downhill" and avoid getting permanently stuck on a single local optimum. Run for millions of generations, discarding an early "burn-in" period before the chain has stabilized (converged) on the high-probability region of tree space, the retained samples approximate the true posterior distribution, and the frequency of each clade across those samples becomes its posterior probability.
Does using likelihood or Bayesian methods guarantee the correct tree?
No — both are powerful, statistically principled frameworks, but neither guarantees correctness. Their accuracy depends on the substitution model being a reasonable approximation of real evolution (model misspecification can still mislead both approaches), on having enough informative sequence data, on adequate taxon sampling (sparse sampling can leave long branches unbroken and still cause attraction artifacts even under ML), and on the actual history matching a single bifurcating tree at all (processes like hybridization, introgression, or horizontal gene transfer violate the very tree model being fit). What likelihood and Bayesian methods do guarantee, under reasonably correct models, is statistical consistency — accuracy that improves as more data are added, unlike parsimony's proven failure mode under long-branch conditions.