Coalescent theory: the genealogy of a sample, TMRCA, and linkage disequilibrium

Intuition Beginner

Instead of tracking an entire population forward in time, coalescent theory looks backward: take a sample of genes alive today and trace their ancestry into the past. Every pair of genes shares a common ancestor at some point in the past -- the point where their lineages coalesce (merge). Continuing backward, all lineages eventually merge into a single ancestral gene, the most recent common ancestor (MRCA) of the entire sample. The time to reach this ancestor is the time to the most recent common ancestor (TMRCA).

The TMRCA depends on population size. In a large population, two randomly chosen genes are unlikely to share a parent in the previous generation, so lineages take a long time to coalesce -- deep genealogies. In a small population, the chance that two genes share a recent parent is high, so lineages coalesce quickly -- shallow genealogies. Roughly, the TMRCA of a sample scales with the effective population size $N_e$.

Genes that sit close together on the same chromosome are linked: they tend to be inherited together and share the same genealogical history. Genes far apart (or on different chromosomes) recombine more often and have independent genealogies. The non-random association between alleles at nearby loci is called linkage disequilibrium (LD). High LD means that knowing the allele at one locus tells you a lot about the allele at the neighbouring locus. LD decays with physical distance because recombination breaks up the associations over time.

Visual Beginner

The left panel shows a Kingman coalescent tree. The first coalescence involves two of six lineages; coalescence rates increase as fewer lineages remain (each pair coalesces at rate $1/(2N_e)$, so more pairs means faster merger). The right panel shows the characteristic decay of LD with distance: nearby loci are strongly correlated, distant loci are nearly independent.

Worked example Beginner

In a population with effective size $N_e=10,000$, what is the expected TMRCA for a sample of $n=2$ genes?

For two genes, there is only one pair of lineages. The coalescence rate is $1/(2N_e)=1/20,000$ per generation. The expected waiting time is:

$$E[T_{\text{MRCA}}]=2N_e=20,000\text{ generations}.$$

For a human population with generation time 25 years, this is $20,000\times 25=500,000$ years. Any two randomly chosen human gene copies share a common ancestor roughly half a million years ago on average.

For a sample of $n=10$ genes, the TMRCA is approximately $4N_e(1-1/n)=4\times 10,000\times 0.9=36,000$ generations, or about 900,000 years. Adding more samples increases the TMRCA only slowly (logarithmically) because the extra lineages coalesce quickly near the tips.

Check your understanding Beginner

Formal definition Intermediate+

The Kingman coalescent

The Kingman coalescent (Kingman, 1982) is a continuous-time Markov process that describes the genealogy of a sample of $n$ genes drawn from a diploid population of effective size $N_e$. The state of the process at any time is a partition of the $n$ genes into blocks of mutually coalesced lineages. When $k$ lineages remain, any specific pair of lineages coalesces at rate $1/(2N_e)$, and there are $\left(\genfrac{}{}{0ex}{}{k}{2}\right)$ pairs. The total coalescence rate is:

$${\lambda }_k=\left(\genfrac{}{}{0ex}{}{k}{2}\right)\cdot \frac{1}{2N_e}=\frac{k(k-1)}{4N_e}.$$

The waiting time $T_k$ until the next coalescence (reducing $k$ lineages to $k-1$) is exponentially distributed with rate ${\lambda }_k$:

$$T_k\sim \text{Exponential}\left(\frac{k(k-1)}{4N_e}\right),E[T_k]=\frac{4N_e}{k(k-1)}.$$

Expected TMRCA

The time to the most recent common ancestor (TMRCA) of the entire sample of $n$ genes is the sum of the waiting times from $k=n$ down to $k=2$:

$$T_{\text{MRCA}}=\sum _{k=2}^n{T}_k.$$

Taking expectations:

$$E[T_{\text{MRCA}}]=\sum _{k=2}^n\frac{4N_e}{k(k-1)}=4N_e\sum _{k=2}^n\left(\frac{1}{k-1}-\frac{1}{k}\right)=4N_e\left(1-\frac{1}{n}\right).$$

For large $n$, $E[T_{\text{MRCA}}]\approx 4N_e$. For $n=2$, $E[T_{\text{MRCA}}]=2N_e$.

Total tree length

The total tree length $L_{\text{total}}$ (the sum of all branch lengths) determines how many mutations the sample is expected to carry. During the interval when $k$ lineages are present, $k$ branches extend simultaneously:

$$E[L_{\text{total}}]=\sum _{k=2}^{n}k\cdot E[T_k]=\sum _{k=2}^n\frac{4N_e}{k-1}=4N_e\sum _{j=1}^{n-1}\frac{1}{j}=4N_e{H}_{n-1},$$

where $H_{n-1}={\sum }_{j=1}^{n-1}1/j$ is the $(n-1)$-th harmonic number. For large $n$, $H_{n-1}\approx \ln (n)+\gamma$ (Euler's constant $\gamma \approx 0.577$).

Tree shape and population history

The shape of the coalescent tree records population history:

• Constant population size. The standard coalescent produces trees where most coalescences occur near the tips (many lineages, high coalescence rate) and the last two lineages take a long time to coalesce. External branches are short relative to internal branches.

• Population growth. In a rapidly expanding population, most coalescences are pushed deep into the past when the population was small. The tree becomes star-like: long external branches and short internal branches. This produces an excess of rare variants (singletons), detected as significantly negative Tajima's D 19.04.02 pending.

• Population bottleneck. A bottleneck (temporary reduction in $N_e$) increases the coalescence rate during the bottleneck period, producing a cluster of coalescences at that time. The tree has long internal branches flanking the bottleneck. After the bottleneck, the surviving lineages take a long time to coalesce in the expanded population. The signature is a deficit of rare variants (significantly positive Tajima's D) and an excess of intermediate-frequency variants.

The structured coalescent

When a population is subdivided into $d$ demes with migration rates $m_{ij}$ between demes $i$ and $j$, the structured coalescent adds migration events to the backward process. A lineage in deme $i$ migrates backward to deme $j$ at rate $m_{ij}$, while two lineages in deme $i$ coalesce at rate $1/(2N_i)$. The probability that two lineages from the same deme coalesce before either migrates determines $F_{ST}$ 19.04.01: lineages from different demes have longer expected coalescence times because they must first migrate into the same deme before coalescing.

Linkage disequilibrium

Linkage disequilibrium (LD) measures the non-random association between alleles at two loci. Consider two biallelic loci with alleles $A/a$ and $B/b$. Let $p_{AB}$, $p_{Ab}$, $p_{aB}$, $p_{ab}$ be the haplotype frequencies. The standard LD measure is:

$$D=p_{AB}\cdot p_{ab}-p_{Ab}\cdot p_{aB}.$$

Under linkage equilibrium (random association), $D=0$. Three normalised measures are in common use:

• $D^{\prime }$ (Lewontin's normalised $D$): $D^{\prime }=D/D_{\max }$, where $D_{\max }=\min (p_A{p}_b,p_a{p}_B)$ if $D>0$ or $\min (p_A{p}_B,p_a{p}_b)$ if $D<0$. $D^{\prime }$ ranges from $-1$ to $+1$ and reaches $\pm 1$ whenever at least one haplotype is absent.

• $r^2$ (squared correlation coefficient): $r^2=D^2/(p_A{p}_a{p}_B{p}_b)$. $r^2$ ranges from 0 to 1 and is the measure most used in GWAS because of its direct relationship to statistical power for tagging SNP selection.

Decay of LD with distance

Recombination breaks down LD between two loci at rate $c$ (the recombination fraction per generation). The expected value of $D$ after $t$ generations satisfies:

$$E[D_t]=D_0(1-c{)}^t.$$

At equilibrium between drift (which creates LD by sampling) and recombination (which destroys it), the expected $r^2$ between two loci at recombination distance $c$ is approximately:

$$E[r^2]\approx \frac{1}{1+4N_{e}c}=\frac{1}{1+\rho },$$

where $\rho =4N_{e}c$ is the population-scaled recombination rate. When $\rho \ll 1$ (tightly linked loci in a small population), $r^2\approx 1$ and LD is high. When $\rho \gg 1$ (loosely linked loci or large population), $r^2\approx 0$ and LD is negligible.

Haplotype blocks

In practice, the genome is organised into haplotype blocks: regions of high LD within which a few common haplotypes account for most of the variation, separated by regions of low LD (recombination hotspots). The size of haplotype blocks depends on $N_e$ and the local recombination rate. In human populations, haplotype blocks in regions of low recombination can span hundreds of kilobases, while blocks in recombination hotspots may be only a few kilobases long. Populations that have experienced recent bottlenecks (e.g., European populations) have longer haplotype blocks than populations with larger long-term $N_e$ (e.g., African populations), because bottlenecks reduce effective recombination by reducing the number of distinct haplotypes.

Key theorem with proof Intermediate+

Theorem (Expected TMRCA of a sample of $n$ genes under the Kingman coalescent). In a diploid population of effective size $N_e$, the expected time to the most recent common ancestor of a sample of $n$ gene copies is $E[T_{\text{MRCA}}]=4N_e(1-1/n)$.

Proof. The Kingman coalescent with $k$ lineages has coalescence rate ${\lambda }_k=\left(\genfrac{}{}{0ex}{}{k}{2}\right)/(2N_e)=k(k-1)/(4N_e)$. The waiting time for the next coalescence is $T_k\sim \text{Exponential}({\lambda }_k)$ with $E[T_k]=1/{\lambda }_k=4N_e/[k(k-1)]$.

The TMRCA is the sum $T_{\text{MRCA}}=T_n+T_{n-1}+\cdots +T_2$. By linearity of expectation:

$$E[T_{\text{MRCA}}]=\sum _{k=2}^n\frac{4N_e}{k(k-1)}=4N_e\sum _{k=2}^n\left(\frac{1}{k-1}-\frac{1}{k}\right).$$

The sum telescopes:

$$\sum _{k=2}^n\left(\frac{1}{k-1}-\frac{1}{k}\right)=\left(\frac{1}{1}-\frac{1}{2}\right)+\left(\frac{1}{2}-\frac{1}{3}\right)+\cdots +\left(\frac{1}{n-1}-\frac{1}{n}\right)=1-\frac{1}{n}.$$

Therefore $E[T_{\text{MRCA}}]=4N_e(1-1/n)$. $\square$

Corollary (Expected total tree length). The expected total tree length is $E[L_{\text{total}}]=4N_e{H}_{n-1}$ where $H_{n-1}={\sum }_{j=1}^{n-1}1/j$ is the $(n-1)$-th harmonic number.

Proof. During the interval $T_k$, exactly $k$ branches extend. So:

$$E[L_{\text{total}}]=\sum _{k=2}^{n}k\cdot E[T_k]=\sum _{k=2}^n\frac{4N_e}{k-1}=4N_e\sum _{j=1}^{n-1}\frac{1}{j}=4N_e{H}_{n-1}.\square$$

Bridge. This result connects directly to the neutral theory 19.04.02 pending: the expected number of segregating sites in a sample is $E[S]=\mu \cdot E[L_{\text{total}}]=\theta H_{n-1}$, which is the basis for Watterson's estimator ${\hat{\theta }}_W=S/H_{n-1}$. The separation of genealogy (coalescent tree) from mutation (Poisson process along branches) is the central insight that makes coalescent-based inference possible.

Exercises Intermediate+

Coalescent-based demographic inference and linkage disequilibrium in practice Master

Skyline plots and PSMC

The skyline plot (Pybus et al., 2000) estimates changes in effective population size over time from a gene genealogy. Given a coalescent tree with coalescence times $t_2,t_3,\dots ,t_n$, the maximum-likelihood estimate of $N_e$ during the interval between coalescences $k$ and $k-1$ is:

$${\hat{N}}_e^{(k)}=\frac{k(k-1)\cdot (t_k-t_{k+1})}{4},$$

in rescaled coalescent units. Stitching together these interval-wise estimates produces a step function tracing $N_e(t)$ through time -- the skyline plot. The classic skyline is piecewise constant with a step at each coalescence; the extended Bayesian skyline (EBSP; Heled and Drummond, 2008) smooths this using a prior on population-size change.

The pairwise sequentially Markovian coalescent (PSMC; Li and Durbin, 2011) extracts demographic history from a single diploid genome. The two haplotypes of a diploid individual are related by a coalescent process that varies along the genome due to recombination. At each position, the local TMRCA of the two haplotypes records the population size at that coalescence time. By scanning along the genome and computing the distribution of local coalescence times, PSMC reconstructs $N_e(t)$ over a time range from roughly $10^4$ to $10^6$ generations ago. The resolution is limited at the recent end (few coalescences in a large population produce short regions of low diversity) and at the ancient end (multiple hits obscure the true coalescence time).

Bayesian skyline and multi-locus coalescent

The Bayesian skyline plot (Drummond et al., 2005) places skyline estimation within a Bayesian framework, jointly inferring the genealogy and population-size trajectory from sequence data. The prior on $N_e(t)$ is a piecewise-constant model with a user-specified number of change points; the posterior distribution provides credible intervals on $N_e(t)$ at each time point. The method is implemented in BEAST (Bayesian Evolutionary Analysis Sampling Trees) and handles uncertainty in both the tree topology and coalescence times.

*BEAST (Heled and Drummond, 2010) extends the coalescent to multiple loci. Under the multispecies coalescent, each locus has its own gene tree, but the gene trees are correlated through a shared species tree. *BEAST jointly estimates the species tree, the gene trees at each locus, the population sizes along each branch of the species tree, and the divergence times -- all within a single Bayesian analysis. The multi-locus approach dramatically improves the precision of species-tree estimation relative to single-locus analyses because each locus provides an independent realisation of the coalescent process conditional on the species tree.

Species delimitation and gene-tree discordance

BPP (Bayesian Phylogenetics and Phylogeography; Yang and Rannala, 2010) uses the multispecies coalescent to delimit species from multi-locus sequence data. The method computes the posterior probability of different species-delimitation models (different assignments of individuals to species), integrating over uncertainty in gene trees and population sizes. BPP is widely used but has been criticised for potentially over-splitting: when population structure is present within a single species, BPP may assign subpopulations to different species with high posterior probability, because the multispecies coalescent attributes coalescence-time differences between subpopulations to separate species rather than to structure within a species.

Incomplete lineage sorting (ILS) is the failure of gene lineages to coalesce within the species branches of the species tree. When two speciation events occur in rapid succession, gene lineages from the three descendant species may coalesce in the ancestral population before the first speciation event. The probability that a random gene tree matches the species tree for three species with short internal branch length $T$ (in coalescent units) is:

$$P(\text{gene tree}=\text{species tree})=1-\frac{2}{3}{e}^{-T}.$$

When $T$ is small (rapid successive speciations), this probability can be well below 1, meaning a substantial fraction of loci have gene trees that disagree with the species tree. This gene tree -- species tree discordance is ubiquitous in rapidly radiating groups (e.g., Darwin's finches, cichlids, Drosophila) and requires the multispecies coalescent for correct inference.

The multispecies coalescent

The multispecies coalescent (Rannala and Yang, 2003) models the joint distribution of gene trees across multiple loci conditional on a species tree with population sizes along each branch. For a species tree with branch lengths ${\tau }_i$ (in coalescent units) and population-size parameters ${\theta }_i$, the probability of a gene tree $G$ at a single locus is:

$$P(G\mid \text{species tree})=\prod _{\text{branches}}P(\text{coalescence times on branch}\mid {\tau }_i,{\theta }_i).$$

On each branch of the species tree, the gene lineages entering that branch coalesce according to the Kingman process with rate scaled by ${\theta }_i$. If $k$ lineages enter a branch of length $\tau$ and $m$ exit (having coalesced to $m$ lineages), the contribution to the likelihood integrates over the coalescence times within the branch. The multispecies coalescent provides the theoretical basis for species-tree inference (ASTRAL, MP-EST, *BEAST), species delimitation (BPP), and testing hypotheses about the timing and order of speciation events.

The ancestral recombination graph

The ancestral recombination graph (ARG; Griffiths and Marjoram, 1997) extends the coalescent to multiple linked loci. In the ARG, each lineage carries a segment of the chromosome. Going backward in time, two types of events occur:

1. Coalescence: two lineages carrying overlapping ancestral material merge at rate $\left(\genfrac{}{}{0ex}{}{k}{2}\right)/(2N_e)$ per pair, same as the standard coalescent.
2. Recombination: a lineage splits at rate $\rho /L$ per generation (where $\rho =4N_{e}r$ is the population-scaled recombination rate and $L$ is the number of loci), with the ancestral material partitioned between the two daughter lineages at the recombination breakpoint.

The ARG is a directed acyclic graph rather than a tree: each recombination event introduces a new node with two parents (the two ancestral segments). The total number of nodes in the ARG grows with both the sample size and $\rho$, making exact computation intractable for whole-genome data. Modern methods (SMC, SMC'; McMillan and Durbin's PSMC; the SMC-based approximations in RELATE and tsinfer) approximate the ARG by assuming that at most one recombination event occurs between coalescences, reducing the graph to a Markov chain along the genome.

The ARG is the full genealogical record of a sample: every coalescence and every recombination event in the history of the sampled chromosomes. A tree sequence (Kelleher et al., 2016) is a compressed encoding of the ARG that stores the local genealogical tree at each position along the genome, sharing tree structure across adjacent positions where the genealogy has not changed. Tree sequences enable exact simulation and efficient inference of the ARG from large genomic datasets (implemented in msprime and tsinfer).

LD mapping and GWAS

Linkage disequilibrium mapping exploits the correlation between nearby loci to map disease genes. If a disease-causing variant exists at one position, it will be in LD with nearby markers. By genotyping a dense set of markers (SNPs) across the genome and testing each for association with a trait, a genome-wide association study (GWAS) identifies regions where marker-trait LD is stronger than expected by chance.

The statistical framework for GWAS is straightforward: for each SNP $j$, regress the phenotype $y$ on genotype $g_j$ (coded as 0, 1, or 2 copies of the minor allele):

$$y_i={\beta }_0+{\beta }_j{g}_{ij}+{ϵ}_i.$$

The test of $H_0:{\beta }_j=0$ produces a $p$-value for each SNP. The threshold for genome-wide significance is approximately $5\times 10^{-8}$, derived from a Bonferroni correction for roughly 1 million independent tests (the approximate number of independent LD blocks in the human genome).

Tagging SNPs

Not every SNP needs to be genotyped directly. A tagging SNP (or tag SNP) is a representative SNP chosen to capture the variation in a region of high LD. If SNP $j$ has $r^2>0.8$ with SNP $k$, then genotyping $j$ provides nearly equivalent statistical power for detecting association with a causal variant at $k$'s position. SNP arrays (e.g., Illumina Global Screening Array, Affymetrix Axiom) exploit this: they genotype 500K--2M carefully chosen tag SNPs that capture common variation across the genome through LD, achieving effective coverage of ~10M common SNPs.

The $r^2$ threshold for tagging depends on the application. For GWAS, $r^2>0.8$ is standard (losing at most 20% power). For imputation (predicting ungenotyped SNPs from a reference panel), the relevant measure is the imputation $r^2$ (squared correlation between true and imputed genotypes), which depends on the frequency of the target SNP and the LD structure of the reference panel.

Phased and unphased data

A haplotype is the combination of alleles at multiple loci on a single chromosome. Most sequencing data is unphased: the genotype at each locus is known (e.g., individual is heterozygous A/a), but the chromosomal assignment (which allele is on which chromosome copy) is not. Phased data assigns each allele to its chromosome of origin.

Phasing matters for LD analysis because LD operates on haplotypes, not genotypes. The measure $D$ and $r^2$ are defined in terms of haplotype frequencies $p_{AB}$, $p_{Ab}$, $p_{aB}$, $p_{ab}$, which require knowing the two haplotypes each individual carries. With unphased data, haplotype frequencies must be estimated, typically via the EM algorithm (Excoffier and Slatkin, 1995): initialise haplotype frequencies, compute expected genotype frequencies under Hardy-Weinberg, update haplotype frequencies from observed genotype counts, and iterate to convergence.

Statistical phasing methods infer haplotypes from unphased genotype data by exploiting LD. The most common approaches are:

• PHASE (Stephens et al., 2001): uses a Bayesian model that builds haplotypes by sampling from a prior distribution favouring haplotypes similar to those already constructed, reflecting the empirical observation that most haplotypes in a population are recombinants of a small number of ancestral haplotypes.
• fastPHASE (Scheet and Stephens, 2006): uses a hidden Markov model where the hidden states are cluster haplotypes and the observed data are genotypes. The HMM models the process of sampling haplotypes from a set of ancestral templates, with transitions corresponding to recombination between templates.
• SHAPEIT (Delaneau et al., 2012): uses a hidden Markov model on a reference panel of known haplotypes, conditioning on the haplotypes of all other individuals in the sample to improve phasing accuracy.
• Eagle (Loh et al., 2016): scales to biobank-sized datasets by using a faster HMM with a positional Burrows-Wheeler transform to efficiently search for matching haplotypes in a reference panel.

Phasing accuracy depends on the density of markers, the complexity of the LD structure, and the availability of a reference panel. For common variants ($\text{MAF}>5\%$), phasing accuracy with a good reference panel exceeds 99%. For rare variants, phasing is less accurate because fewer reference haplotypes carry the variant, reducing the information available for statistical inference.

Gene genealogies, the ARG, and coalescent estimation theory Master

Coalescent estimation of population parameters

The coalescent provides a likelihood framework for estimating population-genetic parameters. The probability of observing a particular genealogy $\mathcal{G}$ with coalescence times $\mathbf{t}=(t_2,t_3,\dots ,t_n)$ under the Kingman coalescent with effective size $N_e$ is:

$$P(\mathcal{G}\mid N_e)=\prod _{k=2}^n\frac{k(k-1)}{4N_e}\exp \left(-\frac{k(k-1)}{4N_e}{t}_k\right).$$

The log-likelihood is:

$$\ell (N_e)=\sum _{k=2}^n\left[\ln \left(\frac{k(k-1)}{4N_e}\right)-\frac{k(k-1)}{4N_e}{t}_k\right].$$

Setting $d\ell /d{N}_e=0$ yields the MLE:

$${\hat{N}}_e=\frac{{\sum }_{k=2}^{n}k(k-1)t_k}{4(n-1)}.$$

This is the coalescent-based estimate of $N_e$ from a known genealogy. In practice, the genealogy is unknown and must be inferred from sequence data, so the full likelihood integrates over genealogies weighted by their sequence likelihood (a computationally demanding problem solved by MCMC in BEAST).

Coalescent hidden Markov models

The sequentially Markovian coalescent (SMC; McVean and Cardin, 2005) approximates the ARG as a Markov process along the genome. At each position, the local genealogy is a coalescent tree. Moving along the genome, recombination events change the local genealogy, but under the SMC approximation, each recombination event modifies the tree by breaking one branch and allowing the displaced lineage to re-coalesce. The SMC' modification (Marjoram and Wall, 2006) allows the displaced lineage to re-coalesce below the break point, improving accuracy.

The SMC framework enables the PSMC and related methods:

• PSMC (Li and Durbin, 2011): applies the SMC to a single diploid genome, reading the local TMRCA along the genome as a proxy for $N_e(t)$.
• MSMC (Schiffels and Durbin, 2014): extends PSMC to multiple genomes, using the coalescence patterns among $k>2$ haplotypes to resolve more recent demographic events.
• SMC++ (Terhorst et al., 2017): scales to hundreds of genomes by using the site frequency spectrum in addition to pairwise coalescence times, extending the time range of inference.
• RELATE (Speidel et al., 2019): estimates genome-wide genealogies at every position, simultaneously inferring the ARG and demographic history.

The variance of coalescent-based estimators

Coalescent-based estimators have high variance because the coalescent process is inherently stochastic: a single genealogy is one realisation of a random process, and the coalescence times are exponentially distributed with high coefficient of variation. For the TMRCA of a sample of $n$ genes:

$$\text{Var}(T_{\text{MRCA}})=\sum _{k=2}^n\text{Var}(T_k)=\sum _{k=2}^n{\left(\frac{4N_e}{k(k-1)}\right)}^2.$$

The dominant contribution comes from $T_2$ (the last coalescence), for which $\text{Var}(T_2)=(2N_e{)}^2$. The standard deviation of the TMRCA is of the same order as its mean, making single-locus estimates unreliable. Multi-locus approaches (*BEAST, MSMC) reduce variance by averaging over independent genealogical realisations, but the reduction scales only as $1/\sqrt{L}$ where $L$ is the number of (effectively independent) loci.

This high variance has practical consequences: demographic inferences from single genes are unreliable, and even genome-wide estimates of $N_e(t)$ have wide credible intervals, particularly for ancient time periods where few coalescence events provide information.

Likelihood-free and ABC methods

When the likelihood is intractable (e.g., for complex demographic models with selection, structure, and admixture), approximate Bayesian computation (ABC) provides a likelihood-free alternative. The procedure is:

1. Sample parameters $\mathbf{\theta }$ from the prior $\pi (\mathbf{\theta })$.
2. Simulate data ${\mathbf{D}}_{\text{sim}}$ under the model with parameters $\mathbf{\theta }$.
3. Compute summary statistics $S({\mathbf{D}}_{\text{sim}})$ and compare to $S({\mathbf{D}}_{\text{obs}})$.
4. Accept parameters for which $\Vert S({\mathbf{D}}_{\text{sim}})-S({\mathbf{D}}_{\text{obs}})\Vert <ϵ$.

ABC methods for coalescent models use summary statistics such as the site frequency spectrum, haplotype diversity, LD decay curves, and the distribution of pairwise differences. The choice of summary statistics critically affects the quality of the posterior approximation: sufficient statistics retain all information, but in practice no finite set of summaries is sufficient for complex demographic models. Regression adjustment (Beaumont et al., 2002) and neural density estimation (Jiang et al., 2017) improve ABC accuracy by correcting for the mismatch between accepted and true posterior.

Full proof set Master

Proposition 1 (Distribution of coalescence times). Under the Kingman coalescent with $k$ lineages in a diploid population of effective size $N_e$, the waiting time $T_k$ until the next coalescence is exponentially distributed with rate $k(k-1)/(4N_e)$, and $T_2,T_3,\dots ,T_n$ are mutually independent.

Proof. In each generation, the probability that a specific pair of $k$ lineages shares a common parent is $1/(2N_e)$. The probability that at least one pair among $\left(\genfrac{}{}{0ex}{}{k}{2}\right)$ pairs coalesces is $p_k=\left(\genfrac{}{}{0ex}{}{k}{2}\right)/(2N_e)+O(1/N_e^2)$, neglecting the probability of two or more simultaneous coalescences (which is $O(1/N_e^2)$). The probability of no coalescence for $t$ consecutive generations is $(1-p_k{)}^t\approx e^{-p_{k}t}$ for large $N_e$. Rescaling time as $\tau =t/(2N_e)$ and taking $N_e\to \infty$, this becomes $e^{-\left(\genfrac{}{}{0ex}{}{k}{2}\right)\tau }$, the survival function of an exponential with rate $\left(\genfrac{}{}{0ex}{}{k}{2}\right)/(2N_e)=k(k-1)/(4N_e)$.

Independence of $T_2,T_3,\dots ,T_n$ follows from the Markov property: the process after each coalescence depends only on the number of remaining lineages, not on the timing of previous coalescences. $\square$

Proposition 2 (Expected TMRCA and total tree length). For a sample of $n$ genes, $E[T_{\text{MRCA}}]=4N_e(1-1/n)$ and $E[L_{\text{total}}]=4N_e{H}_{n-1}$.

Proof. As shown in the Key theorem section. $\square$

Proposition 3 (Variance of TMRCA). The variance of the TMRCA of a sample of $n$ genes is:

$$\text{Var}(T_{\text{MRCA}})=\sum _{k=2}^n{\left(\frac{4N_e}{k(k-1)}\right)}^2.$$

Proof. Since $T_2,T_3,\dots ,T_n$ are independent (Proposition 1):

$$\text{Var}(T_{\text{MRCA}})=\sum _{k=2}^n\text{Var}(T_k)=\sum _{k=2}^n{\left(\frac{1}{{\lambda }_k}\right)}^2=\sum _{k=2}^n{\left(\frac{4N_e}{k(k-1)}\right)}^2.$$

The dominant term is $\text{Var}(T_2)=(2N_e{)}^2=4N_e^2$. For large $n$, the variance converges to:

$$\text{Var}(T_{\text{MRCA}})\to 4N_e^2\sum _{k=2}^{\infty }\frac{16}{k^2(k-1{)}^2}\approx 4N_e^2\cdot \frac{{\pi }^2}{3}-12N_e^2\approx 1.29\times (4N_e{)}^2.$$

The standard deviation is roughly $0.8\times 4N_e$, comparable to the mean of $\approx 4N_e$, confirming the high stochasticity of the coalescent. $\square$

Proposition 4 (Expected LD at drift-recombination equilibrium). At equilibrium between drift and recombination in a diploid population of effective size $N_e$ with recombination rate $c$ between two biallelic loci, the expected squared correlation coefficient is $E[r^2]\approx 1/(1+\rho )$ where $\rho =4N_{e}c$.

Proof (outline). Track the recursion for $E[D^2]$ under drift and recombination. In one generation, recombination reduces $D$ by a factor $(1-c)$: $D_{t+1}=(1-c)D_t+{\delta }_t$, where ${\delta }_t$ is the change due to drift (sampling variance). The drift contribution has $E[{\delta }_t]=0$ and $\text{Var}({\delta }_t)=p_A{p}_a{p}_B{p}_b/(2N_e)$ (the multinomial sampling covariance).

At equilibrium, $E[D_{t+1}^2]=E[D_t^2]=E[D^2]$. Squaring the recursion:

$$E[D^2]=(1-c{)}^{2}E[D^2]+\frac{p_A{p}_a{p}_B{p}_b}{2N_e}.$$

Solving (approximating $(1-c{)}^2\approx 1-2c$ for small $c$):

$$E[D^2]\approx \frac{p_A{p}_a{p}_B{p}_b}{4N_{e}c}=\frac{p_A{p}_a{p}_B{p}_b}{\rho }.$$

Since $r^2=D^2/(p_A{p}_a{p}_B{p}_b)$:

$$E[r^2]\approx \frac{1}{\rho }=\frac{1}{4N_{e}c}.$$

The more precise approximation $E[r^2]\approx 1/(1+\rho )$ (Hill and Robertson, 1968) accounts for the coupling between allele frequencies and $D$ that the simple derivation neglects. $\square$

Proposition 5 (Probability of concordant gene tree and species tree under ILS). For three species (A, B, C) with species tree $((A,B),C)$ and internal branch length $T$ in coalescent units ($T=\tau /(2N_e)$ where $\tau$ is the branch length in generations), the probability that a random gene tree matches the species tree is:

$$P(\text{concordance})=1-\frac{2}{3}{e}^{-T}.$$

Proof. Consider the ancestral population of species A and B, with branch length $T$ in coalescent units. Two gene lineages from A and B enter this ancestral population. They coalesce with probability $1-e^{-T}$ within the branch (at rate 1 in coalescent units). If they coalesce within the branch (probability $1-e^{-T}$), the resulting lineage coalesces with the C lineage in the deeper ancestral population, and the gene tree necessarily matches $((A,B),C)$.

If A and B fail to coalesce within the branch (probability $e^{-T}$), all three lineages (A, B, C) enter the deeper ancestral population simultaneously. By symmetry, each of the three possible resolved topologies has equal probability $1/3$. The topology matching the species tree has probability $1/3$.

Therefore:

$$P(\text{concordance})=(1-e^{-T})\times 1+e^{-T}\times \frac{1}{3}=1-e^{-T}+\frac{e^{-T}}{3}=1-\frac{2}{3}{e}^{-T}.\square$$

When $T=0$ (instantaneous speciation): $P=1/3$ (random, no information). When $T\to \infty$ (long internal branch): $P\to 1$ (no ILS). At $T=1$: $P=1-2/(3e)\approx 0.755$.

Connections Master

• Genetic drift 19.04.01. The coalescent is the retrospective counterpart of the forward-time Wright-Fisher model 19.04.01. The coalescence rate $\left(\genfrac{}{}{0ex}{}{k}{2}\right)/(2N_e)$ derives directly from the Wright-Fisher sampling process: two gene copies coalesce when they share a parent, which happens with probability $1/(2N_e)$ per generation. The coalescent provides the genealogical framework for all drift-based inference, and $N_e$ estimated from coalescence times equals the variance effective size from the forward model.

• Neutral theory 19.04.02 pending. Under neutrality, the genealogical process (coalescent) and the mutational process (Poisson along branches) are independent. This separation is what enables Watterson's estimator, Tajima's D, and all coalescent-based tests of neutrality 19.04.02 pending. Deviations from the standard coalescent (due to selection, population structure, or demography) alter the shape of the genealogy, which in turn alters the distribution of polymorphism -- the basis for detecting selection from sequence data.

• Phylogenetics 19.07.01. The multispecies coalescent bridges population genetics and phylogenetics: gene trees are nested within the species tree, and the coalescent provides the probability model for gene-tree -- species-tree discordance. Methods like ASTRAL and *BEAST use the multispecies coalescent to infer species trees from multi-locus data, accounting for ILS. The coalescent also provides the null model for testing whether discordance among gene trees exceeds the ILS expectation (indicating hybridisation or introgression).

• Hardy-Weinberg extensions 19.02.02 pending. LD is a multilocus extension of the Hardy-Weinberg framework. Where Hardy-Weinberg tests for random mating at a single locus (genotype frequencies = product of allele frequencies), LD tests for random association across loci (haplotype frequencies = product of marginal allele frequencies). Both are tests of independence, and both are violated by evolutionary forces (drift, selection, migration) that create non-random associations.

• Quantitative genetics 19.05.01. LD determines the accuracy of genomic prediction (genomic BLUP, genomic selection). The effective number of independent chromosome segments $M_e$ is approximately $2N_{e}L/\ln (4N_{e}L)$ where $L$ is the genome length in Morgans, and the accuracy of genomic prediction scales as $\sqrt{M_e/(M_e+\lambda )}$ where $\lambda$ is the ratio of residual to genetic variance. LD structure therefore directly controls the feasibility of marker-based breeding.

Historical & philosophical context Master

John Kingman (1982) derived the coalescent in three simultaneous publications, recognising that the genealogy of a sample could be described by a simple continuous-time Markov process. The key insight was that in the large-population limit, only pairwise coalescences matter (the $O(1/N_e^2)$ probability of triple coalescences vanishes), and the process simplifies to a random binary tree with exponentially distributed branch lengths. Kingman showed that this limit is universal: it applies not only to the Wright-Fisher model but to any Cannings model (any exchangeable reproduction scheme with finite offspring variance).

The coalescent perspective reversed the direction of population-genetic reasoning. Before Kingman, population genetics was primarily forward-looking: given allele frequencies today, what will they be tomorrow? After Kingman, the retrospective approach became dominant: given a sample of genes today, what is the distribution of their genealogy? This shift was not merely mathematical -- it reconceptualised what population genetics explains. The forward approach explains the dynamics of populations; the coalescent explains the pattern of genetic variation in samples. Since all empirical data are samples, the coalescent is in many ways the more natural framework.

The application of coalescent theory to LD and recombination was developed by Hudson (1983), who introduced the coalescent with recombination and showed that LD patterns encode information about $\rho =4N_{e}c$. Hudson's work established that the genome carries a record of its own genealogical history, with each region's genealogy correlated with its neighbours through the shared ARG. The practical exploitation of this record for demographic inference (PSMC, MSMC) and for association mapping (GWAS, tagging SNPs) has been one of the most productive applications of population-genetic theory in the genomic era.

The philosophical significance of coalescent theory is that it makes genealogy -- not population -- the primary object of study. In the coalescent framework, the population is a derived concept (the collection of all gene copies), while the genealogy of a sample is fundamental. This genealogical perspective has profoundly influenced how evolutionary biologists think about species: under the multispecies coalescent, a species is a branch of the species tree with an associated population size, and the boundaries between species are defined by the probability that gene lineages cross them without coalescing. The species concept becomes a statement about the genealogical process rather than about phenotypic distinctness or reproductive isolation.

Bibliography Master

@article{Kingman1982,
  author = {Kingman, J. F. C.},
  title = {The coalescent},
  journal = {Stochastic Processes and their Applications},
  volume = {13},
  year = {1982},
  pages = {235--248},
}

@article{Hudson1983,
  author = {Hudson, Richard R.},
  title = {Properties of a neutral allele model with intragenic recombination},
  journal = {Theoretical Population Biology},
  volume = {23},
  year = {1983},
  pages = {183--201},
}

@book{Wakeley2009,
  author = {Wakeley, John},
  title = {Coalescent Theory: An Introduction},
  publisher = {Roberts \& Company},
  year = {2009},
}

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

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

@article{PybusEtAl2000,
  author = {Pybus, Oliver G. and Rambaut, Andrew and Harvey, Paul H.},
  title = {An integrated framework for the inference of viral population history from reconstructed genealogies},
  journal = {Genetics},
  volume = {155},
  year = {2000},
  pages = {1429--1437},
}

@article{DrummondEtAl2005,
  author = {Drummond, Alexei J. and Rambaut, Andrew and Shapiro, Beth and Pybus, Oliver G.},
  title = {Bayesian coalescent inference of past population dynamics from multiple sequences},
  journal = {Molecular Biology and Evolution},
  volume = {22},
  year = {2005},
  pages = {1185--1192},
}

@article{LiDurbin2011,
  author = {Li, Heng and Durbin, Richard},
  title = {Inference of human population history from individual whole-genome sequences},
  journal = {Nature},
  volume = {475},
  year = {2011},
  pages = {493--496},
}

@article{SchiffelsDurbin2014,
  author = {Schiffels, Stephan and Durbin, Richard},
  title = {Inferring human population size and separation history from multiple genome sequences},
  journal = {Nature Genetics},
  volume = {46},
  year = {2014},
  pages = {919--925},
}

@article{HeledDrummond2010,
  author = {Heled, Joseph and Drummond, Alexei J.},
  title = {Bayesian inference of species trees from multilocus data},
  journal = {Molecular Biology and Evolution},
  volume = {27},
  year = {2010},
  pages = {570--580},
}

@article{YangRannala2010,
  author = {Yang, Ziheng and Rannala, Bruce},
  title = {Bayesian species delimitation},
  journal = {Proceedings of the National Academy of Sciences},
  volume = {107},
  year = {2010},
  pages = {15587--15592},
}

@article{RannalaYang2003,
  author = {Rannala, Bruce and Yang, Ziheng},
  title = {Bayes estimation of species divergence times and ancestral population sizes using {DNA} sequences from multiple loci},
  journal = {Genetics},
  volume = {164},
  year = {2003},
  pages = {1645--1656},
}

@article{GriffithsMarjoram1997,
  author = {Griffiths, Robert C. and Marjoram, Paul},
  title = {An ancestral recombination graph},
  journal = {Progress in Population Genetics and Human Evolution},
  volume = {87},
  year = {1997},
  pages = {257--270},
}

@article{McVeanCardin2005,
  author = {McVean, Gilean A. T. and Cardin, Niall J.},
  title = {Approximating the coalescent with recombination},
  journal = {Philosophical Transactions of the Royal Society B},
  volume = {360},
  year = {2005},
  pages = {1387--1393},
}

@article{HillRobertson1968,
  author = {Hill, William G. and Robertson, Alan},
  title = {Linkage disequilibrium in finite populations},
  journal = {Theoretical and Applied Genetics},
  volume = {38},
  year = {1968},
  pages = {226--231},
}

@article{KelleherEtAl2016,
  author = {Kelleher, Jerome and Etheridge, Alison M. and McVean, Gilean},
  title = {Efficient coalescent simulation and genealogical analysis for large sample sizes},
  journal = {PLoS Computational Biology},
  volume = {12},
  year = {2016},
  pages = {e1004842},
}

@article{StephensEtAl2001,
  author = {Stephens, Matthew and Smith, Nicholas J. and Donnelly, Peter},
  title = {A new statistical method for haplotype reconstruction from population data},
  journal = {American Journal of Human Genetics},
  volume = {68},
  year = {2001},
  pages = {978--989},
}

@article{SpeidelEtAl2019,
  author = {Speidel, Leo and Forest, Mathieu and Shi, Sinan and Myers, Simon R.},
  title = {A method for genome-wide genealogy estimation for thousands of samples},
  journal = {Nature Genetics},
  volume = {51},
  year = {2019},
  pages = {1321--1329},
}

@article{TerhorstEtAl2017,
  author = {Terhorst, Jonathan and Kamm, John A. and Song, Yun S.},
  title = {Robust and scalable inference of population history from hundreds of unphased whole genomes},
  journal = {Nature Genetics},
  volume = {49},
  year = {2017},
  pages = {303--309},
}