The multifurcating skyline plot

Abstract

A variety of methods based on coalescent theory have been developed to infer demographic history from gene sequences sampled from natural populations. The ‘skyline plot’ and related approaches are commonly employed as flexible prior distributions for phylogenetic trees in the Bayesian analysis of pathogen gene sequences. In this work we extend the classic and generalized skyline plot methods to phylogenies that contain one or more multifurcations (i.e. hard polytomies). We use the theory of Λ-coalescents (specifically, Beta(2-α,α)-coalescents) to develop the ‘multifurcating skyline plot’, which estimates a piecewise constant function of effective population size through time, conditional on a time-scaled multifurcating phylogeny. We implement a smoothing procedure and extend the method to serially sampled (heterochronous) data, but we do not address here the problem of estimating trees with multifurcations from gene sequence alignments. We validate our estimator on simulated data using maximum likelihood and find that parameters of the Beta(2-α,α) -coalescent process can be estimated accurately. Furthermore, we apply the multifurcating skyline plot to simulated trees generated by tracking transmissions in an individual-based model of epidemic superspreading. We find that high levels of superspreading are consistent with the high-variance assumptions underlying Λ-coalescents and that the estimated parameters of the Λ-coalescent model contain information about the degree of superspreading.

Multiple Merger Genealogies in Outbreaks ofMycobacterium tuberculosis

“…The Beta coalescent, in which the probability of each individual to coalesce in a multiple merger event is regulated by a Beta distribution with parameters α (between 1 and 2) and 2 − α. The Beta coalescent was originally introduced to model populations with sweepstakes reproduction ( Schweinsberg 2003 ), but it was also proposed to capture the genealogies of populations undergoing recurrent bottlenecks and of epidemics characterized by superspreaders ( Tellier and Lemaire 2014 ; Hoscheit and Pybus 2019 ). Lower values of α (closer to one) correspond to larger multiple merger events, and for α = 1 the Beta coalescent corresponds to the Bolthausen–Sznitman (BSZ) coalescent ( Bolthausen and Sznitman 1998 ).…”

“…Thus, to overcome the limitation of assuming synchronous sampling, we encourage future studies to develop MMC models that explicitly consider the time of sampling. Such a model is proposed in Hoscheit and Pybus ( 2019 ) as an extension of the Beta coalescent but without an explicit mechanism to convert real time units into coalescent time units.…”

“…The Beta coalescent has also been proposed to model the genealogies of epidemics with superspreaders. A simulation study showed that the parameter α contains information about the degree of superspreading, with lower values of α corresponding to higher level of superspreading (and larger multiple merger events; Hoscheit and Pybus 2019 ). However, also in this case, an explicit population model is missing.…”

“…Following Hoscheit and Pybus (2019) , we add serial sampling to the MMC and to Kingman’s coalescent with exponential growth simply by stopping the coalescent at times (on the coalescent time scale) where further individuals are sampled. Then, we start a new (independent) coalescent tree that has rates and waiting times as the nonserial coalescent (multiple merger or with growth) started in the last state of the stopped coalescent plus adding one block with a single individual for each individual sampled at this time.…”