19.07.02 · eco-evo-bio / phylogenetics

Molecular clock hypothesis: calibration, rate heterogeneity, and Bayesian divergence dating

stub3 tiersLean: nonepending prereqs

Anchor (Master): Ho, S. Y. W. & Duchene, S. — Mol. Ecol. 23 (2014) 5947-5965

Intuition Beginner

In 1962, Zuckerkandl and Pauling made a striking observation: the number of amino-acid differences between the haemoglobin proteins of different species is roughly proportional to the time since those species shared a common ancestor, as independently estimated from fossils. Humans and gorillas differ at a few positions; humans and horses differ at many more; humans and sharks differ at still more. The pattern looks like a ticking clock — DNA and protein sequences accumulate changes over time at a roughly constant rate.

If you know how fast the clock ticks (the mutation rate , measured in substitutions per site per year) and you count the total number of differences between two species ( substitutions per site), you can estimate when they diverged:

where the factor of 2 accounts for both lineages accumulating substitutions independently since they split. This is the molecular clock.

The clock needs calibration — an independent source of timing information to set . Fossils provide this: if the oldest known fossil of a clade is 50 million years old, the clade's origin must be at least that old. Geological events (the formation of the Isthmus of Panama ~3 million years ago, the separation of New Zealand from Gondwana ~80 million years ago) also provide calibration points.

Two complications immediately arise. First, not all genes tick at the same speed. Mitochondrial genes evolve fast (~10x nuclear genes in animals); ribosomal RNA genes evolve slowly. Even within a gene, some sites are constrained by purifying selection while others are free to change. Second, the clock is not perfectly regular. Generation time, metabolic rate, DNA repair efficiency, and population size all affect substitution rates, causing different lineages to tick at different speeds. These are problems of rate heterogeneity, and modern methods relax the strict clock assumption to accommodate them.

Visual Beginner

The left panel shows an ultrametric tree under a strict clock: all root-to-tip distances are equal because every lineage evolves at the same rate. The right panel shows the same topology under a relaxed clock: branch lengths in substitutions differ across lineages, but the tree is still scaled to time using calibration points (marked with stars). The calibration anchors convert substitutions into millions of years.

Worked example Beginner

Two rodent species differ at 12% of sites in a mitochondrial gene. A fossil calibration indicates that their common ancestor lived 10 million years ago. What is the substitution rate?

Using , rearrange to get :

Now apply this rate to a second pair of species (same gene) that differ at 30% of sites. But 30% is large enough that multiple substitutions at the same site are likely (saturation), so we need a correction. Under the Jukes-Cantor model, the corrected distance is:

The divergence time is:

Without the correction, the naive estimate would be million years — a substantial underestimate. The Jukes-Cantor correction accounts for hidden multiple hits.

Check your understanding Beginner

Formal definition Intermediate+

The strict molecular clock

The strict molecular clock assumes that the substitution rate is constant across all branches of the phylogeny. For two sequences that diverged time units ago, the expected number of substitutions per site is:

where is the evolutionary distance (corrected for multiple hits using an appropriate substitution model) and the factor 2 arises because substitutions accumulate on both lineages independently. Given a calibrated rate , any pairwise distance yields a divergence time .

The relative rate test

The relative rate test (Sarich and Wilson 1967 Science 158, 1200-1203; Wu and Li 1985 Proc. Natl. Acad. Sci. USA 82, 1741-1745) checks whether two lineages evolve at the same rate without knowing . Given three taxa A, B, and an outgroup O, the distances satisfy and , where superscripts denote the lineage along which substitutions accumulate. Under a strict clock, , so . A statistically significant deviation from zero rejects the strict clock for the two lineages.

Calibration points

Divergence-time estimation requires at least one calibration point — a node whose age is independently constrained:

  • Fossil calibrations. The oldest known fossil assignable to a clade provides a minimum age. The prior on the node age is typically a lognormal distribution with the offset (minimum) set to the fossil age, reflecting the expected gap between the true divergence and the oldest known fossil.
  • Biogeographic calibrations. Geological events that create vicariance (e.g., the formation of the Isthmus of Panama ~3 Ma, the breakup of Gondwana) constrain the maximum age of splits between lineages separated by the event.
  • Tip dates. For rapidly evolving viruses sampled over time, each sample has a known collection date. The molecular clock rate is estimated directly from the temporal spread of sampling dates and genetic divergence.

Rate heterogeneity across sites and genes

Not all sites evolve at the same rate. The gamma-distributed rate variation model (Yang 1994 J. Mol. Evol. 39, 105-111) assigns each site a rate drawn from a gamma distribution with shape parameter and mean 1:

Small means strong rate heterogeneity (a few fast-evolving sites, many conserved ones); recovers uniform rates. In practice, a discrete approximation with rate categories is used, and the likelihood at each site is averaged over categories.

Partitioned models assign different substitution models and clock rates to different genes or codon positions, recognising that mitochondrial, nuclear, and chloroplast genes evolve at different tempos.

Relaxed clock models

The uncorrelated lognormal relaxed clock (Drummond et al. 2006 PLoS Biol. 4, e88) assigns each branch an independent rate drawn from a lognormal distribution:

where is the mean rate across branches and is the variance. The key assumption is that rates on adjacent branches are independent draws — no autocorrelation. If the posterior estimate of is near zero, the strict clock is adequate.

The autocorrelated relaxed clock (Kishino et al. 2001 Mol. Biol. Evol. 18, 2050-2060) assumes rates change gradually: , where is the branch duration and is the rate-change parameter.

Bayesian divergence dating

Bayesian divergence-time estimation computes the joint posterior distribution over tree topology , divergence times , substitution rates , and model parameters :

where is the sequence likelihood (computed by Felsenstein's pruning algorithm), is the tree prior (e.g., Yule process with speciation rate ), is the clock prior, and encompasses substitution model priors and calibration distributions. MCMC samples from this posterior, and the output is a set of dated trees summarised by a maximum clade credibility tree with 95% highest posterior density (HPD) intervals on node ages.

Counterexamples to common slips

  • The molecular clock is a universal constant. The clock rate varies across genes, across lineages, and across time. The "clock" is a useful approximation within a limited scope, not a law of nature. Neutral theory predicts rate constancy per year (or per generation) within a gene, but even neutral rates differ among lineages due to variation in generation time, population size, and repair efficiency.

  • More calibration points always improve dating. Poorly chosen calibrations (misidentified fossils, unjustified maximum bounds) introduce systematic error. Cross-validation of calibration points (Near and Sanderson 2004 Syst. Biol. 53, 62-69) can identify outliers that are inconsistent with the remaining data.

  • A relaxed clock is always better than a strict clock. When taxa are closely related and rate variation is minimal, the strict clock is more parameter-parsimonious and yields tighter credible intervals. The relaxed clock adds parameters that can inflate variance unnecessarily when data are limited.

Key theorem with proof Intermediate+

Theorem (Consistency of Bayesian divergence dating under the uncorrelated lognormal relaxed clock). Under the uncorrelated lognormal relaxed clock model with correctly specified substitution model, tree prior, and calibration distributions, the Bayesian posterior on divergence times concentrates around the true times as sequence length .

Proof sketch. The argument proceeds in two stages.

Stage 1: Branch lengths in substitutions per site are consistently estimated. As , the posterior on tree topology and branch lengths (in substitutions per site) concentrates on the true values by the consistency of Bayesian phylogenetic inference under a correctly specified substitution model (Chang 1996 Stat. Sin. 6, 143-164). The pruning algorithm computes site likelihoods that are identifiable — different yield different likelihoods — so the posterior concentrates.

Stage 2: Substitution rates and times are disentangled by calibration. Each branch has length , where is the rate and is the time duration. Under the uncorrelated lognormal prior, independently across branches. The tree prior (Yule or birth-death) constrains the relative times. At least one calibration prior with finite variance anchors the time scale. The joint prior on is proper (integrable) given the calibration, and the likelihood identifies (elementwise product). With consistently estimated and the prior separating rates from times, the posterior on concentrates around the true divergence times.

Bridge. This consistency result connects the substitution models of unit 19.07.01 to absolute time. The pruning algorithm provides the likelihood engine, and the relaxed clock with calibration priors provides the mechanism to decompose branch lengths into rates and times. The bridge is from comparative sequence data to a chronogram — a phylogeny whose nodes are dated in millions of years, enabling all downstream comparative methods.

Exercises Intermediate+

Relaxed clocks, calibration theory, and phylodynamics Master

BEAST2 workflows: tree priors and model selection

BEAST2 (Bouckaert et al. 2019 PLoS Comput. Biol. 15, e1006650) implements Bayesian divergence dating with flexible models for the tree-generating process, clock model, and substitution process. The choice of tree prior depends on the study system:

  • Yule process (pure birth). Nodes are added at rate (speciation rate) with no extinction. Appropriate for higher-level macroevolutionary trees where extinction is negligible relative to the timescale.
  • Birth-death process. Nodes are added at rate and removed at rate (extinction rate). The birth-death prior accounts for lineages that existed but left no living descendants, providing more realistic node-age distributions than Yule for deep trees with significant extinction.
  • Coalescent. For within-species data (population samples), the coalescent prior models the genealogy of sampled individuals as a function of effective population size and its changes over time (constant, exponential growth, Bayesian skyline).

Model selection among clock models (strict vs. relaxed), substitution models, and tree priors uses marginal-likelihood estimation via path sampling or stepping-stone sampling (Baele et al. 2012 Mol. Biol. Evol. 29, 2157-2167), comparing models by Bayes factors.

Fossilized birth-death process and total-evidence dating

The fossilized birth-death (FBD) process (Heath et al. 2014 Proc. Natl. Acad. Sci. USA 111, E2957-E2966; Stadler 2010 J. Theor. Biol. 265, 325-332) jointly models speciation, extinction, and fossil sampling as a single process. Under the FBD, fossils are treated as sampled ancestors (lineages sampled at a time before the present that may or may not be direct ancestors of living taxa) rather than as fixed calibrations on particular nodes. The model parameters are speciation rate , extinction rate , fossil recovery rate , and the proportion of extant taxa sampled .

Total-evidence dating (Ronquist et al. 2012 Syst. Biol. 61, 973-999) combines molecular data from extant taxa with morphological data from both extant and fossil taxa in a single Bayesian analysis under the FBD tree prior. The morphological data inform the phylogenetic placement of fossils; the FBD prior informs the timing. This approach avoids the circularity of fixing a fossil's phylogenetic position in node dating and instead treats it as an inference problem.

Node dating vs tip dating: soft bounds and hard bounds

In node dating, each calibration is a prior distribution on the age of a specific node. Hard bounds set absolute minimum and maximum ages (e.g., uniform(50, 100) Ma); soft bounds allow a small probability that the true age falls outside the specified range, implemented as exponential tails on the calibration distribution. Soft bounds (Yang and Rannala 2006 Syst. Biol. 55, 487-498) are preferred because hard bounds can introduce artefacts when the fossil record is incomplete or the geological age estimate is uncertain.

Cross-validation of calibrations (Near and Sanderson 2004 Syst. Biol. 53, 62-69) removes each calibration in turn, re-estimates the node age from the remaining data and calibrations, and compares the leave-one-out estimate to the calibration prior. Calibrations whose priors are strongly inconsistent with the data-driven estimates are flagged as potential outliers.

The time-dependent rate phenomenon

Empirical studies consistently find that estimated substitution rates decline with the time depth of the analysis: rates estimated from intra-species data are higher than those from inter-species comparisons (Ho et al. 2005 J. Mol. Evol. 61, 532-539; Ho et al. 2011 Mol. Ecol. 20, 3083-3094). Three mechanisms contribute:

  1. Purifying selection on slightly deleterious mutations. Polymorphisms that contribute to short-term rate estimates are removed by selection over longer timescales, reducing the apparent rate.
  2. Sequencing error and damage. In ancient DNA, post-mortem damage (cytosine deamination) inflates apparent substitution rates on short branches.
  3. Saturation and model underfitting. At long timescales, multiple substitutions at the same site obscure the true distance, and substitution models that fail to account for site heterogeneity underestimate the corrected distance, reducing the apparent rate.

Correcting for TDRP requires modelling the contribution of transient polymorphisms to the rate estimate. The gamma-Gaussian mixture approach (Ho et al. 2015) fits a model where the effective rate is a mixture of a long-term rate (fixed substitutions) and a short-term component (polymorphisms decaying at rate inversely proportional to ).

Incomplete lineage sorting and gene tree discordance effects on dating

When gene trees differ from the species tree due to ILS, each gene's divergence time reflects the gene coalescence time, which predates the speciation event. For a gene sampled from two species with divergence time and ancestral effective population size , the expected gene coalescence time is (in generations). Across many loci, the distribution of gene coalescence times is exponential with mean above .

The practical consequence: single-gene molecular clock estimates systematically overestimate species divergence times. Multi-locus approaches address this by: (1) estimating the species tree directly using the multispecies coalescent (StarBEAST2; Ogilvie et al. 2017 Mol. Biol. Evol. 34, 2101-2114), which separates coalescence times from speciation times; or (2) using the pseudo-timeline approach, which adjusts node dates by subtracting the expected ILS bias coalescent units.

Phylodynamics: real-time molecular clocks in epidemiology

Phylodynamics (Grenfell et al. 2004 Science 303, 327-332) applies molecular clock methods to viral sequence data sampled over short timescales (months to decades) to infer epidemiological dynamics. BEAST2 with a coalescent or skygrid tree prior estimates how the effective number of infections (where is generation time) changes over time from the shape of dated viral genealogies.

The molecular clock rate for RNA viruses is typically to substitutions per site per year — fast enough that measurable evolution occurs over epidemiological timescales. For SARS-CoV-2, the rate is approximately substitutions per site per year (~1-2 mutations per genome per month), enabling real-time tracking of variant emergence, geographic spread, and transmission dynamics. The skygrid coalescent prior (Gill et al. 2013 Mol. Biol. Evol. 30, 713-724) estimates a piecewise-constant trajectory, revealing epidemic growth and decline phases directly from sequence data without case-report data.

Full proof set Master

Proposition 1. Under the strict molecular clock with rate and the Jukes-Cantor model, the maximum likelihood estimate of the divergence time for two sequences with observed proportion of differences is .

Proof. The probability of observing a difference at a single site between two sequences that diverged time units ago is under JC69, where the factor 2 arises because substitutions accumulate on both lineages. For sites with differences, the log-likelihood is:

Setting and solving:

Let . Then , so . Substituting and simplifying:

where is the JC-corrected distance.

Proposition 2. For taxa and a strict molecular clock, there are distinct rooted, ranked tree topologies (trees that specify both branching order and the temporal order of nodes).

Proof. A rooted binary tree on taxa has internal nodes. Under a strict clock, the nodes must be ordered in time (no two nodes share the same age, almost surely). A ranked tree assigns a temporal ranking to the coalescence events. The number of ranked genealogies for labelled taxa is:

At the first coalescence event, choose any pair from lineages: choices. At the second event, choose from remaining lineages: choices. Continue until one pair remains. However, some ranked trees correspond to the same unranked topology. The number of distinct unranked rooted topologies is , so the number of rankings per topology varies. The total count of ranked labelled histories is , which is the number of distinct temporal orderings of coalescence events.

Proposition 3 (Effective sample size and MCMC convergence). For a Bayesian phylogenetic analysis with MCMC, the effective sample size (ESS) for a parameter computed from samples with autocorrelation at lag is , and convergence requires .

Proof. The MCMC chain produces samples that are autocorrelated. The variance of the sample mean is:

where is the marginal variance and is the autocorrelation at lag . The "effective" number of independent samples is the value such that , yielding . The threshold is a heuristic ensuring that the posterior mean is estimated with sufficient precision: with 200 independent samples, the standard error of the mean is , giving ~3 significant digits of accuracy.

Connections Master

  • Phylogenetic tree reconstruction 19.07.01. The molecular clock converts branch lengths from substitutions per site (the output of tree reconstruction) into absolute time. Without dating, a phylogeny specifies branching order but not when divergences occurred. The pruning algorithm from 19.07.01 computes the likelihoods that Bayesian dating methods use to weight candidate chronograms.

  • Neutral theory 19.04.02 pending. The molecular clock hypothesis was originally justified by neutral theory: if most substitutions are selectively neutral, their rate of fixation equals the mutation rate independently of population size. This provides the mechanistic basis for rate constancy. Departures from neutrality (positive selection, purifying selection) are the primary causes of rate heterogeneity across genes and lineages.

  • Coalescent theory 19.04.03 pending. The coalescent provides the tree prior for Bayesian divergence dating of within-species data. The coalescent waiting times, parameterised by , set the expected distribution of node ages, and the Bayesian skyline and skygrid models extend the coalescent to allow to vary over time — connecting phylodynamic inference to population history.

  • Wright-Fisher model 19.02.05. The effective population size that enters coalescent-based dating is defined by the Wright-Fisher model. The connection is: determines coalescence rates, which determine the distribution of gene-tree node ages relative to speciation events, which biases molecular clock estimates when ILS is present.

Historical & philosophical context Master

The molecular clock hypothesis originated with Zuckerkandl and Pauling's 1962 paper "Molecular disease, evolution, and genic heterogeneity" (in Horizons in Biochemistry, Academic Press, pp. 189-225) [Zuckerkandl & Pauling 1962]. Comparing haemoglobin sequences across species, they observed that the degree of molecular divergence correlated with the fossil-estimated divergence times — sequences from closely related species were more similar than those from distantly related ones, in rough proportion to the time since their last common ancestor. In a follow-up paper (Zuckerkandl and Pauling 1965, in Evolving Genes and Proteins, Academic Press, pp. 97-166), they explicitly proposed the term "molecular evolutionary clock."

The concept was formalised and extended by Sarich and Wilson (1967 Science 158, 1200-1203), who used immunological distances (a proxy for sequence divergence) to date the human-chimpanzee divergence to approximately 5 million years — far more recent than the 15 million years then favoured by paleoanthropologists. This sparked a major controversy but was ultimately vindicated by DNA sequence data.

The neutral theory of molecular evolution (Kimura 1968 Nature 217, 624-626; King and Jukes 1969 Science 164, 788-798) provided the theoretical justification: if most substitutions are selectively neutral, their fixation rate equals the mutation rate per generation, independent of population size, producing a clock-like accumulation of changes. The generation-time problem — organisms with shorter generations should show faster per-year clocks — was addressed by the nearly neutral theory (Ohta 1972 Nature 242, 194-196), which predicts that slightly deleterious mutations behave neutrally in small populations but are purged in large populations, introducing dependence on both generation time and population size.

The strict clock assumption was progressively relaxed. Sanderson (1997 Mol. Biol. Evol. 14, 1218-1231) introduced nonparametric rate smoothing (NPRS), which penalised rate changes between adjacent branches. Thorne, Kishino, and Painter (1998 Mol. Biol. Evol. 15, 1647-1660) developed the first Bayesian approach to divergence dating with autocorrelated rates. Drummond et al. (2006 PLoS Biol. 4, e88) introduced the uncorrelated lognormal relaxed clock in BEAST, which became the standard tool for Bayesian divergence dating. The fossilized birth-death process (Stadler 2010; Heath et al. 2014) further integrated fossil data into the dating framework through total-evidence approaches.

The philosophical significance of the molecular clock is that it provides an independent temporal framework for evolutionary history — one not dependent on the incompleteness of the fossil record. Molecular dates have pushed back the origins of many groups (e.g., placental mammals to before the K-Pg boundary, contradicting the "explosive model" derived from fossils alone), while also revealing that the fossil record, though incomplete, is more accurate in its relative timing than was once assumed. The tension between molecular and fossil dates remains one of the active frontiers of evolutionary biology, driving methodological improvements in both fields.

Bibliography Master

@incollection{ZuckerkandlPauling1962,
  author = {Zuckerkandl, E. and Pauling, L.},
  title = {Molecular disease, evolution, and genic heterogeneity},
  booktitle = {Horizons in Biochemistry},
  editor = {Kasha, M. and Pullman, B.},
  publisher = {Academic Press},
  pages = {189--225},
  year = {1962}
}

@incollection{ZuckerkandlPauling1965,
  author = {Zuckerkandl, E. and Pauling, L.},
  title = {Evolutionary divergence and convergence in proteins},
  booktitle = {Evolving Genes and Proteins},
  editor = {Bryson, V. and Vogel, H. J.},
  publisher = {Academic Press},
  pages = {97--166},
  year = {1965}
}

@article{SarichWilson1967,
  author = {Sarich, V. M. and Wilson, A. C.},
  title = {Immunological time scale for hominid evolution},
  journal = {Science},
  volume = {158},
  pages = {1200--1203},
  year = {1967}
}

@article{Kimura1968,
  author = {Kimura, M.},
  title = {Evolutionary rate at the molecular level},
  journal = {Nature},
  volume = {217},
  pages = {624--626},
  year = {1968}
}

@article{Ohta1972,
  author = {Ohta, T.},
  title = {Evolutionary rate of cistrons and evolutionary divergence},
  journal = {Nature},
  volume = {242},
  pages = {194--196},
  year = {1972}
}

@article{Yang1994,
  author = {Yang, Z.},
  title = {Maximum likelihood phylogenetic estimation from {DNA} sequences with variable rates over sites: approximate methods},
  journal = {J. Mol. Evol.},
  volume = {39},
  pages = {105--111},
  year = {1994}
}

@article{Kishino2001,
  author = {Kishino, H. and Thorne, J. L. and Bruno, W. J.},
  title = {Performance of a divergence time estimation method under a probabilistic model of rate evolution},
  journal = {Mol. Biol. Evol.},
  volume = {18},
  pages = {2050--2060},
  year = {2001}
}

@article{Drummond2006,
  author = {Drummond, A. J. and Ho, S. Y. W. and Phillips, M. J. and Rambaut, A.},
  title = {Relaxed phylogenetics and dating with confidence},
  journal = {PLoS Biol.},
  volume = {4},
  pages = {e88},
  year = {2006}
}

@article{Sanderson1997,
  author = {Sanderson, M. J.},
  title = {A nonparametric approach to estimating divergence times in the absence of rate constancy},
  journal = {Mol. Biol. Evol.},
  volume = {14},
  pages = {1218--1231},
  year = {1997}
}

@article{Thorne1998,
  author = {Thorne, J. L. and Kishino, H. and Painter, I. S.},
  title = {Estimating the rate of evolution of the rate of molecular evolution},
  journal = {Mol. Biol. Evol.},
  volume = {15},
  pages = {1647--1660},
  year = {1998}
}

@article{Heath2014,
  author = {Heath, T. A. and Huelsenbeck, J. P. and Stadler, T.},
  title = {The fossilized birth-death process for coherent calibration of divergence-time estimates},
  journal = {Proc. Natl. Acad. Sci. USA},
  volume = {111},
  pages = {E2957--E2966},
  year = {2014}
}

@article{Ronquist2012,
  author = {Ronquist, F. and Klopfstein, S. and Vilhelmsen, L. and Schulmeister, S. and Murray, D. L. and Rasnitsyn, A. P.},
  title = {A total-evidence approach to dating with fossils, applied to the early radiation of the {H}ymenoptera},
  journal = {Syst. Biol.},
  volume = {61},
  pages = {973--999},
  year = {2012}
}

@article{YangRannala2006,
  author = {Yang, Z. and Rannala, B.},
  title = {Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds},
  journal = {Syst. Biol.},
  volume = {55},
  pages = {487--498},
  year = {2006}
}

@article{NearSanderson2004,
  author = {Near, T. J. and Sanderson, M. J.},
  title = {Assessing the robustness of divergence time estimates},
  journal = {Syst. Biol.},
  volume = {53},
  pages = {62--69},
  year = {2004}
}

@article{Ho2011,
  author = {Ho, S. Y. W. and Lanfear, R. and Bromham, L. and Phillips, M. J. and Soubrier, J. and Rodrigo, A. G. and Cooper, A.},
  title = {Time-dependent rates of molecular evolution},
  journal = {Mol. Ecol.},
  volume = {20},
  pages = {3083--3094},
  year = {2011}
}

@article{HoDuchene2014,
  author = {Ho, S. Y. W. and Duchene, S.},
  title = {Molecular-clock methods for estimating evolutionary rates and timescales},
  journal = {Mol. Ecol.},
  volume = {23},
  pages = {5947--5965},
  year = {2014}
}

@article{Bouckaert2019,
  author = {Bouckaert, R. and Vaughan, T. G. and Barido-Sottani, J. and Duchene, S. and Fourment, M. and Gavryushkina, A. and Heled, J. and Jones, G. and Kuhnert, D. and De Maio, N. and Matschiner, M. and Mendes, F. K. and Muller, N. F. and Ogilvie, H. A. and du Plessis, L. and Popinga, A. and Rambaut, A. and Rasmussen, D. and Siveroni, I. and Suchard, M. A. and Wu, C.-H. and Xie, D. and Zhang, C. and Stadler, T. and Drummond, A. J.},
  title = {{BEAST} 2.5: an advanced software platform for {B}ayesian evolutionary analysis},
  journal = {PLoS Comput. Biol.},
  volume = {15},
  pages = {e1006650},
  year = {2019}
}

@article{Grenfell2004,
  author = {Grenfell, B. T. and Pybus, O. G. and Gog, J. R. and Wood, J. L. N. and Daly, J. M. and Mumford, J. A. and Holmes, E. C.},
  title = {Unifying the epidemiological and evolutionary dynamics of pathogens},
  journal = {Science},
  volume = {303},
  pages = {327--332},
  year = {2004}
}

@article{Ogilvie2017,
  author = {Ogilvie, H. A. and Heled, J. and Xie, D. and Drummond, A. J.},
  title = {{Computational performance and statistical accuracy of *BEAST} and comparisons with other methods},
  journal = {Mol. Biol. Evol.},
  volume = {34},
  pages = {2101--2114},
  year = {2017}
}

@article{Gill2013,
  author = {Gill, M. S. and Lemey, P. and Faria, N. R. and Rambaut, A. and Shapiro, B. and Suchard, M. A.},
  title = {Improving {B}ayesian population dynamics inference: a coalescent-based model for multiple loci},
  journal = {Mol. Biol. Evol.},
  volume = {30},
  pages = {713--724},
  year = {2013}
}

@book{FutuymaKirkpatrick2017,
  author = {Futuyma, D. J. and Kirkpatrick, M.},
  title = {Evolution},
  edition = {4th},
  publisher = {Sinauer Associates},
  year = {2017}
}

@book{HartlClark2007,
  author = {Hartl, D. L. and Clark, A. G.},
  title = {Principles of Population Genetics},
  edition = {4th},
  publisher = {Sinauer Associates},
  year = {2007}
}