Freshwater mussels house a diverse mussel-associated leech assemblage

Freshwater mussels (Unionida) are one of the most imperiled animal groups worldwide, revealing the fastest rates of extinction. Habitat degradation, river pollution and climate change are the primary causes of global decline. However, biological threats for freshwater mussels are still poorly known. Here, we describe a diverse ecological group of leeches (Hirudinea: Glossiphoniidae) inhabiting the mantle cavity of freshwater mussels. So far, examples of mussel-associated leech species are recorded from East Asia, Southeast Asia, India and Nepal, Africa, and North America. This group comprises a dozen glossiphoniid species with a hidden life style inside the mantle cavity of their hosts largely overlooked by researchers. We show that the association with freshwater mussels evolved independently in three leech clades, i.e. Batracobdelloides, Hemiclepsis, and Placobdella, at least since the Miocene. Seven mussel-associated leech species and two additional free-living taxa are described here as new to science.


Results
Integrative taxonomy of the mussel leech assemblage. We found that leeches are commonly encountered within the mantle cavity of freshwater mussels in East Asia, Southeast Asia, and East Africa ( Fig. 1 and Table 1). In total, the mantle cavities of 3,045 freshwater mussel specimens were examined, with 370 freshwater mussels being infested by 1,334 leeches (except larvae) (Supplementary Dataset 1). The phylogenetic analyses, species delimitation modeling, and morphological assessments reveal that these mussel-associated leeches belong to nine species within two genera, Hemiclepsis and Batracobdelloides (Glossiphoniidae) (  Map of global distribution of leeches associated with freshwater mussels, prevalence of leech infestation of freshwater mussels in the Old World, and living examples of mussel-associated leech species from the mantle cavity of freshwater mussels. (A) Map of global distribution of mussel-associated leeches (species richness is given in open circles). North America: Placobdella montifera and P. parasitica [ 21 ]. The map was created using ESRI ArcGIS 10 software (https://www.esri.com/arcgis); the topographic base of the map was created with Natural Earth Free Vector and Raster Map Data (https://www.naturalearthdata.com) and HydroSHEDS (https://www.hydrosheds.org) (Map: Mikhail Yu. Gofarov). (B) Prevalence of leech infestation recovered in samples of freshwater mussels (Unionida: Unionidae, Margaritiferidae, and Iridinidae) from East Asia, Southeast Asia, and East Africa (N = 3,045 mussels, primary data:   (Kaburaki, 1921) Possible obligate inhabitant of the mantle cavity of freshwater mussels [23][24][25] Jullundur [Jalandhar, Punjab, India] 23 India and Nepal 23,25,28,79 Secondary host and shelter: freshwater mussels Lamellidens (Unionidae: Parreysiinae) 23 Ancestral area reconstruction. Our ancestral area reconstruction combined from three different modeling approaches (S-DIVA + DEC + S-DEC) suggests that the crown group of the Glossiphoniidae had a continuous range throughout East Asia and North America (probability = 99.9%) (Supplementary Table 5 and Supplementary  Fig. 5). The ancestral range of the subfamily Haementeriinae most likely crossed North America and South America with a subsequent vicariance event (probability = 86.7%), while the subfamily Glossiphoniinae originated in East Asia with subsequent dispersal events (probability = 98.5%). The genus Batracobdelloides most likely spread from East Asia to Africa with a subsequent vicariance event (probability = 96.6%). The MRCA of Asian Batracobdelloides mussel-associated leeches most likely had a broad ancestral range in East and Southeast Asia with a subsequent vicariance event associated with the separation of the Asian freshwater drainages from each other (probability = 99.7%). The genus Hemiclepsis, in its turn, originated in East Asia following by an intra-area radiation (probability = 75.8%). Finally, the Asian Hemiclepsis mussel-associated leeches have had a broad ancestral range throughout East and Southeast Asia with a subsequent vicariance event (probability = 99.7%). The patterns outlined above were also supported by each modeling approach separately (Supplementary Table 5 and Supplementary Fig. 5).
Ancestral life style reconstruction. The MRCA of Hemiclepsis was most likely a free-living species (probability = 99.6%), while the MRCA of the clade with mussel-associated leech species belonging to this genus was most likely an obligate mussel inhabitant as its recent descendants (probability = 90.8%) ( Fig. 2 and Supplementary Fig. 6). The ancestral reconstruction for Batracobdelloides suggests that the MRCA of this genus was a free-living species (probability = 60.1%) rather than a free-living species with a facultative hidden stage inside the mantle cavity of freshwater mussels (probability = 37.4%). Reconstruction of the MRCA of a clade containing two African members of Batracobdelloides returns similar probability patterns (63.6 and 36.4%, respectively). In contrast, the MRCA of the clade of Asian Batracobdelloides mussel-associated leeches within this genus was most likely a mussel inhabitant (probability = 94.2%).
Life cycles and feeding of the mussel leeches. Seven stages were observed in the life cycle of Hemiclepsis mussel-associated leeches in the field, most of which occur within the mantle cavity of a host mussel ( Fig. 6 and Supplementary Table 6). The mature leech leaves the mantle cavity, fixes the egg cluster to the dorsal margin of the host shell near the umbo, and covers the brood by its flat body. After hatching, the larvae attach to the abdomen of the parent, which enters the mantle cavity of the host mussel, with a subsequent internal development of the juvenile leeches into adults. In contrast, only six stages were observed in the life cycle of Batracobdelloides mussel-associated leeches ( Fig. 6 and Supplementary Table 7), because they do not need to leave the host mussel for brooding. These leeches place the egg cluster into a tube-like, enclosed cavity in the median section of their abdomen.  Table 1. Mussel-associated leeches of the World (obligate and facultative inhabitants of the mantle cavity of freshwater mussels) with supplement of an overview of free-living species in the genera Batracobdelloides Oosthuizen, 1986 and Hemiclepsis Vejdovsky, 1884 (complete checklist of these genera is presented as Supplementary Note 1). *Molecular data for these nominal taxa is still lacking. n/a -not available.
The multiple records of abundant larvae and young leeches inside the mantle cavity indicate that the mussel leech species probably use freshwater mussels as their secondary host (Supplementary Tables 6 and 7), while direct evidence of this suggestion (e.g. molecular data for the digestive system of the immature leeches) is still to be collected. In contrast, the digestive system (crop) of dissected mature specimens of Batracobdelloides and Hemiclepsis mussel-associated taxa is filled with a dark red or brown blood-like substance that does not resemble the translucent mussel body fluids ( Supplementary Fig. 10). The COI sequences of the crop content of mature leeches reveal that they feed on blood of freshwater fishes that can be considered the primary hosts (Table 1 and Supplementary Table 8). Both groups of the mussel leeches disperse at the adult stage when they leave their mussel hosts feed on fish blood.
Host range and abundance of mussel-associated leeches. Two leech species associated with freshwater mussels, i.e. Hemiclepsis kasmiana comb. rev. and Batracobdelloides tricarinatus, could be considered secondary host generalists, each of which inhabits mussel species from different families ( Table 1). In contrast, other mussel leech species appear to be associated with one or a few genera of freshwater mussels belonging to a single subfamily or even tribe ( Table 1). The proportion of mussels infested by at least one leech (Leech Infestation Prevalence index, LIP; Eq. 1) was significantly influenced by the taxonomic affinities of the host mussels at the tribe level in Southeast Asia In Southeast Asia (Myanmar), mussel-associated leeches infest Lamellidentini, Indochinellini, and Pseudodontini mussels, but they were never recorded in members of the Leoparreysiini, Host reconstructions are shown for the clades of interest. Proposed obligate inhabitants of the mantle cavity of freshwater mussels are colored red. Free-living leech taxa with a hidden stage inside the mantle cavity of freshwater mussels are colored green. Free-living leech species are colored blue. Outgroup taxa are not shown.
The highest mean LIP was found in East Asia (LIP ± s.e.m. = 60.0 ± 6.2%; N = 28), while this parameter was lower in Southeast Asia ( Table 9). The mean intensity of leech infestation was also higher in East Asia (ILI ± s.e.m.  Table 9). High or moderate levels of the mean LIP were observed in several river drainages such as the Amur, Gladkaya, and Razdolnaya (Russia), Seomjin and Geum (South Korea), Hyakuken (Japan), Sittaung, Salween, Ayeyarwady, and Hlaingbwe (Myanmar), and the Albert Nile (Uganda) ( Supplementary Fig. 11).
Taxonomy. Here, we introduce nine leech species new to science based on diagnostic morphological and molecular characters. A complete morphological description of every novel species is given in Supplementary Distribution. East, Southeast and South Asia, with a small radiation in Africa (two species) and Europe (one species) ( Table 1).     Etymology. The name of this species is derived from the name of its type locality, the Hlaingbwe River basin.
Morphological diagnosis. Small leech, which could be distinguished from other congeners by a combination of the following characters: dorsum greenish or white, without clear markings, with one central row of triangular, flattened tubercles with rounded apex, and similar lateral tubercles being broadly scattered (Supplementary Table 4).

Molecular diagnosis.
The new species differs from other congeners by five fixed nucleotide substitutions in the COI gene fragment (Supplementary Table 3). Genetically, it is most closely related to Batracobdelloides yaukthwa sp. nov. (mean COI p-distance = 3.8%).
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1).

Molecular diagnosis. The new species differs from other congeners by seven fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1).
Distribution. Bago, Sittaung, Ayeyarwady, and Salween basins, Myanmar. Etymology. The name of this species means "freshwater bivalve" (yaukthwa) in Burmese reflecting its preference to use freshwater mussels as hosts.

Batracobdelloides yaukthwa
Morphological diagnosis. Small leech, which could be distinguished from other congeners by a combination of the following characters: dorsum reddish or brownish, with one central row of spike-like tubercles and separate spike-like tubercles being scattered laterally (Supplementary Table 4). Table 3). Genetically, it is most closely related to Batracobdelloides hlaingbweensis sp. nov. (mean COI p-distance = 3.8%).

Molecular diagnosis. The new species differs from other congeners by nine fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1). Etymology. The name of this species is derived from the East Asian region, in which it is distributed.

Batracobdelloides koreanus
Morphological diagnosis. Small leech, which could be distinguished from other congeners by a combination of the following characters: dorsum light yellow, with seven longitudinal brown stripes and three rows of weakly developed, almost invisible, rounded tubercles; posterior sucker with radial brownish bands (Supplementary Table 4). Table 3) and one fixed nucleotide substitution in the 18S gene fragment [154T]. Genetically, it is most closely related to Batracobdelloides hlaingbweensis sp. nov. (mean COI p-distance = 5.1%).

Molecular diagnosis. The new species differs from other congeners by 10 fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1).
Distribution. East, Southeast and South Asia, with one species in Siberia and Europe (Table 1). Etymology. The name of this species is derived from the Khanka Lake basin, where the type series was collected.

Hemiclepsis khankiana
Morphological diagnosis. Medium-sized leech, which could be distinguished from other congeners by a combination of the following characters: dorsum yellowish or whitish, with six longitudinal broad, smooth brown stripes; posterior sucker without bands, but with dense, diffuse brown dots (Supplementary Table 4). Table 3). Genetically, it is most closely related to Hemiclepsis kasmiana comb. rev. (mean COI p-distance = 3.8%).

Molecular diagnosis. The new species differs from other congeners by five fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1).
Distribution. Khanka Lake Basin, Russia and probably China. Etymology. The name of this species is derived from the country of Myanmar, where it is distributed.

Hemiclepsis myanmariana
Morphological diagnosis. Small leech, which could be distinguished from other congeners by a combination of the following characters: dorsum yellowish, brownish or whitish, sometimes with unclear longitudinal narrow light brown stripes and rows of light brown dashes; posterior sucker brownish, with large white spots marginally (Supplementary Table 4). Externally, it resembles rather Batracobdelloides than Hemiclepsis by a small, weakly separated head and weakly developed brown markings on the dorsum.

Molecular diagnosis.
The new species differs from other congeners by 11 fixed nucleotide substitutions in the COI gene fragment and four fixed nucleotide substitutions in the 18S gene fragment (Supplementary Table 3). Genetically, it is most closely related to Hemiclepsis khankiana sp. nov. (mean COI p-distance = 7.1%).
Life style. A mussel leech species that seems to be a host mussel specialist (Table 1).
Distribution. Ayeyarwady, Sittaung, Bilin, and Salween basins, Myanmar. www.nature.com/scientificreports www.nature.com/scientificreports/ Morphological diagnosis. Medium-sized leech, which could be distinguished from other congeners by a combination of the following characters: dorsum smooth, orange, with seven rows of ovate yellow spots (spots in the lateral rows are located on the edge of the last annulus of each somite); posterior sucker with large yellow spots (Supplementary Table 4). This novel species resembles Hemiclepsis marginata but could be distinguished from it by larger, ovate yellow spots (vs. smaller, rounded spots) and the presence of the central row of spots (vs. the lack of this feature). Table 3). Genetically, it is most closely related to Hemiclepsis khankiana sp. nov. (mean COI p-distance = 9.1%).

Molecular diagnosis. The new species differs from other congeners by 15 fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A free-living leech species (Table 1).
Distribution. Partizanskaya and Ussuri basins, Primorye Region, Russia. Etymology. The name of this species is derived from the Tumnin River, from which the type series was collected.

Hemiclepsis tumniniana
Morphological diagnosis. Small leech, which could be distinguished from other congeners by a combination of the following characters: dorsum brown or orange, with seven rows of white or yellow spots (spots of the central row are joined to a broad white or yellow stripe, while spots in other rows may disappear in large specimens); seven rows of very low tubercles; posterior sucker whitish, with dense, diffuse brown dots (Supplementary Table 4). Table 3). Genetically, it is most closely related to Hemiclepsis kasmiana comb. rev. (mean COI p-distance = 9.5%).

Molecular diagnosis. The new species differs from other congeners by 27 fixed nucleotide substitutions in the COI gene fragment (Supplementary
Life style. A free-living leech species (Table 1).

Discussion
Species-rich mussel leech assemblage and cryptic diversity of leeches. Here, we report the discovery of a species-rich assemblage of leeches being associated with freshwater mussels (Unionida). In the Old World, this assemblage is known from East and Southeast Asia, India, Nepal, and Africa ( Fig. 1 and Table 1). It includes members of two genera, Batracobdelloides (six species) and Hemiclepsis (three species) that were largely overlooked by researchers. We describe seven species new to science from East and Southeast Asia as a supplement to three mussel-associated leech taxa already known from the Old World, i.e. Batracobdelloides tricarinatus (Africa) 26,27 , B. reticulatus (India and Nepal) [23][24][25] , and Hemiclepsis kasmiana comb. rev. (East Asia) 18,20,21 . Although they might have been overlooked, mussel-associated leeches have not been reported from Europe, Middle East, North and Central Asia, Australia, the Indonesian Archipelago, New Guinea, and the Philippines. In the New World, two leech species, i.e. Placobdella montifera and P. parasitica, were repeatedly recorded in freshwater mussels from North America 1,2,10-12,14 , while such findings from South America are still lacking. In total, 12 glossiphoniid leech species are associated with freshwater mussels globally. The species richness of free-living leeches in Asia also seems to be largely underestimated, because an integrative re-analysis of widespread Palearctic species such as Hemiclepsis marginata reveals the presence of cryptic taxa with restricted ranges, e.g. H. schrencki sp. nov. (a vicariate species replacing H. marginata in the Amur River and smaller freshwater basins of the Japan Sea drainage) and H. tumniniana sp. nov. (a possible endemic lineage to the Tumnin River and, probably, several nearest river systems). Our findings are in agreement with previous researches indicating that nominal species of leeches with broad distribution throughout Eurasia actually represent species complexes containing one or several endemic species-level lineages in various regions of the Asian Subcontinent 31 .
Origin of the mussel leech assemblage. The leech association with freshwater mussels evolved independently in three genera: Batracobdelloides, Hemiclepsis, and Placobdella. The Hemiclepsis mussel-associated leeches seem to be the most ancient clade originated near the Oligocene -Miocene boundary (Fig. 2). Four Batracobdelloides species from Myanmar belong to a clade that is fully supported by all the phylogenetic analyses. This clade originated in the Late Miocene. The phylogenetic position of Batracobdelloides koreanus sp. nov. is still uncertain. This species was recovered as sister to the African B. tricarinatus + B. amnicolus clade by the maximum likelihood and MrBayes phylogenies ( Supplementary Figs 2-3). In contrast, the BEAST phylogeny indicates that it is sister to a clade containing four Batracobdelloides mussel leech species from Myanmar (Fig. 2). Here, we chose the latter hypothesis as corresponding to the biogeographic patterns, and in this case, the Asian Batracobdelloides www.nature.com/scientificreports www.nature.com/scientificreports/ mussel leeches form a monophyletic clade of the mid-Miocene origin. As for the Placobdella taxa being associated with freshwater mussels, they represent two independent lineages within the genus ( Supplementary Figs 1 and 4). In general, two monophyletic clades of leeches evolved in close association with freshwater mussels as early as the Miocene, and several additional lineages have originated independently as facultative inhabitants (parasites or commensals) of the mussel mantle cavity.
Slow evolutionary rates and historical biogeography of the mussel leeches. Our fossil-calibrated phylogenetic model suggests that leeches are characterized by slow rates of molecular evolution, with the mean COI substitution rate of 0.63%/site/Myr (95% HPD 0.51-0.75%/site/Myr). Slow substitution rates are a common feature for several "living fossil" taxa, e.g. freshwater mussels 32,33 , coelacanths 34 , anthozoans 35,36 , sturgeons 37 , and puddle fishes 37 . To the best of our knowledge, we report the first reliable substitution rates for the Hirudinea based on a fossil-calibrated phylogeny, and these values can be applied as external rates to calculate time-calibrated phylogenetic models in the future. Our biogeographic reconstructions strongly support the New World origin for the Haementeriinae and the East Asian origin for the Glossiphoniinae. Two mussel-associated leech clades belonging to the genera Batracobdelloides and Hemiclepsis appear to have originated in East and Southeast Asia with subsequent vicariance and intra-area radiation events. These regions are among the most species-rich hotspots of freshwater bivalve diversity at the global scale [38][39][40][41][42] , and the ancient and diverse freshwater mussel faunas most likely supported the origin and radiation of mussel-associated leeches in Asia. The BEAST phylogeny reveals that the African Batracobdelloides taxa were separated from Asian members of this genus in the Early Miocene that corresponds to similar Miocene vicariance events in other freshwater animals such as the radicine pond snails 43 . These events coincide with the period of direct connection between African and Eurasian plates via Middle East since the closure of the Tethys in the Early Miocene (approximately 20 Myr) 44 .
Association of leeches with freshwater mussels. Previous researchers chiefly advocated in favor of the "clandestine shelter" hypothesis that assumes occasional commensal relationship of leeches with freshwater mussels 1,2,[10][11][12]14 . This assumption was based on field observations in North America, revealing a low infestation rate of freshwater mussels by leeches, i.e. Placobdella montifera and P. parasitica 1,10-12 . In southwestern Louisiana, 21 freshwater mussels were infested with Placobdella montifera among 2,300 mussel specimens examined (LIP = 0.91%), and only 28 adult leeches were found in this sample (ILI = 0.01 l.p.m.) 10 . Our novel data from Asia and Africa reveals that at least two clades of mussel-associated leech species could be considered obligate inhabitants of the mantle cavity of freshwater mussels serving as shelter. Furthermore, we propose that larvae and juvenile mussel-associated leeches could feed on mucus and body fluids of freshwater mussels representing secondary hosts for these leech species 18,20,21 . Molecular sequences of the digestive system content of the adult mussel-associated leeches indicate that they leave their mussel hosts periodically to obtain blood of freshwater fishes serving as the primary hosts. Probably, adult leeches need to use one or several higher-calorie fish blood meals instead of nutritionally sparse mussel haemolymph to complete the life cycle, i.e. to ensure the successful development of eggs. This hypothesis agrees with the idea that evolution of blood-sucking leeches has been moving toward parasitism of animals with a nutrient-rich blood 45 . Such a two-host feeding behavior (both obligate and facultative), when fish blood meals are needed at the final stage of the life cycle just before leech reproduction, appears to be a successful adaptation to freshwater environment, in which availability of vertebrate blood is limited, and many leech species are forced to use nutrient-poor body fluids of invertebrates as the primary feeding source 45 .
It has been suggested that the primary selective pressure driving the evolution of parental care in leeches may have been predation on leech eggs and juvenile stages 46 . From this point of view, a hidden life style of mussel-associated leeches inside the mantle cavity of freshwater mussels could be considered a progressive evolutionary trait in brooding behavior helping to protect juvenile stages from predators. There are a few records of other leech species inside the mantle cavity of other freshwater bivalve groups, e.g. Sphaeriidae 14,47 and Dreissenidae 2,48 . However, none of these bivalve inhabitants were reported in association with the Unionida. Although a snail-associated leech assemblage is also poorly known, it seems to be a species-rich entity, with at least eleven leech taxa using freshwater gastropods as hosts only in North America 14 . Stibarobdella moorei, a unique example of a cephalopod-associated marine leech, uses Octopus bimaculatus (Octopodidae) as the primary host 49 . Association of Alboglossiphonia leeches with freshwater bryozoans recorded in Siberia is another unusual example, illustrating a rather "clandestine shelter" commensalism than a host-parasite relationship 50 . In soft-bottom environments, the marine leech Notostomum cyclostomum (Piscicolidae) uses crab exoskeletons as the hard substrate for its cocoon, but this leech species does not feed on crustaceans but on fish blood 51 . These examples highlight that leeches use various invertebrate taxa as hosts, shelters or brooding substrates and that such hidden associations of Hirudinea with other animal groups may be much more common than it was assumed previously.

Methods
Data sampling. Mussel-associated leeches (N = 1,334 specimens) were collected from the mantle cavity of 3,045 freshwater mussels (Unionida: Unionidae, Iridinidae, and Margaritiferidae) using forceps during our broad-scale survey of freshwater mussels in East Asia (Russian Far East, South Korea, and Japan), Southeast Asia (Myanmar) and East Africa (Uganda). For most samples, the number of adult and juvenile leeches (without larvae) in every mussel specimen was recorded (Supplementary Tables 6-7 and 9 and Supplementary Dataset 1) to estimate the prevalence and intensity of leech infestation in each freshwater mussel sample 21,49 . In order to reveal the life cycle of mussel-associated leeches, from every examined mussel specimen we recorded: (1) the presence and position of leeches brooding on the host shell; (2) the presence of mature leeches carrying eggs www.nature.com/scientificreports www.nature.com/scientificreports/ (Batracobdelloides) and those with larvae attached to their abdomen; (3) the presence of leech larvae attached to the host mussel soft body; (4) the presence of juvenile leeches; and (5) the presence of adult leeches (Supplementary Tables 6 and 7). In the same sampling sites, free-living leeches were collected from stones, wood fragments, tree branches, and macrophytes using forceps. The leeches were fixed in 96% ethanol immediately after collecting. The fixation of leeches was not preceded by anesthesia in 10% ethanol solution to avoid DNA degradation in the hot tropical and subtropical environments where these samples were primarily collected. A few additional mussel leech samples (N = 26 specimens) were obtained from the mantle cavity of freshwater mussel specimens in the collection of the Federal Scientific Center of the East Asia Terrestrial Biodiversity, Far Eastern Branch of the Russian Academy of Sciences, Vladivostok, Russia. Images of live leech specimens and their broods were recorded with a digital camera (Canon EOS 7D, Canon Inc., Japan).

Molecular analyses.
New molecular sequences (COI and partly 18S rRNA) were generated from tissue samples of 118 leech specimens. Additionally, we obtained COI sequences of the crop content that was extracted from 14 mature leech specimens belonging to 9 mussel-associated leech species and 2 free-living taxa (Supplementary Table 8) to estimate taxonomic affinities of the primary host species (freshwater fishes).
Total genomic DNA was extracted from 96% ethanol-preserved samples using the NucleoSpin ® Tissue Kit  Table 2). Additional sequences were obtained from GenBank. Five haplotypes of the Oligochaeta taxa were used as outgroup. The sequence alignment was performed for each gene separately using the MUSCLE algorithm implemented in MEGA7 55 . The 18S rRNA alignment was checked through the Gblocks 0.91b online server 56 58 . This test revealed a little saturation effect in the two partitions even under the assumption of an asymmetrical tree. A partition-homogeneity test with heuristic search through PAUP* v4.0a165 59 shared the significant conflict of phylogenetic signals among the partitions in the dataset (P = 0.01). However, we considered that this conflict does not affect the phylogeny because it seems to reflect a copious homoplasy rather than independent histories of the genes 38,60 . The single gene alignments were joined to a two-locus alignment using FaBox v1.5 (http://users-birc.au.dk/palle/php/fabox) 61 . Maximum likelihood phylogenetic analyses were carried out with an online version of IQ-TREE v1.6.11 62 using an ultrafast bootstrap algorithm 63 and an automatic identification of the most appropriate substitution models 64 . We used IQ-TREE because this software package was found to return phylogenetic reconstructions with best-observed likelihoods compared with other available likelihood-based algorithms 65 . Bayesian phylogenetic reconstructions were performed using MrBayes v3.2.6 66 . Two independent runs, each with one cold and three heated (temperature = 0.1) MCMC chains, were conducted for 25 million generations (sampling every 1,000th generation). Convergence of the MCMC chains to the stationary distribution was assessed visually based on the plotted posterior estimates with Tracer v1.7 67 , and 15% of the sampled trees were discarded as an appropriate burn-in. Bayesian calculations were performed at the San Diego Supercomputer Center through the CIPRES Science Gateway 68 . The best-fit evolutionary models applied to each partition (3 codons of COI+ 18S rRNA) in the MrBayes and IQ-TREE runs are presented in Supplementary Table 11.

Species delimitation and diagnostics of the new taxa.
To diagnose the new species, we used two-step procedure based on the phylogenetic and morphological analyses 39 . First, we applied an automatic species delimitation approach to delimit the Molecular Operational Taxonomic Units (MOTUs) that may correspond to biological species. The Glossiphoniidae COI sequence dataset was compiled using our own data (111 sequences) and available materials obtained from GenBank (421 sequences). These sequences were collapsed to 316 unique haplotypes. The COI haplotype of an unidentified Piscicolidae taxon was used as outgroup. The species delimitation was performed using the Poisson Tree Process (PTP) modeling through the PTP web-service (http:// mptp.h-its.org) 69 . This approach seems to be more appropriate for slowly evolving animal groups such as freshwater mussels 39 and leeches. As an input tree, we used a maximum likelihood consensus phylogeny inferred from an online version of IQ-TREE v1.6.11 62 with an ultrafast bootstrap algorithm 63 . Substitution models applied to each codon position are listed in Supplementary Table 11. The resulting PTP species delimitation model was largely congruent with the modern taxonomy of the Glossiphoniidae, supporting the majority of currently accepted species ( Supplementary Fig. 4). Based on this evidence, we concluded that it is an appropriate model to delimit the species-level units in our dataset. An uncorrected COI mean p-distance to the nearest neighbor of each lineage was calculated in MEGA7 55 . Second, each MOTU within the clades of interest was studied using morphological and biogeographic criteria and was compared with the original descriptions of nominal taxa to link each clade to a biological species. Nine species lacking available names are described under this study as new to science. The molecular diagnosis of every new taxon was designed using fixed nucleotide substitutions, which were estimated for each gene separately using a Toggle Conserved Sites tool of MEGA7 55 at 50% level. For the diagnoses, an alignment of congeneric haplotype sequences was performed using the Muscle algorithm implemented in MEGA7 55 . All deleterious mutations were retained for the analyses. www.nature.com/scientificreports www.nature.com/scientificreports/ Morphological study. The external morphological characters (number and position of eyes, annulation, color, papillation, position of genital pores, and body size) 70 were examined on specimens of the new species and related taxa. Body measurements for the new leech species were performed using a stereomicroscope Leica M165C (Leica Microsystems GmbH, Germany) equipped with an ocular-micrometer as follows: body length (BL), body width (BW), width of anterior sucker (AW), and width of posterior sucker (PW). To study the reproductive and digestive systems, leeches were dissected using a standard approach 70 . The images of specimens and their morphological and anatomical details were taken with stereomicroscopes Leica M165C (Leica Microsystems GmbH, Germany) and Zeiss Axio Zoom.V16 (Carl Zeiss AG, Germany).
Divergence dating and substitution rate estimation. Node ages were estimated with BEAST v1.10.4 71 using the same two-locus dataset as for the IQ-TREE and MrBayes phylogenetic analyses (see above). Substitution models assigned to each partition are listed in Supplementary Table 11. A lognormal relaxed clock and Yule speciation process with continuous quantile parametrization were applied as the priors of the fossil-calibrated Bayesian model. To dating the phylogeny, we used one new crown fossil calibration as follows: †Hirudinea indet. Hard minimum age: 210 Ma (Late Triassic). Diagnosis and phylogenetic placement: This fossil from a fluvio-lacustrine deposit was identified as a leech cocoon 72 . We assume that this freshwater fossil could represent a crown lineage of the clade Glossiphoniiformes + (Hirudiniformes + Erpobdelliformes), because the suborder Oceanobdelliformes contains the primarily marine and brackish water families 73 , earlier members of which were rather marine worms (our unpublished data based on an ancestral area reconstruction analysis). Absolute age estimate: Late Triassic, a ∼80-m-thick succession of coal-bearing fluvio-lacustrine deposits, Section Peak Formation (Victoria Group, Beacon Supergroup), Timber Peak in the Eisenhower Range, north Victoria Land, East Antarctica, ~210 Ma, based on stratigraphy and palynological analyses 72  Three runs, each with 100,000,000 generations, were performed at the San Diego Supercomputer Center (SDSC, University of California, San Diego, USA) through the CIPRES Science Gateway 68 . The resulting log files were checked for convergence of the MCMC chains with Tracer v1.7 67 . All the ESS values were recorded >300. The sets of fossil-calibrated trees inferred from the three runs were joined with LogCombiner v1.10.4 71 applying a 10% burn-in and an additional re-sampling at every 10,000 generation. The resulting set included 27,000 binary fossil-calibrated trees, based on which a consensus fossil-calibrated phylogenetic tree was obtained with TreeAnnotator v1.10.4 71 . The COI and 18S rRNA substitution rates (mean values and 95% HPD) were obtained from the combined log file using Tracer v1.7 67 .
Ancestral area and life style reconstructions. For the ancestral trait and area reconstructions with RASP v3.2 75 , we used the set of 27,000 fossil-calibrated binary trees that were combined from three runs of BEAST v1.10.4 (see above). As a condensed tree, we used the user-specified, fossil-calibrated consensus tree, which was calculated based on this set of trees using TreeAnnotator v1.10.4 (see above). Non-Glossiphoniidae sequences were removed from the tree set with the appropriate option of the software. Ancestral area patterns were reconstructed using three probabilistic algorithms: Statistical Dispersal-Vicariance Analysis (S-DIVA), Dispersal-Extinction Cladogenesis (Lagrange configurator, DEC), and Statistical Dispersal-Extinction Cladogenesis (S-DEC) 75 . Seven distribution areas of the in-group taxa were assigned as follows: (A) Africa; (B) East Asia (in a broad sense with Eastern Siberia); (C) Southeast Asia; (D) North America; (E) Europe; and (F) South America (Supplementary Table 2). Several unlikely range constraints (i.e. AD, AF, BF, CD, CE, CF, and EF) were removed from the prior settings of the analyses. The S-DIVA analyses were calculated with the following parameters: max areas = 2; allow reconstruction with max reconstructions = 100; max reconstructions for final tree = 1,000; and allowing extinctions. The DEC and S-DEC analyses were performed with default settings and max areas = 2. In addition to the reconstructions obtained from each analysis separately, we used summary results of all three kinds of analyses, which were combined with RASP v3.2 75 . To reconstruct ancestral life style patterns, we used a Bayesian MCMC analysis 75 . Three life style types were coded as follows: (A) leeches with a hidden life style within the mantle cavity of freshwater mussels; (B) free-living leeches; and (AB) free-living leeches with a hidden stage inside a mussel. The analysis was computed with 500,000 generations (sampling every 100th generation) and 10 MCMC chains (temp = 0.1). Null distribution was not allowed. To exclude the pre-convergence part of the simulation, a 10% burn-in was applied.
Assessment of leech infestation prevalence and intensity in freshwater mussels. Based on our field data (Supplementary Table 9 and Supplementary Dataset 1), the Leech Infestation Prevalence (LIP, %) index was computed using the following equation: where N IM represents a total number of freshwater mussels infested by at least one leech in a given sample, and N M represents a total number of freshwater mussels in this sample 21,49 .
The intensity of leech infestation (ILI, leeches per mussel [l.p.m.]) was calculated using available field data (Supplementary Table 9 and Supplementary Dataset 1) as follows: L M where ∑N L represents a total number of leeches collected from the mantle cavity of all freshwater mussels in a given sample, and N M represents a total number of freshwater mussels in this sample 21,49 .