BACKGROUND
a common challenge faced by empirical researchers in studies of ecological communities is to identify individuals at the species level from limited information collected from a broad taxonomic range of organisms.
in many cases, useful taxonomic keys for particular groups or regions are not available.
this is because many diverse groups are morphologically cryptic, contain many undescribed taxa, or existing taxonomic literature is conflicting, an issue referred to as the “taxonomic impediment”
<cit> .
in these cases, short dna sequence tags  are frequently surveyed because they can be rapidly and inexpensively collected
<cit> .
dna barcoding initiatives aim to connect these sequence tags to taxa validated by expert taxonomists
<cit> , but at present this is not possible for most groups.
as a result, diversity must frequently be quantified in the absence of a low-level taxonomic framework.
in order to accomplish this, observed dna sequences must be clustered into putative species.
while the delimitation of species is a complex philosophical and biological problem
<cit> , species concepts widely share the idea that species are independently evolving metapopulation lineages
<cit> .
this provides a justification for using genetic data  as the primary data for the diagnosis of these lineages, as they contain the signal of historical processes involved in lineage divergence
<cit> .
as a caveat, lineages identified in this way will not necessarily meet the criteria for species status under any given species concept, such as reproductive isolation from other such lineages, or exhibit morphological, ecological or behavior divergence.
methods used for delimitation of species from barcode data are a subset of those developed for the larger problem of species delimitation.
they can be considered species discovery methods because they must be functional in the absence of good a priori taxonomic information
<cit> .
this contrasts with validation methods , which test specific hypotheses of species status, and assignment methods, which assign unknown individuals to existing species .
of the handful of approaches typically used to discover species limits using genetic data, thresholds based on pairwise sequence distances among individuals are perhaps most commonly applied to cluster sequences into putative species
<cit> .
these methods identify some level of sequence divergence beyond which two individuals cannot be considered the same species.
distance threshold methods have been criticized because they do not account for evolutionary processes
<cit> , and the uncertainty in selecting an appropriate threshold
<cit> , which relies on an observable “barcode gap” between pairwise intraspecific and interspecific dna sequence distances .
pons et al.
<cit>  introduced a model-based alternative to distance threshold methods.
the model, the general mixed yule-coalescent , takes a phylogenetic tree estimated from dna sequence data and assumes that the branching points in the tree correspond to one of two events: divergence events between species level taxa , or coalescent events between lineages sampled from within species .
because the rate of coalescence within species is expected to be dramatically greater than the rate of cladogenesis, the gmyc aims to find the demarcation between these types of branching.
this model has been shown to be useful in several empirical studies
<cit> .
because it is based on a likelihood function that directly models evolutionary processes of interest, it provides a means to ameliorate some of the criticisms leveled at threshold methods.
notably, it has allowed for quantification of uncertainty in delimitation of species
<cit>  and avoids the use of non-independent pairwise sequence distances  as data.
the gmyc model as presently implemented, however, does not account for three potentially large sources of error.
first, it is widely recognized that a variety of factors can cause the genealogy from a particular locus to be discordant with the true history of speciation
<cit> , and the gmyc, like all methods based on a single locus, can be mislead by this discordance.
second, there may be error in the model estimates.
under certain circumstances, the transition from speciation events to coalescent events may be indistinct  causing the model to have a wide confidence interval.
a recent implementation by powell
<cit>  accounts for uncertainty in the threshold parameter and produces model-averaged species limits, but uses point estimates for the other parameters.
finally, phylogenetic error can diminish the accuracy of delimitation results.
the gmyc model requires the user to input a point estimate of the phylogenetic tree and inference is premised on the accuracy of this point estimate.
diversity studies using sequence tags, however, typically use relatively short loci that yield estimates of topology and branch lengths that may have high levels of uncertainty.
this uncertainty could influence the accuracy of the model.
in order to address the second and third potential sources of error, we introduce a bayesian implementation of this model with flexible prior distributions in the statistical scripting language r
<cit> .
it accounts for the error in phylogenetic estimation and uncertainty in model parameters by integrating over uncertainty in tree topology and branch lengths and in the parameters of the model via markov chain monte carlo simulation
<cit> .
it produces marginal posterior probabilities for species that are independent of these factors along with output characterizing the posterior distribution that is suitable for downstream analyses of community structure accounting for uncertainty in species limits and phylogeny using r packages such as picante
<cit> , vegan
<cit> , and ape
<cit> .
we also conduct simulation tests to evaluate the performance of the model and re-analyze a dataset previously analyzed with the likelihood version of the model.
methods
model
given an ultrametric phylogenetic tree estimated from a set of sequences consisting of multiple species and multiple individuals within species, the gmyc model decomposes the tree into its component waiting times between branching events.
these waiting times are the data to be modeled
<cit> .
the parameter of primary interest is a threshold parameter before which the waiting times are modeled according to a yule process and after which, according to k intraspecific coalescent processes, where k is the number of species with more than one unique sequence sampled.
the likelihood of a waiting time under a yule model is:
lxi=λnipe−λnipxi
where the waiting times  are assumed to be exponentially distributed and a function of: the branching rate , the number of lineages in the interval , and a rate change parameter that accounts for the possibility of increasing or decreasing diversification rate with time .
the likelihood  of a waiting time under the coalescent model is:
lxi=λnini−1pe−λnini−1pxi
where the branching rate  can be interpreted as 1/neμ  and the rate change parameter  can be interpreted as population size change over time
<cit> .
the gmyc model combines the above models, and the likelihood of the full model is calculated by assigning lineages in each waiting interval to either the yule process or one of the coalescent processes such that:
b=λk+1ni,k+1pk+1+∑j=l,kλjnijnij−1pj
making the full likelihood of a waiting time:
lxi=be−bxi
where k indexes the branching processes , λk+ <dig> and pk+ <dig> are the branching rate and rate change parameters for the yule process, and λj and pj are the branching rate and population size change parameters for the coalescent process.
following pons et al.
<cit> , we constrain λj and pj to be identical across coalescent processes.
the number of lineages assigned to the yule and coalescent processes in each waiting interval are ni,k+ <dig> and ni,j, respectively.
assignment of lineages in this case is determined by the selection of a threshold.
because the sequence data employed in these analyses are typically from short fragments, and thus likely to yield trees with high levels of uncertainty in topology and branch lengths we implemented this model in a bayesian statistical framework.
it eliminates the reliance on point estimates of the phylogeny and model parameters and by estimating the marginal probabilities of the identity of species, allows one to incorporate that uncertainty in downstream analyses.
our implementation operates as follows.
first, the posterior distribution of trees and branch lengths are characterized using beast
<cit> .
second, a post-burn-in sample of the trees from that posterior distribution is taken, and for each tree, an mcmc simulation is conducted to characterize the posterior distribution of the gmyc model.
following pons et al., we did not allow the λ parameters to be freely estimated, but used a maximum likelihood estimator
<cit> .
the non-normalized posterior density function is as follows:
pt,λj,λk+ <dig> pjpk+1|τ∝pt,λj,λk+ <dig> pjpk+1pτ|t,λjλk+ <dig> pjpk+ <dig>
where t is the threshold.
because each mcmc evaluates the posterior of the gmyc conditioned on the tree, pooled samples from analyses of many trees sampled from the tree posterior result in a posterior probability distribution of species delimitations conditioned only on the sequence data, the phylogenetic model and the gmyc model.
simulation testing
we evaluated the utility of this implementation of the gmyc using three simulation experiments.
in each, we simulated gene trees from species trees using ms
<cit> .
all species trees contained  <dig> species and were generated via a yule process in mesquite
<cit> , randomly sampling  <dig> species from a tree of  <dig> species.
this design was intended to mimic environmental samples of a given taxon, which would not be expected to contain all species in a clade.
in the first experiment we examined the effect of tree depth on model accuracy.
we simulated  <dig> species trees as above and scaled them to four different depths .
when considering how the results translate to haploid, maternally inherited organellar dna, the equivalent tree depths are halved  and n becomes the effective number of females in the population.
we then simulated a single gene tree from each species tree at each depth, sampling five alleles per species.
for each of these trees we sampled from the posterior for  <dig>  generations, discarding the first  <dig>  generations as burn-in and thinning every  <dig> generations, assessed stationarity by examining plots of the parameters by eye, and characterized the posterior distribution of the threshold parameter, which determines the species limits given a tree.
priors on all parameters were uniform distributions; in the case of the threshold parameter, from u and for the p parameters u.
in the second experiment we looked at the influence of sampling.
the species trees with a depth of  <dig> n from the first experiment were used with four different sampling schemes:  <dig> alleles per species,  <dig> alleles per species,  <dig> alleles per species, and a random number of alleles per species, drawn from a lognormal distribution, with a mean and standard deviation of  <dig> .
we used the lognormal distribution because it approximates some real species-abundance distributions
<cit>  that might be observed in actual species delimitation datasets.
we conducted the analyses as in the first experiment.
in the third experiment, we tested the effect of nucleotide sampling and tree estimation on the accuracy of the model .
we selected  <dig> of the simulated gene trees from  <dig> species trees scaled to  <dig> n generations for which the confidence intervals in the analysis overlapped the true value of  <dig> species.
we then simulated dna sequences on those gene trees of  <dig> bp,  <dig> bp,  <dig> bp and  <dig> bp using seq-gen
<cit> .
we assumed θ =  <dig>   and an hky + g model.
we characterized the posterior distribution of trees using the true model of sequence evolution and a strict clock model in beast.
we pruned all identical sequences and ran beast for  <dig> million generations, discarding the first million as burn-in, at which point all parameters for all replicates had effective sample sizes above  <dig> and most above  <dig>  we then ran independent gmyc mcmc analyses on  <dig> trees sampled every  <dig>  generations from the beast posterior distribution of trees, pooled the results and characterized the marginal posterior distribution of the threshold parameter compared to the distribution produced using the true tree.
empirical data analyses
to illustrate how this implementation of the gmyc could be applied; we downloaded from genbank and reanalyzed the dataset from pons et al., the original publication of the gmyc .
we first pruned each alignment to consist only of unique sequences.
since we are not using a true genealogy sampler , identical sequences will result in many zero-length branches at the tip of the tree, and will cause the model to over-partition the dataset.
we then applied a phylogenetic clock model to estimate the posterior distribution of ultrametic trees using beast.
we partitioned models of sequence evolution by codon, and also by gene when multiple genes were used, applying models of sequence evolution selected by dt modsel
<cit> .
we accepted that runs converged when all parameters reached ess values greater than  <dig> and checked that posterior distributions differed from priors.
we explored the use of different tree priors as it is conceivable that in cases where branch-length information is lacking, the prior could strongly influence the posterior.
for bayesian gmyc mcmc analyses, we ran each tree for  <dig>  generations, discarding the first  <dig> as burn-in and sampling every  <dig> generations.
using  <dig> trees sampled from the beast posterior distribution of trees, this resulted in  <dig> samples.
we selected this length of markov-chain because preliminary analyses suggested that stationarity was usually achieved by  <dig> generations.
we compared the posterior distribution from sampling multiple trees to that from the maximum clade credibility tree and examined the effect of changing the prior on the yule rate change parameter .
we compared the posterior distribution to the point estimate produced by the likelihood version of the model, and to the akaike weights
<cit>  of each threshold point.
RESULTS
simulation tests
we first tested the influence of tree depth on model performance.
when deeper trees are simulated, coalescent and yule branching processes are expected to occur on more distinct time scales, and thus in general the model should perform better.
the influence of tree depth is actually confounded by two issues, however.
first, as the tree depth becomes shallower the implied rate of speciation increases because all trees contain  <dig> species.
if the rate of speciation approaches the rate of coalescence within species, then a sharp transition between processes should not be detectable.
second, as the implied rate of speciation increases, more species originate relatively recently.
the expected time to coalescence for a diploid, panmictic population is  <dig> n generations.
cladogenic events occurring more recently than this are expected to be increasingly difficult to delimit for two reasons: they are more likely to yield species that are not monophyletic and thus impossible to accurately identify under this model, and the most recent common ancestor  of the daughter species is more likely to occur more recently than the threshold point.
assuming species monophyly, the expected time to the mrca for two species that diverged  <dig> n generations ago is  <dig> n generations.
therefore all probability should be on thresholds older than  <dig> n generations, and most on thresholds older than  <dig> n generations.
again, when considering maternally inherited, haploid, organellar dna, equivalent times in n generations are halved, and n becomes the effective number of females in the population.
this would give an expected time to mrca of  <dig> n generations.
the results of this test are dramatic .
there is a clear increase in accuracy as well as a decrease in the range of the 95% highest posterior density interval  with increasing tree depth.
at shorter tree depths, the model’s performance diverged from expectations.
when trees are short, larger numbers of species have divergence times younger than  <dig> n generations and thus should not be detectable under the model.
therefore we expected the number of species delimited to be smaller than for deeper trees.
we did not observe this.
instead, replicate posterior means became widely scattered around the true value of  <dig> species.
95% hpds also did not uniformly increase with decreasing tree depth, instead, a wide distribution of hpds were observed.
we also expected threshold times to have a minimum value of approximately  <dig> n generations.
at deeper tree depths this is observed, with  <dig> % and  <dig> % of mcmc steps sampling thresholds younger than  <dig> n generations for tree depths of  <dig> n and  <dig> n, respectively.
however, at  <dig> n and  <dig> n tree depths, 20% and 40% of mcmc steps sampled thresholds younger than  <dig> n generations.
these results indicate that the model performs well under demographic or sampling conditions that result in coalescent and yule processes occurring on very different time scales.
it does not, however, perform optimally when those conditions are not met.
ideally one would hope that as inference of the threshold point became more difficult, that the 95% hpds would increase, but still encompass the true value 95% of the time.
this is not the case at the  <dig> n and  <dig> n tree depths.
hpds generally become broader, but for increasing numbers of simulation replicates, they fail to encompass the true value.
<dig> species arising in  <dig> n generations constitutes a very rapid radiation, with an average of 89% of branches in the species tree shorter than the expected population coalescence time of  <dig> n generations.
failure to accurately assess credibility intervals in this case is likely because in this area of parameter space, the gmyc is no longer an accurate approximation of the real branching process in the gene tree.
rather than there being a threshold between coalescent and speciation branching processes, the two processes are intermixed because there is little time for the independent evolution of lineages prior to speciation.
note that these conditions will cause any dna barcode-based method of species discovery to fail and will also challenge more realistic models utilizing multilocus data and prior information on population assignment.
next we examined the effect of intraspecific sampling.
because the data points used by the model are waiting times between branching events, we expected that with  <dig> species, we would not need extremely high sampling to accurately characterize the model, and that the distribution of samples among species would not be particularly important.
our expectations were met.
we found that sampling of  <dig> individuals per species yielded poor results .
replicate posterior means showed a strong bias towards inference of a large number of species.
sampling of greater than  <dig> individuals per species provided an improvement in the accuracy of the posterior means, but no change in the 95% hpd range.
all sampling schemes greater than  <dig> individuals per species appeared to yield similar results, including the more realistic condition of a lognormal distribution of alleles among species in which a large number  of species are represented by singletons.
while delimiting rare species, particularly from single specimens, is a challenge faced by taxonomists
<cit> , our results suggest that the gmyc method may efficiently delimit species represented by singletons by calibrating the divergence threshold using data from better represented species.
finally, we tested the effects of nucleotide sampling and the incorporation of phylogenetic uncertainty.
we expected to find wider hpds with less sequence data, as uncertainty, particularly in branch lengths should be greater.
we found a mild reduction in accuracy of the posterior means with up to  <dig> bp of sequence, but after that, posterior means converged on those of the true tree.
the 95% hpds improved with the addition of more sequence, but had not quite converged on those estimated from the true tree, even at  <dig> bp .
this suggests that uncertainty in phylogenetic estimation, particularly in typical dna barcode datasets of ~ <dig> bp will contribute substantially to uncertainty in species delimitation.
three factors that could influence the accuracy of the model that were not explored here: migration, population substructure and selection.
papadopoulou et al.
<cit>  examined the effects of migration on the formation of detectable gmyc clusters.
they simulated datasets under an island model and found that even very low levels of migration  caused likelihood ratio tests to fail to reject the null model of a single branching process.
they interpreted this as evidence that the likelihood implementation of the model is conservative and will not infer species at all unless they are strongly isolated.
papadopoulou et al.’s simulations assumed complete demic sampling, but lohse
<cit>  conducted simulations showing that under moderate migration rates  and with a large proportion  of demes unsampled that spurious, significant clusters could be inferred from the true gene genealogies.
in his simulations, lohse showed that when  <dig> demes were sampled from a metapopulation consisting of  <dig> demes, that an average of  <dig> species were inferred, and 80% of replicates rejected the null model.
wakeley
<cit>  described the genealogical pattern resulting from such a process as having two phases that occur on very different time scales: a scattering phase, in which there is rapid coalescence and migration in local demes, and a collecting phase that begins when each remaining lineage is in its own deme and takes a very long time.
in this case the gmyc might see the scattering phase as the “coalescent process” and the collecting phase as the “yule” process.
further exploration of this issue is likely to be important, particularly if the gmyc is applied to phylogenetic samples with deep phylogeographic sampling.
while lohse shows convincingly that this interaction of parameter space with sampling can mislead the gmyc, it is not clear to what extent these problematic areas of parameter space exist in real datasets.
we simulated  <dig> genealogies using ms under the conditions above and observed that the average time to coalescence of all lineages was  <dig>  n generations , with the scattering phase taking the first 4- <dig> n generations.
if we assumed that these  <dig> demes were species level taxa, each with θ =  <dig> , we would expect to observe gmyc clusters with mrcas at a depth of  <dig> - <dig>  substitutions per site and the mrca for all lineages at  <dig>  substitutions per site.
it is unlikely that the collecting phase would have the time to play out under this scenario, as it would take nearly  <dig> million years at a mutation rate of  <dig>  substitutions per site per million years.
if, by contrast, we assume that these demes represent populations at a smaller scale, each with a theta of  <dig> / <dig>  then we would expect mrcas of delimited clusters to be at a depth of  <dig> - <dig>  substitutions per site.
with a typical dna barcode or short mitochondrial dna set of 650- <dig> bp, the scattering phase would be undetectable.
the mrca for all lineages would occur at  <dig>  substitutions per site.
unless this process was considered in the context of a larger species tree, it is unlikely that the gmyc would identify a significant branching threshold.
empirical analyses
we reanalyzed the empirical data used by pons et al.
to illustrate the original formulation of the gmyc so as to provide a direct comparison of the implementations using representative data.
the beast run converged after  <dig> million generations and we discarded  <dig>  million trees as burn-in.
the estimate of the standard deviation of the lognormal distribution of rates did not overlap  <dig>  so we could not use a strict clock with these data.
when using samples of trees from the beast posterior distribution, the mean number of species estimated by the bayesian gmyc was  <dig> and the 95% hpd ranged from  <dig> to  <dig>  the rate change parameter for the yule process ranged as high as  <dig> .
in this model, the fold change in speciation rate from the root to the last speciation event is equal to np/n where n is the estimated number of species and p is the yule rate change parameter.
in this case, given  <dig> species, p =  <dig>  allows for a 30-fold speciation rate increase.
we thought this might be unrealistically high , and thus re-ran the analysis with the prior distribution set to u, or a maximum 2-fold increase.
this minimally influenced the results.
we also analyzed the maximum clade credibility tree under both priors, and using the likelihood implementation of the gmyc.
we compared the results of likelihood and bayesian analyses by calculating akaike weights for each possible threshold in the likelihood analysis.
akaike weights are the relative likelihoods of a set of models and thus suited for qualitatively comparing results among these analyses.
the bayesian gmyc analysis of the maximum clade credibility tree yielded a mean of  <dig> species and a narrower 95% hpd of 38– <dig> species.
when applied to the maximum clade credibility tree, the u prior did change the results, yielding a posterior mean of  <dig> species and a 95% hpd of  <dig> to  <dig>  the maximum likelihood analysis resulted in  <dig> species with a 95% confidence interval of  <dig> to  <dig> species.
the results of the likelihood and bayesian analyses of the maximum clade credibility tree are very similar , particularly when the prior u is placed on the yule parameter, but when sampling trees, the posterior distribution is substantially broader .
at least some of this uncertainty stems from variation in topology.
plots of pairwise probabilities of conspecificity demonstrate this .
the probability distribution of the number of species in the sample is also wider for the posterior distribution than for the akaike weights .
our bayesian approach is similar in spirit to the model-averaging approach of powell
<cit>  in that its goal is to make the inference of species limits and analyses based on them more robust to uncertainty.
there are, however, three major differences.
first, we take into account phylogenetic uncertainty.
as indicated by our results, this is perhaps the most influential difference, although in theory it could also be accounted for using model averaging.
second, our bayesian method requires the specification of prior knowledge.
depending on a researcher’s comfort with assigning prior probability distributions to the rate and threshold parameters, this is either an advantage or a disadvantage.
finally, the treatment of nuisance parameters  is fundamentally different.
while in the model-averaging approach, inferences are conditioned on maximum-likelihood point estimates of nuisance parameters, the bayesian approach integrates out nuisance parameters, giving marginal probabilities of species limits.
these final two differences are intrinsic to bayesian inference, and researchers choosing among methods will need to consider their choice of statistical paradigm.
we note that our confidence intervals in all analyses  are far wider than those of pons et al., which were 46– <dig> species.
this is most likely because of the difference in obtaining ultrametric trees.
pons et al.
used non-parametric rate smoothing
<cit>  on a maximum likelihood estimate of the tree.
while they achieved sensible results with narrower confidence limits than ours, we nevertheless advocate an approach that samples trees and fits them to a clock model using a parametric method.
this allows for a full accounting of uncertainty associated with phylogenetic estimation, albeit at the cost of some precision.
CONCLUSIONS
our results demonstrate that the bayesian implementation of the gmyc model is reasonably reliable given two caveats.
first, the length of the dna sequence is important.
we found that when we sampled only  <dig> bp, or only  <dig> alleles per species, that the performance of the model declined strongly.
second, the model is only useful when the underlying history of the species under consideration lies in particular regions of parameter space.
species that have recently diverged, or clades undergoing rapid radiation are unlikely to be identifiable under the model.
in the latter case, the model may provide misleading estimates and confidence.
cases such as these, however, may be recognizable because the results may be highly unexpected in the context of other sources of data such as morphology or geography.
our implementation of the model provides two main improvements over the original.
first, it allows the specification of prior probabilities on model parameters.
it is our experience that very high values of the yule process rate change parameter sometimes have high likelihood and result in high uncertainty in the threshold parameter .
these high values may be biologically unrealistic, and the specification of an informative prior can reduce the posterior probability of those areas and produce a more accurate estimate of diversity.
second, it allows for the characterization of species limits without use of a point estimate of the phylogeny.
we know that many datasets are associated with substantial uncertainty owing to limited sequence data collection.
the bayesian gmyc method provides marginal probabilities of species identities and will allow downstream estimates of species diversity and community structure  to account for uncertainty underlying species designations.
an important future direction for this work is to implement the multiple-threshold version of the model proposed by monaghan et al.
<cit> , which can account for greater variation in divergence times and effective population sizes than the model implemented here.
it has been shown to provide a better fit to some datasets
<cit> , but will require implementation of a more complex reversible-jump mcmc that allows proposals that change the number of parameters in the model.
it is widely acknowledged that single-locus data are not optimal for the inference of phylogeny, historical demography, or species limits
<cit> .
nevertheless, vast amounts of biological diversity remain undescribed at the level of species, and this limits our ability to understand the evolutionary history of our planet and its current ecological functioning.
available alternative means of describing species diversity, either from molecular or morphological data have major drawbacks in that they are time consuming, expensive, or subject to their own biases and inaccuracies.
single-locus data for many groups are currently being generated on a large scale, and we advocate making the best of this data.
we believe that under certain conditions, the gmyc model can be useful, and that a bayesian framework accounting for uncertainty is most appropriate for these data.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
nr and bc developed the original concept and designed the simulation experiments.
nr wrote all software and conducted the analyses described in the manuscript.
nr and bc interpreted results and wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
figures s <dig>  s <dig>  s <dig>  these figures display the distribution of mcmc samples for each treatment and each replicate within treatments for simulated data.
s <dig> is results from the tree depth simulation, s <dig> is the results from the allele sampling simulation and s <dig> is the results from the nucleotide sampling simulation.
click here for file
acknowledgements
we thank the national science foundation  for funding aspects of this work.
we thank jeremy brown for conversations that initiated this research, and members of the carstens lab  at louisiana state university for discussions related to and comments on this manuscript.
we thank dr. timothy barraclough and two anonymous reviewers for helpful correspondence regarding this work and comments on drafts of the manuscript.

 a common challenge faced by empirical researchers in studies of ecological communities is to identify individuals at the species level from limited information collected from a broad taxonomic range of organisms. in many cases, useful taxonomic keys for particular groups or regions are not available. this is because many diverse groups are morphologically cryptic, contain many undescribed taxa, or existing taxonomic literature is conflicting, an issue referred to as the “taxonomic impediment”. in these cases, short dna sequence tags are frequently surveyed because they can be rapidly and inexpensively collected. dna barcoding initiatives aim to connect these sequence tags to taxa validated by expert taxonomists, but at present this is not possible for most groups. as a result, diversity must frequently be quantified in the absence of a low-level taxonomic framework. in order to accomplish this, observed dna sequences must be clustered into putative species. while the delimitation of species is a complex philosophical and biological problem, species concepts widely share the idea that species are independently evolving metapopulation lineages. this provides a justification for using genetic data as the primary data for the diagnosis of these lineages, as they contain the signal of historical processes involved in lineage divergence. as a caveat, lineages identified in this way will not necessarily meet the criteria for species status under any given species concept, such as reproductive isolation from other such lineages, or exhibit morphological, ecological or behavior divergence. methods used for delimitation of species from barcode data are a subset of those developed for the larger problem of species delimitation. they can be considered species discovery methods because they must be functional in the absence of good a priori taxonomic information. this contrasts with validation methods, which test specific hypotheses of species status, and assignment methods, which assign unknown individuals to existing species. of the handful of approaches typically used to discover species limits using genetic data, thresholds based on pairwise sequence distances among individuals are perhaps most commonly applied to cluster sequences into putative species. these methods identify some level of sequence divergence beyond which two individuals cannot be considered the same species. distance threshold methods have been criticized because they do not account for evolutionary processes, and the uncertainty in selecting an appropriate threshold, which relies on an observable “barcode gap” between pairwise intraspecific and interspecific dna sequence distances. pons et al. introduced a model-based alternative to distance threshold methods. the model, the general mixed yule-coalescent, takes a phylogenetic tree estimated from dna sequence data and assumes that the branching points in the tree correspond to one of two events: divergence events between species level taxa, or coalescent events between lineages sampled from within species. because the rate of coalescence within species is expected to be dramatically greater than the rate of cladogenesis, the gmyc aims to find the demarcation between these types of branching. this model has been shown to be useful in several empirical studies. because it is based on a likelihood function that directly models evolutionary processes of interest, it provides a means to ameliorate some of the criticisms leveled at threshold methods. notably, it has allowed for quantification of uncertainty in delimitation of species and avoids the use of non-independent pairwise sequence distances as data. the gmyc model as presently implemented, however, does not account for three potentially large sources of error. first, it is widely recognized that a variety of factors can cause the genealogy from a particular locus to be discordant with the true history of speciation, and the gmyc, like all methods based on a single locus, can be mislead by this discordance. second, there may be error in the model estimates. under certain circumstances, the transition from speciation events to coalescent events may be indistinct causing the model to have a wide confidence interval. a recent implementation by powell accounts for uncertainty in the threshold parameter and produces model-averaged species limits, but uses point estimates for the other parameters. finally, phylogenetic error can diminish the accuracy of delimitation results. the gmyc model requires the user to input a point estimate of the phylogenetic tree and inference is premised on the accuracy of this point estimate. diversity studies using sequence tags, however, typically use relatively short loci that yield estimates of topology and branch lengths that may have high levels of uncertainty. this uncertainty could influence the accuracy of the model. in order to address the second and third potential sources of error, we introduce a bayesian implementation of this model with flexible prior distributions in the statistical scripting language r. it accounts for the error in phylogenetic estimation and uncertainty in model parameters by integrating over uncertainty in tree topology and branch lengths and in the parameters of the model via markov chain monte carlo simulation. it produces marginal posterior probabilities for species that are independent of these factors along with output characterizing the posterior distribution that is suitable for downstream analyses of community structure accounting for uncertainty in species limits and phylogeny using r packages such as picante, vegan, and ape. we also conduct simulation tests to evaluate the performance of the model and re-analyze a dataset previously analyzed with the likelihood version of the model. methods model given an ultrametric phylogenetic tree estimated from a set of sequences consisting of multiple species and multiple individuals within species, the gmyc model decomposes the tree into its component waiting times between branching events. these waiting times are the data to be modeled. the parameter of primary interest is a threshold parameter before which the waiting times are modeled according to a yule process and after which, according to k intraspecific coalescent processes, where k is the number of species with more than one unique sequence sampled. the likelihood of a waiting time under a yule model is: lxi=λnipe−λnipxi where the waiting times are assumed to be exponentially distributed and a function of: the branching rate, the number of lineages in the interval, and a rate change parameter that accounts for the possibility of increasing or decreasing diversification rate with time. the likelihood of a waiting time under the coalescent model is: lxi=λnini−1pe−λnini−1pxi where the branching rate can be interpreted as 1/neμ and the rate change parameter can be interpreted as population size change over time. the gmyc model combines the above models, and the likelihood of the full model is calculated by assigning lineages in each waiting interval to either the yule process or one of the coalescent processes such that: b=λk+1ni,k+1pk+1+∑j=l,kλjnijnij−1pj making the full likelihood of a waiting time: lxi=be−bxi where k indexes the branching processes, λk+ and pk+ are the branching rate and rate change parameters for the yule process, and λj and pj are the branching rate and population size change parameters for the coalescent process. following pons et al., we constrain λj and pj to be identical across coalescent processes. the number of lineages assigned to the yule and coalescent processes in each waiting interval are ni,k+ and ni,j, respectively. assignment of lineages in this case is determined by the selection of a threshold. because the sequence data employed in these analyses are typically from short fragments, and thus likely to yield trees with high levels of uncertainty in topology and branch lengths we implemented this model in a bayesian statistical framework. it eliminates the reliance on point estimates of the phylogeny and model parameters and by estimating the marginal probabilities of the identity of species, allows one to incorporate that uncertainty in downstream analyses. our implementation operates as follows. first, the posterior distribution of trees and branch lengths are characterized using beast. second, a post-burn-in sample of the trees from that posterior distribution is taken, and for each tree, an mcmc simulation is conducted to characterize the posterior distribution of the gmyc model. following pons et al., we did not allow the λ parameters to be freely estimated, but used a maximum likelihood estimator. the non-normalized posterior density function is as follows: pt,λj,λk+ pjpk+1|τ∝pt,λj,λk+ pjpk+1pτ|t,λjλk+ pjpk+ where t is the threshold. because each mcmc evaluates the posterior of the gmyc conditioned on the tree, pooled samples from analyses of many trees sampled from the tree posterior result in a posterior probability distribution of species delimitations conditioned only on the sequence data, the phylogenetic model and the gmyc model. simulation testing we evaluated the utility of this implementation of the gmyc using three simulation experiments. in each, we simulated gene trees from species trees using ms. all species trees contained species and were generated via a yule process in mesquite, randomly sampling species from a tree of species. this design was intended to mimic environmental samples of a given taxon, which would not be expected to contain all species in a clade. in the first experiment we examined the effect of tree depth on model accuracy. we simulated species trees as above and scaled them to four different depths. when considering how the results translate to haploid, maternally inherited organellar dna, the equivalent tree depths are halved and n becomes the effective number of females in the population. we then simulated a single gene tree from each species tree at each depth, sampling five alleles per species. for each of these trees we sampled from the posterior for generations, discarding the first generations as burn-in and thinning every generations, assessed stationarity by examining plots of the parameters by eye, and characterized the posterior distribution of the threshold parameter, which determines the species limits given a tree. priors on all parameters were uniform distributions; in the case of the threshold parameter, from u and for the p parameters u. in the second experiment we looked at the influence of sampling. the species trees with a depth of n from the first experiment were used with four different sampling schemes: alleles per species, alleles per species, alleles per species, and a random number of alleles per species, drawn from a lognormal distribution, with a mean and standard deviation of. we used the lognormal distribution because it approximates some real species-abundance distributions that might be observed in actual species delimitation datasets. we conducted the analyses as in the first experiment. in the third experiment, we tested the effect of nucleotide sampling and tree estimation on the accuracy of the model. we selected of the simulated gene trees from species trees scaled to n generations for which the confidence intervals in the analysis overlapped the true value of species. we then simulated dna sequences on those gene trees of bp, bp, bp and bp using seq-gen. we assumed θ =  and an hky + g model. we characterized the posterior distribution of trees using the true model of sequence evolution and a strict clock model in beast. we pruned all identical sequences and ran beast for million generations, discarding the first million as burn-in, at which point all parameters for all replicates had effective sample sizes above and most above we then ran independent gmyc mcmc analyses on trees sampled every generations from the beast posterior distribution of trees, pooled the results and characterized the marginal posterior distribution of the threshold parameter compared to the distribution produced using the true tree. empirical data analyses to illustrate how this implementation of the gmyc could be applied; we downloaded from genbank and reanalyzed the dataset from pons et al., the original publication of the gmyc. we first pruned each alignment to consist only of unique sequences. since we are not using a true genealogy sampler, identical sequences will result in many zero-length branches at the tip of the tree, and will cause the model to over-partition the dataset. we then applied a phylogenetic clock model to estimate the posterior distribution of ultrametic trees using beast. we partitioned models of sequence evolution by codon, and also by gene when multiple genes were used, applying models of sequence evolution selected by dt modsel. we accepted that runs converged when all parameters reached ess values greater than and checked that posterior distributions differed from priors. we explored the use of different tree priors as it is conceivable that in cases where branch-length information is lacking, the prior could strongly influence the posterior. for bayesian gmyc mcmc analyses, we ran each tree for generations, discarding the first as burn-in and sampling every generations. using trees sampled from the beast posterior distribution of trees, this resulted in samples. we selected this length of markov-chain because preliminary analyses suggested that stationarity was usually achieved by generations. we compared the posterior distribution from sampling multiple trees to that from the maximum clade credibility tree and examined the effect of changing the prior on the yule rate change parameter. we compared the posterior distribution to the point estimate produced by the likelihood version of the model, and to the akaike weights of each threshold point.  simulation tests we first tested the influence of tree depth on model performance. when deeper trees are simulated, coalescent and yule branching processes are expected to occur on more distinct time scales, and thus in general the model should perform better. the influence of tree depth is actually confounded by two issues, however. first, as the tree depth becomes shallower the implied rate of speciation increases because all trees contain species. if the rate of speciation approaches the rate of coalescence within species, then a sharp transition between processes should not be detectable. second, as the implied rate of speciation increases, more species originate relatively recently. the expected time to coalescence for a diploid, panmictic population is n generations. cladogenic events occurring more recently than this are expected to be increasingly difficult to delimit for two reasons: they are more likely to yield species that are not monophyletic and thus impossible to accurately identify under this model, and the most recent common ancestor of the daughter species is more likely to occur more recently than the threshold point. assuming species monophyly, the expected time to the mrca for two species that diverged n generations ago is n generations. therefore all probability should be on thresholds older than n generations, and most on thresholds older than n generations. again, when considering maternally inherited, haploid, organellar dna, equivalent times in n generations are halved, and n becomes the effective number of females in the population. this would give an expected time to mrca of n generations. the results of this test are dramatic. there is a clear increase in accuracy as well as a decrease in the range of the 95% highest posterior density interval with increasing tree depth. at shorter tree depths, the model’s performance diverged from expectations. when trees are short, larger numbers of species have divergence times younger than n generations and thus should not be detectable under the model. therefore we expected the number of species delimited to be smaller than for deeper trees. we did not observe this. instead, replicate posterior means became widely scattered around the true value of species. 95% hpds also did not uniformly increase with decreasing tree depth, instead, a wide distribution of hpds were observed. we also expected threshold times to have a minimum value of approximately n generations. at deeper tree depths this is observed, with % and % of mcmc steps sampling thresholds younger than n generations for tree depths of n and n, respectively. however, at n and n tree depths, 20% and 40% of mcmc steps sampled thresholds younger than n generations. these results indicate that the model performs well under demographic or sampling conditions that result in coalescent and yule processes occurring on very different time scales. it does not, however, perform optimally when those conditions are not met. ideally one would hope that as inference of the threshold point became more difficult, that the 95% hpds would increase, but still encompass the true value 95% of the time. this is not the case at the n and n tree depths. hpds generally become broader, but for increasing numbers of simulation replicates, they fail to encompass the true value. species arising in n generations constitutes a very rapid radiation, with an average of 89% of branches in the species tree shorter than the expected population coalescence time of n generations. failure to accurately assess credibility intervals in this case is likely because in this area of parameter space, the gmyc is no longer an accurate approximation of the real branching process in the gene tree. rather than there being a threshold between coalescent and speciation branching processes, the two processes are intermixed because there is little time for the independent evolution of lineages prior to speciation. note that these conditions will cause any dna barcode-based method of species discovery to fail and will also challenge more realistic models utilizing multilocus data and prior information on population assignment. next we examined the effect of intraspecific sampling. because the data points used by the model are waiting times between branching events, we expected that with species, we would not need extremely high sampling to accurately characterize the model, and that the distribution of samples among species would not be particularly important. our expectations were met. we found that sampling of individuals per species yielded poor results. replicate posterior means showed a strong bias towards inference of a large number of species. sampling of greater than individuals per species provided an improvement in the accuracy of the posterior means, but no change in the 95% hpd range. all sampling schemes greater than individuals per species appeared to yield similar results, including the more realistic condition of a lognormal distribution of alleles among species in which a large number of species are represented by singletons. while delimiting rare species, particularly from single specimens, is a challenge faced by taxonomists, our results suggest that the gmyc method may efficiently delimit species represented by singletons by calibrating the divergence threshold using data from better represented species. finally, we tested the effects of nucleotide sampling and the incorporation of phylogenetic uncertainty. we expected to find wider hpds with less sequence data, as uncertainty, particularly in branch lengths should be greater. we found a mild reduction in accuracy of the posterior means with up to bp of sequence, but after that, posterior means converged on those of the true tree. the 95% hpds improved with the addition of more sequence, but had not quite converged on those estimated from the true tree, even at bp. this suggests that uncertainty in phylogenetic estimation, particularly in typical dna barcode datasets of ~ bp will contribute substantially to uncertainty in species delimitation. three factors that could influence the accuracy of the model that were not explored here: migration, population substructure and selection. papadopoulou et al. examined the effects of migration on the formation of detectable gmyc clusters. they simulated datasets under an island model and found that even very low levels of migration caused likelihood ratio tests to fail to reject the null model of a single branching process. they interpreted this as evidence that the likelihood implementation of the model is conservative and will not infer species at all unless they are strongly isolated. papadopoulou et al.’s simulations assumed complete demic sampling, but lohse conducted simulations showing that under moderate migration rates and with a large proportion of demes unsampled that spurious, significant clusters could be inferred from the true gene genealogies. in his simulations, lohse showed that when demes were sampled from a metapopulation consisting of demes, that an average of species were inferred, and 80% of replicates rejected the null model. wakeley described the genealogical pattern resulting from such a process as having two phases that occur on very different time scales: a scattering phase, in which there is rapid coalescence and migration in local demes, and a collecting phase that begins when each remaining lineage is in its own deme and takes a very long time. in this case the gmyc might see the scattering phase as the “coalescent process” and the collecting phase as the “yule” process. further exploration of this issue is likely to be important, particularly if the gmyc is applied to phylogenetic samples with deep phylogeographic sampling. while lohse shows convincingly that this interaction of parameter space with sampling can mislead the gmyc, it is not clear to what extent these problematic areas of parameter space exist in real datasets. we simulated genealogies using ms under the conditions above and observed that the average time to coalescence of all lineages was n generations, with the scattering phase taking the first 4- n generations. if we assumed that these demes were species level taxa, each with θ = , we would expect to observe gmyc clusters with mrcas at a depth of - substitutions per site and the mrca for all lineages at substitutions per site. it is unlikely that the collecting phase would have the time to play out under this scenario, as it would take nearly million years at a mutation rate of substitutions per site per million years. if, by contrast, we assume that these demes represent populations at a smaller scale, each with a theta of / then we would expect mrcas of delimited clusters to be at a depth of - substitutions per site. with a typical dna barcode or short mitochondrial dna set of 650- bp, the scattering phase would be undetectable. the mrca for all lineages would occur at substitutions per site. unless this process was considered in the context of a larger species tree, it is unlikely that the gmyc would identify a significant branching threshold. empirical analyses we reanalyzed the empirical data used by pons et al. to illustrate the original formulation of the gmyc so as to provide a direct comparison of the implementations using representative data. the beast run converged after million generations and we discarded million trees as burn-in. the estimate of the standard deviation of the lognormal distribution of rates did not overlap so we could not use a strict clock with these data. when using samples of trees from the beast posterior distribution, the mean number of species estimated by the bayesian gmyc was and the 95% hpd ranged from to the rate change parameter for the yule process ranged as high as. in this model, the fold change in speciation rate from the root to the last speciation event is equal to np/n where n is the estimated number of species and p is the yule rate change parameter. in this case, given species, p =  allows for a 30-fold speciation rate increase. we thought this might be unrealistically high, and thus re-ran the analysis with the prior distribution set to u, or a maximum 2-fold increase. this minimally influenced the results. we also analyzed the maximum clade credibility tree under both priors, and using the likelihood implementation of the gmyc. we compared the results of likelihood and bayesian analyses by calculating akaike weights for each possible threshold in the likelihood analysis. akaike weights are the relative likelihoods of a set of models and thus suited for qualitatively comparing results among these analyses. the bayesian gmyc analysis of the maximum clade credibility tree yielded a mean of species and a narrower 95% hpd of 38– species. when applied to the maximum clade credibility tree, the u prior did change the results, yielding a posterior mean of species and a 95% hpd of to the maximum likelihood analysis resulted in species with a 95% confidence interval of to species. the results of the likelihood and bayesian analyses of the maximum clade credibility tree are very similar, particularly when the prior u is placed on the yule parameter, but when sampling trees, the posterior distribution is substantially broader. at least some of this uncertainty stems from variation in topology. plots of pairwise probabilities of conspecificity demonstrate this. the probability distribution of the number of species in the sample is also wider for the posterior distribution than for the akaike weights. our bayesian approach is similar in spirit to the model-averaging approach of powell in that its goal is to make the inference of species limits and analyses based on them more robust to uncertainty. there are, however, three major differences. first, we take into account phylogenetic uncertainty. as indicated by our results, this is perhaps the most influential difference, although in theory it could also be accounted for using model averaging. second, our bayesian method requires the specification of prior knowledge. depending on a researcher’s comfort with assigning prior probability distributions to the rate and threshold parameters, this is either an advantage or a disadvantage. finally, the treatment of nuisance parameters is fundamentally different. while in the model-averaging approach, inferences are conditioned on maximum-likelihood point estimates of nuisance parameters, the bayesian approach integrates out nuisance parameters, giving marginal probabilities of species limits. these final two differences are intrinsic to bayesian inference, and researchers choosing among methods will need to consider their choice of statistical paradigm. we note that our confidence intervals in all analyses are far wider than those of pons et al., which were 46– species. this is most likely because of the difference in obtaining ultrametric trees. pons et al. used non-parametric rate smoothing on a maximum likelihood estimate of the tree. while they achieved sensible results with narrower confidence limits than ours, we nevertheless advocate an approach that samples trees and fits them to a clock model using a parametric method. this allows for a full accounting of uncertainty associated with phylogenetic estimation, albeit at the cost of some precision.  our results demonstrate that the bayesian implementation of the gmyc model is reasonably reliable given two caveats. first, the length of the dna sequence is important. we found that when we sampled only bp, or only alleles per species, that the performance of the model declined strongly. second, the model is only useful when the underlying history of the species under consideration lies in particular regions of parameter space. species that have recently diverged, or clades undergoing rapid radiation are unlikely to be identifiable under the model. in the latter case, the model may provide misleading estimates and confidence. cases such as these, however, may be recognizable because the results may be highly unexpected in the context of other sources of data such as morphology or geography. our implementation of the model provides two main improvements over the original. first, it allows the specification of prior probabilities on model parameters. it is our experience that very high values of the yule process rate change parameter sometimes have high likelihood and result in high uncertainty in the threshold parameter. these high values may be biologically unrealistic, and the specification of an informative prior can reduce the posterior probability of those areas and produce a more accurate estimate of diversity. second, it allows for the characterization of species limits without use of a point estimate of the phylogeny. we know that many datasets are associated with substantial uncertainty owing to limited sequence data collection. the bayesian gmyc method provides marginal probabilities of species identities and will allow downstream estimates of species diversity and community structure to account for uncertainty underlying species designations. an important future direction for this work is to implement the multiple-threshold version of the model proposed by monaghan et al., which can account for greater variation in divergence times and effective population sizes than the model implemented here. it has been shown to provide a better fit to some datasets, but will require implementation of a more complex reversible-jump mcmc that allows proposals that change the number of parameters in the model. it is widely acknowledged that single-locus data are not optimal for the inference of phylogeny, historical demography, or species limits. nevertheless, vast amounts of biological diversity remain undescribed at the level of species, and this limits our ability to understand the evolutionary history of our planet and its current ecological functioning. available alternative means of describing species diversity, either from molecular or morphological data have major drawbacks in that they are time consuming, expensive, or subject to their own biases and inaccuracies. single-locus data for many groups are currently being generated on a large scale, and we advocate making the best of this data. we believe that under certain conditions, the gmyc model can be useful, and that a bayesian framework accounting for uncertainty is most appropriate for these data. competing interests the authors declare that they have no competing interests. authors’ contributions nr and bc developed the original concept and designed the simulation experiments. nr wrote all software and conducted the analyses described in the manuscript. nr and bc interpreted results and wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 figures s s s these figures display the distribution of mcmc samples for each treatment and each replicate within treatments for simulated data. s is results from the tree depth simulation, s is the results from the allele sampling simulation and s is the results from the nucleotide sampling simulation. click here for file acknowledgements we thank the national science foundation for funding aspects of this work. we thank jeremy brown for conversations that initiated this research, and members of the carstens lab at louisiana state university for discussions related to and comments on this manuscript. we thank dr. timothy barraclough and two anonymous reviewers for helpful correspondence regarding this work and comments on drafts of the manuscript.
here we present a bayesian implementation of an evolutionary model-based method, the general mixed yule-coalescent model. we also re-analyze the dataset used to introduce the original gmyc model. our results suggest that the gmyc model can be useful under a wide variety of circumstances, particularly in cases where divergences are deeper, or taxon sampling is incomplete, as in many studies of ecological communities, but that, in accordance with expectations from coalescent theory, rapid, recent radiations may yield inaccurate results. we conducted simulations testing the effects of species evolutionary history, levels of intraspecific sampling and number of nucleotides sequenced. we found that the model results are improved with addition of dna sequence and increased sampling, although these improvements have limits. recent and rapid divergences result in higher amounts of uncertainty in the model and eventually cause the model to fail to accurately assess uncertainty in species limits. in such cases dna sequence-based species delimitation has been proposed as a means of estimating species boundaries for further analysis. our implementation integrates over the parameters of the model and uncertainty in phylogenetic relationships using the output of widely available phylogenetic models and markov-chain monte carlo simulation in order to produce marginal probabilities of species identities.
 species are considered the fundamental unit in many ecological and evolutionary analyses, yet accurate, complete, accessible taxonomic frameworks with which to identify them are often unavailable to researchers. in such cases dna sequence-based species delimitation has been proposed as a means of estimating species boundaries for further analysis. several methods have been proposed to accomplish this. here we present a bayesian implementation of an evolutionary model-based method, the general mixed yule-coalescent model. our implementation integrates over the parameters of the model and uncertainty in phylogenetic relationships using the output of widely available phylogenetic models and markov-chain monte carlo simulation in order to produce marginal probabilities of species identities.  we conducted simulations testing the effects of species evolutionary history, levels of intraspecific sampling and number of nucleotides sequenced. we also re-analyze the dataset used to introduce the original gmyc model. we found that the model results are improved with addition of dna sequence and increased sampling, although these improvements have limits. the most important factor in the success of the model is the underlying phylogenetic history of the species under consideration. recent and rapid divergences result in higher amounts of uncertainty in the model and eventually cause the model to fail to accurately assess uncertainty in species limits.  our results suggest that the gmyc model can be useful under a wide variety of circumstances, particularly in cases where divergences are deeper, or taxon sampling is incomplete, as in many studies of ecological communities, but that, in accordance with expectations from coalescent theory, rapid, recent radiations may yield inaccurate results. our implementation differs from existing ones in two ways: it allows for the accounting for important sources of uncertainty in the model and in the specification of informative prior distributions that can increase the precision of the model. we have incorporated this model into a user-friendly r package available on the authors’ websites. species delimitationgmycbayesian phylogeneticsdna barcoding 