Leech breach: a first record of the invasive freshwater leech Helobdella europaea (Hirudinea: Glossiphoniidae) in Fiji

Study area and sample collection

A kick-netting technique (Stark et al. 2001) was used to collect leech specimens (n = 50) for genetic analyses from the Koroboya village stream system, Ba Province, Viti Levu, Fiji in 2020 (Fig. 1). The specimens were preserved whole in 80% ethanol.

Fig. 1.

Collection locations and representative specimens of H. cf. europaea from Fiji. (a) Koroboya village stream system and wider Ba Province watershed annotated with coloured sampling locations and dates, (b) the location of Ba Province (red) in Fiji and (c) Fiji relative to the South Pacific region. Fijian specimen features shown include a (d) dorsal view with five rows of distinct black-tipped papillae; (e) cone-shaped dorsal papillae, dark longitudinal stripes and white chromatophores; and (f) ventral view with fully developed young attached – (J), a pair of eyes with (Cs) caudal sucker. Insets (a–c) were generated using QGIS v 3.18.3-Zürich and open source data obtained from the Humanitarian Data Exchange (https://data.humdata.org/dataset/cod-ab-fji) and Natural Earth (https://www.naturalearthdata.com/).

Morphological analyses

Additional specimens obtained from Bavu village in the Nadroga-Navosa province (year 2015) and Ba province (years 2015, 2019, and 2020, n = 52) were morphologically identified as H. europaea following Kutschera (2004), Lai et al. (2009) and Kutschera et al. (2013). Examinations involved measuring specimen sizes (body length at rest), observing form (body shape, presence of a convex dorsum, anterior sucker, posterior (caudal) sucker and pair of eyes with annulation along the axis of the body), colour and pattern.

Phylogenetic analyses

Total DNA was extracted from individual whole leech specimens (n = 32/50 specimens collected) using a CTAB (cetyltrimethylammonium bromide) chloroform/isoamyl alcohol protocol (Adamkewicz and Harasewych 1996), modified as per Lal et al. (2016) at the University of the South Pacific, Suva, Fiji. Purified DNA was quantified with a QUBIT flourometer (Life Technologies, MA, USA) using the Broad Range DNA assay. Subsequently, 16 samples were further processed for Polymerase Chain Reaction (PCR) and Sanger sequencing at the School of Natural Sciences at Massey University in New Zealand.

A 658 bp mitochondrial cytochrome c oxidase subunit I (cox1) DNA fragment was targeted for amplification using the universal primers (Folmer et al. 1994): LCO1490 (5’-GGT CAA CAA ATC ATA AAG ATA TTG G-3’) and HCO2198 (5’-TAA ACT TCA GGG TGA CCA AAA AAT CA-3’). Amplification was carried out on a Whatman T1 (Whatman Biometra, Göttingen, Germany) thermal cycler, with each reaction using EmeraldAmp GT Master Mix (Takara, Kusatsu, Shiga, Japan), 10 pmol of each primer and 2 ng template DNA. Cycling conditions were as follows: 94°C for 3 min; 35 cycles of 94°C for 30 s, 55°C for 30 s, 72°C for 40 s, and a final extension step at 72°C for 5 min with a hold at 10°C. PCR products were assayed using 1.5% agarose gel electrophoresis before clean-up with a SAP-Exo kit (Jena Bioscience, Jena, Germany). Bi-directional sequencing was performed on an ABI 3730 DNA Analyzer (Applied Biosystems, Waltham, MA, USA) following standard instructions.

Chromatograms were edited using CodonCode Aligner ver. 3.7.1–6.0.2 (CodonCode Co., Centerville, MA, USA) in the Geneious v.9.0.4 package (Kearse et al. 2012). Gaps in the sequences were treated as missing data, and final alignment for sequences was by MUltiple Sequence Comparison by Log-Expectation (MUSCLE, (Edgar 2004). The final alignment was subjected to substitution model testing using jModelTest (Posada 2008), with Bayesian Information Criterion (BIC) scores including the negative log likelihood (−lnL) and BIC difference (delta parameter) used to select the optimal model. Edited and aligned forward and reverse sequences were then subjected to a megaBLAST search to retrieve matching cox1 Helobdella sequences from the NCBI GenBank (Table 1). Sequences belonging to Placobdella phalera and Glossiphonia elegans were used as outgroups as per Lai et al. (2009).

Bayesian inference was used for phylogenetic reconstruction in the MrBayes v3.2 package (Huelsenbeck and Ronquist 2001; Ronquist and Huelsenbeck 2003; Ronquist et al. 2012). Parameter settings for the data block were set as follows: lset applyto = (all), nst = 2, rates = gamma; prset applyto = (all), lset applyto = (all). Each analysis incorporated two runs of 100 000 000 generations, with four independent chains for each run, a temperature of 0.10 for the heated chains, sampling frequency of 1000 and burn-in fraction of 20%. The burn-in threshold was selected on the basis that both independent runs had achieved stationarity (i.e. stable log likelihood values reached for all sampled trees, indicated by the average standard deviation of split frequencies <0.01). Convergence was also assessed using Tracer v.1.6 independently (Rambaut et al. 2003). Final trees were retained by selecting only those post-burn in trees producing the highest individual (p) and cumulative posterior (P) probabilities during Markov Chain Monte Carlo (MCMC) computations, of p = 0.000 and P ≥ 0.995 respectively. A consensus tree was then constructed from the final credible sets of trees using the strict consensus method in Dendroscope 3.5.7 (Huson et al. 2007). All phylograms were visualised, inspected and edited in FigTree (https://tree.bio.ed.ac.uk/) v.1.4.2 (Rambaut 2014). Net evolutionary divergence estimates between the different Helobdella sp. taxa, as well as among H. europaea sampled from different collection sites were computed in MEGA6 (Tamura et al. 2013). Standard error estimates were obtained following a bootstrap procedure with replication.

The mitochondrial genome of H. europaea

Ten nanograms of H. europaea total genomic DNA from a Fijian specimen and 62.5 μM of the selective multiple displacement amplification (MDA) primers TATT*A*A and TAAT*T*A (* indicates exonuclease resistant phosphorothioate bonds) were denatured in 5 μL of 1× EquiPhi Reaction Buffer at 95°C for 3 min. The mix was cooled on ice for 3–5 min then added to 20 μL of 1× reaction mix containing EquiPhi29 DNA polymerase Reaction Buffer, DTT, dNTPs and EquiPhi29 DNA Polymerase as per the protocol: https://www.thermofisher.com/order/catalog/product/A39390. The reaction was incubated at 37°C for 3 h followed by an inactivation step at 65°C for 10 min.

After whole genome amplification, impurities were removed from the amplified DNA using a DNA Genomic Clean and Concentrator Column (Zymo Research). The amplified DNA was eluted from the column in 15 μL of water and the concentration of the eluate (290 ng/μL) determined by QUBIT (ThermoFisher). A MinION PCR-free library was then prepared using the Oxford Nanopore Technologies (ONT) rapid barcoding kit (SQK-RBK004) and sequenced on a MinION flongle (R9.4.1) using an ONT Mk1C MinION device running MinKNOW v21.10.8. H. europaea total genomic DNA was also sequenced without MDA on a flongle on the ONT Mk1C MinION device.

Bases were called using Guppy v5.0.17. Adaptors were trimmed using Porechop v0.2.4 (Wick 2018). Assembly of the mitochondrial genome was made using Trycycler (Wick et al. 2021). Short and low-quality reads were removed with NanoFilt https://github.com/wdecoster/nanofilt to produce two filtered read sets: (1) reads of length ≥1000 bases and (2) reads of length ≥5000 bases. BLASTn searches https://anaconda.org/bioconda/blast of read sets against the mitochondrial genome of Placobdella lamothei NC_030269.1 and against the assembled H. europaea mitochondrial genome were used to identify and count mitochondrial reads in MDA enriched and non-enriched DNA. Flye v2.9 was used to assemble contigs from read sets https://anaconda.org/bioconda/flye, Minimap2 v2.24 was used to map reads to linear contigs https://anaconda.org/bioconda/minimap2 and Geneious Prime 2022.2.1 www.geneious.com to map reads to a circular reference. Mapping to a circular reference was necessary as reads that spanned the 5′ and 3′ ends of the assembled mitochondrial genome were not well mapped using a linear mapper.

Comparing the counts of mitochondrial DNA (mtDNA) reads in amplified and unamplified DNA allowed us to calculate the amount of enrichment of mtDNA produced by selective MDA. Annotations to the mitochondrial genome were made using Geneious 9.1.8 https://www.geneious.com and these were compared to mitochondrial genome sequences for P. lamothei NC_030269.1 and Helobdella robusta AF178680.1. The mitochondrial genome sequence of H. europaea is available from the NCBI GenBank under accession number: NC_072606.1.

Ethical standards

Responsibility for the ideas expressed in this paper rest with the authors alone. The authors assert that all procedures contributing to this work comply with applicable national and institutional ethical guidelines on the care and use of laboratory or otherwise regulated animals.

Morphological analyses

Preserved adult specimens are depicted in Fig. 1d–f. Microscopic observations of morphological features revealed visible segments (Fig. 1d), a prominent ‘bumpy’ dorsal surface with distinguishing features (Fig. 1e) and a suction disc at the tail end (Fig. 1f). Other features observed included a moderately flattened ovate-lanceolate body, one pair of punctiform to triangular eyes in segment III (3rd annulus), five rows of distinct, black-tipped papillae of the dorsum, dorsal pigmentation arranged in numerous longitudinal dark grey stripes and the absence of a nuchal scute. All specimens collected from Koroboya village stream had a size range of 3–8 mm (Rashni and Brown 2021) and features were concordant with H. europaea.

Ecological characteristics

The Fijian specimens collected from Ba and Nadroga provinces inhabit slow flowing freshwater streams with heavily modified riparian vegetation. The Koroboya stream is 3–5 m wide and the channel depth range is 0.10–1.0 m with a silted streambed. The streambed was mostly bare during sampling and some areas were covered with Chara sp. (a green macroalga) and leaf litter. Large numbers of leeches were found on silt covered Chara sp. leaves and instream leaf litter dominated by Mangifera indica (mango) leaves. Previously, Helobdella leeches were found on silted Potamogeton crispus (a freshwater aquatic plant) beds in the Qalinabulu creek, Balevuto village, Nadroga.

Phylogenetic analyses

Sixteen samples of leech total DNA from Fiji amplified and were sequenced successfully, returning identical sequences. A megaBLAST search produced 100% identity matches to 12 catalogued H. europaea sequences: AY856049.1, DQ995298.1, AF329052.1, DQ995297.1, MF804537.1, AY856048.1, MG976140.1, FJ000349.1, KC904241.1, KC904242.1, MN335875.1 and DQ995304.1. Given that the 16 Fijian samples contained identical sequences, a consensus sequence (Supplementary Material 3) was generated for further analyses.

The consensus Fiji sequence was pooled with 80 Helobdella and outgroup sequences (Table 1) for alignment and substitution model testing. The Hasegawa–Kishino–Yano model with gamma-shaped rate variation and invariable sites (HKY + G + I) was determined to have the best fit from jModelTest results (−ln L = 2076.5, BIC = 5835.7 and AICc = 4491.6 [corrected Akaike information criterion]), and was subsequently used for Bayesian reconstruction of phylogenetic relationships. Two parallel runs of MrBayes were used to generate the final tree, where the final average standard deviation of split frequencies achieved was 0.0028 with average potential scale reduction factor (PSRF) for parameter values remaining at 1.000. In total, 16 002/20 002 trees were sampled from the parallel MrBayes runs and a burn-in set (20%) discarded to generate an initial subset of 12 802 trees. Subsequently, 89 highly credible trees were retained (posterior probability p ≤ 0.000 and cumulative posterior probability P ≤ 0.995) and a strict consensus final tree generated (Fig. 2).

Fig. 2.

Bayesian phylogenetic reconstruction based on 81 cytochrome c oxidase (cox1) nucleotide sequences. A consensus sequence was generated for the Fijian samples from 16 individual sequences. The tree presented is a strict consensus of 89 highly credible trees (cumulative posterior P ≤ 0.995). Node support values for major clades indicate posterior probabilities while the scale bar represents the number of substitutions per site. The major clade resolved for Helobdella europaea is highlighted in red, while the larger clade for the H. triserialis series sensu stricto as per Morhun et al. (2021) is presented in grey. Taxon images depict H. modesta (a), H. stagnalis (b), H. octatestisaca (c), H. papillata (d), H. robusta (e), H. triserialis (f) and H. europaea (g). Outgroup taxa used were Glossiphonia elegans (h) and Desserobdella phalera syn. Placobdella phalera (i). Images (b, c, f–h) are used with permissions from iNaturalist.org users mentioned in the acknowledgements section. The remaining images are adapted from www.boldsystems.org, with the following attributions: (d) and (i), no rights reserved/public domain; (a) and (e) CC BY CBG Photography Group.

Our reconstruction confirms the assignment of the Fijian leech samples to Helobdella europaea, as the consensus sequence resolved a strong monophyletic clade with H. europaea specimens from Australia, New Zealand, Spain, Taiwan, South Africa, Hungary, Ukraine and the USA. This clade has short internal branch lengths and high posterior probability node support (100). The H. europaea clade is paraphyletic within a larger grouping comprising the sister taxa H. triserialis and H. socimulcensis. Two other large paraphyletic groups were identified, highlighting relationships between H. robusta and H. papillata; and H. modesta with H. stagnalis and H. octatestisaca, respectively. All clades were resolved with high posterior probabilities.

Evolutionary divergence estimates

Relationships identified in the Bayesian tree topology were supported in the pairwise net evolutionary divergence estimates computed between taxa and collection locations (Tables 2 and 3). Pairwise divergence estimates between the Fijian H. europaea sample sequence and all other H. europaea sequences were effectively zero (Table 3), supporting conspecificity of these taxa. A minor degree of divergence was detected between the Fijian and Australian sequences (0.0015 ± 0.0011), however the Australian samples also reflected an identical degree of separation from other conspecific sequences from Spain, Taiwan, the USA, Hungary and South Africa. The largest but still small degree of separation was evident between the Ukrainian sequence (MT258557) and all others (0.0022–0.0037).

The close affinity of H. europaea with H. socimulcensis (0.0297 ± 0.0079) and H. triserialis (0.0304 ± 0.0086) provides support that these three taxa may comprise a species complex (Table 2, values in bold), which is further discussed by Morhun et al. (2021). Similarly, pairwise divergence estimates were low between H. robusta with H. papillata (0.0277 ± 0.0054), and H. stagnalis and H. octatestisaca (0.0084 ± 0.0035), respectively. Estimates were highest between H. octatestisaca and all other taxa with the exception of H. stagnalis, ranging from 0.2027 ± 0.0298 (H. papillata) to 0.2382 ± 0.0349 (H. europaea).

Mitochondrial genome of H. europaea

The genome sequence determined for H. europaea was assembled from high coverage of ONT MinION sequenced MDA mitochondrial DNA reads using Flye v. 2.9 based on the tutorial developed for Trycyler (Wick et al. 2021). Selective MDA using the protocol and primers described enriched mtDNA from 1% to more than 28% reads of the reads ≥1 kb in length. Details on enrichment and genome coverage have been provided in Supplementary Material 1 and 2.

The assembled annotated genome shown in Fig. 3 was found to be very similar in organisation to that of P. lamothei NC_030269.1 and H. robusta AF178680.1. The complete mitochondrial genome sequence is available online under accession NC_072606.1 (https://www.ncbi.nlm.nih.gov/nuccore/NC_072606.1) and as Supplementary Material 4. Unique to H. europaea was the presence of repeats in the control region. These repeats were identified in both amplified and unamplified H. europaea DNA that were sequenced on the MK1C minion device. It was not determined whether the number of repeats were variable within and between individuals of H. europaea.

Fig. 3.

Annotated mitochondrial genome of Helobdella europaea.

Morphological and ecological insights

The external morphology of samples collected in Fiji is consistent with descriptors for H. europaea reported in other studies (Kutschera 2004; Kutschera et al. 2013) and further studies examining specimen gut contents may reveal potential impacts on native aquatic invertebrate fauna. Our observations suggest that H. europaea in Fiji prefers soft-bottomed sections of streams inhabited by macrophytes (P. crispus and Chara sp.) and leaf litter.

Molecular barcoding

Resolution of the Fijian leech samples in a monophyletic clade, strong sequence identity matches and evolutionary divergence computations with catalogued H. europaea specimens confirm this species is now present in Fiji. Results reported here underscore the utility of molecular barcoding tools for identifying invasive or potentially invasive taxa. The cox1 locus has successfully been used to identify H. europaea in various other parts of the world where concerns have been raised on its presence, including Ukraine (Morhun et al. 2021), Spain (Reyes-Prieto et al. 2014) and New Zealand (Siddall and Budinoff 2005).

Molecular barcoding approaches have further characterised evolutionary relationships in the genus Helobdella, which has provided much needed taxonomic clarity. Siddall and Budinoff (2005) validated the current species name of H. europaea, after reporting this leech was originally named H. striata based on specimens collected in Germany and subsequently independently described as Helobdella papillornata from Australia – see Govedich and Davies (1998). Given that H. europaea is endemic to South America and affinities are apparent with other South American Helobdella spp., and also that unfortunately the appropriate name of Helobdella (triserialis) lineata is already assigned to a separate North American species, this name is valid (Siddall and Budinoff 2005).

Several studies have suggested that H. europaea is part of a species complex, which includes H. triserialis and H. socimulcensis (Siddall and Borda 2003; Siddall and Budinoff 2005; Morhun et al. 2021). Reconstructions generated during the current study show support for this H. triserialis complex, given the paraphyly of all H. europaea sequences included with H. triserialis and H. socimulcensis sequences for comparison.

The mitochondrial genome of H. europaea

The genome sequence determined for H. europaea is the first mitochondrial genome reported for the species and genus. Our MDA protocol and primers have also been used to achieve similar levels of mitochondrial (mt) genome enrichment with Papilio (butterfly) DNA, and we anticipate the protocol is applicable for mt genome enrichment in other invertebrates. Assembly and length determination of the control region in H. europaea was challenging, as this region comprises repeat sequences that were a feature of both MDA amplified and unamplified DNA. Dinucleotide repeats also occurred that are also a feature of the control region in Placobdella lamothei (NC_030269). Variation in the length of control region repeats has previously been observed in other animal species (Munwes et al. 2011), and it is possible that variation in the length of the control region also exists within individuals and populations of H. europaea.

The selective multiple displacement amplification and rapid barcoding Oxford Nanopore sequencing protocols used to sequence the mitochondrial genome of H. europaea provide a cost-effective means for sequencing the mitochondrial genomes of invertebrates. These protocols require a heating block and a microcentrifuge, but are PCR-free (i.e. do not require a PCR machine). Thus the protocol has potential applications for biodiversity and biosecurity assessments in a low infrastructure setting.

Implications for watershed management

The presence of H. europaea in the Ba River catchment is concerning given the potential negative impacts (Kutschera 1989; Kutschera 2004; Pfeiffer et al. 2004) this organism might have on the native invertebrate fauna of Fiji. Several area-endemic taxa are present in the Ba River catchment, with aquatic insects, molluscs and worms (Kutschera 2004; Málnás et al. 2016) potentially being vulnerable. A total of 73 freshwater macroinvertebrate taxa have been recorded from this catchment area, of which 14% are confirmed as endemic to Fiji with a further 47% being possibly endemic pending confirmation (Rashni 2020). These endemic taxa include five caddisflies (Abacaria fijiana, Abacaria ruficeps, Anisocentropus fijianus, Goera fijiana and Oxyethira fijiensis), a damselfly, Nesobasis spp., a shrimp (Caridina fijiana), a micro-water strider Fijivelia sp., a water cricket (Hydropedecticus vitiensis) and spring snails Fluviopupa spp. The presence of H. europaea may pose the greatest potential risk to area-endemic spring snails belonging to the genus Fluviopupa (Gastropoda: Tateidae), which are of high conservation significance. Currently a total of 28 Fluviopupa species are recorded from Fiji, all of which are area- and national- endemics (Zielske and Haase 2014), and are included in the Fiji Endangered and Protected Species (EPS) Act 2017.

The likely mode of introduction of H. europaea to Fiji remains unknown at the present time, and a thorough bioassessment of this species across the connected riverine network of both the Ba and Nadroga provinces is strongly recommended to understand the extent of spread and potential impacts on resident aquatic biodiversity. In this regard for Fiji and elsewhere this species may occur, the newly determined mitochondrial genome sequence for H. europaea provides a genome reference for higher resolution comparative analyses and source tracking. Such data will further assist the Fiji Invasive Species Taskforce (FIST) in their decision making for evaluating the status of H. europaea in Fiji.

Although H. europaea has been recognised as an invasive species globally (Kutschera 2004; Málnás et al. 2016; Morhun et al. 2021) and is alien to South Pacific Island countries, at this stage it can only be given the status of an introduced species in Fiji according to the corroborating definitions of Invasive Alien Species (IAS) via the Fiji National Biodiversity Strategies and Action Plans (NBSAPs), Convention on Biological Diversity (CBD) and the IUCN (International Union for Conservation of Nature) (Government of Fiji 2014; IUCN 2018; Department of Environment and Government of Fiji 2020). According to the Fiji NBSAP and CBD definitions, IAS are organisms found outside of their native geographical ranges that have spread and become invasive in their new habitats, and cause harm to biodiversity and other things that humans value (Government of Fiji 2014; Department of Environment and Government of Fiji 2020).

Results presented here offer a foundation on which further research assessing the extent of spread, ecology and ethology of H. europaea in Fiji and elsewhere may be based. In Fiji, outcomes of these future studies will further support the Biosecurity Authority of Fiji, the Ba and Nadroga-Navosa Provincial Administration and the Ba and Nadroga Yaubula Management Support Team to implement containment response protocols and formulate eradication strategies when/if required. Regionally in the Pacific and also globally, through the availability of the complete mitochondrial genome sequence for this species via this study, diagnosing its spread to other regions through traceability investigations is now possible.