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.
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.
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.
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).
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).
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.
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.
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.  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.
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.