MTMLmsBayes: Approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity
 Wen Huang^{1}Email author,
 Naoki Takebayashi^{2},
 Yan Qi^{3} and
 Michael J Hickerson^{1, 3}
DOI: 10.1186/14712105121
© Huang et al; licensee BioMed Central Ltd. 2011
Received: 5 May 2010
Accepted: 3 January 2011
Published: 3 January 2011
Abstract
Background
MTMLmsBayes uses hierarchical approximate Bayesian computation (HABC) under a coalescent model to infer temporal patterns of divergence and gene flow across codistributed taxonpairs. Under a model of multiple codistributed taxa that diverge into taxonpairs with subsequent gene flow or isolation, one can estimate hyperparameters that quantify the mean and variability in divergence times or test models of migration and isolation. The software uses multilocus DNA sequence data collected from multiple taxonpairs and allows variation across taxa in demographic parameters as well as heterogeneity in DNA mutation rates across loci. The method also allows a flexible sampling scheme: different numbers of loci of varying length can be sampled from different taxonpairs.
Results
Simulation tests reveal increasing power with increasing numbers of loci when attempting to distinguish temporal congruence from incongruence in divergence times across taxonpairs. These results are robust to DNA mutation rate heterogeneity. Estimating mean divergence times and testing simultaneous divergence was less accurate with migration, but improved if one specified the correct migration model. Simulation validation tests demonstrated that one can detect the correct migration or isolation model with high probability, and that this HABC model testing procedure was greatly improved by incorporating a summary statistic originally developed for this task (Wakeley's Ψ_{ W } ). The method is applied to an empirical data set of three Australian avian taxonpairs and a result of simultaneous divergence with some subsequent gene flow is inferred.
Conclusions
To retain flexibility and compatibility with existing bioinformatics tools, MTMLmsBayes is a pipeline software package consisting of Perl, C and R programs that are executed via the command line. Source code and binaries are available for download at http://msbayes.sourceforge.net/ under an open source license (GNU Public License).
Background
Comparative phylogeographic inference of multiple codistributed taxa is of central importance in evolutionary biology, biogeography and community ecology [1–5]. Soon it will be routine to use large amounts of genetic data sampled from multiple individuals and multiple nonmodel taxa [6] in combination with other sources of environmental and ecological information to make powerful biogeographic inference such as how climate change affects whole biota or how geographic processes generate and partition patterns of biodiversity across communities [7]. However, simultaneous analysis of data from multiple taxa and multiple unlinked loci presents analytical and computational challenges. By utilizing simulation and summary statistics to avoid the need to calculate an explicit likelihood function, ABC (approximate Bayesian computation) or "likelihoodfree" methods can potentially tackle complex multitaxa demographic models when more exact methods are not efficient [8]. Although some information in the data is sacrificed when only using summary statistics, ABC methods have been shown to compare well against methods that utilize an explicit likelihood function [9–11] by allowing efficient extraction of information from the data under explicit models that can be built from background information [12–14].
Here we present MTMLmsBayes, a computer software pipeline that can be used to test for simultaneous divergence and migration across multiple codistributed taxonpairs given multilocus DNA sequence data. This method uses a coalescentbased model involving multiple taxa that can diverge at different times into taxonpairs and independently experience different magnitudes of migration, population sizechanges, and intragenic recombination. The hierarchical model allows for variation and uncertainty in demographic parameters across taxonpairs while testing specified multiple taxa scenarios of postdivergence migration and estimating hyperparameters that characterize the variability in divergence times across taxonpairs. Uncertainty in mutation rate heterogeneity across loci is also accounted for. For example, this software will allow testing for simultaneous divergence [11] and choosing among alternate multitaxon scenarios such as isolation and migration that can be generated via ecological niche models [15]. Some recent packages have recently made ABC methods accessible to empiricists conducting phylogeographic inference [16–21], and MTMLmsBayes complements these by using HABC for comparative phylogeographic datasets.
Hierarchical Bayesian models
The use of a hierarchical Bayesian framework within the context of ABC has been described elsewhere [10–12, 22, 23]. As in the single locus msBayes [24], our hierarchical Bayesian model includes taxonspecific demographic parameters and locusspecific mutation parameters (Φ) that are conditional on demographic and mutational "hyperparameters" (ϕ) which quantify the variability of Φ among the different taxonpairs and loci. These hyperparameters ϕ, can in turn be conditional on discrete "model indicator parameters" (ϕ^{ Z } ). For example, taxonspecific parameters (Φ) for migration rates can vary across a set of population pairs conditional on either hyperprior distributions ϕ^{ 1 }or ϕ^{ 2 }, which both represent different biogeographic hypotheses about the dynamics of isolation across codistributed taxonpairs. Potentially, Bayesian model choice can first be performed by obtaining the Bayesian posterior probabilities of models ϕ^{ 1 }, ..., ϕ^{ Zmax } and subsequently obtaining the posterior probabilities of other hyperparameters conditional on the model with highest posterior probability or averaged across models conditional on their relative posterior probabilities [22, 25].
Hierarchical ABC
Instead of explicitly calculating the likelihood expression P(Dataϕ^{ Z } , Φ) to get a joint posterior distribution, we sample from the joint posterior distribution P((ϕ^{ Z } , Φ)Data) by simulating the data under a coalescent model using candidate parameters drawn from the joint prior distribution P(ϕ^{ Z } , Φ). A summary statistic vector D of each simulated multitaxon multilocus dataset is then compared to the observed summary statistic vector D* in order to generate random observations from the joint posterior. MTMLmsBayes implements hierarchical ABC by way of a standard rejection/acceptance algorithm [10, 26–30] followed by an optional transformation step.
Specifically, for the simulated i th data set, a set of parameter values and Φ_{ i } are randomly drawn from their joint prior P(ϕ^{ Z } , Φ) and are then used to simulate data and associated D_{ i }. This is repeated until a large number of sample points from the joint prior distribution P(D,ϕ^{ Z } , Φ) have been obtained (typically 10^{6}  10^{7}). The joint posterior distribution for ϕ^{ Z } and Φ is sampled with probabilities proportional to the similarity between D* and each simulated sample of D_{ i }. The summary statistics within each vector D_{ i }are scaled to have unit variance followed by calculating a Euclidian distance between D_{ i }and D*. Subsequently, a userdefined proportion of simulated summary statistic vectors D_{ i }are accepted with their associated parameter values being used to sample the joint posterior. Typically 50010,000 simulated data sets are accepted out of > 10^{6} prior simulations. To improve upon the posterior estimation, an optional transformation step can involve local linear regression for continuous hyperparameters following the scheme of [31] or polychotomous logistical regression for estimating discrete model indicator parameters or discrete integer hyperparameters [25, 32, 33]. Alternatively, one could apply other postacceptance transformation methods [21, 34, 35] such as the general linear model [21].
Models of demography and DNA sequence evolution
As in [24], each taxon consists of an ancestral population of effective size θ_{ A } that splits at time τ into two descendent populations of effective sizes θ_{ A1 } and θ_{ A2 } which then start exponentially growing into sizes θ_{ B1 } and θ_{ B2 } at times τ_{ B1 } and τ_{ B2 } . If there is migration incorporated into the demographic model, each taxonpair has an effective migration rate that occurs after divergence (Nm; where m is the probability of symmetric migration between pairs of descendent sister populations). The parameters Nm, θ_{ A } , θ_{ A1 } , θ_{ A2 } , θ_{ B1 } , θ_{ B2 } , τ_{ B1 } and τ_{ B2 } all independently vary across all codistributed taxonpairs according to uniform prior distributions that are specified by the researcher.
The multiple loci from each taxonpair are assumed to be unlinked and for the mutation model, the JukesCantor [37], equalinput (F81) [38], or HKY model [39] of DNA substitution can be optionally used for all loci [37], with the equalinput model being default. Because the divergence with migration model is generally applied to taxa that have split in the last 10 My, these models should be sufficient [40]. The rate of DNA substitution can vary across unlinked loci such that the rate differences are scaled from the mean of a gamma distribution. Uncertainty in the shape parameter α, is incorporated by randomly drawing α from a uniform hyperprior distribution ranging between 1 and 20 with the scale parameter = 1/α, such that the mean rate scalar is 1.0 across replicate simulations. If there is prior evidence regarding specific patterns in rate variation amongst loci, such as relative rate estimates from samples of outgroup taxa, one can constrain loci to have specific average rate differences. Furthermore, a scalar parameter for each locus can incorporate ploidy differences for loci such as haploid uniparentally inherited mitochondrial and chloroplast DNA, and likewise these scalar parameters can enforce relative differences in generation times across taxonpairs. A uniform prior distribution can be optionally specified for intragenic recombination rates that vary independently across taxa.
Summary Statistics
The summary statistic vector D in MTMLmsBayes can calculate up to 23 summary statistic classes collected from each locus of every taxonpair. These summary statistic classes are of three categories: 1.) whole population summary statistics that treat the taxonpair as a single population sample; 2.) subpopulation summary statistics that are calculated on each of the two descendent population samples, and 3.) summary statistics that quantify differences between the two descendent population samples. Categories 1 and 2 include π, 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, [41], SD(π  θ_{ W } ) the standard deviation in the difference between these two quantities, sH, Shannon's diversity index on allele frequencies [42], and s, Wakeley's population correlation coefficient in the number of pairwise differences [43]. Category 3 includes, π_{ b } and π_{ net } , the total average and net average pairwise differences between two descendent population samples, [44], and Wakeley's s_{ XY } and Ψ_{ W } , two derivations of the interpopulation correlation coefficient in the number of pairwise differences between populations. These latter two summary statistics have been demonstrated as useful for distinguishing migration from isolation models [43, 45].
For every simulated dataset of multiple taxonpairs and multiple loci, a three dimensional vector (D) of these summary statistics is first constructed with dimensions of x summary statistic classes, y taxonpair indicator elements and z loci. Subsequently, a new 3dimensional vector D_{ m }is calculated from D where z_{ m } consists of the first four raw moments across loci [46]. The raw moments are the moments about zero, which can be converted to central moments (the first to forth central moments correspond to mean, variance, skewness, and kurtosis, respectively) using binomial transformation [47]. Moments are used to reduce the dimensionality of summary statistics vector and to capture the distribution of random variables (summary statistics) across loci. To be general, the number of sampled loci can vary amongst taxonpairs such that the length of z varies within D whereas within D_{ m }, z can have up to 4 elements.
When calculating D_{ m }, a final step involves reordering the taxonpair indicator elements of y into descending values of mean π_{ b } across loci such that each of the x columns of summary statistic classes have their y elements tandemly ordered by descending values of π_{ b } . This greatly reduces the combinatorial sample space such that order of the original sampling configuration is not determinant on any corresponding ordered vector of π_{ b } 's (which are predicted to correlate with the corresponding vector of τ's [48]). This strategy takes advantage of the exchangeability of the expected values of π_{ b } across sample sizes and their correlation with each taxonpair's τ (divergence time) [48]. By using this reordering procedure, the Euclidian distance between each simulated (D_{ m })_{ i }and observed D_{ m }* is independent of the ordering of taxonpairs within the sampling configuration and results in a higher correlation between ΔΩ and ΔD_{ m }than when not reordering. Here, ΔΩ is the difference in Ω (dispersion index of divergence times across Y taxon pairs) between pairs of data sets and ΔD_{ m }is the Euclidian distance between their corresponding pairs of summary statistic vectors D_{ m }that are calculated from these corresponding pairs of data sets. This ordering scheme for D_{ m }results in a desired ABC strategy with a higher correlation between summary statistics and estimated parameters (i.e. Ω and D_{ m }). This was confirmed by comparing pairs of simulated data sets and here we verify the improved performance of this sorting procedure via simulations.
Implementation
Running MTMLmsBayes is a four step process that includes: (1) preparation of the input file specifying prior distributions and the sampling configuration from the DNA sequence data; (2) preparation of the observed summary statistic vector, (3) generating a "reference table" of simulated data sets (i.e. coalescent simulations of data sets matching the observed data with regards to the sampling configuration and with parameters drawn from the prior); and (4) performing an acceptance/rejection step to obtain a sample from the posterior distribution. To improve estimation, accepted parameter values sampled from the posterior distribution can be subjected to transformation methods depending on if whether they are continuous parameters (local linear regression) or discrete parameters (polychotomous regression) using R scripts provided by M. Beaumont. Alternatively, one could perform other recently developed methods of postacceptance transformation to improve parameter estimation [21, 34, 35].
Due to the modular pipeline architecture of MTMLmsBayes, users can also opt to use available command line driven R scripts to generate pseudoobserved data sets (i.e. "PODS"; [14]) in order to conduct simulationbased model validation as well as fine tune the ABC conditions with respect to choice of summary statistics and acceptance threshold. In addition, users can use available R scripts to conduct posterior model fitting in order to assess the fit of the simulation models to the observed data [13, 14].
After installing the software via binary installation or compilation of source code, each of the four main steps is performed by executing four corresponding Perl executables on the command line. The data can be formatted as IM files [49], or FASTA files. While the data configuration file now accommodates multiple locus data, MTMLmsBayes can analyze single locus data sets thereby superceding the previous singlelocus msBayes. We distribute MTMLmsBayes as C source code, R scripts and Perl executables to be run on the command line after compiling on Linux, Mac OSX, and most POSIX systems using instructions from the README file. The MTMLmsBayes package is available from http://msbayes.sourceforge.net/ and also includes an online manual with installation/running instructions available from https://docs.google.com/Doc?docid=0AVkCIu87W8ooZGNyc3M2ZDhfNDJkZm5zd3dmcg&hl=en.
Results
For simulationbased testing, we generally compare estimates from PODS with the known hyperparameter values that simulated the PODS [10, 50] and calculate RMSE and RMSPE (root mean square error and root mean square posterior error) using these known values and each posterior mode estimate and the of 500 accepted posterior hyperparameter values in order to gauge the amount of uncertainty and bias associated with posterior estimates. PODS are simulated using random draws from the hyperprior of Ψ, where Ψ ranges from 1 to Y according to a discrete uniform distribution or alternatively are simulated under a history of simultaneous divergence (Ψ = 1; Ω = 0.0). For each set of conditions (i.e. number of loci, migration levels or chosen D_{ m }) we conduct ABC on sets of 100 independently generated PODS and for each we estimate E(τ) and Ω. For each set of 100 PODS and set of conditions we recycle the same 1,500,000 random draws from the prior (reference table), and use 500 accepted draws for ABC posterior estimation. In all cases, the simulated prior and sets of 100 PODS matches exactly with respect to sample size (i.e. number of loci and taxonpairs). Simulated data included of 5 to 20 taxonpairs and 1  64 loci whose length was 1100 basepairs.
Additionally we include an empirical analysis of three Australian avian taxonpairs that are hypothesized to have arisen simultaneously from three codistributed ancestral species due to the emergence of the Carpentarian barrier in northern Australia [51, 52]. Specifically, the three taxonpairs consist of the red backed fairy wren, Malurus melanocephalus melanocephalus and M. m. cruentatus (37 loci of 58  467 base pairs and mean of 27.8 individuals per descendent sister taxon), the blackthroated and longtailed finches, Poephila cincta and P. acuticauda (30 loci of 216  650 base pairs and one individual collected per descendent sister taxon) and the brown and blacktailed treecreepers, Climacteris picumnus and C. melanura (15 loci of 201  358 base pairs and mean of 9.5 individuals per descendent sister taxon).
Summary statistic vector D_{m}
When looking at pairs of PODS, the values of ΔD_{ m }between the pairs of simulated summary statistic vectors will be correlated with ΔΩ under optimal conditions of estimating Ω. Likewise when Ω is fixed at 0.0 (simultaneous divergence), values of ΔD_{ m }should be 0.0 under such optimal conditions for estimating Ω. To verify that ordering elements of y (taxonpair indicators) by the first moment of π_{ b } leads to more accurate estimates of Ω under simultaneous divergence than when ordering y (taxonpairs) by arbitrary order defined in the sampling configuration, we conduct simulation tests and plot frequency histograms of estimates of Ω given that PODS are generated under simultaneous divergence (Figure 2). Not only is the strategy for reordering D_{ m }superior to ordering arbitrarily, this advantage becomes more substantial as the number of taxonpairs increase (Figure 2D and 1F). Due to the exchangeability of π_{ b } across sample sizes, this sorting strategy minimizes ΔD_{ m }between observed and simulated data in cases when ΔΩ = 0.0 and Ω = 0.0 (simultaneous divergence). The increasing advantage of this reordering strategy as the number of sampled taxonpairs increases is expected given that ordering by the magnitude of π_{ b } obviates the need to simulate from the entire combinatorial sample space with regards to all possible orders from which the taxonpairs could be simulated when making the prior. Because this combinatorial sample space quickly increases with number of taxonpairs, ordering by some arbitrary rule such as number of samples per taxonpair results in increasing magnitudes ΔD_{ m }with greater number of taxonpairs given that Ω = 0.0 (Figure 2). Although using only the mean of π_{ b } across loci results in reasonable estimates of Ω, other available summary statistics are available for future applications of MTMLmsBayes that will involve testing complex multispecies histories other than simultaneous divergence. When this software pipeline is expanded to allow data consisting of large numbers of SNPs such that none of the information in the data are lost, we expect that a strategy involving genetic distances instead of Euclidian distances might work well for testing multitaxa hypotheses or alternatively using the first four moments of sets of summary statistics calculated across loci and/or taxa [46, 53]. For further information about the bourgeoning set of methods and strategies being developed for ABC, [12–14] provide thorough overviews.
Number of loci
As expected, increasing numbers of loci lead to more accurate estimates of Ω (Figures 3B, 4B and 5). However, improvement in estimation of E(τ) is not a monotonic increase with the number of loci (Figures 3A and 4A). The performance of estimating E(τ) with 4 loci is worse than a single locus, and the advantage of more loci is not reached until ≥16 loci are used (Figures 3A and 4A). In this case, sufficient sampling of loci is required to overcome the variance introduced by rate heterogeneity across loci. Estimating Ω on the other hand improves with 8 loci and continues to improve. We note that the accuracy of both estimators improves substantially with > 32 loci (Figures 3, 4, and 5).
Rate heterogeneity
Overall, estimating both Ω and E(τ) was relatively insensitive to whether or not the model of acrosslocus rate heterogeneity was correctly specified (Additional file 2). Generally, estimator performance was highest when the PODS had equal rates, but we note that PODS with unequal rates resulted in high accuracy in estimates of Ω and E(τ) irregardless of whether rate heterogeneity or rate uniformity was correctly specified in the prior model.
Postdivergence Migration
Given data sets with 16 loci per each of five taxonpairs, migration had a negative impact on the estimation of Ω and E(τ) but the magnitude of this negative impact depended on whether one used the correct migration model for simulating the prior. As theory predicts [54–56], we generally found that estimates of Ω and E(τ) became less reliable with increasing migration (Figures 6 and 7) even when the correct migration models were used. Migration model misspecification also influenced estimator bias and precision. When the PODS are generated under isolation, the estimators of Ω and E(τ) generally became less accurate with increasing migration model misspecification. Likewise, when PODS were generated under a migration model, model misspecification resulted in higher estimator bias and less precision as quantified by RMSE and RMSPE.
Overall, this simulation analysis demonstrates that quantifying the level of temporal congruence in multitaxa divergence will be augmented if one first tests for migration so that an appropriate hyperprior model can be specified. Therefore it would be wise to test whether a migration or isolation model is justified in the taxonpairs using informational theoretic and MCMC techniques [56, 57] or ABC model choice before quantifying the level of congruence in multitaxa divergence. Overall, this simulation analysis demonstrates that our multilocus test for simultaneous divergence will work better if one is interested in testing simultaneous cessation of all gene flow rather then testing for simultaneous divergence with some postisolation gene flow. However, it remains to be seen whether larger numbers of loci and/or other summary statistics can better infer multitaxa divergence with limited migration or secondary contact.
Hierarchical ABC model choice among competing migration models
Because the accuracy of estimation can depend on assumptions about migration after divergence, one can first use ABC model choice techniques [32, 58] to compare the posterior probability of isolation and postdivergence migration models in the context of our hierarchical multitaxa divergence model. Specifically we treat models of isolation and migration as a set of models specified by a categorical model indicator parameter that can be estimated via ABC. In this case the acceptance rejection step can be followed by a polychotomous regression step that has been shown to improve estimation of discrete categorical parameters [15, 22, 25, 32]. To test the accuracy of this technique, the five taxonpair data was simulated using 3,000,000 random draws from the hyperprior with the three different migration models simulating the data with equal probability (one isolation model and two migration models). For the two migration models, each of the five taxonpair's migration parameter (Nm; number of effective migrants per generation) is independently drawn from a uniform prior distribution (0.0,1.0) or (0.0,10.0) depending on which of the two migration models. Subsequently the posterior for the model indicator parameter conferring to isolation or the two different migration levels (Nm upper bounds of 1.0 and 10.0) was sampled from the 500 closest accepted matches obtained with the ABC algorithm with and without subsequent polychotomous regression.
The accuracy of this ABC model choice procedure was then assessed by conducting this procedure on 100 PODS of five taxonpairs and 16 loci simulated under each of the three different migration models (isolation and Nm upper bounds of 1.0 and 10.0). The probability of choosing the correct migration model ranged from 0.77 to 0.96 when one used a summary statistic vector D_{ m } that included the across loci mean π_{b} and Ψ_{ W } (Additional file 3) whereas using π_{b} only resulted in fewer successful model choices (probability of choosing correct model ranging from 0.570.72). Indeed, Wakeley's Ψ_{ W } was developed as a way to distinguish between migration and isolation models [43] and when harnessed within our hierarchical ABC framework, we demonstrate it to have potential application given a multiple taxonpair model. Additionally, the use of polychotomous regression greatly improved the probability of successful model choice over using direct nontransformed accepted values. Likewise, the Bayes factor support for the correct model was augmented when using π_{b} and Ψ_{ W } as well as polychotomous regression (Additional file 3).
Empirical analysis
To demonstrate how MTMLmsBayes can test for simultaneous divergence given large numbers of loci and postdivergence migration, we used 1537 loci collected from three bird taxonpairs all of which consist of sister taxonpairs that presently span the Carpentarian barrier in northern Australia [51, 52]. This includes the brown and blacktailed treecreepers (Climacteris picumnus and C. melanura), the blackthroated and longtailed finches (Poephila cincta and P. acuticauda) and the eastern and western ssp. of redbacked wrens (Malurus melanocephalus melanocephalus and M. m. cruentatus). Results strongly suggest that all three sister taxonpairs diverged at the same time and hence could have formed by way of the same geoclimatic event that formed the Carpentarian barrier in northern Australia. Furthermore, this result of simultaneous divergence was insensitive to whether or not one incorporated low levels of migration after divergence. The time of this simultaneous divergence was 81,000 y.b.p. and 200,000 y.b.p. under isolation and low migration models respectively.
As expected from theory and shown in our simulation analysis (Figures 6 and 7), tests of simultaneous divergence and divergence time estimation are dependent on model assumptions about postdivergence migration [54–56] and therefore we initially used ABC model choice [32, 58] to compare the posterior probability of two models; complete isolation and postdivergence migration across all three taxonpairs (Figure 1). To generate simulated data for ABC, the three taxonpair data set was simulated 3,000,000 times using random draws from the hyperprior and both isolation and migration models were used to simulate the data with equal probability. Under the migration model, the three values of Nm for each of the three taxonpairs (number of effective migrants per generation) are independently drawn from a uniform prior distribution (0.0,1.0) and assigned to each taxonpair. After conducting the ABC model selection, the low migration model had more support (Pr(migration  D_{ m } ) = 0.65)), yet not enough for "high" or "moderate" Bayes factor support [59]. Alternatively, we estimated E(τ) and Ω from mixed isolation/migration priors such that estimates of E(τ) and Ω are averaged across the relative posterior probability of the isolation model and migration model. In this case, the goal is not to test the models but to obtained estimates of E(τ) and Ω while allowing uncertainty in model selection.
As always, translating scaled divergence time estimates into real time estimates depends on assumptions about DNA mutation rates and here we report real time estimates based on DNA mutation rates reported previously. An assumed mean rate across loci of 5.0 × 10^{9} per site per generation (as reported by [52]) and a two year generation time results in mean divergence time, E(τ) estimates of 81,000 y.b.p. and 200,000 y.b.p. under isolation and low migration models respectively. These estimates are generally consistent with the reported divergence time estimates of the wrens (Lee and Edwards 2008; 270,000 y.b.p) and finches (Jennings and Edwards 2005; 432,000 y.b.p.) using the same rates and a similar coalescentbased isolation with migration model that used Markov chain Monte Carlo [49, 60]. We additionally note that Lee and Edwards [61] estimated low levels of migration (< 1.0 migrants per generation) in the fairy wrens which is also consistent with the higher posterior probability for the low migration model that we found via ABC model selection. Further, the older and less precise estimate of means divergence time under migration than under isolation is expected due to migration breaking up the correlation between coalescent times and divergence times [54–56].
Conclusions
Multispecies comparative phylogeographic inference using genetic data from large numbers of nonmodel taxa will increasingly become a standard tool for understanding the interplay between geography, climate change, speciation, extinction, demographic changes, and species interactions as well as making links between different types of biodiversity, ecological services and ultimately wellinformed conservation policy [62, 63]. Inferring how whole assemblages of species react to putative geographical barriers is central to obtaining these larger goals and MTMLmsBayes will become an important bioinformatics tool for such inference given multilevel models with large amounts of complexity. Phylogeographic data sets with multiple codistributed taxonpairs with genetic data collected from multiple loci are rapidly emerging [64–67], and here we demonstrate that correct inference of simultaneous divergence is somewhat robust against violations in assumptions about among locus rate heterogeneity although incomplete isolation with postdivergence migration can make inference of simultaneous divergence difficult. Furthermore, it is likely that other demographic complexities such as preisolation subdivision, diminishing/accelerating levels of postisolation migration, and recombination are likely to affect inference [68]. Although MTMLmsBayes does optionally allow for intragenic recombination, testing how ignoring this parameter biases inference is beyond the scope of this work and researchers should test for recombination or use nonrecombinant blocks for analysis.
The modular design of MTMLmsBayes further allows simulationbased model validation and posterior predictive model fitting and will be able to interface with other bioinformatics tools developed for ABC [20, 21]. Moreover, the modular design will ultimately allow implementing various constrained analyses for testing an array of multitaxon histories beyond the tests of migration and simultaneous divergence presented here so that researchers will finally be able to make large scale biogeographic inference across whole communities with sufficient demographic complexity.
Availability and requirements
We distribute MTMLmsBayes as C source code, R scripts and Perl executables under opensource, GNU General Public License to be run on the command line after compiling on Linux, Mac OSX, and most POSIX systems using instructions from the README file. The MTMLmsBayes package is available from sourceforge http://msbayes.sourceforge.net/ and also includes an online manual with installation/running instructions available from as well as associated R scripts to conduct simulation testing are available from http://qcpages.qc.cuny.edu/Biology/Hickerlab/Software/Software.html
List of abbreviations
 The abbreviations include ABC:

(Approximate Bayesian Computation)
 HABC:

(Hierarchical Approximate Bayesian computation)
 RMSE:

(Root Mean Square Error)
 RMSPE:

(Root Mean Square Posterior Error)
 and PODS:

(PseudoObserved Data Sets).
Declarations
Acknowledgements
We thank J. Lee, B. Jennings, and S. Edwards for permission to analyze the empirical multilocus data from three avian taxonpairs. We thank the staff of the City University of New York HPCF (High Performance Computing Facility) for computational resources. We thank B. Carstens and the staff of the Southwestern research station for hosting the Statistical Phylogeography course and for the students of this course for making practical suggestions. We thank M. Beaumont for kindly providing R scripts and D. Hudson and E. Stahl for permission to use E. Stahl's finite sites version of his ms coalescent simulator under GNU Public License. Support for N. Takebayashi was supported by National Science Foundation (DEB0640520) and Alaska INBRE Grant Number RR016466 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH). Support for M. Hickerson was provided by National Science Foundation (DEB0743648).
Authors’ Affiliations
References
 Bermingham E, Moritz C: Comparative phylogeography: concepts and applications. Mol Ecol 1998, 7: 367–369. 10.1046/j.1365294x.1998.00424.xView Article
 Arbogast BS, Kenagy GJ: Comparative phylogeography as an integrative approach to historical biogeography. J Biogeogr 2001, 28: 819–825. 10.1046/j.13652699.2001.00594.xView Article
 Coyne JA, Orr HA: Speciation. Sunderland, MA: Sinauer Associates Inc; 2004.
 Avise JC: Phylogeography: The history and formation of species. Cambridge: Harvard University Press; 2000.
 Hubbell SP: The Unified Neutral Theory of Biodiversity and Biogeography. Princeton, NJ: Princeton University Press; 2001.
 Vera C, Wheat CW, Fescemyer HW, Frilander MJ, Crawford DL, Hanski I, Marden JH: Rapid transcriptome characterization for a nonmodel organism using massively parallel 454 pyrosequencing. Mol Ecol 2008., 2371:
 Graham CH, Parra JL, Rahbek C, McGuire JA: Phylogenetic structure in tropical hummingbird communities. Proceedings of the National Academy of Sciences 2009, 106(Supplement 2):19673–19678. 10.1073/pnas.0901649106View Article
 Nielsen R, Beaumont MA: Statistical inferences in phylogeography. Mol Ecol 2009, 18: 1034–1047. 10.1111/j.1365294X.2008.04059.xView ArticlePubMed
 Tallmon DA, Luikart G, Beaumont BA: Comparative evaluation of a new effective population size estimator based on approximate Bayesian computation. Genetics 2004, 167: 977–988. 10.1534/genetics.103.026146PubMed CentralView ArticlePubMed
 Excoffier L, Estoup A, Cornuet JM: Bayesian analysis of an admixture model with mutations and arbitrarily linked markers. Genetics 2005, 169: 1727–1738. 10.1534/genetics.104.036236PubMed CentralView ArticlePubMed
 Hickerson MJ, Stahl E, Lessios HA: Test for simultaneous divergence using approximate Bayesian computation. Evolution 2006, 60: 2435–2453.View ArticlePubMed
 Beaumont MA: Approximate Bayesian Computation in Evolution and Ecology. Annual Review of Ecology, Evolution, and Systematics 2010, 41(1):379–406. 10.1146/annurevecolsys102209144621View Article
 Csilléry K, Blum MGB, Gaggiotti OE, Francois O: Approximate Bayesian Computation (ABC) in practice. Trends Ecol Evol 2010, 25(7):410–418. 10.1016/j.tree.2010.04.001View ArticlePubMed
 Bertorelle G, Benazzo A, S M: ABC as a flexible framework to estimate demography over space and time: some cons, many pros. Mol Ecol 2010, 19(13):2609–2625. 10.1111/j.1365294X.2010.04690.xView ArticlePubMed
 Carnaval A, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C: Stability predicts genetic diversity in the Brazilian Atlantic Forest Hotspot. Science 2009, 323: 785–789. 10.1126/science.1166955View ArticlePubMed
 Lopes JS, Balding D, Beaumont MA: PopABC: a program to infer historical demographic parameters. Bioinformatics 2009, 25(20):2747–2749. 10.1093/bioinformatics/btp487View ArticlePubMed
 Anderson CNK, Ramakrishnan U, Chan YL, Hadly EA: Serial SimCoal: A population genetic model for data from multiple populations and points in time. Bioinformatics 2005, 21: 1733–1734. 10.1093/bioinformatics/bti154View ArticlePubMed
 Cornuet JM, Santos F, Beaumont MA, Robert CP, Marin JM, Balding DJ, Guillemaud T, Estoup A: Inferring population history with DIY ABC: a userfriendly approach to Approximate Bayesian Computation. Bioinformatics 2008, 24(23):2713–2719. 10.1093/bioinformatics/btn514PubMed CentralView ArticlePubMed
 Jobin MJ, Mountain JL: REJECTOR: software for population history inference from genetic data via a rejection algorithm. Bioinformatics 2008, 24: 2936–2937. 10.1093/bioinformatics/btn540PubMed CentralView ArticlePubMed
 Thornton K: Automating approximate Bayesian computation by local linear regression. BMC Genet 2009, 10(1):35. 10.1186/147121561035PubMed CentralView ArticlePubMed
 Wegmann D, Leuenberger C, Neuenschwander S, Excoffier L: ABCtoolbox: a versatile toolkit for approximate Bayesian computations. BMC Bioinformatics 2010, 11(1):116. 10.1186/1471210511116PubMed CentralView ArticlePubMed
 Palero F, Lopes J, Abello P, Macpherson E, Pascual M, Beaumont M: Rapid radiation in spiny lobsters (Palinurus spp) as revealed by classic and ABC methods using mtDNA and microsatellite data. BMC Evol Biol 2009, 9(1):263. 10.1186/147121489263PubMed CentralView ArticlePubMed
 Storz JF, Beaumont BA: Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsattelite DNA variation using a hierarchical Bayesian model. Evolution 2002, 56: 154–166.View ArticlePubMed
 Hickerson MJ, Stahl E, Takebayashi N: msBayes: Pipeline for testing comparative phylogeographic histories using hierarchical approximate Bayesian computation. BMC Bioinformatics 2007, 8: 268. 10.1186/147121058268PubMed CentralView ArticlePubMed
 Fagundes NJR, Ray N, Beaumont M, Neuenschwander S, Salzano FM, Bonatto SL, Excoffier L: Statistical evaluation of alternative models of human evolution. Proceedings of the National Academy of Sciences of the United States of America 2007, 104: 17614–17619. 10.1073/pnas.0708280104PubMed CentralView ArticlePubMed
 Estoup A, Beaumont BA, Sennedot F, Moritz C, Cornuet JM: Genetic analysis of complex demographic scenarios: spatially expanding populations of the cane toad, Bufo marinus . Evolution 2004, 58: 2021–2036.View ArticlePubMed
 Tavaré S, Balding DJ, Griffiths RC, Donnelly P: Inferring coalescence times from DNA sequence data. Genetics 1997, 145: 505–518.PubMed CentralPubMed
 Weiss G, von Haeseler A: Inference of population history using a likelihood approach. Genetics 1998, 149: 1539–1546.PubMed CentralPubMed
 Pritchard JK, Seielstad MT, PL A, Feldman MW: Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Mol Biol Evol 1999, 16: 1791–1798.View ArticlePubMed
 Marjoram PM, Molitor J, Plagnol V, Tavaré S: Markov chain Monte Carlo without likelihoods. Proc Natl Acad Sci USA 2003, 100: 15324–15328. 10.1073/pnas.0306899100PubMed CentralView ArticlePubMed
 Beaumont MA, Zhang W, Balding DJ: Approximate Bayesian computation in population genetics. Genetics 2002, 162: 2025–2035.PubMed CentralPubMed
 Beaumont MA: Joint determination of topology, divergence time and immigration in population trees. In Simulations, Genetics and Human Prehistory. Edited by: Matsumura S, Forster P, Renfrew C. Cambridge: McDonald Institute for Archaeological Research; 2008:135–154.
 François O, Blum MGB, Jakobsson M, Rosenberg NA: Demographic history of european populations of arabidopsis thaliana . PLoS Genet 2008, 4: e1000075.PubMed CentralView ArticlePubMed
 Leuenberger C, Wegmann D: Bayesian Computation and Model Selection Without Likelihoods. Genetics 184(1):243–252. 10.1534/genetics.109.109058
 Blum MGB, François O: Nonlinear regression models for Approximate Bayesian Computation. Statistics and Computing 2010, 20(1):63–73. 10.1007/s1122200991160View Article
 Hudson RR: Generating samples under a WrightFisher neutral model of genetic variation. Bioinformatics 2002, 18: 337–338. 10.1093/bioinformatics/18.2.337View ArticlePubMed
 Jukes TH, Cantor CH: Evolution of protein molecules. In Mammalian protein metabolism. Edited by: Munro HM. New York: Academic Press; 1969:21–123.View Article
 Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Evol 1981, 17: 368–376. 10.1007/BF01734359View ArticlePubMed
 Hasegawa M, Kishino H, Yano TA: Dating of the humanape splitting by a molecular clock of mitochondrial DNA. J Mol Evol 1985, 22: 160–174. 10.1007/BF02101694View ArticlePubMed
 Nei N, Kumar S: Molecular Evolution and Phylogenetics. Oxford: Oxford University Press; 2000.
 Watterson GA: On the number of segregating sites in genetic models without recombination. Theor Popul Biol 1975, 7: 256–276. 10.1016/00405809(75)900209View ArticlePubMed
 Sherwin WB, Jabot F, Rush R, Rossetto M: Measurement of biological information with applications from genes to landscapes. Mol Ecol 2006, 15: 2857–2869. 10.1111/j.1365294X.2006.02992.xView ArticlePubMed
 Wakeley J: Distinguishing migration from isolation using the variance of pairwise differences. Theor Popul Biol 1996, 49: 369–386. 10.1006/tpbi.1996.0018View ArticlePubMed
 Nei M, Li W: Mathematical model for studying variation in terms of restriction endonucleases. Proceedings of the National Academy of Sciences of the United States of America 1979, 76: 5269–5273. 10.1073/pnas.76.10.5269PubMed CentralView ArticlePubMed
 Wakeley J: The variance of pairwise nucleotide differences in two populations with migration. Theor Popul Biol 1996, 49: 39–57. 10.1006/tpbi.1996.0002View ArticlePubMed
 Bazin E, Dawson KJ, Beaumont MA: Likelihoodfree Inference of Population Structure and Local Adaptation in a Bayesian Hierarchical Model. Genetics 2010. genetics. 109.112391 genetics. 109.112391
 Papoulis A: Probability, Random Variables, and Stochastic Processes. 2nd edition. New York: McGrawHill; 1984.
 Takahata N, Nei M: Gene genealogy and variance of interpopulational nucleotide differences. Genetics 1985, 110: 325–344.PubMed CentralPubMed
 Hey J, Nielsen R: Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis . Genetics 2004, 167: 747–760. 10.1534/genetics.103.024182PubMed CentralView ArticlePubMed
 Cook SR, Gelman A, Rubin DB: Validation of Software for Bayesian Models Using Posterior Quantiles. Journal of Computational and Graphical Statistics 2006, 15(3):675–692. 10.1198/106186006X136976View Article
 Jennings WB, Edwards SV: Speciational history of australian grass finches ( poephila ) inferred from thirty gene trees. Evolution 2005, 59: 2033–2047.PubMed
 Lee JY, Edwards SV: Divergence Across Australia's Carpentarian Barrier: Statistical Phylogeography of the RedBacked Fairy Wren ( Malurus melanocephalus ). Evolution 2008, 62(12):3117–3134. 10.1111/j.15585646.2008.00543.xView ArticlePubMed
 Sousa VC, Fritz M, Beaumont MA, Chikhi L: Approximate Bayesian Computation Without Summary Statistics: The Case of Admixture. Genetics 2009, 181(4):1507–1519. 10.1534/genetics.108.098129PubMed CentralView ArticlePubMed
 Rosenberg NA, Feldman MW: The relationship between coalescence times and population divergence times. In Modern Developments in Theoretical Population Genetics. Edited by: Slatkin M, Veuille M. Oxford: University Press; 2002:130–164.
 Arbogast BS, Edwards SV, Wakeley J, Beerli P, Slowinski JB: Estimating divergence times from molecular data on phylogenetic and population genetic timescales. Annu Rev Ecol Syst 2002, 33: 707–740. 10.1146/annurev.ecolsys.33.010802.150500View Article
 Nielsen R, Wakeley J: Distinguishing migration from isolation: A Markov chain Monte Carlo approach. Genetics 2001, 158(2):885–896.PubMed CentralPubMed
 Carstens BC, Stoute HN, Reid NM: An informationtheoretical approach to phylogeography. Mol Ecol 2009, 18(20):4270–4282. 10.1111/j.1365294X.2009.04327.xView ArticlePubMed
 Verdu P, Austerlitz F, Estoup A, Vitalis R, Georges M, ThÈry S, Froment A, Le Bomin S, Gessain A, Hombert J: Origins and Genetic Diversity of Pygmy HunterGatherers from Western Central Africa. Curr Biol 2009, 19(4):312–318. 10.1016/j.cub.2008.12.049View ArticlePubMed
 Kass RE, Raftery A: Bayes factors. Journal of the American Statistical Association 1995, 90: 773–795. 10.2307/2291091View Article
 Rannala B, Yang ZH: Bayes Estimation of Species Divergence Times and Ancestral Population Sizes Using DNA Sequences From Multiple Loci. Genetics 2003, 164: 1645–1656.PubMed CentralPubMed
 Lee JY, Edwards SV: Divergence Across Australia's Carpentarian Barrier: Statistical Phylogeography of the RedBacked Fairy Wren (Malurus melanocephalus). Evolution 2009, 62(12):3117–3134. 10.1111/j.15585646.2008.00543.xView Article
 Knowles LL: Statistical Phylogeography. Annual Review of Ecology, Evolution, and Systematics 2009, 40(1):593–612. 10.1146/annurev.ecolsys.38.091206.095702View Article
 Hickerson MJ, Carstens BC, CavenderBares J, Crandall KA, Graham CH, Johnson JB, Rissler L, Victoriano PF, Yoder AD: Phylogeography's past, present, and future: 10 years after. Mol Phylogen Evol 2010, 54(1):291–301. 10.1016/j.ympev.2009.09.016View Article
 Moyle LC: Ecological and evolutionary genomics in the wild tomatoes (solanum sect. Lycopersicon). Evolution 2008, 62(12):2995–3013. 10.1111/j.15585646.2008.00487.xView ArticlePubMed
 Dolman G, Moritz C: A multilocus perspective on refugial isolation and divergence in rainforest skinks ( carlia ). Evolution 2006, 60(3):573–582.View ArticlePubMed
 Hurt C, Anker A, Knowlton N: A multilocus test of simultaneous divergence across the isthmus of panama using snapping shrimp in the genus alpheus. Evolution 2009, 63(2):514–530. 10.1111/j.15585646.2008.00566.xView ArticlePubMed
 Rogers S, Bernatchez L: The Genetic Architecture of Ecological Speciation and the Association with Signatures of Selection in Natural Lake Whitefish (Coregonus sp. Salmonidae) Species Pairs. Mol Biol Evol 2007, 24(6):1423–1438. 10.1093/molbev/msm066View ArticlePubMed
 Becquet Cl, Przeworski M: Learning about Modes of Speciation by Computational Approaches. Evolution 2009, 63(10):2547–2562. 10.1111/j.15585646.2009.00662.xView ArticlePubMed
Copyright
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/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.