Demographic inference through approximate-Bayesian-computation skyline plots

ABC skyline plot

For a demographic skyline plot analysis within the ABC framework, our model consisted of a single population with constant size that instantaneously changes to a new size n times through time. The parameters (from present to past, as in the coalescent model) are the present scaled population size θ₀ = 4N₀μ (where N₀ is the effective population size in number of diploid individuals and µis the mutation rate per generation) which changes to θ₁ at time τ₁ = T₁μ (where T is the time measured in generations), remains at θ₁ and then it changes to θ₂ at τ₂, and so on, until the last change to θ_n at τ_n. Note that other models and parametrization could have been used for our purpose, as in the alternative model that we present in Section S1.2.

The objective of a standard ABC analysis would be to estimate the posterior distribution for each parameter of the model. In our case, the parameters {(θ_i, τ_i); i ∈ [0, n]} have been treated as nuisance parameters and we focused on inferring from them the trajectory of the scaled effective population size along time, θ(t), as in Drummond et al. (2005). In order to approximate θ(t) we select m times of interest, {t_j; j ∈ [1, m]}. Given a simulation k with parameters {(θ_k,i, τ_k,i); i ∈ [0, n_k]}, derived parameters {θ_k(t_j); j ∈ [1, m]} are obtained as follows: θ_k(t_j) = θ_k,i for i satisfying the condition τ_k,i ≤ t_j < τ_k,i+1 (see Fig. S1 for some examples). For each t_j, inference of the derived parameters θ(t_j) were obtained following standard ABC procedures as described elsewhere  (e.g., Beaumont, Zhang & Balding, 2002). Median and 95% highest posterior density (HPD) intervals of derived parameters θ(t_j) were used to draw ABC skyline plots.

Simulations with different numbers of population size changes can be used for inference because of the use of derived parameters θ(t_j), which are common to all models. We set the prior probability on the number of constant size periods to be Poisson distributed with λ = ln(2) as in Heled & Drummond (2008). This gives equal prior probability to stable populations (a single period of constant size) and changing populations (two or more periods). Thus, posterior probability on the number of periods may be used to discriminate between stable and changing demographies by estimating the Bayes factor of one period (constant population size) versus several demographic periods (variable population size). Posterior probabilities of contrasting models can be obtained by logistic regression as described elsewhere (Beaumont, 2008).

We implemented this approach in a suite of R scripts (R Core Team, 2017) that we named DIYABCskylineplot (Navascués, 2017). For each simulation the number of population size changes is sampled using the prior probabilities. Through a command line version of DIYABC (v2.0, Cornuet et al., 2014), parameter values {(θ_k,i, τ_k,i); i ∈ [0, n_k]} are sampled from the prior distribution, coalescent simulations are performed and summary statistics are calculated (mean across loci of the number of alleles, N_a; heterozygosity, H_e; variance of allele size, V_a, and Garza & Williamson (2001) statistic, M). In addition, the Bottleneck statistic (ΔH; Cornuet & Luikart, 1996), which compares the expected heterozygosity given the allele frequencies with the expected heterozygosity given the observed number of alleles, is calculated in R from the summary statistics provided by DIYABC. Derived parameter values, {θ_k(t_j); j ∈ [1, m]}, are calculated from the reference table (i.e., table of original parameters and summary statistics values for all simulations) produced by DIYABC and their posterior probability distributions are estimated in R using the abc package (Csilléry, François & Blum, 2012).

Simulations

The method described above was evaluated on simulated data (pseudo observed data-set, POD) of contracting and expanding populations. Declining populations had a present effective size of N₀ = 100 diploid individuals that changed exponentially until time T, which had a value of 10, 50, 100 or 500 generations in the past, reaching an ancestral population sizes of N_A, which had a value of 1,000, 10,000 or 100,000 individuals. Expanding populations had a present population size of N₀ with a value of 1,000, 10,000 or 100,000 diploid individuals, which changed exponentially until reaching the size of the ancestral population N_A = 100 at time T, which had a value of 10, 50, 100 or 500 generations in the past. For times older than T, the population size is constant at N_A for all scenarios. In addition, we simulated three constant population size scenarios with N taking a value of 1,000, 10,000 or 100,000. Equivalent scenarios were also evaluated in Girod et al. (2011) and Leblois et al. (2014). PODs were generated for 50 individuals genotyped at 30 microsatellite loci evolving under a generalised stepwise mutation model (GSM, Slatkin, 1995). Additional PODs varying in number of loci (7, 15 or 60 loci) and sample size (6, 12, 25 or 100 diploid individuals) were produced to evaluate the influence of the amount of data in the detection of demographic change. Mutation rate was set to μ = 10⁻³ and P_GSM to 0, 0.22 or 0.74 (P_GSM is the parameter of a geometric distribution determining the mutation size in number of repeats). One hundred replicates (i.e., PODs) were run for each scenario. Therefore, the mutation scaled parameter values are for θ = 4Nμ: 0.4, 4, 40 or 400 and for τ = Tμ: 0.01, 0.05, 0.1 or 0.5. PODs were obtained using the coalescent simulator fastsimcoal (Excoffier & Foll, 2011).

Every POD was analysed with the same set of prior probability distributions that largely includes all parameter values of simulations. Scaled effective size parameters, θ_i, were taken from a log-uniform distribution in the range (10⁻³ , 10⁴) and scaled times, τ_i, from a log-uniform distribution in the range (2.5 × 10⁻⁴,4). A uniform prior in the range (0, 1) was used for mutational parameter P_GSM. For each replicate of each scenario, we obtained the skyline plot (median and 95% HPD intervals of the θ(t_j) posterior distributions) and estimated the Bayes factor between constant size and variable demography by using logistic regression. Estimates of the mutational parameter P_GSM were also obtained for each POD. For each scenario, mean absolute error, bias and proportion of times the true value falls outside the credibility interval were estimated.

Data sets

In addition to PODs, four real data-sets from the literature were re-analysed with the ABC skyline plot described above. The first data-set comes from the whale shark (Rhincodon typus), the largest extant fish. Whale sharks inhabit all tropical and warm temperate seas and are considered an endangered species with a global population decline of more than 50% in the last three generations  (Pierce & Norman, 2016). We have re-analysed a data-set of 478 individuals genotyped at 14 microsatellite loci from Vignaud et al. (2014). The second example is the leatherback turtle (Dermochelys coriacea), the most widely distributed sea turtle found from tropical to sub-polar waters. The global population has been reduced in about 40% in the last three generations. As species, the leatherback turtle is classified as vulnerable, mainly because of the Northwest Atlantic population that shows an increase in number of nests. However, other populations are critically endangered ( Wallace, Tiwari & Girondot, 2013). The data-set re-analysed (215 individuals genotyped at 10 mocrosatellite loci; Molfetti et al., 2013) comes from the Northwest Atlantic population. Last, we re-analysed the data from the populations of two colobus monkeys at the Cantanhez National Park in Guinea Bissau (Minhós et al., 2016). The Western black-and-white colobus (Colobus polykomos, 22 individuals genotyped at 14 loci) and the Temminck’s red colobus (Procolobus badius ssp. temminckii, 23 individuals genotyped at 13 loci) are two sympatric species from the Western African rainforest considered to be vulnerable and endangered respectively ( Oates, Gippoliti & Groves, 2008; Galat-Luong et al., 2016). Data were analysed with the same prior distributions as PODs except for the colobus monkeys datasets, which consist of tetranucleotide markers. Previous evidence suggests that tetranucleotide microsatellite mutations are mainly of only one repeat unit  (e.g., Leopoldino & Pena, 2003; Sun et al., 2012). In order to incorporate this prior knowledge, half of the simulations had P_GSM = 0 (i.e., a strict stepwise mutation model, SMM) and the other half had the parameter sampled from a uniform distribution in the range (0, 1).

For comparison, demographic history of the four real data sets was also explored using the MIGRAINE software (Rousset & Leblois, 2016, http://kimura.univ-montp2.fr/ rousset/Migraine.htm) under the model of a single panmictic population with an exponential change in population size. To infer model parameters, MIGRAINE uses coalescence-based importance sampling algorithms under a maximum likelihood framework (Leblois et al., 2014) using OnePopVarSize model. In this model, MIGRAINE estimates present and ancestral scaled population sizes (θ₀ = 4N₀μ and θ_A = 4N_Aμ) and the scaled time of occurrence of the past change in population size (D = T∕4N, going backward from sampling time, when the population size change began). The past change in population size is deterministic and modelled using an exponential growth or decline that starts at time D. Before time D, scaled population size is stable and equal to θ_A. MIGRAINE allows departure from the strict SMM by using a GSM with parameter P_GSM for the geometric distribution of mutation sizes. Finally, detection of significant past change in population size is based on the ratio of population size (θ_ratio = θ₀∕θ_A). θ_ratio > 1 corresponds to a population expansion and θ_ratio < 1 to a bottleneck. If no significant demographic change is obtained, MIGRAINE is run again under a model of stable demography (a single value of θ) for parameter estimation. For the whale shark data set, MIGRAINE analysis was already done in Vignaud et al. (2014). For the leatherback turtle, MIGRAINE was run using 20,000 trees, 200 points at each iteration and a total of 16 iterations. For the colobus monkeys, we considered 2,000 trees, 400 points at each iteration and a total of 8 iterations.

Simulations

The general behavior of the method can be described from three example scenarios (contraction with θ₀ = 0.4, θ₁ = 40, τ = 0.1, expansion with θ₀ = 40, θ₁ = 0.4, τ = 0.1 and constant size with θ = 40; mutational model with P_GSM = 0.22). These examples correspond to intermediate parameter values. Results for all simulations are available in Supplementary Information 1.

The main output of the analysis is the graphical representation (i.e., the skyline plot) of the inferred demographic trajectory. It consists of a plot with three curves, representing the point estimates (median) and 95% HPD intervals of θ through time. Skyline plots obtained from PODs are congruent with the true underlying demography simulated (Fig. 1), except in the less favorable scenarios with very recent or very small changes in population size (Figs. S2–S8). Although the trajectory of the posterior median of θ and the true trajectory share the same trend (declining, increasing or constant), they sometimes differ in magnitude or time-scale. This disparity is more prominent for bottleneck scenarios.

For a quantitative criterion to assert demographic change we explored the value of posterior probabilities for constant and variable population size models, similar to the scheme proposed by Heled & Drummond (2008). These probabilities (summarised as Bayes factors in Fig. 2) proved to be useful for distinguishing bottleneck and expansion scenarios from demographic stability, although with lower performance for less favorable scenarios (Figs. S9–S15). Constant size scenarios show no evidence for size change. The power to detect demographic change reduces with smaller sample size and lower number of loci (Fig. 2) because summary statistics are estimated with lower precision.

Changes in population size were co-estimated with the mutational model parameter P_GSM. Mean absolute error, bias and proportion of replicates for which the true value was outside the 95% HPD interval are reported in Table 1 for the three example scenarios and in Table 1 for all simulations. Estimates from expanding and stable populations show a relatively low error and bias and a good coverage of the credibility interval (except in the strict SMM case). However, estimates from declining populations show higher error and bias.

Real data

The ABC analyses show evidence of population expansion for the whale shark (BF = 59.62) and the leatherback turtle (BF = 16.65); no evidence for population size changes in the black-and-white colobus (BF = 0.58) and some evidence for a bottleneck in the red colobus (BF = 2.63). Respective skyline plots reflect such trends (Fig. 3). Results from MIGRAINE support the same trends, with θ_ratio significantly larger than one for the whale shark and the leatherback turtle, significantly smaller than one for the red colobus and not significantly different from one for the black-and-white colobus (Table S3). Scaled population size estimates through time are also in agreement, except for the leatherback turtle, where the MIGRAINE result suggests a more ancestral expansion of much greater magnitude.

Regarding the mutational model, a large proportion of multi-step mutations seems to be present in all datasets, with P_GSM estimates: ${\hat{P}}_{GSM}=0.55$ (95% HPD = 0.46–0.62) for the whale shark; ${\hat{P}}_{GSM}=0.50$ (95% HPD = 0.38–0.60) for the leatherback turtle; ${\hat{P}}_{GSM}=0.43$ (95% HPD = 4.05 ×10⁻³–0.53) for the black-and-white colobus; and ${\hat{P}}_{GSM}=0.18$ (95% HPD = 0.02–0.75) for the red colobus (see also Fig. S16). Although very small values of P_GSM are included in the credibility interval from the colobus analyses, the GSM is favoured over the SMM when an ABC model choice analysis is performed (BF = 57.50 for the black-and-white colobus and BF = 10.01 for the red colobus). These results are congruent with estimates of P_GSM by MIGRAINE (Table S3).