DNA Barcode Identification of Freshwater Snails in the Family Bithyniidae from Thailand

Freshwater snails in the family Bithyniidae are the first intermediate host for Southeast Asian liver fluke (Opisthorchis viverrini), the causative agent of opisthorchiasis. Unfortunately, the subtle morphological characters that differentiate species in this group are not easily discerned by non-specialists. This is a serious matter because the identification of bithyniid species is a fundamental prerequisite for better understanding of the epidemiology of this disease. Because DNA barcoding, the analysis of sequence diversity in the 5’ region of the mitochondrial COI gene, has shown strong performance in other taxonomic groups, we decided to test its capacity to resolve 10 species/ subspecies of bithyniids from Thailand. Our analysis of 217 specimens indicated that COI sequences delivered species-level identification for 9 of 10 currently recognized species. The mean intraspecific divergence of COI was 2.3% (range 0-9.2 %), whereas sequence divergences between congeneric species averaged 8.7% (range 0-22.2 %). Although our results indicate that DNA barcoding can differentiate species of these medically-important snails, we also detected evidence for the presence of one overlooked species and one possible case of synonymy.


Introduction
Molecular taxonomic methods have been used extensively to complement morphological approaches for species identification, and for establishing phylogenetic relationships [1][2][3][4][5][6][7][8][9][10]. Particularly, species identification through DNA barcoding has seen rapid adoption over the past decade. Prior DNA barcode studies have clearly established their effectiveness in the delimitation of animal species, and also contributed several advantages [11][12][13]. The ability of DNA barcoding to identify all life stages has particular importance in medical parasitology, where it is not only important to identify the parasite and its final host, but also all its life stages and its intermediate hosts. Thus, a multidisciplinary method of classification that includes morphological, molecular and distributional data is an essential prerequisite for understanding the epidemiology of any parasite-induced disease [7].
Freshwater snails of the family Bithyniidae serve as intermediate hosts for the liver fluke, Opisthorchis viverrini, and related species common in the Greater Mekong subregion (Cambodia, Lao People's Democratic Republic, Vietnam, and Thailand). The infection of this parasite has been associated with several hepatobiliary diseases, including opisthorchiasis, cholangitis, obstructive jaundice, hepatomegaly, cholecystitis, and biliary lithiasis [14][15][16][17][18]. Furthermore, both experimental and epidemiological evidence suggest that liver fluke infections can be an etiological factor of cholangiocarcinoma [19][20][21][22][23][24][25]. Three taxa of Bithynia are involved in the transmission of this parasite [26][27][28] with different species reported as intermediate hosts in different parts of Thailand. B. siamensis goniomphalos is a dominant host in the northeast, while B. funiculata and B. siamensis siamensis serve as hosts in the north and B. siamensis siamensis in the central region [26,29]. Taxonomic keys for differentiation to species in the family Bithyniidae utilized size, shape, color, and sculpture on the shell surface, operculum structure, and shape and arrangement patterns of radular teeth. Because these characters often demonstrate both geographic variation and phenotypic plasticity, morphological characters used to separate species are difficult to score and identifications require expert malacologists [30]. DNA barcoding has effectively identified snail species in other settings [31][32][33][34], therefore we decided to test its effectiveness on Bithyniidae.
The present study is the first to explore the application of DNA barcoding in species identification in the family Bithyniidae. We analyzed variation of the COI barcode region within 10 species/subspecies of Bithyniidae using pairwise sequence comparisons. We then examined the effectiveness of DNA barcoding in differentiating among these species.

Snail collections and preparation
Adult snails of the family Bithyniidae (superfamily Rissoacea) were collected with wire-mesh scoops or by hand in 2009 and 2010 from four regions of Thailand: north, northeast, south, and central ( Figure 1, Table 1). These regions were selected based on results from previous studies [26,28,35]. Each collection site was recorded and its GPS coordinates were determined using a Garmin ® nuvi 203 (Garmin (Asia) Co.,Taiwan). The specimens for this study were collected mostly from public water reservoirs where no permits were required. Owners of the private localities (a rice paddy and a waterfall) were asked for their permission. The owners gave their verbal consent for samples to be collected. All species of those snails are not endangered or protected. The snails were sorted and identified following the protocols in Brandt [26], Chitramvong [36], and Upatham et al. [37]. In addition two subspecies (B. s. siamensis and B. s. goniomphalos) were categorized by geographic distribution.
Each snail was subsequently examined for trematode infections by testing for cercarial shedding twice within a week. Prior to cercarial shedding, the snails were cleaned with dechlorinated tap-water. Shedding was induced under 25 W electric light bulbs for 2 hours at room temperature during the day. For species that shed cercaria at night, black covers were used to achieve total darkness and snails were allowed to shed overnight. Uninfected snails were soaked in phosphate buffered saline (PBS) containing antibiotics (200 unit/ml of penicillin and 100 µg/ml of streptomycin) for 3 to 4 hours before extraction of DNA to ensure that bacterial contamination was minimized.
Each snail was dissected to remove its soft body parts, and kept at -20 °C until further analysis. Each specimen was labeled, databased and imaged. All specimen records are in the project 'JUT-Mitochondrial DNA barcodes identification for snail in family Bithyniidae in Thailand' on BOLD, the Barcode of Life Data Systems [38].

DNA extraction
Total genomic DNA was extracted from whole snail tissue using methods similar to those in Winnepenninckx et al. [39]. Snail tissue was first homogenized in lysis buffer (2% w/v Cetyltri-ammonium bromide; CTAB, 1.4 M NaCl, 0.2% v/v βmercaptoethanol, 20 mM EDTA, 100 mM TrisHCl pH 8, 0.2 mg/ml proteinase K), and then incubated at 55 °C for 6 hours. Subsequently, proteins were precipitated using phenol/ chloroform (1:1) once, followed by phenol/ chloroform/ isoamylalcohol (25:24:1), centrifuged at 13,000 g for 10 min (4°C ) twice, and finally washed with chloroform (1:1). The upper aqueous layer was removed, and DNA was precipitated in isopropanol (2:3 v/v), mixed gently by inverting the tube a few times, put on ice for 15 min, and then spun in a microcentrifuge at 13,000 g for 5 min. After centrifugation, the supernatant was discarded; the DNA pellets were washed in 75% absolute ethanol, and centrifuged at 13,000 g for 5 min. After air-drying, the DNA pellet was re-suspended in TE buffer (10 mM Tris, 1mM EDTA, pH 8.0) and stored at -20 °C until analysis. The DNA concentration and purity were estimated by spectrophotometer (NanoVue, GE Healthcare UK limited,

Amplification and sequencing
PCR protocols followed those used by the Canadian Centre for DNA Barcoding [40], with slight modifications. The PCR reaction was performed on a GeneAmp® PCR System 9700 Thermo Cycler (Applied Biosystem, Foster City, CA). The partial mitochondrial COI gene was amplified using the primers shown in Table 2 [41,42] in a total reaction volume of 50 µl. The amplification reaction consisted of 10xPCR buffer for 5 µl, 10 mM dNTP for 0.25 µl, 50 mM MgCl 2 for 2.5 µl, forward primer for 0.5 µl, reverse primer for 0.5 µl, Platinum Taq polymerase for 0. 24 [43][44][45]. PCR products were visualized on a 1.5% agarose gel and the specific band was cut and its DNA purified and then sequenced in the Biochemistry Department, Faculty of Medicine, Khon Kaen University; Pacific Science Co. LTD (Bangkok, Thailand) and at the Biodiversity Institute of Ontario, Canada.

Data analysis
Forward and reverse DNA sequences were assembled, and edited using Chromas version 2.23 [46], BioEdit v. 5.0.6 [47] and CodonCode v.3.01 (CodonCode Corporation, Dedham, MA). Alignment and homology analysis were performed using CLUSTAL X v. 1.8 [48] and MEGA 4 [49] with pairwise nucleotide sequence divergences calculated using the Kimura 2-parameter (K2P) model [50]. Base composition and distance summaries were obtained using the tools provided on the BOLD workbench (www.boldsystems.org) [38], but only sequences ≥ 350 bp were included in the analysis. A neighbour-joining (NJ) tree was also created using BOLD to provide a preliminary display of the sequence divergences.

Results and Discussion
Ten species/subspecies of Bithyniidae were collected from sites across Thailand (Figure 1 and Figure 2). A total of 217 individuals of these species/subspecies were analyzed for COI, and Neotricula aperta gamma strain (family Hydrobiidae, superfamily Rissoacea) from GenBank (Accession: AF AF188222.1 GI: 11493624 and AF188220.1 GI: 11493620) was used as outgroup. All 217 specimens were identified using morphological characteristics of the adult shells, radular patterns, geographic distribution [35][36][37], and confirmed by a malacologist. From 1-6 individuals of each species/subspecies from each of the five regions were analyzed, as shown in the neighbour-joining tree ( Figure S1). The sequences, and trace files, are available on BOLD (project: JUT).
The pairwise sequence divergences were different among species/subspecies (Table S1). Intraspecific K2P distances averaged 2.3±0.001% (range 0-9.2 %), 4-fold less than the mean congeneric sequence divergence of 8.7±0.002% (range 0-22.2 %). The highest mean intraspecific sequence divergence for an individual species was 4.93±0.22% (range 0-9.2%) for Wattebledia crosseana reflecting the fact that members of this species fell into two distinct sequence clusters ( Table 3). The mean sequence divergence across the family was also high, averaging 17.1% (range 13.0-21.3%).The distributions of intraspecific and interspecific divergences showed limited overlap (Figure 3), because most (65.4%) intraspecific sequences showed less than 2% divergence while 83.4% of the interspecific sequences possessed more than 3%  divergence. As a result, sequence divergences for these snails are similar to those in previous barcoding reports on other organisms [2,12]. Hebert et al. [12] reported that COI sequence divergences among animal species from interspecific COI divergences within the phylum Mollusca averaged 11.1±5.1%. The high intraspecific divergences in W. crosseana and G. wykoffi could indicate the presence of previously unrecognized cryptic species. DNA barcoding has proven invaluable at detecting cryptic species, which in many cases, are subsequently corroborated by life history, morphological or other character sets [51][52][53][54]. For these two snail species, the clusters represent allopatric populations with no apparent morphological differences, so it is currently unclear if they represent merely isolated populations or separate entities with differences yet to be revealed. Conversely, the sharing of identical barcode sequence in G. pygmaea and one northern Thailand population of G. wykoffi may be indicative of introgressive hybridization, incomplete lineage sorting, misidentification, or a previously unrecognized synonymy. Further investigations into these groups are necessary to untangle and confirm these predictions and the use of more holistic approaches to delimit species boundaries will be beneficial.
An important finding in the present study is that the three first intermediate hosts (B. s. siamensis, B. s. goniomphalos and B. funiculata) of Southeast Asian liver fluke can all be distinguished by COI barcodes. All three taxa of Bithynia sp. form monophyletic clusters, with 1.5% divergence between the two subspecies of B. siamensis and both subspecies had 7.1% divergence from B. funiculata (Table 3). Because the two subspecies of B. siamensis are morphologically indistinguishable, the capacity of DNA barcoding to discriminate them is significant. Moreover, morphological similarity has created taxonomic confusion and difficulties in the accurate identification of B. s. siamensis and B. s. goniomphalos which are currently believed to be distributed in the north, central, south and northeast of Thailand [26,29,[36][37][38]. As well, the capacity to rapidly diagnose all stages of the host's life cycle is essential for better understanding of the epidemiology of this parasite-induced disease.
The barcoding success for the Bithyniid species examined in this project was 80%, with nearly all taxa forming discrete monophyletic clusters (Figure 4). The two exceptions are G. pygmaea and one population of G. wykoffi, which share an identical COI sequence (see above). These two taxa might possibly be cryptic species. However, the adult size of G. wykoffi is double that of G. pygmaea. Distinct paraphyly was found in W. crosseana. The results indicated that W. crosseana samples from different localities may well represent cryptic species because they are morphologically similar but genetically distinct. Cryptic species of W. crosseana might be resulted from some factors such as different localities which would develop to different genotype. G. wykoffi was separated into more than one geographically-restricted cluster respectively comprising collection localities from the central, north or northeast regions of Thailand. These clusters might be cryptic species according to this analysis as same as W. crosseana. However, more comprehensive analyses of the systematics of these taxa using more specimens, representing their known geographic distribution, as well as more evidence from independent biological investigations, are required before this hypothesis can be verified. Similar studies which have also been reported in other organisms [52][53][54][55][56][57][58][59], yet over all DNA barcoding has proven reliable in identifying species in more than 90% of the organisms investigated [60]. The neighbour-joining tree and ME analysis also revealed that in general, individuals tended to cluster in accordance with collection localities (Supporting Information, Figure S1, S2). The results from ME analysis were very similar to the neighbor-joining analysis so the latter was used to generate diagrams.
The genera Wattebledia and Bithynia formed monophyletic clusters as well, but Gabbia did not. The selection of Neotricula aperta gamma strain (in the same superfamily) from GenBank as the outgroup appeared legitimate as it clustered separately from other snails in family Bithyniidae. Increased taxon, geographic, and gene sampling would be worthwhile to further explore the two 'barcode outliers' and the ability of COI to infer geographic provenance and phylogenetic affinities in this group.
In summary, the present study has studied genetic-variation in ten species/subspecies of Bithyniidae from Thailand using COI. Sequence divergences were lower for intraspecific than congeneric comparison. Using COI, 80% of the studied snail taxa could accurately identified. In comparison with other methods for identifying snails in this family, DNA barcoding is quicker, easier and more applicable, it is suitable for young snail identification which will be beneficial for understanding the epidemiology of opisthorchiasis transmission.