Comparison of Bayesian Coalescent Skyline Plot Models for Inferring Demographic Histories

Abstract Bayesian coalescent skyline plot models are widely used to infer demographic histories. The first (non-Bayesian) coalescent skyline plot model assumed a known genealogy as data, while subsequent models and implementations jointly inferred the genealogy and demographic history from sequence data, including heterochronous samples. Overall, there exist multiple different Bayesian coalescent skyline plot models which mainly differ in two key aspects: (i) how changes in population size are modeled through independent or autocorrelated prior distributions, and (ii) how many change-points in the demographic history are used, where they occur and if the number is pre-specified or inferred. The specific impact of each of these choices on the inferred demographic history is not known because of two reasons: first, not all models are implemented in the same software, and second, each model implementation makes specific choices that the biologist cannot influence. To facilitate a detailed evaluation of Bayesian coalescent skyline plot models, we implemented all currently described models in a flexible design into the software RevBayes. Furthermore, we evaluated models and choices on an empirical dataset of horses supplemented by a small simulation study. We find that estimated demographic histories can be grouped broadly into two groups depending on how change-points in the demographic history are specified (either independent of or at coalescent events). Our simulations suggest that models using change-points at coalescent events produce spurious variation near the present, while most models using independent change-points tend to over-smooth the inferred demographic history.


S1 Marginal likelihoods of the models
We performed marginal likelihood estimation using stepping-stone sampling [9] using the parallel power posterior sampler in RevBayes [5].Marginal likelihoods were estimated using the same models and data settings as described in the main manuscript.The power posterior MCMC analyses were run for 128 stones with 10,000 MCMC iterations as burnin and then 2,000 iterations for each consecutive stone.

S2 Sequence based analyses
We performed demographic analyses of horse data [8] with nine different models: a Constant model, a Skyline model, a Bayesian Skyline Plot (BSP ) model [1], an Extended BSP (EBSP ) model [4], a Skyride model [6], a Skygrid model [3], a Gaussian Markov Random Field (GMRF ) model [2], a Horseshoe Markov Random Field (HSMRF ) model [2], and our new Skyfish model (similar to [7]).Here, we show the resulting population size trajectories on a timeline from 0 to 1, 200, 000 years ago and the convergence assessment comparing two replicates of each analysis.The convergence assessment spans the time of the original analyses (0 − 500, 000 years ago for the isochronous samples, 0 − 1, 200, 000 years ago for the heterochronous samples).

HSMRF
Figure S1: Population size trajectories estimated with nine different models from isochronous sequence data on a timeline from 0 to 1, 200, 000 years ago.Genealogies and population size were jointly estimated.All analyses considered possible variation in population size between 0 and 500, 000 years ago.The bold line represents the median of the posterior distribution of the population size and the shaded area shows the 95% credible intervals.The models highlighted in violet (BSP, EBSP, Skyline, Skyride) are coalescent event based models, i.e., a change in population size can only occur at the time of a coalescent event.In models highlighted in green, population size changes can happen at specified times, independent from coalescent events.Here, all intervals for these models are equally-sized.In the orange Skyfish model, the number of intervals as well as their duration are estimated.All models except for the Constant model, the EBSP model, and the Skyline model have correlated intervals.MCMCs were run for 100, 000 iterations, sampling every tenth iteration, with a burn-in of 10% and two replicates, yielding 18, 000 samples in total.For plotting, the resulting trajectories were evaluated at 500 exponentially-spaced grid points between 0 and 1, 200, 000 years ago.

HSMRF
Figure S2: Population size trajectories estimated with nine different models from heterochronous sequence data on a timeline from 0 to 1, 200, 000 years ago.Genealogies and population size were jointly estimated.All analyses considered possible variation in population size between 0 and 1, 200, 000 years ago.The bold line represents the median of the posterior distribution of the population size and the shaded area shows the 95% credible intervals.The models highlighted in violet (BSP, EBSP, Skyline, Skyride) are coalescent event based models, i.e., a change in population size can only occur at the time of a coalescent event.In models highlighted in green, population size changes can happen at specified times, independent from coalescent events.Here, all intervals for these models are equally-sized.In the orange Skyfish model, the number of intervals as well as their duration are estimated.All models except for the Constant model, the EBSP model, and the Skyline model have correlated intervals.MCMCs were run for 100, 000 iterations, sampling every tenth iteration, with a burn-in of 10% and two replicates, yielding 18, 000 samples in total.For plotting, the resulting trajectories were evaluated at 500 exponentially-spaced grid points between 0 and 1, 200, 000 years ago.

S3 Simulation Study
S3.1 Simulations using the results of the BSP analysis to simulate

Figure S3 :Figure S4 :
FigureS3: Convergence assessment of analyses with isochronous sequence data.Kolmogorov-Smirnov test statistic (D) was calculated for the posterior distributions of two independent MCMC runs at 500 exponentially spaced grid points.The solid blue line depicts the threshold of 0.0921.

Figure S6 :Figure S7 :Figure S8 :
FigureS6: Simulation study results, simulated under the population size trajectory from the Constant analysis.10 simulations of sequence data with 36 tips were performed under the resulting population size trajectory from the Constant analysis with isochronous data.Analyses of each of the simulations were performed with the BSP model (top two rows) and the Skyfish model (bottom two rows).The bold green line represents the median of the posterior distribution of the population size and the shaded area shows the 95% credible intervals.The true trajectory used for simulations is depicted in orange.

Figure S9 :Figure S10 :Figure S11 :Figure S12 :
Figure S9: Comparing constant and linearly changing population size within intervals with the Skyfish model from heterochronous data.The bold line represents the median of the posterior distribution of the population size and the shaded area shows the 95% credible intervals.

Figure S14 :Figure S15 :Figure S17 :Figure S18 :
Figure S14: Kolmogorov-Smirnov test statistic comparing the posterior distributions of BSP analyses with BEAST1 and RevBayes at 500 equally spaced grid points.The solid blue line depicts the threshold of 0.0921.

Table S1 :
Overview of log-transformed marginal likelihoods estimated with a stepping-stone algorithm.Highest support is highlighted with bold font.
Simulation study results, simulated under the population size trajectory from the BSP analysis.simulations of sequence data with 36 tips were performed under the resulting population size trajectory from the BSP analysis with isochronous data.Analyses of each of the simulations were performed with the BSP model (top two rows) and the Skyfish model (bottom two rows).The bold green line represents the median of the posterior distribution of the population size and the shaded area shows the 95% credible intervals.The true trajectory used for simulations is depicted in orange.S3.2 Simulations using the results of the Constant analysis to simulate