 Methodology article
 Open access
 Published:
An improved approximateBayesian modelchoice method for estimating shared evolutionary history
BMC Evolutionary Biology volume 14, Article number: 150 (2014)
Abstract
Background
To understand biological diversification, it is important to account for largescale processes that affect the evolutionary history of groups of codistributed populations of organisms. Such events predict temporally clustered divergences times, a pattern that can be estimated using genetic data from codistributed species. I introduce a new approximateBayesian method for comparative phylogeographical modelchoice that estimates the temporal distribution of divergences across taxa from multilocus DNA sequence data. The model is an extension of that implemented in msBayes.
Results
By reparameterizing the model, introducing more flexible priors on demographic and divergencetime parameters, and implementing a nonparametric Dirichletprocess prior over divergence models, I improved the robustness, accuracy, and power of the method for estimating shared evolutionary history across taxa.
Conclusions
The results demonstrate the improved performance of the new method is due to (1) more appropriate priors on divergencetime and demographic parameters that avoid prohibitively small marginal likelihoods for models with more divergence events, and (2) the Dirichletprocess providing a flexible prior on divergence histories that does not strongly disfavor models with intermediate numbers of divergence events. The new method yields more robust estimates of posterior uncertainty, and thus greatly reduces the tendency to incorrectly estimate models of shared evolutionary history with strong support.
Background
Understanding the processes that generate biodiversity and regulate community assembly is a major goal of evolutionary biology. Largescale changes to the environment, including geological and climatic events, can affect the evolutionary history of entire communities of codistributed species and their associated microbiota. For example, by partitioning communities, such an event can isolate groups of populations and cause a temporal cluster of speciation events across codistributed taxa. Given the dynamic nature of our planet, such biogeographical processes likely play a significant role in determining diversification rates and patterns. At recent timescales, temporal clusters of diversification caused by biogeographical events can leave a signature in the genetic variation within and among the affected lineages. Thus, methods for accurately estimating models of shared evolutionary events across codistributed taxa from genetic data are important for better understanding how regional and global biogeographical processes affect biodiversity.
This inference problem is challenging due to the stochastic nature by which mutations occur in populations and how they are inherited over generations [1, 2]. Thus, a method for estimating historical patterns of divergences across taxa should explicitly model the stochastic mutational and ancestral processes that generate and filter the genetic variation we observe in presentday genetic data. An appealing approach would be a comparative, Bayesian modelchoice method for inferring the probability of competing divergence histories while integrating over uncertainty in mutational and ancestral processes via models of nucleotide substitution and lineage coalescence. The sample space of such a modelchoice procedure would include all models ranging from a single divergencetime parameter (i.e., simultaneous divergence of all codistributed taxa) to the fully generalized model in which each taxon diverged at a unique time.
The software package msBayes implements such an approach in an approximateBayesian modelchoice framework [3, 4]. The method models temporally clustered divergences across taxa caused by a biogeographical event (or a “divergence event”) as a single, instantaneous occurrence. In other words, a divergence event causes a set of taxa to share the same moment of divergence along a continuous time scale (i.e., simultaneous divergence). Given aligned sequence data for Y pairs of populations, msBayes estimates the number of divergence events shared among the pairs, the timing of the events, and the assignment of pairs to the events, while integrating out uncertainty in demographic parameters and the genealogical histories of the sequences. Thus, the method samples over all possible divergence models of differing dimensionality (i.e., all the possible partitions of Y pairs to 1,2,…,Y divergencetime parameters), and, in so doing, estimates the posterior probability of each model.
msBayes has been used to address biogeographical questions in a variety of empirical systems. Some examples include (1) whether the rise of the Isthmus of Panama caused codivergence among species of echinoids codistributed across the Pacific and Atlantic sides of the isthmus [3], (2) if an historical seaway across the Baja Peninsula caused codivergence across species of squamates and mammals codistributed both north and south of the putative seaway [5], (3) if species of gallwasps and their associated parasitoids share divergences across putative glacial refugia [6], and (4) whether repeated fragmentation of the oceanic Islands of the Philippines during Pleistocene sealevel fluctuations caused diversification of vertebrate taxa distributed across the islands [7]. Such applications of the method often result in strong posterior support for codivergence among all or subsets of the taxa investigated (e.g., [3, 5–12]).
For priors on divergencetime and demographic parameters, msBayes uses continuous uniform probability distributions. This causes divergence models with more divergencetime parameters to integrate over a much greater parameter space with low likelihood yet high prior density, which can result in small marginal likelihoods relative to models with fewer divergencetime parameters [13, 14]. Given that the marginal likelihood of a model weighted by its prior is what determines its posterior probability, this can cause support for models with fewer divergence events [7, 15]. This is not a critique of Bayesian model choice in general; comparing models by their marginal likelihoods provides a “natural” penalty for overparameterization and can be a great strength of the Bayesian approach. However, given the sensitivity of marginal likelihoods to the prior, care is needed when selecting prior distributions [14]. Selecting distributions that will often place high prior density in large regions of parameter space with low likelihood can lead to small marginal likelihoods of parameterrich models even if they are correct.
Furthermore, msBayes uses a discrete uniform prior over the number of divergence events 1,2,…,Y. Because there are many more possible assignments of population pairs to intermediate numbers of divergence events, this imposes a prior on divergence models that puts most of the prior mass on models with either very few or very many divergencetime parameters (see Figure five of [7]; for brevity I will refer to this prior as “Ushaped”). Given that models with many divergence events can have small marginal likelihoods due to the uniform priors on divergencetime parameters, the Ushaped prior will effectively create a strong prior preference for models with very few divergence events.
Recently, Oaks et al. [7, 15] found via simulation that msBayes will often strongly support models with a small number of divergence events shared among taxa, even when divergences were random over broad timescales. They suggested this behavior was due to the combination of uniform priors on parameters causing small marginal likelihoods of richer models and the Ushaped prior on divergence models. Hickerson et al. [16] suggested the problem was caused by sampling error, and proposed as a solution an approximateBayesian model averaging approach that samples over empirically informed uniform priors. However, Oaks et al. [15] evaluated the approach proposed by Hickerson et al. [16] using simulations and found that it did not mitigate the method’s propensity to incorrectly infer clustered divergences, and often preferred priors that excluded the true values of the model’s parameters. Here, I describe a new approach that successfully mitigates spurious inference of codivergence while avoiding negative side effects of empirically informed uniform priors.
In this study, I introduce a new method, implemented in the software dppmsbayes, that extends the model of msBayes. I use this method to test whether alternative parameterizations and priors improve the behavior of the approximateBayesian modelchoice approach to estimating shared divergence events. The new approach uses a Dirichletprocess prior (DPP) over all possible models of divergence, and gamma and beta probability distributions in place of uniform priors on many of the model’s parameters. Using simulations, I show that the new implementation has improved robustness, accuracy, and power compared to the original model. The results confirm that the improved performance of the new model is due to a combination of (1) more flexible priors on divergencetime and demographic parameters that avoid placing high prior density in improbable regions of parameter space, and (2) a diffuse Dirichletprocess prior that does not strongly disfavor divergence models with intermediate numbers of divergence events. After reanalyzing sequence data from 22 pairs of taxa from the Philippines [7] under the new model, I find a large amount of posterior uncertainty in the number of divergence events shared among the taxa; a result in contrast with the original msBayes model and congruent with intuition given the richness of the model and the relatively small amount of information in the data.
Methods
The model
In this section, I describe the model, which is a modification of the model implemented in msBayes[4, 7]. The code implementing the new model is freely available in the opensource software package dppmsBayes (https://github.com/joaks1/dppmsbayes). To perform the analyses described below, I used the freely avaliable, opensource software package PyMsBayes (https://github.com/joaks1/PyMsBayes), which provides a multiprocessing interface to msBayes and dppmsBayes. I performed the work described below following the principles of Open Notebook Science. Using versioncontrol software, I make progress in all aspects of the work freely and publicly available in realtime at https://github.com/joaks1/msbayesexperiments. All information necessary to reproduce my results is provided there. I follow much of the notation of Oaks et al. [7], but modify it to aid in the description of the new model. A summary of my notation can be found in Table 1.
I assume an investigator is interested in inferring the distribution of divergence times among Y pairs of populations. For each pair i, n_{ i } genome copies have been sampled, with n_{1,i} copies sampled from population 1, and n_{2,i} sampled from population 2. From these genomes, let k_{ i } be the number of DNA sequence loci collected for population pair i, and K be the total number of unique loci sampled across the Y pairs of populations. I use X_{i, j} to represent the multiple sequence alignment of locus j for population pair i. \mathbf{X}=({X}_{1,1},\dots ,{X}_{Y,{k}_{Y}}) is the full dataset, i.e., a vector of sequence alignments for all pairs and loci. Let G_{i, j} represent the gene tree upon which X_{i, j} evolved according to fixed HKY85 substitution model parameters ϕ_{i,j}. The investigator must specify the parameters of all \mathit{\varphi}=({\varphi}_{1,1},\dots ,{\varphi}_{\mathit{\text{Y}},{k}_{Y}}) substitution models by which the alignments evolved along the \mathbf{G}=({G}_{1,1},\dots ,{G}_{\mathit{\text{Y}},{k}_{Y}}) gene trees. Furthermore, the investigator must specify a vector of fixed constants \mathit{\rho}=({\rho}_{1,1},\dots ,{\rho}_{\mathit{\text{Y}},{k}_{Y}}) that scale the populationsize parameters for known differences in ploidy among loci and/or differences in generation times among population pairs. Lastly, the investigator must also specify a vector of fixed constants \mathit{\nu}=({\nu}_{1,1},\dots ,{\nu}_{\mathit{\text{Y}},{k}_{Y}}) that scale the populationsize parameters for known differences in mutation rates among loci and/or among taxa.
With X,ϕ,ρ, and ν in hand, the joint posterior distribution of the model is given by Bayes’ rule as
which can be expanded using the chain rule of probability into components that are assumed to be independent to get
where T=(T_{1},…,T_{ Y }) is a vector of population divergence times for each of the Y pairs of populations, Θ=(Θ_{1},…,Θ_{ Y }) is a vector of the demographic parameters for each of the Y population pairs, υ=(υ_{1},…υ_{ K }) is a vector of locusspecific mutationrate multipliers for each of the K loci, α is the shape parameter of a gammadistributed prior on υ, and p(Xϕ,ρ,ν), is the probability of the data (or the marginal likelihood of the model) given the fixed constants provided by the investigator.
To avoid calculating the likelihood terms of Equation 2, I distill each sequence alignment X into a vector of insufficient summary statistics S, thus replacing the full dataset \mathbf{X}=({X}_{1,1},\dots ,{X}_{Y,{k}_{Y}}) with vectors of summary statistics for each alignment {\mathbf{S}}^{\ast}=\left({\mathit{\text{S}}}_{1,1}^{\ast},\dots ,{S}_{Y,{k}_{Y}}^{\ast}\right) Optionally, for each population pair, the means of the summary statistics can be calculated across the k loci, and the vector can be further reduced to {\mathbf{S}}^{\ast}=\left({S}_{1}^{\ast},\dots ,{S}_{Y}^{\ast}\right). With S^{∗} in hand, we can estimate the approximate joint posterior distribution
where B_{ ε }(S^{∗}) is the multidimensional Euclidean space around the vector of summary statistics, the radius of which is the tolerance ε. The sources of approximation are the insufficiency of the statistics and the ε being greater than zero. I describe the full model in detail before delving into the numerical method of estimating the approximate model.
Likelihood and genetree prior terms of Equation 2
The likelihood and genetree prior terms of Equation 2 can be expanded out as a product over population pairs and loci
The first term, p(X_{i, j}G_{i, j},ϕ_{i, j}), is the probability of the sequence alignment of locus j for population pair i given the gene tree and HKY85 [17] substitution model parameters [18, i.e., the “Felsenstein likelihood”]. The model allows for an intralocus recombination rate r, which, for simplicity, is assumed to be zero in Equation 2. If r is nonzero, this term requires an additional product over the columns (sites) of each sequence alignment to allow sites to have different genealogies. The second term, p(G_{i, j}T_{ i },Θ_{ i },υ_{ j },ρ_{i, j},ν_{i, j}), is the probability of the gene tree under a multipopulation coalescent model (i.e., species tree) where the ancestral population of pair i diverges and gives rise to the two sampled descendant populations. Each Θ contains the following demographic parameters: The mutationratescaled effective sizes (θ = 4N μ) of the ancestral, θ_{ A }, and descendant populations, θ_{D1} and θ_{D2}; the proportion of the first, ζ_{D1}, and second population, ζ_{D2}, that persist during bottlenecks that begin immediately after divergence in forwardtime; the proportion of time between present and divergence when the bottlenecks end for both populations, τ_{ B }; and the symmetric migration rate between the descendant populations, m. Thus, the probability of the n_{ i }−1 coalescence times (node heights) of gene tree G_{i, j} is given by a multipopulation Kingmancoalescent model [19] where the ancestral population of size θ_{A,i}ρ_{i, j}ν_{i, j}υ_{ j } diverges at time T_{ i } into two descendant populations of constant size θ_{D1,i}ρ_{i, j}ν_{i, j}υ_{ j }ζ_{D1,i} and θ_{D2,i}ρ_{i, j}ν_{i, j}υ_{ j }ζ_{D2,i}, which, after time T_{ i }τ_{B,i}, grow exponentially to their present size θ_{D1,i}ρ_{i, j}ν_{i, j}υ_{ j } and θ_{D2,i}ρ_{i, j}ν_{i, j}υ_{ j }, respectively. Following divergence, the descendant populations of pair i exchange migrants at a symmetric rate of m_{ i }.
Additional prior terms of Equation 2
The term p(α) is the prior density function for the shape parameter of the gammadistributed prior on rate heterogeneity among loci. This prior is α∼U(1,20). The prior probability of the vector of locusspecific mutationrate multipliers given α then expands out as a product over the loci
where each υ is independently and identically distributed (iid) as υ∼G a m m a(α,1/α). If the recombination rate r is allowed to be nonzero, the prior term p(r) would be added to Equation 2, and the prior would be r∼G a m m a(a_{ r },b_{ r }), where a_{ r } and b_{ r } are specified by the investigator.
The prior term for the demographic parameters, p(Θ), expands out into its components and as a product over the Y pairs of populations
The priors for the demographic parameters are {\theta}_{A}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{A}},{b}_{{\theta}_{A}}),{\theta}_{D1}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}),{\theta}_{D2}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}),{\zeta}_{D1}\sim \mathit{\text{Beta}}({a}_{{\zeta}_{D}},{b}_{{\zeta}_{D}}),{\zeta}_{D2}\sim \mathit{\text{Beta}}({a}_{{\zeta}_{D}},{b}_{{\zeta}_{D}}),{\tau}_{B}\sim U(0,1), and m∼G a m m a(a_{ m },b_{ m }), where the hyperparameters of each prior distribution can be specified by the investigator. By default, θ_{ A }, θ_{D1}, and θ_{D2} share the same prior (i.e., {a}_{{\theta}_{A}}={a}_{{\theta}_{D}} and {b}_{{\theta}_{A}}={b}_{{\theta}_{D}}), but a separate gammadistributed prior can be assigned to θ_{ A }. Also, the ζ_{D1},ζ_{D2}, and m parameters are optional (i.e., the investigator can assume that there has been no migration between populations of each pair and/or the population size of each descendant population has been constant through time).
Priors on divergence models
The prior term for the vector of divergence times for each of the Y pairs of populations, T, can be expanded as
where τ is an ordered set of divergencetime parameters {τ_{1},…,τ_{τ}} whose length τ can range from 1 to Y, and t is a vector of indices (t_{1},…,t_{ Y }), where t_{ i }∈{1,…,τ}. These indices map each of the Y pairs of populations to a divergencetime parameter in τ. Thus, T is the result of applying the mapping function
to each population pair i, such that T = (T_{1} = f(τ,t,1),…,T_{ Y }=f(τ,t,Y)).
Biologically speaking, τ contains the times of divergence events, the length of which τ is the number of divergence events shared across the Y pairs of populations. For example, if τ contains a single divergencetime parameter τ_{1}, all Y pairs of populations are constrained to diverge at this time (i.e., t would contain the index 1 repeated Y times, and T would contain the value τ_{1} repeated Y times), whereas if it contains Y divergencetime parameters, the model is fully generalized to allow all of the pairs to diverge at unique times.
Unlike the model implemented in msBayes, here I place priors on t and τ, rather than τ and τ. As a result, t determines the number of divergencetime parameters (τ) in the model. Below, I first describe the prior used for τ and the timescale it imposes on the model before discussing the priors implemented for t.
Each τ within τ is iid as τ∼G a m m a(a_{ τ },b_{ τ }), where a_{ τ } and b_{ τ } are specified by the investigator. Thus, given the number of unique divergencetime classes in t, this determines the probability of prior term p(τt). The divergence times are in coalescent units relative to the size of a constant reference population, which I denote θ_{ C }, that is equal to the expectation of the prior on the size of the descendant populations
Given the size of the descendant populations are iid as {\theta}_{D1}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}) and {\theta}_{D2}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}), this becomes
More specifically, the τ parameters are in units of θ_{ C }/μ generations, which I denote as 4N_{ C } generations. Thus, each τ within τ is proportional to time and can be converted to the number of generations of the reference population, which I denote {\tau}_{{G}_{C}}, by assuming a mutation rate and multiplying by the effective size of the reference population
Thus, for each of the divergence times in τ to be on the same scale, the relative mutation rates among the pairs of populations are assumed to be known and fixed according to the userprovided values in ν.
As described by Oaks et al. [7], to get the divergence times in units proportional to the expected number of mutations, they must be scaled by the realized population size for locus j of populationpair i
where {\stackrel{\u0304}{\theta}}_{D,i} is the mean of θ_{D1} and θ_{D2} for pair i. This gives us the vector of scaled divergence times \mathcal{T}=({\mathcal{T}}_{1,1},\dots ,{\mathcal{T}}_{\mathit{\text{Y}},{k}_{Y}}).
As for the prior term p(t), the total sample space of t is all the possible partitionings of the Y pairs of populations into 1 to Y divergencetime classes, where each partitioning consists of nonoverlapping and nonempty subsets whose union is the Y pairs. Hereinafter, I refer to these partitionings as “ordered” divergence models or partitions. The total number of possible partitions is a sum of the Stirling numbers of the second kind over all possible numbers of categories τ
which is the Bell number [20]. The original msBayes model samples over the unordered realizations of t, such that the sample space is reduced to all the possible integer partitions of Y[4, 7, 21–23] (Additional file 1: Table S1). I denote the set of all possible integer partitions of the Y pairs of populations as a(Y) and the length of that set as a(Y), and I hereinafter refer to these integer partitions as “unordered” divergence models or partitions. The advantages, disadvantages, and justification of ignoring the order of t is discussed in detail below.
I implement two prior probability distributions over the space of all possible divergence models (t). The first simply gives all possible unordered partitions of Y elements equal probability
i.e., a discrete uniform prior over all the integer partitions of Y (unordered divergence models). I denote this prior as t∼D U{a(Y)}.
The second prior is based on the Dirichlet process, which is a stochastic process that groups random variables into an unknown number of discrete parameter classes [24, 25]. The Dirichlet process has been used as a nonparametric Bayesian approach to many inference problems in evolutionary biology [26–31]. Here, I use the Dirichlet process to place a prior over all possible ordered partitions of Y population pairs into divergencetime parameter classes (i.e., “divergence events”). As discussed above, the time of each divergencetime parameter is drawn from the base distribution τ∼G a m m a(a_{ τ },b_{ τ }). The partitioning of the population pairs to divergencetime classes is controlled by the concentration parameter χ, which determines how clustered the process will be. I take a hierarchical approach and use a prior probability distribution (i.e., hyperprior) for χ[32]. More specifically, I use a gammadistributed prior χ∼G a m m a(a_{ χ },b_{ χ }), where a_{ χ } and b_{ χ } are specified by the investigator. I use t∼D P(χ) to denote this Dirichletprocess prior.
This provides a great deal of flexibility for specifying the prior uncertainty regarding divergence models. The concentration parameter χ determines the prior probability that any two pairs of populations i and j will be assigned to the same divergencetime parameter
and also the prior probability of the number of divergencetime parameters
where c(·,·) are the unsigned Stirling numbers of the first kind. Equations 15 and 16 show that smaller values of χ will favor fewer divergencetime parameters, and thus more clustered models of divergence, whereas larger values will favor more divergencetime parameters, and thus less clustered models of divergence.
Differences between this model and the original msBayesmodel
The prior on divergence models
One of the key differences between this model and that of msBayes[4] is the prior distribution on divergence models. As discussed in Oaks et al. [7], in msBayes the prior used for t is a combination of a discrete uniform prior over the possible number of divergence events τ from 1 to Y with a multinomial distribution on the number of times each index of τ appears in t, with the constraint that all τ parameters are represented at least once (see Equation two of [7]). I denote this prior used by msBayes as t∼D U{1,…,Y}. Oaks et al. [7] discuss how placing a uniform prior over the number of divergence parameters (denoted τ here, and as Ψ in [4]) imposes an “Ushaped” prior over divergence models (t; see Figure five(B) of [7]). To avoid this, I place priors directly on the sample space of divergence models, thus eliminating the parameter Ψ from the model. I introduce two priors on divergence models: (1) a prior that is uniform over all unordered divergence models, and (2) a Dirichletprocess prior on all ordered divergence models. The latter provides an investigator with a great deal of flexibility in expressing their prior beliefs about models of divergence.
Estimating ordered divergence models
As mentioned above, msBayes samples over unordered divergence models (i.e., unordered partitions of the Y pairs of populations). That is, the identity of each population pair, and all the information associated with it, is discarded. In my implementation, inference can be done on either unordered or ordered models of divergence. This is discussed in more detail in the description of the ABC implementation below.
The priors on nuisance parameters
I have replaced the use of continuous uniform distributions for priors on many of the model’s parameters (τ,θ_{ A },θ_{D1},θ_{D2},ζ_{D1},ζ_{D2},r,m) with more flexible parametric distributions from the exponential family. I introduce gammadistributed priors for rate parameters that have a sample space of all positive real numbers (τ,θ_{ A },θ_{D1},θ_{D2},r,m), and betadistributed priors for parameters that are proportions bounded by zero and one (ζ_{D1} and ζ_{D2}). These priors provide an investigator with much greater flexibility in expressing prior uncertainty regarding the parameters of the model.
In addition, I have modified the prior on the sizes of the descendant populations of each pair. As described by Oaks et al. [7], msBayes uses the joint prior
such that the userspecified uniform prior on descendant population size is a prior on the mean size of the two descendant populations of each pair. Under my model, the sizes of the descendant populations of each pair are iid as {\theta}_{D1}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}) and {\theta}_{D2}\sim \mathit{\text{Gamma}}({a}_{{\theta}_{D}},{b}_{{\theta}_{D}}). This relaxes the assumption that the sizes of the two descendant populations are interdependent and negatively correlated.
Flexibility in parameterizing the model
In the new implementation, I provide the ability to control the richness of the model. For the θ parameters, by default, the model is fully generalized to allow each population pair to have three parameters: θ_{ A },θ_{D1}, and θ_{D2}. Furthermore, if an investigator prefers to reduce the number of parameters, any model of θ parameters nested within this general model can also be specified, including the most restricted model where the ancestral and descendant populations of each pair share a single θ parameter.
I also provide the option of eliminating the parameters associated with the postdivergence bottlenecks in the descendant populations of each pair (τ_{ B }, ζ_{D1}, and ζ_{D2}), which constrains the descendant populations to be of constant size from present back to the divergence event. Also, rather than eliminate the bottleneck parameters, I allow ζ_{D1} and ζ_{D2} to be constrained to be equal, which removes one free parameter from the model for each of the population pairs.
Overall, my implementation allows an investigator to specify a model that has as many as seven parameters per population pair (θ_{ A },θ_{D1},θ_{D2},τ_{ B },ζ_{D1},ζ_{D2}, and m) or as few as one parameter per pair (θ), in addition to the n_{ i }−1 coalescencetime parameters (i.e., the node heights of the gene tree).
Time scale
As described above, divergence times are in units of θ_{ C }/μ generations, where θ_{ C } is the expectation of the prior on descendantpopulation size. As described by Oaks et al. [7], in msBayes, θ_{ C } is half of the upper limit of the continuous uniform prior on the mean of the descendant population sizes. This is only equal to the expectation of the prior if the lower limit of the prior is zero.
ABC estimation of the posterior of the model
Sampling from the prior
To estimate the approximate posterior of Equation 3, I use an ABC rejection algorithm. The first step of this algorithm entails collecting a random sample of parameter values from the joint prior and their associated summary statistics. Each sample is generated by (1) drawing values of all the model’s parameters, which I denote Λ, from their respective prior distributions; (2) simulating gene trees \mathbf{G}=({G}_{1,1},\dots ,{G}_{Y,{k}_{Y}}) for each locus of each population pair by drawing coalescent times from a multipopulation Kingmancoalescent model given the demographic parameters; (3) simulating sequence alignments \mathbf{X}=({X}_{1,1},\dots ,{X}_{Y,{k}_{Y}}) along the gene trees under the HKY85 substitution parameters \mathit{\varphi}=({\varphi}_{1,1},\dots ,{\varphi}_{Y,{k}_{Y}}) that have the same number of sequences and sequence lengths as the observed dataset; and (4) calculating population genetic summary statistics \mathbf{S}=({S}_{1,1},\dots ,{S}_{Y,{k}_{Y}}) from the simulated sequence alignments. Optionally, an additional step can be performed to reduce the summary statistics to the means across loci for each population pair to get S=(S_{1},…,S_{ Y }). Either way, S contains the same summary statistics as those estimated from the observed data S^{∗}. After repeating this procedure n times, we have a random sample of parameter vectors Λ=(Λ_{1},…,Λ_{ n }) from the model prior and their associated vectors of summary statistics \mathbb{S}=({\mathbf{S}}_{1},\dots ,{S}_{n}).
For all of the analyses below, I use four summary statistics for each pair of populations: π[33], θ_{ W }[34], π_{ n e t }[35], and S D(π−θ_{ W }) [36]. Furthermore, in addition to model parameters, each sample Λ also contains four statistics that summarize T: the mean (\stackrel{\u0304}{T}), variance \left({s}_{T}^{2}\right), dispersion index ({D}_{T}={s}_{T}^{2}/\stackrel{\u0304}{T}), and the number of divergence time parameters (τ). Previously, these have been denoted as E(τ)V a r(τ),Ω, and Ψ, respectively [3, 4, 7]. I use \stackrel{\u0304}{T} and {s}_{T}^{2} in place E(τ) and V a r(τ) to make clear that these values do not represent the prior or posterior expectation/variance of divergence times. I use D_{ T } in place of Ω to clarify that this is a statistic rather than a parameter of the model. Lastly, I use τ in place of Ψ, because the number of divergencetime parameters is no longer a parameter in the new implementation.
Obtaining an approximate posterior from the prior samples
I use a rejection algorithm to retain an approximate posterior sample of Λ from the prior sample Λ=(Λ_{1},…,Λ_{ n }). First, the observed summary statistics S^{∗}, and the summary statistics of the prior samples \mathbb{S}=({\mathbf{S}}_{1},\dots ,{S}_{\mathbf{n}}), are standardized using the means and standard deviations of the statistics from the prior sample (i.e., the prior mean is subtracted from each statistic, and the difference is divided by the prior standard deviation). After all statistics are standardized, the Euclidean distance between S^{∗} and each S within is calculated. The samples that fall within a range of tolerance ε around S^{∗} are retained. The range of tolerance is determined by specifying the desired number of posterior samples to be retained. Posthoc adjustment of the posterior sample can also be performed with a number of regression techniques [37–39]. For analyses below, I use the general linear model (GLM) regression adjustment [39] as implemented in ABCtoolbox v1.1 [40], which Oaks et al. [7] showed performs very similarly to weighted locallinear regression and multinomial logistic regression adjustments [37] for msBayes posteriors.
Ordering of taxonspecific summary statistics
As alluded to in the model description, msBayes does not maintain the order of the taxonspecific summary statistics S within each S. Rather, the summary statistics are reordered by descending values of average pairwise differences between the descendant populations (π_{ b }) [4, 41]. This has the advantage of reducing the sample space of possible divergence models t, but there are at least two disadvantages. First, additional information in the data is lost. By discarding the identity of the Y pairs of populations, all pairspecific information about the amount of data (e.g., the number of gene copies collected from each of the populations [n_{1} and n_{2}], the number of loci, and the length of the loci), and the taxon and locusspecific parameters (ϕ,ν,ρ, and υ) is lost. Second, the results are more difficult to interpret, because divergence models and parameter estimates cannot be directly associated to the taxa under study.
The reordering of the summary statistic vectors also has an important implication for the ABC algorithm. When calculating the Euclidean distance between the observed data and each simulated dataset, the summary statistics being compared often represent sequence alignments of different taxon pairs and/or loci. More specifically, the summary statistics calculated from the observed sequence alignments are being compared to summary statistics calculated from datasets simulated with potentially different (1) numbers of sequences (n_{1} and n_{2}), (2) length of alignments, (3) numbers of loci (k), (4) HKY85 model parameters (ϕ), (5) mutationrate multipliers (ν), and (6) ploidy multipliers (ρ).
In the original descriptions of the msBayes method [3, 4], this reordering is justified by the fact that the expected value of π_{ b } is unrelated to sample size n_{1} and n_{2} and thus exchangeable among pairs. This is incorrect for two reasons. First, the entire vector of summary statistics S for each pair of populations is reordered across pairs, which implies that the justification for reordering π_{ b } applies to all the statistics within each S. However, the expectations for statistics that estimate gross diversity (e.g., π and θ_{ W }) are not independent of sample size for structured populations (e.g., the divergent pairs of populations modeled by msBayes), and other statistics are not independent of sample size in general (e.g., S D(π−θ_{ W })). Second, and more importantly, having the same expectation does not ensure random variables are exchangeable. Rather, for variables to be exchangeable their marginal distributions must be the same (i.e., they must be identically distributed). None of the summary statistics used by msBayes, including π_{ b }, have this property when there is any variation among taxa or loci in the (1) numbers of sequences (n_{1} and n_{2}), (2) length of alignments, (3) numbers of loci (k), (4) HKY85 model parameters (ϕ), (5) mutationrate multipliers (ν), or (6) ploidy multipliers (ρ). Whenever such variation is present (i.e., nearly all empirical applications), the taxonspecific summary statistics S are not exchangeable, and the reshuffling of the summary statistic vectors is not mathematically valid.
The magnitude of the affect of this violation of exchangeability is not known. Huang et al. [4] demonstrated that the reordering of the summary statistic vectors can greatly increase the method’s tendency to infer a single divergence event. By definition, if the summary statistic vectors were exchangeable, the reordering would not change the likelihood or posterior (barring sampling error). Thus, the results of Huang et al. [4] suggest the reordering of the statistics is potentially introducing sizeable error to the analysis.
For comparability with msBayes, I maintain the option for reordering taxonspecific summary statistics by π_{ b }. However, by default, the order is preserved, and ordered divergence models are estimated. In all of the simulationbased analyses described below, the summary statistic vectors are exchangeable, because the simulated datasets have the same (1) numbers of sequences, (2) length of sequences, (3) numbers of loci, (4) HKY85 model parameters, (5) mutationrate multipliers, and (6) ploidy multipliers.
Assessing modelchoice behavior and robustness
Following the simulationbased approach of Oaks et al. [7], I characterize the behavior of several models under the ideal conditions where the data are generated from parameters drawn from the same prior distributions used for analysis (i.e., the prior is correct). I selected the following four model priors for these analyses (Table 2).

1.
The M _{ m s B a y e s } model represents the original msBayes implementation with the Ushaped prior on unordered divergence models and uniform priors on divergencetime and demographic parameters; t∼D U{1,…,Y},τ∼U(0,10),θ _{ A }∼U(0,0.05), and {\stackrel{\u0304}{\theta}}_{D}\sim U(0,0.05).

2.
The M _{ Ushaped } model with the Ushaped prior of msBayes on unordered divergence models, but with exponential priors on divergencetime and demographic parameters; t∼D U{1,…,Y},τ∼E x p(m e a n=2.887),θ _{ A }∼E x p(m e a n=0.025),θ _{D1}∼E x p(m e a n=0.025), and θ _{D2}∼E x p(m e a n=0.025).

3.
The M _{ U n i f o r m } model with a uniform prior over unordered divergence models and exponential priors on divergencetime and demographic parameters; t∼D U{a(Y)},τ∼E x p(m e a n=2.887),θ _{ A }∼E x p(m e a n=0.025),θ _{D1}∼E x p(m e a n=0.025), and θ _{D2}∼E x p(m e a n=0.025).

4.
The M _{ D P P } model with a Dirichletprocess prior on ordered divergence models and exponential priors on divergencetime and demographic parameters; t ∼ D P(χ∼G a m m a(2,2)),τ ∼ E x p(m e a n=2.887), θ _{ A } ∼ E x p(m e a n=0.025),θ _{D1} ∼ E x p(m e a n=0.025), and θ _{D2}∼E x p(m e a n=0.025).
I selected the exponential prior on divergence time used in models M_{ D P P },M_{ U n i f o r m }, and M_{ U s h a p e d } to have the same variance as the uniform prior in model M_{ m s B a y e s }. I selected the exponential prior on population size used in models M_{ D P P }, M_{ U n i f o r m }, and M_{ U n i f o r m } to have the same mean as the uniform prior in model M_{ m s B a y e s }, so that all four models have the same θ_{ C } and thus the same units of time. All of the models were the same in other respects, with three free θ parameters for each population pair, two uniformly distributed (b e t a(1,1)) ζ_{ D } parameters per pair, no migration, no recombination, and resorting of taxonspecific summary statistics by π_{ b } (i.e., sampling unordered divergence models). For all simulations, I used a data structure of eight population pairs, with a single 1000 basepair locus sampled from 10 individuals from each population.
For each of the four models, I simulated 1×10^{6} samples from the prior and 50,000 datasets, also drawn from the prior. I then analyzed each of the simulated datasets, retaining a posterior of 1000 samples from the respective prior. A GLMregression adjusted posterior was also estimated from each of the posterior samples [39]. To assess the robustness of each of the four models, I also analyzed the datasets simulated under the other three models. Overall, for each model, I produced 200,000 posterior estimates, 50,000 from the datasets simulated under that model, and 150,000 from the datasets simulated under the other three models.
For each set of 50,000 simulated datasets, I used the posterior estimates to assess the modelchoice behavior of each model. I did this by assigning the 50,000 estimates of the posterior probability of onedivergence event to 20 bins of width 0.05, and plotted the estimated posterior probability of each bin against the proportion of replicates in that bin with a true value consistent with one divergence event [7, 42]. Ideally, the estimated posterior probability of the onedivergence model should estimate the probability that the onedivergence model is correct. For large numbers of simulation replicates, the proportion of the replicates in each bin for which the onedivergence model is true will approximate the probability that the onedivergence model is the correct model. Thus, if the method has the desirable behavior such that the estimated posterior probability of the onedivergence model is an unbiased estimate of the probability that the onedivergence model is correct, the points should fall near the identity line. For example, let us say the method estimates a posterior probability of 0.90 for 1000 datasets simulated from the prior. If the method is accurately estimating the probability that the onedivergence model is correct given the data, then the onedivergence model should be the true model in approximately 900 of the 1000 replicates. Any trend away from the identity line indicates the method is biased in the sense that it is not accurately estimating the probability that the onedivergence model is the correct model.
I constructed these plots using two criteria for the onedivergence model: (1) the number of divergencetime parameters (τ=1) and (2) the dispersion index of divergence times (D_{ T }<0.01). For the latter, D_{ T }<0.01 has been commonly used as an arbitrary criterion for a single “simultaneous” divergence event (e.g., [3, 5, 6]). I focused on the onedivergence model to assess modelchoice behavior, because it is often of biogeographic interest and is easily comparable among the three different priors used on divergence models.
In addition to the four models above, I also assessed the behavior of a model that samples over ordered divergence models (i.e., the order of the taxonspecific summary statistic vectors were maintained for the observed and simulated datasets); all other settings were identical to the M_{ D P P } model. I denote this model as M°_{ D P P }. I simulated 1×10^{6} prior samples and 50,000 datasets, and analyzed them as above. I was not able to analyze the simulated datasets of the other models under the ordered model, because the identity of the population pairs is not contained in the simulations of the other models.
Assessing power
I evaluated the power of the same four models (Table 2) to detect random variation in divergence times using methods similar to Oaks et al. [7]. For all power simulations, I used a data structure identical to that of the empirical dataset of Philippine vertebrates analyzed by Oaks et al. [7], which consists of 22 pairs of populations. Due to the larger number of pairs, I used a different hyperprior on the concentration parameter for the M_{ D P P } model; I used a prior of t∼D P(χ∼G a m m a(1.5,18.1)) over divergence models for the model M_{ D P P }. All other aspects of the four models in Table 2 were identical to those used in the validation analyses described above. For each of the four models, I generated 2 × 10^{6} samples from the prior.
Next, I simulated datasets from three series of models in which the divergence times of the 22 pairs were random (i.e., no clustering; τ=22). The models comprising each series differ in the variance of the distribution from which the divergence times are randomly drawn. When the variance of random divergence times is small, all of the models in Table 2 are expected to struggle to detect this variation and will often incorrectly estimate highly clustered models of divergence (i.e., few divergence events). The goal is to assess how much temporal variation in random divergence times is necessary before the behavior of the models of Table 2 begins to improve. This will determine the timescales over which the models can reliably detect random variation in divergence times and avoid spurious inference of clustered divergence models.
Specifically, I simulated datasets from the following three series of six models (Table 3).

1.
The {\mathcal{\mathcal{M}}}_{\mathit{\text{msBayes}}} models are identically distributed as M _{ m s B a y e s } except the divergence times for each of the 22 pairs of populations are randomly drawn from a series of uniform distributions, U(0,τ _{ m a x }), where τ _{ m a x } was set to: 0.2, 0.4, 0.6, 0.8, 1.0, and 2.0, in 4N _{ C } generations.

2.
The {\mathcal{\mathcal{M}}}_{\mathit{\text{Uniform}}} models are identically distributed as M _{ U n i f o r m } and M _{ D P P } except the 22 divergence times are randomly drawn from the same series of uniform priors as above.

3.
The {\mathcal{\mathcal{M}}}_{\mathit{\text{Exp}}} models are also identically distributed as M _{ U n i f o r m } and M _{ D P P } except the 22 divergence times are randomly drawn from a series of of exponential distributions: E x p(m e a n=0.058), E x p(m e a n=0.115), E x p(m e a n=0.173), E x p(m e a n=0.231), E x p(m e a n=0.289), and E x p(m e a n=0.577). These exponential distributions have the same variance as their uniform counterparts in the first two series of models.
For each of the six models in each of the three series of models, I simulated 1000 datasets (18,000 datasets in total). I then analyzed each simulated dataset under all four prior models (Table 2), producing 72,000 posterior estimates, each with 1000 samples. I also estimated a GLMregression adjusted posterior from each of the posterior samples [39].
An empirical application
I also assessed the behavior of the newly implemented models when applied to the empirical dataset of Oaks et al. [7], which is comprised of sequence data from 22 pairs of taxa from the Philippine Islands ([43]; Dryad DOI: 10.5061/dryad.5s07m). I analyzed these data under five different models, which are detailed in Table 4. All of these models except one ({\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{simple}}}) have six free demographic parameters per pair of taxa (θ_{ A },θ_{D1},θ_{D2},τ_{ B },ζ_{D1},and ζ_{D2}), in addition to the n_{ i }−1 coalescent times. Three of these models use a Dirichletprocess prior on divergence models: {\mathbf{M}}_{\mathit{\text{DPP}}},{\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{inform}}},\text{and}\phantom{\rule{0.3em}{0ex}}{\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{simple}}}. The M_{ D P P } model represents the priors that Oaks et al. [7] would have selected to reflect their prior uncertainty about the parameters of the model if provided the more flexible distributions that are now implemented. To assess prior sensitivity, the {\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{inform}}} model uses a more informative exponentially distributed prior on divergence times, but otherwise is identical to M_{ D P P }. To assess sensitivity to parameterization, I also applied the simplest possible model under the new implementation \left({\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{simple}}}\right) with only a single demographic parameter (θ) per taxon pair, in addition to the n_{ i }−1 coalescent times. I also applied the original msBayes model (M_{ m s B a y e s }) with priors selected to make the results directly comparable to those of the M_{ D P P } model; the uniform prior on divergence times was selected to have the same variance as the exponential prior of the M_{ D P P } model, and the prior on population size was selected to have the same mean so that the models are on the same timescale. I also applied a model with a uniform distribution over divergence models (M_{ Uniform }). For each of these models, I simulated 2 × 10^{7} samples from the prior, and retained an approximate posterior of the 10,000 samples with the smallest Euclidean distance from the summary statistics calculated from the empirical sequence alignments.
To compare models that sample over ordered versus unordered models of divergence, I also analyzed the data from the subset of ninetaxon pairs that are sampled from the Islands of Negros and Panay in the Philippines. The model I used for these analyses had a Dirichletprocess prior over divergence models and two demographic parameters (θ_{ A } and (θ_{ D }) for each pair of taxa, in addition to the n_{ i }−1 coalescent times (see Table 4 for details). One of the models, which I denote \mathbb{M}{\text{\xb0}}_{\mathit{\text{DPP}}}, maintained the identity of the taxon pairs and sampled over ordered models of divergence, while the other ({\mathbb{M}}_{\mathit{\text{DPP}}}) resorted the summary statistics of the pairs by π_{ b }, losing the identity of the taxa and thus sampled over unordered models of divergence. For both analyses, I simulated 5 ×10^{7} samples from the prior and retained an approximate posterior of 10,000 samples.
Results
Validation analyses: Estimation accuracy
In terms of estimating the variance of divergence times (D_{ T }), the models with exponentially distributed priors (M_{ U s h a p e d }, M_{ U n i f o r m }, and M_{ D P P }) perform similarly when applied to datasets generated under all four of the models in Table 2 (Additional file 1: Figure S1). The M_{ m s B a y e s } model performs similarly to these models when applied to its own datasets, but is sensitive to model violations and is more biased when applied to data generated under the other three models (Additional file 1: Figure S1). Results are similar for the GLMadjusted estimates of D_{ T }, albeit the regression adjustment tends to improve estimates of this continuous statistic for all the models (Additional file 1: Figure S2).
The same general pattern is seen for estimates of \stackrel{\u0304}{T}, with (1) all four models performing similarly when applied to the data generated under the M_{ m s B a y e s } model, (2) the models with exponentially distributed priors performing similarly when applied to data generated under the other three models, and (3) the M_{ m s B a y e s } model is sensitive to model violations and is more biased whenever applied to data generated under other models (Additional file 1: Figure S3). Also, the regression adjustment tends to slightly improve estimates of this continuous statistic for all of the models (Additional file 1: Figure S4).
In terms of estimating the number of divergence events (τ), the M_{ D P P } model has the lowest root mean square error (RMSE) when applied to data generated under most of the models of Table 2 (Additional file 1: Figure S5). The M_{ m s B a y e s } model performs slightly better when applied to its own data, but is the worst performer when applied to data generated under other models (Additional file 1: Figure S5). There is a trend of M_{ D P P }>M_{ U n i f o r m }>M_{ U s h a p e d }>M_{ m s B a y e s } in terms of estimation accuracy as measured by RMSE when the models are applied to data generated under most of the models (Additional file 1: Figure S5). Unlike for the continuous statistics, regression adjustment of this discrete statistic tends to increase estimation bias; all of the models tend to underestimate τ after the GLMadjustment (Additional file 1: Figure S6).
Validation analyses: Modelchoice accuracy
The msBayes model, and my modification of it, is a modelchoice method with the primary purpose of estimating the probabilities of models of divergence across taxa. Thus, it is critical to assess the method’s ability to accurately estimate the posterior probabilities of divergence models. Consistent with the findings of Oaks et al. [7], my results demonstrate that the unadjusted estimates of divergencemodel posterior probabilities are generally more accurate than regressionadjusted estimates (compare the plots along the upperleft to lowerright diagonal for Figure 1 versus Additional file 1: Figure S7 and Figure 2 versus Additional file 1: Figure S8). Regression adjustment results in biased estimates of the posterior probability of the onedivergence model when all model assumptions are satisfied, which is well illustrated in Additional file 1: Figure S8. As a result, I will focus my discussion of the results on the unadjusted estimates.
I find that all four models accurately estimate the posterior probability of the onedivergence model when applied to their own datasets (i.e., when the prior is correct; see diagonal of Figures 1 and 2). The M_{ U n i f o r m } and M_{ D P P } models show robustness to prior violations and perform well when applied to data generated under other models (Figures 1 and 2). However, both are less accurate and tend to underestimate the probability of the onedivergence model when applied to the data generated under M_{ U s h a p e d } (Figures 1 and 2). In contrast, the M_{ m s B a y e s } model is biased toward overestimating the posterior probability of the onedivergence model when applied to data generated under the other three models (Figures 1 and 2). This bias is particularly strong whenever divergence models are not distributed under its Ushaped prior (Figure 1C–D). The other model with the Ushaped prior on divergence models, but exponential priors on parameters (M_{ U s h a p e d }), performs similarly to the M_{ m s B a y e s } model in that it performs well when applied to its own data, but overestimates the probability of the onedivergence model when applied to data generated by the other models (Figures 1 and 2). However, the bias is stronger in the M_{ m s B a y e s } model than M_{ U s h a p e d }.
Overall, the results suggest that the M_{ D P P } and M_{ U n i f o r m } models are relatively robust in terms of modelchoice accuracy, and when model violations do cause them to be biased, they tend to underestimate the probability of the model with a single, shared divergence event. In contrast, the M_{ m s B a y e s } model is very sensitive to model violations, and strongly overestimates the probability of the onedivergence model whenever the model is misspecified. Furthermore, the results suggest that using exponentially distributed priors on nuisance parameters rather than uniform priors helps the M_{ U s h a p e d } model perform better than M_{ m s B a y e s }, but it is still hindered by the Ushaped prior on divergence models and tends to overestimate the probability of the onedivergence model whenever there are violations of the model.
Validation analyses: Ordered divergence models
The results show that the method performs similarly when sampling over ordered models of divergence (Additional file 1: Figures S9 and S10). This suggests that the method is not adversely affected by the increase in the number of possible discrete models (from 22 unordered to 4140 ordered models) when there are eight pairs of populations. This is encouraging, because, as discussed above, estimating unordered models of divergence by shuffling the summary statistic vectors calculated from the sequence alignments is not valid for most empirical datasets. Given these results, estimation of unordered divergence models should be avoided for empirical applications of the method.
Power analyses: Estimation accuracy
All of the models I evaluated (Table 2) struggle to estimate the variance of divergence times D_{ T } regardless of which of the three series of models (Table 3) the data were generated under (Additional file 1: Figures S11–13). The models with the Ushaped prior on divergence models (M_{ m s B a y e s } and M_{ U s h a p e d }) tend to underestimate the variance in divergence times (Plots A–L of Additional file 1: Figures S11–13). whereas the models with Uniform or Dirichletprocess priors over divergence models tend to overestimate variance in divergence times (Plots M–X of Additional file 1: Figures S11–13).
When the divergence times of the 22 population pairs are randomly drawn from a series of exponential priors ({\mathcal{\mathcal{M}}}_{\mathit{\text{Exp}}}), the M_{ D P P } model is the best estimator of D_{ T }, followed by M_{ U n i f o r m } (Additional file 1: Figure S11). The M_{ m s B a y e s } model is strongly biased toward underestimating D_{ T }, estimating values of zero for most of the replicates across all the data models of {\mathcal{\mathcal{M}}}_{\mathit{\text{Exp}}} (Additional file 1: Figure S11). The results of the M_{ U s h a p e d } model are intermediate between those of M_{ m s B a y e s } and the new models M_{ D P P } and M_{ U n i f o r m } (Additional file 1: Figure S11).
Similarly, when the true divergence times are randomly drawn from a series of uniform priors ({\mathcal{\mathcal{M}}}_{\mathit{\text{Uniform}}}), the M_{ D P P } and M_{ U n i f o r m } models tend to overestimate the variance in divergence times, whereas the M_{ m s B a y e s } model underestimates D_{ T }, estimating values of zero for most replicates across all the data models of {\mathcal{\mathcal{M}}}_{\mathit{\text{Uniform}}} (Additional file 1: Figure S12). Again, the performance of the M_{ U s h a p e d } model is intermediate between the M_{ m s B a y e s } and M_{ D P P }/ M_{ U n i f o r m } models (Additional file 1: Figure S12). The results are very similar when the four models are applied to the data simulated under the {\mathcal{\mathcal{M}}}_{\mathit{\text{msBayes}}} series of models (Additional file 1: Figure S13).
Power analyses: Model choice
The modifications of the msBayes model decrease the method’s bias toward clustered divergences when applied to data generated under random divergence times (Figure 3 and Additional file 1: Figures S14–16). The M_{ m s B a y e s } model performs the worst of the four models across all three series of datagenerating models, inferring a single divergence event across most of the 18,000 simulations (Figure 3A–D and plots A–F of Additional file 1: Figures S14–16). Importantly, the M_{ m s B a y e s } model tends to strongly support these estimates of one divergence across most of the simulations (Figure 4A–D and plots A–F of Additional file 1: Figures S17–19). The M_{ D P P } model also prefers the onedivergence model when divergences are random over narrow windows of time, but performs much better when divergences are random over a timescale of 1–2 coalescent units (Figure 3M–P and plots S–X of Additional file 1: Figures S14–16). However, even when M_{ D P P } infers the onedivergence model over narrow timescales, the posterior probability support is always low (Figure 4M–P and plots S–X of Additional file 1: Figures S17–19). The M_{ U n i f o r m } model never infers the onedivergence model in any of the simulation replicates but still tends to infer relatively few (4–6) divergence events when divergences are random over longer periods (Figure 3I–L and plots M–R of Additional file 1: Figures S14–16). Using exponential priors on divergencetime and demographic parameters does increase the power of the M_{ U s h a p e d } model compared to M_{ m s B a y e s } across all three series of data models, but the Ushaped prior still prevents the model from performing as well as the M_{ D P P } and M_{ U n i f o r m } models (Figure 3 and Additional file 1: Figures S14–16).
The improved power of the new models is even more pronounced when looking at estimates of the variance of divergence times (D_{ T }) across the simulations (Figure 5 and Additional file 1: Figures S20–22). The performance among the models is so different, that the histograms of D_{ T } estimates cannot be plotted along a shared xaxis. The M_{ D P P } and M_{ U n i f o r m } models perform similarly across all three series of data models, inferring values of D_{ T } consistent with one divergence event (D_{ T }<0.01) in almost none of the replicates across all the simulations. In contrast, the M_{ m s B a y e s } model infers values consistent with a single divergence event in most of the replicates across all the simulations. Using exponential priors on divergencetime and demographic parameters greatly increases the power of the M_{ U s h a p e d } model to detect variation in divergence times relative to M_{ m s B a y e s }, but it still has less power than the models with Dirichletprocess or uniform priors across divergence models (Figure 5 and Additional file 1: Figure S20–22). Although the D_{ T } threshold of 0.01 is arbitrary, Oaks et al. [7] did show via simulation that the true value of D_{ T } will almost always be greater than 0.01 when divergences are random over periods of 0.1 coalescent units or more (see Figure Sfour of [7]).
As mentioned above, the increased power of the new models is also evident when looking at the estimated posterior probability of the onedivergence model across the power analyses (Figure 4 and Additional file 1: Figures S17–19). The M_{ D P P } and M_{ U n i f o r m } models estimate low posterior probability of τ=1 across all of the simulations. This is in contrast to the M_{ m s B a y e s } model, which infers high posterior probabilities of a single divergence for most replicates across all simulations (Figure 4 and Additional file 1: Figures S17–19). The exponential priors on divergencetime and demographic parameters (model M_{ U s h a p e d }) result in lower estimates of the probability of one divergence when compared to M_{ m s B a y e s }, but higher estimates when compared to M_{ U n i f o r m } and M_{ D P P } (Figure 4 and Additional file 1: Figures S17–19). The M_{ D P P } and M_{ U n i f o r m } models do frequently support the onedivergence model according to a Bayes factor criterion of greater than 10, but still less frequently than the M_{ m s B a y e s } model. This result is not surprising given the extremely small prior probability of the onedivergence model under the M_{ D P P } and M_{ U n i f o r m } models (i.e., very few posterior samples of the onedivergence model will result in a large Bayes factor under these models). However, the small posterior probability of the onedivergence model estimated under M_{ D P P } and M_{ U n i f o r m } should prevent an investigator from overinterpreting the Bayes factor as strong support for clustered divergences.
Lastly, when looking at the estimated posterior probability of D_{ T } being consistent with one shared divergence (p(D_{ T }<0.01B_{ ε }(S^{∗}))), I find the same pattern of model behavior, with M_{ D P P } and M_{ U n i f o r m } inferring low probabilities across all simulations, M_{ m s B a y e s } inferring high probabilities, and M_{ U s h a p e d } inferring intermediate values (Figure 6 and Additional file 1: Figures S23–25).
Empirical results
As expected based on the results of Oaks et al. [7], when the Philippines data are analyzed under the M_{ m s B a y e s } model, there is strong support for very few divergence events shared among all 22 pairs of taxa, with a maximum a posteriori (MAP) estimate of oneshared divergence (Figure 7A). When these data are analyzed using models allowed by the new implementation, there is much less support for highly clustered models and much greater uncertainty regarding the number of divergence events shared among the taxa, especially under the DPP models (Figure 7B–E). Figure 7 also shows the prior distribution across the number of divergence events (τ) for each model, as well as the average prior probability of an unordered and ordered model of divergence (t) across τ. Estimates under the new models tend to be similar to the prior, which is expected under such a parameterrich model when there is limited information from the data (four summary statistics from a single locus for each pair of taxa).
The disparity between the results of the M_{ m s B a y e s } model and the new models is even more pronounced when looking at the 10 divergence models (t) estimated to have the highest probability under each of the models (Additional file 1: Figures S26–30). Again, the new models estimate more divergences, a large amount of posterior uncertainty, and an order of magnitude smaller probability for their respective MAPdivergence model when compared to the M_{ m s B a y e s } model (Additional file 1: Figures S26–30).
Figure 8 shows the estimated posterior probability distribution over the number of divergence events when the data from the ninetaxon pairs from the Islands of Negros and Panay are analyzed under DPP models that sample over unordered ({\mathbb{M}}_{\mathit{\text{DPP}}}) and ordered (\mathbb{M}{\text{\xb0}}_{\mathit{\text{DPP}}}) models of divergence. The results are similar under both models and, again, yield a large amount of uncertainty about the number of divergence events that is similar to the prior uncertainty.
The small difference between the results of the {\mathbb{M}}_{\mathit{\text{DPP}}} and \mathbb{M}{\text{\xb0}}_{\mathit{\text{DPP}}} models is consistent across multiple analyses, and thus could be due to error introduced to the {\mathbb{M}}_{\mathit{\text{DPP}}} model by the invalid shuffling of the summary statistic vectors. Both models estimate a similar set of 10 unordered divergence models with the highest posterior probability (Additional file 1: Figures S31 and S32).
The main advantages of the \mathbb{M}{\text{\xb0}}_{\mathit{\text{DPP}}} model over the {\mathbb{M}}_{\mathit{\text{DPP}}} are that (1) the incorrect shuffling of the summary statistic vectors is avoided, (2) the identity of the taxa is maintained, and thus a fully marginalized estimate of divergence times across the taxa can be obtained (Additional file 1: Figure S33), and (3) the probability of codivergence among any set of taxa can be estimated from the posterior sample.
Discussion
My results demonstrate that using alternative priors on parameters and divergence models improved the behavior of the msBayes model. In the new implementation, modelchoice estimation is more accurate and shows greater robustness to model violations (Figures 1 and 2). The original model is very sensitive to violations and, when present, strongly overestimates the probability of onedivergence event shared across all taxa (Figures 1 and 2). When more appropriate priors are used for divergencetime and demographic parameters, and either a Dirichletprocess or uniform prior applied across divergence models, the model is less sensitive to violations, and, when violations do cause bias, the method tends to underestimate the probability of models with temporally clustered divergences (Figures 1 and 2). Given that clustered models are often of particular interest to biogeographers, this behavior of the new method can be considered conservative.
The modifications also improve the method’s power to detect random variation in divergence times, reducing the tendency to estimate clustered divergences (Figures 3, 4, 5 and 6). My results are similar to those of Oaks et al. [7] in that I find msBayes will often infer strong support for clustered divergences when divergences are random over quite broad timescales (Figures 3, 4, 5 and 6). My results expand on this by showing that this behavior is consistent across a range of conditions underlying the data. The new method, dppmsBayes, has greater power to detect random temporal variation in divergences, is less prone to spurious inference of clustered divergence models, and much less likely to incorrectly infer such models with strong support (Figures 3, 4, 5 and 6).
By evaluating a model intermediate between the old and new implementation (M_{ U s h a p e d }), I was able to determine the relative affects of my modifications to the model. Across all of the analyses, the results show that using better priors on divergencetime and demographic parameters alone does improve the performance of the method. The magnitude of the bias toward inferring support for the onedivergence model when there are model violations is reduced when the exponential priors are used in place of the uniform priors (Figures 1 and 2). Furthermore, using exponential priors improves the method’s power to detect temporally random divergences (Figures 3, 4, 5 and 6). Throughout the analyses, the intermediate model (M_{ U s h a p e d }) performs better than the msBayes model, but not as well as the models with alternative priors on divergence models. This suggests, as predicted by Oaks et al. [7, 15], that the tendency of msBayes to erroneously support models of temporally clustered divergences is caused by a combination of (1) small marginal likelihoods of models with more τ parameters due to uniform priors on divergencetime and demographic parameters and (2) the Ushaped prior on divergence models giving low prior density to models with intermediate numbers of divergence parameters. The former essentially rules out models with many τ parameters, which causes the latter to act like an "Lshaped" prior with a spike of prior density on the onedivergence model. Given the parameter richness of the model and the relatively small amount of information in the summary statistics, it is not surprising that the combination of these two factors can create a strong tendency to infer clustered models of divergence.
While the modifications improve the behavior of the model, I urge caution when using the method and interpreting its results. The method attempts to approximate the posterior of a very parameterrich model using relatively little information from the data. For example, when applied to the dataset of 22 taxon pairs from the Philippines [7], the model has as many as 604–625 free parameters (depending on τ), and samples over 1002 unordered divergence models. Even under the simplest possible model allowed under the new implementation, the model still has 471–492 free parameters. Furthermore, the stochastic coalescent and mutational processes being modeled predict a large amount of variation in possible datasets even when the parameter values are known. The richness and stochastic nature of the model makes for a difficult inference problem, especially when using a small number of summary statistics calculated from the sequence alignments of each taxon pair. The populationgenetic summary statistics used by the method contain little information about many of the free parameters in the model. Thus, I expect the improved method will still be sensitive to priors, and the power, while improved, may still be low. While there is much less prior sensitivity under the new model compared to those observed by Oaks et al. [7], there is still an effect when comparing the results of the empirical data analyzed under a diffuse (M_{ D P P }) and informative \left({\mathbf{M}}_{\mathit{\text{DPP}}}^{\mathit{\text{inform}}}\right) divergencetime prior (Figure 7C versus D). The fact that the posterior shifts toward the prior under the informative prior suggests that the shift away from the prior toward fewer divergence events under the diffuse prior might still be caused by small marginal likelihoods of models with more divergencetime parameters (Figure 7).
Nonetheless, it is reassuring to see a large amount of posterior uncertainty when the new implementation is applied to the empirical datasets (Figures 7 and 8). Applications of the msBayes model often result in strong posterior support for estimated scenarios (e.g., [3, 5–12]), as I found here (Figure 7). Given the richness of the model, the variance of the processes being modeled, and the relatively small amount of information in the summary statistics calculated from the sequence data, finding strong posterior support for any scenario is unexpected. Based on results of the empirical and power analyses (Figures 4, 6, 7 and 8), the new implementation more accurately reflects posterior uncertainty and avoids spurious support for biogeographical scenarios.
I also urge caution when using dppmsBayes due to the lack of theoretical validation of Bayesian model choice when the full data are replaced by summary statistics that are insufficient for discriminating across models under comparison [44], which is certainly the case here. Robert et al. [44] demonstrated that ABC estimates of model posterior probabilities can be inaccurate when such acrossmodel insufficient statistics are used.
Given all of these caveats, I encourage investigators to view this method as a means of exploring their data for general temporal patterns of divergences across taxa, rather than a rigorous means of evaluating hypotheses. As recommended by Oaks et al. [7], any results from the method should be accompanied by (1) analyses under a variety of priors to assess the assumptions underlying model inference and the prior sensitivity of the results, and (2) simulationbased power analyses to provide insight into the temporal resolution of the method. Both approaches are important to help guide the interpretation of results.
Given the difficulty of this estimation problem, I anticipate that fulllikelihood methods that can leverage all of the information present in the sequence data will become increasingly important for robustly estimating shared evolutionary history across taxa [45]. With improving numerical methods for sampling over models of differing dimensionality [46, 47], advances in Monte Carlo techniques [48], and increasing efficiency of likelihood calculations [49], analyzing rich comparative phylogeograpical models in a fulllikelihood Bayesian framework is becoming computationally practical, especially when considering that simulating millions of random datasets from the prior under the simple ABC rejection approach is inefficient and computationally nontrivial.
Conclusions
I introduced a new model for estimating shared divergence histories across taxa from DNA sequence data within an approximateBayesian modelchoice framework. The new method, dppmsBayes, takes a nonparametric approach to the problem by using a Dirichletprocess prior on the temporal distribution of divergences across taxa. The new method shows improved robustness, accuracy, and power compared to the existing method, msBayes. Compared to msBayes, the new approach better estimates posterior uncertainty, which greatly reduces the chances of incorrectly estimating biogeographical scenarios of shared divergence events. This is important, because models of shared divergence events are often ofparticular interest to researchers who employ these methods. This new tool will allow evolutionary biologists to better leverage comparative genetic data to assess the affects of regional and global biogeographical processes on biodiversity.
References
Hudson RR: Gene genealogies and the coalescent process. Oxf Surv Evol Biol. 1990, 7 (1): 144.
Wakeley J: Coalescent Theory: An Introduction. 2009, Greenwood Village, Colorado, USA: Roberts and Company Publishers
Hickerson MJ, Stahl EA, Lessios HA: Test for simultaneous divergence using approximate Bayesian computation. Evolution. 2006, 60 (12): 24352453.
Huang W, Takebayashi N, Qi Y, Hickerson MJ: MTMLmsBayes: Approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity. BMC Bioinformatics. 2011, 12: 1doi:10.1186/14712105121
Crews SC, Hickerson MJ, Leaché AD: Two waves of diversification in mammals and reptiles of Baja California revealed by hierarchical Bayesian analysis. Biol Lett. 2007, 3 (6): 646650. doi:10.1098/rsbl.2007.0368
Stone GN, Lohse K, Nicholls JA, FuentesUtrilla P, Sinclair F, Schönrogge K, Csóka G, Melika G, NievesAldrey JL, PujadeVillar J, Tavakoli M, Askew RR, Hickerson MJ: Reconstructing community assembly in time and space reveals enemy escape in a Western Palearctic insect community. Curr Biol. 2012, 22 (6): 532537.
Oaks JR, Sukumaran J, Esselstyn JA, Linkem CW, Siler CD, Holder MT, Brown RM: Evidence for climatedriven diversification? a caution for interpreting ABC inferences of simultaneous historical events. Evolution. 2013, 67 (4): 9911010. doi:10.1111/j.15585646.2012.01840.x
Barber BR, Klicka J: Two pulses of diversification across the Isthmus of Tehuantepec in a montane Mexican bird fauna. Proc R Soc BBiol Sci. 1694, 277: 26752681. doi:10.1098/rspb.2010.0343
Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C: Stability predicts genetic diversity in the Brazilian Atlantic forest Hotspot. Science. 2009, 323 (5915): 785789. doi:10.1126/science.1166955
Chan LM, Brown JL, Yoder AD: Integrating statistical genetic and geospatial methods brings new power to phylogeography. Mol Phylogenet Evol. 2011, 59 (2): 523537. doi:10.1016/j.ympev.2011.01.020
Plouviez S, Shank TM, Faure B, DaguinThiebaut C, Viard F, Lallier FH, Jollivet D: Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South. Mol Ecol. 2009, 18 (18): 39033917. doi:10.1111/j.1365294X.2009.04325.x
Voje KL, Hemp C, Saetre GP, Stenseth NC, Flagstad Ø: Climatic change as an engine for speciation in flightless Orthoptera species inhabiting African mountains. Mol Ecol. 2009, 18 (1): 93108. doi:10.1111/j.1365294X.2008.04002.x
Jeffreys H: Theory of Probability. 1939, Oxford, UK: Clarendon Press
Lindley DV: A statistical paradox. Biometrika. 1957, 44: 187192.
Oaks JR, Linkem CW, Sukumaran J: Implications of uniformly distributed, empirically informed priors for phylogeographical model selection: A reply to Hickerson et al. 2014. arXiv:1402.6397 [qbio.PE] http://arxiv.org/abs/1402.6397.
Hickerson MJ, Stone GN, Lohse K, Demos TC, Xie X, Landerer C, Takebayashi N: Recommendations for using msbayes to incorporate uncertainty in selecting an ABC model prior: A response to Oaks et al. Evolution. 2014, 68 (1): 284294. doi:10.1111/evo.12241
Hasegawa M, Kishino H, Yano TA: Dating of the humanape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985, 22 (2): 160174.
Felsenstein J: Evolutionary trees from DNA sequences: A maximum likelihood approach. J Mol Evol. 1981, 17: 368376.
Kingman JFC: The coalescent. Stochastic Process Appl. 1982, 13: 235248.
Bell ET: Exponential numbers. Am Math Month. 1934, 41: 411419.
Sloan NJA: The online encyclopedia of integer sequences, sequence A000041. http://oeis.org/A000041.
Sloan NJA: The online encyclopedia of integer sequences, Sequence A008284. http://oeis.org/A008284.
Malenfant J: Finite, closedform expressions for the partition function and for Euler, Bernoulli, and Stirling numbers. 2011, arXiv:1103.1585v6 [math.NT] http://arxiv.org/abs/1103.1585.
Ferguson TS: A Bayesian analysis of some nonparametric problems. Ann Stat. 1973, 1 (2): 209230.
Antoniak CE: Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Ann Stat. 1974, 2 (6): 11521174.
Lartillot N, Philippe H: A Bayesian mixture model for acrosssite heterogeneities in the aminoacid replacement process. Mol Biol Evol. 2004, 21 (6): 10951109. doi:10.1093/molbev/msh112
Huelsenbeck JP, Andolfatto P: Inference of population structure under a Dirichlet process model. Genetics. 2007, 175 (4): 17871802. doi:10.1534/genetics.106.061317
Huelsenbeck JP, Suchard MA: A nonparametric method for accomodating and testing acrosssite rate variation. Syst Biol. 2007, 56 (6): 975987. doi:10.1080/10635150701670569
Larget B, Baum DA, Smith SD, Rokas A, Ané C: Bayesian estimation of concordance among gene trees. Mol Biol Evol. 2007, 24 (2): 412426. doi:10.1093/molbev/msl170
Heath TA, Holder MT, Huelsenbeck JP: A Dirichlet process prior for estimating lineagespecific substitution rates. Mol Biol Evol. 2011, 29 (3): 939955. doi:10.1093/molbev/msr255
Heath TA: A hierarchical Bayesian model for calibrating estimates of species divergence times. Syst Biol. 2012, 61 (5): 793809. doi:10.1093/sysbio/sys032
Escobar MD, West M: Bayesian density estimation and inference using mixtures. J Am Stat Assoc. 1995, 90 (430): 577588.
Tajima F: Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983, 105 (2): 437460.
Watterson GA: On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975, 7 (2): 256276.
Takahata N, Nei M: Gene genealogy and variance of interpopulational nucleotide differences. Genetics. 1985, 110 (2): 325344.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123 (3): 585595.
Beaumont M, Zhang W, Balding DJ: Approximate Bayesian computation in population genetics. Genetics. 2002, 162: 20252035.
Blum MGB, François O: Nonlinear regression models for Approximate Bayesian Computation. Stat Comput. 2009, 20 (1): 6373.
Leuenberger C, Wegmann D: Bayesian computation and model selection without likelihoods. Genetics. 2010, 184: 243252. doi:10.1534/genetics.109.109058
Wegmann D, Leuenberger C, Neuenschwander S, Excoffier L: ABCtoolbox: a versatile toolkit for approximate bayesian computations. BMC Bioinformatics. 2010, 11: 116doi:10.1186/1471210511116
Nei M, Li WH: Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc Nat Acad Sci. 1979, 76 (10): 52695273. doi:10.1073/pnas.76.10.5269
Huelsenbeck JP, Rannala B: Frequentist properties of Bayesian posterior probabilities of phylogenetic trees under simple and complex substitution models. Syst Biol. 2004, 53 (6): 904913.
Oaks JR, Sukumaran J, Esselstyn JA, Linkem CW, Siler CD, Holder MT, Brown RM: Evidence for climatedriven diversification? a caution for interpreting ABC inferences of simultaneous historical events. Dryad Digital Repository. 2012, doi:10.5061/dryad.5s07m
Robert CP, Cornuet JM, Marin JM, Pillai NS: Lack of confidence in approximate Bayesian computation model choice. Proc Nat Acad Sci. 2011, 108 (37): 1511215117.
Sukumaran J: Geographies and genealogies: Phylogeographic simulation and Bayesian approaches to statistical phylogeographic model selection. PhD thesis,. 2012, University of KansasLawrence, Kansas, USA
Green PJ: Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995, 82 (4): 711732.
Lemey P, Rambaut A, Drummond AJ, Suchard MA: Bayesian phylogeography finds its roots. PLoS Comput Biol. 2009, 5 (9): 1000520
Sankararaman S, Jordan MI, BouchardCôté A: Phylogenetic inference via sequential Monte Carlo. Syst Biol. 2012, 61 (4): 579593. doi:10.1093/sysbio/syr131
Ayres DL, Darling A, Zwickl DJ, Beerli P, Holder MT, Lewis PO, Huelsenbeck JP, Ronquist F, Swofford DL, Cummings MP, Rambaut A, Suchard MA: BEAGLE: an application programming interface and highperformance computing library for statistical phylogenetics. Syst Biol. 2012, 61 (1): 170173.
Acknowledgements
I am greatful to Melissa Callahan, Mark Holder, Emily McTavish, Daniel Money, Jordan Koch, Adam Leaché, Peter Foster, two anonymous reviewers, and Christian Robert and his blog (http://xianblog.wordpress.com/) for insightful comments that greatly improved this work. I thank the National Science Foundation for supporting this work (DEB 1011423 and DBI 1308885). This work was also supported by the University of Kansas (KU) Office of Graduate Studies, Society of Systematic Biologists, Sigma Xi Scientific Research Society, KU Department of Ecology and Evolutionary Biology, and the KU Biodiversity Institute. I also thank Mark Holder, the KU Information and Telecommunication Technology Center, KU Computing Center, and the iPlant Collaborative for the computational support necessary to conduct the analyses presented herein.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The author declare that he has no competing interests.
Authors’ contributions
All aspects of this work were done by JRO.
Electronic supplementary material
12862_2014_2640_MOESM1_ESM.pdf
Additional file 1: Supporting table and figures. PDF of supporting Table S1 and Figures S1S33. As referenced in the main text. (PDF 8 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Oaks, J.R. An improved approximateBayesian modelchoice method for estimating shared evolutionary history. BMC Evol Biol 14, 150 (2014). https://doi.org/10.1186/1471214814150
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1471214814150