MOL E C UL AR A PPROACHE S I N F RE S HWATER ECOLOGY

Comparison of four species-delimitation methods

applied to a DNA barcode data set of insect larvae
for use in routine bioassessment
Bryan P. White1,4, Erik M. Pilgrim2,5, Laura M. Boykin3,6, Eric D. Stein1,7, and Raphael D. Mazor1,8
1
Southern California Coastal Water Research Project, Costa Mesa, California 92626 USA
2
US Environmental Protection Agency, Ecological Exposure Research Division, National Exposure Research Laboratory, Cincinnati,
Ohio 45268 USA
3
ARC Centre of Excellence, Plant Energy Biology, M316, University of Western Australia, 35 Stirling Highway, Crawley 6009,
Western Australia, Australia

Abstract: Species delimitation (grouping individuals into distinct taxonomic groups) is an essential part of evo-
lutionary, conservation, and molecular ecology. Deoxyribonucleic acid (DNA) barcodes, short fragments of the
cytochrome c oxidase subunit I (COI) gene, are being used in environmental bioassessments to assign speci-
mens to putative species, but no method for delimiting DNA barcodes into species-level entities is universally
accepted. We investigated the effect of delimitation methods on outcomes of bioassessments based on DNA
barcodes. We used 2 tree-construction methods (neighbor joining [NJ], maximum likelihood [ML]) and 4 clas-
ses of species-delimitation criteria (distance-based, bootstrap support, reciprocal monophyly, and coalescent-
based) with a DNA barcode data set consisting of 3 genera and 2202 COI sequences. We compared species de-
limitations for Baetis (Ephemeroptera:Baetidae), Eukiefferiella (Diptera:Chironomidae), and Simulium (Diptera:
Simuliidae) from different streams. We assessed congruence among trees and compared species abundances
and estimated species richness among methods. NJ followed by use of a standard barcoding distance cutoff
(2%) yielded the greatest number of putative species. All other delimitation methods yielded similar, but lower,
richness. Differences in species delimitations produced by various methods might have been caused by con-
founding factors, such as possible parthenogenesis in Baetis and rare haplotypes in abundant species of Baetis
and Simulium. Eukiefferiella presented the fewest discrepancies among delimitations. Each method can be re-
garded as producing a separate line of evidence contributing to the delimitation of separately evolving lineages.
The increased resolution offered by DNA barcoding can yield important insights into the natural history of or-
ganisms, but the power of these observations is limited without the use of multigene and multilocus data sets.
Key words: DNA barcoding, bioassessment, biomonitoring, species delimitation, coalescent theory, cytochrome
oxidase I, COI, benthic macroinvertebrates, GMYC

Deoxyribonucleic acid (DNA) barcoding is a molecular 2011, Stein et al. 2013, 2014). The conclusion is that mo-
taxonomic method wherein an ∼650-base pair (bp) region lecular methods might improve our view of stream diver-
of the cytochrome c oxidase subunit I (COI) mitochondrial sity because they provide increased taxonomic resolution.
gene is used as a species-level molecular identification In most bioassessment programs, estimates of species
marker in animals (Hebert et al. 2003). DNA barcoding richness are obtained via morphological identification.
has the potential to affect the field of bioassessment (also However, this method can be difficult and time consum-
called biomonitoring) where biological measures, such as ing, and morphology cannot be used to differentiate cryp-
species richness and taxonomic composition, are used to tic species (Bickford et al. 2007) or species complexes
draw conclusions about the health of an ecological system. (Hajibabaei et al. 2006). The ability to differentiate such
Incorporation of DNA barcode methods into bioassess- species is important because presently unknown differ-
ment programs has been evaluated theoretically (Jones ences in their tolerance to pollution, reproductive timing,
2008) and empirically (Pilgrim et al. 2011, Sweeney et al. feeding mechanisms, or other ecological traits may provide

E-mail addresses: 4bryanw@sccwrp.org; 5pilgrim.erik@epa.gov; 6lboykin@mac.com; 7erics@sccwrp.org; 8raphaelm@sccwrp.org

DOI: 10.1086/674982. Received 11 March 2013; Accepted 13 September 2013; Published online 18 December 2013.
Freshwater Science. 2014. 33(1):338–348. © 2014 by The Society for Freshwater Science.
 Volume 33 March 2014 | 339

clues regarding the health of a stream system (Verberk trees to find the highest log-likelihood (lnL) tree. We fo-
et al. 2013). Moreover, species-delimitation methods are cused on 3 insect genera that are widely encountered in
relevant to more than the outcomes of routine bioassess- freshwater bioassessment and whose species are difficult
ment programs. They may have utility in the fields of to identify morphologically: Baetis (Ephemeroptera:Bae-
criminal wildlife forensics (Dawnay et al. 2007), biodiver- tidae), Eukiefferiella (Diptera:Chironomidae), and Simu-
sity indexing (Janzen et al. 2009), detection of fish-market lium (Diptera:Simuliidae). Decisions regarding species de-
replacements (Maralit et al. 2013), ecology (Valentini et al. limitation may affect bioassessment metrics. Therefore,
2009), and biosecurity (Boykin et al. 2012). In each applica- we addressed the following questions: 1) Do different de-
tion, routine, repeated sampling of organisms should yield limitation methods yield different numbers of detectable
consistent sets of species designations across different lab- species? 2) Are differences among methods in the number
oratories. of species detected associated with the number of haplo-
Use of molecular taxonomic methods, such as DNA bar- types (i.e., are haplotype-rich species more difficult to de-
coding, to identify unknown organisms or organisms im- limit, or vice versa)? 3) Do estimates of species richness
practical to identify to species level on a routine basis has produced by different delimitation methods differ among
led to questions related to how to delimit species based sampling sites?
only on a DNA barcode. Rapid decreases in the cost of
high-throughput Sanger sequencing and the advent of M E T H O DS
Next Generation Sequencing (NGS; Shendure and Ji 2008)
Study site and genera
have led to an exponential increase in the rate at which DNA
We obtained a subset (2202) of COI sequences from
barcodes are being uploaded to public databases. Thus far,
Baetis, Eukiefferiella, and Simulium from a bioassess-
1.2 million COI sequences have been released on the Bar-
ment study of 5 streams in the Los Angeles, California
code of Life Data System (BOLD; Ratnasingham and Hebert
(USA), area (Stein et al. 2014). Benthic macroinverte-
2007; www.barcodinglife.org), as part of the International
brate samples were taken from 2 reaches at each stream.
Barcode of Life Project (iBOL, http://ibol.org/), and another
Expert taxonomists identified specimens morphologically
1 million remain uploaded but unreleased. This vast amount
using a standard level of taxonomic effort. These experts
of data has created a pressing need for efficient, objective,
identified 3 distinct species of Baetis (Baetis tricaudatus,
and readily reproducible algorithms for species delimitation.
Baetis adonis, and a 3rd unknown, but recognizably dis-
Many proponents of DNA barcoding have suggested the use
tinct, species Baetis sp. CA), whereas they identified Sim-
of measures of genetic distance to designate species. Inves-
ulium and Eukiefferiella species only to genus. We treated
tigators have recommended species limits based on an aver-
data for each genus separately. The sequences used in
age genetic distance of ≥2% (see below) among individuals
our study are publicly available under the BOLD projects
in different putative species (Ball et al. 2005, Zhou et al.
CFWIA through CFWIJ (see Table S1 for a complete list
2009) or on a level of interspecific variation that is 10× the
of BOLD sample identification codes and GenBank ac-
intraspecific variation (Hebert et al. 2004). Application of
cessions).
coalescent-theoretic methods (Knowles and Carstens 2007,
Rodrigo et al. 2008, Zaldívar-Riverón et al. 2010, Zhang et al.
2011, Nuñez et al. 2012, Vuataz et al. 2012), the principal of Sequence data and haplotype collapsing
genealogical sorting (Cummings et al. 2008), machine learn- We selected closely related genera as outgroups for
ing methods (Bertolazzi et al. 2009, Weitschek et al. 2013), a each data set (Table 1). We used a minimum sequence
heuristic-search-strategy method (O’Meara 2010), Bayesian length requirement of 500 base pairs (bp) to reduce un-
statistical methods (Yang and Rannala 2010, Zhang et al. certainty during NJ and phylogenetic analyses. We trans-
2011), and a multimethod ‘tip to root’ approach (Boykin et al. lated sequences to amino acids in MEGA (version 5.1;
2012) have been proposed as ways to increase objectivity Tamura et al. 2011) and aligned them in MUSCLE (ver-
and biological relevance of species delimitation. We eval- sion 3.8.31; Edgar 2004). The final alignment was cor-
uated 4 major classes of species-delimitation criteria (the rected manually so that sequences lacked gaps and con-
criteria by which a ‘haplotype cluster’ is granted a species- sisted of an uninterrupted open-reading-frame, which led
level status): genetic distance-based (DB), bootstrap sup- us to conclude that no insertions, deletions, or pseudo-
port (BSS), reciprocal monophyly (RM), and coalescent- genes were present in the data sets. Following alignment,
based (CB). We applied these criteria by constructing we used an open-source, custom Perl script developed for
2 types of phylogenetic trees. First, we constructed neigh- this analysis (dnab_collapse.pl, https://github.com/bpwhite
bor joining trees (NJ), which implement a clustering algo- /bioinformatics-toolbox) to reduce the number of individ-
rithm that always finds the ‘first best’ (balanced minimum ual sequences and, thus, the computational requirements
evolution) tree given the data set. Second, we constructed for each analysis. Important information, such as location,
maximum likelihood trees (ML), which implement an al- haplotype identification, and taxonomic identification were
gorithm that heuristically searches a subset of all possible preserved. The end result of this process was sequences
340 | Comparison of four species-delimitation methods B. P. White et al.

Table 1. Number of species identified in each genus with each et al. 2003, 2004, Zhou et al. 2009, 2011, Oceguera-Figueroa
species-delimitation method. ML = maximum likelihood, NJ = et al. 2010, Sweeney et al. 2011). However, Srivathsan and
neighbor-joining, CB = coalescent-based, RM = reciprocal Meier (2011) and Collins et al. (2012) recently suggested
monophyly, BSS = bootstrap support. that K2P is rarely the best nucleotide model for COI-only
Taxon ML+CB ML+RM NJ+BSS NJ+DB data sets. We applied the BSS criterion to the bootstrapped
NJ tree to define putative species based on a bootstrap sup-
Baetis 3 4 4 5 port cutoff of 95%.
Eukiefferiella 10 10 10 9
RM is a statistical approach based on the principal that
Simulium 8 7 7 10
individuals from different species will separate consistently
into distinct monophyletic clades with >95% statistical sup-
that were either unique haplotypes or haplotypes that were port. RM can be applied to either maximum parsimony
present at ≥2 sites. Each remaining sequence was automat- (MP) or ML trees, but a more rigorous test of monophyly
ically annotated with the abundance of that haplotype at (not done here) requires that the observed branching pat-
each location so that information about the diversity and tern be tested against a random branching pattern (Ro-
abundance of haplotype clusters could be garnered quickly. senberg 2007).
We implemented the ML+RM method by first identi-
fying the optimal nucleotide model for each data set with
Intraspecific-variation method (NJ+DB)
jModelTest (version 0.1.1; Posada 2008). Following nucle-
The DB criterion is based on use of an a priori genetic
otide model selection, we constructed ML phylogenetic
distance threshold as the cutoff for deciding whether 2 in-
trees using a Bayesian phylogenetic program, BEAST (ver-
dividuals are members of the same species. This criterion
sion 1.7.4; Drummond et al. 2012) with a coalescent-tree
is predicated on the idea that intraspecific genetic varia-
prior and 3 different molecular clock models: strict, relaxed
tion is small relative to interspecific variation. For exam-
lognormal, and relaxed exponential. Each clock model be-
ple, if the distance cutoff is 2% (Herbert et al. 2003, Meyer
gan with a normally distributed clock rate with a mean of
and Paulay 2005, Rivera and Currie 2009, Sweeney et al.
0.02 substitutions/million y (s/my) (Brown et al. 1979) and
2011), and the calculated genetic distance between indi-
a standard deviation of 0.005 s/my. We ran a Monte Carlo
vidual A and B is 2.5%, then the 2 individuals are assigned
Markov Chain (MCMC) simulation for 10 million steps
to different species. Variations on this method include use
and sampled trees from the MCMC at 1000-step intervals.
of average intraspecific distances (Hebert et al. 2004, Zhou
We checked parameter values for effective sample sizes
et al. 2009, 2011) or variable thresholds depending on the
(ESS) >200 and convergence by plotting marginal probabil-
taxa (Sweeney et al. 2011). The DB method typically has
ities in Tracer (version 1.5; Drummond et al. 2012). We
been applied to NJ trees computed with the algorithm of
discarded the first 20% of trees sampled as burn-ins. We
Saitou and Nei (1987) and the Kimura-2-parameter (K2P;
loaded the remaining 8001 trees into TreeAnnotator (ver-
Kimura 1980) measure of genetic distance.
sion 1.7.4; Drummond et al. 2012) for construction of the
We applied a DB criterion by calculating the nearest-
maximum clade credibility (MCC) tree and calculation of
neighbor distance (smallest interspecific distance) between
posterior probabilities and node ages. After trees were an-
haplotype clusters using the Species Delimitation (ver-
notated, we used a Bayes factor (BF) analysis to test whether
sion 1.04; Masters et al. 2011) plugin for Geneious (ver-
the data were clock-like (Drummond and Rambaut 2007).
sion 5.6.5; Biomatters, http://www.geneious.com/). Two
In this analysis, the marginal likelihoods of each tree are es-
haplotype clusters that contained a pair of nearest neigh-
timated using the harmonic mean, and the possible im-
bors with >2% K2P distance from each other were consid-
provement of one model over another is assessed by divid-
ered different putative species.
ing their marginal likelihoods when those models differ by
only 1 parameter (in this case, the clock model). The ratio
Statistical methods (NJ+BSS and ML+RM) of this division is the BF. An improvement of one model
BSS is the proportion of bootstrap replicates in which over another is considered significant if BF > 2. Following
particular sequences clustered together when the NJ al- selection of the best tree clock model, the MCC was anno-
gorithm is applied (Felsenstein 1985). For example, if a tated in FigTree (version 1.4; http://tree.bio.ed.ac.uk/soft
node achieves 95% bootstrap support, then that node ware/figtree/).
and all of its children were grouped together in 95% of
the bootstrap replicates. This method has been used in
large-scale DNA barcoding studies by Zhou et al. (2009, Modeling method (ML+CB)
2011), Lakra et al. (2011), and Mecklenburg et al. (2011). CB is a modeling approach derived from population
We implemented the NJ+BSS method using 1000 boot- genetics and is based on the principal that individuals
strap replicates in MEGA. We used K2P distance because that possess different species-level coalescent points (hypo-
it is considered ‘standard’ in DNA barcode studies (Hebert thetical ancestors of haplotypes, alleles, or species from
 Volume 33 March 2014 | 341

which point all current members of a population were de- Intraspecific variation method (NJ+DB)
scended) are members of different species. Here we con- The NJ+DB method delimited more putative Baetis
sider the COI gene tree to be analogous to the species tree, species than all other methods (Table 1) because it split
but in many cases, gene trees do not match species trees Baetis 1 into 2 species (1 and 2) (Fig. 1A, B). The NJ+DB
(Liu and Pearl 2007). Analysis of multiple genes typically is method delimited fewer putative Eukiefferiella species
required to obtain an accurate species coalescent point. than all other methods (Table 1) because it lumped Eu-
A custom R script was created (dnab_coalesce.r, also kiefferiella 5 and 6 into 1 species. The genetic distance
available from: https://github.com/bpwhite/bioinformatics between Eukiefferiella 5 and 6 was 1.9%, thus species 6
-toolbox) to run the CB species-delimitation analysis. This missed the cutoff by 0.1% (Fig. 2A, B). The NJ+DB method
script makes use of the splits package (http://r-forge.r delimited more putative Simulium species than all other
-project.org/projects/splits/), and imports the resultant methods (Table 1) because it split Simulium 1 into 2 spe-
MCC trees for each data set into the General Mixed Yule cies (1 and 2) and Simulium 9 into 2 species (9 and 10)
Coalescent (GMYC) function (gmyc). The gmyc function (Fig. 3A, B).
finds the ML threshold for the transition threshold from
a Yule process (interspecific branching rates) to a coales-
cent process (intraspecific branching rates) (Pons et al. Statistical methods (NJ+BSS and ML+RM)
2006, Fontaneto et al. 2007, Knowles and Carstens 2007, Both statistics-based methods resulted in identical spe-
Monaghan et al. 2009, Nuñez et al. 2012, Vuataz et al. cies delimitations in all 3 genera, but the use of ML tree
2012). A likelihood ratio test is automatically performed construction methods increased support values for many
to compare the coalescent model to a Yule model of evolu- nodes over the bootstrap support values (Figs 1–3).
tion, and if the ratio results in a p-value <0.05, the coales-
cent model is accepted over the Yule model and the Modeling method (ML+CB)
putative species entities can be considered statistically sig- Each data set had a different nucleotide model.
nificant. We ran the GMYC model test for a single ML Hasegawa-Kishino-Yano + gamma [HKY+G] was selected
threshold and multiple ML thresholds, wherein the thresh- for Baetis, general time reversible + gamma (GTR+G) was
old was allowed to vary across lineages. We compared the selected for Simulium, and general time reversible + in-
results of the 2 models with a χ2 goodness-of-fit test (also variant + gamma (GTR+I+G) was selected for Eukieffe-
available in the splits package under the function com- riella (Table 2). The lognormal and exponential relaxed
pare). The multiple threshold test was considered an im- clock models were not significant improvements over the
provement over the single threshold test if the data fit that strict clock model for any genus (BF < 2 in all cases), so
model significantly better (p < 0.05). After model selection, the strict clock model was used for both Baetis and Euki-
the dnab_coalesce.r script outputs the resultant species efferiella. In the case of Simulium, negative branch lengths
delimitation using the spec.list function into a comma- in the strict clock tree made the application of the GMYC
separated-value (CSV) format for import into other pro- model impossible. We used the lognormal clock tree in-
grams. stead because it had only a slightly faster mean rate than
the strict clock tree (strict: 0.199 vs lognormal: 0.244). The
Species richness and abundance calculations Simulium MCMC may have been undersampled because
We assigned 4 putative species identifications (one for the ESS for the likelihood parameter was <200, whereas the
each method) to each individual sequence and used a χ2 coalescent parameter was >200 (Table 2). For all 3 genera,
goodness-of-fit test to assess whether species richness the GMYC model was selected over the null model of evo-
was affected by the method used. We assessed the effects lution (Table 2), and the multiple ML threshold model
of species-delimitation methods on stream species rich- was not a significant improvement over the single thresh-
ness by summing the number of putative species en- old model (Table 2). The ML+CB method produced the
countered in each stream for each method. We assessed same number of species as the ML+RM and NJ+BSS
where shifts in the abundance of species might occur by methods, but the designations of those species were dif-
summing the number of individuals given a particular ferent, for example, Baetis 3 was not split into Baetis 3
species identification for each method. and 4 (Fig. 1A, B), whereas it was in the other 3 meth-
ods. Moreover, Simulium 8 and 9 were identified by the
ML+CB method, but Simulium 10 was not (Fig. 3A, B).
RESULTS
Data reduction of the 3 data sets decreased the num-
bers of sequences from 951 to 201 for Baetis, 906 to 389 Shifts in species abundance and richness
for Simulium, and 345 to 32 for Eukiefferiella. The num- The only difference in species abundances among de-
ber of putative species did not differ among delimitation limitation methods was for Baetis 1 and 2. The NJ+DB
methods for any genus (χ2Table 2 = 1.327, p > 0.05). method yielded 2 species consisting of 364 and 429 indi-
342 | Comparison of four species-delimitation methods B. P. White et al.

Figure 1. A.—Maximum likelihood, strict clock tree of 201 Baetis cytochrome c oxidase subunit I (COI) sequences computed using
BEAST with a coalescent tree prior and the Hasegawa-Kishino-Yano + gamma (HKY+G) nucleotide model. Species entities were
delimited using the gmyc function of the splits package 1.0-14. B.—Neighbor-joining tree computed using MEGA 5.05 and the Kimura
2-parameter distance nucleotide model. Putative species are numbered in order of appearance on the tree and highlighted in shaded
grey boxes. Putative species that undergo splitting are indicated by dashed boxes. A = ML+CB, B = ML+RM, C = NJ+BSS, D = NJ+DB,
* = 0.80–0.94 node support, ** = 0.95–0.99 node support, *** = 1.00 node support.

viduals, whereas the other methods yielded 1 species with crease in resolution over the standard level of identifi-
793 individuals. Minor differences in species abundances cation. The differences among species delimitations, al-
among methods were present in all 3 genera, but these though not statistically significant, tended to be associated
differences were limited mostly to the presence or absence with abundant and diverse taxa. This result suggests that
of a few rare haplotypes. Species richness differed consis- the uncertainty associated with species delimitations de-
tently among methods and sites (Fig. 4). When estimates rived from DNA barcoding does not arise from the algo-
differed among methods, the NJ+DB method typically rithm used, but is a byproduct of the inherent limitations
produced higher species richness than the other methods. of COI as a species-level phylogenetic marker. Differences
In all cases the RM and BSS method produced identical among delimitation methods are not likely to result in large
abundance and richness estimates. changes in bioassessment metric scores based on taxon rich-
ness.
The absence of noticeable differences in species abun-
DISCUSSION dances under different delimitation methods (except Bae-
Four species-delimitation methods applied to a data tis 1 and 2) suggests that richness metrics that take
set of 2202 COI sequences from 3 genera of insect larvae abundances into account might be unaffected by different
from southern California yielded similar estimates of delimitation methods, whereas presence–absence-type rich-
species richness. Where standard morphological identifi- ness metrics may be more directly influenced by even sub-
cation effort yielded 5 distinct taxa (Baetis adonis, Baetis tle shifts in species designations (e.g., Simulium 2 and 10
tricaudatus, Baetis sp. CA, Simulium, and Eukiefferiella), exist only under the NJ+DB method and have extremely
DNA barcodes yielded 19 to 25 putative species, a 4× in- low abundances [1 and 2 individuals, respectively]). How-
 Volume 33 March 2014 | 343

Figure 2. A.—Maximum likelihood, strict clock tree of 32 Eukiefferiella cytochrome c oxidase subunit I (COI) sequences
computed with the GTR+I+G nucleotide model. B.—Neighbor-joining tree. See Fig. 1A, B for details of species delimitation, tree
construction, and figure conventions.

ever, even those small differences may not be large enough uals), of increased diversity with increased intraspecific
to significantly affect bioassessment metrics beyond the sampling effort, a result that further supports the hypothe-
greater taxonomic resolution already provided by DNA bar- sis of Bergsten et al. (2012). In our data set, Baetis 1 and
coding. 2 had an average intraspecific K2P distance of 1.6% and
Species delimitations that differed among methods a maximum pairwise K2P distance of 5.1% when lumped
tended to be associated with high-abundance species with together according to the NJ+BS, ML+RM, and ML+CB
many different COI haplotypes, a potentially important criteria. The less frequently encountered Baetis species
trend warranting further examination. For example, pu- (3, 4, and 5), which made up only 17% of Baetis encoun-
tative species that experienced splits under the NJ+DB tered (158 individuals), did not exhibit intraspecific pair-
method (Baetis 1–2, 3–4, and Simulium 1–2, 9–10) were wise distances >1%. These differences in sampling effort
usually very abundant (>30 individuals encountered). Large reflect natural abundances of these species and not targeted
populations tend to be less prone to chance events like ge- sampling effort toward one species or another. In contrast,
netic drift that eliminate rare haplotypes (given enough we also found many individuals in Baetis 1 and 2 that were
time), and it seems plausible that very abundant species separated by great geographic distances but shared iden-
will maintain more rare haplotypes in their gene pools than tical haplotypes. For example, 8 individuals found in Co-
will less abundant species. This idea received support nejo Creek had haplotypes identical to 50 individuals in
from Bergsten et al. (2012), who concluded that the uncer- Big Tujunga Wash. The 2 streams are separated by a geo-
tainty of species identifications for Agabini diving beetles graphic distance of 64 km and an elevation difference of
increased significantly when sampled over increasing geo- 335 m (see Stein et al. 2014 for stream descriptions). We
graphic distances. They found that a sample size of 70 in- did not see more than 1 or 2 instances of shared haplo-
dividuals was necessary to capture 95% of intraspecific types over great distances in Simulium or Eukiefferiella.
diversity. We observed a pattern in Baetis sp. 1 and 2, The multimodal pattern of genetic variation within
which made up 83% of Baetis encountered (793 individ- the Baetis 1 and 2 complex consists of a mixture of low
344 | Comparison of four species-delimitation methods B. P. White et al.

Figure 3. A.—Maximum likelihood, lognormal clock tree of 389 Simulium cytochrome c oxidase subunit I (COI) sequences
computed with the GTR+G nucleotide model. B.—Neighbor-joining tree. See Fig. 1A, B for details of species delimitation, tree
construction, and figure conventions.

and high diversity over broad distances. Such a pattern zoan mitochondrial inheritance and that shared haplotypes
might be explained if these species were parthenogenetic. are a consequence of evolutionarily recent dispersal events.
Many mayfly species exhibit parthenogenesis (Bergman The increased taxonomic resolution offered by DNA
and Hilsenhoff 1978, McCafferty and Morihara 1979, Funk barcoding could be used to understand better the life histo-
et al. 2006, 2008), and Funk et al. (2010) suggested that ries of the seemingly cryptic Baetis observed in our study,
most, if not all, mayflies may be facultatively parthenoge- which might in turn help researchers use traits-based ap-
netic. Females may reproduce parthenogenetically, but proaches to ecology and improve future bioassessment
offspring can be either male or female and can readily re- tools (Verberk et al. 2013). In the case of Baetis spp. in
vert to sexual reproduction. Vanoverbeke and De Meester southern California, the addition of nuclear loci to DNA
(1997) found no relationship between geographic distance barcode data sets might help researchers distinguish be-
and genetic distance in parthenogenetically reproducing tween normal and parthenogenetic modes of inheritance
populations of cladoceran branchiopods (Daphnia magna), (Buckley et al. 2008). When systematists undertake taxo-
a result similar to our finding of shared haplotypes over nomic revisions of morphologically cryptic species in light
large distances (64 km). However, we cannot rule out the of molecular data (e.g., in the Atyaephyra genus of fresh-
possibility that this pattern is a by-product of normal meta- water shrimp; Christodoulou et al. 2012), specific life-
Table 2. Taxa used, number of cytochrome c oxidase subunit I (COI) sequences, number of haplotypes for each taxa, outgroups used for each taxa, Barcode of Life Database
(BOLD) sample identification numbers (IDs), GenBank accession numbers for each outgroup, nucleotide model used, tree log likelihood (lnL), effective sample size (ESS) of tree lnL
and coalescent BEAST parameters, null (Yule) model lnL from gmyc function, General Mixed Yule Coalescent (GMYC) model lnL, p-value of the χ2 goodness-of-fit test between
the null and GMYC model, p-value of the comparison between the single and multiple (multi) GMYC threshold, and number of putative species detected under the GMYC model
(ML entities). Mya = million years ago.
GenBank Null GMYC
COI BOLD accession Nucleotide ESS of ESS of model model Multi Threshold ML
Taxon sequences Haplotypes Outgroups sample ID number model Tree lnL tree lnL coalescent lnL lnL GMYC p p (mya) entities

Baetis 906 104 Fallceon sp. 10-SCCWRP-6041 JN297857 HKY+G −3987.374 203.5 963.9 1195.81 1220.916 <<0.0001 0.884 0.973 3
Centroptilum 09LJ17 HM423544
triangulifer
Callibaetis 09NBMAY-0053 JQ663249
ferrugineus
Simulium 951 304 Prosimulium 08-SWRC-1082 JF287830 GTR+G –6018.239 22.5 368.4 2258.584 2312.053 <<0.0001 0.917 1.063 7
mixtum
Prosimulium FJ524559 FJ524559
travisi
Eukiefferiella 345 32 Chironomus ChiBlank_F05 HM137935 GTR+I+G −3826.534 857.1 3744.5 128.73 145.469 <<0.0001 0.999 0.333 10
riparius
Chironomus JQ350720 JQ350720
kiiensis
Chironomus CAUS016-09 –
jacksoni
346 | Comparison of four species-delimitation methods B. P. White et al.

multigene data sets in routine bioassessments might be

cost prohibitive because linking multiple genes to a single
voucher specimen requires sorting and molecular tagging
of each individual voucher specimen. As metabarcoding
of environmental samples becomes more prevalent (Carew
et al. 2013, Ji et al. 2013), sorting of individual organisms
might become less problematic. One solution might be to
maintain multigene reference libraries of local taxa identi-
fied with rigorous species-delimitation methods, but to use
single-gene methods during routine field sampling (e.g.,
COI- or 16S-based DNA barcoding coupled with distance-
based algorithms). Multigene reference libraries would al-
low researchers to know a priori when instances of mito-
chondrial introgression (incorrect lumping) or ancestral
polymorphisms (incorrect splitting) might produce spuri-
ous single-gene delimitations so that degrees of confidence
could be assigned to molecularly derived species delimi-
tations.

AC KNOW LE DGEMENTS
Figure 4. Differences in species richness observed at
We thank Tomochika Fujisawa for providing an advance re-
different southern California stream reaches.
lease copy of the Species List script essential to the coalescent
delimitation analysis. We are also thankful to Daniel Pickard for
providing taxonomic expertise and valuable feedback on Baetis
history traits, such as parthenogenesis, could be included
species identifications. We are grateful to the Canadian Centre for
in descriptions, transferred to traits databases, and associ- DNA Barcoding, which provided sequencing with support from
ated with DNA barcodes. the Life Technologies Corporation. We also thank the 2 anony-
Coalescent models take into account the natural birth mous referees for their feedback which greatly improved the qual-
and death processes of populations, avoid the use of a ity of the manuscript.
priori distance cutoffs (e.g., the 2% cutoff ), and provide a
statistical framework for testing species delimitations. Fur-
thermore, the GMYC model provides estimates of the L I T E R ATU R E C I T E D
times of transitions from inter- to intraspecies branching Ball, S. L., P. D. N. Hebert, S. K. Burian, and J. M. Webb. 2005.
patterns and effectively links the fields of population genet- Biological identifications of mayflies (Ephemeroptera) using
ics and phylogenetics. However, multiple genes from both DNA barcodes. Journal of the North American Benthological
nuclear and mitochondrial DNA are required to obtain a Society 24:508–524.
Bergman, E. A., and W. L. Hilsenhoff. 1978. Parthenogenesis in
robust estimate of the species coalescence (Heled and
the mayfly genus Baetis (Ephemeroptera: Baetidae). Annals
Drummond 2010, Fujita et al. 2012). If we regard a species
of the Entomological Society of America71:167–168.
as any separately evolving metapopulation of lineages (de Bergsten, J., D. T. Bilton, T. Fujisawa, M. Elliott, M. T. Monaghan,
Queiroz 2007), the ML+CB method stands beside distance M. Balke, and A. P. Vogler. 2012. The effect of geographical
and monophyly-based methods as a separate and unique scale of sampling on DNA barcoding. Systematic Biology 61:
line of evidence in the diagnosis of distinct lineages. Thus, 851–869.
combining the results of several species-delimitation meth- Bertolazzi, P., G. Felici, and E. Weitschek. 2009. Learning to
ods might allow researchers to draw confident conclusions classify species with barcodes. BMC Bioinformatics10:S7.
when diagnosing a lineage (as in Boykin et al. 2012). Bickford, D., D. J. Lohman, N. S. Sodhi, P. K. Ng, R. Meier,
Differentiating between cases of rare, divergent haplo- K. Winker, K. K. Ingram, and I. Das. 2007. Cryptic species
types within a species and cryptic species is one of the as a window on diversity and conservation. Trends in Ecol-
ogy and Evolution 22:148–155.
greatest challenges for users of DNA barcodes as species-
Boykin, L. M., K. F. Armstrong, L. Kubatko, and P. D. Barro.
level markers. This challenge probably will be overcome
2012. Species delimitation and global biosecurity. Evolution-
only by using multigene and multilocus reference librar- ary Bioinformatics 8:1–37.
ies. A combination of mitochondrial genes (e.g., COI, Brown, W. M., M. George, and A. C. Wilson. 1979. Rapid evo-
CYTB, and 16S) and nuclear genes (e.g., 18S, 28S, ITS1, lution of animal mitochondrial DNA. Proceedings of the Na-
and ITS2) would allow robust estimates of species coa- tional Academy of Sciences of the Unites States of America
lescence and improved phylogenetic resolution. Use of 76:1967–1971.
 Volume 33 March 2014 | 347

Buckley, T. R., D. Attanayake, D. Park, S. Ravindran, T. R. Jewell, Hajibabaei, M., D. H. Janzen, J. M. Burns, W. Hallwachs, and
and B. B. Normark. 2008. Investigating hybridization in the P. D. N. Hebert. 2006. DNA barcodes distinguish species of
parthenogenetic New Zealand stick insect Acanthoxyla (Phas- tropical Lepidoptera. Proceedings of the National Academy
matodea) using single-copy nuclear loci. Molecular Phylo- of Sciences of the United States of America 103:968–971.
genetics and Evolution 48:335–349. Hebert, P. D. N., A. Cywinska, S. L. Ball, J. R. deWaard. 2003.
Carew, M. E., V. J. Pettigrove, L. Metzeling, and A. A. Hoff- Biological identifications through DNA barcodes. Proceedings
mann. 2013. Environmental monitoring using next genera- of the Royal Society of London Series B: Biological Sciences
tion sequencing: rapid identification of macroinvertebrate 270:313–321.
bioindicator species. Frontiers in Zoology 10:1–15. Hebert, P. D. N., M. Y. Stoeckle, T. S. Zemlak, and C. M.
Christodoulou, M., A. Antoniou, A. Magoulas, and A. Koukouras. Francis. 2004. Identification of birds through DNA barcodes.
2012. Revision of the freshwater genus Atyaephyra (Crustacea, PLoS Biology 2:e312.
Decapoda, Atyidae) based on morphological and molecular Heled, J., and A. J. Drummond. 2010. Bayesian inference of
data. ZooKeys 229:53–110. species trees from multilocus data. Molecular Biology and
Collins, R. A., L. M. Boykin, R. H. Cruickshank, and K. F. Evolution 27:570–580.
Armstrong. 2012. Barcoding’s next top model: an evaluation Janzen, D. H., W. Hallwachs, P. Blandin, J. M. Burns, J. Cadiou, I.
of nucleotide substitution models for specimen identifica- Chacon, T. Dapkey, A. R. Deans, M. E. Epstein, B. Espinoza,
tion. Methods in Ecology and Evolution 3:457–465. J. G. Franclemont, W. A. Haber, M. Hajibabaei, J. P. W. Hall,
Cummings, M. P., M. C. Neel, and K. L. Shaw. 2008. A genea- P. D. N. Hebert, I. D. Gauld, D. J. Harvey, A. Hausmann, I. J.
logical approach to quantifying lineage divergence. Evolu- Kitching, D. Lafontaine, J.-F. Landry, C. Lemaire, J. Y. Miller,
tion 62:2411–2422. J. S. Miller, L. Miller, S. E. Miller, J. Montero, E. Munroe,
Dawnay, N., R. Ogden, R. McEwing, G. R. Carvalho, and R. S. S. R. Green, S. Ratnasingham, J. E. Rawlins, R. K. Robbins, J. J.
Thorpe. 2007. Validation of the barcoding gene COI for use Rodriguez, R. Rougerie, M. J. Sharkey, M. A. Smith, M. A.
in forensic genetic species identification. Forensic Science Solis, J. B. Sullivan, P. Thiaucourt, D. B. Wahl, S. J. Weller,
International 173:1–6. J. B. Whitfield, K. R. Willmott, D. M. Wood, N. E. Woodley,
de Queiroz, K. 2007. Species concepts and species delimitation. and J. J. Wilson. 2009. Integration of DNA barcoding into an
Systematic Biology 56:879–886. ongoing inventory of complex tropical biodiversity. Molecu-
Drummond, A. J., and A. Rambaut. 2007. BEAST: Bayesian lar Ecology Resources 9:1–26.
evolutionary analysis by sampling trees. BMC Evolutionary Ji, Y., L. Ashton, S. M. Pedley, D. P. Edwards, Y. Tang, A. Nakamura,
Biology 7:214. R. Kitching, P. M. Dolman, P. Woodcock, F. A. Edwards, T. H.
Drummond, A. J., M. A. Suchard, D. Xie, and A. Rambaut. Larsen, W. W. Hsu, S. Benedick, K. C. Hamer, D. S. Wilcove, C.
2012. Bayesian phylogenetics with BEAUti and the BEAST Bruce, X. Wang, T. Levi, M. Lott, B. C. Emerson, and D. W. Yu.
1.7. Molecular Biology and Evolution 29:1969–1973. 2013. Reliable, verifiable and efficient monitoring of biodiver-
Edgar, R. C. 2004. MUSCLE: multiple sequence alignment with sity via metabarcoding. Ecology Letters 16:1245–1257.
high accuracy and high throughput. Nucleic Acids Research Jones, F. C. 2008. Taxonomic sufficiency: the influence of taxo-
32:1792–1797. nomic resolution on freshwater bioassessments using ben-
Felsenstein, J. 1985. Confidence limits on phylogenies: an ap- thic macroinvertebrates. Environmental Reviews 16:45–69.
proach using the bootstrap. Evolution 39:783–791. Kimura, M. 1980. A simple method for estimating evolutionary
Fontaneto, D., E. A. Herniou, C. Boschetti, M. Caprioli, rate of base substitutions through comparative studies of nuc-
G. Melone, C. Ricci, and T. G. Barraclough. 2007. Indepen- leotide sequences. Journal of Molecular Evolution 16:111–120.
dently evolving species in asexual bdelloid rotifers. PloS Bi- Knowles, L. L., and B. C. Carstens. 2007. Delimiting species with-
ology 5:e87. out monophyletic gene trees. Systematic Biology 56:887–895.
Fujita, M. K., A. D. Leache, F. T. Burbrink, J. A. McGuire, and Lakra, W. S., M. S. Verma, M. Goswami, K. K. Lal, V.
C. Moritz. 2012. Coalescent-based species delimitation in an Mohindra, P. Punia, A. Gopalakrishnan, K. V. Singh, R. D.
integrative taxonomy. Trends in Ecology and Evolution 27: Ward, and P. Hebert. 2011. DNA barcoding Indian marine
480–488. fishes. Molecular Ecology Resources 11:60–71.
Funk, D. H., J. K. Jackson, and B. W. Sweeney. 2006. Taxonomy Liu, L., and D. K. Pearl. 2007. Species trees from gene trees:
and genetics of the parthenogenetic mayfly Centroptilum tri- reconstructing Bayesian posterior distributions of a species
angulifer and its sexual sister Centroptilum alamance (Ephe- phylogeny using estimated gene tree distributions. Systematic
meroptera:Baetidae). Journal of the North American Bentho- Biology 56:504–514.
logical Society 25:417–429. Maralit, B. A., R. D. Aguila, M. F. H. Ventolero, S. K. L. Perez,
Funk, D. H., J. K. Jackson, and B. W. Sweeney. 2008. A new par- and M. D. Santos. 2013. Detection of mislabeled commer-
thenogenetic mayfly (Ephemeroptera: Ephemerellidae: Eurylo- cial fishery by-products in the Philippines using DNA bar-
phella Tiensuu) oviposits by abdominal bursting in the sub- codes and its implications to food traceability and safety. Food
imago. Journal of the North American Benthological Society Control (in press).
27:269–279. Masters, B. C., V. Fan, and H. A. Ross. 2011. Species delimitation–
Funk, D. H., B. W. Sweeney, and J. K. Jackson. 2010. Why a Geneious plugin for the exploration of species boundaries.
stream mayflies can reproduce without males but remain bi- Molecular Ecology Resources 11:154–157.
sexual: a case of lost genetic variation. Journal of the North McCafferty, W. P., and D. K. Morihara. 1979. The male of
American Benthological Society 29:1258–1266. Baetis macdunnoughi Ide and notes on parthenogenetic pop-
348 | Comparison of four species-delimitation methods B. P. White et al.

ulations within Baetis (Ephemeroptera: Baetidae). Entomo- Stein, E. D., B. P. White, R. D. Mazor, P. E. Miller, and E. M.
logical News 90:26–28. Pilgrim. 2013. Evaluating ethanol-based sample preservation
Mecklenburg, C. W., P. R. Møller, and D. Steinke. 2011. Biodi- to facilitate use of DNA barcoding in routine freshwater
versity of arctic marine fishes: taxonomy and zoogeography. biomonitoring programs using benthic macroinvertebrates.
Marine Biodiversity 41:109–140. PloS ONE 8:e51273.
Meyer, C. P., and G. Paulay. 2005. DNA barcoding: error rates Sweeney, B. W., J. M. Battle, J. K. Jackson, and T. Dapkey.
based on comprehensive sampling. PLoS Biology 3:e422. 2011. Can DNA barcodes of stream macroinvertebrates im-
Monaghan, M. T., R. Wild, M. Elliot, T. Fujisawa, M. Balke, D. J. prove descriptions of community structure and water qual-
Inward, D. C. Lees, R. Ranaivosolo, P. Eggleton, T. G. Bar- ity? Journal of the North American Benthological Society 30:
raclough, and A. P. Vogler. 2009. Accelerated species inven- 195–216.
tory on Madagascar using coalescent-based models of species Tamura, K., D. Peterson, N. Peterson, G. Stecher, M. Nei, and
delineation. Systematic Biology 58:298–311. S. Kumar. 2011. MEGA5: Molecular Evolutionary Genetics
Nuñez, J. J., A. Vejar-Pardo, B. E. Guzmán, E. H. Barriga, and Analysis using maximum likelihood, evolutionary distance,
C. S. Gallardo. 2012. Phylogenetic and mixed Yule-coalescent and maximum parsimony methods. Molecular Biology and
analyses reveal cryptic lineages within two South American Evolution 28:2731–2739.
marine snails of the genus Crepipatella (Gastropoda: Calyp- Valentini, A., F. Pompanon, and P. Taberlet. 2009. DNA bar-
traeidae). Invertebrate Biology 131:301–311. coding for ecologists. Trends in Ecology and Evolution 24:
Oceguera-Figueroa, A., V. León-Règagnon, and M. E. Siddall. 110–117.
2010. DNA barcoding reveals Mexican diversity within the Vanoverbeke, J., and L. De Meester. 1997. Among-populational
freshwater leech genus Helobdella (Annelida: Glossiphoni- genetic differentiation in the cyclical parthenogen Daph-
idae). Mitochondrial DNA 21:24–29. nia magna (Crustacea, Anomopoda) and its relation to geo-
O’Meara, B. C. 2010. New heuristic methods for joint species graphic distance and clonal diversity. Hydrobiologia 360:
delimitation and species tree inference. Systematic Biology 135–142.
59:59–73. Verberk, W. C. E. P., C. G. E. van Noordwijk, and A. G.
Pilgrim, E. M., S. A. Jackson, S. Swenson, I. Turcsanyi, E. Hildrew. 2013. Delivering on a promise: integrating species
Friedman, L. Weigt, and M. J. Bagley. 2011. Incorporation of traits to transform descriptive community ecology into a
DNA barcoding into a large-scale biomonitoring program: predictive science. Freshwater Science 32:531–547.
opportunities and pitfalls. Journal of the North American Vuataz, L., M. Sartori, J. L. Gattolliat, and M. T. Monaghan.
Benthological Society 30:217–231. 2012. Endemism and diversification in freshwater insects of
Pons, J., T. G. Barraclough, J. Gomez-Zurita, A. Cardoso, D. P. Madagascar revealed by coalescent and phylogenetic analysis
Duran, S. Hazell, and A. P. Vogler. 2006. Sequence-based of museum and field collections. Molecular Phylogenetics
species delimitation for the DNA taxonomy of undescribed and Evolution 66:979–991.
insects. Systematic Biology 55:595–609. Weitschek, E., R. Velzen, G. Felici, and P. Bertolazzi. 2013.
Posada, D. 2008. jModelTest: phylogenetic model averaging. BLOG 2.0: a software system for character-based species
Molecular Biology and Evolution 25:1253–1256. classification with DNA Barcode sequences. What it does,
Ratnasingham, S., and P. D. Hebert. 2007. BOLD: the Barcode how to use it. Molecular Ecology Resources. doi:10.1111
of Life Data System (http://www.barcodinglife.org). Molecu- /1755-0998.12073
lar Ecology Notes 7:355–364. Yang, Z., and B. Rannala. 2010. Bayesian species delimitation
Rivera, J., and D. C. Currie. 2009. Identification of Nearctic using multilocus sequence data. Proceedings of the National
black flies using DNA barcodes (Diptera: Simuliidae). Mo- Academy of Sciences of the United States of America 107:
lecular Ecology Resources 9:224–236. 9264–9269.
Rodrigo, A., F. Bertels, J. Heled, R. Noder, H. Shearman, and P. Zaldívar-Riverón, A., J. J. Martínez, F. S. Ceccarelli, V. S. De
Tsai. 2008. The perils of plenty: what are we going to do with Jesús-Bonilla, A. C. Rodríguez-Pérez, A. Reséndiz-Flores,
all these genes? Philosophical Transactions of the Royal So- and M. A. Smith. 2010. DNA barcoding a highly diverse
ciety of London Series B: Biological Sciences 363:3893–3902. group of parasitoid wasps (Braconidae: Doryctinae) from a
Rosenberg, N. A. 2007. Statistical tests for taxonomic distinctive- Mexican nature reserve. Mitochondrial DNA 21:18–23.
ness from observations of monophyly. Evolution 61:317–323. Zhang, C., D. X. Zhang, T. Zhu, and Z. Yang. 2011. Evaluation
Saitou, N., and Nei M. 1987. The neighbor-joining method: a of a Bayesian coalescent method of species delimitation. Sys-
new method for reconstructing phylogenetic trees. Molecu- tematic Biology 60:747–761.
lar Biology and Evolution 4:406–425. Zhou, X., S. J. Adamowicz, L. M. Jacobus, E. DeWalt, and P. D. N.
Shendure, J., and H. Ji. 2008. Next-generation DNA sequenc- Hebert. 2009. Toward a comprehensive barcode library for
ing. Nature Biotechnology 26:1135–1145. Arctic life–Ephemeroptera, Plecoptera, and Trichoptera of
Srivathsan, A., and R. Meier. 2011. On the inappropriate use of Churchill, Manitoba, Canada. Frontiers in Zoology 6:30.
Kimura-2-parameter (K2P) divergences in the DNA-barcoding Zhou, X., J. L. Robinson, C. J. Geraci, C. R. Parker, O. S. Flint, D.
literature. Cladistics 28:190–194. Etnier, D. Ruiter, R. E. DeWalt, L. M. Jacobus, and P. D. N.
Stein, E. D., B. P. White, R. D. Mazor, J. K. Jackson, J. M. Battle, Hebert. 2011. Accelerated construction of a regional DNA
P. E. Miller, E. M. Pilgrim, and B. W. Sweeney. 2014. Does barcode reference library: caddisflies (Trichoptera) in the
DNA barcoding improve performance of traditional stream Great Smoky Mountains National Park. Journal of the North
bioassessment metrics? Freshwater Science 33:302–311. American Benthological Society 30:131–162.