MTMLmsBayes: approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity.  
Jump to Full Text  
MedLine Citation:

PMID: 21199577 Owner: NLM Status: MEDLINE 
Abstract/OtherAbstract:

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). 
Authors:

Wen Huang; Naoki Takebayashi; Yan Qi; Michael J Hickerson 
Related Documents
:

21095857  Toward unsupervised adaptation of lda for braincomputer interfaces. 20831547  An analysis and comparison of commonly available united kingdom prescribing resources. 21354417  The hormonal and circadian basis for insect photoperiodic timing. 21459847  The intfold server: an integrated web resource for protein fold recognition, 3d model q... 17758897  Greenhouse warming: dueling models: future u.s. climate uncertain. 8805757  Random effects models in latent class analysis for evaluating accuracy of diagnostic te... 
Publication Detail:

Type: Journal Article; Research Support, N.I.H., Extramural; Research Support, U.S. Gov't, NonP.H.S. Date: 20110103 
Journal Detail:

Title: BMC bioinformatics Volume: 12 ISSN: 14712105 ISO Abbreviation: BMC Bioinformatics Publication Date: 2011 
Date Detail:

Created Date: 20110201 Completed Date: 20110711 Revised Date: 20130703 
Medline Journal Info:

Nlm Unique ID: 100965194 Medline TA: BMC Bioinformatics Country: England 
Other Details:

Languages: eng Pagination: 1 Citation Subset: IM 
Affiliation:

Biology Department, City University of New York, Queens College, 6530 Kissena Blvd, Flushing, NY 113671597, USA. wenhuang19@yahoo.com 
Export Citation:

APA/MLA Format Download EndNote Download BibTex 
MeSH Terms  
Descriptor/Qualifier:

Algorithms Animals Bayes Theorem Birds / genetics Demography Evolution, Molecular Gene Flow* Genetic Loci Genetic Variation Models, Statistical* Phylogeography Sequence Analysis, DNA Software* 
Grant Support  
ID/Acronym/Agency:

RR016466/RR/NCRR NIH HHS 
Comments/Corrections 
Full Text  
Journal Information Journal ID (nlmta): BMC Bioinformatics ISSN: 14712105 Publisher: BioMed Central 
Article Information Download PDF Copyright ©2011 Huang et al; licensee BioMed Central Ltd. openaccess: Received Day: 5 Month: 5 Year: 2010 Accepted Day: 3 Month: 1 Year: 2011 collection publication date: Year: 2011 Electronic publication date: Day: 3 Month: 1 Year: 2011 Volume: 12First Page: 1 Last Page: 1 ID: 3031198 Publisher Id: 14712105121 PubMed Id: 21199577 DOI: 10.1186/14712105121 
MTMLmsBayes: Approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity  
Wen Huang1  Email: wenhuang19@yahoo.com 
Naoki Takebayashi2  Email: ntakebayashi@alaska.edu 
Yan Qi3  Email: YQi@gc.cuny.edu 
Michael J Hickerson13  Email: michael.hickerson@qc.cuny.edu 
1Biology Department, City University of New York, Queens College, 6530 Kissena Blvd, Flushing, NY 113671597, USA 

2Institute of Arctic Biology and Department of Biology and Wildlife, 311 Irving I Bldg, University of Alaska, Fairbanks, AK 99775, USA 

3The Graduate Center of the City University of New York, 365 5th Ave, New York, NY 10016, USA 
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.
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}].
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 ith 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}].
MTMLmsBayes generates finite sites DNA sequence data simulated under a coalescent demographic model to perform ABC. The data generation step is accomplished by msDQH which is a version of Hudson's classic coalescent simulator which simulates finite sites DNA sequence data under specified demographic model [^{36}]. The general multiple taxonpair hierarchical ABC model of divergence with migration and size change that can be implemented in MTMLmsBayes is presented elsewhere [^{11},^{24}] and generally involves using the multiple taxonpair data to estimate hyperparameters and parameter summaries that quantify the variability in divergence times across Y taxonpairs (Additional File 1; Figure 1). This includes Ψ, the number of different divergence times across Y taxonpairs, which follows a discrete uniform prior distribution ranging from 1 to Y. Additionally one can estimate the mean divergence time, E(τ), where τ is the time of divergence between a pair of descendent population pairs (in coalescent time units of 4N generations, where N is the sum of current effective population sizes of the two descendent sister populations), as well as estimate Ω, the dispersion index of τ, (Var(τ)/E(τ)). If one conducts the analysis using a partially constrained model where the number of divergence times (Ψ) is held to a single value across the Y taxonpairs, one can subsequently estimate each of the Ψ divergence times (τ_{1}, ..., τ_{Ψ}), as well as the number of taxa that split at each of the Ψ times (Ψ_{1}, ...,Ψ_{Ψ}). Following the algorithm detailed in [^{11},^{24}], the Ψ divergence times τ_{1}, ..., τ_{Ψ }are randomly drawn from a userspecified uniform prior distribution and these Ψ divergence times are subsequently assigned randomly to Ψ taxonpairs of the Y taxonpairs with the remaining Y  Ψ taxonpairs randomly picking divergence times from τ_{1}, ..., τ_{Ψ }with replacement.
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.
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.
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.
To ascertain how well MTMLmsBayes quantifies the congruence of divergence times under a number of different conditions, we conduct an extensive simulation analysis by generating PODS (pseudoobserved data sets; [^{14}]) and quantifying the accuracy and precision of estimates on the known parameter values used to generate the PODS. Specifically, we assessed: 1.) the advantage of reordering elements of y (taxonpair indicators) within D_{m }by descending magnitude of π_{b }averaged across loci with respect to estimating Ω as a function of number of taxonpairs (Y) within the sample (Figure 2); 2.) the effect of increasing numbers of loci (1, 4, 8, 16, 32, and 64 loci) when estimating E(τ) and Ω (Figures 3, 4, and 5); 3.) the consequences of allowing for and ignoring rate heterogeneity across loci (Additional file 2); and 4.) how different levels of postdivergence migration influence estimates of E(τ) and Ω and how this is influenced by migration/isolation model misspecification (Figures 6 and 7).
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).
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.
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).
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.
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.
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).
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.
Hyperparameter estimates of Ψ and Ω indicate an inference of simultaneous divergence, with Ψ = 1 having the highest probability irregardless of which migration/isolation model is used, (Pr(Ψ = 1 D_{m}, isolation) = 0.67; Pr(Ψ = 1 D_{m}, migration) = 0.58; and Pr(Ψ = 1D_{m}, mixed model) = 0.60). Likewise, Ω (the dispersion index characterizing the variability in divergence times) indicated synchronous divergence irregardless of migration model with mode estimates of Ω = Var(τ)/E(τ)) = 0.0 across all three migration/isolation models (Figure 8). The resulting Bayes factor comparing models of simultaneous divergence (Ψ = 1) and nonsimultaneous divergence (Ψ>1) did depend on whether migration was assumed with moderate support for simultaneous divergence given isolation (B(Ψ = 1,Ψ>1) = 4.05), weak support for simultaneous divergence given migration (B(Ψ = 1,Ψ>1) = 2.76) and the mixed model (B(Ψ = 1,Ψ>1) = 2.92). Consistent with our expectations that allowing migration will result in divergence time estimates with more uncertainty, the posterior estimates of mean divergence time and tests of simultaneous divergence are less precise under the low migration model than under a pure isolation model, and the posterior estimates of mean divergence time, (E(τ)), are older under migration than under isolation (Figure 8).
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}].
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.
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
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).
WH, NT, YQ, and MJH developed C, R, and Perl routines for the multitaxa/multiloci model with rate heterogeneity for ABC estimation and model choice. WH and MJH did the extensive simulation testing and MJH conducted the empirical analysis. MJH and NT maintain MTMLmsBayes and NT developed the installation configurations. All authors read and approved the final version of the manuscript.
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).
References
Bermingham E,Moritz C,Comparative phylogeography: concepts and applicationsMol EcolYear: 1998736736910.1046/j.1365294x.1998.00424.x  
Arbogast BS,Kenagy GJ,Comparative phylogeography as an integrative approach to historical biogeographyJ BiogeogrYear: 20012881982510.1046/j.13652699.2001.00594.x  
Coyne JA,Orr HA,SpeciationYear: 2004Sunderland, MA: Sinauer Associates Inc  
Avise JC,Phylogeography: The history and formation of speciesYear: 2000Cambridge: Harvard University Press  
Hubbell SP,The Unified Neutral Theory of Biodiversity and BiogeographyYear: 2001Princeton, NJ: Princeton University Press  
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 pyrosequencingMol EcolYear: 20082371  
Graham CH,Parra JL,Rahbek C,McGuire JA,Phylogenetic structure in tropical hummingbird communitiesProceedings of the National Academy of SciencesYear: 2009106Supplement 2196731967810.1073/pnas.0901649106  
Nielsen R,Beaumont MA,Statistical inferences in phylogeographyMol EcolYear: 2009181034104710.1111/j.1365294X.2008.04059.x19207258  
Tallmon DA,Luikart G,Beaumont BA,Comparative evaluation of a new effective population size estimator based on approximate Bayesian computationGeneticsYear: 200416797798810.1534/genetics.103.02614615238546  
Excoffier L,Estoup A,Cornuet JM,Bayesian analysis of an admixture model with mutations and arbitrarily linked markersGeneticsYear: 20051691727173810.1534/genetics.104.03623615654099  
Hickerson MJ,Stahl E,Lessios HA,Test for simultaneous divergence using approximate Bayesian computationEvolutionYear: 2006602435245317263107  
Beaumont MA,Approximate Bayesian Computation in Evolution and EcologyAnnual Review of Ecology, Evolution, and SystematicsYear: 201041137940610.1146/annurevecolsys102209144621  
Csilléry K,Blum MGB,Gaggiotti OE,Francois O,Approximate Bayesian Computation (ABC) in practiceTrends Ecol EvolYear: 201025741041810.1016/j.tree.2010.04.00120488578  
Bertorelle G,Benazzo A,S M,ABC as a flexible framework to estimate demography over space and time: some cons, many prosMol EcolYear: 201019132609262510.1111/j.1365294X.2010.04690.x20561199  
Carnaval A,Hickerson MJ,Haddad CFB,Rodrigues MT,Moritz C,Stability predicts genetic diversity in the Brazilian Atlantic Forest HotspotScienceYear: 200932378578910.1126/science.116695519197066  
Lopes JS,Balding D,Beaumont MA,PopABC: a program to infer historical demographic parametersBioinformaticsYear: 200925202747274910.1093/bioinformatics/btp48719679678  
Anderson CNK,Ramakrishnan U,Chan YL,Hadly EA,Serial SimCoal: A population genetic model for data from multiple populations and points in timeBioinformaticsYear: 2005211733173410.1093/bioinformatics/bti15415564305  
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 ComputationBioinformaticsYear: 200824232713271910.1093/bioinformatics/btn51418842597  
Jobin MJ,Mountain JL,REJECTOR: software for population history inference from genetic data via a rejection algorithmBioinformaticsYear: 2008242936293710.1093/bioinformatics/btn54018936052  
Thornton K,Automating approximate Bayesian computation by local linear regressionBMC GenetYear: 20091013510.1186/14712156103519583871  
Wegmann D,Leuenberger C,Neuenschwander S,Excoffier L,ABCtoolbox: a versatile toolkit for approximate Bayesian computationsBMC BioinformaticsYear: 201011111610.1186/147121051111620202215  
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 dataBMC Evol BiolYear: 20099126310.1186/14712148926319900277  
Storz JF,Beaumont BA,Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsattelite DNA variation using a hierarchical Bayesian modelEvolutionYear: 20025615416611913661  
Hickerson MJ,Stahl E,Takebayashi N,msBayes: Pipeline for testing comparative phylogeographic histories using hierarchical approximate Bayesian computationBMC BioinformaticsYear: 2007826810.1186/14712105826817655753  
Fagundes NJR,Ray N,Beaumont M,Neuenschwander S,Salzano FM,Bonatto SL,Excoffier L,Statistical evaluation of alternative models of human evolutionProceedings of the National Academy of Sciences of the United States of AmericaYear: 2007104176141761910.1073/pnas.070828010417978179  
Estoup A,Beaumont BA,Sennedot F,Moritz C,Cornuet JM,Genetic analysis of complex demographic scenarios: spatially expanding populations of the cane toad,Bufo marinusEvolutionYear: 2004582021203615521459  
Tavaré S,Balding DJ,Griffiths RC,Donnelly P,Inferring coalescence times from DNA sequence dataGeneticsYear: 19971455055189071603  
Weiss G,von Haeseler A,Inference of population history using a likelihood approachGeneticsYear: 1998149153915469649540  
Pritchard JK,Seielstad MT,PL A,Feldman MW,Population growth of human Y chromosomes: a study of Y chromosome microsatellitesMol Biol EvolYear: 1999161791179810605120  
Marjoram PM,Molitor J,Plagnol V,Tavaré S,Markov chain Monte Carlo without likelihoodsProc Natl Acad Sci USAYear: 2003100153241532810.1073/pnas.030689910014663152  
Beaumont MA,Zhang W,Balding DJ,Approximate Bayesian computation in population geneticsGeneticsYear: 20021622025203512524368  
Beaumont MA,Matsumura S, Forster P, Renfrew CJoint determination of topology, divergence time and immigration in population treesSimulations, Genetics and Human PrehistoryYear: 2008Cambridge: McDonald Institute for Archaeological Research135154  
François O,Blum MGB,Jakobsson M,Rosenberg NA,Demographic history of european populations of arabidopsis thalianaPLoS GenetYear: 20084e100007518483550  
Leuenberger C,Wegmann D,Bayesian Computation and Model Selection Without LikelihoodsGenetics184124325210.1534/genetics.109.10905819786619  
Blum MGB,François O,Nonlinear regression models for Approximate Bayesian ComputationStatistics and ComputingYear: 2010201637310.1007/s1122200991160  
Hudson RR,Generating samples under a WrightFisher neutral model of genetic variationBioinformaticsYear: 20021833733810.1093/bioinformatics/18.2.33711847089  
Jukes TH,Cantor CH,Munro HMEvolution of protein moleculesMammalian protein metabolismYear: 1969New York: Academic Press21123  
Felsenstein J,Evolutionary trees from DNA sequences: a maximum likelihood approachJ Mol EvolYear: 19811736837610.1007/BF017343597288891  
Hasegawa M,Kishino H,Yano TA,Dating of the humanape splitting by a molecular clock of mitochondrial DNAJ Mol EvolYear: 19852216017410.1007/BF021016943934395  
Nei N,Kumar S,Molecular Evolution and PhylogeneticsYear: 2000Oxford: Oxford University Press  
Watterson GA,On the number of segregating sites in genetic models without recombinationTheor Popul BiolYear: 1975725627610.1016/00405809(75)9002091145509  
Sherwin WB,Jabot F,Rush R,Rossetto M,Measurement of biological information with applications from genes to landscapesMol EcolYear: 2006152857286910.1111/j.1365294X.2006.02992.x16911206  
Wakeley J,Distinguishing migration from isolation using the variance of pairwise differencesTheor Popul BiolYear: 19964936938610.1006/tpbi.1996.00188693431  
Nei M,Li W,Mathematical model for studying variation in terms of restriction endonucleasesProceedings of the National Academy of Sciences of the United States of AmericaYear: 1979765269527310.1073/pnas.76.10.5269291943  
Wakeley J,The variance of pairwise nucleotide differences in two populations with migrationTheor Popul BiolYear: 199649395710.1006/tpbi.1996.00028813013  
Bazin E,Dawson KJ,Beaumont MA,Likelihoodfree Inference of Population Structure and Local Adaptation in a Bayesian Hierarchical ModelGeneticsYear: 2010 genetics. 109.112391. 20382835  
Papoulis A,Probability, Random Variables, and Stochastic ProcessesYear: 19842New York: McGrawHill  
Takahata N,Nei M,Gene genealogy and variance of interpopulational nucleotide differencesGeneticsYear: 19851103253444007484  
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. persimilisGeneticsYear: 200416774776010.1534/genetics.103.02418215238526  
Cook SR,Gelman A,Rubin DB,Validation of Software for Bayesian Models Using Posterior QuantilesJournal of Computational and Graphical StatisticsYear: 200615367569210.1198/106186006X136976  
Jennings WB,Edwards SV,Speciational history of australian grass finches (poephila) inferred from thirty gene treesEvolutionYear: 2005592033204716261740  
Lee JY,Edwards SV,Divergence Across Australia's Carpentarian Barrier: Statistical Phylogeography of the RedBacked Fairy Wren (Malurus melanocephalus)EvolutionYear: 200862123117313410.1111/j.15585646.2008.00543.x19087188  
Sousa VC,Fritz M,Beaumont MA,Chikhi L,Approximate Bayesian Computation Without Summary Statistics: The Case of AdmixtureGeneticsYear: 200918141507151910.1534/genetics.108.09812919189952  
Rosenberg NA,Feldman MW,Slatkin M, Veuille MThe relationship between coalescence times and population divergence timesModern Developments in Theoretical Population GeneticsYear: 2002Oxford: University Press130164  
Arbogast BS,Edwards SV,Wakeley J,Beerli P,Slowinski JB,Estimating divergence times from molecular data on phylogenetic and population genetic timescalesAnnu Rev Ecol SystYear: 20023370774010.1146/annurev.ecolsys.33.010802.150500  
Nielsen R,Wakeley J,Distinguishing migration from isolation: A Markov chain Monte Carlo approachGeneticsYear: 2001158288589611404349  
Carstens BC,Stoute HN,Reid NM,An informationtheoretical approach to phylogeographyMol EcolYear: 200918204270428210.1111/j.1365294X.2009.04327.x19765225  
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 AfricaCurr BiolYear: 200919431231810.1016/j.cub.2008.12.04919200724  
Kass RE,Raftery A,Bayes factorsJournal of the American Statistical AssociationYear: 19959077379510.2307/2291091  
Rannala B,Yang ZH,Bayes Estimation of Species Divergence Times and Ancestral Population Sizes Using DNA Sequences From Multiple LociGeneticsYear: 20031641645165612930768  
Lee JY,Edwards SV,Divergence Across Australia's Carpentarian Barrier: Statistical Phylogeography of the RedBacked Fairy Wren (Malurus melanocephalus)EvolutionYear: 200962123117313410.1111/j.15585646.2008.00543.x  
Knowles LL,Statistical PhylogeographyAnnual Review of Ecology, Evolution, and SystematicsYear: 200940159361210.1146/annurev.ecolsys.38.091206.095702  
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 afterMol Phylogen EvolYear: 201054129130110.1016/j.ympev.2009.09.016  
Moyle LC,Ecological and evolutionary genomics in the wild tomatoes (solanum sect. Lycopersicon)EvolutionYear: 200862122995301310.1111/j.15585646.2008.00487.x18752600  
Dolman G,Moritz C,A multilocus perspective on refugial isolation and divergence in rainforest skinks (carlia)EvolutionYear: 200660357358216637502  
Hurt C,Anker A,Knowlton N,A multilocus test of simultaneous divergence across the isthmus of panama using snapping shrimp in the genus alpheusEvolutionYear: 200963251453010.1111/j.15585646.2008.00566.x19154357  
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 PairsMol Biol EvolYear: 20072461423143810.1093/molbev/msm06617404398  
Becquet Cl,Przeworski M,Learning about Modes of Speciation by Computational ApproachesEvolutionYear: 200963102547256210.1111/j.15585646.2009.00662.x19228187 
Figures
Article Categories:

Previous Document: Does osteoporosis predispose falls? A study on obstacle avoidance and balance confidence.
Next Document: The impact of audit and feedback on nodal harvest in colorectal cancer.