Abstract

Leeches (Annelida: Hirudinea) possess powerful salivary anticoagulants and, accordingly, are frequently employed in modern, authoritative medicine. Members of the almost exclusively marine family Piscicolidae account for 20% of leech species diversity, and they feed on host groups (e.g., sharks) not encountered by their freshwater and terrestrial counterparts. Moreover, some species of Ozobranchidae feed on endangered marine turtles and have been implicated as potential vectors for the tumor-associated turtle herpesvirus. In spite of their ecological importance and unique host associations, there is a distinct paucity of data regarding the salivary transcriptomes of either of these families. Using next-generation sequencing, we profiled transcribed, putative anticoagulants and other salivary bioactive compounds that have previously been linked to blood feeding from 7 piscicolid species (3 elasmobranch feeders; 4 non-cartilaginous fish feeders) and 1 ozobranchid species (2 samples). In total, 149 putative anticoagulants and bioactive loci were discovered in varying constellations throughout the different samples. The putative anticoagulants showed a broad spectrum of described antagonistic pathways, such as inhibition of factor Xa and platelet aggregation, which likely have similar bioactive roles in marine fish and turtles. A transcript with homology to ohanin, originally isolated from king cobras, was found in Cystobranchus vividus but is otherwise unknown from leeches. Estimation of selection pressures for the putative anticoagulants recovered evidence for both positive and purifying selection along several isolated branches in the gene trees, and positive selection was also estimated for a few select codons in a variety of marine species. Similarly, phylogenetic analyses of the amino acid sequences for several anticoagulants indicated divergent evolution.

Salivary glands are integral organs for blood-feeding (hematophagous) animals, in that they secrete hemotoxic proteins, such as anticoagulants, into feeding sites to counter normal hemostatic processes (Munro et al., 1992; Fry et al., 2009; Kvist et al., 2013; Lemke et al., 2013). Blood-feeding leeches must counteract host hemostatic processes during their extended feeding period (approximately 25 min in Hirudo medicinalis, where they often consume over 10 times their weight in blood) and multi-month digestion period (Lent et al., 1988). Due to the potency of these salivary gland secretions in preventing clotting of blood (Greinacher and Warkentin, 2008), leeches are frequently used in modern medicine as post-surgical instruments following reconstructive and plastic surgery (Whitaker et al., 2004) and more generally as relievers of venous congestion (Conforti et al., 2002; Rajic and Deleyiannis, 2016; Shammas et al., 2016). In addition to their enduring value in medicine, anticoagulants more broadly have been critically important for pest control (e.g., rodenticides), the understanding of disease vectors, and deciphering complex evolutionary relationships (Mans et al., 2002a; Rokyta et al., 2011; Siddall et al., 2011; Song et al., 2011; Low et al., 2013).

Salivary transcriptome profiles (sialomes) increasingly have proven to be a rich source of comparative information from unrelated blood feeders that rely on similar vertebrate hosts (see Fry et al., 2009). Myriad marine species of Piscicolidae and Ozobranchidae (Annelida: Hirudinea) are clear candidates for salivary transcriptomic research, because they may possess unique adaptations for feeding in saltwater settings and on hosts that physiologically and anatomically differ from their freshwater and terrestrial counterparts. Specifically, piscicolids typically feed on fishes (cartilaginous and non-cartilaginous), whereas ozobranchids feed on turtles (Williams and Burreson, 2006; Utevsky et al., 2007; McGowin et al., 2011). Whereas Kvist et al. (2016) surveyed the anticoagulants found in the salivary transcriptome of the marine piscicolid Pontobdella macrothela (to date, the only survey of a piscicolid), Ozobranchids, which appear to be sister to piscicolids (Apakupakul et al., 1999), have not been surveyed for these loci. Piscicolids and ozobranchids also may differ in the number of anticoagulants they transcribe, given that they show differences in their host preference and that preliminary work has suggested that leeches lacking host specificity transcribe a greater diversity of these loci (Kvist et al., 2014, 2016).

Anticoagulants derived from piscicolid and ozobranchid taxa also merit further study due to these families' impacts in aquaria and on endangered host species, respectively. Inadvertent introduction of piscicolids into closed aquarium settings results in external ulcerations and, in severe infestations, host death (Marancik et al., 2012). Ozobranchus species cause turtle leech erosion disease (Bunkley-Williams et al., 2008), and are speculated to be vectors of fibropapilloma-associated herpesvirus in some sea turtles (Greenblatt et al., 2004). Elucidation of bioactive salivary loci from piscicolids and ozobranchids may serve as an initial framework in which to study leech-induced disease in marine fish and turtles.

The blood of elasmobranchs (sharks, skates, and rays) and of non-cartilaginous fishes represents highly divergent food sources—elasmobranch blood often contains 2 to 3 orders of magnitude more urea than that of non-cartilaginous fishes (Wells et al., 1986; Hammerschlag, 2006). Nonetheless, these strongly denaturing conditions (Zou et al., 1998; Bennion and Daggett, 2003) seem to present no impediment to those piscicolid species that effectively parasitize elasmobranchs.

Here, we sought to characterize bioactive loci in the salivary gland of marine leeches. We examined the diversity of anticoagulants, how that diversity related to the variety of hosts infested by the leeches, and whether or not anticoagulants from elasmobranch-feeding piscicolids were under specific selection pressures due to the varying host preference and the denaturing environment of shark blood.

MATERIALS AND METHODS

Specimen collection and transcriptome sequencing

Piscicolids and ozobranchids new to this study were collected in the wild and from aquaria (Table I). Specimens were preserved in RNAlater (Thermo Fisher Scientific, Waltham, Massachusetts) and then identified to species morphologically with specific identities later corroborated by transcribed mitochondrial cytochrome c oxidase subunit I sequences. Salivary glands were dissected aseptically when leeches were of sufficient size. For specimens with <2 mm anterior width, the anterior portion of the leech, where salivary tissue is contained, was used. Where possible, tissues from 3 specimens were pooled for each species in order to achieve a more representative transcriptomic profile. Specimens were generally considered to be fed, as most were found on hosts or recently had been on hosts. Thereafter, RNA was extracted from the dissected tissue, quantified, and sequenced using transcriptomic approaches.

Table I

Marine and brackish Piscicolidae and Ozobranchidae specimen and sequencing information.

Marine and brackish Piscicolidae and Ozobranchidae specimen and sequencing information.
Marine and brackish Piscicolidae and Ozobranchidae specimen and sequencing information.

For all specimens except Branchellion torpedinis (detailed below), RNA was extracted using the protocol of Kvist et al. (2014); tissue was placed in TRIzol and macerated by pestle or by ZR BashingBeads (Zymo Research, Irvine, California) before the addition of, and centrifugation with, chloroform. The resulting upper aqueous RNA layer was then removed and purified using an RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quality control and quantification of total RNA were conducted using a Qubit Fluorometer (Thermo Fisher Scientific) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California). Sample libraries were prepared using either TruSeq mRNA v2 or TruSeq stranded mRNA kits (Illumina, San Diego, California). Paired-end (2 × 125 base pairs) sequencing used Illumina HiSeq 2500 sequencers at the New York Genome Center.

For B. torpedinis, the salivary gland complex was frozen in liquid nitrogen and ground to a fine powder, and 5 mg portions of this powder were used for mRNA isolation using the mRNA Direct Micro kit (Thermo Fisher Scientific). Thereafter, cDNA was synthesized using SMART (Clontech, Mountain View, California) for library construction, with reverse transcription using Superscript III (Invitrogen, Waltham, Massachusetts) and PCR amplification using Platinum Taq DNA Polymerase High Fidelity (Invitrogen). PCR cycles consisted of denaturation at 94 C for 15 sec, annealing/extension at 68 C for 10 min, and, on the last cycle, a final extension for 12 min. After an initial round of 20 cycles, the PCR amplification was monitored by agarose gel electrophoresis of a 5-μl aliquot every 3–4 cycles, to prevent over-amplification and loss of longer cDNA transcripts. The amplified cDNA was size-selected on a Chroma-Spin 400 column (Clontech), and fractions containing cDNA > 250 base pairs were pooled into <1 kb and >1 kb fractions, to facilitate subsequent re-amplification and obtain cDNA for pyrosequencing. The small and large fractions were combined in equal proportions prior to sequencing. The cDNA sample was normalized with the Trimmer protocol (Evrogen, Moscow, Russia) prior to shearing into ∼650 base pair fragments and preparation of the emulsion phase library using Roche reagents. The library was sequenced by 454 pyrosequencing (Roche, Basel, Switzerland) using the ½ plate scale. Normalization, preparation of the emulsion phase PCR library, and 454 pyrosequencing were carried out at the University of Georgia Genomics Facility (http://dna.uga.edu/).

Transcriptome assembly and BLAST queries

Quality of transcriptomic sequences, both for newly acquired data as well as for those published previously (Kvist et al., 2016), was reviewed in FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Illumina sequences were trimmed using Trimmomatic (Bolger et al., 2014) with suggested settings for paired-end data. These sequences were then assembled using Trinity (Grabherr et al., 2011) for Illumina data and Newbler V2.30 for 454 pyrosequencing. Assembled contigs were then converted to open reading frames (ORFs) of 30 amino acids or longer using TransDecoder (Haas et al., 2013).

Identities of putative anticoagulants and bioactive salivary components were determined from ORFs using BLASTp against a local database consisting of known anticoagulants derived from leeches and other principally blood-feeding organisms with a minimum matching e-value of e−5. The local database was updated from prior studies (Kvist et al., 2014, 2016). The best-scoring matches were then reciprocally BLASTed against the UniProtKB/Swiss-Prot database using BLASTp in order to identify potential better matches against non-anticoagulant loci. In cases where better matches were found, these transcriptome sequences were disregarded. Better matches were considered as those with highest bit scores, given that e-values (among other things) are dependent on database size, which varied substantially between local searches and those against UniProtKB/Swiss-Prot. The same procedure was conducted both for ureases (using a local database created from metazoans ureases in UniProtKB/Swiss-Prot) and for post-translational modification enzymes (using a local database from loci found in Siddall et al., 2016).

BLAST matches that survived reciprocal searches were then both assessed for putative signal peptides with SignalP ver. 4.1 (Petersen et al., 2011) and aligned by codon positions using TranslatorX (Abascal et al., 2010) via MUSCLE (Edgar, 2004a, 2004b), applying the default strategy. Alignments were then used for all phylogenetic and evolutionary analyses.

Evolutionary and selection analyses

The ratio of non-synonymous to synonymous mutations was used to look for natural selection across loci of interest. When this ratio has more non-synonymous substitutions, then positive (Darwinian) selection is suspected. In contrast, when this ratio has more synonymous ones, then negative (purifying) selection is presumed to be the acting force. These analyses only focused on putative bioactive salivary loci that were found with signal peptides, as these have higher probability of being functional within the leeches. For the 2 transcriptomes of Ozobranchus margoi, in order to avoid pseudo-replication, if both had a given locus, only the closer match was used for analyses. Using default settings in Datamonkey (Pond and Frost, 2005; Murrell et al., 2012), aligned sequences were fitted with a model, and tests of positive selection along branches (BREL: branch-site random effects likelihood) and individual codons (MEME: mixed effects model of evolution) were performed. These comparisons were restricted to alignments that contained amino acids from at least 5 different species, with a minimum of 1 locus from species feeding on elasmobranchs or non-cartilaginous fishes. This threshold was applied to increase the robustness of the analyses.

To further compare anticoagulant evolution between these groups and as a second tier of orthology determination, phylogenetic reconstructions of each anticoagulant were then produced using parsimony, employing TNT (Goloboff et al., 2008) with 1,000 search replicates (each with 5 rounds of ratcheting, drifting, and fusing). Given the lack of objectively appropriate outgroup sequences for the anticoagulants, the unrooted strict consensus trees are presented. Support was calculated using 1,000 bootstrap pseudo-replicates.

Due to the potential for disulfide bonds to increase urea resistance, the numbers of cysteines also were compared. Specifically, the percentage of cysteines for each bioactive salivary locus was pooled by prey source (elasmobranchs or non-cartilaginous fishes). Differences between these paired groups were then tested using Wilcoxon signed-rank tests in R (R Core Team, 2014). For the cysteine comparison, each locus had to have a minimum of 1 matching sequence from an elasmobranch feeder and 1 from a non-cartilaginous fish feeder.

RESULTS

Transcriptomic assembly and BLAST queries for bioactive salivary

Raw sequencing reads were deposited in NCBI's Short Read Archive (SRA) as study SRP133366. Assembled transcriptomes were composed of between 54,830 and 220,149 contigs for Illumina data; the 454 assembly for B. torpedinis had 4,734 contigs (Table I). Mean contig length for each transcriptome ranged from 408 to 1,507 base pairs; median contig length ranged from 308 to 668 base pairs (Table I). In total, 345 assembled transcripts from all species studied matched 149 reference bioactive salivary components (Table II); all of these survived the reciprocal BLAST search against UniProtKB/Swiss-Prot.

Table II

Anticoagulant and other bioactive salivary loci found in marine leech transcriptomes. The 149 loci identified with BLAST searches are summarized by major groups. The number of loci found per group is listed for each species, with the number having signal peptides listed in parentheses. All loci matched had an e-value of e−5 or better.

Anticoagulant and other bioactive salivary loci found in marine leech transcriptomes. The 149 loci identified with BLAST searches are summarized by major groups. The number of loci found per group is listed for each species, with the number having signal peptides listed in parentheses. All loci matched had an e-value of e−5 or better.
Anticoagulant and other bioactive salivary loci found in marine leech transcriptomes. The 149 loci identified with BLAST searches are summarized by major groups. The number of loci found per group is listed for each species, with the number having signal peptides listed in parentheses. All loci matched had an e-value of e−5 or better.

The number of ORFs matching putative anticoagulants varied from 15 to 71 in the different transcriptomes (mean = 43). Of these, transcripts with verifiable signal peptides numbered from 5 to 23 (mean = 15). Signal peptides were recovered for many bioactive loci, including antistasin, a disintegrin and metalloprotease with thrombospondin motifs (ADAM[TS]), bdellin, C-type lectin, cystatin, destabilase I, therostasin, manillase, snaclec, leech-derived tryptase inhibitor (LDTI), infestin, leech antiplatelet protein (LAPP), eglin C, cystatin, ghilanten, piguamerin, reprolysin, rhodniin, and several others (see Table II). Additional bioactive loci of interest that were recovered but lacked signal peptides include acocostatin, apyrases, and saratin.

Many of these transcripts matched local database sequences from leeches, but a large number matched transcripts from other, mainly blood-feeding, invertebrates. The bulk of these loci were isolated from blood-feeding arthropods, such as lice (Pediculus), assassin bugs (Triatoma), mosquitoes (Aedes, Anopheles, and Culex), ticks (Argas, Amblyomma, Haemaphysalis, Ixodes, Ornithodoros, Rhipicephalus), biting flies (Lutzomyia and Tabanus), and copepods (Lepeophtheirus). These loci included ADAM (some with thrombospondin motifs; hereafter referred to as ADAM[TS]), c-type lectin, cystatin, apyrase, 5′ nucleotidase, rhodniin, infestin, and serine protease inhibitors. Additionally, BLAST matches were found for ADAM and c-type lectin from non-blood-feeding ants (Camponotus) and flesh flies (Sarcophaga), respectively. Matches to other invertebrate loci included nematodes (Brugia, Haemonchus, Loa, Wuchereria, Trichinella) and trematodes (Schistosoma).

Leech-derived transcripts also matched vertebrate sequences from venomous snakes, including Elapidae (Austrelaps, Demansia, and Ophiophagus), kraits (Bungarus), Colubridae (Philodryas), and vipers (Agkistrodon, Crotalus, Daboia, Deinagkistrodon, Echis, Gloydius, Ovophis, Sistrurus, and Trimeresurus). The matching loci were generally assigned to the ADAM(TS) family or other metalloproteases. Various c-type lectins (snaclecs) from vipers also were recovered.

Only a single putative urease (from Cystobranchus vividus) survived reciprocal BLAST searches. Post-translational modification loci, which likely play important roles relating to leech bioactive salivary compounds (Siddall et al., 2016), were found in all the reviewed species. These included carbonic anhydrase, glycosyltransferase, protein disulfide-isomerase, and oligosaccharyltransferase. The loci found varied between leech species but did not correspond to elasmobranch-feeding vs. non-cartilaginous-feeding species (or any other notable pattern).

Comparisons based on elasmobranch vs. non-cartilaginous fish prey

Five putative anticoagulants (therostasin, antistasin, destabilase I, ADAM[TS], and manillase) remained after filtering out loci that did not possess a signal peptide and for which less than 5 comparative sequences were recovered. Although some amino acid sites were found to be under positive selection for 4 of the 5 anticoagulants reviewed, purifying selection seems to be acting at a broader scale across the amino acid positions. On average, 2 sites were under positive selection for each bioactive salivary locus. Several additional codons showed P-values only slightly higher than 0.05 for several loci (e.g., destabilase has 3 additional codons with P-values <0.10). Positively selected codons were also found along piscicolid branches (Fig. 1), including those leading to species that feed both on elasmobranch (Oxytonostoma typica) and non-cartilaginous fishes (Pterobdella cf. abditovesiculata and the branch uniting the Heptacyclus species).

Figure 1

Phenograms for bioactive salivary compounds found to have codons under positive selection. For each compound, the phenogram topology is consistent, while selection pressures are displayed for each codon found to have some positive selection occurring. Branches are colored according to selection pressure for a given codon. Light gray or blue indicates purifying selection, while dark gray or red indicates positive selection; black represents neutral selection pressures. Thicker branches indicate statistically significant positive selection, i.e., more non-synonymous changes than expected by chance alone. Branch lengths represent the predicted nucleotide substitutions per site. Color version available online.

Figure 1

Phenograms for bioactive salivary compounds found to have codons under positive selection. For each compound, the phenogram topology is consistent, while selection pressures are displayed for each codon found to have some positive selection occurring. Branches are colored according to selection pressure for a given codon. Light gray or blue indicates purifying selection, while dark gray or red indicates positive selection; black represents neutral selection pressures. Thicker branches indicate statistically significant positive selection, i.e., more non-synonymous changes than expected by chance alone. Branch lengths represent the predicted nucleotide substitutions per site. Color version available online.

Three of the 5 putative anticoagulants that were reviewed appeared to have overall positive selection on at least 1 given branch (Fig. 2). As with analyses of codon selection, evidence of positive selection was found on branches leading to both Ozobranchidae (for therostasin and ADAM[TS]) and several representatives of Piscicolidae (for 3 anticoagulants).

Figure 2

Phenograms for bioactive salivary compounds focusing on overall selection by branch. In the color version (available online), blue indicates purifying selection, while red indicates positive selection; black represents neutral selection pressures. Thicker branches indicate statistically significant positive selection. Branch lengths represent the predicted nucleotide substitutions per site.

Figure 2

Phenograms for bioactive salivary compounds focusing on overall selection by branch. In the color version (available online), blue indicates purifying selection, while red indicates positive selection; black represents neutral selection pressures. Thicker branches indicate statistically significant positive selection. Branch lengths represent the predicted nucleotide substitutions per site.

Topologies from phylogenetic reconstructions of the bioactive salivary loci were best resolved for manillase, ADAM(TS), and therostasin (Fig. 3). Heptacyclus species had particularly long branches in the reconstruction of the ADAM(TS) and manillase trees; P. macrothela had a modestly long branch in the therostasin tree. Anticoagulants from Heptacyclus species were always adjacent in the phylogenies.

Figure 3

Unrooted phylogenetic reconstructions of bioactive salivary compounds. Support values are summaries of 1,000 bootstrap pseudo-replicates. Branch lengths represent the number of nucleotide changes.

Figure 3

Unrooted phylogenetic reconstructions of bioactive salivary compounds. Support values are summaries of 1,000 bootstrap pseudo-replicates. Branch lengths represent the number of nucleotide changes.

Leeches feeding on elasmobranchs had almost identical (marginally higher on average) percentages of cysteines in their putative bioactive salivary components when compared to those of leeches feeding on non-cartilaginous fishes. The marginal difference was not statistically significant when including (P = 0.19) or excluding (P = 0.24) the turtle-feeding Ozobranchus with the non-cartilaginous fish feeders. Accordingly, cysteines in these loci may not provide protection against the denaturing effects of urea, which is found in high quantities in elasmobranch blood.

DISCUSSION

The present study represents the largest comparison of leech anticoagulants to date. The marine and brackish leeches examined transcribed a pharmacopeia of 149 putative anticoagulants and other bioactive salivary loci that have previously been linked to anticoagulation. These loci disrupt a variety of steps in the blood coagulation cascade by inhibiting the formation of various loci, such as factor Xa, thrombin, and serine proteases. The loci varied substantially regarding amino acid composition, molecular selection pressure, and phylogenetic topology.

It is clear from these findings that, when screening for anticoagulants, reviewing a broad suite of closely related blood-feeding organisms is critical for uncovering the full array of loci transcribed. Each species had at least 1 transcript that matched a locus not found in any of the other species reviewed. A third of the species possessed 10 or more putative anticoagulation loci that were not recovered from any other species reviewed. Such substantial differences were even found for congeneric species of Heptacyclus.

Numerous loci found in this study were originally isolated from leeches (e.g., antistasin, LDTI, and LAPP). However, many of the transcripts best matched putative orthologs isolated from a variety of non-leech animal taxa, including invertebrates (e.g., blood-feeding arthropods) and vertebrates (e.g., venomous snakes). The majority of these BLAST matches were for ADAM(TS), but a broad diversity of loci were found. The ADAM(TS) family includes secretory loci that are divided into 20 families and impart diverse biological functions through functional disintegrin, metalloprotease, and thrombospondin domains (Tang, 2001). These loci have numerous described anticoagulant roles that likely aid in blood feeding, including inhibition of fibrinogen binding and platelet aggregation, cleavage of the Von Willebrand factor, and endothelial cell disruption (Siigur et al., 2001; Tang, 2001; Mans et al., 2002b; Francischetti, 2010). Additionally, a transcript with putative orthology to ohanin, described from the king cobra and known to slow locomotion in mice (Pung et al., 2005), was recovered in the salivary transcriptomes of C. vividus and also possessed a signal peptide. Several other loci that have never before been recovered in leeches (mostly variants of ADAM[TS] and zinc metalloproteases) were also found, but they often lacked signal peptides, suggesting that they may not be secreted by the leeches or that assemblies of these loci may have been incomplete. For example, a putative ortholog of acocostatin, a disintegrin-like peptide from copperhead snakes that inhibits platelet aggregation, amongst other functions (Teklemariam et al., 2011), was recovered. Shared characteristics between mammalian platelets and thrombocytes in other vertebrates (Pica et al., 1990; Soslau et al., 2005) suggest similar anti-hemostatic roles during leech feeding in fish and turtles.

These results and similar findings from other studies focusing on different animals (Fry et al., 2009; Macrander et al., 2015; Rosenfeld et al., 2016) confirm the importance of reviewing venom, poison, and anticoagulant loci from a broad diversity of organisms in comparative transcriptomic and genomic studies. The possible connection between venoms and leech anticoagulants warrants further studies into the similarities of bioactive pathways of the loci, as well as their phylogeny, in order to infer orthology and evolutionary history. Although these loci have yet to be functionally tested in leeches, their ubiquity in salivary (and venom) glands and the fact that many have signal peptides suggest that these loci retain similar functionality.

Phylogenetic and anticoagulant studies alike suggest that leeches developed from a single sanguivorous ancestor (Apakupakul et al., 1999; Siddall et al., 2016). These transcriptomic data, along with the work by Kvist et al. (2016), suggest that destabilase in particular may have arisen earlier in leech evolution than previously believed; putative orthologs for this locus were recovered from both O. margoi and a variety of piscicolids. Destabilase has demonstrated activity in breaking isopeptide linking within fibrin clots and inhibition of platelet aggregation, which may have roles in disintegrating blood clots at leech feeding sites (Baskova and Nikonov, 1991). While formerly considered to be phylogenetically restricted to Arhynchobdellida (Min et al., 2010), destabilase most likely originated in the common ancestor of arhynchobdellids, piscicolid, and ozobranchids—these groups form 1 of the 2 principal clades in leeches (Apakupakul et al., 1999). In contrast, some factors that are present in both glossiphoniid and arhynchobdellid leeches, such as hirudin (Müller et al., 2016, 2017; Siddall et al., 2016), were not found in any of the marine leech transcriptomes.

In addition to the variety of loci recovered from marine leeches, several of these bioactive salivary components displayed heterogeneous degrees of selection pressures, both positive and purifying, and phylogenetic topologies. These findings highlight suitable protein sequences for further functional investigation and hint at loci that putatively represent unique ways in which leeches overcome host coagulation. For example, Heptacyclus species had an unusually long branch for the phylogenetic reconstruction of manillase and ADAM(TS), and the branch to O. typica was under strong positive selection for an ADAM(TS) locus. However, there was no clear ecological or taxonomic pattern in selection tests and phylogenetic reconstructions, as is often the case for leech anticoagulants (see Kvist et al. [2013, 2014, 2016] regarding selection on leech anticoagulants and Arcà et al. [2014] regarding mosquito anticoagulants).

Although future studies of piscicolid salivary transcriptomics should also include freshwater species, the principally marine focus (C. vividus and P. cf. abditovesiculata were from brackish water) of the present study provided a platform for reviewing anticoagulant diversity of each species with regard to its host specificity, which had not been fully possible in prior studies. No discernible pattern was found relating the number of anticoagulants to the number of known hosts. Elasmobranch-feeding piscicolids had somewhat fewer anticoagulant transcripts found, but B. torpedinis and P. macrothela species feed on a diversity of elasmobranchs (Sawyer et al., 1975). In contrast, prior work had hinted that leeches feeding on a broad array of hosts might have greater anticoagulant diversity. This is highlighted by the particularly large number of anticoagulants found in the terrestrial leech Haemadipsa interrupta (Kvist et al., 2014), which shows a wider host array (see Tessler et al., 2018), as compared to the more limited host diversity found for the elasmobranch-feeding P. macrothela (Kvist et al., 2016) and select turtle- and frog-feeding species of Placobdella (Siddall et al., 2016).

In terms of the denaturing aspect of elasmobranch blood, no notable differences were found between the anticoagulant (or urease) repertoires of elasmobranch-feeding and non-cartilaginous-fish–feeding leeches or the selection pressures on these loci. Differences between elasmobranch-feeding and non-cartilaginous-fish–feeding piscicolids may yet be found among the bacterial communities in their gut, nephridia, and bladders. Distinctive bacterial communities are found in nephridia and bladders—areas associated with nitrogenous wastes, such as urea—and are thought to detoxify and degrade these wastes (Dev, 1964; Sawyer, 1986; Kikuchi et al., 2009; Nelson and Graf, 2012). Highly efficient degrading and detoxifying bacterial species may therefore represent a path for adaptation to elasmobranch feeding in piscicolids.

The present study provides insight into the bioactive salivary components of 2 families of marine leeches that have otherwise mainly been the focus of systematics studies (Williams and Burreson, 2006; Utevsky et al., 2007) and studies on their potential to transmit viruses and parasites to hosts (Siddall and Desser, 1992, 1993; Greenblatt et al., 2004). The large numbers of anticoagulants found here help to fill gaps in knowledge regarding the diversity of anticoagulants across leeches. Elucidation of these loci provides clues into the evolution of blood feeding in leeches and the potential role of these factors in the host–parasite relationship. This study also identifies outlier loci that possibly represent novel function and are therefore worthy of further exploration. Furthermore, our results detail the importance of intense sampling of both close and distant relatives for salivary transcriptomic studies.

ACKNOWLEDGMENTS

We thank Eugene Burreson, Troy Tuckey, and Emily Loose (Virginia Institute of Marine Science); Julio Lorda (Universidad Autónoma de Baja California); Taylor White (Sitka Sound Science Center); Jeff Schwenter and Mike Arendt (South Carolina Department of Natural Resources); Shane Boylan (South Carolina Aquarium); Huntsman Marine Science Center; and the veterinary and husbandry staff at the Georgia Aquarium for their assistance with this project and their aid in our procurement and collection of the specimens utilized here. Thanks go to Rob DeSalle, Danielle de Carle, and Spencer Galen for helping refine early versions of this manuscript. Furthermore, we thank AMNH for a Lerner-Gray Marine Research Grant awarded to M.T. for helping to fund fieldwork in Canada.

LITERATURE CITED

LITERATURE CITED
Abascal,
F.,
R.
Zardoya,
and
M. J.
Telford.
2010
.
TranslatorX: Multiple alignment of nucleotide sequences guided by amino acid translations
.
Nucleic Acids Research
38
:
W7
13
.
Apakupakul,
K.,
M. E.
Siddall,
and
E. M.
Burreson.
1999
.
Higher level relationships of leeches (Annelida: Clitellata: Euhirudinea) based on morphology and gene sequences
.
Molecular Phylogenetics and Evolution
12
:
350
359
.
Arcà,
B.,
C. J.
Struchiner,
V. M.
Pham,
G.
Sferra,
F.
Lombardo,
M.
Pombi,
and
J. M.
Ribeiro.
2014
.
Positive selection drives accelerated evolution of mosquito salivary genes associated with blood-feeding
.
Insect Molecular Biology
23
:
122
131
.
Baskova,
I. P.,
and
G. I.
Nikonov.
1991
.
Destabilase, the novel [varepsilon]-([gamma]-Glu)-Lys isopeptidase with thrombolytic activity
.
Blood Coagulation and Fibrinolysis
2
:
167
172
.
Bennion,
B. J.,
and
V.
Daggett.
2003
.
The molecular basis for the chemical denaturation of proteins by urea
.
Proceedings of the National Academy of Sciences of the United States of America
100
:
5142
5147
.
Bolger,
A. M.,
M.
Lohse,
and
B.
Usadel.
2014
.
Trimmomatic: A flexible trimmer for Illumina sequence data
.
Bioinformatics
30
:
2114
2120
.
Bunkley-Williams,
L.,
E. H.
Williams,
J. A.
Horrocks,
H. C.
Horta,
A. A.
Mignucci-Giannoni,
and
A. C.
Poponi.
2008
.
New leeches and diseases for the hawksbill sea turtle and the West Indies
.
Comparative Parasitology
75
:
263
270
.
Conforti,
M. L.,
N. P.
Connor,
D. M.
Heisey,
R.
Vanderby,
D.
Kunz,
and
G. K.
Hartig.
2002
.
Development of a mechanical device to replace medicinal leech (Hirudo medicinalis) for treatment of venous congestion
.
Journal of Rehabilitation Research and Development
39
:
497
504
.
Dev,
B.
1964
.
Excretion and osmoregulation in the leech, Hirudinaria granulosa (Savigny)
.
Nature
202
:
414
.
Edgar,
R. C.
2004
a
.
MUSCLE: A multiple sequence alignment method with reduced time and space complexity
.
BMC Bioinformatics
5
:
113
. doi:.
Edgar,
R. C.
2004
b
.
MUSCLE: Multiple sequence alignment with high accuracy and high throughput
.
Nucleic Acids Research
32
:
1792
1797
.
Francischetti,
I. M. B.
2010
.
Platelet aggregation inhibitors from hematophagous animals
.
Toxicon
56
:
1130
1144
.
Fry,
B. G.,
K.
Roelants,
D. E.
Champagne,
H.
Scheib,
J. D. A.
Tyndall,
G. F.
King,
T. J.
Nevalainen,
J. A.
Norman,
R. J.
Lewis,
R. S.
Norton,
et al.
2009
.
The toxicogenomic multiverse: Convergent recruitment of proteins into animal venoms
.
Annual Review of Genomics and Human Genetics
10
:
483
511
.
Goloboff,
P. A.,
J. S.
Farris,
and
K. C.
Nixon.
2008
.
TNT, a free program for phylogenetic analysis
.
Cladistics
24
:
774
786
.
Grabherr,
M. G.,
B. J.
Haas,
M.
Yassour,
J. Z.
Levin,
D. A.
Thompson,
I.
Amit,
X.
Adiconis,
L.
Fan,
R.
Raychowdhury,
Q.
Zeng,
et al.
2011
.
Full-length transcriptome assembly from RNA-Seq data without a reference genome
.
Nature Biotechnology
29
:
644
652
.
Greenblatt,
R. J.,
T. M.
Work,
G. H.
Balazs,
C. A.
Sutton,
R. N.
Casey,
and
J. W.
Casey.
2004
.
The Ozobranchus leech is a candidate mechanical vector for the fibropapilloma-associated turtle herpesvirus found latently infecting skin tumors on Hawaiian green turtles (Chelonia mydas)
.
Virology
321
:
101
110
.
Greinacher,
A.,
and
T. E.
Warkentin.
2008
.
The direct thrombin inhibitor hirudin
.
Thrombosis and Haemostasis
99
:
819
829
.
Haas,
B. J.,
A.
Papanicolaou,
M.
Yassour,
M.
Grabherr,
P. D.
Blood,
J.
Bowden,
M. B.
Couger,
D.
Eccles,
B.
Li,
M.
Lieber,
et al.
2013
.
De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis
.
Nature Protocols
8
:
1494
1512
.
Hammerschlag,
N.
2006
.
Osmoregulation in elasmobranchs: A review for fish biologists, behaviourists and ecologists
.
Marine and Freshwater Behaviour and Physiology
39
:
209
228
.
Kikuchi,
Y.,
L.
Bomar,
and
J.
Graf.
2009
.
Stratified bacterial community in the bladder of the medicinal leech, Hirudo verbana
.
Environmental Microbiology
11
:
2758
2770
.
Kvist,
S.,
M. R.
Brugler,
T. G.
Goh,
G.
Giribet,
and
M. E.
Siddall.
2014
.
Pyrosequencing the salivary transcriptome of Haemadipsa interrupta (Annelida: Clitellata: Haemadipsidae): Anticoagulant diversity and insight into the evolution of anticoagulation capabilities in leeches
.
Invertebrate Biology
133
:
74
98
.
Kvist,
S.,
G.-S.
Min,
and
M. E.
Siddall.
2013
.
Diversity and selective pressures of anticoagulants in three medicinal leeches (Hirudinida: Hirudinidae, Macrobdellidae)
.
Ecology and Evolution
3
:
918
933
.
Kvist,
S.,
A.
Oceguera-Figueroa,
M.
Tessler,
J.
Jiménez-Armenta,
R. M.
Freeman,
G.
Giribet,
and
M. E.
Siddall.
2016
.
When predator becomes prey: Investigating the salivary transcriptome of the shark-feeding leech Pontobdella macrothela (Hirudinea: Piscicolidae)
.
Zoological Journal of the Linnean Society
179
:
725
737
.
Lemke,
S.,
C.
Müller,
E.
Lipke,
G.
Uhl,
and
J.-P.
Hildebrandt.
2013
.
May salivary gland secretory proteins from hematophagous leeches (Hirudo verbana) reach pharmacologically relevant concentrations in the vertebrate host?
PLoS One
8
:
e73809
. doi:.
Lent,
C. M.,
K. H.
Fliegner,
E.
Freedman,
and
M. H.
Dickinson.
1988
.
Ingestive behaviour and physiology of the medicinal leech
.
Journal of Experimental Biology
137
:
513
527
.
Low,
D. H. W.,
K.
Sunagar,
E. A. B.
Undheim,
S. A.
Ali,
A. C.
Alagon,
T.
Ruder,
T.N.W.
Jackson,
S. P.
Gonzalez,
G. F.
King,
A.
Jones,
et al.
2013
.
Dracula's children: Molecular evolution of vampire bat venom
.
Journal of Proteomics
89
:
95
111
.
Macrander,
J.,
M. R.
Brugler,
and
M.
Daly.
2015
.
A RNA-seq approach to identify putative toxins from acrorhagi in aggressive and non-aggressive Anthopleura elegantissima polyps
.
BMC Genomics
16
:
221
. doi:.
Mans,
B. J.,
A. I.
Louw,
and
A. W. H.
Neitz.
2002
a
.
Evolution of hematophagy in ticks: Common origins for blood coagulation and platelet aggregation inhibitors from soft ticks of the genus Ornithodoros
.
Molecular Biology and Evolution
19
:
1695
1705
.
Mans,
B. J.,
A. I.
Louw,
and
A. W. H.
Neitz.
2002
b
.
Savignygrin, a platelet aggregation inhibitor from the soft tick Ornithodoros savignyi, presents the RGD integrin recognition motif on the Kunitz-BPTI fold
.
The Journal of Biological Chemistry
277
:
21371
21378
.
Marancik,
D. P.,
A. D.
Dove,
and
A. C.
Camus.
2012
.
Experimental infection of yellow stingrays Urobatis jamaicensis with the marine leech Branchellion torpedinis
.
Diseases of Aquatic Organisms
101
:
51
60
.
McGowin,
A. E.,
T. M.
Truong,
A. M.
Corbett,
D. A.
Bagley,
L. M.
Ehrhart,
M. J.
Bresette,
S. T.
Weege,
and
D.
Clark.
2011
.
Genetic barcoding of marine leeches (Ozobranchus spp.) from Florida sea turtles and their divergence in host specificity
.
Molecular Ecology Resources
11
:
271
278
.
Min,
G.-S.,
I. N.
Sarkar,
and
M. E.
Siddall.
2010
.
Salivary transcriptome of the North American medicinal leech, Macrobdella decora
.
Journal of Parasitology
96
:
1211
1221
.
Müller,
C.,
M.
Haase,
S.
Lemke,
and
J.-P.
Hildebrandt.
2017
.
Hirudins and hirudin-like factors in Hirudinidae: Implications for function and phylogenetic relationships
.
Parasitology Research
116
:
313
325
.
Müller,
C.,
K.
Mescke,
S.
Liebig,
H.
Mahfoud,
S.
Lemke,
and
J.-P.
Hildebrandt.
2016
.
More than just one: Multiplicity of hirudins and hirudin-like factors in the medicinal leech, Hirudo medicinalis
.
Molecular Genetics and Genomics
291
:
227
240
.
Munro,
R.,
M.
Siddall,
S. S.
Desser,
and
R. T.
Sawyer.
1992
.
The leech as a tool for studying comparative haematology
.
Comparative Haematology International
2
:
75
78
.
Murrell,
B.,
J. O.
Wertheim,
S.
Moola,
T.
Weighill,
K.
Scheffler,
and
S. L.
Kosakovsky Pond.
2012
.
Detecting individual sites subject to episodic diversifying selection
.
PLoS Genetics
8
:
e1002764
. doi:.
Nelson,
M. C.,
and
J.
Graf.
2012
.
Bacterial symbioses of the medicinal leech Hirudo verbana
.
Gut Microbes
3
:
322
331
.
Petersen,
T. N.,
S.
Brunak,
G.
von Heijne,
and
H.
Nielsen.
2011
.
SignalP 4.0: Discriminating signal peptides from transmembrane regions
.
Nature Methods
8
:
785
786
.
Pica,
A.,
A.
Lodato,
M. C.
Grimaldi,
and
F.
Della Corte.
1990
.
Morphology, origin and functions of the thrombocytes of Elasmobranchs
.
Archivio Italiano di Anatomia e di Embriologia
95
:
187
207
.
Pond,
S. L. K.,
and
S. D. W.
Frost.
2005
.
Datamonkey: Rapid detection of selective pressure on individual sites of codon alignments
.
Bioinformatics
21
:
2531
2533
.
Pung,
Y. F.,
P. T. H.
Wong,
P. P.
Kumar,
W. C.
Hodgson,
and
R. M.
Kini.
2005
.
Ohanin, a novel protein from king cobra venom, induces hypolocomotion and hyperalgesia in mice
.
Journal of Biological Chemistry
280
:
13137
13147
.
Rajic,
A. J.,
and
F. W.-B.
Deleyiannis.
2016
.
Determining the duration of leech therapy in the treatment of acute venous congestion in prefabricated free flaps
.
Plastic and Reconstructive Surgery
137
:
495e
496e
. doi:.
R Core Team
.
2014
.
R: A language and environment for statistical computing
.
R Foundation for Statistical Computing
,
Vienna, Austria
.
Rokyta,
D. R.,
K. P.
Wray,
A. R.
Lemmon,
E. M.
Lemmon,
and
S. B.
Caudle.
2011
.
A high-throughput venom-gland transcriptome for the Eastern Diamondback Rattlesnake (Crotalus adamanteus) and evidence for pervasive positive selection across toxin classes
.
Toxicon
57
:
657
671
.
Rosenfeld,
J. A.,
D.
Reeves,
M. R.
Brugler,
A.
Narechania,
S.
Simon,
R.
Durrett,
J.
Foox,
K.
Shianna,
M. C.
Schatz,
J.
Gandara,
et al.
2016
.
Genome assembly and geospatial phylogenomics of the bed bug Cimex lectularius
.
Nature Communications
7
:
10164
.
Sawyer,
R. T.
1986
.
Leech biology and behaviour: Feeding biology, ecology, and systematics
.
Oxford University Press
,
New York, New York
,
1065
p
.
Sawyer,
R. T.,
A. R.
Lawler,
and
R. M.
Overstreet.
1975
.
Marine leeches of the eastern United States and the Gulf of Mexico with a key to the species
.
Journal of Natural History
9
:
633
667
.
Shammas,
R. L.,
A.
Cornejo,
L. P.
Poveromo,
A. D.
Glener,
and
S. T.
Hollenbeck.
2016
.
Outcomes of leech therapy for venous congestion after reconstructive surgery
.
Journal of the American College of Surgeons
223
(
Suppl. 1
):
S100
. doi:.
Siddall,
M. E.,
M. R.
Brugler,
and
S.
Kvist.
2016
.
Comparative transcriptomic analyses of three species of Placobdella (Rhynchobdellida: Glossiphoniidae) confirms a single origin of blood feeding in leeches
.
Journal of Parasitology
102
:
143
150
.
Siddall,
M. E.,
and
S. S.
Desser.
1992
.
Ultrastructure of gametogenesis and sporogony of Haemogregarina (sensu lato) myoxocephali (Apicomplexa: Adeleina) in the marine leech Malmiana scorpii
.
The Journal of Protozoology
39
:
545
554
.
Siddall,
M. E.,
and
S. S.
Desser.
1993
.
Ultrastructure of merogonic development of Haemogregarina (sensu lato) myoxocephali (Apicomplexa: Adeleina) in the marine leech Malmiana scorpii and localization of infective stages in the salivary cells
.
European Journal of Protistology
29
:
191
201
.
Siddall,
M. E.,
G.-S.
Min,
F. M.
Fontanella,
A. J.
Phillips,
and
S. C.
Watson.
2011
.
Bacterial symbiont and salivary peptide evolution in the context of leech phylogeny
.
Parasitology
138
:
1815
1827
.
Siigur,
J.,
A.
Aaspõllu,
K.
Tõnismägi,
K.
Trummal,
M.
Samel,
H.
Vija,
J.
Subbi,
and
E.
Siigur.
2001
.
Proteases from Vipera lebetina venom affecting coagulation and fibrinolysis
.
Haemostasis
31
:
123
132
.
Song,
Y.,
S.
Endepols,
N.
Klemann,
D.
Richter,
F.-R.
Matuschka,
C.-H.
Shih,
M. W.
Nachman,
and
M. H.
Kohn.
2011
.
Adaptive introgression of anticoagulant rodent poison resistance by hybridization between old world mice
.
Current Biology
21
:
1296
1301
.
Soslau,
G.,
P. J.
Prest,
R.
Class,
R.
George,
F.
Paladino,
and
G.
Violetta.
2005
.
Comparison of sea turtle thrombocyte aggregation to human platelet aggregation in whole blood
.
Comparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology
142
:
353
360
.
Tang,
B. L.
2001
.
ADAMTS: A novel family of extracellular matrix proteases
.
International Journal of Biochemistry and Cell Biology
33
:
33
44
.
Teklemariam,
T.,
A. I.
Seoane,
C. J.
Ramos,
E. E.
Sanchez,
S. E.
Lucena,
J. C.
Perez,
S. A.
Mandal,
and
J. G.
Soto.
2011
.
Functional analysis of a recombinant PIII-SVMP, GST-acocostatin; an apoptotic inducer of HUVEC and HeLa, but not SK-Mel-28 cells
.
Toxicon
57
:
646
656
.
Tessler,
M.,
S. R.
Weiskopf,
L.
Berniker,
R.
Hersch,
K. P.
McCarthy,
D. W.
Yu, and
M. E.
Siddall.
2018
.
Bloodlines: Mammals, leeches, and conservation in southern Asia
.
Systematics and Biodiversity
,
p
.
1
9
. doi: .
Utevsky,
S. Y.,
A. Y.
Utevsky,
S.
Stefano,
and
T.
Peter.
2007
.
Molecular phylogeny of pontobdelline leeches and their place in the descent of fish leeches (Hirudinea, Piscicolidae)
.
Zoologica Scripta
36
:
271
280
.
Wells,
R. M.,
R. H.
McIntyre,
A. K.
Morgan,
and
P. S.
Davie.
1986
.
Physiological stress responses in big gamefish after capture: Observations on plasma chemistry and blood factors
.
Comparative Biochemistry and Physiology Part A: Physiology
84
:
565
571
.
Whitaker,
I. S.,
D.
Izadi,
D. W.
Oliver,
G.
Monteath,
and
P. E.
Butler.
2004
.
Hirudo medicinalis and the plastic surgeon
.
British Journal of Plastic Surgery
57
:
348
353
.
Williams,
J. I.,
and
E. M.
Burreson.
2006
.
Phylogeny of the fish leeches (Oligochaeta, Hirudinida, Piscicolidae) based on nuclear and mitochondrial genes and morphology
.
Zoologica Scripta
35
:
627
639
.
Zou,
Q.,
S. M.
Habermann-Rottinghaus,
and
K. P.
Murphy.
1998
.
Urea effects on protein stability: Hydrogen bonding and the hydrophobic effect
.
Proteins: Structure, Function, and Genetics
31
:
107
115
.