Testing for simultaneous divergence (vicariance) across different population-pairs that span the same historical barrier to gene flow is of central importance to evolutionary biology, biogeography and community ecology [

1-

4]. Such inferences inform processes underlying speciation, community composition, range delineation, and the ecological consequences of climatic changes. Estimating a population divergence time with an appropriate statistical model [

5] can be accomplished in a variety of ways [

6-

8], yet analyzing comparative phylogeographic data with multiple co-occurring species pairs that vary with respect to demographic parameters and pairwise coalescent times is less straightforward.

Instead of conducting an independent analysis on every population-pair and testing the hypothesis of temporal concordance based on this set of independent parameter estimates of divergence time, the hierarchical model employed by msBayes follows the suggestion of [

9] by concurrently estimating three hyper-parameters that characterize the mean, variability and number of different divergence events across a set of population-pairs. The model employed in msBayes allows estimation of these hyper-parameters across a multi-species data set while explicitly incorporating uncertainty and variation in the sub-parameters that independently describe each population-pair's demographic history (divergence time, current, ancestral and founding effective population sizes), post-divergence migration rate and recombination rate. The msBayes software pipeline is based on the introduction of the approximate Bayesian computation (ABC) method for sampling from the hyper-posterior distribution for testing for simultaneous divergence [

10]. We review the important features here. Although the current implementation is for a single locus per species-pair, future implementations will allow analysis of multi-loci data per species/population pair.

In contrast to previous ABC-like models [

11-

15], our TSD is accomplished by implementing a hierarchical Bayesian model in which the sub-parameters (Φ; within population-pair parameters) are conditional on "hyper-parameters" (

*ϕ*) that describe the variability of Φ among the

*Y *population-pairs. For example, divergence times (Φ) can vary across a set of population pairs conditional on the set of hyper-parameters (

*ϕ*) that varies according to their hyper-prior distribution. Instead of explicitly calculating the likelihood expression P(Data |

*ϕ*,Φ) to get a posterior distribution, we sample from the posterior distribution P((

*ϕ*,Φ) | Data) by simulating the data

*K *times under the coalescent model using candidate parameters drawn from the prior distribution P(

*ϕ*,Φ). A summary statistic vector

**D **for each simulated dataset is then compared to the observed summary statistic vector in order to generate random observations from the joint posterior distribution

*f*(

*ϕ*_{i},Φ

_{i}|

**D**_{i}) by way of a rejection/acceptance algorithm [

16] followed by an optional weighted local regression step [

15]. Loosely speaking, hyper-parameter values are accepted and used to construct the posterior distribution with probabilities proportional to the similarity between the summary statistic vector from the observed data and the summary statistic vector calculated from simulated data.

The hierarchical model consists of ancestral populations that split at divergence times *T*_{Y }= {*τ*_{1}...*τ*_{Y}} in the past. The hyper-parameter set, *ϕ *quantifies the degree of variability in these *Y *divergence times across the *Y *ancestral populations and their *Y *descendent population pairs: (1) Ψ, the number of possible divergence times (1 ≤ Ψ ≤ *Y*); (2) E(*τ*), the mean divergence time; and (3) Ω, the ratio of variance to the mean in these *Y *divergence times, Var(*τ*)/E(*τ*). The sub-parameters for the *i*-th population-pair (Φ_{i}) are allowed to vary independently across *Y *population pairs and include divergence time (*τ*_{i}), current population sizes, ancestral population sizes, post-divergence founding population sizes, durations of post-divergence population growth, recombination rates, and post-divergence migration rates. The multiple population-pair splitting model is depicted in Figure . Each divergence time parameter *τ *is scaled by 2*N*_{AVE }generations, where *N*_{AVE }is the parametric expectation of *N *(effective population size) across *Y *population pairs given the prior distribution.

The summary statistic vector

**D **employed in msBayes currently consists of up to six summary statistics collected from each of the

*Y *population pairs (

*π*,

*θ*_{W}, Var(

*π *-

*θ*_{W}),

*π*_{net},

*π*_{b}, and

*π*_{w}). This includes

*π*, the average number of pairwise differences among all sequences within each population pair,

*θ*_{W }the number of segregating sites within each population pair normalized for sample size, [

17], Var(

*π *-

*θ*_{W}) in each population pair, and

*π*_{net}, Nei and Li's net nucleotide divergence between each pair of populations [

18]. This last summary statistic is the difference (

*π*_{b }-

*π*_{w}) where

*π*_{b }is the average pairwise differences between each population pair and

*π*_{w }is the average pairwise differences within a sister pair of descendent populations. The default setting includes the first four aforementioned summary statistics because they were found to be a least correlated subset of a larger group [

19], however, future versions of msBayes will allow users to choose other summary statistics.

An extensive simulation study was conducted in [

10] to evaluate the performance of our hierarchical ABC model. Because comparative phylogeographic studies are often conducted on multi-species data sets that include rare taxa from which it is difficult to obtain samples from many individuals, we extend the previous evaluation to explore the effectiveness of msBayes in conducting a TSD given small sample sizes (≤ 5 individuals per population pair).