Module 6: Phylogeography & Molecular Clocks

John C. Avise and colleagues coined the term phylogeography in a landmark review (Avise, Arnold, Ball et al. 1987, Annual Review of Ecology and Systematics 18, 489–522) that synthesised a decade of mitochondrial-DNA work on vertebrate populations into an explicit discipline. The central insight was that intraspecific genealogies, when mapped back onto geography, reveal the historical demographic and biogeographic processes that shaped a species — glacial refugia, range expansions, vicariant events, secondary contact — at a resolution far below that of traditional species- level biogeography.

Avise’s “phylogeographic concordance” principle argues that when multiple unrelated taxa share the same phylogeographic break, the break reflects a shared biogeographic history rather than idiosyncratic species biology. Classic examples include the southeastern-versus-southwestern US Pleistocene discontinuity documented across > 30 vertebrate taxa (Soltis et al. 2006 Molecular Ecology), and the east–west Indo-Pacific break at the Sunda-Sahul shelf documented in reef fishes, corals and crustaceans (Avise 2009 Journal of Biogeography).

The initial molecular toolkit was narrow: animal mitochondrial DNA (mtDNA), especially the rapidly evolving control region (d-loop) for recent events and cytochrome oxidase I (COI) or cytochrome b for older events; chloroplast DNA (cpDNA) in plants; and allozymes. Over the past two decades this has broadened dramatically: nuclear single-copy loci like RAG1 or beta-fibrinogen intron; multi-locus SNP panels; reduced-representation RAD-seq (Miller et al. 2007, Andrews et al. 2016 Nat. Rev. Genet.); and whole-genome resequencing from both living and ancient specimens. The field’s statistical core, however, remained the same: coalescent theory.

Sir John Kingman (1982, Journal of Applied Probability 19A, 27–43 and related papers) constructed the mathematical foundation of modern phylogeography. The coalescent is the backwards-in-time process that traces the ancestry of a sample of size \(n\) drawn from a Wright–Fisher population of effective size \(N_e\). Kingman proved that under weak selection and in the limit of large \(N_e\), the lineages coalesce pairwise, and the waiting time \(T_k\) for the \(k\)th coalescence (from \(k\) to \(k-1\) lineages) is exponentially distributed.

Summing over all coalescences gives the expected time to most recent common ancestor:

Under the infinite-sites mutation model, the number of segregating sites \(S\) in a sample is Poisson with mean \(\theta\,H_{n-1}\) where \(\theta = 4N_e\mu L\) is Watterson’s theta. This directly links the observable polymorphism (SNPs, haplotype diversity) to the underlying demographic parameter \(N_e\), giving population geneticists a calibrated window onto population history.

Coalescent simulators like ms (Hudson 2002), msprime(Kelleher et al. 2016 PLoS Comp. Biol.), and SLiM (Haller & Messer 2019) implement the coalescent with arbitrarily complex demographic histories, structured populations, selection, and recombination. They are the workhorse of modern statistical phylogeography.

Zuckerkandl and Pauling (1965 Journal of Theoretical Biology) proposed the molecular-clock hypothesis: that if mutations accumulate roughly linearly with time, genetic divergence between species is a calibrated clock. The hypothesis was controversial for decades (neutral vs selectionist debates of the 1970s–80s) but is now accepted as approximately true for putatively neutral markers over sufficiently long time scales, with the important caveat that clock rates vary both among lineages and over time.

Animal mtDNA became the phylogeographic workhorse because of its surprisingly regular clock. Brown, George and Wilson (1979 PNAS) estimated ~2% sequence divergence per million years for mammalian mtDNA — the famous “2% rule.” Subsequent work has shown this applies best to protein-coding sites in mammals and birds; other taxa and loci have very different rates (insect COI ~1.15–3.5% per Myr, plant cpDNA ~10× slower than animal mtDNA).

Clock calibration comes from three sources: (1) fossil-calibrated node ages on phylogenies (e.g. crown Primates at 62–68 Ma), (2) biogeographically-dated vicariant events (Panamanian closure at ~3 Ma for marine geminate species; Messinian salinity crisis at 5.96–5.33 Ma), and (3) direct per-generation mutation-rate estimates from parent-offspring trios (Kong et al. 2012 Nature 488, 471–475, giving human de novo mutation rate ~1.2 × 10\(^{-8}\) per bp per generation).

Strict clocks frequently fail. Relaxed-clock methods allow rate variation across branches. Drummond & Rambaut (2007 BMC Evolutionary Biology) implemented this in BEAST using either an uncorrelated log-normal relaxed clock (branch rates drawn from a log-normal distribution) or an autocorrelated relaxed clock (rate on a branch correlated with the rate on its parent). BEAST2 (Bouckaert et al. 2014, 2019) and RevBayes (Höhna et al. 2016) now form the standard Bayesian toolkit. Fossilised birth-death models (Stadler 2010) integrate fossils as genealogically-tip-dated sequences, greatly reducing calibration bias.

Godfrey Hewitt’s synthesis (Hewitt 2000, Nature 405, 907–913; and 2004 Phil. Trans. Roy. Soc. B) established the Pleistocene refugial paradigm as a central organising principle for temperate phylogeography. The alternating glacial–interglacial cycles of the last 2.6 Myr, with a dominant 100 kyr cycle over the past ~800 kyr, repeatedly pushed mid-latitude biotas into southern refugia during glacial maxima and released them northward during interglacials.

For Europe, Hewitt identified three principal glacial refugia — Iberia, Italy, and the Balkans — each separated from the others by the Alpine–Pyrenean–Dinaric mountain barrier. After the Last Glacial Maximum (LGM, ~21 ka) each refugium seeded a re-colonisation of the European continent, producing three recurrent phylogeographic patterns:

• Grasshopper pattern (\(Chorthippus\,parallelus\)): only the Balkan refugium recolonised northern Europe; Iberian and Italian lineages remained trapped south of the Pyrenees and Alps.
• Bear pattern (\(Ursus\,arctos\)): Iberian and Balkan lineages met in central Europe, producing an east–west suture zone along the Rhine.
• Hedgehog pattern (\(Erinaceus\)): all three refugia contributed, producing a tri-partite suture in central Europe.

North America shows parallel but distinct refugia: the Beringian refugium in far-northwestern Alaska/Yukon (source of circumpolar boreal taxa), the Great Basin / southwestern refugia, Florida and the Gulf Coast, and isolated nunatak and fjord refugia along the Pacific coast (Shafer et al. 2010 Molecular Ecology). The “suture zones” where diverging lineages met during the Holocene correspond to Remington’s (1968) Pleistocene hybrid zones — the Appalachians, the Rockies–Great Plains boundary, and the Pacific Northwest.

Haplotype-diversity gradients are the phylogeographic fingerprint of post-glacial expansion: southern refugia retain high diversity, while northern leading edges show “allele-surfing” founder effects and low diversity (Edmonds et al. 2004 PNAS). The South-to-North diversity gradient is now a textbook diagnostic pattern for recent range expansion.

Barton and Hewitt (1985, Annual Review of Ecology and Systematics 16, 113–148) laid the theoretical foundation for hybrid-zone analysis with the tension-zone model: narrow clines maintained by a balance between dispersal pulling alleles across the zone and selection against hybrids purifying the zone. At steady state the cline width

Hybrid zones are natural laboratories for speciation. Classic systems include the European house-mouse (Mus musculus) hybrid zone at the domesticus/musculus boundary in central Europe; the European \(Bombina\) toad hybrid zone; the \(Anopheles\,gambiae\)/\(A.\,coluzzii\) contact in west Africa; and the recently documented songbird hybrid zones across North America (Toews et al. 2016 Current Biology).

Suture zones are geographic regions where many independent hybrid zones coincide, indicating a shared historical contact event. Remington (1968) catalogued five major North American suture zones; subsequent genomic work has validated most as Pleistocene re-contact zones. The Tibetan plateau and Hengduan Mountains are a major suture zone for Asian birds and mammals (Qu et al. 2014 Molecular Ecology).

The phylogeographic concordance principle (Avise 1992, 2000) predicts that different taxa, having experienced the same historical biogeographic event, should show concordant phylogenetic breaks at congruent ages. When observed, concordance is strong evidence for a shared biogeographic history; absence argues for idiosyncratic species-specific processes.

Classic successes include:

• Indo-Pacific Sunda–Sahul break: convergent east–west mtDNA splits in reef fishes, corals, gastropods and echinoderms with coalescent-dated splits clustering around the LGM Sundaland lowstand (Crandall et al. 2012 Nature Comm.).
• Central American Isthmus: geminate species pairs across the 3.1–3.5 Ma Panamanian closure in > 40 marine taxa (Knowlton & Weigt 1998 Proc. Roy. Soc. B), used as the canonical vicariant calibration for marine molecular clocks.
• Australian Monsoonal Tropics: convergent refugial breaks across frogs, snakes, geckos and birds of the Kimberley–Arnhem land barrier (Moritz et al. 2009 TREE).
• Neotropical Andean vicariance: shared east–west Andean breaks across rodents, frogs, birds and insects, dated variously to the late Miocene final uplift (Hoorn et al. 2010 Science).

Modern statistical comparative phylogeography tests concordance formally. Hickerson et al. (2006) developed msBayes for hierarchical ABC estimation of shared divergence times. The newer PhyloGeoCoal(Gehara 2017) and ABLE (Beeravolu 2018) frameworks provide genome-scale inference of co-divergence parameters.

Sewall Wright’s F-statistics (1951) quantify the partitioning of genetic variation across hierarchical population structure. At the most useful level:

\(F_{ST}\) of ~0.05 corresponds to ~5 migrants per generation, sufficient to homogenise allele frequencies; values > 0.15 indicate substantial structure. Weir & Cockerham (1984 Evolution) developed the standard estimator. Isolation-by-distance (Wright 1943; Slatkin 1993 Evolution) adds a spatial dimension: \(F_{ST}/(1-F_{ST}) \propto \log(\mathrm{geographic\,distance})\) in 2-D continua.

Model-based assignment methods — STRUCTURE (Pritchard, Stephens & Donnelly 2000), ADMIXTURE (Alexander et al. 2009), and sparse non-negative matrix factorisation in sNMF — infer individual admixture proportions directly from SNP genotypes without a priori population definitions. They have revolutionised fine-scale phylogeography by revealing cryptic population structure and historical admixture (Patterson, Price & Reich 2006).

Heng Li and Richard Durbin (2011, Nature 475, 493–496) introduced the Pairwise Sequentially Markovian Coalescent (PSMC), a transformative method that infers the history of effective population size from a single diploid genome. The two chromosome copies in an individual provide an implicit sample of size two; the distribution of TMRCAs along the genome, estimated from the density of heterozygous sites, encodes \(N_e(t)\) over many orders of magnitude of time.

PSMC and its successors — MSMC/MSMC2 (Schiffels & Durbin 2014 Nature Genetics) for multiple samples, SMC++ (Terhorst, Kamm & Song 2017 Nature Genetics) for thousands of genomes, and the deep-learning-based dnadna (Sanchez 2021) — form the modern toolkit. PSMC traces population history from roughly 20 ka to 5–10 Ma and has been applied to every vertebrate with a sequenced genome, from giant pandas (Zhao et al. 2013) to polar bears (Miller et al. 2012) to the woolly mammoth (Palkopoulou et al. 2015).

The most famous PSMC results include: (1) the human out-of-Africa bottleneck at ~60–70 ka with \(N_e\) dropping to ~1500 in non-African populations; (2) the broader Eurasian/Asian expansion after ~40 ka; (3) the independent population histories of great apes showing chimpanzees with consistently larger \(N_e\) than humans, and gorillas with a major bottleneck at ~1–2 Ma; and (4) the deep archaic-human structure revealed by comparison with Neanderthal and Denisovan genomes.

Direct sequencing of DNA from museum specimens, archaeological bones, and permafrost remains transformed phylogeography by providing actual temporal samples instead of inferred genealogies. The breakthrough paper by Pääbo (1985 Nature) demonstrated PCR amplification from Egyptian mummies; the 1990s saw recovery of mtDNA from Pleistocene permafrost-preserved mammoths and ground sloths. The ancient-DNA revolution accelerated dramatically with Illumina short-read sequencing after 2006 and the development of specialised damage-aware library-preparation protocols (Rohland et al. 2015 Phil. Trans.).

Temporal calibration of substitution rates has been the single most important methodological advance. Tip-dated phylogenetic inference using serially sampled ancient sequences (Rambaut 2000 Bioinformatics) reveals that within-species rates are typically several times faster than the long-term between-species rate — the “time-dependent rate phenomenon” (Ho, Phillips, Cooper & Drummond 2005 Molecular Biology and Evolution). This has led to major revisions of both inferred demographic histories and biogeographic ages.

Landmark ancient-DNA phylogeography projects include the woolly mammoth demographic history (Palkopoulou et al. 2015 Current Biology), the cave-bear genome (Bar-Gal 2013), horse domestication (Fages 2019 Cell), dog origins (Frantz et al. 2016 Science), and the peopling of the Americas (Raghavan et al. 2015 Science, Moreno-Mayar 2018 Science). In every case the ancient DNA overturned or refined a phylogeographic hypothesis previously built on modern samples alone.

The “next-gen” revolution (Illumina 2008; PacBio, Nanopore long-read 2015+) reduced per-genome sequencing costs from $1 billion (Human Genome Project) to under $1000 in 2015 and under $200 today. Whole-genome phylogeography is now routine for vertebrates, and increasingly for invertebrates, plants, and microbes. Projects such as Earth BioGenome (Lewin 2018 PNAS), Bat1K, B10K birds, and the Vertebrate Genomes Project aim at reference genomes for all vertebrate species.

Whole-genome data enable inference of: (1) recombination-rate variation across the genome, which shapes local coalescence times; (2) introgression and hybridisation detectable via D-statistics (Green et al. 2010 Science) and f-statistics (Patterson 2012); (3) selective sweeps differentiating populations at specific loci (Voight et al. 2006 PLoS Biology); (4) linkage-disequilibrium decay as an estimator of local \(N_e\) (McVean 2007); and (5) local ancestry tracts that reveal migration history at sub-chromosomal resolution (Hellenthal et al. 2014 Science).

Landscape genetics (Manel et al. 2003 TREE) combines phylogeography with spatially explicit landscape ecology, treating gene flow as a flux through a cost surface rather than through a featureless plane. The key statistical tools are isolation-by-resistance (McRae 2006 Evolution), which models gene flow as electrical conductance through a heterogeneous resistance grid, and circuit-theoretic connectivity metrics operationalised in the Circuitscape software.

Landscape-genetic studies have demonstrated that roads, rivers, urban areas, mountain ridges, and even linear agricultural features impose population- genetic structure at surprisingly small scales. Cushman et al. (2006 American Naturalist) showed that black-bear gene flow in the US Northwest is better predicted by terrain ruggedness than by Euclidean distance. Coulon et al. (2004 Molecular Ecology) mapped roe-deer gene flow through agricultural hedgerow networks. These studies feed directly into conservation corridor design (see Module 8).

Approximate Bayesian Computation (ABC; Beaumont, Zhang & Balding 2002 Genetics) provides the statistical backbone for modern comparative phylogeography. ABC dispenses with analytical likelihoods by simulating under a demographic model, computing summary statistics from the simulations, and retaining those close to the observed statistics. Posterior distributions of demographic parameters emerge directly from the retained simulations.

BEAST (Bayesian Evolutionary Analysis by Sampling Trees; Drummond & Rambaut 2007, 2012) implements full likelihood inference with relaxed clocks, substitution models, demographic priors, and discrete or continuous phylogeographic diffusion models (Lemey, Rambaut, Welch & Suchard 2010 Mol. Biol. Evol.). The Bayesian skyline/skyride plots (Drummond, Rambaut, Shapiro & Pybus 2005) reconstruct non-parametric \(N_e(t)\) from serially-sampled or calibrated trees, and have become a standard figure in phylogeographic publications.

Machine-learning approaches are rapidly taking over. Deep-learning inference tools (Flagel, Brandvain & Schrider 2019 MBE; ReLERNN for recombination) take raw genome summaries as input and infer demographic parameters via convolutional neural nets. Simulation-based likelihood-free inference will likely replace most ABC workflows over the next decade.