Biodiversity of the Buffalo Leeches Genus Hirudinaria (Arhynchobdellida, Hirudinidae) in Southern Thailand Revealed from DNA Barcoding

© 2022 Academia Sinica, Taiwan Open Access Biodiversity of the Buffalo Leeches Genus Hirudinaria (Arhynchobdellida, Hirudinidae) in Southern Thailand Revealed from DNA Barcoding Putita Jiranuntskul1, Ausadang Boonporn1, Chalita Kongrit1, Somsak Panha2, and Ekgachai Jeratthitikul1,* 1Animal Systematics and Molecular Ecology Laboratory, Department of Biology, Faculty of Science, Mahidol University, Bangkok 10400, Thailand. *Correspondence: E-mail: [email protected] (Jeratthitikul). E-mail: [email protected] (Jiranuntskul); [email protected] (Boonporn); [email protected] (Kongrit) 2Animal Systematics Research Unit, Department of Biology, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand. E-mail: [email protected] (Panha) Received 9 April 2022 / Accepted 8 October 2022 / Published 26 December 2022 Communicated by Benny K.K. Chan Leeches in the genus Hirudinaria Whitman, 1886, also known as buffalo leeches, are blood-sucking ectoparasites of vertebrates. Although they are widely distributed in Asia and had been highly abundant in the past, studies on diversity and taxonomy of this genus are still scarce. There is probably a large amount of cryptic diversity yet to be discovered, particularly from mainland Southeast Asia. In this study, we used morphology and DNA barcoding with a COI gene fragment to explored the diversity of Hirudinaria leeches in the southern region of Thailand, where a unique geographic feature could have led to the diversification of freshwater biota. Molecular phylogenetic analyses and species delimitation approaches (ABGD, bPTP, GMYC, and BOLD) revealed the presence of four putative species of Hirudinaria leeches from southern Thailand, including H. bpling, H. thailandica, and two morphologically cryptic lineages of H. manillensis. Compared to other leech genera, genetic distances of Hirudinaria leeches were relatively low (0.11–0.65% within species; 3.72–14.36% between species) and barcoding gaps were very narrow (1.54–2.88%). The species diversity, distribution pattern, and a phenomenon of low genetic divergence of Hirudinaria leeches in southern Thailand could be explained by an ancient seaway, paleo-drainage, and anthropogenic activities. Key words: Hirudinea, Species delimitation, Cryptic species, Systematics, Molecular phylogeny. BACKGROUND Leeches in the genus Hirudinaria Whitman, 1886, also known as buffalo leeches, are blood-sucking ectoparasites that feed on many vertebrates (Phillips and Siddall 2009). Their main hosts are mammals, particularly cattle, but they are also found feeding on other vertebrates, such as reptiles, amphibians, and fishes (Lai and Chen 2010). Hirudinaria leeches have been recorded to occasionally feed on human blood (Lai and Chen 2010); however, they also benefit humans in medical applications. For example, dry leech powder is used as traditional medicine and is noted for having a wide range of beneficial pharmaceutical properties (Enguang 2008). Living leeches are also used as an alternative therapy to remove excess fluid from patients (Walsmann and Markwardt 1985; Lent 1986). Moreover, they are used as a model organism in laboratories to Citation: Jiranuntskul P, Boonporn A, Kongrit C, Panha S, Jeratthitikul E. 2022. Biodiversity of the buffalo leeches genus Hirudinaria (Arhynchobdellida, Hirudinidae) in southern Thailand revealed from DNA barcoding. Zool Stud 61:84. doi:10.6620/ZS.2022.61-84. Zoological Studies 61:84 (2022) doi:10.6620/ZS.2022.61-84 1 © 2022 Academia Sinica, Taiwan study the nervous system and symbiosis (Graf et al. 2006; Elliott and Kutschera 2011). Hirudinaria species are widely distributed in Asian countries, including Thailand, the Philippines, Vietnam, Malaysia, Indonesia, China, Taiwan, Bangladesh, Sri Lanka, India, and Nepal (Nesemann and Sharma 2001; Zhang et al. 2008; Lai and Chen 2010). Currently, there are four valid species in this genus, namely Hirudinaria manillensis (Lesson, 1842), Hirudinaria javanica (Wahlberg, 1856), Hirudinaria bpling Phillips, 2012, and Hirudinaria thailandica Jeratthitikul & Panha in Jeratthitikul et al., 2020. The latter two species were recently described from Thailand (Phillips 2012; Jeratthitikul et al. 2020). Although they have wide distribution ranges today and had high abundance in the past, research on the diversity and taxonomy of Hirudinaria leeches is still scarce, with only a few published studies available (e.g., Phillips 2012; Chong et al. 2014; Tubtimon et al. 2014; Zulhisyam et al. 2014; Jeratthitikul et al. 2020). In brief, Tubtimon et al. (2014) found evidence of genetic variation in Hirudinaria species in northeastern Thailand. Chong et al. (2014) reported on two types of H. manillensis in Malaysia, which were placed in different clusters of a COI gene tree and showed variation in the position of the anus. Based on these findings, there are probably a number of Hirudinaria species yet to be discovered, particularly from mainland Southeast Asia. Meanwhile, their natural populations in some areas are severely reduced because of the over-harvesting of leeches, habitat destruction, and water pollution (Lai and Chen 2010). A DNA barcode is a short sequence from a standardized gene region that can be used to discriminate species (Allendorf and Luikart 2007). A DNA barcode is usually conserved among conspecific organisms but has enough genetic divergence to separate two different species (Hebert et al. 2003a). For animals, the region most widely used as a DNA barcode is a part of cytochrome oxidase subunit I (COI) in mitochondrial DNA, because its evolution rate can appropriately reveal both recent and deep phylogenetic relationships (Hebert et al. 2003b). Furthermore, universal primers for this region are available, and have been successfully used for a wide variety of animal taxa (Folmer et al. 1994). In animals with high levels of diversity, DNA barcoding based on the COI gene has been successfully used to resolve taxonomic problems, such as revealing intraspecific variation, morphological crypticity, or phenotypic plasticity. DNA barcoding thus can overcome the limitations of relying solely on morphological data (Hebert et al. 2003b 2004; Hebert and Gregory 2005). Other applications of DNA barcoding include rapid species identification, matching males and females, or individuals at different life stages, and assigning unidentified specimens to a species (Allendorf and Luikart 2007; Dela Cruz et al. 2018). Leeches are a group of animals with high morphological variability, and cryptic genetic diversity has been found in many cases (Siddall et al. 2007). DNA barcoding based on the COI gene has been successfully applied along with other datasets to resolve the taxonomy of different groups of leeches with a high rate of monophyletic recovery and classification accuracy, such as in the class Hirudinea (Siddall and Burreson 1998), subclass Euhirudinea (Apakupakul et al. 1999), order Arhynchobdellida (Borda and Siddall 2004), and family Hirudinidae (Phillips and Siddall 2009). However, classifications at the family and generic level need to be revised because some traditional characteristics, such as mode of cocoon deposition, habitat, and feeding habit have led to nonmonophyletic grouping (Apakupakul et al. 1999; Borda and Siddall 2004; Phillips and Siddall 2009). For species identification, the COI gene has been used to identify medicinal leeches (Hirudo verbena and Hirudo medicinalis), which usually have high morphological variation (Siddall et al. 2007). Furthermore, COI gene fragments have been used to support the validity of new leech species, along with morphological evidence (Wang et al. 2022). Hirudinaria leeches are widely distributed throughout Thailand. Aside from Tubtimon et al. (2014), who surveyed buffalo leeches in the northeastern part of the country, there have been no other extensive surveys of these leeches in Thailand. The southern region of Thailand is of particular interest, because its unique geographic features could lead to diversification of freshwater biota (de Bruyn et al. 2005; Grismer et al. 2016; Bohlen et al. 2020). In this study, we explored the diversity of buffalo leeches (Hirudinaria spp.) in southern Thailand based on morphology and DNA barcoding of a COI gene fragment. MATERIALS AND METHODS Specimen collection, preparation, and morphological identification Animal use protocols in this study were approved by the Mahidol University-Institute Animal Care and Use Committee (MU-IACUC) under the approval number MU-IACUC 2018/003. Buffalo leeches were collected from 22 locations in southern Thailand (Fig. 1, Table S1). Additional samples from other regions were also included in this study. For newly obtained specimens, leeches were euthanized by the two-step method as suggested by AVMA Guidelines for the page 2 of 14Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan Euthanasia of Animals (AVMA 2020). Living specimens were placed in a container filled with fresh water. Then, 95% (v/v) ethanol was gradually added to the container, starting from approximately 5% (v/v) concentration until the specimens were fully anesthetized, which could be observed by fully relaxed body muscles. Specimens were moved to 70% (v/v) ethanol to complete the euthanasia and for tissue fixation. Fixed specimens were preserved in 95% (v/v) ethanol. Some of the specimens were photographed for their live coloration before euthanization. Approximately 1 cm2 of the dorsal muscle was cut and kept in 95% (v/v) ethanol for DNA extraction. Other body parts were kept as voucher specimens, and were deposited into the Mahidol University Museum of Natural History, Bangkok, Thailand (MUMNH). Morphological identification was made based on descriptions in previous publications (Sawyer 1986; Sawyer et al. 1998; Nesemann and Sharma 2001; Lai and Chen 2010; Phillips 2012; Tubtimon et al. 2014; Jeratthitikul et al. 2020). DNA extraction, amplification, and sequencing Genomic DNA was extracted from the muscle tissue using a NucleoSpin tissue kit (MACHERYNAGEL, Germany). The DNA barcoding region of the COI gene was amplified by LCO1490Hiru (5'- ATT CTA CTA GTC ATA AAG ATA TTG G-3') and HCO2198Hiru (5'-AAA ATC AAA ATA TAT ACT TCT GGA TG-3') (Jeratthitikul et al. 2020). The PCR mixture contained approximately 10 ng of genomic DNA, 15 µl of EmeraldAmp PCR Master Mix (TAKARA BIO, Japan), 0.9 µl of forward and reverse primers, and 10.2 µl of ddH2O. The T100 thermal cycler (BIO-RAD, United States) was used to carry out the PCR. The thermal cycling was set as follows: initial step at 94°C for 3 min, followed by 35 cycles of denaturation at 94°C for 45 s, annealing at 40°C for 45 s, Fig. 1. Map showing sampling localities of buffalo leeches genus Hirudinaria in (A) Asia and (B) southern Thailand. Dotted lines indicate hypothetical fauna transition zones in southern Thailand: Isthmus of Kra and Surat Thani-Krabi Line. page 3 of 14 Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan extension at 72°C for 1 min, and the final extension at 72°C for 5 min. The PCR products were checked by 1.5% agarose gel electrophoresis under UV light. Later, the products were purified using the PEG precipitation method. The products were sequenced on both strands using the same primer pair on the ABI prism 3730XL automated sequencer (BIONEER, Republic of Korea). The obtained sequences in both directions were assembled using MUSCLE as implemented in MEGA X v.10.2 (Kumar et al. 2018). All sequences obtained in this study were deposited into the Barcode of Life Data System (BOLD) and the GenBank database under accession numbers OM415425–OM415511 (Table S1). Sequence alignments were generated using the MUSCLE option as implemented in MEGA X. The aligned sequences included 87 sequences of leech samples newly obtained in this study, 22 sequences of Hirudinaria species from previous studies, and sequences of four outgroups (Hirudo medicinalis, Hirudo orientalis, Haemopis sanguisuga, and Poecilobdella nanjingensis). The number of unique haplotypes was calculated using “haplotypes” package (Aktas 2020) in R studio (RStudio Team 2020). Phylogenetic analyses The best fit nucleotide substitution models for phylogenetic tree reconstruction were selected by jModelTest v.2.1.10 (Guindon and Gascuel 2003; Darriba et al. 2012). The TPM3uf+I+G and HKY+G models were selected for maximum likelihood (ML) and Bayesian inference (BI) trees, respectively. The ML gene tree was reconstructed in IQ-TREE v.1.6.12 (Nguyen et al. 2015) with 1,000 bootstrapping tests. For BI analysis, a tree was conducted in MrBayes v.3.2.6 (Ronquist et al. 2012) via CIPRES portal (Miller et al. 2015). Ten million MCMC generations were run and sampled every 1,000 generations. The first 25 percent of sampled trees were discarded as burn-in. The ultrametric tree was reconstructed in BEAST v.2.6.4 (Bouckaert et al. 2014) under the Yule speciation model. The tree was run for ten million generations and sampled every 1,000 generations. The effective sample size (ESS) was checked in Tracer v.1.7.2 (Rambaut et al. 2018). The results from two runs were combined in LogCombiner v.2.6.4 and summarized in TreeAnnotator v.2.6.4. Final trees from all analyses were visualized and edited in FigTree v.1.4.4 (Rambaut 2018). Species delimitations Species delimitation was performed using four approaches, namely automated barcode gap (ABGD), Bayesian implementation of Poisson Tree Processes model (bPTP), the Generalized Mixed Yule Coalescent model (GMYC), and cluster analysis in the BOLD system. Firstly, the ABGD method was performed at https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb. html, using the distance matrix based on the KimuraTwo-Parameter (K2P) model exported from MEGA X. All settings were left as default except the relative gap width, which was set as 1.0. Secondly, the bPTP analysis was conducted on the bPTP online server (https://species.h-its.org). The ML tree from IQ-TREE was used as the input tree. The analysis was run for 500,000 MCMC generations, 500 of thinning, and burn-in as 0.1. The run was carried out using the rooted tree with the outgroup removed. Thirdly, the GMYC was performed using the ultrametric tree from BEAST v.2.6.4 as an input. The analysis was conducted using the “splits” package in the R program (Ezard et al. 2021). For the fourth analysis, the sequences deposited into BOLD were analyzed in cluster sequence analysis using the BOLD alignment option and pairwise distance model. In cluster sequence analysis, Refined Single Linkage algorithm (RESL) compared new sequences with sequences already in the database in order to identify operational taxonomic units (OTUs). It initially applied 2% sequence divergence as a threshold, then used a pattern of divergence to refine OTUs. Each OTU was then assigned a Barcode Index Number (BIN). The final consensus species delimitation was assigned based on the results of six methods, including the OTUs from the four species delimitation approaches, phylogenetic clades with ≥ 70 bootstrap support in ML tree and ≥ 0.95 posterior probability in BI trees (Huelsenbeck and Hillis 1993; Larget and Simon 1999), and morphological identification. The consensus results of four out of the six methods were used as criteria to determine final species boundaries (García-Melo et al. 2019). DNA barcode analysis Newly obtained sequences in this study and sequences from previous studies available in the BOLD were analyzed by several tools in the BOLD system: distance analysis, barcode gap analysis (BGA), diagnostic character, and sequence composition. The alignment option used in each tool was the MUSCLE alignment algorithm. Both Kimura Two-Parameter (K2P) and pairwise distance models were applied in each analysis. Unfortunately, some sequences from previous studies were unlabeled or were labeled under other genera. Those sequences could have biased the results, so they were excluded from some analyses. Alternatively, intraand interspecific genetic distances page 4 of 14Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan were calculated in MEGA X using Kimura TwoParameter (K2P) and p-distance models. RESULTS Morphological study Leech specimens from southern Thailand that were available for morphological identification were classified into three morphological groups. (1) Hirudinaria thailandica: this species had a green ventral surface, ejaculatory ducts that were inserted medially into the atrium, and common oviduct that was near the female gonopore. (2) Hirudinaria bpling (Fig. 2A, B): this species showed a green ventral surface similar to H. thailandica, but differed in internal reproductive organs. Its ejaculatory ducts inserted antroventrally into the atrium and its common oviduct opened near the vaginal caecum. (3) The red ventral group consisted of species in the H. manillensis complex. The H. manillensis complex could be clearly identified by its red ventral side with two lateral stripes (Fig. 2C, D). In this study, we failed to obtain specimens of another species in the genus, Hirudinaria javanica. Based on our knowledge, this species is distributed in the northeastern region of Thailand, and its morphology is very unique. Whereas male and female gonopores are separated by five annuli in their congeners, gonopores of H. javanica are separated by seven annuli. Furthermore, the ventral color of this species is yellowish brown. Phylogenetic analyses Eighty-seven DNA sequences were newly obtained in this study. When including 22 sequences from previous studies, the total length of sequences was 657 bp, and contained 47 unique haplotypes. The final aligned matrix consisted of 499 conserved sites, 158 variable sites, and 146 parsimony-informative sites. Phylogenetic trees from three approaches showed similar topologies in major clades (Figs. 3, S1, S2); therefore, only the BI tree is shown in figure 3. All four nominal species were monophyletic and had high support values (ML bootstrap > 70; Bayesian posterior probability of BI and BEAST > 0.95). The genus Hirudinaria was divided into two principal clades. The first clade consisted of H. thailandica and H. bpling. Within the H. thailandica clade, there were specimens collected from locations above approximately 9 degrees north latitude (square symbols in Fig. 1), together with leeches from the Chao Phraya and Khorat Plateau basins. One sequence (JQ738405) that was previously identified as ‘H. bpling’ by Chong et al. (2014) from Malaysia was also included in this clade. The H. bpling clade was divided into two subclades, although this subdivision was supported solely by the BEAST analysis. One H. bpling subclade included samples from Trang (A134, A135), Satun (A72, A74), Chumphon (A78), and Nakhonsithammarat (A68) provinces, and specimens from Malaysia. Another subclade consisted of specimens from Suratthani (A131), Nakhonsithammarat (A132), Krabi (A76), Phuket (PK), and Phang-Nga (A136, HIR435) provinces. The latter subclade also included sequences of the type series of the species from Phan-Nga Province (JQ646012 and JQ646013; Phillips 2012). The second clade consisted of H. javanica and four deeply divergent clades of leeches that were morphologically identified as H. manillensis. The clade of H. javanica was placed at the base of the tree. We failed to obtain H. javanica from the south of Thailand. All specimens in this study were collected from northeastern Thailand except for one sequence Fig. 2. External morphology of living Hirudinaria leeches from southern Thailand. (A) dorsal and (B) ventral sides of H. bpling from Satun Province (C) dorsal and (D) ventral sides of H. manillensis 3 from Songkhla Province. Scale bar = 1 cm. page 5 of 14Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan Fig. 3. BEAST ultrametric tree of buffalo leeches genus Hirudinaria. Leeches from southern Thailand are highlighted in bold. Numbers on nodes are bootstrap values from ML tree generated by IQ-TREE, Bayesian posterior probability from BI tree generated by MrBayes, and from ultrametric tree generated by BEAST, respectively. Black bars indicate morphological identification (MORPHO) and delineated OTUs suggested by four species delimitation approaches (GMYC, bPTP, BIN, and ABGD). Grey bars indicate samples that were not available for morphological identification. page 6 of 14 Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan from Malaysia. Four lineages of H. manillensis were grouped together with high support values. 1) The H. manillensis 1 clade included specimens from northeast Thailand, four provinces in southern Thailand (Suratthani [A131], Nakhonsithammarat [A68], Phuket [PK], and Udonthani [A41]), and sequences from previous studies, including one sequence of H. manillensis (GQ368746) from Thailand, two sequences (JQ738400-1) that had been morphologically identified as Poecilobdella granulosa by Chong et al. (2014), and two sequences (SK3-018, SR2-018) from Sakon Nakhon and Surin provinces in northeast Thailand. 2) The H. manillensis 2 clade was a singleton from Puerto Rico (AY425449). 3) The H. manillensis 3 clade consisted of specimens from Phatthalung (A69) and Songkhla (A71) provinces of southern Thailand, a sequence from Malaysia (FJ610329), and two sequences from Vietnam (GQ368747-8). 4) Finally, the H. manillensis 4 clade included sequences of H. manillensis from India (KT693106-8), which we did not examine morphologically in this study. Species delimitations The ABGD method, BOLD system, and bPTP method showed consistent results of seven OTUs, including H. thailandica, H. bpling, four lineages of H. manillensis, and H. javanica (Figs. 3, 4). On the other hand, the GMYC approach revealed two additional OTUs, one in the H. thailandica clade and another in the H. bpling clade (Fig. 3). These additional OTUs correspond to the subclades in our phylogenetic trees. However, the morphological methods could successfully distinguish only four groups, which were H. thailandica, H. bpling, the red ventral group (H. manillensis 1 and H. manillensis 3), and H. javanica. Across the multiple approaches, the consensus regarding these OTUs revealed four species of Hirudinaria leeches from southern Thailand, namely H. thailandica, H. bpling, and two putative species of H. manillensis (H. manillensis 1 and H. manillensis 3). DNA barcode analysis DNA barcode analysis was done using seven OTUs, as suggested by the species delimitation analyses in the previous section (Fig. 3). The BIN discordance analysis showed five clusters containing more than one species (BOLD:AAM3860, BOLD:AAX1213, BOLD:AAM3859, BOLD:ACQ2297, BOLD:ACB6877). However, all discordances were due to the species labels in the BOLD system. Some samples were labeled as Hirudinaria sp. and some were under the genus Poecilobdella. No insertions, deletions, or stop codons Fig. 4. Result from ABGD analysis showing the stability at 7 OTUs. 10 20 30 40 50 0.0010 0.0017 0.0028 0.0046 0.0077 0.0129 0.0215 0.0359 0.0599 Number of OTUs Recursive Partition Initial Partition 0 Prior intraspecific divergence (P) page 7 of 14Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan were detected, indicating no pseudogenes in our dataset. Average genetic distances between species are shown in table 1. Genetic distance ranged from 0.11 to 0.65% within species (average of 0.36%), and 3.72 to 14.36% between species (average of 10.41%). The highest genetic divergence was that between H. bpling and H. manillensis 1, and the lowest divergence was between H. manillensis 1 and H. manillensis 2. The mean genetic distance among the four H. manillensis lineages was 7.28%. For the BOLD analysis, the warning signal, which indicates the excess of intraspecific distance over the genetic distance between nearest neighbors (NN), was not detected. All pairwise comparisons had both maximum and mean intraspecific distances less than the distance to NN (Fig. 5). The barcoding gap between intraand interspecific individuals was between 1.54 and 2.88% (Fig. 6). In the barcode gap analysis, results based on pairwise distance and K2P model were similar, therefore only results based on K2P model are presented in table 1. DNA diagnostic characters Average nucleotide base compositions were T = 41.4%, C = 15.5%, G = 16.2%, A = 26.9%. The average GC contents were 41.1%, 42.8%, and 11.0% in codon position 1, 2, and 3, respectively. The DNA diagnostic character analysis classifies nucleotide characters into five categories: diagnostic, diagnostic or partial, partial, partial or uninformative, and invalid characters. Among these categories, the number of diagnostic characters was highest and found in all species. Partial characters were found in four species (H. bpling, H. javanica, H. thailandica, and H. manillensis 1). Diagnostic or partial, partial or uninformative, and invalid characters were not identified in any species. Samples of H. manillensis 2 and H. manillensis 4 were not included in the DNA diagnostic character analysis due to insufficient sampling and mislabeling, respectively. The summarized nucleotide characters are presented in table S2. DISCUSSION Genetic divergence between lineages In this study, species delimitation approaches identified seven species of Hirudinaria leeches. Based on genetic distance, these seven species had intraand interspecific divergence comparable to previous studies on Hirudinaria. The average COI intraspecific divergence in this study (0.36%) was similar to Chong et al. (2014) and Jeratthitikul et al. (2020) (0.2–0.8%). The range of intraspecific genetic divergence in this study (0.11–0.65%) was also similar. The genetic divergence between congeners in this study (3.71–14.36%) was comparable to that between Hirudinaria species in previous studies (4.51–14.4%) (Phillips 2012; Tubtimon Fig. 5. Genetic distance based on K2P model from BOLD. (A) comparison of maximum intraspecific distance of each species and distance to its nearest neighbor (B) comparison of mean intraspecific distance of each species and distance to its nearest neighbor. Red diagonals indicate where intraspecific distance equals distance to nearest neighbor. page 8 of 14 Zoological Studies 61:84 (2022) © 2022 Academia Sinica, Taiwan et al. 2014; Chong et al. 2014; Jeratthitikul et al. 2020). It was noted that all lineages of H. manillensis in the previous studies were considered as a single species. On the other hand, Hirudo leeches, a closely related taxon, had relatively high intraand interspecific genetic distance compared to Hirudinaria leeches in this study. Possible explanations are that Hirudo species have a wide distribution range across continents, and specimens were collected throughout their ranges (Siddall et al. 2007; Trontelj and Utevsky 2012; Kutschera and Elliott 2014; Wang et al. 2022). In addition, there is evidence that Hirudinaria leeches in Thailand and neighboring countries are directly transported by the leech trade (personal communication with local traders) and/or by attaching to their mammal host during animal trading (GIAHS 2021). High dispersal ability thus leads to low intraspecific genetic diversity in some Hirudinaria species such as H. manillensis clades 1 and 3. Barcoding gaps between intraand interspecific genetic distances were 1.54–2.88% in this study. This was very narrow compared to other leeches and annelid species (Siddall et al. 2007; Martinsson et al. 2017; Ye et al. 2017; Prantoni et al. 2018). Although several studies suggested a 10% threshold of genetic divergence to Fig. 6. K2P pairwise comparisons generated from MEGA X showing barcoding gaps of 1.54 to 2.88% between intraand interspecific genetic distances of buffalo leeches genus Hirudinaria. Table 1. Average genetic divergence matrix based on uncorrected K2P distance of 657-bp COI gene fragment sequences (% ± S.E.) of buffalo leeches Hirudinaria spp. Interspecific genetic divergences are shown below the diagonal, and intraspecific distances within each taxon are shown in bold Taxa 1. 2. 3. 4. 5. 6. 7. 1. H. bpling 0.43 ± 0.01 2. H. thailandica 4.36 ± 0.01 0.65 ± 0.03 3. H. manillensis 114.36 ± 0.01 13.20 ± 0.03 0.25 ± 0.02 4. H. manillensis 2 13.63 ± 0.03 12.43 ± 0.04 3.72 ± 0.07 N/A 5. H. manillensis 3 12.84 ± 0.02 12.41 ± 0.02 9.29 ± 0.03 8.65 ± 0.06 0.32 ± 0.02 6. H. manillensis 4 11.55 ± 0.02 10.10 ± 0.02 8.01 ± 0.04 7.69 ± 0.02 6.31 ± 0.02 0.11 ± 0.06 7. H. javanica 12.99 ± 0.02 12.38 ± 0.04 10.612 ± 0.08 10.51 ± 0.11 12.20 ± 0.09 10.46 ± 0.05 0.40 ± 0.07 page 9 of 14Zoological Studies 61:84 (2022)