ASAP: assemble species by automatic partitioning

Here, we describe Assemble Species by Automatic Partitioning (ASAP), a new method to build species partitions from single locus sequence alignments (i.e., barcode data sets). ASAP is efficient enough to split data sets as large 104 sequences into putative species in several minutes. Although grounded in evolutionary theory, ASAP is the implementation of a hierarchical clustering algorithm that only uses pairwise genetic distances, avoiding the computational burden of phylogenetic reconstruction. Importantly, ASAP proposes species partitions ranked by a new scoring system that uses no biological prior insight of intraspecific diversity. ASAP is a stand‐alone program that can be used either through a graphical web‐interface or that can be downloaded and compiled for local usage. We have assessed its power along with three others programs (ABGD, PTP and GMYC) on 10 real COI barcode data sets representing various degrees of challenge (from small and easy cases to large and complicated data sets). We also used Monte‐Carlo simulations of a multispecies coalescent framework to assess the strengths and weaknesses of ASAP and the other programs. Through these analyses, we demonstrate that ASAP has the potential to become a major tool for taxonomists as it proposes rapidly in a full graphical exploratory interface relevant species hypothesis as a first step of the integrative taxonomy process.

are widely used, as they are easy to apply on DNA-barcoding data sets, even large, and precisely because they do not necessitate predefined species hypotheses. The most popular are general mixed Yule-Coalescent model (GMYC) (Pons et al., 2006), Poisson tree process (PTP) (Zhang et al., 2013), both first developed in a maximum likelihood framework, and later extended to a Bayesian framework (Reid & Carstens, 2012), and automatic barcode gap discovery (ABGD) . GMYC and PTP take a phylogenetic tree as input and estimate rates of branching events to infer which part of the tree more likely follows a speciation model (the deepest part) and which part follows a coalescent model (subtrees of the shallowest part). The species partition is found by maximizing the likelihood of the transition between these two branching rates, GMYC in absolute time (hence the need for an ultrametric tree), PTP in mutational time at different nodes of the tree. GMYC first inferred a single transition event between the two rates (speciation versus coalescent); PTP first had two rates (speciation and coalescent).
Both were later expanded to infer "multiple thresholds", allowing several transitions to occur in different subtrees (Kapli et al., 2017;Monaghan et al., 2009).
Contrary to the two previous methods, ABGD uses only pairwise genetic distances (no tree is inferred) and automatically identifies in their distribution the so-called "barcode gap". This gap marks the limit between the smaller intraspecific distances and the larger interspecific distances. From the gap, a distance threshold is estimated and used to partition the samples into putative species. A coalescent model is used to identify the position of the most likely barcode gap, based on a maximal genetic intraspecific divergence P defined a priori by the user. Consequently, users must provide a range of P in which ABGD identifies one or several barcode gaps and the method outputs the corresponding species partitions. For a single data set, ABGD thus eventually proposes several partitions that correspond to different prior values P. In its recursive version, ABGD is applied on each group of the initial partition, and eventually splits them when internal barcode gaps are detected.
Although the results obtained with the various methods often vary depending on data set characteristics (e.g., Blair & Bryson, 2017), the main conclusions of these studies are: (a) All methods generally perform well (but see e.g., Dellicour & Flot, 2018) being mostly congruent (i.e., providing similar species partitions) with each other and with the species partitions inferred from independent data (e.g., other molecular markers, morphological data, ecological data); (b) All of them perform poorly when the number of sampled individuals per species is too low (Ahrens et al., 2016), or when the contrast of intra-versus interspecific divergences is mild. This contrast varies with species ages, mutation rates, population sizes, strengths of the selection and degrees of within-species population structure (Pante et al., 2015;Pentinsaari et al., 2017;Ritchie et al., 2016); mPTP was in particular developed to overcome this issue (Kapli et al., 2017); (c) Partitions proposed by the three methods sometimes differ, each of them being able to infer the "correct" species when the two others fail. This led some authors to propose that all three methods (among with eventually others) should be applied jointly and compared (Ducasse et al., 2020); and (d) Although there are several exceptions (e.g., Blair & Bryson, 2017), ABGD in particular, and PTP to a lesser extent, tend to lump species more than GMYC (Pentinsaari et al., 2017). Conversely, the multiple-threshold version of GMYC is particularly prone to oversplitting (Fujisawa & Barraclough, 2013;Kekkonen & Hebert, 2014).
In comparison with GMYC and PTP, ABGD has the advantage of being very fast, mainly because it bypasses the phylogenetic reconstruction. Furthermore, because ABGD identifies a species partition for each value of P defined a priori, several partitions may be proposed, reflecting the uncertainty stemming from the data and encouraging the user to evaluate the relevance of the ABGD partitions in the light of other data, as it is recommended in an "integrative taxonomy" approach. However, ABGD does not provide a score for each partition that would help the user to identify the "best" partition(s), and this probably constitutes the main drawback of ABGD (judging from the numerous comments and questions the authors of ABGD have received from the users).
In this article, we describe a new method of species delimitation, still based on pairwise genetic distances, but which implementation provides a score for each defined partition and overcomes the challenge of a priori defining P. Our new algorithm, ASAP, still provides several partitions, more or less fine-grained, but ranked using a new scoring system. Importantly, we also develop a full graphical web-interface to ease its usage. However, ASAP, like any other method, must not replace the taxonomist work, as any partition of species must be subsequently tested against other evidences in an integrative taxonomy framework. This is especially crucial as ASAP uses single-locus data that are known to bear weaknesses.

| Overview of the ASAP software
ASAP is a self-contained program written in C. Users can use ASAP either through a full graphical web-interface (https://bioin fo.mnhn. fr/abi/publi c/asap), or download and compile the sources for local usage (same url).
Our algorithm is an ascending hierarchical clustering, merging sequences into "groups" that are successively further merged until all sequences form a single group. At each merging step, the assignment of all sequences into groups is named a "partition". The first partition contains as many groups as sequences (no grouping was yet done) whereas the last partition is a single group with all sequences inside. Larger groups are created by merging groups of the previous partition together. We characterize all newly created partition in two complementary ways. First, we assign to it a probability that quantifies the chances that each of its new groups is a single species.
Second, we compute the width of the barcode gap between the previous and this new partition. Both metrics (probability and barcode gap width) are combined into a single asap-score that is used to rank the partitions.

| Ranked distances
We first start by computing, when not provided, all pairwise distances between the n sequences of the alignment. Distances are then ranked by increasing values. The efficiency of the algorithm stems from the fact that each distance is only considered once in increasing order for clustering purposes.

| Hierarchical clustering
The clustering process starts with a first partition where each sequence belongs to a different group. ASAP then treats each of the ranked distances one by one in increasing order (equal distances are treated together) as a threshold value for delimiting groups: sequences separated by a distance equal to the current value d C are clustered into the same group. Consequently, when sequences that were in different groups are clustered together, the previous groups are merged into a new larger group, and is associated to the current clustering distance, d C . Importantly, a new partition can have a single new group or several new ones when several sequences from different groups are merged independently into different groups for the same distance d C . When a new partition is built, the clustering process pauses. ASAP then scores all new groups with a probability of panmixia. It also scores the new partition using an ad hoc score computed from both the barcode gap width and probabilities of panmixia. After the group(s) and partition scoring, ASAP then continues the clustering by looking after the next distances until another partition is built. The algorithm stops when all sequences are merged into a single final group.

| Computing p-values For each group
We aim at computing a p-value for a newly created group that is a merge of two or more subgroups. We compute Π intra the average pairwise distance between sequences within the subgroups and Π inter the average pairwise distance among sequences of different subgroups ( Figure 1). We then compare Π intra to its theoretical distribution, computed by Monte-Carlo simulations of a neutral coalescent model assuming a single panmictic species with a sample size m and a coalescent mutation rate θ = Π inter /[2 × (1 − 1/m)]. The value of θ is set so that in the simulations the distance between sequences connected by the most recent common ancestor (MRCA) of the group (π inter ) is equal, on average, to the observed one: E[π inter ] = Π inter . This relates to the average time to the MRCA that is 2 × (1 − 1/m), expressed in coalescent time (Wakeley, 2009). We compute the p-value as the fraction of replicates where the simulated π intra is equal or lower than the observed Π intra . The number of replicates is updated on the fly to have correct estimations of low p-values. Put differently, it quantifies under H0 (one single species) the probability of observing a diversity Π intra or less within the subgroups given that the divergence between the subgroups is on average Π inter .

For partitions
We compute the probability to observe π intra or less diversity within all subgroups of the current partition (that are groups of previous partition before the merge) assuming that all new groups of the current partition are independent coalescent realizations with θ estimated for each group independently.

| Recursive splits
Once a new partition is built, ASAP tests for each of the groups of the partition whether its p-value is lower than a given risk (by default 1%) and consequently should be split. When a group is split, ASAP recursively descends to all its subgroups and assesses whether they should be split as well.

| Relative barcode gap width
ASAP also computes a relative barcode gap width associated to the current partition (Appendix S1). The partition is associated to a threshold distance d T that is the mid-point between the current distance, d C (with rank r C ), that triggered the merging and the previous distance in the list d C-1 (with rank r C − 1). A barcode gap corresponds to a "jump" in the distance values in only few ranks.
While increasing only few ranks in the list, the distance will "jump" from a value that is (much) less than d T to a value that is (much) higher than d T . To quantify the barcode gap width, ASAP scans downward the distance list from d C-1 until it finds the first distance smaller than 0.9d C-1 : this is d L which rank is r L in the list. It then scans from d C the distance list upward until it finds the first distance above 1.1d C : this is d H which rank is r H . The relative gap width W is defined as: to compute the "relative" width of the gap; the "+1" only prevents the ratio to be very high when distance values are very small. The higher the W, the larger the barcode gap.

| Outputs
At the end of the clustering, ASAP scores and sorts all the different partitions using two criteria: their p-value sorted (see Section 2.2.3.2) by increasing order (the smallest p-value has rank 1) and their rank of relative barcode gap width (see Section 2.2.5) sorted by decreasing order (the largest gap has rank 1). The asap-score is the average of both ranks: the smaller, the better. Furthermore, ASAP produces a graphical output where each node of the hierarchical clustering is colour-coded depending on its probability of being a panmictic species (see Section 2.2.3.1). Thus, the colour guides the user finding which nodes may be split into smaller groups. Several other graphical options are provided to help the user navigate among partitions and choose the "most relevant" partition, beyond a simple naive use of the asap-score (Appendix S2).

| Tests on empirical data
To compare the results obtained by four methods (ASAP, (m)PTP, (m) GMYC and ABGD), we selected 10 empirical COI data sets covering various taxa (birds, mammals, amphibians, insects, crustaceans and molluscs) and including 44 to 2,574 specimens that belong to 5 to 643 species (Table 1) (Borisenko et al., 2008; Elias-Gutierrez F I G U R E 1 An illustration of the clustering algorithm on a small data set of nine sequences. On the lower part, we report how ASAP proceeds (downward in the figure) through the list of ranked distances (on the left), merging successively sequences into groups (highlighted in coloured blocks). For each new group, ASAP computes a p-value that this new group is a panmictic species (values reported on the right part) based on pairwise differences within (intra) and between (inter) subgroups. Furthermore, each time a new group is created, a new partition is built (a sequence of blocks in the central part) that is associated to the current distance d C . The distances d C at which the partitions are instantiated are represented in a phenetic tree (top part). Each node is a group, each horizontal dashed line is a partition. For each newly created partition, ASAP also computes a probability of panmixia (p-val) and a relative gap width metrics (W). Then using their respective ranks (given in parenthesis), ASAP computes an ad hoc ASAP-score: the lower the score, the better the partition [Colour figure can be viewed at wileyonlinelibrary.com] S1 S3 S4 S2 S5 S7 S9 S8 S6 1 S1,S3 0.0029 2 S1,S4 0.0029 3 S3,S4 0.0029 4 S1,S2 0.0058 5 S2,S3 0.0058 6 S2,S4 0.0058 7 S7,S9 0.0231 8 S7,S8 0.0289 9 S1,S5 0.  (1)  For all empirical data sets, we used the web version of ABGD, with default parameters. Only the initial partitions were considered, and only the more stable partition(s) (i.e., the partition(s) found with several P in the vicinity of the barcode gap) was (were) reported.
For ASAP, we used a recursive split probability of 0.01 (see Section 2.2.4), and report (a) the partition with the best asap-score as well as (b) the partition that is closest to the "correct" one among the two best partitions, according to their asap-scores. For GMYC and mGMYC, ultrametric trees were reconstructed using BEAST 2 , with an independent GTR substitution model for each codon position. Relative divergence times were estimated using a relaxed log-normal clock with a coalescent prior and a constant population size, following the recommendations of Monaghan et al. (2009)

| Simulations
We measured the power of ABGD, GMYC, (m)PTP and ASAP to retrieve the correct species partition in various scenarios using Monte Carlo simulations. We used a "multispecies coalescent" framework (Rannala & Yang, 2003) with different options and parameters using Monte-Carlo simulations, as described previously . Note that contrarily to the standard multispecies coalescent, the species tree is here drawn from a probability distribution. The home-made C simulator is available upon request.
Briefly, for each replicate, we generate a species tree using either Note: Each line represents a data set which numbers of sequences (#seq) and species (#spec) are reported in the provided reference. We compare the "true" number of species to the predictions made by the partition ranked first by ASAP (ASAP 1st), by the "best" partition among the two first predicted by ASAP (ASAP 1st-2nd), the "best" partition by ABGD and the unique partition predicted by PTP, mPTP, GMYC and mGMYC. There is no partition for Birds by GMYC and mGMYC as we were not able to obtain a Bayesian tree given the large number of sequences. Cells were coloured in dark grey when predictions were very accurate (at most 5% different from the referenced number of species) and with light grey when accurate (between 5% and 10%). model (all species arose at the same time). Radiation (hard polytomy) models cases where all speciation events follow each other quickly and where no mutations have occurred between the first (the root) and the last speciation event. We used a backward coalescent version of these models that we have previously used for ABGD evaluation . For the radiation model a unique speciation event, exponentially distributed with rate r, is drawn. For the Yule model (n sp −1) speciation events are drawn with identical rate (Lambert & Stadler, 2013).
Once the species tree is obtained, we assign sequences to species uniformly, with at least one sequence per species. All species (current and ancestral) are assumed to be of equal effective size (i.e., PTP, we used as input the "true" gene genealogy (the one simulated for the replicates) not only to fasten the simulation (i.e., skipping the phylogenetic reconstruction) but also to assess their power when the phylogeny is perfectly reconstructed. We would like to emphasize that only ASAP used unprocessed data (polymorphic sites) without any biological insights (no prior, no phylogeny reconstruction nor calibration).

| Empirical data sets
We first assessed the ability of ASAP through a proxy that is its ability to retrieve the "correct" number of species in 10 empirical data sets (Table 1). The data sets were selected to represent test cases of different sizes (from 44 sequences/5 species to 2,574 sequences/643 species). We first report the number of species predicted in the partition with the best asap-score (ASAP first): we found that in 4/10 of the data sets, the partition with the best asap-score is "very close" to the reference one (less than 5% difference in terms of species numbers) and that 8/10 is "close" (<10% difference). If we also consider the partition with the second best asap-score (ASAP first and second), the degree of accuracy increases to 6/10 for the very close ones and 9/10 for the close ones. This is a good indication that ASAP users should consider not only the partition with the best asap-score but also few subsequent ones. It is important to report that here no extra biological knowledge was considered for ASAP predictions.
One could for example use threshold distances (e.g., d T or d C ) to prefer one partition over another despite a poorer asap-score (e.g., in most clades intraspecific diversity is typically on the order of 1%, not on the order of 10%). Obviously, other criteria and characters should also be used to choose a final species partition, in an integrative taxonomy context.
One of the ASAP main qualities is that it is extremely fast compared to any method that relies on tree reconstruction. The online version takes 45 s for the largest data set of Table 1  On a curated unpublished moth data set, it took 6 min 35 on the website to delimit 2,466 species (best asap-score) or 2,067 (second best asap-score) from 9,396 sequences. Subsequent partitions with lower asap-scores are close to one or the other of these two first partitions. Because of its rapidity, ASAP web server accepts up to 10 4 sequences (unlike the ABGD server).
We also took the opportunity of analysing the 10 data sets to assess the performance of other methods: ABGD which is solely based on pairwise distances, PTP and mPTP that were run on an ML trees (i.e., RaxML) and GMYC and mGMYC on an ultrametric trees estimated by a Bayesian MCMC method (i.e., BEAST). Results (Table 1) show that ABGD performance is similar to ASAP first-second, that PTP and mPTP tend to not perform very well, that GMYC performs very well provided that the number of species is not too large and that, as previously reported in the literature, mGMYC generally oversplits (Fujisawa & Barraclough, 2013;Kekkonen & Hebert, 2014).
Note that ABGD performances are somehow overestimated as we report the partition that is the closest to the reference one over the whole range of P. We could not use GMYC for the largest data set as the Bayesian tree reconstruction did not converge after several weeks of computation.

| Simulated data sets
We then assess the theoretical performance of ASAP using Monte-Carlo simulations of a multispecies coalescent framework. In brief, a random species tree is generated using either a Radiation model, where all species arose in single event, or a Yule model, where the speciation events occur at constant rate independently in all branches. In both model, we tune the separation of time scales (speciation versus intraspecific coalescent events) using a speciation rate that is expressed in coalescent time (i.e., N generations per unit of time). The lower the speciation rate, the better the separation of time scales. For example, when the speciation rate is 0.1, speciation events are 10 times slower than pairwise coalescent events within species.

| The impact of speciation rate on ASAP
We first examine the ability of ASAP to correctly retrieve four species in both speciation models as a function of the speciation rate (from 0.001 to 1). We report in Figure 3 the fraction of runs where ASAP was able to correctly retrieve the four species (top panel) and the average number of predicted species, regardless of their composition (bottom panel). We assess the quality of the partition with the best asap-score (ASAP first) as well as the quality of the partition that is the closest to the truth among the two best partitions (ASAP first-second).
We observe that for low rates of speciation, the best partition proposed by ASAP correspond exactly to the four species. This is an "easy" case where the two timescales are well separated. As the speciation rate increases, both time scales overlap and it becomes harder to delineate species using pairwise genetic differences at a single locus. When the speciation rate is larger than 1, speciation events are more recent than intraspecific divergence so that individuals within species are no more different than individuals between species.
ASAP performs usually better with the Radiation than with the Yule model. This is especially striking for moderate speciation rate F I G U R E 3 Performance of ASAP as a function of the speciation rate. For two alternative models of speciation (Radiation and Yule), we report the fraction of replicates where ASAP find the four correct species (top panels). We considered either only the partition with the best asap-score (ASAP-1) or the partitions ranked first and second (ASAP-1/2). Obviously, the latter has better performance. We also report the average number of predicted species, regardless they are correct or not (bottom panels). Each point is evaluated on 500 replicates (e.g., 0.03). For radiations, most of the errors correspond to oversplit, as illustrated by the average number of predicted species that is larger than four. Under the Yule model with four species, there are three independent speciation events and consequently there is a higher chance to generate at least one very recent speciation event that would be invisible in regard of sequence divergence. Indeed, the most recent event is exponentially distributed with rate 3r. As a consequence, contrarily to the radiation model, ASAP failures correspond for this rate to cases where it lumps the two closest species into a single one.

| The impact of the number of species on ASAP
Second, we explore the impact of the number of species for a fixed sample size of 200 sequences, with r = 0.01, a moderately challenging speciation rate. We report the average number of predicted species regardless of their composition for both the radiation and the Yule models. Results (Figure 4) show (

| The impact of the number of species on ABGD, PTP and GMYC
We apply the same analysis to ABGD, (m)PTP and GMYC. We would like to emphasize again that we assessed their power under optimal conditions: a single "excellent" prior for ABGD representing a perfect knowledge of intraspecific diversity and the "true" simulated tree for (m)PTP and GMYC, bypassing their main limitations, that is having a correctly reconstructed phylogenetic tree. As a consequence, we here overestimate their power for realistic biological situations where only a set of sequences is available (neither the true tree nor prior knowledge of intraspecific diversity is known). ASAP, on the contrary, directly uses the sequences and needs no prior biological insight or phylogenetic reconstruction.
The power assessments of the methods (Figure 4) show that ABGD retrieves well the correct partition when speciation occur as a single radiation but has a limited power when speciations follow a Yule model. On the contrary, we found that GMYC performs very well for the Yule model but is less efficient for a radiation model.
Interestingly mPTP consistently split a constant small number of species. It thus performs poorly when the number of species is low but quite well when the number of species is 50 or more.

| D ISCUSS I ON
We introduced a new species delimitation program, ASAP, fully exploratory, in the sense that it does not require any a priori knowledge, neither on the number of species, the species composition, or any biological information, such as a phylogenetic tree or a prioridefined intraspecific genetic distances. Only pairwise genetic distances are used to build a list of partitions ranked by a score. This composite score is computed using the probabilities of groups to be panmictic species and the barcode gap widths. ASAP overcomes the two mains limitations of ABGD, namely (a) the need for an a priori defined P; and (b) the lack of a scoring system.
However, and contrary to some other methods, ASAP still outputs several partitions, ranked by their asap-scores. A list of the "best" partitions (10 by default) is provided in the output together with their gap-width score, their p-value, their threshold distance d T and the number of species they correspond to.
The graphical output of ASAP has four main components (Appendix S2): 1. A list of partitions ranked by their asap-score that putatively correspond to species hypothesis, 2. A plot of the asap-score as a function of d C . We report the asapscore of all partitions (not only the best ones) as a function of the clustering distance d C to appreciate whether all good partitions have similar d C or whether "potentially good" partitions can drastically differ in size.
3. An ultrametric clustering tree of all sequences, where the distance to the leaves lengths correspond to the distance d C at which these sequences were clustered in the same group. All nodes of this tree are colour-coded depending on their p-value (the darker the more it differs from a panmictic species). 4. A "boxed-species" graph, where species hypotheses in the different partitions are represented as vertical boxes in front of the ultrametric tree.
When a partition is selected by a click in any of the three panels, it is automatically highlighted in the two other components.
We also propose a complementary representation, where we display the hierarchical tree with, at its leaves, the 10 best ASAP partitions where their groups are depicted as boxes (that are similar to the boxes of Figure 1).
We have evaluated ASAP strengths and weaknesses using both real and simulated data. Our benchmark shows that ASAP performs well delivering partitions in a matter of minutes even for data sets as large as 10 4 sequences. ASAP is thus meant to be applied on large single-locus data sets when no species hypothesis is available, as typically produced in DNA-barcoding projects. Although the web version limits the input to 10 4 sequences, more sequences can be analysed using a local command-line version of ASAP (sources are available on the webserver). and applying several methods to a given data set is a strategy commonly applied to maximize the probability to detect species complexes, identified as groups of species whose limits vary depending on the method.
Importantly, several other methods can also be used to delimit species, such as BINs (Ratnasingham & Hebert, 2013), Jmotu (Jones et al., 2011) or VSEARCH (Rognes et al., 2016), among others (e.g., Rannala & Yang, 2020). We are also aware that the number of predicted species is only a proxy to assess the performance of the different methods. Indeed, other metrics such as the F-measure (Larsen & Aone, 1999) or the number of splits or merges (Ratnasingham & Hebert, 2013) give also insightful information. Some of them are even implemented in meta-analysis software such as LIMES (Ducasse et al., 2020), which could be used to perform a more extensive benchmark of all existing methods using a wider spectrum of metrics.
More generally, and as advocated by the proponents of the integrative approach in taxonomy, the use of a single marker with a single method of species delimitation should be avoided, precisely because each method has its own limitations. Some methods are based on a phenetic criteria (e.g., ASAP and ABGD) while others on phylogenetic criteria (e.g., (m)GMYC and (m)PTP). Furthermore a single locus may not follow the species history, because of introgression and incomplete lineage sorting. This is particularly true for species in the grey zone, in which the gene tree may differ from the species tree, and the coalescent events may be older than the speciation events (De Queiroz, 2005). For this reason, we recommend that single-locus methods are to be used as a first step of the species delimitation process that is to propose primary species hypotheses.
This is for example useful in groups for which there is no pre-existing hypotheses to test, or for which unknown/incorrectly delimited species represent the majority of the diversity (e.g., microbial communities or hyperdiverse groups of eukaryotes, such as insects, spiders, nematodes, molluscs). Furthermore, DNA barcodes are now routinely produced using NGS approaches, providing large numbers of sequences often not assignable to known and sequenced species (Kennedy et al., 2020), and for which methods such as ASAP are welcome to e.g., compare species diversity among sites.
In a second step it is then the responsibility of the taxonomist to evaluate with other methods (in particular, methods that will evaluate alternative partitions of species) and/or lines of evidence (such as other genetic markers, morphology or ecology) whether the proposed hypotheses are robust, or not. In this context, methods such as ASAP, ABGD, (m)PTP and (m)GMYC should thus be seen as a formalized and reproducible way to propose species hypotheses in groups where no such hypotheses exist, or, if they do exist, that are better to be ignored.

AUTH O R CO NTR I B UTI O N S
S.B., G.A., and N.P. designed the method; G.A. developed the algorithm and tested it on simulated data sets; S.B. wrote the program and created the web-interface; N.P. performed the tests on real data sets; G.A., and N.P. wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
ASAP is available at https://bioin fo.mnhn.fr/abi/publi c/asap. Data sharing is not applicable to this article as no new data were created or analysed in this study.
The software used to simulate multispecies coalescent with random speciation time was written in C and is available upon request, as well as the simulated data sets. All real data sets are directly accessible from the ASAP website.