Demographic inference through approximate-Bayesian-computation skyline plots

The skyline plot is a graphical representation of historical effective population sizes as a function of time. Past population sizes for these plots are estimated from genetic data, without a priori assumptions on the mathematical function defining the shape of the demographic trajectory. Because of this flexibility in shape, skyline plots can, in principle, provide realistic descriptions of the complex demographic scenarios that occur in natural populations. Currently, demographic estimates needed for skyline plots are estimated using coalescent samplers or a composite likelihood approach. Here, we provide a way to estimate historical effective population sizes using an Approximate Bayesian Computation (ABC) framework. We assess its performance using simulated and actual microsatellite datasets. Our method correctly retrieves the signal of contracting, constant and expanding populations, although the graphical shape of the plot is not always an accurate representation of the true demographic trajectory, particularly for recent changes in size and contracting populations. Because of the flexibility of ABC, similar approaches can be extended to other types of data, to multiple populations, or to other parameters that can change through time, such as the migration rate.

S1 Supplementary Methods S1.1 Approximating θ(t) from parameters θ 0 , θ 1 , ..., θ n In order to work with a set of common demographic parameters that will represent the demographic trajectory θ(t), the parameters of the constant piecewise demography (θ 0 , θ 1 , ..., θ n ) are transformed to parameters {θ(t 1 ), θ(t 2 ), ..., θ(t m )} as described in the main text. A graphical representation of this transformation is provided in Supplementary Figure S1, which shows the equivalence between the original and the derived parameters for three example simulations. Figure S1. Approximation of θ(t) from original parameters. Three simulations of demographic histories are represented, with one, none and two population size changes, respectively. Simulation 1 Simulation 2 Simulation 3 Approximation of θ(t) sim θ(t = 0) θ(t = 1) θ(t = 2) θ(t = 3) θ(t = 4) 1 4 4 40 40 40 2 1.5 1.5 1.5 1.5 1.5 3 0.04 11 11 178 178 S1.2 Alternative ABC skyline plot In addition to the procedure described in the main text, other ways to obtain skyline plots in approximate Bayesian Computation (ABC) can be devised, particularly by changing the number and type of models and the prior probabilities of parameters. For comparison, we present one of these alternatives. In this ap-proach we consider a single piecewise model with a single demographic change.
The parameters of the model are θ 0 , θ 1 and τ 1 , following the nomenclature presented in the main text. Demographic inference is performed on the derived parameters θ(t j ) in the same way as described in the main text. Evidence for population size change will be considered by evaluating the ratio θ 0 /θ 1 in a similar way as in MIGRAINE. This alternative approach was evaluated on the three example scenarios (contraction with θ 0 = 0.4, θ 1 = 40, τ = 0.1, expansion with θ 0 = 40, θ 1 = 0.4, τ = 0.1 and constant size with θ = 40; mutational model with P GSM = 0.22). Simulated data sets were generated for 50 individuals genotyped at 30 microsatellite loci.

S1.3 Comparison with alternative methods
We present a comparison of results between methods implemented in four programmes: DIYABCskylineplot (Navascués, 2017), MIGRAINE (Leblois et al., 2014), VarEff (Nikolic and Chevalet, 2014) and BEAST (Drummond et al., 2012), which can make inferences of demographic change from microsatellites following a generalized stepwise mutation model (GSM). This is not intended to be a thorough cross comparison, but to provide some intuition about how different their results can be, their limitations and, thus, the potential interest of their simultaneous use. We knowingly exclude MSVAR (Beaumont, 1999;Storz and Beaumont, 2002) from this comparison because it only implements a strict stepwise mutation model (SMM). We also excluded methods based on summary statistics (Cornuet and Luikart, 1996;Kimmel et al., 1998;Garza and Williamson, 2001) as we considered them to be superseded by ABC inference. However, Girod et al. (2011) evaluate those methods on equivalent simulated scenarios. A more thorough evaluation of the methods implemented in MIGRAINE and VarEff under other scenarios can be found in the original publications (Leblois et al., 2014;Nikolic and Chevalet, 2014). To our knowledge, there is no work addressing the performance of the extended Bayesian skyline plot method included in BEAST on microsatellite data. All four methods have been used to analyse seven data sets, three simulations (PODs, pseudo-observed data set) and four real data sets presented in the main text (leatherback turtle, whale shark, black-and-white colobus and red colobus). The three PODs that we evaluated correspond to the example scenarios used through this work (contraction with θ 0 = 0.4, θ 1 = 40, τ = 0.1, expansion with θ 0 = 40, θ 1 = 0.4, τ = 0.1 and constant size with θ = 40; with a generalised stepwise mutational model (GSM) with P GSM = 0.22). Simulated data sets were generated for 50 individuals genotyped at 10 microsatellite loci. The reduction of the number of loci respect to the simulations presented in the main text was decided to reduce the computing time of the analysis with BEAST.
Models and assumption of the methods implemented in DIYABCskylineplot and MIGRAINE are described in the main text. The analyses of the four real data case studies correspond to those presented in the main text. For the analysis of the three PODs, DIYABCskylineplot was run with the same priors described in the main text; and MIGRAINE was run using 2000 trees, 400 points at each iteration and a total of 12 iterations.
We used BEAST v1.8.0 (Drummond et al., 2012) to perform an extended Bayesian skyline plot analysis (Heled and Drummond, 2008). The general approach is, in several ways, similar to the one used in our ABC analysis but implemented in a MCMC coalescent sampler. The main difference is in the demographic model, which is a linear piecewise model (instead of constant piecewise model). That is, the model is defined by the scaled effective population size values θ i at n times points (i ∈ [0, n]). Effective population size changes linearly from θ i to θ i+1 between consecutive time points i and i + 1. The number of time points defining the model is explored using a Poisson distribution, which allows discriminating between stable and changing demographies as described in the main text. Regarding the microsatellite mutational model, BEASTS v1.8.0 allows for complex features, such as mutational bias between insertion and deletions of repeat units, allele size dependency of mutation rate and presence of multi-step mutations (Wu and Drummond, 2011). We have limited our analysis to the model equivalent to the GSM (called two-phase model with a single parameter in BEAST). For each data set, we run three independent chains of 10 9 steps, sampling parameters every 10 5 steps and discarding the first 10 8 steps as a burn-in, except for the red colobus for which three additional chains were run to achieve effective sample sizes larger than 100. Due to the computational burden of the analysis, we opted to reduce the sample size to 50 random individuals for the whale shark and leatherback turtle (original sample sizes were 478 and 215 individuals respectively) and to fix the GSM parameter to P GSM = 0.55 for the whale shark and P GSM = 0.56 for the leatherback turtle. These values were based on the estimates obtained with DIYABCskylineplot and MIGRAINE (see Table S3 and Figure S16). For the POD of a contracting population, due to problems of convergence and computing time, the P GSM was also fixed to the true value 0.22. For the rest of the analysis (colobus species and PODs from expanding and constant population) we used a uniform prior between zero and one for the P GSM parameter. All other options and priors were left to its default values.
The forth method is the composite-likelihood approach implemented in the R package VarEff (Nikolic and Chevalet, 2012). In this method, the data are reduced to the distribution of pairwise differences (difference in number of motifs) among gene copies for the sample. The probability that two gene copies differ by x number of motifs given a constant piecewise demographic model with parameters {(θ i , τ i ); i ∈ [0, n]} and a GSM with parameter P GSM can by calculated through numerical integration (Chevalet and Nikolic, 2010). Nikolic and Chevalet (2012) exploit this result to estimate the parameters of the demographic model with an MCMC algorithm. However, the mutational parameter P GSM is considered to be known and it is not estimated but fixed a priori. We run this analysis considering a model with three population size changes, uniform prior for the time τ = tµ between 0 and 4 (to explore the same time interval as with the ABC method), mutational parameter fixed to P GSM = 0.37 for the black-and-white colobus, P GSM = 0.215 for the red colobus, P GSM = 0.55 for the whale shark, P GSM = 0.56 for the leatherback turtle and P GSM = 0.22 (i.e. the true value) for the three PODs. These values were based on the estimates obtained with DIYABCskylineplot and MIGRAINE (see Table S3 and Figure S16). We run a chain with 10 5 samples, separated of 100 steps and taken after a burn-in period of 10 5 steps. All other parameters were set as recommended by authors (Nikolic and Chevalet, 2014).

S2.1 Alternative ABC skyline plot
The alternative parametrization to obtain a skyline plot on the ABC framework described in Supplementary section S1.2 yielded very similar results to the approach described in the main text (Supplementary Figure S17). Credibility intervals are narrower than in Figure 1, reflecting the stronger autocorrelation of effective population size through time when only one population size change is allowed. This should not be interpreted as a better performance of the inference as it reflects the differences in the prior setting. The ability to detect population size change is similar using both sets of priors (Supplementary Figure S18).

S2.2 Comparison with alternative methods
Graphical results of the four demographic inference methods are presented in Supplementary Figure S19 for PODS and Supplementary Figure S20 for real data. Demographic parameters estimated by MIGRAINE are presented in Supplementary Table S5. Evidence for variable demography from DIYABCskylineplot and BEAST is presented in Supplementary Table S7.
Regarding the analysis of PODs, were the true values are known, DIYABCskylineplot and MIGRAINE inferences described the three different demographies correctly co-estimating the mutational parameter P GSM . Credibility and confidence intervals contained the true values of parameters and were fairly congruent between both methods. Results from the extended Bayesian skyline plot from BEAST were less satisfactory. First, for the contracting scenario, co-estimation of demography and mutational model was computationally too expensive. Fixing the mutational parameter P GSM to its true value allowed finishing the computations in a reasonable time (i.e. few days) but is not a solution available for real data for which the parameter value will be unknown. For constant and expansion scenarios, demography and mutational model were co-estimated with mixed success. The expansion was well detected, though true values were consistently outside the credibility interval. In contrast, BEAST gave very strong evidence of variable demography for the constant model, with a skyline plot suggesting a demographic expansion and an estimate of P GSM far from the true value. Finally, VarEff performed well for characterizing the contraction and the constant size, but no the expansion. Note that original evaluation of the VarEff method already shows some limits in the detection of recent demographic increase, but it is able to characterize older expansions (Nikolic and Chevalet, 2014). Nevertheless, given that VarEff requires fixing the value of P GSM and we used the true value, these results are not of much relevance to the application of the method to real data.
The analyses of real data with the extended Bayesian skyline plot and VarEff highlight some of the problems described above with the additional difficulty of the lack of knowledge on the true demography and mutational process. Analyses with BEAST required reducing the sample size and fixing the mutational parameter to reduce the computation for two of the datasets (whale shark and leatherback turtle). In order to fix the values of P GSM for the analysis in BEAST and VarEff, we based our choice on the results from DIYABCskylineplot and MIGRAINE analysis of the same datasets. This is a practice that should be avoided in the analysis of real data and we did it as the only way to provide some comparison of results among methods. VarEff analysis did not reveal any striking demographic change in the history of those populations, suggesting a constant population size. Extended Bayesian skyline plots from BEAST were fairly flat for all four study cases and inferred systematically lower effective population sizes than the other methods. Considering the Bayes factor, evidence was either negative or weak for population size change: 1.36 for the whale shark, 1.75 for the leatherback turtle, 0.78 for the black-and-white colobus and 1.68 for the red colobus. Although the amount of data was much lower than that used by the other three methods for the whale shark and the leatherback turtle, the reduced datasets are similar in size to the PODs for which it was able to retrieve signal of demographic change.
The limited number of simulations does not allow us to conclude much about the performance of the different methods and the computation cost of BEAST precludes increasing the number of simulations. Nevertheless, this comparison is informative about the applicability of the methods to real data and their limitations. In our opinion, the use of VarEff on real data is difficult to justify, as it requires some knowledge on the mutational process that is rarely available. We want to stress that using the estimates of P GSM from the same dataset by an alternative method (such as DIYABCskylineplot or MIGRAINE) is not a good statistical practice and should not be done for the analysis of real data. The extended Bayesian skyline plot implemented in BEAST can potentially co-estimate the demography and mutational process if computation resources and time are not a constraint. However, the potential biases and pitfalls of the method have not been characterized. The failure to infer the correct demography for the simulated constant size scenario is anecdotal (many replicates would be necessary for a proper evaluation), but does not give confidence in the methods. Given that alternative methods with much less computation cost and thoroughly validated exist (i.e. DIYABCskylineplot and MIGRAINE), there is no reason to invest the computational resources required for the analysis in BEAST. Note that these remarks are valid in the context of the demographic inference from microsatellite data and should not be extrapolated to other analysis that are implemented in BEAST.     (Cornuet and Luikart, 1996). Estimates of genetic diversity are mean values from the 100 simulated data sets for each combination of demographic parameters (θ 0 , current scaled effective population size; θ 1 , ancestral scaled population size; τ , time since population size change). Standard deviations are shown in parentheses.  (Cornuet and Luikart, 1996).   Figure S17. Alternative ABC Skyline plots. Superimposed skyline plots (median in black, and 95%HPD interval in grey of the posterior probability distribution for θ(t)) from 100 replicates for example contraction (θ 0 = 0.4, θ 1 = 40, τ = 0.1), expansion (θ 0 = 40, θ 1 = 0.4, τ = 0.1) and constant size (θ = 40) scenarios with mutational model P GSM = 0.22. True demography is shown in orange. Note that present is at τ = 0 (left). These estimates are obtained through the alternative ABC skyline plot approach described in the Supplementary Methods S1.2.  Figure S18. Detection of demographic change with the ratio θ 0 /θ 1 in ABC. Point estimate (median, represented with black dots) and 95% credibility intervals (grey bars) for the ratio θ 0 /θ 1 , from 100 replicates for example contraction (θ 0 = 0.4, θ 1 = 40, τ = 0.1), expansion (θ 0 = 40, θ 1 = 0.4, τ = 0.1) and constant size (θ = 40) scenarios with mutational model P GSM = 0.22. Position of replicates in the abscissa is arbitrary. Credibility intervals not containing the value θ 0 = θ 1 (black horizontal line) will be considered as evidence for expansion or contraction. These estimates are obtained through the alternative ABC skyline plot approach described in the Supplementary Methods S1.2.  Figure S19. Comparison of methods on three simulated datasets. Black lines: skyline plots (median and 95%HPD interval of the posterior probability distribution) obtained with DIYABCskylineplot. Blue lines: skyline plots (median and 90%HPD interval of the posterior probability distribution) obtained with VarEff. Pink lines: skyline plots (median and 95%HPD interval of the posterior probability distribution) obtained with BEAST. Demographic trajectories based on parameters point estimates from MIGRAINE analysis are shown with a green line for reference. True demography is shown in orange. Note that present is at τ = 0 (left).