Exploring Barbronia species diversity and phylogenetic relationship within Suborder Erpobdelliformes (Clitellata: Annelida)

Introduction

Barbronia is a genus of predatory freshwater leeches belonging to the family Salifidae (Erpobdelliformes, Clitellata, Annelida) with only three species have previously been reported in China: Barbronia weberi, Barbronia zhejiangica Yang, 1996, and Barbronia yunnanensis Yang, Wang & Zhang, 1997 (Yang, 1996; Yang, Wang & Zhang, 1997). The former one is a well-known introduced leech species, however, the latter two are rarely known endemic Barbronia species. Barbronia weberi was originally described by Blanchard (1897) in Indonesia, but now it is broadly found in China (Ta-Hsiang, 1974; Zhao et al., 2020), India (Bandyopadhyay & Mandal, 2005; Chandra & Mahajan, 1971; Ghate, 1991), Myanmar (Eriksen et al., 2022), United States (Rutter & Klemm, 2001; Sawyer & Sawyer, 2018), Brazil (Pamplin & Rocha, 2000), Mexico (Garduño-Montes de Oca et al., 2016; Oceguera-Figueroa, León-Règagnon & Siddall, 2005), United Kingdom (Sawyer, 1986; Sawyer & Sawyer, 2018), Germany (Kutschera, 2004; Nehring, 2006), Hungary (Ludanyi et al., 2019), Italy (Genoni & Fazzone, 2008), Spain (Pavluk, Pavluk & Rasines, 2011), Netherlands (Van Haaren et al., 2004), Australia (Govedich et al., 2003; Govedich, Bain & Davies, 2002), New Zealand (Mason, 1976), and South Africa (Nakano & Nguyen, 2015). Barbronia weberi has achieved its current wide distribution not only through human activities (Govedich et al., 2003), but also through its ability to reproduce cocoons asexually, without the need for cross-fertilization (Sawyer, 2020). Moreover, B. weberi is the second intermediate host of the parasitic trematode Australapatemon niewiadomski (Blasco-Costa, Poulin & Presswell, 2016).

The taxonomic status of the widely distributed invasive species B. weberi needs to be evaluated on a global scale (Oceguera-Figueroa et al., 2011), which involved B. weberi formosana, B. wuttkei, B. arcana, and Barbronia sp. Molecular identification has become a valuable complement to morphological taxonomy over the past two decades. A partial mitochondrial gene COI sequence has been used to recognize Barbronia species. Barbronia wuttkei was originally described as Erpobdella wuttkei based solely on morphological traits (Kutschera, 2004), however, it was later transferred to the genus Barbronia but retained as a valid name by Grosser & Trontelj (2008) based on COI phylogenetic evidence. The genetic analysis using the COI sequence-based species delimitation approach, Poisson Tree Process (PTP), has led to the initial assumption that Barbronia wuttkei and B. arcana should be classified under the taxon B. weberi (Klass et al., 2021). One Barbronia sp. (GenBank accession: MN503261), as well as another specimen (GenBank accession: MF458701, marked as Erpobdella sp.), are recently recognized as B. gwalagwalensis on basis of COI data (Klass et al., 2021).

The phylogeny of Barbronia was initially inferred solely based on 18S ribosomal DNA sequences alone, B. weberi was found to be closely related to a clade of two erpobdellids: Dina lineata (synonym: Erpobdella lineata) and Erpobdella octoculata (Trontelj, Sket & Steinbrück, 1999). Although an attempt was made to amplify a fragment of the mitochondrial marker 12S rDNA from the B. weberi specimen by Trontelj, Sket & Steinbrück (1999), it failed using the primer pair 12S-A/12S-B (Siddall, 2002). Consequently, only the unique 18S rDNA sequence of B. weberi was used in subsequent phylogenetic study (Borda & Siddall, 2004a; Oceguera-Figueroa, León-Règagnon & Siddall, 2005). A comprehensive phylogenetic analysis was conducted using combined data from partial sequences of 18S rDNA, 28S rDNA, 12S rDNA, and cytochrome c oxidase subunit I (COI) published by Borda & Siddall (2004b). This study confirmed that the Barbronia species (B. weberi, B. weberi formosana, B. gwalagwalensis, and an undescribed Barbronia species) formed a monophyletic group closely related to Linta be. Another phylogenetic result supported that B. weberi and B. arcana were part of the Salifidae clade, a sister group of Erpobdellidae (Oceguera-Figueroa, León-Règagnon & Siddall, 2005). An expanded dataset including six Barbronia specimens represented all mentioned Barbronia species and B. wuttkei, forming a monophyletic group that was a sister group to a clade of Salifa and Linta species (Oceguera-Figueroa et al., 2011). Recently, Barbronia represented by three taxa (Anderson, Braoudakis & Kvist, 2020), and Barbronia weberi only (Tessler, Siddall & Oceguera-Figueroa, 2018), are closely related to a clade that includes Odontobdella blancharidi and Mimobdella japonica.

The implementation of next-generation sequencing (NGS) technologies has led to a substantial increase in the number of whole mitochondrial genomes in the GenBank database (Kuang & Yu, 2019), and has gradually become a mainstay of species identification and phylogeny research (Franco-Sierra & Diaz-Nieto, 2020; Kortsinoglou et al., 2020). To date, however, only about 20 sequences of various gene fragments from Barbronia specimens are available in the GenBank database. Eight of these sequences are the mitochondrial COI barcodes, six partial 18S rDNA sequences, four short 28S rDNA sequences, and only two 12S rDNA sequences, representatively. Neither a complete mitochondrial genome of Barbronia species nor even the family Salifidae has been published. This gap significantly limits our understanding of either the species identification or phylogenetic relationships of Barbronia within Salifidae.

In the present study, the whole mitochondrial genome, as well as the 18S rDNA and 28S rDNA of Barbronia cf. gwalagwalensis, are assembled and annotated; the species boundary of Barbronia species is evaluated using two kinds of single-locus COI based species delimitation methods; the phylogenetic relationships of Barbronia within the Suborder Erpobdelliformes are conducted based on both mitochondrial and nuclear data.

Materials and Methods

Results

General features of mitochondrial genomes and nuclear rDNAs

The newly assembled mitochondrial genome of B. cf. gwalagwalensis is AT-rich (71.9%) circular mapping molecule (total length, 14,847 bp), with an average coverage of 209-fold. The components of B. cf. gwalagwalensis mitogenome consist of 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), 2 rDNAs (12S and 16S rDNA), and a possible control region (451 bp, the longest non-coding reign is located between tRNA R and tRNA H). All mitochondrial genes of Barbronia are encoded on the same strand, and gene order are generally consistent with previously published erpobdellid mitogenomes (Fig. 1 and Table 3). Among 13 Barbronia PCGs, both NAD2 and NAD5 genes are inferred to use ATT as an initiation codon, the COX3 gene is initiated with ATA, and the remaining genes used ATG as a start codon. The predicted secondary structures of 22 tRNAs had a similar clover leaf shape (Fig. 2), and these tRNAs range in size from 60 to 69 bp, with the shortest tRNA gene being tRNA L and the longest one being tRNA Q. The complete length of nuclear 18S rDNA and 28S rDNAs are 1,876 and 2,863 bp, respectively, the former one is identical but longer than the partial 18S rDNA sequence (AY786462) of B. gwalagwalensis.

Table 3 :

General features of the mitochondrial genomes of B. cf. gwalagwalensis and other specimens.

Taxa	Size (bp)	%GC	CDSs	tRNAs	rRNAs	Accession
B. cf. gwalagwalensis	14,847	28.1	13	22	2	OQ339201*
Erpobdella testacea	14,495	26.9	13	20	2	MT584166
Erpobdella japonica	14,725	27.9	13	22	2	MF358688
Erpobdella octoculata	13,035	27.9	12	20	2	MT410851
Erpobdella octoculata	15,580	27.5	13	22	2	OM257408
Erpobdellidae sp.	14,746	30.2	13	22	2	MW43152
“Erpobdella octoculata”	14,407	28.4	13	22	2	KC688270

DOI: 10.7717/peerj.17480/table-3

Note:

The specimen (KC688270 or NC_023927, named as “Erpobdella octoculata”) is likely to be a misidentified taxon or the outcome of sequencing contamination. The accession numbers highlighted with asterisks in bold indicate newly added data.

The outcome of single-locus-based species delimitation

The B. weberi groups and the B. gwalagwalensis group listed in Table 1 correspond to two main well-supported (PP > 0.99) clades in COI phylogenetic analysis (Fig. 3). All members of the B. weberi group are consistently recognized as a well-supported PSH, including the previous described taxa B. weberi, B. weberi formosana, B. arcana, Barbronia wuttkei (clade A in Fig. 3). The primary species hypothesis (PSH) of B. weberi groups proposed by bPTP are consistent with bGMYC, but the difference between two results is the species delimitation of individuals listed in the B. gwalagwalensis group (clade B in Fig. 3). In the bPTP results, the five individuals of the B. gwalagwalensis group in clade B are divided into two PSHs. One well supported PSH comprise only B. cf. gwalagwalensis, and the other moderately supported PSH encompass the paratype of B. gwalagwalensis (GenBank accession number: AY786455) and the remaining specimens (MF458701, MN295405 and MN503261, Fig. 3 and Table 1). However, the five individuals of the B. gwalagwalensis group are either split into four well-supported PSHs or collectively classified as one moderate PSH in the bGMYC analysis (Fig. 3). In the four well-supported PSHs, one PSH contains MF458701 and AY786455, and each of the remaining PSHs is formed separately from MN295405, MN503261and OQ339201. Barbronia cf. gwalagwalensis (OQ339201) alone is consistently recognized as a separately well-supported PSH in both bGMYC and bPTP analyses (threshold 0.95–1), although it is also somewhat classified with the paratype B. gwalagwalensis (AY786455) in a moderately supported PSH (threshold 0.5–0.9) by bGMYC. Moreover, the uncorrected p-distance of COI sequences between B. cf. gwalagwalensis (OQ339201) and the paratype B. gwalagwalensis (AY786455) is 4.4% (Fig. 4), which is substantially higher compared to the p-distance (ranged from 0.2% to 1.2%) among four individuals in clade B (Fig. 4) collected from Myanmar (MN295405), South Korea (MN503261), South Africa (AY786455) and France (MF458701) (Fig. 4 and Table 1). The p-distance between above four individuals and B. cf. gwalagwalensis ranged from 3.8% to 4.5%. In contrast, the p-distances of seven B. weberi individuals in clade A (Fig. 3) sampled from South Africa, the United States, Mexico, Costa Rica, and Germany are less than 0.8%, and the interspecific uncorrected p-distance between the B. weberi group and the B. gwalagwalensis group (clade A and clade B, respectively, Fig. 3) ranged from 6.4% to 8.7%. Considering the result of species delimitation and the sufficiently large p-distance mention above, the current specimen (OQ339201) is treated as B. cf. gwalagwalensis.

Phylogenetic relationships within Erpobdelliformes

The best-fit model for each partition, GTR + I + G, was selected using PartitionFinder 2. Three Haemopis leeches belonging to Haemopidae (Hirudiniformes) were used as outgroups both in Bayesian Inference (BI) analysis and the maximum likelihood (ML) analysis. Both the BI tree and ML tree topologies show generally well-supported for four main monophyletic groups (PP = 1.00, and BS > 95%), despite differences in the placement of the four clades (Fig. 5). The first clade (PP = 1.00, and BS = 99%) is Salifidae represented by nine specimens from five genera (Barbronia, Mimobdella, Odontobdella, Linta, and Salifa), the second clade (PP = 1.00, and BS = 99%) is Orobdellidae consisted of 18 Orobdella species, the third one (PP = 1.00, and BS = 95%) is Gastrostomobdellidae sensu stricto on behalf of three Gastrostomobdella species, and the last one is a well-supported Erpobdellidae clade (PP = 1.00, and BS = 100%) of 18 Erpobdella species. Barbronia is a well-supported (PP = 1.00, BS = 100%) monophyletic group, and its closely related clade consisted of Mimobdella and Odontobdella. The genus Salifa is recovered as paraphyletic with the inclusion of Linta be. Gastrostomobdellidae sensu stricto (containing genus Orobdella only), but not Gastrostomobdellidae sensu lato (Orobdella + Gastrostomobdella), is a well-supported monophyletic group. The monophyly of Erpobdelliformes was well supported in Bayesian Inference (BI) analysis, with a posterior probability (PP) of 1, but it was not corroborated by the maximum likelihood (ML) analysis.

Discussion

To classify Barbronia cf. gwalagwalensis as a distinct species, additional evidence is needed since cryptic speciation is common among clitellates (Erséus & Gustafsson, 2009; Liu et al., 2017; Martinsson & Erséus, 2021). Van Haaren et al. (2004) identified one Dutch specimen as belonging to the B. assiuti/weberi complex, acknowledging that distinguishing the species or identifying it based solely on the morphological characteristics of this single specimen is also not feasible, which highlights the need for careful distinction between individual sample identities and broader species classification. The evidence for new species should encompass both multiple loci-based genetic approaches and morphological analyses from numerous specimens (Dellicour & Flot, 2015; Jorger & Schrodl, 2013; Stengel et al., 2022; Sukumaran, Holder & Knowles, 2021). Initially, this specimen was morphologically identified as the common species B. weberi due to the observation of accessory gonopores (the anterior and posterior ones separately close to the male and female gonophores). The presence of two accessory gonopores in this specimen rules out the possibility of it being B. zhejiangica or B. yunnanensis (Yang, Wang & Zhang, 1997). Molecular results show that this specimen is more closely related to B. gwalagwalensis and the Korean Barbronia specimen (KF966549), the latter of which has been tentatively assigned to Barbronia cf. zhejiangica by Klass et al. (2021). Interestingly, accessory gonopores have also been observed in B. gwalagwalensis (Westergren & Siddall, 2004). However, there are no previous records of B. gwalagwalensis in China. Unfortunately, a thorough morphological examination of this small specimen is unfeasible due to the complete utilization of its body for DNA extraction, a necessary step for constructing a high-quality library in next-generation sequencing. Therefore, the current classification for this specimen is as the taxon B. cf. gwalagwalensis, rather than describing it as a new species. This situation highlights the potential for B. gwalagwalensis to have a broader distribution than previously thought. Initially identified in South Africa, this species has also been found in Myanmar, as noted by Klass et al. (2021). Further studies, involving more samples from across Asia and Europe and detailed morphological analysis, are required to draw definitive conclusions.

Over last two decades, species identification based on molecular data has indeed provided a valuable complement to morphological taxonomy (Mahadani et al., 2022), aided by the increasing availability of genetic techniques (Dellicour & Flot, 2015). In the current study, seven specimens of B. weberi constituted a well-supported clade in the COI phylogenetic results, this clade is also recognized as one valid species both in bGMYC and bPTP analyses (Fig. 3). The current species delimitation results support that the previously proposed taxa (B. arcana, B. weberi formosana and B. wuttkei) are synonyms of B. weberi rather than valid species, which is generally consistent with previous studies (Klass et al., 2021; Oceguera-Figueroa et al., 2011). The COI sequence similarity among the seven individuals within the B. weberi group (clade A in Fig. 3), collected from South Africa, the United States, Mexico, Costa Rica, and Germany are nearly identical. The lack of molecular differentiation (small p-distance) between them can be explained by the large effective population size but low rate of mutation, caused by both the relatively large distribution of invasive B. weberi and the potentially low speciation rate (i.e., reproduce cocoon without cross-fertilization as mentioned in the introduction). However, it is important to note that a gene tree is not always consistent with a species tree when population sizes are large and speciation rates are high (Dellicour & Flot, 2015), single locus-based methods may not be sufficient for delimiting specimens listed in the group of B. gwalagwalensis, especially considering the relatively high genetic diversity between specimens sampled from southern Africa, France, Korea, and China (Fig. 4). On the basis of all available Barbronia data listed in Table 1, all Barbronia specimens tend to lump into one taxon by applying multiple loci-based species delimitation method. This is likely due to insufficient data. In this scenario, utilizing data from both the mitochondrial genome and nuclear loci collected from each specific B. gwalagwalensis specimen (listed in Table 1) is expected to provide a more promising solution for resolving the species delimitation issue in B. gwalagwalensis, rather than relying solely on the partial mitochondrial gene COI.

Mitochondrial genomes and full-length nuclear rDNA have been now routinely applied to assess species boundaries and deep relationships in many phylogenetic studies due to obtaining sufficient data efficiently through NGS (Jia et al., 2023; Moreno-Carmona et al., 2023; Prada et al., 2023). However, a complete mitochondrial genome of any Barbronia species or of the family Salifidae remains unpublished, highlighting a significant gap in our understanding of their systematics. This lack of such data presents a challenge in accurately determining species delimitation and phylogenetic relationships within these groups, underscoring the need for more comprehensive genomic research. In the current study, the first complete mitochondrial genome of salifid leeches is assembled, and compared with whole or incomplete mitochondrial genome sequences of erpobdellids (Fig. 1 and Table 3), representing three nominal species and one unrecognized species in the GenBank. The specimen named Erpobdella octoculata (KC688270 and NC_023927) is likely to be a misidentified specimen or the outcome of sequencing contamination (Oceguera-Figueroa et al., 2016). Clearly, it is far from insufficient using these data to estimate either the species delimitation of Barbronia species or the phylogenetic relationships of Barbronia within Erpobdelliformes (Arhynchobdellida, Hirudinea). Therefore, four conventional molecular markers used in many previous studies of Erpobdelliformes are collected and analyzed, including full-length 18S and 28S rDNA sequences and COI and 12S rDNA sequences extracted from the newly acquired mitochondrial genome. The four well-supported clades, i.e., Salifidae, Orobdellidae, Gastrostomobdellidae sensu stricto, and Erpobdella, are clearly distinguished in the current molecular phylogenetic results. The Barbronia clade is sister to a clade consisting of Odontobdella blanchardi and Mimobdella japonica with well supported in the current study (Fig. 5), and the monophyly of the Barbronia is congruent with previous research (Anderson, Braoudakis & Kvist, 2020; Klass et al., 2021; Nakano et al., 2018; Nakano & Nguyen, 2015; Oceguera-Figueroa et al., 2011). However, the Salifa genus represented by two species is not monophyletic in the current study and in Nakano et al. (2018) result, which is not consistent with previous studies (Klass et al., 2021; Nakano & Nguyen, 2015). The Orobdellidae clade is well supported, which is congruent with previous studies (Nakano, 2016a, 2016b, 2022; Nakano et al., 2018; Nakano & Nguyen, 2015; Nakano, Ramlah & Hikida, 2012). Additionally, previous studies have shown that Erpobdella forms a well-supported clade within the Erpobdellidae family. However, incorporating more species from Gastrostomobdella, rather than just one, provides stronger support for a closer relationship between the Erpobdella clade and the Gastrostomobdellidae family (Nakano et al., 2018; Nakano & Nguyen, 2015). Despite these insights, the robust phylogenetic positions of four families within Erpobdellidae remains uncertain in the current study, these findings underscore the intricate evolutionary relationships within the Erpobdellidae family, indicating the need for further research using more data to fully elucidate the phylogeny of these leech genera. To solve these problems, further studies involving more sampled specimens, or at least simultaneously obtaining more molecular data (i.e., mitochondrial genomes and rDNAs) from a few valuable specimens using NGS are urgently needed. The forthcoming systematics study of B. gwalagwalensis and related species is expected to benefit significantly from the methods employed in our current research, particularly in acquiring mitochondrial genomes and full-length rDNA sequence. This approach effectively bridges gaps between sequences amplified with different primers, notably improving results in cases where conventional primers are less effective (Siddall, 2002; Trontelj, Sket & Steinbrück, 1999).

Conclusions

The first complete mitochondrial genome, full-length 18S and 28S rDNAs of salifid leeches within the family Salifidae are provided. The species delimitation results based on COI data supported the view that the previously proposed taxa (Barbronia arcana, B. weberi formosana, and B. wuttkei) are synonyms of B. weberi. However, the taxonomic status of B. gwalagwalensis and B. cf. gwalagwalensis need to be further studied, necessitating extensive fieldwork across Asia and Europe. Barbronia formed a well-supported clade, including B. weberi, B. gwalagwalensis, and B. cf. gwalagwalensis, and Barbronia is sister to a well-supported clade comprising Odontobdella and Mimobdella. They constitute a strongly supported monophyletic group and are incorporated into the Salifidae family. Additionally, the families Salifidae, Gastrostomobdellidae sensu stricto (containing only the genus Orobdella) and Orobdellidae are all well-supported monophyletic clades. According to the results of this study, the strategy of obtaining both whole mitochondria and nuclear markers from extensively sampled salifid species using NGS is expected to elucidate the species diversity of B. gwalagwalensis and the evolutionary relationship of Salifidae.