BMC Evolutionary Biology BioMed Central Methodology article

Background Marine allopatric speciation is an enigma because pelagic larval dispersal can potentially connect disjunct populations thereby preventing reproductive and morphological divergence. Here we present a new hierarchical approximate Bayesian computation model (HABC) that tests two hypotheses of marine allopatric speciation: 1.) "soft vicariance", where a speciation involves fragmentation of a large widespread ancestral species range that was previously connected by long distance gene flow; and 2.) peripatric colonization, where speciations in peripheral archipelagos emerge from sweepstakes colonizations from central source regions. The HABC approach analyzes all the phylogeographic datasets at once in order to make across taxon-pair inferences about biogeographic processes while explicitly allowing for uncertainty in the demographic differences within each taxon-pair. Our method uses comparative phylogeographic data that consists of single locus mtDNA sequences from multiple co-distributed taxa containing pairs of central and peripheral populations. We use the method on two comparative phylogeographic data sets consisting of cowrie gastropod endemics co-distributed in the Hawaiian (11 taxon-pairs) and Marquesan archipelagos (7 taxon-pairs). Results Given the Marquesan data, we find strong evidence of simultaneous colonization across all seven cowrie gastropod endemics co-distributed in the Marquesas. In contrast, the lower sample sizes in the Hawaiian data lead to greater uncertainty associated with the Hawaiian estimates. Although, the hyper-parameter estimates point to soft vicariance in a subset of the 11 Hawaiian taxon-pairs, the hyper-prior and hyper-posterior are too similar to make a definitive conclusion. Both results are not inconsistent with what is known about the geologic history of the archipelagos. Simulations verify that our method can successfully distinguish these two histories across a wide range of conditions given sufficient sampling. Conclusion Although soft vicariance and colonization are likely to produce similar genetic patterns when a single taxon-pair is used, our hierarchical Bayesian model can potentially detect if either history is a dominant process across co-distributed taxon-pairs. As comparative phylogeographic datasets grow to include > 100 co-distributed taxon-pairs, the HABC approach will be well suited to dissect temporal patterns in community assembly and evolution, thereby providing a bridge linking comparative phylogeography with community ecology.


Background
Allopatric speciation is an enigma in many marine organisms because larval dispersal can potentially connect disjoint populations and thereby prevent the reproductive and morphological divergence that arises from prolonged isolation [1][2][3][4][5][6][7]. This is especially enigmatic in the Indo-Pacific region where many species range freely across this expanse without evidence for barriers to genetic exchange. This marine region harbors the planet's highest species diversity and endemism of marine fauna, and with the absence of explicit barriers, some have pushed controversial models of sympatric speciation to explain this elevated diversity [8,9]. Even the proposed competing models of geographic speciation in the Indo-Pacific remain contentious and generally involve different facets of the classic dispersal or vicariance models for speciation. Under the one model (the "soft vicariance" model), speciations in peripheral archipelagoes result from a large widespread patchy ancestral species range connected by long distance gene flow that is eventually interrupted by oceanographic changes in temperature, sea level and/or currents [10][11][12] which leads to peripheral isolation and endemism. Under a second model (the "colonization" model), speciations in peripheral (i.e peripatric) archipelagoes emerge from sweepstakes centrifugal colonizations from high-diversity central areas followed by prolonged periods of isolation with potential inward range shifts towards the central region [8,9,13,14]. As in the case of terrestrial systems, using genetic data to distinguish these two scenarios is difficult because their expected genetic signatures are often similar, a situation that is exacerbated if demographic changes such as expansions and bottlenecks occur after an isolating event [10,[15][16][17].
Likewise, discerning the modes of isolation and speciation using phylogenetic and phylogeographic data is often fraught with uncertainty because species can potentially shift their ranges [18] or loose the population genetic patterns associated with colonization >> 2N generations subsequent to isolation. Traditionally, vicariance and dispersal histories have been tested using phylogenetic approaches that use area cladograms [19,20], consensus methods [21], or parsimony [22] in combination with some method of ancestral character state reconstruction. Although many of these classic methods were biased to find vicariance, recent methods incorporate more complex biogeography histories [23,24] such as maximum likelihood [25,26] and Bayesian methods that use both distributional and phylogenetic data [27]. Likewise, empirical studies have increasingly found dispersal/colonization to be a more common force behind allopatric speciation [15,16,[28][29][30][31]. Regardless, such analyses are often circular or ambiguous [32], and ancestral character reconstruction methods are always going to be hindered when elevated homoplasy in biogeographic patterns obscures the inferences in the older parts of a phylogeny [18,31,33].
Here we present an entirely different approach to testing for vicariance and dispersal histories. Instead of using phylogenetic comparative methods, we use coalescent population genetics to estimate ancestral demographic patterns across co-distributed taxa within a community. While this is in the spirit of previous suggested approaches that blend systematics and population genetics [34], the hierarchical Bayesian approach presented here tests vicariance and dispersal across taxon-pairs instead of doing so one at a time. Specifically we extend a hierarchical approximate Bayesian model (HABC) [35,36] in order to quantify the strength of these two alternative models of allopatric isolation across marine endemic taxa that are co-distributed in peripheral archipelagoes.
Using HABC allows sidestepping the requirement of an explicit likelihood function. Instead, it uses a probabilistic simulation model to generate data sets to compare with the empirical data. By using summary statistics, one can easily compare the simulated and empirical data in order to estimate parameters of the simulation model via an approximate sample of the posterior distribution. In HABC we use hyper-parameters that describe processes across co-distributed taxon-pairs as well as sub-parameters that describe the demographic history of each taxonpair.
First we describe the population genetic models whose hyper-parameters we want to estimate, and then we describe HABC. After detailing the HABC model, the summary statistics and the HABC implementation, we test these two biogeographic hypotheses given two comparative phylogeographic datasets: mtDNA CO1 data collected from multiple cowrie gastropod species that are endemic to the Hawaiian and Marquesan archipelagos. Specifically we use HABC to test whether marine vicariance ("soft vicariance") (H 1 ) or colonization (H 2 ) is the dominant isolating mechanism in either of these two marine communities. After using HABC to choose the best model of community isolation, we then use HABC to estimate temporal congruence in soft vicariance and/or colonization. We specifically use mtDNA sequence data collected from each co-distributed peripheral endemic taxon as well as each of the respective sister species which are usually more geographically widespread (Additional file 1).
Although this method specifically addresses questions relevant to the species diversity and patterns of endemism in the Indo-Pacific, it will be broadly applicable to many comparative phylogeographic datasets and biogeographic settings.

Soft vicariance and colonization
Rather than classical terrestrial vicariance where a large ancestral population is broken up into two isolated sister populations ( Figure 1A), our "soft vicariance" scenario (H 1 ; Figure 1B) has two ancestral populations with effective sizes (θ τ ) 1 and (θ τ ) 2 that are connected by high to moderate gene flow (M 1 = 1.0 to 100.0 migrants per generation) until τ V , when M 1 decreases to 0.0 -1.0 migrants per generation (M 2 ). If this second period is prolonged, then effective isolation and divergence can occur [37]. At τ V , the sizes of the two sister populations ((θ τ ) 1 and (θ τ ) 2 ) remain the same size or begin to grow exponentially until they reach their present sizes (θ 1 and θ 2 ) at τ = 0 depending on the draw from the prior ( Figure 1B). As in [16,17], time of vicariance, population sizes and migration rates are all free to vary across taxon-pairs according to their prior distributions.
Under the colonization scenario (H 2 ; Figure 1C), one of the sister populations is founded by a very small number of individuals (θ τ ) 2 that come from a larger source population (θ τ ) 1 at the time of colonization, τ C with subsequent isolation. The small colonizing population (θ τ ) 2 , then grows exponentially until it reaches its present effective size of θ 2 at τ = 0. The primary parametric expectation that distinguishes marine vicariance (H 1 ) from colonization (H 2 ) is the relatively small effective population size of the colonized population ((θ τ ) 2 ) at the putative time of colonization (τ C ). Secondarily, the possibility of gene flow subsequent to the vicariance event τ V further distinguishes H 1 and H 2 , although this assumption can potentially be relaxed. With regards to different patterns in the molecular genetic data under H 1 and H 2 , under colonization (H 2 ) samples from peripheral populations will likely accumulate a surplus of rare alleles due to having a current effective population size that greatly expanded from a small size after the colonization time τ C . In addition, there is likely to be generally more genetic diversity under soft vicariance (H 1 ) due to there being two ancestral populations rather than just one.
To statistically quantify the relative support of these two hypotheses (H 1 and H 2 ) across Y co-distributed peripheral endemics and their sister taxa given DNA sequence data, we extend and modify the hierarchical approximate Bayesian computation (HABC) framework of [35,36]. We also use this framework to estimate the temporal congruence of both vicariance and colonization across the Y phylogeographic data sets.
Although one could independently estimate (θ τ ) 2 in each of the Y data sets and use all of the independent posterior densities of (θ τ ) 2 to measure the support of H 1 and H 2 across the Y pairs, implementation of a hierarchical model accomplishes this from a single analysis and uses more information from the data via "borrowing strength" [38][39][40].

Hierarchical approximate Bayesian computation
Our implementation of HABC is based on the framework presented in [35,36,41], and we review the important features here. In HABC, sub-parameters (Φ; within taxonpair parameters) are conditional on "hyper-parameters" (φ) that quantify the variability of Φ among the Y taxonpairs. 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 a coalescent model using candidate parameters randomly drawn from the joint hyper-prior and sub-prior distribution P(φ, Φ). A summary statistic vector D i for each simulated dataset is then compared to the observed summary statistic vector D* in order to generate random observations from the joint posterior distribution f(φ i , Φ i |D i ) by way of a rejection/acceptance algorithm followed by a weighted local linear regression step [42].
The rejection/acceptance algorithm involves calculating a summary statistic vector from the observed data and each of the K simulated data sets. Each simulated data set is generated using parameters that are randomly drawn from the joint prior. Following [35,42], K Euclidian distances between the normalized observed summary statistic vector D* and each of the K normalized summary statistic vectors are then calculated (||D i -D*|| = d). An arbitrary proportion (tolerance) of the K simulations with the lowest d values are then used to obtain an approximate sample from the joint posterior after weighting and transforming the accepted parameter values using local linear regression [35,42]. After the local linear regression step, accepted and transformed parameter values that fall outside their respective prior bounds are subsequently transformed to have the values of their respective prior boundaries. For example, if an accepted parameter value with a uniform prior of [0.0,1.0] is transformed by local linear regression to a negative number, it is subsequently transformed to 0.0.

Hierarchical Model of Community Colonization and Vicariance
Model hyper-parameters and sub-parameters are listed and described in Additional file 2. Three hyper-parameters are drawn from their respective hyper-prior distributions (Additional file 2A), and these include: 1.) Z, the number of descendent populations per Y taxon-pairs that arise by colonization at times ; 2.) the number of different vicariance times Three models of allopatric isolation In addition to hyper-parameter estimation, we also use the HABC algorithm to sample from the posterior distributions of sub-parameter summaries (Additional file 2E) in order to quantify the support for H 1 and H 2 and secondarily estimate levels of temporal congruence in colonization and/or soft vicariance. Namely, we obtain estimates of the arithmetic means of three sub-parameters (E((θ τ )2), E(τ C ), E(τ V )), as well as the dispersion indexes of τ C and τ V (Ω C and Ω V respectively). The sub-parameter summary E((θ τ ) 2 ) is expected to be < 0.05 if H 2 is dominate across the Y taxon-pairs (Z = Y). The dispersion indexes Ω C and Ω V of the Z colonization times ( ) and the (Y -Z) vicariance times ( ) measure the ratio of the variance to the mean of these two sets of times and are therefore expected to be ≈ 0.0 when if there is temporal congruence in colonization or soft vicariance.

Two Stage Model Implementation
We implement the hierarchical analysis in two stages. In stage 1, an unconstrained general model is used to quantify the support for H 1 and H 2 across the Y taxon pairs. This first stage is accomplished by simulating K random draws from a general hyper-prior and using the HABC algorithm to sample from the posteriors of? Z and E((θ τ ) 2 ). Under our hierarchical model, H 1 and H 2 are equally probable because Z is a hyper-parameter drawn from the discrete uniform hyper-prior distribution [0, Y].
In the stage 2 analysis, K random draws are taken from a constrained hyper-prior where Z is fixed to be the mode of its posterior distribution obtained in stage 1. There are equal numbers of hyper-prior draws (K) obtained in both

Summary Statistic Vector for HABC Acceptance/Rejection Algorithm
In order to implement the HABC procedure, we use two modified versions of the summary statistic vector D used in [35]. For the Marquesas summary statistic vector (D Marquesas ), we calculate eight summary statistics collected from each taxon-pair (56 total). This includes π b (average pairwise differences between each central and peripheral Marquasan taxon-pair), π (average pairwise differences among all individuals within each taxon-pair), π w (average pairwise differences within descendent populations of each taxon-pair), θ W (Watterson's estimator of θ of each taxon-pair), and Var(π -θ W ). For π, θ W and Var(π -θ W ), each of the Y taxon-pairs are treated as a single population sample ( ). The Marquesas summary statistic vector, D Marquesas also includes calculating π, θ W , and Var(πθ W ) in each of the Y peripheral Marquasan samples ( ) that are putatively colonized and we denote these as π 2 , (θ W ) 2 , and Var(π -θ W ) 2 . Under this scheme, the vector D Marquesas is where each of the Y rows correspond to the Y taxon-pairs (Y = 7) and the eight columns correspond to the eight summary statistic classes. After these 8 × Y summary statistics are calculated, we must choose a way to consistently order the Y rows within D Marquesas . Instead of consistently ordering the rows by each taxon-pair's sample size, we increase the efficiency of our HABC estimator by ordering the rows based on the taxon-pair's Tajima's D [43] calculated from the taxon-pair's peripheral population sample ( ). For example, row 1 would contain the eight summary statistics collected from the taxon-pair with the lowest Tajima's D, and the Y th row would contain the eight summary statistics collected from the taxon-pair with the highest Tajima's D.
The motive for this ordering procedure is to extract more information from the data with respect to the estimated hyper-parameters than would be accomplished by ordering consistently by sample size [35]. For an efficient HABC estimator, there should be a strong correlation between pair-wise differences in hyper-parameter values (i.e. E(θ τ ) 2 or Z) and Euclidian distances between corresponding pairs of summary statistic vectors from corresponding pairs of simulated data sets. If ordering by sample size rather than ranked values of an informative summary statistic, pair-wise values of Z or E(θ τ ) 2 are not predicted to correlate with Euclidian distances of D calculated from corresponding pairwise simulated data sets. This is because sample size has no bearing on how each of the Y taxon-pairs are assigned to histories H 1 and H 2 when drawing values of Z from the hyper-prior. On the other hand, ordering by an informative summary statistic will minimize Euclidian distances among data sets with equal or similar values of Z regardless of which of the Y taxon pairs were assigned histories H 1 and H 2 . The consequent improved accuracy in HABC estimation that results from this ordering procedure is based on the exchangeability of the Y rows within D Marquesa (D 1 ,...,D Y ). If Φ i and D i are invariant to the permutations of the indexes (1 ,..., Y) and the i th taxon-pair's sample size is unrelated to the expectation of its Φ i or D i , there is exchangeability in the model [44]. We order by the peripheral Tajima's D because it is a summary statistic that is predicted to be informative with respect to demographic parameter differences between histories H 1 and H 2 (i.e. the ratio of each taxon-pair's (θ τ ) 2 and θ 2 ).
For the analysis of the 11 Hawaiian taxon-pairs, we use a reduced summary statistic vector D Hawaii to avoid null values that would arise in population samples that only included one individual (Additional file 1B). The vector D Hawaii in this case is , and the rows were likewise sorted by Tajima's D calculated from the peripheral populations.

Two Comparative Phylogeographic Implementations
We use this HABC method on two comparative phylogeographic data sets of Pacific cowrie gastropods (Cypraeidae). The first dataset consists of seven sister taxon-pairs of cowries that each consist of a descendent endemic species or sub-species that is distributed within the peripherally located Marquesas archipelago as well each sister taxon that is more geographically widespread (Additional file 1A). The second dataset consists of eleven sister taxonpairs of cowries co-distributed within the peripherally located Hawaiian archipelago. Like in the former data set, each pair consists of a Hawaiian endemic and a more widespread sister (Additional file 1B). Both data sets consisted of 614 base pairs of the CO1 mtDNA locus collected from 2-93 individuals per taxon-pair (Additional file 1). The HABC procedure was implemented using a modified version of the MSBAYES comparative phylogeographic software pipeline [35,36] consisting of several C and R programs that are run with a Perl "front-end" and utilizes a finite sites version of Hudson's classic coalescent simulator [45]. For both analyses, 2,000,000 random draws were sampled from the hyper-prior and 1,000 -2,000 accepted draws were used to construct hyper-posterior samples (tolerance of 0.0005 and 0.001 respectively) using the HABC acceptance/rejection algorithm.
The prior bounds for hyper-parameters and sub-parameters are given in Additional file 2. To explore the sensitivity of using different prior assumptions, we used two We calculate Bayes factors to compare the relative hyperposterior support of either history (soft vicariance and colonization) being dominate across all Y taxon pairs (e.g Z = 0 or Z = Y) against all other scenarios including mixed scenarios. We accomplish this by comparing relative hyper-posterior support of these two scenarios while accounting for the relative hyper-prior support for these two scenarios [47]. To calculate this Bayes factor, we use an arbitrary partition of hyper-parameter space to deline-ate where H 1 or H 2 is dominate across all Y taxon-pairs. For example, to evaluate the evidence of colonization being dominate across all Y taxon pairs (Z = Y) against all other scenarios (Z <Y), the approximate Bayes factor B(Z = Y, Z <Y) is the ratio of the two approximate hyper-posteriors of these two scenarios divided by the ratio of the two hyper-priors of these two scenarios,

Marquesas Results
The HABC analysis yielded a hyper-posterior that strongly supports H 2 , where colonization histories ( In stage 2, where we constrained the hyper-parameter Z to be 7, the hyper-posterior best supported temporal concordance in colonization across all seven Marquesas endemics. In this case, estimates of Ψ C and Ω C were 1.00 (95% quantiles: 1.00 -2.14) and 0.00 (95% quantiles: 0.00 -0.17) respectively (Additional file 3A; Figure 2). The mean time of colonization E(τ C ) was 1.58 My ago if we assume a 1% divergence rate per My ( Figure 2D). The

Hawaii Results
In contrast to the Marquesas analysis, the hyper-posteriors were much more similar to the hyper-priors ( Figure 3). The less informative posteriors are consistent with the lower Hawaiian sample sizes (Additional file 1B) and consequence of using of a reduced number of summary statistics (D Hawaii ) for the HABC acceptance/rejection algorithm. Although the hyper-posterior of Z given the Hawaiian data suggests a mixed history of both colonization and soft vicariance (Figure 3), larger sample sizes will be required to verify this. This weak inference is demonstrated by Bayes factors giving weak support for vicariance or colonization across all 11 Hawaiian endemics. Specifically, the calculated Bayes factors B(Z = 0, Z > 0), B(Z = 11, Z < 11), B(E((θ τ ) 2 ) < 0.05, (E((θ τ ) 2 ) = 0.05) < 1.0) and B(E((θ τ ) 2 ) = 0.05, E((θ τ ) 2 ) < 0.05) < 1.0) where all weak (< 1.0). Although the Hawaiian data yielded weak inference, the credibility intervals for Z ranged from 0.00 to 9.22 (Additional file 3C and 3D), suggesting that the history of isolation likely involved both soft vicariance and colonization. Likewise, estimates of Ω C and Ω V did not suggest simultaneous vicariance or colonization, and likewise yielded less informative posteriors than obtained in the Marquesas analysis (Figures 3C &3D). As was the case of the Marquesas analysis, hyper-posterior estimates were not sensitive to tolerance, prior assumptions on the upper

Simulated Data Sets
One of the chief advantages of HABC and ABC methods is the ease at which one can evaluate the performance, bias, and precision of the estimator via simulations. Specifically, we can simulate pseudo-observed data sets with known hyper-parameter values and compare estimates with their true values. Even though the most time-consuming task is to simulate a large enough prior and/or hyper-prior in ABC and HABC, once it is produced for the analysis of the real observed data set, it can subsequently be used to quantify bias and precision on the pseudoobserved data sets.
To this end, we simulate pseudo-observed data sets with known values of Z and E((θ τ ) 2 ) under the general model (stage 1), and Ψ C , Ψ V , E(τ C ), E(τ V ), Ω C , and Ω V under the constrained model (stage 2). We repeat this for both sample sizes used in the empirical implementation (D Marquesas and D Hawaii ). To simulate a data set (pseudo-observed data), all hyper-parameters and sub-parameters were randomly drawn from the prior and a corresponding D Marquesas or D Hawaii was subsequently calculated. For each corresponding pseudo-observed D Marquesas and D Hawaii , a hyper-posterior sample was obtained from the HABC rejection-sampling algorithm given 2,000,000 random draws from the hyper-prior. For every pseudo-observed (simulated) data set, an estimate is obtained from the posterior mode. We report estimates using a tolerance of 0.001 corresponding to 2,000 accepted draws from the hyper-prior. For both D Marquesas and D Hawaii and stage 1 and two, 250 estimates were made from 250 simulated pseudo-observed datasets.
In addition to these evaluations, we assess the ability to distinguish H 1 and H 2 when the true history is a special asymmetrical case of soft vicariance (H 1 ). To this end, we obtain HABC estimates of Z and E((θ τ ) 2 ) on 250 pseudoobserved datasets of size D Marquesas simulated under H 1 where true values of Z and E((θ τ ) 2 ) are 0 and 0.0 -0.05 respectively. Estimates for Z and E((θ τ ) 2 ) were obtained using a general prior (stage 1).
Another factor we explored was three other methods of post-acceptance transformation other than local linear regression (LLR) for obtaining a sample of the hyper-posterior distribution of the Z hyper-parameter. Because Z is a discrete integer (ranging from 0 to Y), it might be most appropriate to preserve Z as a discrete integer when using an ABC regression technique that implements polychotomous logit regression (PLR) [48][49][50]. We therefore used simulations to compare the effectiveness of LLR with: 1. PLR; 2. using the raw accepted values (RAW); and 3. a cumulative logit regression model (CLR). To compare the estimator bias and precision of these four methods, we calculated the root mean square error (RMSE) from 1000 estimates using each of these four methods on 1000 simulated pseudo-observed data sets with parameters randomly drawn from the priors. The methods for CLR and For each estimate, tolerance was 0.001 (2,000 accepted draws) using the local regression algorithm.

Results of Simulation Testing
Our simulations identified conditions under which the HABC estimator is reliable as well as conditions under which it is less reliable. At stage 1 of the HABC analysis, the estimates of E((θ τ ) 2 ) were consistently close to the cor-responding true values of E((θ τ ) 2 ), although the larger sample sizes matching the Marquesas data set (D Marquesas ) yielded more accurate estimates of E((θ τ ) 2 ) than using sample sizes matching the smaller Hawaii data set (D Hawaii ; Figures 4A &4B). On the other hand, estimates of Z were less accurate, yet using D Marquesas resulted in more accurate estimates of Z than when using D Hawaii ( Figures  4C &4D). Additionally, estimates of Z representing colonization across all Y taxon-pairs never resulted in false positives given D Marquesas ( Figure 4C). Specifically, when estimates of Z were 6 or 7, true Z values were always 6 or 7. Conversely, lower true values of Z yielded less accurate estimates of Z (Figures 4C &4D). For example, when estimates of Z were 0, true values ranged from 0 -4 given D Marquesas (Figure 4C). In general, the simulation analysis verified the statistical confidence in the empirical inference of colonization across all seven Marquesas endemics, yet demonstrated there to be more uncertainty in the inferred history of the 11 Hawaiian endemics.
Under the constrained stage 2 model, the simulations revealed a strategy for estimating the variability in the Z colonization times. Specifically, the best strategy would be to use estimates of Ω C rather than Ψ C ( Figure 5). However, there was less accuracy and precision in the Ω C estimates given D Hawaii ( Figure 5E). Unlike estimates of Ω C , estimates of Ω V or E(τ V ) were not as reliable ( Figure 5H), perhaps due to the very small amount of migration (0 to 1.0 migrants per generation) after each "soft vicariance" event (τ V ) under H 1 .
Estimator Performance: 250 true hyper-parameter values plotted against their posterior mode estimates (stage 2 model) Estimator performance under constrained history of asymmetrical soft vicariance (True Z = 0; True E((θ τ ) 2 ) < 0.05; stage 1 model) Figure 6 Estimator The simulations also revealed that we have the ability to distinguish H 1 and H 2 when the true history is a special asymmetrical case of vicariance (H 1 ) where Z = 0 and each (θ τ ) 2 ranges from 0.0 to 0.05 under the D Marquesas sample size configuration ( Figure 6). In this case, 63% of the Z estimates were ≤ 1 (Figure 6A &6B). Likewise, estimates of E((θ τ ) 2 ) were a reliable indicator of detecting H 1 under this special asymmetrical case of vicariance (H 1 ). Even though true values of E((θ τ ) 2 ) ranged from 0.0 to 0.05, estimates of E((θ τ ) 2 ) in this case were upwardly bias, with 90% of the E((θ τ ) 2 ) estimates ranging from 0.21 to 0.48 ( Figure 6C). Even though the E((θ τ ) 2 ) estimator is upwardly biased given this special case, it is upwardly biased in the direction of correct

Model Assumptions and Robustness
While previous studies have used multi-locus population genetic data to reconstruct the demography and geography of speciation [51][52][53][54][55], here we use single locus mtDNA data to look at patterns across multiple co-distributed taxa. Although single locus inference can be hazardous in the face of coalescent variance and the possibility of selection, our approach offers the possibility to look at patterns of community assembly when the community consists of many non-model organisms where only "barcode" DNA sequence data can be feasibly collected. Not only does our model incorporate the stochasticity of single-locus coalescent variance across taxa, by combining datasets into a hierarchical Bayesian analysis we gain statistical "borrowing strength" [38]. The "borrowing strength" of HABC is achieved by making inferences across groups (i.e. co-distributed taxa) by pooling information across the groups without assuming the groups are from the same population [39,40]. This allows estimating congruence across groups in sub-parameters while borrowing strength from the full comparative phylogeographic sample. This "borrowing strength" translates into higher sample size depending on the magnitude of the "pooling factor" which represents the degree to which sub-parameter estimates (Φ) are pooled together from hyper-parameter estimates of φ, rather than estimated independently from each phylogeographic dataset [56].
The possibility of selection at the tightly linked mtDNA genome could bias results of our analytical method [57], especially if balancing selection occurred in the mtDNA genome in ancestral or descendent taxa such that coalescent events were much older than neutral expectations. Likewise, if positive selective sweeps occurred on the mtDNA genome after colonization or vicariance, estimates of Z and E((θ τ ) 2 ) could be biased to reflect colonization [58,59]. However, if positive selection only occurred at the mtDNA genome before colonization or vicariance, then the timing of these isolation events could be better estimated due to reduction of ancestral polymorphism. Furthermore, because selection and demography are ultimately confounding, our method will be less reliable if mitochondrial positive selection is prevalent in the comparative phylogeographic sample. This could be especially troublesome if peripatric speciation by colonization or vicariance involves positive selection at mtDNA genes that allow adaptive divergence in novel peripheral habitats [60][61][62]. Nonetheless, results from our HABC method can be considered conservative with respect to inferring the geographic and demographic history of isolation across a peripheral community. For example, a strong inference of colonization across an entire data set could result from both strong positive selection and/or small effective colonizing population sizes, but it would be extraordinary if strong positive selection occurred at the mtDNA genome of all Y co-distributed taxon-pairs.
A strong result of temporal concordance in isolation is also conservative with respect to violations from a uniform molecular clock model. If rate variation occurred across the Y taxa, then we would expect an inference in temporal discordance unless rates were inversely proportional with actual isolation times. Nevertheless, this HABC method will be most useful when applied to co-distributed taxon-pairs that are closely related. Our empirical application was restricted to cowrie gastropods (Cypraeidae) and the COI loci we used only marginally rejected rate constancy [63].
Although Bayesian methods are less robust if results are heavily dependent on prior bounds, results from the empirical data were not sensitive to our exploration of different prior assumptions. The overall inferences and hyper-posterior estimates from the empirical data were not sensitive to model assumptions regarding the priors of θ (Additional file 3) or post-isolation migration. Additionally, all four methods of post-acceptance transformation (LLR, PLR, CLR and RAW) yielded identical mode estimates of Z given the Marquesas data and similar mode estimates of Z given the Hawaii data (ranging from 4 -5).
Another consideration is how deviations from a panmictic Wright-Fisher model could have affected our HABC estimates. Although the sampling scale is large in some of the source species (Additional file 1), the cowrie gastropod taxa we included have high dispersal capabilities and are therefore not likely to have elevated within species subdivision [63]. This is confirmed by the relatively low levels of within species average pair-wise differences (π 1 and π 2 ) in both data sets (Additional file 1). If intra-species migration rates are > 1, our idealized coalescent model assuming intra-species panmixia is somewhat appropriate (Slatkin 1985). Even with some population some structure, a standard coalescent model can suffice if a species consists of many small demes with at least moderate migration such that number of demes is approximated by a scaled effective population size [64][65][66][67]. If ancestral population structure is approximately scaled by ancestral effective population size, then our chosen priors are conservative because we allow for ancestral population sizes that are two to four times as large as the observed pair-wise distances of extant population sizes (π 1 and π 2 ; Additional file 1). However, our method should not be applied to populations that are heavily structured over large geographic scales.

Vicariance and Dispersal in Marine Communities
Vicariance and dispersal speciation could be hugely relevant in the marine realm, especially within the highly diverse Indo-Pacific region that is dominated by islands rather than long continuous coastlines. Explanations for elevated Indo-Pacific diversity in the centrally located "coral triangle" portion of this Indo-Pacific region (Philippines, Malay Peninsula, and New Guinea) usually revolve around sympatric speciation followed by outward range shifts [68] or peripatric speciation followed by inward range shifts [69,70]. However, the plausibility of the first hypothesis of sympatric speciation is very controversial on theoretical grounds [71,72] as well as being very difficult to test empirically [18]. On the other hand, the second hypothesis of peripatric speciation is a much more likely force behind Indo-Pacific diversification if long distance oceanic dispersal to peripheral populations is sufficiently low for isolation and subsequent reproductive or ecological divergence to emerge between central and peripheral archipelagoes [1][2][3][4][5][6][7].
Instead of the classic vicariance model, under our marine vicariance model (or "soft vicariance") an ancestral range is inter-connected by long distance gene flow that is interrupted by oceanographic changes in temperature, sea level and/or currents [10][11][12]. In this case we might predict that co-distributed peripheral endemics became isolated simultaneously in taxa with lower dispersal capability. On the other hand, our marine colonization model is more similar to the classic dispersal or peripatric model. Here, allopatric isolation arises via sweepstakes colonization where the timing of colonization could be predicted to occur randomly across the co-distributed endemics after an archipelago emerges from geological processes.

Hawaiian vs Marquesas Dynamics
Both Hawaii and the Marquesas have some of the highest levels of endemism in all of Oceania [73][74][75], but HABC analyses support different histories for their endemic species. The Hawaiian archipelago is the most isolated island chain in the Indo-West Pacific and has existed for > 70 My with coral reefs since at least 35 My [76]. It is now well established that terrestrial endemics can be older than the oldest emergent island in the archipelago, a pattern resulting from initial colonization of older, now subsided islands, followed by dispersal to new islands after emergence [77]. Thus, there has been ample opportunity for isolation and speciation, perhaps even more so for marine taxa. Although the lower sample sizes lead to greater uncertainty associated with the Hawaiian estimates (Figure 3), the hyper-parameter 95% credibility intervals suggest a strong inference of isolation via soft vicariance in at least a subset of the Hawaiian taxon-pairs (Z = 0.0 -9.22; E((θ τ ) 2 ) = 0.30-0.95). If this is the case, then occupation of the Hawaiian archipelago was much older than the Marquesas archipelago, consistent with geologic evidence. Moreover, the inference of soft vicariance in a number of the taxon pairs suggests that there was greater potential for migrants between this archipelago and the central Pacific and Indo-Pacific triangle regions during older periods. Such connectivity of the Hawaiian chain to the remaining Indo-West Pacific via Johnston Atoll has been suggested in other studies [78][79][80].
In contrast, we find strong inference of isolation via colonization across all seven cowrie gastropod endemics codistributed in the Marquesas, as well as a strong inference of temporal concordance in this colonization. Unlike other island chains in the Pacific Ocean that have older seamounts and atolls trending away from most recent island (e.g. Hawaii), the Marquesan hotspot is unusual in being quite young. The ages of Marquesan islands range from 1.3 My (Fatu Hiva) to 6.0 My (Eiao) [81]. If we apply a molecular clock of 1% divergence per My [63], the inferred timing of simultaneous colonization of the Marquesas archipelago is from 0.84 -1.90 My (Additional file 3A), and is consistent with the young origins of the islands. If we accept our strong inference of temporal concordance, it could be argued that this assemblage of cowries colonized via an episodic oceanographic event that caused a surge in gene flow from the central Pacific region. Given that the Marquesas has one of the highest levels of marine endemism in Oceania [73,74], it will be interesting to see if HABC analyses on other taxon-pairs show similarly young divergences and temporal congruence.

Conclusion
Although soft vicariance and colonization are likely to produce relatively similar genetic patterns when only a single taxon-pair is considered, our simulation analysis shows that our hierarchical Bayesian model can potentially detect if either history is a dominant process across a marine community. The empirical implementation of our method yields a strong inference of isolation via simultaneous colonization across all seven cowrie gastropod endemics co-distributed in the Marquesas. In contrast, our method shows a strong inference of isolation via soft vicariance in at least a subset of the 11 Hawaiian taxon-pairs, although the smaller sample size resulted in less certainty in our estimates.
Our HABC method exemplifies the utility in "statistical phylogeographic" approaches [82,83] rather then qualitative and descriptive approaches that make large inferences from the small details observed from gene trees [84]. The HABC approach accomplishes this by analyzing all the phylogeographic datasets at once in order to make across taxon-pair inferences about biogeographic processes while explicitly allowing for uncertainty in the demographic differences within each taxon-pair.
Although the approach described here uses HABC to test for two particular biogeographic explanations of allopatric diversification across co-distributed taxa, the HABC framework is flexible and therefore can provide a skeleton for testing other biogeographic models from comparative phylogeographic data. Indeed, one of the original objectives of phylogeography comparative was to resolve deepseated questions about how climate change drives community assembly and evolution of whole biotas [85]. However, this goal has so far been unrealized [86][87][88] because comparative phylogeographic studies rarely involve more than a handful of co-distributed species. Comparative phylogeographic datasets are bound to have explosive growth as collecting DNA sequence data across a wide diversity of co-distributed taxa scales up to the level of comprehensive ecosystem sampling. Such "community-scale" comparative phylogeographic data sets could potentially test classic biogeographic hypotheses (e.g. vicariance versus dispersal) at the community level [89,90], as well as test controversial and fundamental hypotheses in community ecology such as Hubbell's Neutral theory [91], Tillman's stochastic competitive assembly model [92], and Diamond's niche assembly rules [93,94]. As comparative phylogeographic datasets grow to include > 100 co-distributed taxon-pairs, the HABC approach will be well suited to dissect temporal patterns in community assembly and thereby provide a bridge linking comparative phylogeography with community ecology.
anonymous reviewers and the communicating editor H. Zauner for useful comments and suggestions to improve the manuscript.