Integrative taxonomy by molecular species delimitation: multi-locus data corroborate a new species of Balkan Drusinae micro-endemics

Taxonomy offers precise species identification and delimitation and thus provides basic information for biological research, e.g. through assessment of species richness. The importance of molecular taxonomy, i.e., the identification and delimitation of taxa based on molecular markers, has increased in the past decade. Recently developed exploratory tools now allow estimating species-level diversity in multi-locus molecular datasets. Here we use molecular species delimitation tools that either quantify differences in intra- and interspecific variability of loci, or divergence times within and between species, or perform coalescent species tree inference to estimate species-level entities in molecular genetic datasets. We benchmark results from these methods against 14 morphologically readily differentiable species of a well-defined subgroup of the diverse Drusinae subfamily (Trichoptera, Limnephilidae). Using a 3798 bp (6 loci) molecular data set we aim to corroborate a geographically isolated new species by integrating comparative morphological studies and molecular taxonomy. Our results indicate that only multi-locus species delimitation provides taxonomically relevant information. The data further corroborate the new species Drusus zivici sp. nov. We provide differential diagnostic characters and describe the male, female and larva of this new species and discuss diversity patterns of Drusinae in the Balkans. We further discuss potential and significance of molecular species delimitation. Finally we argue that enhancing collaborative integrative taxonomy will accelerate assessment of global diversity and completion of reference libraries for applied fields, e.g., conservation and biomonitoring.


Background
Species represent a fundamental information unit in biological research [1][2][3]. Species-specific abundance fluctuations integrated with autecological attributes are used to assess stream health and ecological water quality, evaluate the potential of disease and parasite vectors, and were found to be highly informative in species distribution modelling approaches [4][5][6][7][8][9][10][11][12][13][14][15]. Aggregate taxa, i.e., taxonomic entities comprising more than one species, often do not provide sufficient resolution to reap the power of ecological analysis [16][17][18]. Thus, estimation of, e.g., ecological water quality in compliance with the EU Water Framework Directive crucially depends on precise taxonomy to delineate and define species [16].
Given the taxonomic impedimentthe worldwide decline of taxonomic competence to identify species based on morphological charactersbiological sciences and policymaking are severely hampered by difficulties in compiling relevant and up-to-date baseline diversity data [19][20][21]. Indeed, taxonomy and assessment of eukaryotes remains primarily reliant on identification and comparison of morphological characters to define and address species [22][23][24][25]. However, characterization of species ideally uses different sources of information (ecological, morphological, anatomical, physiological, genomic, geographical or others) in an integrative taxonomic approach for the identification, delineation and description of taxa [2,3,26,27]. Advances in molecular genetic methods recently promoted molecular taxonomy: species recognition and delineation based on unique genomic characters [28][29][30].
Also, despite the demonstrated informativeness of molecular taxon delimitation to test species hypotheses, not all studies that successfully employ molecular taxonomic taxa delimitation tools follow through to describe these new taxonomic entities [3,45,46]. Interestingly, this has been related to the complex taxonomic procedures associated with the formal description of new species [45]. However it is more likely that molecular species hypotheses are ignored in integrative approaches due to insufficient morphological or ecological support (e.g., [47]).
In this contribution we benchmark results of molecular species delimitation against morphologically well-defined taxa in a highly diverse group of caddisflies. The subfamily Drusinae (Insecta, Trichoptera) constitutes an ideal model taxon to assess potential congruence of traditional and molecular taxonomic methods. This group of mostly coldstenotopic species inhabiting Eurasian mountain ranges exhibits highly disjunct distribution patterns and high levels of micro-endemism, indicative of small-scale allopatric diversification and persistence of isolated lineages over geological time [48][49][50][51][52]. However, historic introgressiona process that complicates species delimitationwas demonstrated in some species of Drusinae [49]. In the Western Balkans, taxonomic richness of Drusinae is particularly high [50,51,[53][54][55][56][57][58] and was presumably shaped by multiple glacial cycles and karstification (cf. [48,51,59]). Additionally, Western Balkan Drusinae are morphologically distinct with multi-locus molecular data showing minimal differences which potentially result from recent speciation [51,59]. Together, these conditions make Drusinae a good model for testing the suitability and precision of species delimitation methods.
Here, we assess taxonomic informativeness of four recently developed exploratory molecular species delimitation tools by inferring entity richness hypotheses on a 3798 bp, 6 loci (mtCOI5-P, mtCOI3-P, 16S mrDNA, CADH, WG, 28S nrDNA), 14 morphospecies dataset comprising 1 new and 1 recently described species. As a test case, we aim to clarify the systematic status of a morphologically distinct potential new species and corroborate its distinctiveness in an integrative taxonomic approach. We predict that only methods directly integrating information from several loci will provide taxonomically conclusive results.

Collection and taxonomic methods
Adult specimens were collected using sweep nets, larvae were collected by hand-picking. Collected specimens were stored in 96% EtOH. Specimens were cleared for genitalic dissections and examinations using either the Qiagen Blood and Tissue Kit for DNA-extraction according to the manufacturer's recommendation and subsequent KOHtreatment [60], or KOH-treatment. Nomenclature of male terminalia follows [61] (for Limnephilus flavicornis Fabricius) using the simplifying terms "superior appendages" for the lateral processes of segment X (cerci sensu [62]), and "intermediate appendages" for the sclerite and the anterior process of segment X (paraproct sensu [62]). Illustrations were prepared according to [63] in which pencil drawings made with a camera lucida are digitized, edited and inked in Adobe Illustrator (v. 16.0.4, Adobe Systems Inc.).
Morphology-based delimitation of species-level taxa was achieved in a classical comparative taxonomic approach: we scrutinized as many specimens as possible from as many populations as possible from as many Drusinae species as possible to discriminate intraspecific from interspecific variation in more than 300 adult and larval features (e.g., structure of male and female terminalia, wing venation, larval feeding apparatus, and larval pronotum shape; [64]) and unequivocally assigned specimens to species-level taxa (cf. [50,[55][56][57][58]64]).

Phylogenetic inference and molecular species delimitation
Nucleotide substitution models for each partition were selected according to the Bayesian Information Criterion in the model test module of Mega v6 [71] (Table 1). For phylogenetic analysis, the mt16S and nu28S fragments were not partitioned; protein coding genes (mtCOI5-P, mtCOI3-P, nuCADH, nuWnt1) were partitioned by codon position.
Single gene phylogenies were estimated by Bayesian Inference through BEAST 2 [72] (5 × 10 9 generations, sampling every 10,000th generation). Analyses were run 4× independently to assure topological convergence. BEAST log files were examined in Tracer v1.6 [73] to assess when runs had reached a stationary phase. A maximum clade credibility tree was estimated via TreeAnnotator v1.8.1 [74] based on the sampled trees after discarding the first 30% as burn-in. Also, congruence of phylogenetic signal among data partitions was assessed by examining ≥0.95 posterior probability topologies of single gene analyses.
A species trees was estimated using *BEAST [75] as implemented in BEAST 2 using unique specimen identifiers as the species trait, i.e. without a priori species definitions. We ran species tree analysis assuming a Yule speciation tree prior for 5 × 10 8 generations 4× independently, sampling every 10,000th generation. *BEAST log files were examined in Tracer v1.6 to assess if runs had reached a stationary phase and converged on model parameters; maximum clade credibility trees were estimated as described above.
We then performed molecular species delimitation using tools complying with the following set of criteria: (1) the method was designed as naive exploratory tool without a priori assignment of specimens to groups or assumptions about relationships between specimens (we consequently excluded BPP [76], as this method requires at least an a priori clustering of specimens into groups and a guide tree); (2) the method was originally designed for molecular species delimitation (we consequently excluded applications like Brownie [32,77], or SpeDeSTEM [78][79][80]); (3) the method was designed to exploit 'standard' partial genetic sequence data.
To assess potential congruence of molecular species delimitation methods, we performed a series of analyses employing several analytical resources selected as described above and report results against the benchmark of 14 morphologically identifiable species: (I) Automatic Barcode Gap Detection [ABGD] was performed for each locus and a concatenated sequence data set via the ABGD webmask [81]. ABGD is a tool designed to infer species hypotheses based on automatized identification of barcode gaps between inter-and intraspecific pairwise distances in partial sequence data sets. The method does not make assumptions about data structure or evolutionary history, and only requires input data (a single locus alignment) to be sufficiently variable. Pairwise distances are computed either as simple p-distances, or as substitution-corrected distances (via either JC69 [82], or K2P [83] models). Barcode gaps are discovered as slope maxima of a function describing the relation of pairwise distance ranks and pairwise distances. The method partitions specimens into groups in a recursive manner in which every group is split again, until no further splits are possible, while integrating user-provided priors on maximum and minimum intraspecific differentiation and barcode gap width. The prior on intraspecific divergence (denoted P in the original publication and the software [81]) defines the threshold between intra-and interspecific pairwise distances, and is iterated from minimum to maximum  [81]) defines sensitivity of the algorithm as scaling factor of empiric maximum intraspecific divergence to estimate minimum barcode gap width between intra-and interspecific distances [by default X = 1.5]. From the thus provided set of possible partitions (up to 10 using default settings) Puillandre et al. [81] suggest to select the most plausible one(s) and assess their potential informativeness in an integrative approach. To avoid oversplitting of single species, we used the default distance metric (JC69) as this was previously found to produce more conservative species hypotheses [84,85]. Likewise, we used the default minimum barcode gap width prior to derive relatively conservative species hypothesesmodifications (i.e., increasing or decreasing the default value) thereof were previously reported to either increase or decrease numbers of species hypotheses proposed, with smaller X values generally leading to more delimited entities and vice versa [84,85]. Assuming that the potential barcode gap space would be satisfyingly described, we used default intraspecific divergence minima and maxima, and report delimitations at partition maximum as identified through recursive partitioning. Following the suggestions by Puillandre et al. [81], we only present species hypotheses at partition maxima as these correspond most closely to expected numbers of taxa.
(II) We used the GMYC [Generalized Mixed Yule Coalescent] model [86,87], implementing single and multiple thresholds via the 'splits' package in R 3.2.1 [88,89] on single gene trees and a *BEAST species tree to infer GMYC species. The GMYC model aims to discern stochastic birth-death processes (effectively a pure-birth Yule model) between species from neutral coalescent processes within species by analysis of time intervals between branching events (which, in turn, can be summarized as combination of independent Poisson processes) in time-calibrated single gene trees. Input prerequisites require a well-sampled, well-estimated, ultrametric single neutral locus tree that ideally represents the true species genealogy in absence of population structure and population size fluctuation. The method defines sets of species hypotheses based on single or multiple threshold times that potentially distinguish coalescent events from speciation events, and searches for a single maximum likelihood model of mixed speciation and diversification processes across the search space, i.e., sets of species hypotheses. Species hypotheses thus delineated correspond to the phylogenetic species concept.
(III) We analyzed single gene trees and a *BEAST species tree using the PTP [Poisson Tree Processes] model via the PTP webmask [76] using both heuristic ML and Bayesian implementations of the PTP algorithm. Somewhat similar to the GMYC model, the PTP model aims to discern speciation processes among species from diversification processes within species, but analyses numbers of substitutions between branching events instead of time intervals. Input prerequisites enforce the same assumptions as the GMYC model, but this method does not require an ultrametric input tree to delineate entities corresponding to the phylogenetic species concept. This delineation is instead achieved by heuristically inferring species delimitations and searching for a delimitation pattern that maximizes likelihood of a mixed model describing speciation and diversification processes as two independent Poisson process classes across the search space, i.e., sets of species hypotheses. Removing the outgroup in initial runs did not affect delimitation results; we consequently did not use this option.
(IV) We performed combined species tree estimation and species delimitation analysis as available via Species Tree And Classification Estimation, Yarely [STACEY] in BEAST 2 [90]. Simplified, this method (as its predecessor, DISSECT [91]) is an extension of the multispecies coalescence model used in *BEAST [75], in which a birth-death-collapse model is used to estimate the species tree [90,91]. Further, specialized operators are included that model population sizes along branches, prune and regraft subtrees, modify node heights, and merge tips to minimal clusters. The method aims to maximize tree likelihood over a Bayesian tree space by using a MCMC model in which single tips can be merged to minimal clusters to estimate a species or minimal cluster (SMC) tree, while specific priors ensure compatibility between species and gene trees [90]. User-supplied priors define behaviour of MCMC moves, and provide a probability space for the expected number of species primarily through the Collapse Weight prior in combination with a Collapse Height parameter [90,91]. While not extensively tested, a wide range of values [1e-4-1e-7] for Collapse Height was found to provide similar species delimitation results [91]. Contrastingly, the Collapse Weight prior was found to confound delimitation results if fixed, which can be circumvented if estimated during the MCMC process [90,91]. Here, we assess importance of prior space settings for growth rate and population size scaling parameters in STACEY analysis, and test influence of ploidy settings on species delimitation results.
Initially, we assumed a birth-death speciation tree prior while using a Collapse Height of 0.0003, and estimated Collapse Weight with an initial value of 0.5 using a beta prior (1,1) around [0,1]; following suggestions for prior choice in species tree analysis using *BEAST [92], Jeffreys prior was used for growth rate and population scaling factor; the relative death was estimated using a beta prior (1,1) around [0,1]. We used equal ploidy settings, following results and arguments presented in [93][94][95]. We chose this approach to avoid disproportionate influence from mitochondrial partial sequence data and, consequently, treat each gene tree as likely as any other to diverge from the species tree. We did not modify the Collapse Height as preliminary experimental runs confirmed the patterns described in Jones et al. [91]: we found that values larger than or equal to 1e-3 lead to merging of all included specimens to very few (1-2) groups. Further, the NodeReheight operator was set to 3× its value as suggested by [90]. Genealogical relationships were estimated via STACEY 4× independently (1 × 10 7 generations, sampling every 5,000th generation) after incorporating suggestions obtained from an initial run. STACEY log files were examined as stated above. Support for tree topologies estimated by STACEY were examined by constructing a maximum clade credibility tree running TreeAnnotator v1.8 after discarding the first third of all estimated trees. Species delimitations based on trees estimated by STACEY were assessed using speciesDA ( [96]), using the same burn-in, a collapse height of a tenth of the average branch length (corresponding to a value of 0.0005), and default similarity cut-off. We conducted additional analyses to explore the sensitivity of the method to ploidy settings, and prior space for growth rate and population scaling factor. We used different combinations of ploidy settings for mitochondrial loci (using either the same value as for nuclear loci, or ¼ of the value used for nuclear loci which is commonly used in species tree analysis), Jeffreys priors, and logarithmic normal priors for growth rate and population scaling factor. Logarithmic normal priors were used to estimate both parameters, where one set mimicked empirical posterior distributions of growth rate and population scaling factor (where growth rate M = 2.5, S = 1.2, and population size M = −8.5, S = 2, respectively) while the others covered a prior space around M∈ [3,5,7], S∈[1.5, 2.5, 3.5] (resulting in nine possible prior combinations). Further, we edited ploidy settings in the original setup file and re-ran the original STACEY analysis to check for congruence with the new version as we noticed version-dependent differences between setup file templates. These analyses were run for 10 × 10 9 generations, sampling every 10,000th generation, and analysed as described above. In total we thus tested the same data set with 24 prior setting combinations (cf. Table 3.) Both GMYC and PTP methods were run on single gene trees estimated via BEAST exclusively, as taxa delimitations using both methods on BEAST trees were found to be consistent [97]. Further, both methods were used to estimate taxa delimitations on a *BEAST species tree. The practice of inferring species hypotheses on a species tree estimated through *BEAST represents a violation of many of the assumptions that lay the base for both the GMYC and the PTP model and we strictly advise against this approach. Here, however, we took this measure to allow for a more comprehensive comparison of species delimitation methods.

Performance of molecular species delimitation methods
Results from species-delimitation methods were found to be incongruent ( Fig. 1, Table 2). Benchmarked against morphologically differentiable taxa, the majority of methods did not recover relevant species hypotheses. Only STACEY produced estimates corresponding to morphologically differentiable entities.   Specimens referred to by unique identifiers are grouped by species; results are displayed relative to morphologically distinguishable taxa where identical minuscules indicate group membership as inferred through molecular species delimitation methods, where each column represents an independent analysis. Grey cells indicate conflicts of species hypotheses and the taxonomic benchmark, white cells indicate congruence of morphological delimitation and molecular species delimitation. Results of single-locus species delimitation methods are grouped by locus. Labels: A, results of STACEY analysis; B, results of ABGD analysis of a concatenated sequence dataset; C, results of PTP analysis of a *BEAST species tree; D, results of bPTP analysis of a *BEAST species tree; E, results of single threshold GMYC analysis of a *BEAST species tree; F, results of multiple threshold GMYC analysis of a *BEAST species tree; 1, results from ABGD analysis; 2, results from PTP analysis; 3, results from bPTP analysis; 4, results from single threshold GMYC analysis; 5, results from multiple threshold GMYC analysis. For analytical details, see text. Order of specimens corresponds to order of specimens in Fig. 1 Table 2).

Effect of prior choice on STACEY delimitation results
Topology of final maximum clade credibility trees did not differ between runs while node support and branch lengths differed between single runs. STACEY recovered the same species hypotheses in runs with congruent ploidy settings, independent of growth rate and population scaling factor prior choice (Table 3). Ploidy settings considerably affected species delimitation results: When using the same ploidy settings for all loci, STACEY recovers 14 species hypotheses corresponding to the 14 morphological species included in the majority of analyses (Fig. 1., Tables 2 and 3).
However, using one quarter of the nuclear ploidy value in STACEY analyses results in a higher number of species hypotheses (Table 3). These analyses concordantly suggest 17 groups [12 morphological species] (overpartitioning D. osogovicus and D. zivici sp. nov.) and support D. popovi and three distinct D. zivici sp. nov. groups as different entities using a collapse height corresponding to one tenth of the average branch length (oscillating between 0.0006 and 0.0007 between different analyses) in species delimitation analysis. Females of the new species are most similar to females of D. popovi but exhibit (1) in dorsal view distinct, rounded lateral shoulders of segment X, (2) in dorsal Table 3 Tabular summary of results of molecular species delimitation via STACEY under different parameter settings Specimens referred to by unique identifiers are grouped by species; results are displayed relative to morphologically distinguishable taxa and identical minuscules indicate group membership as inferred through STACEY using different prior settings. Grey cells indicate conflicts of species hypotheses and taxonomic benchmark, white cells indicate congruence of morphological delimitation and molecular species delimitation. Compound labels separated by slashes indicate prior settings for growth rate and population scale factor. Labels: J, Jeffreys prior; E/E, empiric priors corresponding to a logarithmic normal growth rate prior around (2.5, 1.2) and a logarithmic normal population scale factor prior around (−8.5, 2); S, logarithmic normal prior (M = 3, S = 1.5); M, logarithmic normal prior (M = 5, S = 2.5); L, logarithmic normal prior (M = 7, S = 3.5); mth, ploidy settings of mitochondrial loci corresponding to ¼ of the value used for nuclear loci; mtd, equal ploidy settings for all loci; asterisk, original analysis. For analytical details, see text. Order of specimens corresponds to order of specimens in Fig. 1 view a ragged outline of the lateral lobes of segment X. Drusus popovi females have an evenly rounded lateral outline of segment X in dorsal view and evenly rounded lateral lobes of segment X.

Species description
Larvae of the new species are most similar to D. serbicus Marinković-Gospodnetić as larvae of both species have an intermittent lateral line ( [98,99]), but exhibit a pronotum with a distinct, rounded pronotal ridge (type B sensu [98]). Larvae of D. serbicus have an annular pronotal ridge (type E sensu [98]).
Male genitalia (Fig. 2a-e) Tergite VIII brown, in dorsal view cranially broadly incised; lacking setation; spinate area as two suboval laterocaudal lobes medially connected by a broad band of spines, embracing a medial wide indentation, flanked by membraneous less sclerotized areas. Ninth abdominal segment (IX) in caudal view ventrally as wide as dorsally; in lateral view medially with a distinct, rounded caudad protrusion and a ventral protrusion, embracing the base of the inferior appendages. Superior appendages in lateral view subcircular, dorsally elongate; in caudal and dorsal view compressed, medially concave. Intermediate appendages in lateral view rounded, rough, high (higher than the superior appendages); in dorsal view the tips medially more proximal, extending widely laterally: bar-shaped, laterally rounded, rough; in caudal view rectangular, tips as wide as the base of the intermediate appendages. Inferior appendages (gonopods sensu [62]) in lateral view subovate with a slight triangular dorsomedial protrusion, straight; in caudal, dorsal and ventral view proximal part broad, distal part slender, straight. Parameres simple, with a distinct medial thorn-like spine preceeded by 2 smaller spines.
Female genitalia (Fig. 3a-d) Segment IX setation abundant, concentrated in the dorsal half; lateral lobe of segment IX membraneous, in lateral view oblique triangular, the ventral edge ventrally protruding and approximately twice as long as the dorsal edge, with a distinctly sclerotized setose dorsal part protruding caudally, in dorsal and ventral view slender, projecting laterally, in caudal view dorsal sclerotized setose part somewhat triangular. Segment X in lateral view oblique rectangular, with a distal, triangular caudal protrusion; in dorsal view hexagonal, medially widest, with rounded shoulders, 2 small median caudal lobes, and distally with 2 semicircular, irregularly-edged lateral lobes, each with a lateral setose area; ventrally unsclerotized, open. Supragenital plate in lateral view quadrangular, ventrally longer than dorsally, with a small dorsal protrusion, caudal line slightly indent; in ventral view quadrangular with a medioventral indentation forming 2 rounded lateroventral lobes; in caudal view wide rectangular, dorsally slightly wider than ventrally, with a dorsomedial protrusion; in dorsal view with a medial dorsal protrusion. Vulvar scale in lateral view triangular, rather straight, approximately as long as the supragenital plate; in ventral view slender with 3 lobes: 2 lateral lobes, digitiform, roundly oval, straight; 1 median, digitiform, of lesser width than lateral lobes, length approximately 2/3rds of that of lateral lobes.
Legs brown, distally slightly lighter (Fig. 4e-g). Abdomen white (Fig. 4h), dorsal gills from II praesegmental position to V praesegmental position, lateral gills absent, ventral gills from II postsegmental position to VI postsegmental position; lateral line as a short segment in the first quarter of II, continuously from last quarter of II to first half of VIII (Fig. 4j); abdomen I with 1 dorsal and 2 lateral protuberances, posterior sclerites absent on lateral protuberances, setal areas sa1-3 fused dorsally and ventrally, sternum bearing 2 setae with distinct basal plates; abdomen VIII with 2 long and 2 short posterodorsal setae on either side; abdomen IX with 1 posterodorsal seta on either side, dorsal sclerite IX semicircular, pale brown with 7 long and several shorter setae (Fig. 4h-i). Case simple, constructed of mineral particles.

Etymology
Named for Miroslav Živić, biophysicist, for his continuous support of faunistic surveys.

Drusinae diversity in the Balkans
Integrated assessment of morphology and results from multi-locus molecular species delimitation corroborate the new species D. zivici and indicate a sister species relationship with D. popovi. The description of this species increases total Drusinae diversity in the Balkans to 43 species [50][51][52][53][54][55][56][57][58]64], demonstrating the significance of this region comprising a highly endemic fauna and flora for European biodiversity [100]. High diversity of aquatic taxa in the Balkans can be attributed to species communities of isolated habitats [101][102][103][104], including numerous Drusinae species. Intriguingly, Drusinae distribution patterns in the Balkans indicate continuance of segregated populations potentially dissociated by climatic or orogenetic processes that in combination with limited dispersal potential [52,105] resulted in small-scale allopatric speciation and extant disjunct distribution of species. As several other taxonomic groups exhibit combinations of ecological traits similar to those of Drusinae (e.g., some Leuctridae, Consorophylax spp.; [15,[106][107][108]), Drusinae may serve as model taxon to estimate unexplored biodiversity. Indeed, recent taxonomic surveys recovered several new Drusinae species in the Western Balkans, demonstrating the necessity and urgency of organismic studies in Europe [109]. Additionally, the Balkan aquatic biodiversity is imminently threatened by anthropogenic habitat modification and climate change [110,111]. Endemic taxa such as the majority of Drusinae species are particularly vulnerable to both climate change and habitat degradation [112][113][114][115]. Socio-economic change and the resulting acceleration of natural resource exploitation [110,[116][117][118][119][120][121] threaten natural habitats world-wide with deleterious effects on biodiversity.
Conservation measures to counter such development crucially depend on comprehensive faunistic data to identify the most critically imperiled habitats and implement adequate protection measures. While such comprehensive data is currently not available, the highly diversified and ecologically specialized Drusinae could serve as umbrella taxon for highland stream and spring habitats in the Western Balkans.

Model taxon, limitations of the molecular dataset and molecular species delimitation
The Drusinae are morphologically well-defined, yet high diversity and intricate morphological characters make identification of several species challenging. Recognition and delineation of new species thus requires high taxonomic skill, and is best supported by additional, independent information such as molecular sequence data. We therefore assume that the Drusinae are representative of taxa that are likely to be investigated using molecular species delimitation tools. The molecular markers used in this study were successfully used to infer species-level relationships in Drusinae and were also found to corroborate taxonomic status of various new species within Drusinae [55,56,64,68]. We thus assume they are sufficiently informative to support molecular taxon delimitation. Indeed, combined species tree estimation and taxa delimitation inferred on a multi-locus data set via STACEY delineated 14 entities corresponding to the 14 morphologically distinguishable species included.
However, the majority of molecular species delimitation tools exploiting single-locus data only did not produce conclusive taxon hypotheses when benchmarked against morphologically distinct entities. Likewise, automated barcode gap analyses on a concatenated molecular data set and GMYC and PTP analyses of a species tree estimated through *BEAST did not recover species hypotheses corresponding to morphologically identifiable species.
Taxonomic estimates based on different loci differ distinctly, indicating a certain necessity to select loci based on their informativeness [122]. We found delimitation results inferred on single locus data to reflect locus variability/proportion of parsimony-informative sites, as more variable loci led to a higher number of proposed species hypotheses. However, the majority of these delimitations are taxonomically not informative. We consequently suggest that, if only single locus data be used, a locus should be selected (ideally through benchmarking locus informativeness based on a known set of species-level taxa) that is conservative for single species in the focal group (e.g., [67]). Similar to our results, other studies report overpartitioning or overaggregation of morphological species by single-locus molecular taxa delimitation tools. Differential variability of molecular markers and differences in phylogenetic inference method used as well as effects of population sizes and speciation rates, but also low traditional taxonomic resolution were previously found to affect significance of molecular species delimitation [97,[122][123][124][125]. Altogether, single locus molecular species delimitation and identification tools currently seem limited in their general applicability. Results obtained here corroborate overaggregation and overpartitioning of morphological species by several single-locus molecular species delimitation tools, indicating a moderate potential as additional information source in integrative taxonomy [97,[122][123][124][125]. Ultimately, this problem is also related to deviation of single gene trees from the true species tree [126,127]. STACEY likely outperforms single-locus approaches by directly integrating available information in a multispecies coalescent model to estimate a SMC tree. Further, we found the STACEY algorithm to be rather impervious to deviations in prior space controlling model speciation rate and population sizes. However, substantial variation of delimitation results can be induced when modifying ploidy settings. In line with [93][94][95] we argue that using the same ploidy settings for all loci represents a more robust approach to estimate species or minimal cluster trees by equally weighting variability in each locus. Our results strongly support promoting STACEY as the method of choice for integrative taxonomy, at least in caddisflies. Yet, further studies are necessary to gauge how well this method performs in other groups.

Integrating molecular species delimitation and traditional taxonomy
Molecular species delimitation and identification tools can be used to infer operational taxon units for evolutionary or ecological analysis, or to obtain alternative taxa hypotheses or discover distinct evolutionary lineages in morphologically uniform taxa [128][129][130][131][132][133][134]. Using these tools is particularly beneficial for studying evolutionary or ecological patterns in taxonomically understudied areas [128,134]. When integrated into larger databases these data on evolutionary lineages and species can be further used for large-scale evolutionary or phylogenetic studies [66]. Also, taxonomic resolution in ecological assessment of water quality could be greatly increased through molecular species identification [131,135]. Such an approach would impart increased and standardized resolution of existing multi-metric indices and thus harbours the potential to enhance existing water resource management schemes.
However, while integrative taxonomy based on (initial) molecular genetic species identification has been proposed as remedy to the taxonomic impediment [136][137][138], the majority of taxonomic studies do not exploit this opportunity. This is likely due to financial or temporal limitations, and because species delimitation studiesincluding some of our owndo not follow through with formal species description or implementation of other taxonomic consequences.
Nevertheless, the potential to develop and test species hypotheses using molecular data is likely to act as an incentive for accelerated taxonomic and ecological research. Collaborative integrative taxonomic projects comprising molecular taxonomists and classical taxonomists are particularly likely to expedite discovery of global biodiversity. Increased collaborative taxonomic efforts will further provide a wealth of information for other disciplines of biological sciences, like biomonitoring of aquatic ecosystems [132,133,[139][140][141], or conservation ecology [142][143][144][145][146]. Currently, constraints on integrative taxonomy, community meta-barcoding, or eDNA approaches in environmental monitoring assays are imposed by inadequate completeness and precision of molecular databases [19,[147][148][149][150]. We anticipate collaborative integrative taxonomic approaches will accelerate alpha-taxonomy and thus provide reference data in addition to on-going international efforts to develop reference libraries. Further, development of conservation management plans or novel tools for aquatic ecosystem assessments will benefit from data thus compiled.

Conclusions
Only multi-locus molecular species delimitation via STACEY reliably delineated molecular species corresponding to morphologically identifiable taxa, confirming a priori expectations on taxonomic significance of different molecular species delimitation tools. We assume that taxonomically relevant molecular species delimitation tools hold potential to accelerate identification of new species (cf. [137,138]), local and global biodiversity estimation and thus enforcement of conservation policies [151] by providing a meaningful assessment of biodiversity richness. However, the capacities of purely molecular species identification based on single-locus data seem to have been overestimated. This clearly demonstrates that while molecular genetic procedures will likely be of relevance in routine monitoring applications, they are currently not fit to serve as surrogate for 'classical' explorative taxonomy [152][153][154].
Nevertheless, we expect multispecies-coalescent-based molecular species delimitation to mitigate the taxonomic impediment and accelerate taxonomic and ecological studies. Under the prevailing biodiversity crisis we direly need to uncover what we are losing fastestlife's uncharted diversity [150,154]. Thus, we advocate a truly holistic integrative taxonomy in order to comprehensively scrutinize our planets declining biodiversity, and so provide essential information for applied biologists, ecologists, conservationists and policy-makers.

Additional file
Additional file 1: Specimen and data set details. Description of data: Summary of collection details, molecular sequence data used and voucher specimen storage. Abbreviations: Lat, latitude; Lon, longitude. (PDF 415 kb)