Summary

• We can estimate phylogenetic trees from DNA sequence data using models of sequence evolution.
• Model-based phylogenies have likelihoods.
• Likelihoods can be used to test multiple competing hypotheses.

Learning objectives

• Understand biological sources of discordance between gene trees and species trees.
• Recognize cytonuclear discordance in plant phylogenetics.
• Be familiar with some species delimitation techniques.
• Apply the basic test of historical introgression (D-statistic / ABBA-BABA) between species.

Gene trees and species trees

Each locus has its own genealogy; these gene trees can be different from the species tree.

Three biological sources of gene tree discordance:

1. Incomplete lineage sorting (ILS) - example: rapid speciation in Andean tomato relatives Pease et al. (2016)
2. Gene duplication and loss - example: polyploidy leads to different trees due to differential homeolog retention Jiao et al. (2011)
3. Hybridization / introgression - example: hybrid sunflower species in extreme habitats Rieseberg et al. (2003)

Cytonuclear discordance

Cytonuclear discordance is common in plants

Chloroplasts can move between species through hybridization. Sometimes there is little evidence of the historical gene flow on the nuclear genome, but the chloroplast remains captured.

Chloroplast capture among annual sunflowers (Helianthus). Fig. 2 from Lee-Yaw et al. (2019)

When concatenated ML fails

For a four-taxon species tree \((((A,B)x,C)y,D)\), the most probable gene tree is not the species tree when \(x\) and \(y\) are very short.

The anomaly zone was given by Degnan and Rosenberg (2006) as:

\[ y < a(x) = \log\!\left[\frac23 + \frac{3e^{2x}-2}{18\,(e^{3x}-e^{2x})}\right] \]

For example, at \(x=y=0.1\): \(P[(((A,B),C),D)] = 0.10\), but \(P[((A,B),(C,D))] = 0.13\).

When concatenated ML fails

The anomaly zone boundary a(x) plotted against x. Based on analytical results from Degnan and Rosenberg (2006).

The coalescent for a single population

\[ P(\text{Two sequences do not coalesce in i generations}) \approx 1 - e^{-\frac{-i}{2N}} \approx e^{-T} \]

\(T\) is our coalescent unit and \(T = \frac{\text{# of generations}}{2N}\).

A Wright-Fisher process for 10 individuals.

Deep coalescence

ILS: How much discordance does the coalescent predict?

For a rooted triple with internal branch length \(T\), the probability that a gene tree matches the species tree is

\[ P(\text{gene tree matches species tree}) = 1 - \tfrac{2}{3} e^{-T} \]

and each of the two discordant gene trees occur with probability

\[ P(\text{gene tree is discordant}) = \frac{1}{3} e^{-T} \]

The coalescent predicts that the discordant gene trees are equally probable Hudson (1983).

Deep coalescence

ILS is reduced as the internal branch length \(T\) increases.

Predictions of gene tree discordance. Fig 2. of Tiley et al. (2020)

Four taxa, three site patterns

• Take the species tree \((((P1,P2),P3),O)\). The outgroup carries the ancestral allele (A); a mutation on an internal branch makes the tips below it derived (B) (Green et al. 2010).
• BBAA comes from the matching gene tree. ABBA and BABA need a discordant one.
• Under ILS the two discordant trees are equally likely, so ABBA and BABA SNPs should be equally common — the same symmetry as the three gene-tree probabilities.

Patterson’s \(D\)

\[ D = \frac{n_{\text{ABBA}} - n_{\text{BABA}}}{n_{\text{ABBA}} + n_{\text{BABA}}} \]

• \(D = 0\): ILS alone. \(D>0\): excess ABBA, so P3 shares alleles with P2. \(D<0\): excess BABA, so P3 shares alleles with P1 (Green et al. 2010; Durand et al. 2011).
• Needs a biallelic SNP at every site, with one allele per taxon (or a frequency per population); sites with any other pattern (BBAA, AABA, …) are ignored.
• Assumes the outgroup gives the ancestral state and the same substitution rate in P1 and P2.
• A test of the null (no gene flow), not an estimate of how much gene flow.

Worked example: counting patterns

Worked example: genome-wide \(D\)

Suppose a whole-genome SNP set gives \(n_{\text{ABBA}} = 1{,}250\) and \(n_{\text{BABA}} = 1{,}000\):

\[ D = \frac{1250 - 1000}{1250 + 1000} = \frac{250}{2250} = 0.111 \]

• Is 0.111 different from 0? A binomial test (ABBA vs BABA = 50:50) treats every SNP as independent: \(Z = 250/\sqrt{2250} = 5.3\).
• But neighbouring SNPs are linked and share a genealogy, so the effective sample size is much smaller and that \(Z\) is anti-conservative.
• Fix: resample the genome in blocks, not the SNPs — the block jackknife (Green et al. 2010; Patterson et al. 2012).

Worked example: block jackknife

The limitations of \(D\)

• Significant \(D\) ≠ proof of recent hybridization. Ancestral population structure in the P3 + P2 ancestor can also give ABBA excess (Durand et al. 2011).
• Direction is ambiguous: \(D>0\) can be P3→P2 or P2→P3, and gene flow from a ghost lineage related to P3 gives the same signal.
• Does not estimate the amount of gene flow \(f_4\)-ratio, \(f_d\) and \(f_{\text{dM}}\) estimate the introgressed fraction, and work better in genomic windows than \(D\) (Martin, Davey, and Jiggins 2015)

In plants, these statistics are widely used on target-capture or other reduced representation data. We have to be mindful of ploidy and duplications though!

The MSC: a censored coalescent

The gene-tree density, one population at a time

For a population of duration \(\tau\) that \(m\) lineages enter and \(n\) leave (Rannala and Yang 2003, eq. 8):

\[ \underbrace{\prod_{j=n+1}^{m}\frac{2}{\theta}\exp\left\{-\frac{j(j-1)}{\theta}\,t_j\right\}}_{\text{each coalescence, and which pair}} \;\times\; \underbrace{\exp\left\{-\frac{n(n-1)}{\theta}\Big[\tau - \textstyle\sum_j t_j\Big]\right\}}_{\text{no more coalescences before } \tau} \]

• \(t_j\) = waiting time while \(j\) lineages remain; the total rate is \(\binom{j}{2}\cdot\frac{2}{\theta} = \frac{j(j-1)}{\theta}\).
• This is the popgen result \(f(T_j)=\binom{j}{2}e^{-\binom{j}{2}T_j}\) with time rescaled, \(t = T\theta/2\), then censored at \(\tau\).
• \(f(G,t\mid\Theta)\) for the whole locus = the product over populations.

The MSC likelihood

\[ \ell(\Theta) = \sum_{i=1}^{L}\log f(X_i\mid\Theta) = \sum_{i=1}^{L}\log\left\{\sum_{G_i}\int_{t_i} \bbox[#DCEBFA,3px]{f(G_i,t_i\mid\Theta)}\; \bbox[#DDF3EA,3px]{f(X_i\mid G_i,t_i)}\,dt_i\right\} \]

• \(\bbox[#DCEBFA,3px]{\text{blue}}\): the coalescent density of gene tree \(G_i\) and its times — the previous slide.
• \(\bbox[#DDF3EA,3px]{\text{green}}\): Felsenstein’s likelihood of locus \(i\)’s alignment on that gene tree — the first half of this deck.
• Unknown gene trees are summed and integrated over, just as unknown ancestral states were.

Full-likelihood methods

\[ \ell(\Theta) = \sum_{i=1}^{L}\log\left\{\sum_{G_i}\int_{t_i} \bbox[#DCEBFA,3px]{f(G_i,t_i\mid\Theta)}\; \bbox[#DDF3EA,3px]{f(X_i\mid G_i,t_i)}\,dt_i\right\} \]

• You supply only the sequences; MCMC jointly samples the gene trees and the species tree — so this is only practical with Bayesian methods.
• Can be very slow, or impossible, for many loci and moderate numbers of tips.
• Software: StarBEAST2 (Ogilvie, Bouckaert, and Drummond 2017), BPP (Rannala and Yang 2017; Flouri et al. 2018).

Summary methods: ASTRAL

\[ \ell(\Theta) \approx \sum_{i=1}^{L}\log\left\{\sum_{G_i}\int_{t_i} \bbox[#DCEBFA,3px]{f(G_i,t_i\mid\Theta)}\; \cancel{f(X_i\mid G_i,t_i)}\,dt_i\right\} \]

• Estimate each gene tree first (e.g. IQ-TREE), then treat it as known data.
• ASTRAL finds the species tree sharing the most quartets with the gene trees: consistent under the MSC — if the gene trees are right (Zhang et al. 2018).

Summary methods: ASTRAL

Large species tree estimations with ASTRAL. Fig 2. of One Thousand Plant Transcriptomes Initiative (2019)

When the gene trees are wrong

Error in gene-tree estimation can mislead summary methods, but not full-likelihood methods.

Simulated species trees estimation with summary (MP-EST) and full-likelihood (BPP) methods. From Xu and Yang (2016)

How much error is in empirical data?

• A non-trivial example: In Malpighiales, gene-tree estimation error explains more gene-tree variation (52%) than gene flow (32%) or ILS (15%) (Cai et al. 2021).

Single-locus delimitation: barcoding

Before more easily obtainable data, barcoding studies were a common approach for species identification.

• One standard locus that could be easily amplified from anything was sequenced, such as rbcL or matK for land plants (CBOL Plant Working Group 2009).

Conceptual illustration of the barcoding gap

Testing species hypotheses with the MSC

Using the MSC and model selection techniques to accept or reject species hypotheses.

\[\mathrm{LRT} = 2\,(\ln L_{1} - \ln L_{0})\]

Bayesian approaches: delimiting species on a guide tree

• Each delimitation is a different MSC model; reversible-jump MCMC moves among them (Yang and Rannala 2010).
• Output: the posterior probability of each delimitation, and of the number of species.

Using reverse jump moves to estimate posterior probabilities for species hypotheses.

Bayes-factor delimitation

\[ \mathrm{BF}_{10} = \frac{p(D\mid M_1)}{p(D\mid M_0)}, \qquad p(D\mid M) = \int \underbrace{p(D\mid\theta, M)}_{\text{likelihood}}\;\underbrace{p(\theta\mid M)}_{\text{prior}}\,d\theta \]

• \(p(D\mid M)\) is the marginal likelihood: the likelihood averaged over the prior.
• It can’t be read off an MCMC run. It needs a sampling technique across multiple MCMC runs to remove effects of the priors (Xie et al. 2011).
• Feasible for small samples, similar to how it was applied to sub-clades across Nastanthus–Gamocarpha (Denham et al. 2019).