Effects of phylogenetic reconstruction method on the robustness of species delimitation using single‐locus data

1Species are a fundamental unit for many fields of biology, yet their identification and delimitation are rarely straightforward (Hebert. ). Molecular techniques allow for rapid and broad assessment of diversity of poorly known groups or where traditional techniques are difficult (Blaxter; Tang.; Fontaneto ). Well-established metrics for species delimitation exist (see Sites & Marshall; Birky.; Flot, Couloux & Tillier; Puillandre. ) but only a few are grounded in evolutionary theory and do not require a priori hypotheses regarding species entities (O'Meara; Yang & Rannala ). Fewer still are designed for large-scale single-locus marker surveys (Fujisawa & Barraclough; Zhang. ). With the increased frequency of DNA taxonomy studies and their potential marriage with next generation sequencing technologies (NGS – Creer. ), there is a need to determine potential sources of bias on diversity estimates. Here, we evaluate robustness of the generalised mixed Yule coalescent model (GMYC) and the Poisson tree process (PTP) species delimitation methods to different approaches of phylogenetic reconstruction of the gene trees. Robustness was assessed by how topological and branch length variation introduced by phylogenetic methods influences delimitation estimates in terms of species richness and identity.

2A special branch of phylogenetic species delimitation (see Sites & Marshall ) is coalescent-based species delimitation methods (Pons.; Fontaneto.; Zhang. ), which combine coalescent theory with diversification models to infer the transition point between population and species-level processes on a gene tree. These approaches provide objective, clade-specific threshold(s) with which to delimit evolutionarily significant units (ESUs) of diversity (akin to species, as defined by the Evolutionary Species Concept – Simpson ). These methods provide an alternative to operational taxonomic unit (OTU) picking methods, which rely on arbitrary, clade-specific sequence similarity thresholds (Barraclough. ).

3The GMYC is one of the most popular coalescent-based species delimitation methods and is designed for single-locus data (Fujisawa & Barraclough; although it can be used with concatenated-loci data, Pons.; Fontaneto. ) and has been used to describe new species (Birky. ). The method separately models the fit of Yule (pure birth; Yule ) and coalescent processes (Hudson ) to an ultrametric tree to define the transition from species-level to population-level processes, used to delimit ESUs. The PTP (Zhang. ) is a recently developed method that models speciation and coalescent events relative to numbers of substitutions rather than time, and uses heuristic algorithms to identify the most likely classification of branches into population and species-level processes, used to delimit ESUs. This approach assumes either that substitutions are clocklike or, if substitution rates vary across the tree, that coalescent and speciation events occur at a constant rate per substitution event, rather than per unit of time. The key advantage of the PTP, however, is that it is devised for non-ultrametric trees.

4Several studies have evaluated factors that could bias accuracy of the GMYC and PTP. For the GMYC, simulation studies have addressed the effects of various aspects of sampling (Papadopoulou.; Bergsten.; Reid & Carstens; Talavera, Dinca & Vila ), population size and speciation rates (Esselstyn.; Fujisawa & Barraclough ). For the PTP, simulations have been used to evaluate the effect of birth rates (i.e. evolutionary distances between species) and sampling unevenness (Zhang. ). Less attention has been paid to the influence of different phylogenetic methods for reconstructing the underlying gene tree. For coalescent-based species delimitation, phylogenetic and branch smoothing (defined as methods that correct rate heterogeneity to make ultrametric, clocklike trees) methodology are potentially large sources of bias if branch length and topological variation is introduced by different phylogenetic methods, for example, by different treatment of unresolved nodes and rate heterogeneity. Zero-length branches introduce infinite (logarithmic) branching rate artefacts that might bias species delimitation and underestimate (early placement of the threshold) or overestimate (recent placement) species diversity (GMYC), and heterogeneity in the rate of molecular evolution among lineages would violate the assumption that branching events can be modelled against substitutions directly (PTP). It is well known that different methods of rate smoothing introduce variability in branch lengths (Drummond & Suchard ) that can ultimately affect inferences made from the tree (Rutschmann ); artificially variable branch lengths might therefore result in variable diversity estimates with the GMYC. A previous assessment of the effect of phylogenetic methods on GMYC ESU estimates showed that certain method combinations perform poorly (Talavera, Dinca & Vila ), but is not clear whether this is generally true.

5The GMYC, in combination with at least 11 different phylogenetic and 9 smoothing methods (Table S1), has been used in over 150 studies. BEAST (Drummond & Rambaut ) is the most popular software for obtaining gene trees (48·9%), followed by MrBayes (25% – Ronquist. ) and RAxML (8·3% – Stamatakis ). BEAST is also the most popular software for rate smoothing (53·3%), followed by r8s (28·5% – Sanderson ), PATHd8 (6·7% – Britton. ) and chronopl (5·5% – Paradis, Claude & Strimmer ). It is not clear from the literature why one particular phylogenetic method is favoured. Is BEAST chosen (i) due to historical preference, (ii) because a posterior sample of trees is desired, (iii) because it does not require a post hoc rate-smoothing step, or (iv) because it provides more accurate species hypotheses than other methods? We address the latter issue for both the GMYC and PTP by systematically evaluating their performance given different phylogenetic methods across several data sets.

6We evaluate the GMYC and PTP methods using cytochrome c oxidase subunit 1 (COI) data sets, first, where the species boundaries and diversity are relatively well known: cowries (Meyer & Paulay ), Drosophila spp. and Romanian butterflies (Dinca. ). Secondly, we compare the methods using 13 COI data sets of Rotifera, a phylum where the taxonomy is much less resolved, the sampling not as comprehensive, and where the benefit of DNA taxonomy is expected to be the greatest. We provide guidelines for maximising the robustness of species hypotheses based on single-locus data with respect to phylogenetic method.

7We compiled over 12 000 COI sequences forming 16 data sets (Table S2), corresponding mostly to genera (Rotifera + Drosophila) but also a family (Cypraeidae [cowries]; Meyer & Paulay ) and a comprehensive geographical sample comprising several families (99% of Romanian butterfly species; Dinca. ). Tree reconstruction followed standard protocols (Data S1; Fig. S1): (i) align sequences with out-groups (Table S3) using MAFFT v6.814b (Katoh, Asimenos & Toh ), (ii) remove non-unique haplotypes (for comparability the same matrix was used for all analyses, although this step is not necessary prior to generation of BEAST trees, see Talavera, Dinca & Vila ), (iii) reconstruct gene trees and (iv) make gene trees ultrametric. Gene trees were generated using distance (UPGMA – Sokal & Michener; neighbour joining – Saitou & Nei ), maximum likelihood (GARLI – Zwickl; RAxML – Stamatakis; PhyML – Guindon. ) and Bayesian inference (MrBayes – Huelsenbeck & Ronquist; BEAST – Drummond & Rambaut ). Post hoc branch smoothing (not necessary for BEAST and UPGMA trees) was performed using the R 2.15.2 (R Core Team ) package ape 3.0.7 functions (chronopl and chronos – Paradis, Claude & Strimmer ), PATHd8 (Britton. ) and r8s (Sanderson ).

8The presence of unresolved nodes and rate heterogeneity was measured directly from the trees. For each non-ultrametric gene tree, rate heterogeneity was measured as the standard deviation of the root to tip distances, where a greater standard deviation signifies greater rate heterogeneity. Analysis of BEAST trees was used to quantify whether the different species delimitation methods lead to different diversity estimates also where there are no unresolved nodes.

9The GMYC method with a single threshold (ST-GMYC), multiple thresholds (MT-GMYC; Monaghan.; Fujisawa & Barraclough ) and a multimodel approach (MM-GMYC; Powell; Fujisawa & Barraclough ) was applied to each ultrametric gene tree using the splits 1·0–11 (Ezard, Fujisawa & Barraclough ) R package. PTP analyses were performed using its webserver (http://species.h-its.org/). For each clade, up to 25 different GMYC and 30 PTP estimates were made. Primarily, the PTP analysis was used with non-ultrametric gene trees (PTP-raw: trees without post hoc smoothing, as intended by Zhang. ), but smoothed trees were also used (PTP-all: all trees) for a direct comparison with the GMYC input trees.

10The deviance of each ESU estimate from the expected diversity was gauged using the absolute difference between observed (ESUX) and expected (ESUexpected) diversity, standardised among data sets by dividing by the average diversity of the focal data set (ESUmeanA: including the focal ESU estimate). ESUexpected was either obtained from the morphological species count (ESUmorph) or the average of the species counts from across all trees (ESUmeanB: excluding the focal ESU estimate). For the three data sets where the species boundaries have been better evaluated, the morphological species count was determined using either the GenBank species name (Drosophila and Romanian butterflies) or expert advice (cowries; C. Meyer pers. comm.). In the absence of a reliable taxonomic species count for the 13 Rotifera clades, ESUmeanB was used as a conservative estimate of species richness. The use of ESUmeanB as a proxy for ESUexpected was validated by the relationship between the residual variation derived from ESUmorph and from ESUmean for the non-Rotifera data sets (File S1). Residual variation was determined for each gene tree and species delimitation method (see File S2 for examples of the calculations).

11Correspondence between ESUs and ESUmorph, in terms of species identity, was evaluated for the three non-Rotifera data sets. For each species delimitation estimate, the number of morphospecies that were split, lumped, or an exact match to an ESUmorph was counted. Exact matches are where an ESU contains all species from a single morphospecies and no other. Morphospecies are split if they are found in more than one ESU and lumped if multiple morphospecies are present within a single ESU. These counts were performed for ST-GMYC, MT-GMYC, PTP-all and PTP-raw, but not for MM-GMYC because the method returns non-integers.

12Generalised linear mixed models (GLMM; Bates. ) with a Poisson error structure were used to ascertain how residual variation varies with species delimitation methods, combination of phylogenetic and smoothing methods, rate heterogeneity, and presence of unresolved nodes. Three GLMMs were used to ask: For (i) all trees and (ii) BEAST trees, how does residual variation vary with species delimitation and phylogenetic methods? (iii) For trees with post hoc smoothing, how does residual variation vary with species delimitation and phylogenetic methods, presence of unresolved nodes and rate heterogeneity? For each of the models, residual variation was used as the response variable and clade was blocked out as a random effect.

13A generalised linear model (GLM) with a quasibinomial error structure was used to assess if the proportion of morphospecies that are an exact match to an ESU (response variable) differed among species delimitation methods, clades and combination of phylogenetic and smoothing methods (explanatory variables).

14For each of the models, significant differences among the levels were identified using post hoc Tukey HSD tests (multcomp 1.3-1 R package – Hothorn, Bretz & Westfall ). All analyses were performed in R.

15For each of the 16 clades, ten gene trees (4x BEAST, MrBayes, GARLI, PhyML, RAxML, NJ and UPGMA) and 25 ultrametric trees were generated. MrBayes analysis of the cowrie data set (1459 tips) failed to converge, leading to a total of 159 gene trees and 396 ultrametric trees. These were analysed using the ST-GMYC (396 analyses), MT-GMYC (374 analyses) and MM-GMYC (286 analyses). For the cowries, only the ST-GMYC was performed owing to computational demands; MT-GMYC ran for over a week on a 3GHz processor with 8GB RAM without reaching a local likelihood optimum. The reduced number of analyses for the MM-GMYC is due to the method not accommodating trees with unresolved nodes without manual input (the logarithm of zero-length branches produces an infinite branching rate), which would not have been achievable within the scope of the present study. In total, 475 PTP and 1056 GMYC analyses were performed (Table S3).

16For each data set, the number of ESUs estimated (non-Rotifera, Fig.1a–c; Rotifera, Fig. S2) and their residual variation (Table1; Fig.1d–g) varied among species delimitation methods. For all data sets, the PTP estimates, especially PTP-raw, best matched the expected diversity (Table1; Fig.1d–g). GMYC estimates varied depending on whether one or multiple thresholds or a multimodel approach was used (non-Rotifera, Fig.1a–c; Rotifera, Fig. S2); the MM-GMYC and MT-GMYC inferences were the most consistent (Fig.1d–g; Table1). The ST-GMYC estimates were more variable (non-Rotifera, Fig.1a–c; Rotifera, Fig. S2) and differed more from the expected diversity (Fig.1d–g; Table1). Reanalysis of these data without BEAST and UPGMA trees removes some of the significant differences associated with delimitation methods (Table1). For the non-Rotifera data set, there are no differences among the species delimitation methods. For the Rotifera data set, ST-GMYC is the only significantly different method (Table1). When only BEAST trees were analysed, there were no significant differences in residual variation among delimitation methods (Fig.2d–e; Table S4).

17Different combinations of phylogenetic and smoothing methods resulted in varied ESU estimates (non-Rotifera, Fig. S3; Rotifera, Fig. S2) and residual variation (Fig.3; Fig. S4; Table S5). The tendency to under- or overestimate diversity relative to the mean varied randomly among trees smoothed with different methods (Table S6). Gene trees smoothed with the chronopl and chronos functions typically led to highly variable ESU estimates (non-Rotifera, Fig. S3; Rotifera, Fig. S2) that differed from the expected diversity (Fig.3; Fig. S4; Table S6).

18The proportion of unresolved nodes differed among gene trees (Fig.2c), from none for BEAST trees to 24% for NJ and 43·8% for MrBayes trees. Increased residual variation is related to the presence of unresolved nodes for Rotifera (GLMM: t = 5·92, df = 1046, P < 0·0001; Fig.2b) but not non-Rotifera (GLMM: t = 0, df = 81, P = 1; Fig.2a) clades and interacts with rate heterogeneity for both Rotifera (GLMM: t = −3·27, df = 1046, P = 0·0011) and non-Rotifera (GLMM: t = −3·28, df = 81, P = 0·0015) data sets.

19The proportion of morphospecies that were inferred as an ESU did not differ significantly among the different species delimitation methods (Fig.4; Table1). Most of the combinations of phylogenetic and smoothing methods produced similar proportions of exact matches (Table S6; Fig.4), but those smoothed with chronopl or chronos produced significantly lower proportions of exact matches, resulting from either higher levels of lumping or splitting. Differences were the largest among the data sets, with a significantly higher proportion of exact matches for the Romanian butterflies than for the cowries (GLMbinomial: Z = 3·39, P < 0·001; Fig.4), but no differences when compared to the Drosophila. The proportion of exact matches was on average 63 ± 2% and was highest for the Romanian butterflies (70·4 ± 3·5%), followed by Drosophila (59·1 ± 2·6), and cowries (58 ± 4·4%; Fig.4; Fig. S5). The type of mismatches differed in proportion between the three data sets (Fig. S5): cowries were typically split, Drosophila were lumped (ST-GMYC, PTP) and split (MT-GMYC), and the Romanian butterflies were lumped.

20Good taxonomy is central to any discipline using species as a fundamental unit. Coalescent-based, phylogenetic species delimitation clusters sequences into evolutionarily significant units. This approach relies heavily on the underlying tree and is affected by the choice of phylogenetic methods (Talavera, Dinca & Vila ). Our results indicate that the PTP method produces ESU estimates that are more robust to phylogenetic reconstruction methods than the GMYC method, except when BEAST trees are used.

21Specifically, residual variation in ESU estimates was lowest for PTP-raw. The three implementations of the GMYC method differed in how robust they were to phylogenetic methods (MM-GMYC>MT-GMYC>ST-GMYC). As expected, the MM-GMYC produced ESU estimates that were more robust to different phylogenetic methods, although the MM-GMYC estimate is typically an average. Species delimitation using both the PTP and GMYC methods was consistent (lower residual variation) for BEAST trees, possibly because they require no post hoc smoothing step and contain no unresolved nodes. In contrast, analysis of NJ trees resulted in particularly large deviations in ESU estimates, irrespective of smoothing method. This is not surprising given that NJ is a clustering method that does not rely on an evolutionary model (Saitou & Nei ), known to underperform if the distance measure is not a correct estimate of nucleotide substitutions (Tateno, Takezaki & Nei ). There is also a large increase in residual variation associated with chronopl and chronos branch smoothing, which are particularly prone to haphazard lumping and splitting ESUs irrespective of the degree of between-branch smoothing (λ) chosen (File S3). This finding concurs with that of Talavera, Dinca & Vila who found that GMYC analyses of NJ trees smoothed with PATHd8, chronopl or chronos produced aberrant ESU counts.

22To quantify parameters that differ among trees and may affect ESU estimates, we analysed the effect of rate heterogeneity and unresolved nodes, which are either characteristics of poor tree reconstruction (methodological or sample issues) or real features of the data. We found a significant effect of both these parameters on the GMYC and PTP output: diversity estimates for trees with highly variable rates and/or unresolved nodes deviated more widely from the expected diversity than clocklike, resolved trees (e.g. BEAST trees). Branch smoothing of trees with highly variable substitution rates can lead to exaggerated stretching of branches (Drummond & Suchard ), which will detriment all coalescent-based species delimitation methods that use branch lengths as an input. Unresolved nodes in the tree impinge on correct diversity estimates because their resolution can lead to artefacts in the branch length data (e.g. infinite branching rates) that could result in misplaced coalescent thresholds used for delimitation. For the GMYC, splitting might occur if infinite branching rates are found closer to the tips, while for the PTP, it might result from increased average, observed intraspecific cohesiveness resulting from no increase in branch lengths with more tips. Contrarily, the diversity could be underestimated if the unresolved nodes are closer to the root for the reciprocal reasons. Whether unresolved nodes and rate heterogeneity in the data are correlates or causes of incorrect diversity estimates remains to be tested. Encouragingly, their effect is alleviated when BEAST trees are used as input.

23As a measure of how species identity differed among the methods, we assessed the proportion of ESUs that were exact matches to traditional species (morphospecies). We found similar levels of species richness to the traditional taxonomy but varying levels of discordance in identity between the traditional and DNA taxonomy. The proportion of exact matches was on average 63%, and variation in this was associated primarily, with combination of phylogenetic and smoothing methods and taxonomic group, but less with species delimitation method. We found no significant differences in the proportion of exact matches between the species delimitation methods, although the PTP method was qualitatively higher. The largest differences were between the three clades, and potentially points to the varying levels of taxonomic work in these groups. These differences seem to be driven by aberrantly deviant ESU estimates (in terms of richness and identity) associated with the use of chronopl and chronos smoothing methods, which typically split the cowrie morphospecies and lumped the Drosophila and Romanian butterfly species.

24Traditional species of the Romanian butterflies appear to be supported by DNA taxonomy perhaps because the data set represents a geographical (rather than taxonomic) sample. Species are expected to appear more distinct in such a sample because the closest relatives of most sampled species will not be sampled (Bergsten. ). Although the proportion of morphospecies that were lumped, relative to split, indicates that lower intraspecific sampling in this clade is over-representing the Yule process in the tree and thus missing some of the ESUs. The higher intraspecific sampling for the cowries and Drosophila indicates that the splitting of these species could be associated with unresolved taxonomy (Packer. ) or overlapping intra- and interspecific variation (Meyer & Paulay; Wiemers & Fiedler ). While efforts have been made to resolve the taxonomy of these groups (Meyer & Paulay; O'Grady & Markow ), a more concerted effort is required to address the gap between DNA and traditional taxonomy across the entirety of these clades (C. Meyer pers. comm).

25By assessing ESU counts across 16 data sets with over 1500 separate species delimitation analyses, we have shown that the PTP-raw model with any robust gene tree and the GMYC used on BEAST trees produce consistently robust and, on average, accurate species estimates. These findings can probably be extrapolated to other genetic markers: COI and 18S are typically used for animals (Tang., ), multiple markers (e.g. 16S) for bacteria (Barraclough.; Morlon. ), ITS for fungi (Powell. ) and multiple markers (e.g. matK and rbcL) for plants (Hollingsworth & CBOL Plant Working Group ). Although the variability of these markers will likely yield different degrees of coalescent clustering and species separation (Tang. ) that warrants a more thorough evaluation.

26Coalescent-based species delimitation is likely to gain in popularity: either to facilitate the description of biodiversity in an integrative, iterative way as a tool to tackle the burgeoning taxonomic crisis (Puillandre. ), or to cluster sequences from NGS studies (Creer.; Chariton. ). The latter would benefit from evolutionary approaches that provide a deeper understanding of the nature and extent of diversity (Barraclough. ). Applying coalescent-based species delimitation to NGS is currently limited by the amount of variability, the short (but ever increasing) read lengths, the amplification success of the markers used and the computational expense of the coalescent-based metrics. As with all DNA taxonomy studies, primers need to be designed to combat the low amplification success of certain primers (Zhan. ), robust bioinformatics pipelines need to be developed (S. Creer pers. comm.), and sampling regimes that are representative of intra- and interspecific variability and geographical range should be considered (Papadopoulou.; Lohse; Bergsten.; Talavera, Dinca & Vila ).

27The PTP method is appealing when speed is essential because ultrametric trees are not required (Zhang. ), meaning that some of the problems encountered and the additional computation required with branch smoothing may be circumvented. However, the PTP makes the assumption that branching events scale with substitutions rather than time, which might be violated when substitution rates are heterogeneous. The GMYC with a BEAST tree provided equally consistent estimates but obtaining BEAST trees is computationally expensive. However, when rate heterogeneity is high and can be adjusted across the tree estimation using models, perhaps by use of well-informed internal calibration priors, then diversity estimation might benefit from sophisticated dating and diversity estimation procedures. We feel that the GMYC is more true to the speciation process, in that speciation and coalescence happen over time and not necessarily in relation to how many substitutions occur in marker genes. While the transition between speciation and coalescent processes, used to delimit species, might be mirrored by differences in the number of species- and population-level substitutions (PTP), time is a more direct expression of the process; therefore, methods that separately model this transition (time; GMYC) and phylogenetic methods that formally correct for substitution rate variation among species (e.g. BEAST) are conceptually more appropriate. We recommend use of both PTP and GMYC methods with the appropriate phylogenetic tree or choosing between them on a case-by-case basis, bearing in mind the differences in speed and underlying theory inherent in the two methods. The PTP method with non-ultrametric trees is currently quicker to implement than the GMYC, especially the MM-GMYC, although the speed of the GMYC could be increased with parallelisation. Both the phylogenetic and species delimitation steps become computationally demanding for larger data sets (e.g. NGS studies). Such data sets, which are often taxonomically broad, are likely to violate use of a single substitution rate and so increased parameterisation and prior information is more likely to yield trees that better reflect the data and thus provide more realistic diversity estimates. We envisage that better phylogenetic handling of substitution rate heterogeneity within the samples, irrespective of delimitation method, and the use of ESU nodal support as a proxy for species identity confidence, would further improve the delimitation of primary species hypotheses from single-locus marker surveys.