Extensive cryptic diversity in the cosmopolitan sludge worm Limnodrilus hoffmeisteri (Clitellata, Naididae)

Limnodrilus hoffmeisteri Claparède, 1862 is a common freshwater worm, often regarded as an indicator of organic pollution. The taxonomic status of this species is controversial due to great variation in morphological features. Numerous morphological forms of L. hoffmeisteri are recorded in the literature, especially from Europe and North America. Today, DNA-based species delimitation assumes that species boundaries can be more objectively and effectively estimated using genetic data rather than with morphological data alone. To investigate if L. hoffmeisteri is a single species, 295 worms identified as either L. hoffmeisteri or other similar (congeneric) morphospecies, using currently accepted morphological criteria, were collected from 82 locations in the northern hemisphere. The number of primary species hypotheses (PSHs) was first explored with cytochrome oxidase subunit I (COI), the proposed DNA barcode for animal species, and with data for all specimens. Both automatic barcoding gap discovery (ABGD) and the Bayesian general mixed Yule coalescent (bGMYC) model revealed the existence of ≥25 distinct PSHs (COI lineages) in our dataset. Then, smaller samples of individuals representing major COI lineages were used for exploration of a nuclear locus, the internal transcribed spacer (ITS) region. In the ITS gene tree (81 sequences), generated by BEAST, 16 well-supported terminal groups were found, but not all of these groups were congruent with the PSHs found in the COI tree. As results across these different analyses were inconsistent, we resorted to analyzing reciprocal monophyly between gene trees and used a minimum consensus of all evidence, suggesting that there are 13 separately evolving lineages (=13 species) within our sample. The smallest uncorrected COI p-distance between these species is 12.1%, and the largest intraspecific p-distance is 16.4%, illustrating the problem of species delimitation with a DNA-barcoding gap as the sole criterion. Ten of these species are morphologically identified as “L. hoffmeisteri,” the remaining three can be attributed to morphologically distinct congeneric species. An individual from the type locality in Switzerland was designated as a neotype of L. hoffmeisteri sensu stricto. This worm belongs to one of the ten species, and this lineage is widely distributed in Europe, Asia, and North America. The remaining nine species show a mixed distribution pattern; some appear to be endemic to a restricted area, others are Holarctic. Our results provide clues to the future revalidation of some of the nominal species today placed in synonymy with L. hoffmeisteri. A BEAST analysis, based on previously published and newly generated 16S data, suggested that this complex contains also other species than those studied by us. By integrating additional genetic data, it will be possible to identify these and additional specimens in future studies of Limnodrilus, and the neotype provides a baseline for further revisions of the taxonomy of the L. hoffmeisteri complex.


Introduction
The common oligochaetous clitellate Limnodrilus hoffmeisteri Claparède, 1862, in the family Naididae, subfamily Tubificinae (sensu Erséus et al. 2008), has long been arbitrarily treated as a single and widely distributed species in sediments of various freshwater habitats (Kennedy 1965), where it plays an important role in benthic food webs (Maciorowski et al. 1977;Woodward et al. 2008). By mixing deposited matter through burrowing, feeding, and respiration, it has a profound influence on organic components/contaminants in sediments and is even used in the decontamination of sludge containing hazardous organic compounds and heavy metals (Fischer and Beeton 1975;Martinez and Levinton 1996;Matisoff et al. 1999;Vorob'ev et al. 2010;Zhang et al. 2012). L. hoffmeisteri is also regarded as a biological indicator of eutrophication due to its high abundance in organically enriched sediments and tolerance to hypoxia (Uzunov et al. 1988;Burton 1992). Furthermore, it is one of many clitellates involved as definitive hosts in the life cycles of myxozoan parasites that are severely affecting fish populations (Xiao and Desser 1998;Kent et al. 2001;Atkinson et al. 2007;Marton and Eszterbauer 2012). Not surprisingly, the ecologically important taxon L. hoffmeisteri is frequently used as a model organism in both basic and applied science, despite the fact that morphological variability has prompted debate and controversy over its taxonomic status for more than a century.
L. hoffmeisteri and its congener Limnodrilus udekemianus Claparède, 1862 were originally described from a stream near Geneva in Switzerland (Claparède 1862). After Claparède's work, numerous new species were attributed to the genus (Ratzel 1868;Eisen 1879;Vejdovský 1884;Michaelsen 1900;Southern 1909;Nomura 1913;Chen 1940). All Limnodrilus Claparède, 1862 lack (dorsal) hair chaetae, which are common in most other genera of the Tubificinae, and a sexually mature Limnodrilus individual bears a pair of conspicuous cuticular penis sheaths. In some cases, the penes are several times longer than the diameter of the worm itself, but they are normally withdrawn into deep invaginations of the body wall. The sheaths are tube-like and variably slender, with the distal end of the tube modified, showing great morphological variation both within and among taxa (Brinkhurst and Jamieson 1971;Kathman and Brinkhurst 1998).
Escalating during the twentieth century, the taxonomic literature of Limnodrilus was loaded with discussions on the morphology of the hard structures, i.e., the chaetae and the penis sheaths. In particular, attention was given to the shape of the distal part and length/width ratios of the penis sheaths. Much effort was made to discriminate the interspecific differences from intra-specific variation of these characters within the genus (Kennedy 1969; Barbour et al. 1980;Ladle and Bird 1980;Dzwillo 1984;Ohtaka 1985;Steinlechner 1988;Ohtaka et al. 1990). Chen (1940) pointed out that some alleged differences, e.g., in penis length/width ratios, are due to stage of development, measurement errors, and optical artifacts. Due to inconsistencies in published descriptions and inconclusive results in the above comparative and morphometric studies, the taxonomy of Limnodrilus, with time, has progressed from splitting to lumping. In a global revision of aquatic oligochaetes (Brinkhurst and Jamieson 1971), more than 40 synonyms, plus three intra-specific nominal forms/varieties, were all regarded as a single species, and from an ecological point of view, L. hoffmeisteri became Bthe commonest, most widely distributed tubificid^in the world (ibid., p. 467). Still, however, the recognition of two main morphotypes within L. hoffmeisteri has persisted in the recent (post-1971) literature; in one type of worms, the penis sheath ends in a circular plate, perpendicular to the tube, and in the other type, the sheath bears an asymmetrical, hood-like lobe at its distal end (Brinkhurst and Jamieson 1971).
In recent years, the notion that the taxon L. hoffmeisteri represents a species complex, rather than one species, has gained credibility, along with the development of molecular systematics. Preliminary molecular studies, based on analyses of small mitochondrial 16S rDNA data sets, suggested that the morpho-taxon L. hoffmeisteri is likely to encompass cryptic lineages (Beauchamp et al. 2001;Erséus and Gustafsson 2009;Marton and Eszterbauer 2012). It is also known that several other, common and widely distributed clitellate morphospecies are complexes of cryptic species (i.e., morphologically more or less indistinguishable, but genetically distinct species) (Erséus and Gustafsson 2009;James et al. 2010;Erséus 2014, 2017).
Mitochondrial markers, especially the cytochrome oxidase subunit I (COI) gene, have been favored for initial phylogenetic analyses in a wide range of animal groups, to estimate lineage divergence within a large sample of individuals (Hebert et al. 2004;Bickford et al. 2007;Paz and Crawford 2012;Geiger et al. 2014). However, species delimitation based on mitochondrial genes only may not reflect species boundaries as accurately as an integrative approach exploring a combination of various data sources (Will et al. 2005;Spooner 2009;Dupuis et al. 2012;Carstens et al. 2013). Essentially, a multilocus strategy provides independent estimates of both mitochondrial and nuclear genealogical histories, and congruence among estimates provides strong evidence of actual species divergence. The aim of this study is to explore species boundaries among worms referred to as L. hoffmeisteri. We investigate a large sample of specimens from the northern hemisphere, using a combination of mitochondrial and nuclear markers. In addition, to establish the identity of L. hoffmeisteri Claparède, 1862, sensu stricto, a neotype is selected and described from the type locality near Geneva.

Material and methods
Specimen collection, preservation, and morphological identification Specimens morphologically attributable to L. hoffmeisteri and a few other species of the genus Limnodrilus were collected from 82 sampling sites in 19 countries (Fig. 1) and initially preserved in 80-95% ethanol (details of specimens, collection sites, and sequence accession numbers are provided in Supplementary Table S1). The anterior part of each sampled worm was stained in alcoholic paracarmine solution and mounted in Canada balsam on a microscope slide, following Erséus (1994), to serve as a physical voucher used for morphological examination. The vouchers are deposited in the Swedish Museum of Natural History (SMNH), Stockholm, and (one specimen) in the University Museum of Bergen (ZMBN), Norway. The posterior parts of the worms were used for DNA extraction. The mounted specimens were examined under an Olympus BX60 compound microscope equipped with a digital camera DXM 1200, using species identification keys by Kathman and Brinkhurst (1998) and van Haaren and Soors (2013).
The sample collected in Seymaz River, the type locality of L. hoffmeisteri, was neither morphologically nor genetically homogeneous, but one sexually mature individual from this sample was designated to serve as neotype for this taxon. It represents the morphotype (at this site) showing the best fit with Claparède's (1862) original description.
DNA extraction, amplification, sequencing, and alignment DNA was extracted from the posterior part of each individual using the DNeasy Blood & Tissue kit or the EZNA Tissue DNA kit (Omega Bio-Tek, Norcross, GA, USA). A part of the mitochondrial gene COI gene was amplified and sequenced for all specimens, and for subsets of samples, partial mitochondrial 16S rDNA and the nuclear ITS region were amplified and sequenced. Amplifications were carried out in a 25-μl volume reaction with 1 μl of each primer, 1 μl of template DNA, and 6 μl of water mixed with 15 μl Red Taq DNA Polymerase Master Mix (VWR, Haasrode, Belgium). The thermo-cycling procedures for the three markers are given in Supplementary Table S2. DNA sequencing was performed either by Macrogen Ltd. (Seoul, Korea), Eurofins MWG Operon (Ebersberg, Germany), or by the Beijing Genomics Institute (Beijing, China). Sequences were assembled from the chromatogram files and were manually checked for inconsistencies and polymorphisms in Geneious Pro 6.1.7.
COI and 16S sequences of BLimnodrilus hoffmeisterip ublished in the NCBI database (GenBank) were downloaded and compared with our new sample, to find out which of them were relevant (i.e., genetically close) to our data. Two COI records from NCBI were ignored: one (Acc. No. AF534865) was likely from a species of the naidid genus Bothrioneurum (Erséus, unpublished data), and the other (EU311396) appeared to be from a vertebrate, as suggested by a BLAST search (http://blast.ncbi.nlm.nih.gov/Blast.cgi).
New sequences are deposited in NCBI; for accession nos., see Supplementary Table S1. The previously published NCBI Fig. 1 Sampling locations represented by triangles sequences, 43 of COI and 51 of 16S, and all identified as L. hoffmeisteri were incorporated into the respective datasets (see Supplementary Table S1 for details). The alignment of the protein-coding COI sequences was produced with no gaps in Geneious V6.1.8. The haplotypes of all COI sequences were generated from a COI alignment with ambiguous codes treated as N using DnaSP V5; invariable sites were included and sites with missing data were not considered (Rozas and Rozas 1995). The alignment of haplotypes is available in TreeBASE (accession 20383). The partitions of the ITS sequences (i.e., ITS1, 5.8S and ITS2) were first annotated by the software ITSx (Bengtsson-Palme et al. 2013) and manually checked. Alignments of each partition of ITS and the partial 16S gene were performed with the L-INS-i algorithm in the MAFFT (Katoh and Standley 2013) plugin for Geneious, and the scoring matrix was chosen as 1PAM/k = 2 with a gap opening penalty of 1.53.

Gene tree inferences
Four ultrametric trees of the three genes (i.e., one for each of 16S and ITS, COI with both an Ball sequences^tree and one based on Bhaplotypes^) were separately generated using a strict clock model with a Yule speciation prior in BEAST 1.8.2 (Drummond et al. 2012). Posterior distributions were estimated twice by sampling every 1000 generations in 10 million Markov chain Monte Carlo (MCMC) steps, under the GTR substitution model suggested by the PartitionFinder V1.1.1 (Lanfear et al. 2012), all other parameters were used as default. Two independent outputs were combined by LogCombiner (from the BEAST package; Drummond et al. 2012) and then inspected for convergence (ESS > 200) in Tracer 1.6, discarding the first 1000 trees of each run as burn-in. In addition, uncorrected pairwise genetic distances (p-distances) for recognized species were calculated in MEGA V6 (Tamura et al. 2013).

Single locus and multi-locus species delimitation
Under the hypothesis that the nominal L. hoffmeisteri is a complex of cryptic species, the number of primary species hypotheses (PSHs) were first explored using a larger initial COI alignment (all specimens) utilizing automatic barcode gap discovery (ABGD, http://wwwabi.snv.jus-sieu. fr/public/abgd/). ABGD partitions sequences into PSHs by inferring a range of prior intra-specific divergence from the data itself, under the assumption of a barcoding gap between the intra-species and the inter-species variation (Puillandre et al. 2012).
The COI (all) sequences and the COI haplotype data sets were, respectively, also analyzed to test species boundaries using the Bayesian implementation of the generalized mixed Yule coalescent (GMYC) approach. With the assumption of a constant speciation rate (no extinction) among speciation events, GMYC classifies species by maximizing the likelihood to optimize shifts in the branching patterns of an ultrametric tree (a tree in which all the lengths of paths from the root to the tips are equal) (Pons et al. 2006). This process relies on the prediction that independent evolution (neutral coalescent process) leads to the appearance of distinct genetic clusters, which are separated by longer internal branches (Fujisawa and Barraclough 2013). The GMYC method may introduce phylogenetic error if analyzed on the basis of a single tree only, but bGMYC, which is based on Bayesian methodology, can account for this uncertainty by sampling multiple ultrametric trees (Reid and Carstens 2012). In the bGMYC analysis, 100 trees were sampled from the posterior distribution of trees obtained in the BEAST analyses of COI and its unique haplotype data set. Under each tree topology, 2 × 10 4 generations were ran with a thinning interval of 100 using the bGMYC package 1.02 in R. The first 5000 generations were discarded as burn-in. The convergence of MCMC chains in the bGMYC analysis was assessed by checking the evolution graph of the posterior probability against the number of generations.
To test whether PSHs derived from the ABGD and bGMYC analyses are likely to represent real species, a multilocus Bayesian method (Bayesian phylogenetics and phylogeography; BPP V3.2 (Yang and Rannala 2010), which accommodates the species (population) phylogeny and coalescent processes in extant and ancestral species, was applied to examine the PSHs derived from the ABGD analyses. We tested PSHs under three different prior probability distributions: large θ (2, 10) and moderate τ (2, 100); medium θ (2, 100) and moderate τ (2, 100); and small θ (2, 1000) and moderated τ (2, 100). Each analysis was ran (MCMC 100,000 generations with 8000 treated as burn-in) twice to confirm consistency between runs. However, the BPP analyses failed to provide a fully resolved set of species hypotheses, probably due to the small population size and deep divergences within some of our lineages.

Species delimitation by congruence and reciprocal monophyly
As accounted for in the BResults^section, the number of PSHs s varied widely between the single-locus methods and datasets. We therefore also approached our species delimitation in a stepwise manner based on a number of criteria of congruence. First, we only considered PSHs that had been delimited in at least one of the above analyses as possible species hypotheses. We then lumped PSHs that were sister groups in the respective trees until we had well-supported (>0.99 pp) groups that are reciprocally monophyletic across all COI and ITS (gene) trees.

Morphological analysis
Penis sheaths (length, basal width, shape of distal end) were morphologically examined in adult specimens, and chaetal features (relative thicknesses, and lengths, of distal teeth) were noted for all individuals (Table 1). The variation in the length/ basal-width ratio of the sheaths was quantitatively calculated by cluster analysis in the IBM SPSS Statistics. The software BayesTraits V2.0 was used to search the Pagel's lambda value (range 0 to 1) for the best predictive distribution of the given traits on transformation of the COI phylogeny under a Brownian motion model of trait evolution (Pagel and Meade 2013). The lambda values close to 1 indicate significant phylogenetic signal. The shapes of both penis sheaths and anterior chaetae were used as supplementary information when evaluating gene trees.

General information about the molecular data sets
Two hundred and ninety-five new COI sequences, with NCBI accession numbers XX000000-XX000000 (see Supplementary  Table S1), were obtained in this study. Therefore, together with 43 public sequences, the whole COI data set contains 337 sequences. The COI alignment is 658 bp long, including 293 variable sites, of which 273 are informative. Sixty-six haplotypes were generated from 470 sites in the COI alignment in DnaSP, excluding sites with missing data or ambiguities. New 16S and ITS sequences were obtained from 81 sampled individuals representing the major COI lineages. The 16S sequences, including 51 from NCBI, range between 335 and 514 bp. The 16S alignment is 526 bp long, with 168 of 202 variable sites being informative. The ITS alignment is 1099 bp long, of which 332 are variable, and 260 are informative sites.

Single locus gene trees
BEAST analyses of all COI sequences and COI haplotype produced largely congruent trees, but they are different in some terminal branches ( Supplementary Fig. S1). The ITS BEAST tree is less congruent with the COI trees, with differences in both terminal and deeper branches ( Fig. 2 and Supplementary  Fig. S2).
The 16S gene tree ( Supplementary Fig. S3) showed general agreement with the COI tree and was used primarily for matching previously published 16S sequences with the corresponding COI clusters, and later on, species hypotheses (see Supplementary Table S1). Most of the 51 downloaded NCBI 16S sequences were nested within our own lineages, but five North American ones (AF325976, AF325977, AF325978, and AF325981 from Colorado, and EU160488 from New York) were not. They all appear, however, to be closely related to the species that we refer to as BL. hoffmeisteri IV^below. Finally, a single 16S sequence from New York (EU160493) came out as the sister to 16S sequence of our BL. hoffmeisteri V,^a singleton in our sample (see Supplementary Fig. S3 and  Table S1).
Due to the great similarity between the COI and 16S gene trees, the 16S dataset was not used for primary species delimitation.

Primary species delimitation
The ABGD analyses of the COI dataset resulted in 31 primary species hypotheses (PSHs), when the initial partitions were used ( Fig. 2 and Supplementary Fig. S4). The number of recursive partitions varied wildly under different settings of the prior threshold, but the 31 PSHs from the initial partitioning were the same groups as those found monophyletic in the COI maximum clade credibility trees derived from BEAST, both in the tree based in the whole dataset and in the tree based on unique haplotypes.
The results from the bGMYC analyses based on trees estimated from the two COI datasets were not identical ( Fig. 2 and Supplementary Fig. S1): the analysis of all COI sequences divided the data set into 32 PSHs, whereof 24 were well supported (pp > 0.95) and the remaining eight PSHs got lower support, whereas in the analysis of the haplotype trees, the dataset was divided into 25 PSHs, whereof 15 were well-supported (pp > 0.95). The outcome of the bGMYC analysis based on whole COI sequences not only contained a higher number of well supported PSHs than those based on COI haplotypes but it was also generally congruent with COI PSHs suggested by ABGD, except that the PSHs GA31 was split into PSHs Gb31 and Gb32 (Supplementary Figs. S1-2). The 32 PSHs found in the bGMYC analysis of all sequences are illustrated as one of the columns to the right of the COI haplotype tree (Fig. 2).
The ITS dataset was divided into 16 PSHs by ABGD (Fig. 2), using a relative gap width of 1. If regarded as groups of specimens, these 16 groups were also found as wellsupported clades in the ITS gene tree ( Fig. 2 and Supplementary Fig. S2). Only five groups of these are identical to those found by all analyses of COI, but the other ITS groups were either lumped, or incongruent, with at least some of the COI PSHs, as shown by the vertical bars in Fig. 2. That is, these groups need to be lumped to more inclusive groups to be reciprocally monophyletic for both markers (see also below).
A consensus among all evidence is that a minimum of 13 species exist in our sample, whereof ten are morphologically identified as L. hoffmeisteri, and in Fig. 2, these are denoted as lineages I-X. Additionally, four COI lineages (Hap1, Hap8, Hap37, Hap52 in Fig. 2 and Supplementary Table S1) may represent additional hypothetical species within L. hoffmeisteri, but they were represented by specimens for which ITS sequencing failed, i.e., they are not yet supported by any nuclear data. The three other species were identified as Limnodrilus claparedianus (LC), Limnodrilus maumeensis Brinkhurst and Cook, 1966 (LM), and BL. claparedianus-cervix^(CC), i.e., a form morphologically intermediate between L. claparedianus and Limnodrilus cervix Brinkhurst, 1963(see Hiltunen 1967Ohtaka et al. 2006).
There is strong nodal support (PP ≥ 0.99 in COI, 16S and ITS) for each of the currently recognized ten species of L. hoffmeisteri ( Fig. 2 and Supplementary Figs. S2-3, species I-X), as well as the three congeneric lineages CC, LC, and LM, but also full statistical support (PP ≥ 0.99) for nonmonophyly of the L. hoffmeisteri complex as a whole; at least lineages CC and LM are nested within the latter.

Pairwise genetic distances
The species suggested by all data were inputted for an uncorrected pairwise distance calculation. Then, the overall mean uncorrected p-distances for COI is 15.4% and for 16S about 9.7%. The mean distances between species (COI 16.9%; 16S 11.1%) were distinctly higher than the mean divergences within them (COI 2.8%; 16S 1.6%). The smallest COI interspecific p-distance within our dataset was 12.1% (between species VII and species VIII), whereas the largest mean value (16.4%) of intra-specific variation for COI was observed within species IX ( Table 2).

Morphology of male genitalia and chaetae and attribution to species names
We examined all mature worms within our species of the L. hoffmeisteri complex and paid particular attention to the cuticular penis sheaths. We confirmed the existence of two main morphotypes based on the distal end of the sheaths, described previously by, e.g., Brinkhurst and Jamieson (1971): the so-called typical form, where the distal part of the sheath forms an asymmetrical and strongly tilted hood ( Fig. 3 and Supplementary Figs. S6, species VII-VIII) and the plate-topped or spiralis form, where the distal end forms a more or less orthogonal plate, which may be circular (Fig. 3, species I, Figs. 4b and 5j and Supplementary Fig. S5, species I) or somewhat asymmetrical (Fig. 3, species IX, Figs. 4a and 5i and Supplementary Fig. S6, species IX). In many cases, however, individual penis sheaths were difficult to classify using these two basic shapes, and only a few representative penis sheaths of the ten L. hoffmeisteri species and two of the other species are illustrated here ( Fig. 3 and Supplementary  Figs. S5-6). Among the other (congeneric) species, one mature specimen was identified as L. maumeensis (Fig. 3, species LM), as its penis sheaths have a relatively thick wall with an obvious neck (see original description by Brinkhurst and Cook 1966). Two other species in our material had their own distinct penis types: one easily identified as L. claparedianus (cf. Hiltunen 1967, fig . 22), the other as BL. claparedianus-cervix,^which appears to be intermediate between those two species (Fig. 3 species CC; cf. Hiltunen 1967, fig. 23; Ohtaka et al. 2006;Brinkhurst and Cook 1966, fig . 7C).
The length and width of the penis sheaths were measured in 91 sexually mature specimens, selected from nine of the ten species in the L. hoffmeisteri complex, and the other three recognizable congeneric morpho-species mentioned above. There was only one, sexually immature, specimen of species V, so its penis morphology is unknown. The longest penis sheath (1012 μm) was observed in the BL. claparedianus-cervix^lineage (species CC; specimen CNK47), and the shortest one, belonging to species VIII (specimen CE3139), was only about 270 μm long. Basal width ranged from 21 μm (CE8604, species I) to 67 μm (CE1178, species LM). Penis sheath length had a somewhat bimodal distribution, indicating two morpho-groups in our material: one group consists of Fig. 2 The maximum clade credibility tree for mitochondrial COI haplotypes (left) and groups from the corresponding nuclear ITS tree (right). Vertical columns of bars in the middle denote (starting from left) primary species hypotheses estimated by (i) ABGD on COI sequences, (ii) bGMYC on COI sequences, (iii) bGMYC on COI haplotypes and (iv) ABGD on ITS sequences, and (v) final species delimited by reciprocal monophyly in all COI and ITS trees. Species labeled I-X are identified as members of the L. hoffmeisteri complex, CC as BL. claparedianus-cervix,^LM as L. maumeensis, and LC as L. claparedianus. Nodal support of the tree corresponds to posterior supports from BEAST analysis, with a black dot signifying ≥0.99 values. Stars indicate Bayesian posterior probabilities of ≥0.95 in the bGMYC analyses. Hap numbers at terminals refers to COI haplotype of the specimens sequenced L. claparedianus, BL. claparedianus-cervix,^and species III of the L. hoffmeisteri complex, all with penis sheaths longer than 600 μm. The other group contains L. maumeensis and the remaining eight species of L. hoffmeisteri, which generally have shorter penis sheaths (Fig. 6). The three species/ lineages (species III, CC, and LC) with penes >600 μm are not found close to each other on the trees. In the BayesTraits analyses, the lambda values for penis length and the penis sheath length/width ratio were 0.77 and 0.64, respectively, suggesting that these characters carry a phylogenetic signal. However, because the sheath ratio variation is continuous and overlapping, we failed to unequivocally distinguish between nine of the ten species (I-II, IV-X) within the L. hoffmeisteri complex, as well as species III from BL. claparedianus-cervix- (Figs. 3 and 6). The number of anterior chaetae varied considerably, typically from four to eight per bundle throughout the material studied, and the ranges overlapped between species. The shape of chaetae located in the anterior segments (at least in segments II-VI) did not vary much among conspecific individuals, but differences in relative thicknesses and lengths of distal teeth were recognized between some species (Table 1 and Fig. 2). This pattern was also found in immature individuals, although the chaetae were smaller.

Species diversity and distribution
The geographical distribution of the species recognized in the present study is summarized in Fig. 7. Species VII, IX and X, L. claparedianus, and BL. claparedianus-cervix^are truly Holarctic; they occur in all northern continents. Of these, species IX was the most frequent group in our material, with sampling sites in ten different countries. Species I and III were both collected in Europe and the USA. The remaining species (II, IV-VI, L. maumeensis), plus the four unnamed haplotypes 1, 8, 37, and 52 (not yet assigned to species), are each only known from one country or continent. For the time being, they are putatively endemic to their respective areas.
Four countries (Switzerland, Sweden, China, and the USA) were more intensively sampled than others (Fig. 7). Species diversity within them varies greatly, with different combinations of species. The dominant species of each country is not always cosmopolitan. For instance, species VI is strongly dominant in China, and species II is the second most dominant species in Switzerland; these species are so far only known from these countries.
Neotype designation for L. hoffmeisteri sensu stricto With the current delimitation of several cryptic species previously included in this taxon, the inevitable question is the identity of L. hoffmeisteri sensu stricto. Therefore, we here designate a neotype, using a genetically typified specimen from the type locality. Claparède's (1862) original description focused on other anatomical characters, but his illustration of a penis sheath and one anterior chaeta allowed for a good morphological match with only one of the four BL. hoffmeisteriŝ pecies found at this locality by us (species IX, see below, Description of neotype). In particular, Claparède's (1862) illustration (plate 1, fig. 1) showed a relatively straight penis sheath with a circular and orthogonal distal plate; the single chaeta illustrated (plate 3, fig. 12) had short teeth, with the upper tooth only slightly longer than the lower. The other species (II, VII, and X) now occurring at the locality have Specimen from Sweden (CE3093). e Chaetae, somewhat distorted by pressure in slide-mounting. f Atrium, vas deferens and penis sheath. at atrium, ed ejaculatory duct, pr prostate, ps penis sheath, sa spermathecal ampulla, sd spermathecal duct, sf sperm funnel, sz spermatozeugma, vd vas deferens penis sheaths with more angled or hooded distal ends ( Fig. 3 and Supplementary Figs. S5-6); species II has chaetae with the upper tooth much longer than the lower, and species VII has a shorter upper tooth (Table 1 and Fig. 2). The designation of a neotype with molecular information and physical deposition is an essential prerequisite for a future formal taxonomic revision of the L. hoffmeisteri complex and Limnodrilus as a whole.
Summary information about the remaining (still unidentified/unnamed) species, and four additional, but unclassified COI haplotypes studied herein, is provided as Supplementary Document S1. Several junior synonyms are likely to occur in the literature, but they are not yet recognized.
Type locality According to Claparède (1862), the original material was collected in Switzerland, BDans le lit de la Seime près Villette, canton de Genève^(the type locality). Neotype collected by Yingkui Liu in the same river (today called Seymaz), at 46.199°N, 6.194°E, on 24 Aug 2014.
Other genetic information This species also contains COI Hap43 (CE4924; COI barcode in GenBank KY369547) and Fig. 6 A scatter plot of the length of the penis sheath against the ratio between length and basal width (right corner) of species. Species (SHSs) labeled I-IV and VI-X are identified as members of the L. hoffmeisteri complex, CC as BL. claparedianus-cervix,^LM as L. maumeensis, and LC as L. claparedianus Fig. 5 Micrographs of slide-mounted specimens of the investigated Limnodrilus hoffmeisteri sensu stricto (species IX), specimen IDs and 50 μm scale on lower right. a Neotype. Ventral chaetae in an anterior segment. b A topotypic specimen (CE22807). Chaetae in III. c Specimen from Sweden (CE3093). Chaetae in VIII. d Specimen from China (CE12493). Chaetae in VI. e A topotypic specimen (CE22811). Head of one spermatozeugma. f A topotypic specimen (CE22811). Vas deferens near atrium (left) and near sperm funnel (right). g Specimen from China (CE12493). Distal end of penis sheath. h Specimen from China (CE12498). Distal end of penis sheath. i Neotype. Distal end of penis sheath. j A topotypic specimen (CE22811). Distal end of penis sheath. k Neotype. A whole penis sheath Hap47 (CE1840; barcode KY369478). For further information, see Supplementary Table S1.
Geographical distribution Species is apparently widely distributed in the world; identified by us using combination of molecular and morphological data in material from Europe (Switzerland, Germany, The Netherlands, and Sweden); one GenBank record from Italy (EU117546), Asia (China, Japan, Malaysia, and Thailand), and North America (Canada and USA).
Description of neotype Length about 20 mm (estimated from photo of complete worm taken before the DNA sample was removed), width 0.34 mm in segment VIII, to 0.36 mm in XI. Pharynx in segments II-III, chloragogen cells dense on gut beginning in V. Segments II-VIII with 5-8 (median 7) fully developed chaetae per bundle (Figs. 4c and 5a), and often a partially developed replacement chaeta; post-clitellar segments mostly with 4-5 chaetae per bundle; dorsal and ventral bundles similar. All chaetae bifid, sigmoid, with nodulus at about the distal 1/3; teeth short and moderately divergent, approximately equal (4-5 μm long), or with upper tooth slightly longer; upper tooth usually slightly thinner than lower. Ventral chaetae absent in segment XI; no modified spermathecal chaetae in X. Chaetal length 100-132 μm in segments II-IX; 100-110 μm in middle segments, similar in dorsal and ventral bundles.
Clitellum weak, segments X-XII. Male pores in line with ventral chaetae, about level with dorsal chaetae in segment XI; pore opens into a simple penial bursa about 170 μm deep. Spermathecal pores in line with ventral chaetae, about midway between chaetae in segment X and anterior septum (9/10). Sperm funnel conical, lips about 100 μm long, projecting into segment X from 10/11. Vas deferens difficult to follow, but at least 2600 μm long; diameter to 42 μm near atrium, somewhat thinner (26-36 μm) in middle portion, and narrower (26 μm) near sperm funnel. Vas deferens composed of a single layer of ciliated epithelium to 10 μm thick (thinner near sperm funnel), without obvious muscle layer. Atrium mostly in segment XII within gonadal sacs; narrowly fusiform, narrowing gradually to ejaculatory duct, and more abruptly to vas deferens; length 620 μm, diameter to 100 μm ental to prostate attachment, about 60 μm in ectal part (Fig. 4a). Atrium with a very thin (2 μm) muscle layer, and a thick epithelium; epithelial cells granular, irregular, with indistinct boundaries. Atrium joined near midpoint by a single, stalked prostate gland; prostate a fan-like cluster of irregular lobes; entire mass to 250 μm high; prostatic cells granular and indistinct. Ectal (ejaculatory) duct of atrium about 300 μm long by up to 30 μm wide, with a thin, non-ciliated epithelium; appearing slightly wrinkled and less evenly cylindrical than vas deferens.
Shaft of cuticular penis sheath nearly straight (Figs. 4a, b, and 5k); wall thickness 2-3 μm. Sheath length 530 μm; diameter 48 μm at basal (proximal) end, tapering to 28 μm near middle, and only slightly narrower (26 μm) at distal end (below head); ratio of length to basal width 11. Distal end (Bhead^) of one penis sheath forms a plate 85 μm wide (Figs. 4a and 5i, k), at most slightly wider on one side, one side bent upward; plate at a slight angle from shaft. Head of second penis sheath nearly circular, appearing more strongly angled from shaft (Fig. 4b), possibly from compression. Penis sheath surrounded by a close-fitting layer of epithelium (to 12 μm thick), continuous with lining of penial bursa and continuing as a thin layer within the sheath; entire structure surrounded by a thick layer (to 30 μm) of loosely arranged, transverse-spiral muscle fibers; near distal end, the muscles detach from sheath and join the body wall behind the male pore.
Pharynx in segments II-III, about equally developed dorsally and ventrally; chloragogen cells dense on gut, always beginning in segment V. Long, convoluted, lateral blood vessels visible in segments II-VII and in IX; in VIII, lateral blood vessels shorter and greatly dilated, modified as Bhearts.D orsal vessel ventrolaterally displaced in middle-posterior segments, beginning in segment X. In CE12498, one pair lateral blood vessels in posterior segments; network of capillary blood vessels not visible in epidermal layer.

Species
Europe Asia  North America  I  AT , CH, DE, DK, SE  US  II  CH  III  BE  US  IV  US  V  US  VI  CN  VII  BE, CH, DE, FI, IT, SE  CN  CA, US  VIII  FI, NL, SE, UK, IT  IX  CH, DE, NL, SE, IT  CN, JP, MY, TH  CA, US  X  CH, EE, FI, IT ,NL, SE  CN  US  LC  CH, NL Male and spermathecal pores as above; paired, inconspicuous transverse openings in line with ventral chaetae. Female pores inconspicuous, on 11/12. Clitellum weak, segments X-XII. Both testes and ovaries large, usually extending through X and XI, respectively.
Remarks This morphological description was intended to verify consistency with Claparède's (1862) original L. hoffmeisteri description. It is not an attempt to review or evaluate the vast morphological literature that has followed the original description. No original types are known to exist (Reynolds and Wetzel 2015). When revisiting Claparède's original locality (Seymaz River, Canton Geneva) in August, 2014, four of the species now recognized as members of the L. hoffmeisteri complex (species II, VII, IX, and X), as well as L. claparedianus, were collected. In addition to this, Vivien et al. (2015), in a barcoding study of Swiss Naididae found species I and L. udekemianus in streams of the neighboring Canton of Vaud. Of these, species IX is the one that best conforms to Claparède's description of L. hoffmeisteri, and therefore, we selected a representative of this lineage from River Seymaz as a neotype of this taxon.
All specimens used in the description were sexually mature, with well-developed testes and ovaries, mature eggs in the egg sacs, sperm on the male funnels, and spermatozeugmata in the large spermathecal ampullae. Female funnels were not obvious in any specimen (which seems typical for the genus). Penis sheath form and proportions are perhaps the most consistently described characters in Limnodrilus literature; however, we note that the commonly reported length/width ratio varies greatly, depending on which part of the sheath is measured. Here, we use the basal (proximal) end of the sheath for width measurements, in keeping with recent morphological literature. Proportions of penes in the new material, representing group IX specimens from three localities, are consistent with Claparède's illustration (Claparède 1862, plate I, fig. 1). Claparède gave a Bcopulatory organ^length/width ratio (possibly based on the terminal width near the distal part) of around 5-6, but measuring his drawing gives a penis sheath length/basal-width of about 11, as noted by Piguet (1913). Claparède's (possibly schematic) illustration shows a round (symmetrical) head, orthogonal to the shaft; in contrast, the head is asymmetrical in some of the present specimens, and the longer side may be somewhat reflexed. The apparent angle of the head may be related to either viewing aspect or compression in slide-mounting. However, penis sheaths are not strongly narrowed and curved distally in our material, and the distal ends do not widen abruptly to form a strongly angled hood as in the Btypical^form of L. hoffmeisteri shown in the recent guides (Brinkhurst 1960, fig. 4G-H;Brinkhurst 1965, fig. 4a;Ohtaka 1985, fig. 5B; Ohtaka et al. 1990, fig. 3A-F;Hiltunen 1967, fig. 20;van Haaren and Soors 2013, fig. 230). Instead, they seem more similar to the Bplate-topped^or intermediate forms, e.g., fig.  5A in Ohtaka (1985), fig. 4 in Ohtaka et al.(1990).
Other characters related to morphology of reproductive organs are not emphasized in recent literature and have not been consistently described. The new material has a relatively elongate atrium, probably more similar to Claparède's illustration than to the atria of some other species (e.g., Limnodrilus profundicola [cf. figure 4 in (Fend et al. 2016)]). Penial musculature is often prominent in the genus, with a characteristic spiral arrangement in some worms that have been attributed to L. hoffmeisteri (Eisen 1885;Chen 1940;Moore 1905). This is clearly visible on the neotype and other material examined here, but Claparède's statement Bcomposée de fibres longitudinales el de fibres circulairesd oes not clearly describe the same structure. Although the form of the spermathecal duct has been used as a taxonomic character within the genus in the past, the original description is rather vague. The neotype and associated specimens have a relatively elongate (tubular) spermathecal duct, compared with, e.g., L. profundicola (Cui et al. 2015, fig. 7;Fend et al. 2016).
Although subsequent comparisons rely heavily on chaetal morphology, chaetae were barely mentioned by Claparède (1862): relatively slender, 6-8 anteriorly, 4-6 in mid-body, 2-3 posteriorly; according to the single illustration (plate 3, fig. 16), the upper tooth is slightly longer than the lower. Later works on species that have been attributed to L. hoffmeisteri have focused on length (particularly relative length) of the distal teeth: most commonly, teeth have been considered subequal, but with the lower tooth sometimes thicker (e.g., Southern 1909;Chen 1940;Brinkhurst 1965;Ohtaka 1985). Chaetal morphology seems rather consistent in the present material; although the upper tooth may be slightly longer than the lower, it is never distinctly longer or shorter; teeth are moderately divergent and not strongly curved.
Claparède's description included internal somatic characters that are difficult to observe in our alcohol-preserved material. There is no evidence of a highly vascularized body wall in posterior segments (as in L. udekemianus) in the single specimen with some remaining posterior segments. The hearts in segment VIII only seem typical for the genus (Chen 1940). Nephridia could not be seen clearly. In some guides (e.g., Timm 2009, van Haaren andSoors 2013), the dense, continuous layer of chloragogen tissue on the gut, beginning on segment V seems an important specific character of L. hoffmeisteri when separating this taxon from L. udekemianus, although Claparède (1862, p. 226) stated that the coating of chloragogen cells (cellules pigmentaires) begins on segment V in both species.

Species delimitation
Single locus data are widely and conveniently used to estimate species boundaries in large samples of specimens of morphology-based taxa (Nicolas et al. 2012;Lv et al. 2014;Trebitz et al. 2015). We tried to disentangle the mitochondrial diversity of the highly variable L. hoffmeisteri, by examining specimens from all continents of the northern hemisphere with ABGD and bGMYC. The parallel analysis of a nuclear marker (ITS) suggested that these approaches overestimated the number of species. ABGD may over-split the hypothetical species groups, as the population mutation rate is hard to estimate from genetic data without an a priori known number of species (Puillandre et al. 2012). The two barcoding gaps observed in our COI ABGD analysis ( Supplementary Fig. S4) may indicate that the number of sampled individuals is insufficient.
In our bGMYC analyses, we used COI sequences of all specimens, as well as unique haplotypes, when we generated the input trees. Unique haplotypes are often used to reduce computational time (Reid and Carstens 2012;Talavera et al. 2013;Lemer et al. 2014), although removing identical sequences may overestimate population size when inferring gene trees (Drummond and Bouckaert 2015). On the other hand, if identical sequences are not removed, zero length branches increase the risk of over-splitting (Reid and Carstens 2012).
Although mitochondrial DNA (COI) generally has a higher evolutionary rate and a smaller effective population size than that of nuclear markers (Rautenberg et al. 2012), by itself, it is generally not sufficient for the accurate assessment of species boundaries (Knowles and Carstens 2007;Esselstyn et al. 2012;Dellicour and Flot 2015). Thus, the primary species hypotheses delimited from a large mitochondrial single-locus data set need to be corroborated by other data. In particular, nuclear markers have been suggested for this purpose (e.g., Elias et al. 2007;Knowles and Carstens 2007;Achurra and Erséus 2013), and we used ITS, together with COI, to delimit our species. The congruence among independent estimates of genealogical history provides evidence of ten actual species in the L. hoffmeisteri species complex. We also conclude that the lineages identified as L. claparedianus, L. maumeensis, and BL. claparedianus-cervix^are species.
In the L. hoffmeisteri complex, there are several examples of deep intraspecific coalescence in COI, with up to 16.4% pairwise distances (for species IX; Table 2). This is particularly noteworthy, as the minimum distance (between species VII and VIII) in our studied sample is only 12.1%. Thus, it cannot be automatically assumed that COI lineages, which are that far apart, represent different species, without testing this with nuclear data. This again illustrates how problematic it is to delimit species with a DNA-barcoding gap as the sole criterion, in Clitellata (e.g., Achurra and Erséus 2013;Martinsson et al. 2013Martinsson et al. , 2015Martinsson and Erséus 2017) as well as in other animal groups (e.g., Munoz et al. 2011;Webb et al. 2012;Hogner et al. 2012)

Morphological variation
In our sample of the L. hoffmeisteri complex, we observed great morphological variation in the length/width ratio and the distal shape of the penis sheaths, with a large overlap across species (Fig. 6). The two main morphological types, i.e., the Btypical^and Bplate-topped^ends of the sheath, were frequently encountered, but there is also considerable variation within some species of the complex, including apparently intermediate forms (see, e.g., Supplementary  Fig. S6, species IX). Therefore, the morphology of the penis sheath alone may not be sufficient for circumscribing species and identifying specimens in this group. Additionally, as the lack of consistent penial features in our evolutionary lineages makes it difficult to match them with taxonomic entities described in the literature, resolving nomenclatural issues within this complex will be a difficult task.
On the other hand, our study indicates that chaetal characters may be correlated with genetic variation within the L. hoffmeisteri complex. Some of our recognized species showed apparent diagnostic features in the anterior chaetae (Table 1 and Fig. 2). For example, specimens with the upper chaetal teeth almost equal to the lower teeth are clustered together in our gene trees. We found that individuals with the lower teeth being slightly longer and thicker than the upper ones constitute one large clade (Fig. 2, species VII-VIII and Hap 52), which is consistent with the description of L. parvus Southern, 1909, to date mostly regarded as a synonym of L. hoffmeisteri (see Brinkhurst and Jamieson 1971). However, Southern (1909, fig. 5) claimed that the penis sheath of his taxon Bbecomes very narrow distally, before expanding into a funnel-like mouth^ (Southern 1909, fig. 5), which seems to fit the species VII and VIII well ( Supplementary Fig. S6).

Geographical distribution
Several species of the L. hoffmeisteri complex are widely distributed, which may have been facilitated by, e.g., human activities or passive transport by birds (Milbrink and Timm 2001;Koel et al. 2010). As many as seven of the ten species in the complex, plus one of the unnamed haplotypes (Hap1), were collected in North America; five of these species are also present in Europe, and three also in China. Similarly, there are seven species in Europe, of which five are also North American. Moreover, it is likely that several of our species are among the many published records of BL. hoffmeisteri^(sensu lato) from the southern parts of the world. Among the other species studied by us, L. claparedianus is represented in all northern continents, and this taxon is also known from all parts of the southern hemisphere except Antarctica (Brinkhurst 1966, Brinkhurst 1971, Brinkhurst and Marchese 1988, Pinder and Brinkhurst 2000. It is noteworthy that Sweden, with its northern location (>55°N) and recent glacial history, is the home of at least five of the ten recognized species of the L. hoffmeisteri complex. However, this high diversity may be explained by a sampling bias.

Conclusion
This study has confirmed that the common taxon L. hoffmeisteri, as previously viewed, is not a single cosmopolitan species. Instead, it is a group of at least ten more or less cryptic species. We recognized some morphological features (penial and chaetal shapes) that may discriminate some lineages within the L. hoffmeisteri complex, but they do not, at this stage, appear sufficient for a full taxonomic resolution of all species. The group as such has been considered to be cosmopolitan, and our study supports this inasmuch as several of the cryptic species are indeed widely distributed in the world, at least in the northern hemisphere.
An ultimate goal of future work should be to classify all different species closely associated with this iconic taxon, and the designation of the L. hoffmeisteri neotype is meant to be an important first step towards a complete revision of this large complex on a global scale. An objective, DNAbased species circumscription of Limnodrilus species will benefit the interpretation of their distributional and Table 2 The maximum interspecific and the minimum intraspecific pairwise distances between COI sequences in the study samples. Species labeled I-X are identified as members of the L. hoffmeisteri complex, CC as BL. claparedianus-cervix,^LM as L. maumeensis, and LC as L. claparedianus   I  II  III  IV  V  VI  VII  VIII  IX  X  CC  LC  LM   I  ecological patterns and will enable a more accurate monitoring of environmental disturbance wherever these abundant and ubiquitous organisms are present.