Multiple merger coalescent inference of effective population size

Abstract

Variation in a sample of molecular sequence data informs about the past evolutionary history of the sample's population. Traditionally, Bayesian modelling coupled with the standard coalescent is used to infer the sample's bifurcating genealogy and demographic and evolutionary parameters such as effective population size and mutation rates. However, there are many situations where binary coalescent models do not accurately reflect the true underlying ancestral processes. Here, we propose a Bayesian non-parametric method for inferring effective population size trajectories from a multifurcating genealogy under the [Formula: see text]-coalescent. In particular, we jointly estimate the effective population size and the model parameter for the Beta-coalescent model, a special type of [Formula: see text]-coalescent. Finally, we test our methods on simulations and apply them to study various viral dynamics as well as Japanese sardine population size changes over time. The code and vignettes can be found in the phylodyn package.This article is part of the theme issue '"A mathematical theory of evolution": phylogenetic models dating back 100 years'.

Keywords: Beta-coalescent; Gaussian processes; Lambda-coalescent; Multiple mergers coalescent.

Figure 1.

The effect of $\alpha$ in tree topology and coalescent times. (a) Average block size for different values of $\alpha$ and number of tips. (b) Total coalescent rates ${\lambda }_b={\sum }_{k=2}^b\left(\genfrac{}{}{0px}{}{b}{k}\right){\lambda }_{b,k}$ when there are $b$ lineages for different values of $\alpha$, scaled by $\left(\genfrac{}{}{0px}{}{b}{2}\right)$, the rate under Kingman's coalescent.

Figure 2.

A multifurcating tree with eight total lineages and six coalescent events labelled with the coalescent times ${\displaystyle \mathit{t}}$, sampling times ${\displaystyle \mathit{s}}$, block sizes ${\displaystyle \mathit{m}}$ and extant lineages $A(t)$. Here, there are $L=3$ sampling times.

Figure 3.

Comparison of different methods estimating exponential $N_e(t)$. (a) A simulated genealogy under the Beta-coalescent with $\alpha =1.5$ and the exponential growth $N_e(t)$ trajectory. The genealogy has 39 internal nodes and 50 tips. (b) The effective population size trajectories reconstructed by different methods. Note the y-axis is plotted in log-scale. Solid lines are the median trajectories, and dotted lines are the 95% credible bands. The true trajectory is plotted in purple, and the solid black line and shaded band are the median and 95% credible band of the inferred $N_e(t)$ using the true $\alpha =1.5$.

Figure 4.

(a) A dated genealogy generated using 266 RSV-A sequences from Nextstrain. (b) The reconstructed effective population size trajectories for RSV-A sequences in south and southeast Asia. Note the $y$-axis is plotted in log-scale. Solid lines are the median trajectory, and dotted lines are the 95% credible band.

Figure 5.

(a) A dated phylogeny generated using 171 Enterovirus D68 sequences from Nextstrain. (b) The reconstructed effective population size trajectories for Enterovirus D68 sequences in Europe. Note the $y$-axis is plotted in log-scale. Solid lines are the median trajectory, and dotted lines are the 95% credible band.

Figure 6.

(a) A dated phylogeny generated using all Sardinops melanostictus sequences. There are 156 tips and 95 internal nodes. (b) The subtree constructed of only sequences sampled in 1990: 106 tips and 61 internal nodes.

Figure 7.

The reconstructed effective population size trajectories for (a) the 106 sequences sampled in 1990 and (b) all Sardinops melanostictus sequences. Note the $y$-axis is plotted in log-scale. Solid lines are the median trajectory, and dotted lines are the 95% credible band.

Figure 8.

(a) The serially sampled UPGMA binary genealogy of all Japanese sardine sequences. (b) Comparison of the reconstructed $N_e(t)$ using BNPR for the binary UPGMA and the MCMC $N_e(t)$ for the multifurcating tree in figure 6a. Note that the $T_{MRCA}$ for the multifurcating tree in figure 6a is around 29 000 years versus 23 000 years for the binary genealogy shown in (a), which is why the inferred $N_e(t)$ profiles have different ranges for $t$.

Figure 9.

An example of inference of effective population size of a phylogenetic tree using BNPR under various coalescent models.

Figure 10.

Log total coalescent rate when there are $b$ lineages. Solid lines show the true value for different values of $b$ and dashed lines show our approximation.

Figure 11.

Posterior distribution of $\alpha$ for the tree in figure 3. The blue and red lines are the mean and median values, respectively.

Figure 12.

Trace plots for the MCMC chain for the tree in figure 3. The marginal trace plot for $\log N_e(t)$ is taken at $t\approx 5.23$, which is the 70th grid point.

Figure 13.

Boxplots of inferred $\alpha$ values from simulations with trees generated from (a): uniform trajectory, (b): exponential trajectory, (c): boombust trajectory. Each row is a fixed $N$; each column is a fixed $\alpha$; and each colour is a sampling time. The black line is the true $\alpha$.

Figure 14.

The three ’true’ genealogies used for simulation. Each was generated under a Beta-coalescent with $\alpha =1.5$ and an exponential effective population size trajectory.

Figure 15.

Results of inferred $\alpha$ and $N_e(t)$ values under the block-size MLE, hybrid and MCMC methods based on reconstructed genealogies. (a) Scatterplots of inferred $\alpha$ values: the horizontal line is the true $\alpha =1.5$. (b) Scatterplots of the coverage of the 95% credible intervals of the inferred $N_e(t)$ values. (c) Scatterplots of the ${\log }_{10}$ MSE of the inferred $N_e(t)$ values. Note the red triangles in each figure represent the value from applying our methods on the true genealogy in figure 14. The green and blue colours represent the fast mutation scenario, and the original mutation scenario, respectively.