Introduction
Bloodfeeding requires several adaptations to overcome biotic barriers imposed by hosts and to enable efficient feeding. Such adaptations include morphological, behavioral and molecular features that allow the parasites to create an incision and inhibit hemostasis, inflammation and other responses from the host immune system (Ribeiro and Francischetti, Reference Ribeiro and Francischetti2003; Ribeiro and Arcà, Reference Ribeiro and Arcà2009; Min et al., Reference Min, Sarkar and Siddall2010). In this regard, leeches have attracted attention, especially due to their intense use in medicine in the 18th and 19th centuries and the potent anticoagulants isolated from salivary secretions of different species, especially from the European medicinal leech, Hirudo medicinalis Linnaeus, 1758 (Elliott and Kutschera, Reference Elliott and Kutschera2011).
Traditionally, leeches were subdivided into two orders based on the morphology of their feeding apparatus. Members of Rhynchobdellida are proboscis-bearing and insert their proboscis subcutaneously in order to feed on blood. By contrast, arhynchobdelid leeches are non-proboscis-bearing and the bloodfeeding taxa within this group use denticle-bearing jaws to create wounds from which they feed (Sawyer, Reference Sawyer1986). Studies based on molecular data have suggested that Rhyncobdellida is paraphyletic (e.g. Apakupakul et al., Reference Apakupakul, Siddall and Burreson1999) and, as a result, Tessler et al. (Reference Tessler, de Carle, Voiklis, Gresham, Neumann, Cios and Siddall2018a) erected five orders to accommodate known leech diversity (Oceanobdelliformes, Glossiphoniformes, Erpobdelliformes, Americobdelliformes and Hirudiniformes). Additionally, phylogenetic studies have contributed to clarifying the evolution of bloodfeeding in leeches. Siddall and Burreson (Reference Siddall and Burreson1995) were first to propose bloodfeeding as the behavior of the last common ancestor of all leeches, with subsequent independent transformations to macrophagy within Hirudinida. This notion is also supported by subsequent phylogenetic analyses (Apakupakul et al., Reference Apakupakul, Siddall and Burreson1999; Borda and Siddall, Reference Borda and Siddall2004; Tessler et al., Reference Tessler, de Carle, Voiklis, Gresham, Neumann, Cios and Siddall2018a).
Salivary transcriptomic data and the presence of putative anticoagulants in both non-proboscis and proboscis-bearing leeches have been used as evidence for an ancestral bloodfeeding state in Hirudinida (see Siddall et al., Reference Siddall, Brugler and Kvist2016). However, some of these studies have focused on BLAST determination for the majority of the putative anticoagulants found, instead of tree-based orthology determination, as in other studies on non-proboscis-bearing leeches (Kvist et al., Reference Kvist, Min and Siddall2013, Reference Kvist, Brugler, Goh, Giribet and Siddall2014, Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017; Müller et al., Reference Müller, Haase, Lemke and Hildebrandt2017). The lack of a proper phylogenetic hypothesis of the major leech anticoagulants questions previous conclusions that the presence of such molecules across the leech tree is evidence of bloodfeeding behavior in the last common ancestor. Therefore, an analysis including putative anticoagulants of a wide phylogenetic diversity of leeches is needed to determine which molecules are orthologous to proteins that were present in the last common ancestor of leeches.
Limnobdella mexicana Blanchard, 1893 is a member of Praobdellidae, a family of jawed leeches. Although little is known about the host diversity and bloodfeeding behavior of this species, members of this family are known to cause internal hirudiniasis – leech infestations in orifices such as the buccal mucosa, nasal cavity and vesicles – and most have been found to feed indiscriminately on mammals (including humans) and birds (Oka, Reference Oka1934; Phillips et al., Reference Phillips, Arauco-Brown, Oceguera-Figueroa, Gomez, Beltrán, Lai and Siddall2010; Siddall et al., Reference Siddall, Rood-goldman, Barriot and Barboutist2013); and recently a member of this family has been recorded as feeding on crustaceans (Nakano et al., Reference Nakano, Tomikawa, Sakono and Yoshikawa2017). Previous leech salivary transcriptomic analyses have reported a positive relationship between host-generality and diversity of anticoagulants (Kvist et al., Reference Kvist, Brugler, Goh, Giribet and Siddall2014). Thus, it is predicted that leeches with higher levels of host-generality, such as that recorded for members of the family Praobdellidae, possess a proportionately higher diversity of anticoagulants. Despite the medical and veterinary importance of praobdellid leeches, there are no studies that describe the anticoagulant repertoire in members of this family. To remedy this, herein we explore the diversity of anticoagulants and other bioactive proteins secreted by the salivary glands of L. mexicana.
Methods
Specimen collection and identification
Leech specimens were collected in Río Pilón, Rayones, Nuevo León, Mexico (25°01′31.39″N, 99°59′38.7″W) on 11–12 July 2015. Salivary gland tissue was dissected from a single specimen using sterile tools while the leech was submerged in RNAlater (Thermo Fisher Scientific). Tissue masses were stored in RNAlater at −20 °C until RNA extraction. The specimen was identified based on morphology using specialized literature (Sawyer, Reference Sawyer1986).
RNA extraction, cDNA library construction and sequencing
Total RNA was extracted using a standard TRIzol-based (Life sciences) protocol, as described below. Salivary tissue mass was frozen in liquid nitrogen, ground to powder with a drill and homogenized in 500 μL of TRIzol and the sample was incubated for 5 min at room temperature (RT). Then, 100 mL of bromochloropropane was added and the sample was mixed by vortexing. After incubation for 10 min at RT, the sample was centrifugated at 16 000 rpm for 15 min at 4 °C. The aqueous layer was recovered, mixed with 500 µL of isopropanol and incubated overnight at −20 °C. The sample was then centrifugated for 15 min at 16 000 rpm at 4 °C to allow the formation of an RNA pellet, which was then washed twice with 1 mL of 75% isopropanol and centrifuged at 16 000 rpm at 4 °C for 15 and 5 min, respectively, in the first and second washing step, air dried and eluted in 100 µL of RNA Storage solution (Ambion, Austin, TX, USA). Dynabeads mRNA purification kit (Invitrogen, Carlsbad, CA, USA) was used following the manufacturer's protocol for total RNA purification. After incubation of total RNA at 65 °C for 5 min, the sample was incubated on a rocker for 30 min with 200 µL of magnetic beads and then washed twice with washing buffer. This step was repeated after an incubation for 5 min. Thereafter, mRNA was eluted in 15 µL of Tris-HCL buffer. Quality and quantity of mRNA as measured with a pico RNA assay in an Agilent 21 000 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and quantity was also measured using a RNA assay on a Qubit fluorometer (Life Technologies, Gaithersburg, MD, USA).
The cDNA library was constructed on the Apollo 324 automated system using the PrepX mRNA kit (Wafergen, Fremont, CA, USA) following the manufacturer's protocol. Full-length cDNAs were fragmented to a mean insert size of 500 bp. The sample was inline barcoded with Illumina indices to allow for pooling with other libraries prior to sequencing. Library quality and size selection were checked on the Agilent 2100 Bioanalyzer (Agilent Technologies) with the HS DNA assay. Sequencing was performed in 1/6th of a lane on the Illumina HiSeq 2500 platform (paired end, 150 bp) at the FAS Center for Systems Biology at Harvard University, Cambridge, MA, USA.
Data sanitation, assembly and gene prediction
Before assembly, raw data were filtered and trimmed with Trim Galore! ver. 0.5 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore) applying a PHRED score of 30. Assembly was conducted with Trinity on Galaxy ver. 2.2.0 with pair distance set to 800. Open Reading Frames (ORF) among the contigs were predicted with Transdecoder ver. 5.2.0.
Annotation
Nucleotide and amino acid ORFs were BLASTed against a local anticoagulant database adapted from Kvist et al. (Reference Kvist, Brugler, Goh, Giribet and Siddall2014, Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017) and Tessler et al. (Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b) and the Swiss-Prot database, using BLASTp and BLASTx with a cut-off e-value of 1 × 10−10. In order to identify protein domains within the predicted ORFs, we searched against the Pfam database (Finn et al., Reference Finn, Coggill, Eberhardt, Eddy, Mistry, Mitchell, Potter, Punta, Qureshi, Sangrador-Vegas, Salazar, Tate and Bateman2016) using HMMER ver. 3.2.1 (Eddy, Reference Eddy2011). In addition, SignalP ver. 4.1 (Petersen et al., Reference Petersen, Brunak, Von Heijne and Nielsen2011), with the sensitive D-cutoff value, was used to determine putative signal peptides and to select ORFSs with higher e-values when compared to hits for the same ORF against the local anticoagulant database. This was done in order to ensure that hits against anticoagulants in our local database were better than any potential hit against a non-anticoagulant in any other database. The remaining ORFs were manually examined and significant hits with a higher e-value against Swiss-Prot were still reported when the hits against both databases belonged to the same protein family. ORFs were then aligned with putative anticoagulant orthologues for further phylogenetic analyses. All aligments were performed using the web version of MAFFT ver. 7 (Katoh and Standley, Reference Katoh and Standley2013) with automatic assigment of alignment strategy, the BLOSUM62 scoring matrix and gap opening penalty of 5.00. Gene trees were obtained employing the parsimony criterion in TNT ver. 1.5 (Goloboff and Catalano, Reference Goloboff and Catalano2016) with 1000 replicates, each with 5 rounds of ratcheting, drifting and tree fusing. Trees were left unrooted due to the lack of obvious outgroups for leech anticoagulants. Putative anticoagulants were then aligned with their respective archetypal anticoagulant sequence and visualized with JalView ver. 2 (Waterhouse et al., Reference Waterhouse, Procter, Martin, Clamp and Barton2009).
Results
Sequencing, assembly and BLAST searches
In total, 11 029 677 raw sequence reads were produced by the Illumina run and 65 833 contigs were assembled with Trinity – the data are deposited in the Short Read Archive at NCBI under accession number PRJNA 532 931. Transdecoder predicted a total of 1,019 ORFs with significant matches (e-value of 1 × 10−10 or below) to known anticoagulants and other salivary proteins from leeches, as well as other bloodfeeding organisms. A total of 137 hits were recovered against 15 widely known anticoagulation factors from different inhibitor families, suggesting the presence of homologues to several anticoagulants. Putative orthologues to thrombin inhibitors included hirudin and hirustasin from Hirudo medicinalis Linnaeus, 1758, bufrudin from Hirudinaria manillensis (Lesson, 1842), and an antithrombin from the king cobra Ophiophagus hannah (Cantor, 1836). Further, 27 hits against the factor Xa inhibitor antistasin from Haementeria officinalis de Phillipi, 1849 were recovered, 6 hits against therostasin from Theromyzon tessulatum (Müller, 1774) and 4 hits against ghilanten from Haementeria ghilianii de Phillipi, 1849. Fully 35 hits against inhibitors of platelet aggregation were recovered, including destabilase, saratin, bdellin-B, and disintegrins and metalloproteases with thrombospondin motifs (ADAM[TS]), as well as peptides inhibiting other pathways, such as leech derived tryptase inhibitor (LDTI). Several ORFs matched non-leech putative anticoagulants, including Kazal-type serpins and 5′-apyrase. Many putative anticoagulants possessed signal peptides at the 5′-end, suggesting their secretion out of the salivary cells – these included manillase, destabilase, antistasin, bdellin, LDTI, 5′-apyrase and therostasin. Table 1 lists all of the recovered anticoagulants and other bioactive salivary proteins related to bloodfeeding, together with their respective Pfam domains and indication of presence or absence of signal peptides.
Table 1. Putative anticoagulant proteins with high scoring matches in salivary transcripts from Limnobdella mexicana
MF, molecular function; CC, cellular component; BP, biological process.
Gene trees
We describe the unrooted tree topologies using the terminology proposed by Wilkinson et al. (Reference Wilkinson, McInerney, Hirt, Foster and Embley2007), in which a ‘clan’ in an unrooted tree is equivalent to a ‘monophyletic group’ in a rooted tree and ‘adjacent group’ is equivalent to ‘sister group’. Gene trees for multiple sequence alignments of putative antistasin-like factor Xa inhibitors, putative hirudin-like thrombin inhibitors and hirudin-like factors (HLFs), ADAM[TS], bdellin, destabilase, eglin c, ficolin, manillase and saratin are shown in Fig. 1 and Supplementary Fig. S1. Gene trees for the remaining putative anticoagulants were not inferred due to lack of comparative sequences. Hirudin-like and antistasin-like proteins showed a high degree of molecular diversity (Fig. 1A and B).
Fig. 1. Consensus tree of the most parsimonious trees of leech-derived putative antistasin-like proteins (A), hirudin-like proteins (B) and saratin/LAPP (C).
In the antistasin family tree (Fig. 1A), sequences from Limnobdella mexicana were nested in each of the three main clans that were recovered in our analysis. These represent the main lineages of putative antistasin-like factor Xa inhibitors: antistasin, therostasin and piguamerin/guamerin/bdellastasin. Within the antistasin clade, putative antistasin-like inhibitors from Piscicolidae form a clan adjacent to sequences from Glossiphoniiformes and Hirudiniformes. The archetypal antistasin and the archetypal ghilanten form a well-supported clan (99% parsimony bootstrap support [PBS]), adjacent to a sequence from Ozobranchus margoi (Apathy, 1890). Putative antistasins (DN22758 and DN19809) recovered in the L. mexicana salivary transcriptome nest as the adjacent group of sequences from arhynchobdellid species. The therostasin clan exhibits a high degree of molecular diversity and is formed by sequences recovered in previous studies from members of Oceanobdelliformes, Glossiphoniiformes, as well as one sequence from L. mexicana (DN457). The third main clan is formed by guamerin, piguamerin, bdellastasin and other sequences that presented BLAST hits against antistasin and ghilanten in previous studies. A clan, adjacent to therostasins, formed exclusively by antistasins and ghilantens from Placobdella species is also present in our analysis.
Trees for hirudin-like thrombin inhibitors (Fig. 1B) showed that the archetypal hirudin, hirulin and bufrudin form a clan that is separated from the remaining sequences by a proportionally long branch. Several hirudin-like factors form a clade, in which HLF1 cluster as adjacent to the remaining HLFs. The remaining clan is composed of hirudin-like antithrombins from a wide range of leech taxa, including both rhynchobdellids and arhynchobdellids, adjacent to L. mexicana hits against hirudin and bufrudin.
In the antiplatelet tree (Fig. 1C), there is no clear dichotomy between archetypal saratin and LAPP, suggesting that these nomenclaturally distinct anticoagulation factors may indeed be orthologous. The L. mexicana sequence is positioned on a long branch adjacent to saratin from Aliolimnatis fenestrata (Moore, 1939) and two different species of the genus Haementeria.
ADAM[TS] sequences recovered in the L. mexicana transcriptome (DN24510 and DN13743) form a clan adjacent to ADAM[TS] from Oxytonostoma typica Malm, 1863 (Supplementary Fig. S1A). In the manillase gene tree (Supplementary Fig. S1B), predicted ORFs from L. mexicana form a clan with archetypal manillase and an additional sequence from Haemadipsa interrupta (Moore, 1935). In the bdellin gene tree (Supplementary Fig. S1C), the L. mexicana (DN2497) sequence falls within a clade formed by H. medicinalis, Hirudo nipponia Whitman 1886 and the archetypal sequence of bdellin. The L. mexicana-derived eglin c putative orthologue also sits on a long branch when compared to other branches of the tree (Supplementary Fig. S1D). The hit against destabilase I falls within a clan formed by the archetypal sequence for destabilase and sequences belonging to Hirudo medicinalis and Macrobdella decora (Say, 1824) (Supplementary Fig. S1E). For ficolin (Supplementary Fig. S1F), the L. mexicana sequence and a H. medicinalis sequence form an adjacent group to another ficolin-like molecule from H. medicinalis.
Sequence similarity and pairwise alignments
Sequence similarity varied considerably between the newly-sequenced, putative anticoagulants and their respective archetypal variants. Destabilase I showed amino acid sequence similarity of 76.47%, manillase showed 59%, eglin c showed 51.43%, ficolin showed 51.43%, Kazal type serpin showed 49.12%, DN24002_c0_g1_i1 showed 39.58% against hirudin, DN24002_c0_g1_i2 showed 41.18% against hirudin, cystatin showed 43.14%, therostasin showed 34.44%, hirustasin showed 32.73%, 5′-apyrase showed 29.25%, bdellin showed 28.22%, DN19809 showed 28.16% against antistasin, DN22758 showed 25.74% against antistasin, the snake antithrombin showed 23.33%, saratin (orthologous to LAPP; see Kvist et al., Reference Kvist, Sarkar and Siddall2011) showed 19.89%, DN24901 showed 12.5% against bdellastasin, ADAM showed 10.76%, and ADAM[TS] showed 8,88%.
Pairwise alignments of selected putative anticoagulants are shown in Fig. 2 and Supplementary Fig. S2. Analysis of sequences with hits against antistasin showed important differences when compared to the archetypal anticoagulant. For example, the number of cysteines, involved in disulfide bonds, and the number of antistasin-like domains differ between the comparisons. The alignment of the best hit against antistasin-like proteins is shown in Fig. 2A. DN22758 from L. mexicana is shorter than the archetypal sequence and 12 out of the 20 cysteines present in the archetypal sequence are conserved in the L. mexicana-derived putative orthologue. One of the other hits against antistasin, ORF DN19809, shows full conservation of all 20 cysteine residues when aligned with the archetypal sequence (Fig. 2B). However, the newly generated sequence also shows the presence of a C-terminal region that is rich in cysteines and that is not shared by the archetypal antistasin. ORFs DN457 (best hit against therostasin) and DN24901 (best hit against bdellastasin) are formed by a single antistasin-like domain (Fig. 2C and D), suggesting orthology with therostasin (DN457), and piguamerin, guamerin or bdellastasin (DN24901).
Fig. 2. MAFFT-based alignment of antistasin-like proteins from Limnobdella mexicana and their respective archetypal anticoagulant. The archetypal sequence is shown on top of each alignment. A: putative antistasin (DN22758) from Limnobdella mexicana aligned with the know sequence of antistasin (AAA29193). B: putative bdellastesin/piguamerin/guamerin (DN24901) from Limnobdella mexicana aligned with the know sequence of bdellastesin (P82107). C: putative therostasin (DN457) from Limnobdella mexicana aligned with the know sequence of therostasin (Q9NBW4). D: putative antistasin (DN19809) from Limnobdella mexicana aligned with the know sequence of antistasin (AAA29193). Red boxes indicate conserved cysteines.
Pairwise alignment of the putative hirudin and the archetypal sequence (Supplementary Fig. S2A and S2B) evinced full conservation of the six cysteine residues present in both the newly generated sequence and the archetypal variant. However, the critical regions for anticoagulation activities in both N-terminal and C-terminal ends show high variation between the aligned sequences.
Discussion
The present study details, for the first time, the salivary transcriptome of a praobdellid leech. Species of this family are known for infestations in body orifices of a wide diversity of hosts. Beyond the 15 BLAST hits against well-known anticoagulants, we also found proteins that have previously been associated with anticoagulation processes, but that are not well characterized in leeches. Some of these hits include a thrombin inhibitor from the king cobra O. hannah, a Kazal type serpin and a 5′-apyrase. Similar hits have been recovered in previous studies of leech salivary transcriptomes, yet the phylogenetic comparison between the leech-derived variants and those of other organisms, coupled with low sequence similarities, suggests that these are not orthologous molecules (Kvist et al., Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017; Tessler et al., Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b). This is further exacerbated by the fact that functional characterization of these molecules is lacking for leeches. The significant BLAST-hits detailed here included inhibitors of platelet aggregation, factor Xa inhibitors and thrombin inhibitors. With regard to factor Xa inhibitors and thrombin inhibitors, whereas the initial BLAST searches of our newly generated sequences recovered crude identities, the phylogenetic analyses much refined their IDs. For example, one of the initial BLAST hits against antistasin was revealed to be a different factor Xa inhibitor (bdellastasin) by the phylogenetic analyses. This underscores the utility of phylogenetic analysis for orthology determination and the importance of including as much diversity as possible in the data matrices.
Orthology determination
Given the low sequence similarity between some of the putative anticoagulants recovered in the L. mexicana transcriptome and their respective archetypal sequences, the detailed examination of gene trees for each of the putative anticoagulants became of increasing importance. Indeed, the phylogenetic trees showed some evidence of paralogy within L. mexicana-derived putative factor Xa inhibitors (Fig. 1A), as well as within putative hirudin-like thrombin inhibitors (Fig. 1B), when compared to their highest scoring BLAST hits (Table 1).
Antistasin-like factor Xa inhibitors canonically inhibit the activation of prothrombin by factor Xa (Krowarsch et al., Reference Krowarsch, Cierpicki, Jelen and Otlewski2003; Koh and Kini, Reference Koh and Kini2009). They are formed by domains that are rich in cysteine and, whereas therostasin, guamerin and piguamerin are formed by a single domain, antistasin and ghilanten are formed by two antistasin-like domains repeated in tandem (Salzet, Reference Salzet2001; Koh and Kini, Reference Koh and Kini2009). We found, in our analyses, BLAST hits against three antistasin-like factor Xa inhibitors: ghilanten, therostasin and antistasin. Consideration of the details of the pairwise alignments of factor Xa inhibitors disclosed low sequence similarity between these ORFs, with regards to both the number of cysteines and the number of antistasin-like domains. In the antistasin family tree (Fig. 1A), ORF DN457 forms a clan together with sequences annotated as both antistasin and therostasin, and this clan also includes the archetypal sequence of therostasin. DN457 is formed by a high number of cysteines distributed in a single domain, 17 of which are shared with the archetypal sequence for therostasin (Fig. 2D). The placement of DN457 with a clan including the archetypal sequence of therostasin, coupled with the fact that our sequence seems to only include a single antistasin domain, strongly suggests orthology with therostasin.
ORF DN24901 presented a BLAST hit against antistasin, but the alignment against the archetypal sequence shows important sequence differences, including the number of cysteines involved in disulfide bonds, suggesting structural divergence. In the antisitasin family tree, DN24901 falls within a clade formed by three archetypal anticoagulants: piguamerin, guamerin and bdellastasin (Fig. 1A). The phylogenetic position of these three archetypes and the presence of a single antistasin-like domain (Fig. 2B) suggest an orthologous relationship with one of these anticoagulants. Following the results of the phylogenetic tree, we hypothesize that this sequence is orthologous with bdellastasin, but detailed functional characterization of these poorly understood proteins is needed prior to any definitive statement.
As in previous studies (e.g. Kvist et al., Reference Kvist, Brugler, Goh, Giribet and Siddall2014), the archetypal sequences for antistasin and ghilanten form a well-supported clan (PBS = 100%), suggesting that these are orthologous sequences that could be circumscribed by a single name (see Fig. 1A). Antistasin was first isolated from the glossiphoniid leech Haementeria officinalis and, later, ghilanten was isolated from Haementeria ghilianii and showed higher anticoagulant potency (Dunwiddie et al., Reference Dunwiddie, Thornberry, Bull, Sardana, Friedman, Jacobs and Simpson1989; Blankenship et al., Reference Blankenship, Brankamp, Manley and Cardin1990). The alignment of ghilanten and antistasin shows a high degree of sequence similarity between the C-terminal region of both sequences, with most variations lying on the N-terminal region (data not shown); the active site of antistasin is located in the N-terminal region, and the amino acid differences found between both anticoagulants were previously linked to differences in anticoagulatory power (Dunwiddie et al., Reference Dunwiddie, Thornberry, Bull, Sardana, Friedman, Jacobs and Simpson1989; Blankenship et al., Reference Blankenship, Brankamp, Manley and Cardin1990). Given the phylogenetic position of antistasin and ghilanten and their high sequence similarity, coupled with their identical function, we suggest that ghilanten should be recognized as an orthologue of antistasin. Two putative antistasin sequences (DN22758 and DN19809) were recovered in L. mexicana (Fig. 2A and D). These two sequences are adjacent to putative antistasin recovered from H. medicinalis, Hirudo verbana Carena, 1820, M. decora, A. fenestrata and H. interrupta, as well as the archetypal sequences for ghilanten and antistasin (Fig. 1A).
Two putative hirudin-like thrombin inhibitors were recovered in our BLAST search. ORFs DN24002_c0_g1_i1 and DN24002_c0_g1_i2 recovered BLAST hits against bufrudin and hirudin, respectively. Two regions are critical for the inhibition of thrombin by hirudin-like thrombin inhibitors. The first three amino acids of the N-terminal regions bind to thrombin, penetrating into its active site. In addition, the last 10 amino acids of the C-terminal region of hirudin are important residues involved in the specific binding with the exosite. As a result, changes to those amino acids cause critical alteration to its binding efficiency (Sharp, Reference Sharp1996; Salzet, Reference Salzet2002). Recently, the discovery of hirudin-like factors, proteins similar to hirudin but that lack anticoagulant capabilities, has raised questions about the capacity of hirudin-like antithrombin proteins for inhibiting coagulation (Müller et al., Reference Müller, Mescke, Liebig, Mahfoud, Lemke and Hildebrandt2016, Reference Müller, Haase, Lemke and Hildebrandt2017). Our alignments show a low degree of amino acid residue conservation in the critical C-terminal region in our hits against hirudin and bufrudin, when compared to the archetypal hirudin. Previous studies have demonstrated that a high level of variation within this region is common among hirudin-like proteins, even among proteins that show anticoagulatory activities (Scacheri et al., Reference Scacheri, Nitti, Valsasina, Orsini, Visco, Ferrera, Sawyer and Sarmientos1993; Müller et al., Reference Müller, Mescke, Liebig, Mahfoud, Lemke and Hildebrandt2016, Reference Müller, Haase, Lemke and Hildebrandt2017; Siddall et al., Reference Siddall, Brugler and Kvist2016). Our recovered hits are highly similar (84.54%) to each other, and they fall within a clan formed by a polytomy of putative antithrombin inhibitors and hirudin-like factors. Further investigations are needed to clarify whether sequences within this polytomy are paralogous or orthologous and if they possess any anticoagulatory activity. Furthermore, hirudins (including the archetypal sequence) from H. verbana, and H. medicinalis, as well as the archetypal bufrudin from H. manillensis form a clan (Fig. 1B), indicating paralogy between hirudins and hirudin-like factors, as suggested by Müller et al. (Reference Müller, Haase, Lemke and Hildebrandt2017).
BLAST searches also recovered hits against numerous putative inhibitors of platelet aggregation. In our bdellin, destabilase I and manillase phylogenetic reconstructions (Supplementary Figs S1B, S1C and S1E), the ORFs with significant hits against those anticoagulants nested within clades formed by species of Hirudiniformes, always including the archetypal sequences, suggesting orthology between our newly generated sequences and the archetypal anticoagulants. Saratin, ADAM[TS], eglin c and ficolin showed low sequence similarity in pairwise alignments with their respective archetypal anticoagulant (Supplementary Figs 2H, 2I). Saratin, is an inhibitor of the von Willebrand factor-mediated platelet aggregator that has been proposed to be orthologous to leech antiplatelet protein (LAPP) (Barnes et al., Reference Barnes, Krafft, Frech, Hofmann, Papendieck, Dahlems, Gellissen and Hoylaerts2001; Cruz et al., Reference Cruz, Eidt, Drouilhet, Brown, Wang, Barnes and Moursi2001; Kvist et al., Reference Kvist, Sarkar and Siddall2011). The placement in our analysis of saratin from L. mexicana as adjacent to LAPP sequences from glossiphoniid species and the presence of a clade formed by the archetypal saratin, saratin sequences from H. medicinalis and M. decora, and one LAPP sequence from H. medicinalis, suggest orthology between LAPP and saratin (Fig. 1C).
ADAM and ADAM[TS] are members of the metalloproteinase superfamily. Both of these proteins have been linked to biological processes, including the cleavage of the von Willebrand factor (ADAM[TS]13) that plays a major role in blood coagulation (Tang, Reference Tang2001; Furie and Furie, Reference Furie and Furie2005; Porter et al., Reference Porter, Clark, Kevorkian and Edwards2005; Brocker et al., Reference Brocker, Vasiliou and Nebert2009). Putative metalloproteinase proteins have been recorded in some leech salivary transcriptomes (Kvist et al., Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017; Tessler et al., Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b), but no functional characterizations have been performed and their role in leech salivary secretions is still unknown. Hits against both families of metalloproteinases were recovered with low sequence similarity against their respective archetypal anticoagulants (ADAM = 10.76% and ADAM[TS] = 8.88%). In our phylogenetic reconstruction, putative ADAM and ADAM[TS] from L. mexicana form a clan, adjacent to a sequence from Oxytonostoma typica (Supplementary Fig. S1A).
Ficolins are reported in numerous mammals, in which they play an important role as pattern-recognition receptors of pathogens in the immune system (Fujita et al., Reference Fujita, Matsushita and Endo2004). Non-immune functions of ficolin-like proteins have been reported in snake venoms, where they seem to act as a trigger for platelet aggregation (Ompraba et al., Reference Ompraba, Chapeaurouge, Doley, Devi, Padmanaban, Venkatraman, Velmurugan, Lin and Kini2010). Regarding bloodfeeding animals, ficolin has been reported in several leech and arthropod salivary transcriptome studies, which suggests a regulatory role in coagulation during bloodfeeding (Ribeiro et al., Reference Ribeiro, Arcà, Lombardo, Calvo, Phan, Chandra and Wikel2007; Ribeiro and Arcà, Reference Ribeiro and Arcà2009; Min et al., Reference Min, Sarkar and Siddall2010; Kvist et al., Reference Kvist, Brugler, Goh, Giribet and Siddall2014; Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017). As in those studies, we suggest that L. mexicana possesses proteins orthologous to leech ficolins, based on the fact that our best hit nests within a clade formed by two sequences from Hirudo medicinalis (Supplementary Fig. S1F).
Eglin c is an inhibitor of elastase that has been extensively considered as a candidate for drug development (Salzet, Reference Salzet2001). Previous studies have demonstrated low sequence similarity between putative eglin c, sometimes to the inclusion of various cysteine residues and sometimes in their complete absence (Min et al., Reference Min, Sarkar and Siddall2010; Kvist et al., Reference Kvist, Min and Siddall2013). In L. mexicana, the best hit against eglin c possesses a cysteine in position 47 that is not shared with any other eglin c from other leech species. The phylogenetic reconstruction of leech eglin c showed that the L. mexicana hit forms a clan with a sequence from H. interrupta that also possesses a single cysteine, but in a different position. The archetypal eglin c forms a separate clan with sequences from A. fenestrata and H. medicinalis (Supplementary Fig. S1D).
Although our analyses clarify the phylogenetic relationships of some major leech anticoagulants, especially antistasin-like and hirudin-like anticoagulants, we recognize that the position of the root in a tree may affect orthology determination (see Ballesteros and Hormiga, Reference Ballesteros and Hormiga2016); if the root of a tree is placed within a monophyletic assemblage this assemblage will be rendered paraphyletic after rooting. Since there is no objective way to determine the root for the phylogeny of these major putative leech anticoagulants, coupled with the fact that no study has investigated the presence of putative orthologues within non-hirudinean clitellates, we left the trees unrooted. Given the importance of the root for estimations of orthology using phylogenies, future efforts should focus on finding a logical root for the trees of major leech anticoagulants.
Association of host specificity and anticoagulant diversity
Previous studies on leech salivary transcriptomes have suggested a positive association between the number of anticoagulants and host-generality (Kvist et al., Reference Kvist, Brugler, Goh, Giribet and Siddall2014). Therefore, it is hypothesized that host-generalist species have a wider anticoagulant repertoire compared to host-specialists. However, Tessler et al. (Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b) found no significant pattern between piscicolid leeches restricted to elasmobranch hosts and generalist leeches of the same family, though the number of anticoagulants found in elasmobranch-restricted parasitic leeches was slightly lower. Limnobdella mexicana is a member of Praobdellidae, a family that is known to be associated with nasopharyngeal infections; host records for members of this family include mammals, reptiles, amphibians and invertebrates (Phillips et al., Reference Phillips, Arauco-Brown, Oceguera-Figueroa, Gomez, Beltrán, Lai and Siddall2010; Siddall et al., Reference Siddall, Rood-goldman, Barriot and Barboutist2013; Nakano et al., Reference Nakano, Tomikawa, Sakono and Yoshikawa2017). The number of anticoagulants found in L. mexicana was slightly higher than most of the previous studies of leech salivary transcriptomes, except for that of H. interrupta. However, comparison between different available datasets is exacerbated by the fact that sequencing efforts and technologies vary among different studies. Moreover, previous studies have focused on the presence of homologues of archetypal anticoagulants, not accounting for the presence of multiple duplicates (Min et al., Reference Min, Sarkar and Siddall2010; Kvist et al., Reference Kvist, Brugler, Goh, Giribet and Siddall2014; Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017; Siddall et al., Reference Siddall, Brugler and Kvist2016) that might possess minute sequence and structural differences such as those recorded for hirudin variants, specifically hirudin-like factors (Scacheri et al., Reference Scacheri, Nitti, Valsasina, Orsini, Visco, Ferrera, Sawyer and Sarmientos1993; Müller et al., Reference Müller, Haase, Lemke and Hildebrandt2017). Nonetheless, the sequencing effort and determinations of putative anticoagulants applied in the present study is exactly the same as in Kvist et al. (Reference Kvist, Oceguera-Figueroa, Tessler, Jiménez-Armenta, Freeman, Giribet and Siddall2017). Thus, the higher number of putative anticoagulants found might be an indication of a positive association between host diversity and diversity of anticoagulants in leeches, but further investigations, putatively involving full genome sequencing, are still needed to confirm such a relationship.
Implication for evolution of bloodfeeding in leeches
The origin of bloodfeeding in leeches has been widely debated. Sawyer (Reference Sawyer1986) suggested a bloodfeeding ancestor for proboscis-bearing leeches, but a macrophagous ancestor for Arhyncobdellida. However, more recent phylogenetic studies have suggested bloodfeeding as the ancestral behaviour of Hirudinida, with independent losses within proboscis and non-proboscis-bearing lineages (Apakupakul et al., Reference Apakupakul, Siddall and Burreson1999; Tessler et al., Reference Tessler, de Carle, Voiklis, Gresham, Neumann, Cios and Siddall2018a). Transcriptomic data have been used to test the hypothesis of a bloodfeeding ancestor of leeches, under the perception that non-bloodfeeding leeches will maintain anticoagulation factors based on this ancestral donation, as well as the fact that bloodfeeding leeches across the leech phylogeny should share orthologous proteins if these were derived from a common ancestor. The presence of orthologues in Glossiphoniidae and Piscicolidae, when compared to Hirudiniformes, for each of hirudin, mannilase, piguamerin, bdellastasin, destabilase and LDTI, recorded in previous studies, suggest that these were derived from an ancestor of modern leeches that approaches the time frame of the ancestor of all leeches (Siddall et al., Reference Siddall, Brugler and Kvist2016; Tessler et al., Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b). In other words, it has been suggested that if (e.g.) Piscicolidae shares orthologous anticoagulant proteins with Hirudinidae, the ancestor that provided these proteins existed early in leech evolution, approaching the base of the leech tree. However, previous studies that recorded the presence of these anticoagulants in many Piscicolidae and Glossiphoniidae species have approached the problem of orthology determination using a simple BLAST search or performed phylogenetic analysis that did not include the majority of the sequences available for the anticoagulant families studied (Siddall et al., Reference Siddall, Brugler and Kvist2016; Tessler et al., Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b).
Comparison between BLAST hits and gene trees of their respective putative orthologous proteins showed incongruence of antistasin-like and hirudin-like putative anticoagulants. Therefore, previous records of these anticoagulants can be questioned. Our phylogenetic analysis of putative anticoagulants found in L. mexicana contributes to the orthology determination of several families of leech anticoagulants and, consequently, serves to clarify the origin of these major protein lineages. Our antistasin family tree shows that putative antistasin, bdellastasin and therostasin are found in both proboscis- and non-proboscis-bearing leeches, suggesting that the ancestral leech possessed the molecular armoury to inhibit factor Xa, even if this inhibition is not directly associated with bloodfeeding.
On the other hand, hirudin derived from hirudiniform species form a clan in the hirudin gene tree, and this clan is clearly separated from hirudin-like factor (HLFs) and other putative thrombin inhibitors from proboscis- and non-proboscis-bearing leeches. Since HLFs do not possess anticoagulatory capabilities and their function is still unknown (Müller et al., Reference Müller, Mescke, Liebig, Mahfoud, Lemke and Hildebrandt2016, Reference Müller, Haase, Lemke and Hildebrandt2017), not much can be inferred about the anticoagulatory activity of hirudin-like proteins in the ancestral leech, as suggested by Siddall et al. (Reference Siddall, Brugler and Kvist2016) and Tessler et al. (Reference Tessler, Marancik, Champagne, Dove, Camus, Siddall and Kvist2018b). Our analyses also present evidence of orthology between hits for manillase and LAPP/saratin recovered from early and late branching lineages of leeches, pointing to a possible single origin of these proteins. Therefore, with evidence detailed above, the sole presence of a BLAST hit against a putative anticoagulant should not be interpreted as evidence for the origin of certain molecular function. BLAST homology determination should be followed by robust orthology tests.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182019000593
Author ORCIDs
Rafael Eiji Iwama, 0000-0002-2430-4178
Acknowledgements
The authors acknowledge Alejandro Manzano-Marín for assistance with software usage and Michael Tessler for the curation of the anticoagulant local data used on BLAST searches and for creating the script used to run TNT. We thank Rosa Fernández for technical assistance during RNA extraction and we thank the Willi Hennig Society for making TNT freely available. Two anonymous reviewers provided comments that help to make this article clearer.
Financial support
This project was funded by a NSERC Discovery Grant and an Olle Engkvist Byggmästare Foundation scholarship (SK), as well as internal funds from the Faculty of Arts and Sciences, Harvard University (GG).
Conflicts of interest
None.
Ethical standards
Not applicable.
