Introduction

Inventorying strategies are key to successful conservation programs. Overall species composition in protected areas helps to identify management priorities and informs the implementation of usage programs [1]. Furthermore, the success of a particular conservation effort [2] relies on effective longitudinal biodiversity monitoring in ways that also allow for evaluation of environmental change, ecosystem health, and the emergence of pathogens [3, 4].

Live trapping, line transects, and track surveys are common inventory methods for terrestrial tetrapods, each with its own inherent limitations and biases [5, 6]. Passive camera trapping is often considered more efficient especially in terms of human resources and expertise [7, 8]. Insofar as the likelihood of a species triggering a sensor is inversely proportional to body size [9], camera trapping the majority of mammals requires specific modifications [10] or the use of intermittent video capture instead [11]. Meanwhile reptiles and amphibians, which are more imperiled than mammals [12, 13], are represented by comparatively few camera trap studies [14]. Of the foregoing, only live-trapping and active capture allow for comprehensive genetic biodiversity monitoring of the fullest range of vertebrate diversity in protected areas while allowing for demographic studies and the detection of cryptic species [3, 15, 16].

In aquatic habitats, passive environmental DNA (eDNA) sampling has proven robust when employed in a longitudinal design [17] and useful in detecting invasive and rare aquatic species [18, 19]. Though eDNA is ill-suited to terrestrial monitoring, a burgeoning use of an invertebrate ingested DNA (iDNA) approach shows promise [20, 21, 22, 23]

Sampling iDNA from leeches encompasses some of the benefits of other survey strategies [24] for those leech groups that have both sanguivorous and terrestrial life-history strategies [25]. Much like investigators in track and transect surveys, leeches are actively sampling their environment, but unlike biting flies leeches do not travel large distances between blood meals [26]. Residual mitochondrial DNA remains detectable in the gut of leeches for months post-feeding [27] and it is becoming clear that terrestrial leeches are particularly adept at revealing small animal diversity in ways otherwise only possible with trapping and active collecting [26].

Each of cytochrome b, 12S rDNA, 16S rDNA and cytochrome c oxidase subunit 2 gene have been explored as a barcode locus in iDNA isolates from mosquitoes, carrion flies, biting midges, ticks and leeches [21, 27, 28, 29, 30]. The mitochondrial ribosomal genes are among the more taxonomically universal markers [31], though rarely are they used in combination. Here we explore several iDNA strategies toward the recovery iDNA in recent as well as >10 year old museum collections of terrestrial leeches while also assessing leech gut associated bacterial symbionts and vertebrate trypanosome parasites detectable in blood meals.

Materials and methods

Results

Among the 761 individual leeches examined (S1 Table) with conventional sequencing methods, 437 did not yield an amplification product either for 16S rDNA or for 12S rDNA. Excluding highest-scoring sequence matches to humans, cattle, water buffalo, pig, chicken and house-rats, there were 39 confirmed high-scoring (better than 1e^-20) sequences from 12S rDNA. Twenty-one (54%) of those samples also previously yielded matches to non-human and non-domesticated animals on the basis of 16S rDNA. Altogether, 12S rDNA and 16S rDNA data combined were dominate by carnivores, primates and artiodactylids (Fig 1B). The proportion of leeches from which DNA was extracted without regard to the visual presence or absence of bloodmeals that yielded sequence information relative to animals of interest (15%) well-exceeds the 5% of leeches observed with visible blood meals in legacy collections examined here with next-gen sequencing (below). Average sequence identity to high scoring matches in nr/nt was 97.7% for 12S rDNA and 95.9% for 16S rDNA.

Among those isolates conventionally sequenced for 12S rDNA, tetrapod species identified came from three classes (S1 Table): Mammalia, represented by six orders (Artiodactyla, Carnivora, Chiroptera, Eulipotyphla, Primates, Rodentia); Aves, represented by two orders (Galliformes; Strigiformes); and Amphibia, represented by one order (Anura). Among those isolates previously conventionally sequenced for 16S rDNA, tetrapod species identified came from two classes: Mammalia, represented by five orders (Artiodactyla, Carnivora, Chiroptera, Primates, Rodentia); and Aves, represented by one order (Galliformes). In Bangladesh, most iDNA isolates were primate (38%, esp. macaque), carnivore (31%, esp. civets) or ungulate (21%, esp. muntjak). Similarly, across Cambodia, most isolates were carnivore (33%, esp. civet), primate (22%, esp. macaque) or ungulate (22%, esp. red deer and mouse deer). Isolates from China were predominantly ungulate (30%, esp. muntjak and red deer), rodent (24%, esp. white bellied rat) or carnivore (18%, esp. civet).

In a manner not apparent with the mammal specific 16S rDNA locus, individual sample amplifications and conventional sequencing of the 12S rDNA locus provided results that included amphibians like a microhylid from Cambodia (e-value = 1e^-14, 98% identity) and a species of Fejervarya from China (e-value = 6e^-52, 99% identity). While it also found matches to birds, none of these presented convincing scores both for e-value and sequence identity. Rodents also were better represented among 12S rDNA sequences (24%) than among 16S rDNA data (2%) from China, Bangladesh and Cambodia combined. In the case of both loci, many species were recovered only once: a Northern treeshrew (Tupaia belangeri, e-value = 8e^-114, 96% identity) from Bangladesh, each of two squirrels (Menetes berdmorei, e-value = 5e^-29, 96% identity; Proscirillus murinus, e-value = 8e^-22, 100% identity), a gaur (Bos gaurus, e-value = 5e^-133, 99% identity), and a roundleaf bat (Hipposideros bicolor, e-value = 4e^-73, 93% identity) from Cambodia, as well as one squirrel (Dremomys rufigenis, e-value = 3e^-106, 93% identity), a mole (Scapanulus oweni, e-value = 2e^-20, 97% identity), and a silver pheasant (Lophura nycthemera, evalue = 1e^-125, 96% identity) from China.

It was apparent that method of field preservation influenced amplification success. All samples from China and Bangladesh had been fixed in RNAlater; samples from Cambodia had been preserved in either RNAlater or in 95% ethanol. Of the 55 samples representing target wild taxa, 39 had been fixed in RNAlater and 16 in 95% ethanol, and of the 23 Cambodian samples providing vertebrate 12S rDNA sequences, only 2 had been preserved in ethanol. Regardless, amplification and sequencing of 12S rDNA was confirmatory of prior 16S rDNA results in 7 samples from Bangladesh, 11 samples from Cambodia, and 3 samples from China. It also provided additional data for 3 samples from Bangladesh, 12 samples from Cambodia, and 4 samples from China, none of which previously yielded data for 16S rDNA.

The vertebrate 12S locus could not discriminate among species of macaque from Bangladesh or Cambodia (see S1 Table), whereas the 16S locus was able to distinguish Macaca mulatta and Macaca sinensis with sequence identities of 97% or higher. Conversely, scores of 87% to 95% from 16S rDNA suggesting African civets (Civettictis civetta) in Bangladesh which were properly identified as the Large Indian civet (Vivera zibetha) at 100% sequence identity with 12S rDNA. Nonetheless, there were also disagreements between loci. A herpestid from Bangladesh was determined by 12S rDNA to be Herpestes javanicus (evalue = 1e^-43, 99% identity), not an African Crossarchus obscurus (evalue = 1e^-107, 93% identity) as suggested by 16S rDNA. In contrast, a rodent from Cambodia was more plausibly identified with 16S rDNA as Niviventer excelsior (evalue = 4e^-41, 92% identity) as opposed to being an Amazonian Oecomys auyantepui (evalue = 1e^-24, 95% identity) as suggested by 12S rDNA. From China, porcupine iDNA from one sample was either Athurus macrourus (evalue = 1e^-30, 100% identity on 12S rDNA) or Hystrix indica (evalue = 1e^-105, 93% identity on 16S rDNA), both of which are geographically plausible.

Genomic iDNA next-gen sequences from the Australian leech C. tanae yielded significant scores against a pademelon (Thylogale billardierii: tRNA-Leu+nd1, evalue = 3e^-142), humans (Homo sapiens: 12S, evalue 0.0; cox3, evalue = 0.0; cytb+tRNA-Thr, evalue = 0.0), and a domestic dog (Canis lupus familiaris: cox2, evalue = 0.0; nd4, evalue = 2e^-159; nd5, evalue = 0.0). Reads from the Malaysian leech H. interrupta matched those of a mouse deer (Tragulus javanicus: tRNA-Phe, 6e^-101; tRNA-Trp+tRNA-Ala+tRNA-Asn+tRNA-Cys+tRNA-Tyr+cox1, evalue = 0.0; cox1, evalue = 0.0; cox2, evalue = 0.0; nd3+tRNA-Arg+nd4L, evalue = 0.0; nd2 evalue = 0.0; nd4 evalue = 0.0; nd4+tRNA-His+tRNA-Ser+tRNA-Leu+nd5 evalue = 0.0) with an overall 97.9% sequence identity. Those from the Madagascan C. fallax returned significant matches for the lesser hedgehog tenrec (Echinops telfairi: cox1, evalue = 1e^-123), for the ring-tailed mongoose (Galidia elegans: nd2, evalue = 0.0; cytB, evalue = 0.0; tRNA-Thr+tRNA-Pro+control region, evalue = 0.0) with an overall 98.8% sequence identity, and in particular for the narrow-striped mongoose (Mungotictis decemlineata: nd5, evalue = 0.0; nd4, evalue = 0.0; atp6+cox3, evalue = 0.0; tRNA-Phe+12S evalue = 1e^-159; nd1, evalue = 2e^-167; cox1, evalue = 0.0; cox1+tRNA-Ser+tRNA-Asp+cox2, evalue = 0.0; nd5+nd6, evalue = 0.0) with an overall 91.5% sequence identity.

Thirteen isolates amplified and sequenced conventionally yielded 18S rDNA reads for trypanosomes. One resulted from a specimen of Haemadipsa rjukjuana (specimen 5B3) from Goaligong, China, that had human 16S rDNA in its gut. The remaining 12 were specimens of Haemadipsa ornata from Bangladesh. Two evidenced muntjak iDNA (3C1, 4H10), three had cow iDNA (3E5, 3B2, 4C2), two had pig iDNA (4F1, 4C2), one had macaque iDNA (3A4), and one evidenced mongoose iDNA (4H4). Next-generation sequence reads resulting from amplifications of trypanosome 18S rDNA only yielded best matches from nr/nt for pooled samples of C. tanae: 72 reads assembling into 2 contigs and 11 singletons. Phylogenetic analysis of all of the foregoing conventionally and next-generation sequenced 18S rDNA reads in combination with those of putatively related species of Trypanosoma revealed four principal clades (Fig 2), all of which were relatively closely related to Trypanosoma theileri and Trypanosoma cyclops (clade C in Fig 2). Conventional reads from leeches that had fed on a macaque (3A4) grouped with Trypanosoma cyclops (clade C in Fig 2) previously only known from macaques, as did isolate 4H4 that had fed on a mongoose. Conventional reads from leeches that had fed on muntjak (3C1 and 4H10) grouped together in a clade (clade A in Fig 2) with others from leeches that had cow iDNA and with reads obtained from next generation sequencing of C. tanae. However, the clade with the most reads obtained from next generation sequencing of C. tanae grouped with an unnamed trypanosome isolate previously known from a wallaby (clade D in Fig 2).

Next-generation sequence reads resulting from amplifications of tetrapod 12S rDNA from C. tanae numbered 77 that assembled into 9 contigs, of which 7 were human and one was a dog. Next-generation sequence reads resulting from amplifications of bacterial 16S rDNA from C. tanae numbered 143 that assembled into 7 contigs, leaving 9 singletons, and yielding best scoring matches (evalue < 2e^-21) to cultured and uncultured microbial taxa: Corynebacteriales, Rikenellaceae, beta-proteobacteria, alpha-proteobacteria, Prevotellaceae, and Firmicutes. The contig representing the largest proportion (45%) of reads from amplifications, and an additional single read left out of all contigs, were both closely related to other Rikenellaceae previously recovered from the digestive tract of Australian and African aquatic leeches in the family Hirudindae (Fig 3).

In terms of best scoring genome-wide matches by genomic iDNA contigs to available microbial sequences in nr/nt, the most frequently recovered (Fig 3) were species of Mesorhizobium, Aminobacter, and other Rhizobiales (44%, median evalue = 2e^-113), as well as Delftia species and other Burkholderiales (23%, median evalue = 4e^-26). The remainder was distributed across alpha-proteobacteria (8%), beta-proteobacteria (8%), gamma-proteobacteria (8%), Actinobacteria (4%), and nitrogen fixing Rhodospirillaceae (4%).

Discussion

Our findings reinforce the notion that amplifiable DNA in legacy museum specimens can reveal aspects of an organisms’ biology and ecological interactions with the surrounding community [3]. Specifically, these terrestrial leech specimens were able to provide information regarding geographically co-occurring tetrapods and their blood parasites as well as for the leeches themselves and their associated microbiomes. Significantly, this study provides further insight into best practices regarding standardizing iDNA protocols [24, 26, 27, 32]. Unlike studies leveraging blood-feeding dipterans [23] visible blood meals in the crop of leeches appear to be unnecessary for the success of iDNA biodiversity surveys. And while new studies might consider the superiority of RNAlater over ethanol preservation in the field for de novo comprehensive surveys, it is clear that greater-than decade-old ethanol-preserved samples decade-old retain sufficient mitochondrial DNA to be of complementary value.

As inferred elsewhere [23] amplification and sequencing of both 12S and 16S loci is preferable to a single-locus iDNA approach however, study designs for indirect assessments of biodiversity, be that through iDNA surveys of leeches or eDNA from water sources [42], still need to account for the relative unevenness of coverage of reference databases. From the next-generation shotgun genomic sequencing approach there was a precise lack of overlap among mitochondrial loci obtained and identified as Malagasy mongoose species: Galidia elagans versus Mungotictis decemlineata. Where loci were represented for Galidia elegans in public reference databases, those received strong (>98%) matches from next-gen data; and where those loci did not exist for G. elegans, next-gen data matched loci in the mitochondrial genome of M. decemlineata exceedingly well (>91%). Thus, we conclude that the leech had fed on a G. elegans not on a M. decemlineata notwithstanding a larger number of sequences matching a larger number of loci for M. decemlineata, because the reference database includes fewer loci for the former. Had we only sequenced a single target locus, including the cox1 barcoding locus, we would have readily concluded otherwise. Nor for that matter does a lower e-value on one locus necessarily imply higher confidence in taxonomic identification relative to a worse e-value on a different locus especially without reference to how close a match the next-best taxon might be. Without thorough knowledge of the structure and completeness of the reference databases against which environmental samples are to be compared, blithe and uncritical acceptance of high-scoring matches is bound to be misleading in ways that are worse than merely uninformative. This should also discourage any clustering of sequencing reads into MOTUs prior to taxonomic assignment [43] with an arbitrarily chosen level of similarity. The mean sequence identity recovered against Tragulus javanicus for more than half a dozen loci exceeded the 97% identity measure often used for clustering MOTUs regardless of the fact that this species is not known to exist in Malaysia. Yet scores against available data for the Malaysian Tragulus kanchil were worse. We take this to be indicative of a different Malaysian taxon (e.g., Tragulus napu) for which no comparable data exist. The possibility of acquiring sequence information for undescribed taxa should not be ignored [44].

Any source of direct or indirect data on forest biodiversity entails unique biases (Fig 1). Leech iDNA, while more adept at recovering smaller ground-dwelling tetrapods than camera traps [32], are presumed to be insensitive to arboreal-dwelling species (e.g, lemurs) where unlike dipterans [23] leeches spend little time. Nonetheless, it is apparent from the combined 12S and 16S data that exploiting the gut contents of tropical terrestrial haemadipsid leeches has the potential to sample non-mammalian and mammalian diversity alike. Next-generation amplicon studies employing additional primer sets in combination that are exclusive and specific to amphibians, squamates or birds might yet recover more diversity than was apparent here [18].

Natural history collections serve not only as fundamental libraries of spatiotemporally dynamic patterns of biodiversity, for example of vertebrate taxa in regions exposed to climate change, but also serve as an invaluable resource relating to their associated parasites, pathogens and microbiotas [45]. In many cases, first-appearances, historical distributions, and range-shifts of parasitic helminths and protists, have only be recovered from tissues and related information originally collected with the intent of the host alone [46, 47]. Understanding parasite diversity has itself been driven by retrospective assessments of existing collections of their hosts [48, 49, 50]. Trypanosomes are especially well retained in and amplifiable from bloodfeeding vectors like tsetse flies [51]. Fig 2, based on isolates from leeches, suggests that there is a greater diversity yet to be characterized and taxonomically described from terrestrial Australian tetrapods, particularly as it concerns clades A and B. Similarly, amplifiable trypanosome DNA was recovered from 25 of 44 terrestrial leeches sampled in Australia and New Guinea [52], much like finding trypanosomes in the digestive tracts of leeches up to 44 days after feeding [53]. None of this, however, should be taken as evidence of leeches playing a functional role in vectoring trypanosomes among terrestrial tetrapods. Regrettably, ozobranchid leeches have been implicated as a vector of marine turtle herpesvirus [54, 55] in light of finding viral DNA in their gut contents yet without any positive evidence of transmission by those leeches. Likewise, trypanosomatids are demonstrably viable even in vitro with minimal nutritional provision [56, 57]. As such their presence in the crop of any invertebrate should not be taken as evidence of the latter’s status as a viable infection vector. Nevertheless, whatever diversity can be characterized from iDNA is meaningfully representative of the hematozoa resident in geographically coincident vertebrate hosts.

Outside of serious infection from resident gut flora of aquatic leeches [58, 59, 60] there has been no prior evaluation of the microbiome of terrestrial haemadipsid leeches notwithstanding their frequent ectoparasitism of humans [61, 62, 63]. Metagenomic data obtained here for pooled samples of C. tanae from Australia suggest that species of Aeromonas uniformly found in the crop of “medicinal” leeches [64, 65, 66, 67] and implicated in human infections are not well-represented in the crop of terrestrial leeches. Other stereotypical leech gut related symbionts from Rikenellaceae [39, 68, 69] are present in terrestrial leeches, with isolates grouping biogeographically (Fig 3). Lineages obtained here from C. tanae group with other Gondwanan (African and Australian) taxa. There is also an apparent unambiguous symbiotic association for Rhizobiaceae with heamadipsid leech crop contents in a manner that is strikingly similar to what was discovered for obligate intracellular mycetome endosymbionts of freshwater leech species in the genus Placobdella [28, 38]. Presently, unlike the intracellular bacterial symbionts of Placobdella species there does not appear to be a tissue-specific relationship for these microbes that otherwise are best known for symbioses in leguminous plants. Further careful morphological work will be required to assess whether haemadipsid crop-associated “Lambertion organs” [70] are the location of these bacteria and whether, like the esophageal bacteria of the genus Placobdella, the bacteria retain all of the functional iron-dependent nitrogen-fixing genetic loci.

Supporting information