Delineating Species with DNA Barcodes: A Case of Taxon Dependent Method Performance in Moths

The accelerating loss of biodiversity has created a need for more effective ways to discover species. Novel algorithmic approaches for analyzing sequence data combined with rapidly expanding DNA barcode libraries provide a potential solution. While several analytical methods are available for the delineation of operational taxonomic units (OTUs), few studies have compared their performance. This study compares the performance of one morphology-based and four DNA-based (BIN, parsimony networks, ABGD, GMYC) methods on two groups of gelechioid moths. It examines 92 species of Finnish Gelechiinae and 103 species of Australian Elachistinae which were delineated by traditional taxonomy. The results reveal a striking difference in performance between the two taxa with all four DNA-based methods. OTU counts in the Elachistinae showed a wider range and a relatively low (ca. 65%) OTU match with reference species while OTU counts were more congruent and performance was higher (ca. 90%) in the Gelechiinae. Performance rose when only monophyletic species were compared, but the taxon-dependence remained. None of the DNA-based methods produced a correct match with non-monophyletic species, but singletons were handled well. A simulated test of morphospecies-grouping performed very poorly in revealing taxon diversity in these small, dull-colored moths. Despite the strong performance of analyses based on DNA barcodes, species delineated using single-locus mtDNA data are best viewed as OTUs that require validation by subsequent integrative taxonomic work.


Introduction
After little progress over a long interval, the past decade has seen the development of several analytical methods which employ DNA sequences to delimit species boundaries [1][2][3][4][5][6][7][8][9][10]. Another innovation, DNA barcoding [11,12], was originally developed for specimen identification using a standardized segment of the mitochondrial genome (a 648bp region of the cytochrome c oxidase subunit I, COI), but its utility for species discovery was soon recognized [13][14][15][16]. The results showed some differences among higher taxa (e.g., North American plusiine moths), but much more detailed investigation of these effects is merited. This study extends past work by comparing the performance of five methods with two groups of moths, Finnish gelechiines and Australian elachistines. The Finnish gelechiine fauna includes 180 species belonging to 49 genera ( [51], J. Kullberg, pers. comm.), while the Australian elachistines include 148 species in three genera [52]. Both subfamilies are members of the Gelechioidea, one of the largest radiations within the Lepidoptera, and provide a well-defined set of reference species established following detailed morphological and ecological studies [52][53][54]. Because these moths are generally small and dull-colored, they are very difficult for taxonomic studies so there are many undescribed species [55][56][57]. The two groups do have one difference; the Australian elachistines present a challenge for DNA barcoding due to their close affinities and supposed recent origins [58], while barcodes have successfully discriminated many European gelechiines (e.g., [14,59,60]). The adoption of DNA barcode-based methods for the delineation of species in the Gelechioidea and other taxa sharing their biological attributes has great potential to accelerate species delineation. Executing identical analyses with two datasets, one more challenging than the other, provides an opportunity to evaluate the performance of the differing methods in this context.
This study employs four commonly used methods for DNA-based species delineation, including one older but still widely used approach, statistical parsimony networks (TCS), and three recent methods: Barcode Index Numbers (BIN), Automatic Barcode Gap Discovery (ABGD) and Generalized Mixed Yule Coalescent (GMYC). These methods were selected for inclusion based on their general popularity and their strong performance in a previous study [9]. BIN analysis always generates only one number of OTUs for each set of DNA sequences, while the other approaches do not because key parameter values (TCS and ABGD) or input trees (GMYC) can vary. In addition to comparing the performance of these DNA-based approaches, we obtained results from a morphology-based analysis using external characters. We study differences between the outcomes for two datasets, considering both the count and composition of the putative species (OTUs) produced by each analysis. Finally, we evaluate the performance of the methods with singletons, as well as monophyletic and non-monophyletic (i.e., para-and polyphyletic) species.

Taxon sampling
Specimens of 92 species of Australian Elachistinae were sampled from the Australian National Insect Collection (ANIC) and the Finnish Museum of Natural History (MZH). Specimens of 103 species of Finnish Gelechiinae (tribes Teleiodini, Gelechiini and Gnorimoschemini) were sampled from the private collection of M.M. during 2008-2012 (17 of the latter specimens were collected from Denmark, Estonia, France, Latvia, Russia and Sweden). One additional elachistine specimen was analyzed from the Agricultural Scientific Collections Unit (ASCU), and three specimens of gelechiines from the private collection of Erkki Laasonen. Three to five specimens per species were usually sampled when available, targeting recently collected individuals from diverse geographic localities. Relatively few specimens of each taxon were analysed to maximize the species coverage. Larger sample sizes were examined for a few species whose taxonomic status is controversial. One or two legs were removed from dry pinned specimens for DNA extraction. Specimens were identified by L.K. (Elachistinae) and M.M. (Gelechiinae) following the taxonomy of Kaila [52] and Huemer and Karsholt [53,54], respectively. BOLD Sample and Process IDs, GenBank accession numbers and other details of our sequence data can be retrieved from S1 Table. DNA extraction, PCR amplification and sequencing DNA extraction, PCR, and sequencing were performed at the Canadian Centre for DNA Barcoding following standard high-throughput protocols [61]. The first round of PCR employed the primers LepF1 and LepR1 [11] which generate a 658bp amplicon that spans the barcode region of COI. In cases of failure, two additional PCR reactions were carried out to recover 306bp and 407bp amplicons using a standard primer set [62]. If one of these reactions was successful, an effort was made to obtain a barcode compliant record (>497bp) by amplifying shorter regions of COI using the primer sets described in Hebert et al. [63]. All sequences were aligned using the BOLD Aligner in the Barcode of Life Data Systems (BOLD [64]) and then inspected visually for stop codons and frameshift mutations in MEGA5 [65].

Comparison between datasets
Several attributes were studied to expose differences between the two datasets. Intra-and interspecific pairwise distances were calculated in the BOLD workbench employing the "Barcode Gap Analysis" tool, and visualized using the "sppDist" function in SPIDER [66] available in R [67]. The incidence of monophyly was quantified using the "monophyly" function of SPIDER. Pairwise distances for all sequences included in the analysis were calculated using a K2P distance model in MEGA5.

Morphological sorting
In order to simulate the process of species recognition through morphological sorting, we recruited an experienced lepidopterist (M.N.) without previous knowledge of Australian elachistines or Finnish gelechiines to sort pinned specimens into OTUs using external morphology, mainly wing patterns, an approach similar to that employed in previous studies (e.g., [43]). The test collection included from one to five specimens of 96 species of Elachistinae and 83 species of Gelechiinae with representatives of most species that were used for DNA barcoding and a few additional taxa ( Table 1). The individuals included were not DNA barcode vouchers, but other similarly identified museum specimens.

OTU delineation based on DNA barcodes
The Barcode Index Number System (BIN [9]), statistical parsimony networks (TCS [19,20]) and Automatic Barcode Gap Discovery (ABGD [8]) rely on pairwise sequence distances between specimens to determine the number of OTUs within a dataset. The RESL algorithm, which forms the basis of the BIN system, employs a three-stage procedure which starts with single linkage clustering using a fixed 2.2% threshold. This phase is followed by Markov clustering, which aims to improve the accuracy of the OTUs, and finally the Silhouette criterion compares the different clustering schemes from Markov clustering and chooses the option with the highest Silhouette score. ABGD employs a two-phase system which initially divides sequences into OTUs based on a statistically inferred barcode gap (i.e., initial partitioning), and subsequently conducts a second round of splitting (i.e., recursive partitioning). ABGD has three key parameters: (i) X, which is an estimate of relative gap width, and (ii) minimum and (iii) maximum values of prior intraspecific divergence (P), which are used to detect the barcode gap. The default P-values typically produce a range of OTU counts. TCS produces the most parsimonious solution for a particular cut-off value (90-99% cut-off values are available) by combining pairs of specimens with the lowest genetic distances. The procedure continues until the cut-off value is exceeded. The higher the cut-off, the smaller the number of steps needed to exceed it and the greater the count of unconnected networks recognized. In other words, selecting a high cut-off value produces a high species count and vice versa. The Generalized Mixed Yule Coalescent (GMYC [3,4,21]) differs strongly from the other methods because it is a model-based approach, aiming to discover the maximum likelihood solution for the threshold between the branching rates of speciation and coalescent processes on a tree. The number and composition of OTUs is inferred by counting the lineages crossing the threshold. The BIN analysis was done using a stand-alone version of RESL which is scheduled for public release in the near future. Standard BIN assignments are available on BOLD v3.6 (http:// www.boldsystems.org), but they are generated through the analysis of all barcode sequences on BOLD, meaning that the results are not strictly comparable with those obtained with other methods (because they are based on a more inclusive dataset). Statistical parsimony networks were calculated using software TCS v1.21 [20] with separate analyses for ten cut-off values (90%-99%). ABGD analyses were performed on 24-25 March 2013 on the web interface (http://wwwabi.snv.jussieu.fr/public/abgd/). Because the default value for relative gap width (X = 1.5) did not produce a result for either dataset, two lower values (X = 0.8, 1.0) were used. 1.0 was the highest value that could be applied as ABGD did not produce results for the Gelechiinae dataset with X = 1.1. ABGD provides the option of using three distance metrics: Jukes-Cantor (JC [68]), Kimura 2 parameter (K2P [69]) and simple p-distances. We conducted analyses using all three metrics with both values (0.8, 1.0) of X, resulting in six analyses per dataset. All results using prior limits for intraspecific divergence ranging from P = 0.001-0.1 were recorded. Defaults were employed for all other parameter values.
GMYC requires a fully-resolved ultrametric chronogram as input. In order to test the effect of different input trees on GMYC, we calculated chronograms using three approaches: unweighted pair group method with arithmetic means (UPGMA [70]) and two Bayesian inference gene trees constructed with a Yule pure birth model [71,72] and constant size coalescent [73] tree priors. UPGMA trees have rarely been used with GMYC analyses [25,74], likely reflecting concerns with their effectiveness in phylogeny estimation (e.g., [75]), but they are an attractive option because of their speed and simplicity. The UPGMA trees used in this study were constructed in MEGA5 using a K2P distance model. Model selection for Bayesian analyses was performed a priori with jModeltest v.0.1.1 using the Akaike information criterion (AIC) [76] and a posteriori with Bayes factors implemented in Tracer v.1.5 [77]. GTR+G+I was the preferred model for Elachistinae with both methods, and for the Gelechiinae with Bayes factors, but not with jModeltest where HKY+G+I was preferred. As GTR+G+I was the second ranked option for jModeltest, the same model was employed for both datasets. The fit of clock models and tree priors were also estimated using Bayes factors with preference for the uncorrelated relaxed lognormal clock model over a strict clock and coalescent tree prior over Yule prior. Bayesian inference trees were constructed using BEAST [78,79]. XML files (S1-S4 Appendices) were made with the BEAUti v1.7.1 interface with the following settings: GTR+G+I substitution model; empirical base frequencies; 4 gamma categories; all codon positions partitioned with unlinked base frequencies and substitution rates. An uncorrelated relaxed lognormal clock model was used with rate estimated from the data and ucld.mean parameter with uniform prior employing 0 as the lower and 10 as the upper boundary. All other settings employed defaults. The length of the MCMC chain was 40 000 000 sampling every 4000. All BEAST runs were executed in BioPortal [80] and the resultant ESS values and trace files of runs were evaluated in Tracer. Two independent runs were combined using LogCombiner v.1.7.1 with 20% burn-in. Maximum clade credibility trees with 0.5 posterior probability limit and node heights of target tree were constructed in TreeAnnotator v1.7.1. Both single-and multiple-threshold GMYC analyses were conducted in R with the packages APE [81] and SPLITS [82]. All analyses related to GMYC were performed with haplotype data collapsed in ALTER [83].

Direct examination of OTU composition
A simple comparison of OTU counts to the number of reference species can be misleading because similar results can be produced by varying levels of congruence between species and OTU boundaries if splits and merges are counterbalanced. In order to acquire a deeper insight, we estimated the correspondence between the boundaries of OTUs and reference species by assigning each OTU as a MATCH, SPLIT, MERGE or MIXTURE [9]. A MATCH results when the specimens assigned to an OTU include all those assigned to a reference species. By contrast, a SPLIT represents the case where members of a reference species are divided into two or more OTUs, while a MERGE represents the case where two or more reference species are assigned to a single OTU. MIXTURES involve more complex cases where members of two or more reference species are involved in both merger and splitting. Each OTU can only be assigned to one of these categories. The performance of each method with singletons as well as with monophyletic and nonmonophyletic species was studied by dividing datasets according to the results of the monophyly analysis by SPIDER (see Comparison between Datasets), and conducting a direct examination of congruence as described above. No additional OTU delineation analyses were executed with partitioned data.

Results
A total of 562 full-length sequences (654bp; the original length 658bp is reduced by the BOLD aligner as it removes the first and three last bases) were recovered. These included 307 sequences (187 haplotypes) from 92 species in 25 genera of Gelechiinae and 255 sequences (178 haplotypes) from 103 species (including 6 undescribed species) in two genera of Elachistinae (Table 1). These datasets provide coverage for all Finnish gelechiine species in the tribes Teleiodini, Gelechiini and Gnorimoschemini, and for 65.5% of all Australian elachistines. We only included full-length sequences to minimize the possible effects of missing bases on the outcomes of subsequent analyses. All sequences are available in public databases (for GenBank accessions see S1 Table; BOLD

Dataset comparison
Intraspecific distances in the Gelechiinae varied from 0.00% to 2.94% (mean = 0.39%, SE = 0.01) while distance to the nearest neighbor (NN) species ranged from 0.92% to 11.25% (mean = 6.33%, SE = 0.02) so there were few cases of overlap between intra-and interspecific distances (Fig 1). A similar pattern was observed in the Elachistinae with intra-specific distances ranging from 0.00% to 2.3% (mean = 0.28%, SE = 0.01), while NN distances varied from 0.00% to 11.02% (mean = 3.48%, SE = 0.03) (Fig 1). Because these distance measures reflect past decisions on species boundaries (which may be incorrect), they may be biased, but this effect can be reduced by calculating pairwise distances without a priori grouping. This analysis confirmed that the average distance among all sequences was lower for Australian elachistines (mean = 0.099) than for Finnish gelechiines (mean = 0.13) (Fig 2). The proportion of monophyletic groups was also very different: 99% of the Gelechiinae species were monophyletic (85 monophyletic, 1 non-monophyletic, 6 singletons), but only 75% of the Elachistinae (52 monophyletic, 17 non-monophyletic, 34 singletons).

Morphological sorting
The number of putative species resulting from morphological sorting was close to the reference count for both the Gelechiinae (91 vs. 83 reference species) and the Elachistinae (97 vs. 96). However, the composition of the OTUs showed a poor match with accepted taxonomy (Fig 3). Only 29% of the Gelechiinae species were correctly assigned, and 17% of the Elachistinae. As well, a very high proportion (40%) of the OTUs in both subfamilies represented MIXTURES of two or more species.

OTU counts
OTU counts produced by the DNA-based delimitation methods ranged from 90 to 122 for the Gelechiinae (Figs 4a, 5a and 5b) and from 27 to 159 for the Elachistinae (Figs 4b, 5a and 5b). Only one method generated the same OTU count as the number of reference species (92) for the Gelechiinae: Automatic Barcode Gap Discovery with relative gap width (X = 1.0) and prior intraspecific divergence (P = 0.0215). None of the methods generated the same OTU count as the number of reference species (103) for the Elachistinae. The relative number of OTUs versus   Taxon Dependent Performance of DNA-Based Delineation Methods the reference species count varied between the two datasets: most OTU counts for the Gelechiinae were higher than the reference count of 92 species (Fig 6a), whereas most for the Elachistinae were lower than the reference count of 103 species (Fig 6b).
BIN: BIN generated a single outcome for each dataset, recognizing 97 OTUs for the Gelechiinae, and 84 for the Elachistinae. (Fig 4a and 4b). In general, the BIN results followed the same pattern as the other methods with a higher OTU count than the reference sequence count for the Gelechiinae and a lower OTU count for the Elachistinae (Fig 4a and 4b). When compared with the other methods, the BIN results were approximately in the middle of the performance plots (orange bars in Fig 6a and 6b).
TCS: TCS was used with all available cut-off values (90-99%, resulting in ten outcomes. As expected, the lowest cutoff (90%) always generated the fewest clusters while the highest (99%) generated the most (Gelechiinae: 95-113, Elachistinae: 72-94) (Fig 4a and 4b). The relative position of the OTU counts differed between the two datasets: the TCS results were scattered among those for the other methods in the Gelechiinae, while those for the Elachistinae had low values (green bars in Fig 6a and 6b). All results from TCS were higher than the reference species count for the Gelechiinae, but lower for the Elachistinae (black vs. green bars in Fig 6a and  6b).  ABGD: ABGD was used with two values of relative gap width (X) and three distance metrics (p, JC, K2P). All OTU counts resulting from varying values of prior intraspecific divergence (P) were recorded (Fig 5a and 5b). All analyses produced zero OTUs for the Gelechiinae when P = 0.0599 and the initial partition with 90 OTUs was reached when P = 0.0359. All distance metrics behaved similarly with the Gelechiinae generating OTU counts ranging from 90 to 122 (X = 0.8) and from 90 to 115 (X = 1.0) (Figs 4a, 5a and 5b). The pattern changed with the Elachistinae dataset (Figs 4b, 5a and 5b). The two values of X produced differences for the initial partitions: either one (X = 1.0) or two (X = 0.8) OTU counts were generated by the initial partition (Fig 5a). ABGD behaved similarly with recursive partitions (Fig 5b), but there were differences between the three distance metrics. P-distance produced OTU counts ranging from 27 to 138 (X = 0.8) and from 110 to 125 (X = 1.0), while JC (X = 0.8: 60-122; X = 1.0: 88-117) and K2P (X = 0.8: 76-122; X = 1.0: 94-118) generated more constrained counts when X = 0.8 (Fig  5a). ABGD was the only method to produce the same OTU count as the number of reference species for the Gelechiinae with high prior intraspecific divergence value (P = 0.0215). By contrast, the closest match (100 OTUs) for the Elachistinae was generated by two low P-values (0.00278 and 0.00167) (Fig 6a and 6b).
GMYC: GMYC was used with three input trees: UPGMA and two Bayesian chronograms constructed in BEAST with Yule and coalescent tree priors. The results of both single-and multiple-threshold models were recorded (Table 2), although only one analysis indicated a better fit for the multiple-threshold model (UPGMA starting tree with Elachistinae; χ 2 = 93.22, d. f. = 21, P<<0.001). The likelihood ratio test was highly significant for all analyses, indicating rejection of the null model (OTU count = 1). The single-threshold model generally produced lower cluster counts than the multiple-threshold for the Gelechiinae and the starting trees constructed in BEAST resulted in lower counts than the UPGMA trees (Fig 4a). All GMYC analyses recognized more OTUs than the reference species count for the Gelechiinae (97-111 OTUs vs. 92 species) (purple bars in Fig 6a). GMYC behaved differently with the Elachistinae dataset. The single-threshold analysis based on an UPGMA starting tree recognized 159 OTUs, which was far more than the other analyses (Fig 4b). The multiple-threshold model based on the UPGMA tree also generated a high OTU count (110), whereas the GMYC analyses with Bayesian trees were more stable, recognizing from 93 to 96 OTUs (Fig 4b). The reference species count was between UPGMA and BEAST results (black vs. purple bars in Fig 6a and 6b).

Method performance based on OTU composition
The percentage of MATCHES was generally higher for the Gelechiinae (63-96%) than for the Elachistinae (46-77%; logistic regression: estimate = -1.29, n = 67, P<0.0001), reflecting the increased proportion of MERGES and MIXTURES for the Elachistinae (Fig 7a and 7b). The results generated by ABGD with p-distance were excluded from these statistical tests due to the strongly discordant results generated by this method for the Elachistinae (MATCH: 18-70%). The highest percentage of MATCHES was produced by ABGD for both datasets, but with very different P-values (Gelechiinae: p-distance, JC, K2P, X = 0.8, P = 0.0215; Elachistinae: JC, X = 1, P = 0.001) ( Table 3, Fig 7a and 7b). ABGD also generated the lowest percentage of MATCHES for the Elachistinae (Fig 7b), while the poorest result for the Gelechiinae was produced by the multiple-threshold GMYC with UPGMA input tree (Fig 7a). BIN: BIN produced a high percentage of MATCHES for the Gelechiinae (90%), but substantially less for the Elachistinae (67%).
TCS: TCS generated a varying proportion of MATCHES depending on the cut-off value. The highest percentage (91%) of MATCHES for the Gelechiinae was obtained with 92% and 93% cut-off values, while the best result (70%) for the Elachistinae was obtained with a cut-off value of 98%.
ABGD: The initial partition produced 89% MATCHES for the Gelechiinae irrespective of distance metric or value of relative gap width. The highest percentage of MATCHES was generated by P = 0.0215 (X = 0.8, 96%; X = 1.0, 93%) and the lowest by P = 0.001 (X = 0.8, 83%; X = 1.0, 84%). By contrast, the two values of relative gap width and the different distance metrics had clear effects on the performance of ABGD for the Elachistinae. The two initial partitions produced by X = 0.8 differed in their percentage of MATCHES (P = 0.001-0.00278: 64-71%; P = 0.00464-0.00215: 18-57%). X = 1.0 generated only one initial partition, which was the same as the partition with lower P-values of X = 0.8. JC (46-77%) and K2P (57-75%) generated a similar percentage of MATCHES, although JC was more variable. P-distance performed very poorly (18-70%), especially with P-values from 0.00464 to 0.0215. The most congruent outcome with the reference species (77% MATCHES) was produced by P = 0.001 with JC and X = 1, but the same P-value generated also the second highest percentage of MATCHES with K2P and both values of X.

Singletons, mono-and non-monophyletic species
Most singletons in the Gelechiinae dataset matched with their corresponding reference species (BIN, all TCS, all GMYC with single-threshold, most ABGD with X = 0.8: 100%), whereas the     percentage of MATCHES varied much more for the Elachistinae (Fig 8a and 8b). BIN and TCS with 95% cut-off produced a similar percentage of MATCHES for the Elachistinae (74% and 76%, respectively). Other results of TCS varied from 68% (90% cut-off) to 88% (99% cut-off). The highest percentage of MATCHES for the Elachistinae was produced by GMYC with UPGMA starting tree and single-threshold model (94%) (S2 Table); the rest varied from 76% (UPGMA, multiple-threshold) to 88% (BEAST with coalescent prior, single-threshold) (Fig 8a  and 8b). The results of ABGD spanned a wide range from 26% (p-distance, X = 0.8, P = 0.0215 Initial) to 94% (JC, X = 0.8, P = 0.001; K2P, X = 0.8, P = 0.001) (Fig 8a and 8b). The examination of mono-and non-monophyletic (i.e., para-or polyphyletic) species revealed that none of the methods was effective in delimiting non-monophyletic taxa. The only exception was ABGD with small P-values, which managed to deliver one or two MATCHES. On the other hand, ABGD also produced a high number of MIXTURES together with these MATCHES. As only one species in the Gelechiinae dataset was non-monophyletic and six were singletons, the difference was minor between the analyses including all data and only monophyletic species (all data: mean = 86.5%; monophyletic: mean = 87.3%; logistic regression: estimate = 0.07, n = 58, P = 0.4207). However, the performance was significantly improved for the Elachistinae (all data: mean = 63.9%; monophyletic: mean = 73.9%; logistic regression: estimate = 0.47, n = 76, P<0.0001) (Fig 9a and 9b). The taxon-dependent pattern in general performance remained even after removing all non-monophyletic species and singletons as the percentage of MATCHES was still significantly higher for the Gelechiinae than for the Elachistinae dataset (Fig 9a and 9b; logistic regression: estimate = -0.89, n = 67, P<0.0001). The performance of individual methods when non-monophyletic species and singletons were excluded was rather similar to the performance revealed by all data (Fig 9a, 9b and S3 Table).

Discussion
This study has compared the performance of five species delineation methods (BIN, TCS, ABGD, GMYC, and external morphology) with two groups of Lepidoptera: Finnish Gelechiinae and Australian Elachistinae. The difference between the groups was evident as the Gelechiinae had a wider barcode gap (Fig 2) and included more monophyletic species than the Elachistinae. The Gelechiinae also seem more morphologically diverse as indicated by their assignment to 25 genera, while all but one of the 92 species of Australian elachistines are in a single genus. Our results reveal a striking difference between the two taxa in the effectiveness of varied delineation methods in recovery of current species boundaries. Performance was generally higher for all methods with the Gelechiinae than the Elachistinae. The range between the   lowest and highest OTU counts was smaller, and the percentage of MATCHES was also higher for the Gelechiinae than for the Elachistinae. The higher proportion of MATCHES for the Gelechiinae remained evident when monophyletic species and singletons were studied separately, although the performance improved in all methods.
Explaining the differential success in species delineation between Gelechiinae and Elachistinae The direct examination method depends upon the reliability of the species boundaries in the reference species. Cases of discordance between the boundaries of reference species and OTUs can arise in two ways. They can reflect errors between true species boundaries and those recognized by current taxonomy. In these cases, the reference species have been wrongly delineated, and a discordant OTU might reveal the true species boundary. Mistakes caused by taxonomic errors can be corrected when discovered. However, cases of discordance can also arise due to biological factors when a reference species corresponds with the true species boundary, but OTU delineation methods cannot recover it because of weaknesses in the analytical method or in the underlying data. Taxonomic errors reflect the subjectivity that is often involved in drawing the line between two species. In addition to errors caused by insufficient knowledge of the focal taxa, taxonomic errors can arise through the use of unsuitable characters. For instance, an inappropriate reliance on wing venation led to oversplitting in one species complex of European elachistids [84]. As a general rule, the accuracy of species delineation improves as a particular group experiences recurrent study.
One major biological factor is the age of species. Young species tend to have unclear boundaries because processes such as hybridization and introgression are ongoing. As ecological differences can arise very rapidly via divergent selection [85,86], their use in species diagnosis can enhance the discovery of young species. Because DNA barcodes are single-locus data from the mitochondrial genome, they usually cannot recover OTUs which follow true species boundaries in cases where introgression is prevalent [87]. In general, the delineation of such species is challenging for all methods and a combination of several types of characters coupled with a careful sampling scheme is required (e.g., [88,89]). However, this approach is not an option for the first phase of species discovery when synoptic methods, such as DNA barcodes, are the best tools.
In addition to taxonomic and biological factors which influence the precision of the reference species, the sampling scheme can also impact the performance of the delineation methods. Optimally, a sufficiently large number of specimens covering the whole geographic distribution of each species would be included, and the study would include all known species of the focal clade. Unfortunately, this optimal scenario is rarely feasible, especially when a high number of poorly-known species is included. However, it should be noted that including a small number of specimens per species and/or sampling only distantly related species may artificially improve the congruence between species and OTUs. Restricted geographic sampling can have a similar effect, although the impact of geographical distance on the intraspecific variance has shown taxon-dependence [90][91][92]. Hausmann et al. [93] studied the performance of the BIN system in the well-studied geometrid moths with a large-scale sampling scheme, covering most parts of Europe. They reported rather poor performance (67% MATCHES), but concluded that many cases of discordance reflected flaws in the current taxonomy rather than problems with the method. Restricted intraspecific sampling can also raise the number of singletons in the data, which might affect the method performance. However, no effect of this type was discovered here (see Fig 8) or in a previous study evaluating the performance of GMYC [48].
Eight species of Gelechiinae and 32 species of Elachistinae were delineated differently in this study than in current taxonomy (i.e., three or four out of four methods produced discordant results, see species marked with asterisks in Table 1). To evaluate the effect of the sampling scheme, these species were studied for their number of BINs (i.e., OTUs delineated by the RESL algorithm) on BOLD (results in S4 Table). As these BINs are based on all sequence data on BOLD, the sampling effort was increased for most gelechiine species. No conflicts between the results in this study and the BINs were revealed. As additional specimens were available for seven of eight species from various parts of Europe and North America, sampling-based error is an unlikely explanation. Instead, it is possible that these reference species, especially the five species which were SPLIT in this study, each reflect a case of overlooked species and, thus, the discordance observed between current species boundaries and OTUs reflects a taxonomic error. However, two gelechiine species (Scrobipalpa artemisiella and S. stangei) were MERGED in the same OTU here as well as in the same BIN on BOLD. This case of discordance might be due to a biological factor as these species can easily be separated by morphology and have different life histories. Sampling-based error could not be evaluated for one gelechiine species (Scrobipalpa bryophiloides) as no additional sequences were available.
The examination of the BIN records on BOLD provided only a few additional records for the Australian Elachistinae so the effect of sampling could not be evaluated. However, many discordant results for the Elachistinae species were MERGES, leading to a higher number of specimens per OTU. As a consequence, the overall sample size per OTU was larger for the Elachistinae than the Gelechiinae. Seventeen elachistine species formed three groups with highly discordant results between the reference species and the delineated OTUs (see asterisks with numbers in Table 1). E. lachnella was included to this comparison, because it was MERGED in the same BIN with E. nodosae. This discordance between the result of this study and the BIN was due to one intermediate E. nodosae specimen on BOLD. As these three groups included 36, 12, and 10 sequences, sampling effort was not low so sampling-based error is unlikely. Instead, both biological and taxonomic factors may explain the observed discrepancy. These species were originally delineated based on ecological traits (in particular phenology, host plant selection and the shape of leaf mines that their larvae produce) which were correlated with small, but consistent morphological differences, a pattern compatible with their recent origin. Their young age is further supported by the distribution of pairwise genetic distances in Fig 2. However, as these species were described very recently (2011), they have not yet experienced critical re-examination so the species hypotheses cannot be considered as fully robust. As well, their young age might reflect taxonomical complications such as introgression which would make species boundaries difficult to interpret.
Fourteen Elachista species, which were originally delineated based on morphological differences, were MERGED with their sister species. In one case (E. ophelma and E. catagma), the representatives of each species formed a distinct subcluster within the shared OTU. Two other species (E. gerasmia and E. physalodes) showed a similar pattern, but one specimen was grouped with the subcluster otherwise comprised solely of its sister species. In three cases, the subcluster of one species was nested within its sister species (E. anolba and E. averta; E. zophosema and E. litharga; E. stictifica and E. platysma). Two pairs (E. ophthalma and E. nr. ophthalma; E. sp. ANICLK1 and E. sp. ANICLK3) included undescribed, but morphologically distinct species, which were MERGED with their sister species. These cases may also reflect both biological and taxonomic factors associated with recently diverged taxa. Only two Elachista species were designated as SPLITS (E. carcharota and E. discina). The split within E. carcharota was associated with large geographical distance (western vs. eastern Australia), the specimens were originally deemed conspecific because morphological differences were not apparent. E. discina was divided into two OTUs from sites in proximity. Both cases of SPLITS may reflect problems introduced by the small number of samples per species, but the possibility of overlooked cryptic species cannot be excluded.

Method performance
Some general patterns in method performance were present regardless of the taxon. BIN, TCS with cut-off value 95%, and GMYC with Bayesian input trees and the single-threshold model produced similar results for both datasets (97-102 OTUs, 89-91% MATCHES for Gelechiinae; 81-96 OTUs, 62-67% MATCHES for Elachistinae). GMYC analyses based on Bayesian trees and BIN performed slightly better than TCS (95% cut-off), especially with the Elachistinae. As the performance of GMYC and BIN was similar for the Elachistinae, there was no evidence to support Zhang et al.'s [10] contention that tree-based approaches are superior for taxa with a narrow barcode gap.
The performance of GMYC was found to be very sensitive to the starting tree. UPGMA trees produced poor results with regard to both OTU count and composition (Figs 4-6, Table 3). Similar sensitivity has been reported in previous studies, which have tested this method with different trees [39,40,48,94]. Tang et al. [94] noted that starting trees transformed to ultrametric by post hoc branch smoothing (e.g., by employing function 'chronopl' in R) perform especially poorly. This feature complicates the use of GMYC with large datasets, because computationally expensive BEAST trees seem to be the only reliable option. Another noteworthy feature of GMYC is the weak performance of the multiple-threshold model which was also detected in a previous simulation study [21].
ABGD produced both the highest and some of the lowest percentages of MATCHES for both datasets. As different prior intraspecific divergence (P) values (when used with default parameters) lead ABGD to generate variable OTU counts, it would be optimal to choose a fixed P-value to gain just one result. P = 0.01 has been proposed [8], and this value performed well for the Gelechiinae (90-92% MATCHES), but poorly for the Elachistinae (50-59% MATCHES). We conclude that the adoption of one P-value can result in either high or low performance, depending on the focal group. Without a fixed P, ABGD generates a range of outcomes, meaning that the user must choose the 'correct' result, compromising the objectivity of this DNA-based method.
ABGD also showed considerable sensitivity to the distance metric adopted. The results with p-distance were most discordant, but K2P and JC also produced variable results for the Elachistinae. Similar discordance was observed in a study on Australian hypertrophine moths [17]. As the effect of distance metric on barcoding results is minor [95], the divergent OTU counts arising from the use of different distance metrics in ABGD seems to reflect a methodological problem. Puillandre et al. [8] have noted that ABGD requires 3-5 specimens per species for ideal performance, but this criterion is difficult to meet, especially if the number of taxa is unknown. As this issue has been discussed earlier [17], we only point out that our intraspecific sampling was mostly below this minimum level, but the general performance was still mainly congruent with the other tested methods.

External morphology vs. DNA barcodes
Groupings based on external morphology has been the primary method for species delineation for centuries. The results from the morphological sorting in this study generated a low percentage of MATCHES for both subfamilies and a high proportion of MIXTURES, a particular challenge for subsequent taxonomic work. Although this test involved low sample sizes, it still provides an estimate of the relative efficacies of OTU designation via external morphology versus DNA-based methods. As with many other gelechioid moths, both elachistines and gelechiines are small and dull-colored, often lacking clearcut differences in external morphology. Furthermore, many elachistine species are sexually dimorphic [52], which might have contributed to the lower number of MATCHES for the Elachistinae. As the superfamily Gelechioidea includes many undescribed species, the need for efficient tools to expedite taxonomic workflows is of high importance. The present study reveals that the sole reliance on external morphology for the initial phase of taxonomic work will slow progress. We do expect that performance would have been improved if sorting had been done by an experienced gelechioid taxonomist, but this is not a general solution because many insect groups lack experts.

Conclusions
Our results affirm the general effectiveness of current algorithmic methods for species delineation together with DNA barcodes as a tool for initial species discovery. Such analyses will be particularly useful for poorly-known groups with constrained external phenotypic variation. However, we urge careful selection of methods and parameters (and starting trees where applicable) as the same approach can produce results whose quality varies depending on the focal taxon, parameter values, and distance metrics. Furthermore, a parameter value which provides a high-quality outcome in one group can generate poor results for another. The combined use of several methods provides one way to obtain a more robust estimate of species boundaries [96]. Because the focal taxa are generally poorly known in studies aiming to delineate putative species, little information on evolutionary history is usually available. Examination of the width of the barcode gap with pairwise distances without a priori grouping does provide a preliminary estimate of the levels of divergence for the group under study, a potential aid to the interpretation of results.
Some authors have indicated that DNA barcode-based methods are not useful in the absence of prior knowledge on the focal group (e.g., [91]), but we disagree. Due to its speed, simplicity and objectivity, the analysis of DNA barcode data with species delineation methods is the most feasible tool for large-scale screening of poorly-known biodiversity. It provides an accelerated start for subsequent studies which can employ broader sampling and examine more characters.
Supporting Information S1 Table. A table including