Molecular and morphological identification of Biomphalaria species from the state of São Paulo, Brazil

Abstract DNA barcoding and morphological characters were used to identify adult snails belonging to the genus Biomphalaria from 17 municipalities in the state of São Paulo, Brazil. The DNA barcode analysis also included twenty-nine sequences retrieved from GenBank. The final data set of 104 sequences of the mitochondrial cytochrome oxidase I (COI) gene was analyzed for K2P intraspecific and interspecific divergences, through tree-reconstruction methods (Neighbor-Joining, Maximum Likelihood and Bayesian inference), and by applying different models (ABGD, bPTP, GMYC) to partition the sequences according to the pattern of genetic variation. Twenty-seven morphological parameters of internal organs were used to identify specimens. The molecular taxonomy of Biomphalaria agreed with the morphological identification of specimens from the same collection locality. DNA barcoding may therefore be a useful supporting tool for identifying Biomphalaria snails in areas at risk for schistosomiasis.

Identification of Biomphalaria specimens to the species level and analysis of infection by S. mansoni are key elements of surveillance strategies for schistosomiasis control and elimination (PAHO 1968, WHO 2013).Shell morphology is of limited use for identifying different species of snails in this genus (Paraense 1966;Jarne and Théron 2001), and therefore the anatomical characters described by Paraense (1961Paraense ( , 1975Paraense ( , 2001) ) are used instead.However, identification of Biomphalaria solely based on morphological characters is constrained by phenotypic plasticity, the limited descriptions of cryptic species, and the difficulty in applying species-diagnostic characters to juvenile specimens (Carvalho et al. 2008;Teodoro et al. 2010).The issue of how useful molecular tools may be in the identification of Biomphalaria snails has become particularly important in recent years as there is consensus among malacologists that morphological identification using internal anatomical parameters is susceptible to error, especially when the snails being analyzed belong to complexes of morphologically similar species (Paraense 1972(Paraense , 1974(Paraense , 1988;;Spatz et al. 1999;Vidigal et al. 2000).To overcome these limitations and difficulties associated with traditional taxonomy, various methodologies based on molecular markers have been developed.
When used in conjunction with bioinformatics tools and sequence databases, DNA barcoding routinely facilitates the identification of biological species (Ratnasingham and Hebert 2007;Casiraghi et al. 2010).This technique is based on the polymorphism of a short region (approximately 600 bp long) of the mitochondrial cytochrome c oxidase 1 (COI) gene (Hebert et al. 2003).DNA barcode includes a series of strategies for delimiting species into molecular operational taxonomic units (MOTUs) using a combination of laboratory and bioinformatics methods (Fontaneto et al. 2013).The most important strategies for identifying MOTUs include analysis of intraspecific and interspecific genetic distances, and analyses based on population and phylogenetic models.These approaches include (ABGD) (Puillandre et al. 2012) and the barcode index number (BIN) system (Ratnasingham and Hebert 2013), which use algorithms based on the partition of molecular data according distance methods, and the generalized mixed Yule coalescent (GMYC) method (Fujisawa and Barraclough 2013) and Bayesian Poisson Tree Processes (bPTP) method (Zhang et al. 2013).
DNA barcoding has been used to augment morphological identification of Bulinus in Africa (Kane at al. 2008;Stothard et al. 2013;Standley et al. 2014), and yielded better results than identifications based on shell characters.Although there are over 500 COI sequences in GenBank from snails of the genus Biomphalaria found in African and Neotropical regions, most DNA barcoding studies use African species.There is therefore a dearth of knowledge about the effectiveness of DNA barcoding in taxonomic identification of Neotropical species of Biomphalaria (Standley et al. 2011;Tuan et al. 2012).
Here, we investigate the utility of analysis of distributions of intraspecific and interspecific COI divergences based on genetic distances, tree reconstruction methods based on Bayesian inference, Maximum Likelihood (ML), and K2P-Neighbor-Joining (NJ) grouping of sequences, and the ABGD, GMYC and bPTP methods for delimitation of Biomphalaria species in conjunction with schistosomiasis field surveys.

Experimental design
Planorbids were collected in 17 municipalities in the state of São Paulo, Brazil between May 2012 and January 2013 (Fig. 1).The collection points were georeferenced with a Garmin ETrex Summit® GPS (Table 1).
Samples were collected from freshwater habitats in the Paranapanema, Tietê, Ribeira do Iguape and Paraíba do Sul River basins and the northern coast of São Paulo that had been previously surveyed and classified according to the risk for schistosomiasis transmission as part of a program to monitor snails that are intermediate hosts of S. mansoni (Biomphalaria).
In accordance with the methods described in the Brazilian Ministry of Health Schistosomiasis Surveillance and Control Program (Ministry of Health 2008), snails were collected at sampling stations in each freshwater body and grouped into batches according to their origin.Most of the snails in each batch were then exposed to artificial light in the laboratory to determine whether they were infected with cercariae.At least two specimens from each batch were used for morphological analysis and at least two for the DNA barcode analysis.
DNA barcoding was applied to 75 adult snails taken from samples collected in the field.Only snails that did not have any parasite larvae in their digestive gland and ovotestis were used for molecular identification.Shells were removed by compress-  1).
ing each snail between two slides.After removing the shell fragments, each crushed snail was transferred to a clean Petri dish.The portion of the cephalopodal mass corresponding to the foot was excised under a stereo microscope with forceps and scissors and used as starting material for isolation of total DNA.To maximize the efficiency of genomic DNA purification we used fresh material that had not been fixed.Each specimen was then dissected and identified to the species level based on the presence or absence of the renal ridge and the most informative characters of the male and female copulatory organs.DNA barcoding was carried out in a blind fashion, i.e., without prior knowledge of the general morphological characteristics identified in the animal.
An additional 118 adult specimens were taken from the same field samples (at least two per batch) and scored for 27 morphological characters used by Paraense (1975Paraense ( , 1981Paraense ( , 1984Paraense ( , 2001) ) in his descriptions of Neotropical species of the genus Biomphalaria.The soft parts were removed from the shell after placing the snails in 70°C for 40 seconds and then fixing them in Railliet-Henry's solution (distilled water 930 mL, sodium chloride 6 g, formalin 50 mL and glacial acetic acid 20 mL).After at least 24 hours of fixation, the specimens were dissected following Deslandes' (1951) protocols to examined the renal tube and reproductive system.Specimens were not anesthetized  (Mulvey and Bandoni 1994).
prior to fixation to ensure that the procedure followed was the same as that used in our malacology laboratories.
The longitudinal renal ridge is considered the gold-standard character for differentiating B. glabrata (Paraense and Deslandes 1959) from other species in the genus, in which the ridge is absent.The anterior and posterior regions of the vagina were examined.The proportions for the diameter and length of the oviduct were based on the nidamental gland; for the diameter of the uterus, the cephalic portion of the nidamental gland was used, for the length of the uterus, the posterior region of the vagina; for the length of the spermathecal duct, the body of the spermatheca; and for the length of the anterior region of the vagina, the posterior region of the vagina.The relative proportions of the organs or structures were used for comparisons together with the shell and mantle pigmentation pattern.

DNA extraction, amplification and sequencing
DNA isolation was carried out with the DNeasy Tissue Kit (Qiagen®).A fragment of the COI gene (~600 bp) was amplified with the LCO/HCO primers (Folmer et al. 1994).Polymerase chain reaction (PCR) was carried out in a total volume of 50 µL and the following reaction mixture: 10-100 ng of DNA, 0.2 mM of each dNTP, 0.10 µM of each primer and 1 U of Taq DNA polymerase in the supplied reaction buffer.The cycling conditions consisted of an initial 3 min step at 95°C for denaturation; 25 cycles of 1 min at 95°C, 1 min at 47°C and 1 min 30 s at 72°C and a final extension step of 7 min at 72°C (Tuan et al. 2012).PCR products were purified with a Qiagen purification kit and then sequenced in the Biotechnology Center at the Butantan Institute in an ABI3100 automated sequencer (Applied Biosystems®).

Molecular data analysis
The electropherograms obtained from forward and reverse sequencing of each specimen were corrected using CHROMAS (Technelysium Pty Ltd.) and then aligned with CLUSTALX version 1.8 (Thompson et al. 1997).The aligned sequences were edited with BIOEDIT version 7.0 (Hall 1999), and the general polymorphism of the sequences was calculated in DNAsp version 5 (Librado and Rozas 2009).
The final alignment consisted of a matrix of 75 COI sequences from the collected specimens (36 B. tenagophila, 12 B. occidentalis, 10 B. glabrata, 9 B. straminea, 1 B. intermedia, 7 B. peregrina) and 29 COI sequences of Biomphalaria from other Neotropical areas that were retrieved from GenBank (Table 1).
Intraspecific and interspecific genetic distances (Kimura 1980) were calculated by pairwise comparison of the sequences of all the individuals using the Kimura 2-parameter (K2P) method with the MEGA 6 (Molecular Evolutionary Genetics Analysis) package (Tamura et al. 2013).Three tree-based methods were performed for phyloge-netic reconstructions.The K2P distance matrix was used to reconstruct a Neighbor-Joining (NJ) tree.MEGA 6 was also used to perform Maximum Likelihood analysis.In the ML analysis, the GTR+I+G model of sequence evolution was chosen using the Akaike information criterion as implemented in MODELTEST 2.3 (Nylander 2004).The reliability of NJ and ML topologies was evaluated using bootstrap support with 1000 replicates.The parameters estimated by MODELTEST were also used in a Bayesian Markov-Chain Monte Carlo (MCMC) analysis in MRBAYES 3.1 (Huelsenbeck and Ronquist 2001;Ronquist and Huelsenbeck 2003).Two simultaneous independent searches were run for 1.5 x 106 generations, with trees saved every 100 generations, and the first 1.500 sampled trees of each search discarded as "burn-in".

Morphological analysis
The morphological identifications of the 118 adult snails that were studied are presented in Table 2.The results of morphological analysis revealed the following: Shell: the presence of a carina, the shape of the whorls and the shape of the shell aperture distinguished B. tenagophila and B. occidentalis from the other species in the group.Renal tube: The presence of renal ridge was observed in all the B. glabrata specimens studied.Pigmentation of the mantle: adult specimens of B. tenagophila, B. glabrata and B. occidentalis had more uniform pigmentation than the four other species studied, which had blotchy pigmentation.Reproductive system: the presence of a vaginal pouch in B. tenagophila and its absence in B. occidentalis differentiates these two species.Biomphalaria straminea and B. intermedia had marked variation in the posterior region of the vagina; in the former, the corrugation in this region was markedly wavy while in the latter it was swollen.
Biomphalaria peregrina differed from the species in the B. straminea complex (B.straminea and B. intermedia) in the width of the oviduct, the length of the uterus, the length of the spermathecal duct and the length of the anterior region of the vagina.Biomphalaria intermedia differed from B. straminea in the number of ovotestis diverticula, the length of the oviduct, the presence of an oviduct pouch, the number of prostate diverticula and the width of the uterus.Biomphalaria oligoza, B. peregrina, B. intermedia and B. straminea are differentiated by the number and shape of the ovotestis diverticula, appearence and size of seminal vesicle, the number and shape of the prostate diverticula, and the shape of the prepuce.The 27 morphological characters used to identify Biomphalaria are detailed in Table 2.

Molecular analysis
The final alignment matrix for the 104 sequences consisted of 549 characters including 25% polymorphic, 21% parsimony-informative and 12 unique sites (Table 3).
The K2P sequence divergence for intraspecific comparisons ranged from 0.0% to 4.0%, while for interspecific comparisons the corresponding figure varied from 4.0% to 12% (Table 4).The greatest intraspecific genetic distances were observed between specimens of B. peregrina from SP and Rio Grande do Sul (southern Brazil) (4.0%) and specimens of B. glabrata from Rio Grande do Sul and Puerto Rico (3.9%).
The frequency distribution of the 104 analyzed sequences indicates that although there were some extreme pairwise distances (>3%) in B. glabrata, B. tenagophila, B. peregrina and B. straminea; intraspecific and interspecific divergences did not overlap (Fig. 2A).Nevertheless, a typical barcode gap was not observed in this dataset.A closer inspection of the distances for each taxonomic group shows that there is a clear barcode gap between B. glabrata, B. straminea, B. peregrina and B. intermedia.There was no clear barcode gap between closely related B. tenagophila, B. t. guaibensis and B. occidentalis (interspecific distance 3-4%) (Fig. 2 C, D, E, F).
When run using the default settings, ABGD recovered five different subunits of B. glabrata.This result may be explained by the pronounced genetic variation in this species, but the possibility that these subgroups represent cryptic taxa cannot be ruled out.
The trees generated by the Bayesian, ML and NJ methods (Fig. 3) delineated six well supported groups (posterior probabilities and bootstrap values ≥90) congruent with the current classification of Biomphalaria.The only B. intermedia sequence appeared in a distinct branch supported by low Bayesian and bootstrap values.

Discussion
This study sought to determine the utility of DNA barcoding in delimiting species in freshwater snails of the genus Biomphalaria.The Bayesian, ML and NJ analyses (Fig. 3, Suppl.material 1) yielded trees with well-supported internal branches (≥90), resolving six out of the seven taxa as monophyletic groups.
The assessment of the potential of DNA barcode for species differentiation in Biomphalaria essentially revolves around the comparison of results of the morphological and molecular analysis of closely similar or taxonomically ambiguous species.In the case of the three taxa in the B. tenagophila complex, one character that is normally effective for specific identification is the vaginal pouch, which is present in B. tenagophila and B. t. guaibensis but not in B. occidentalis.(The anatomical features of these three taxa were illustrated by Tuan et al. 2012).Although we did not observe this in our material, in some specimens of B. occidentalis there is a slight projection of the ventral wall of the vagina (Paraense 1981), which raises questions regarding the distinctness of this taxon.
The intraspecific genetic distance within B. tenagophila showed values with a range from 0 to 3% (Table 4, Figs 2, 3).A high level of genetic divergence within this species was obtained for sequences associated with specimens collected in Juquiá (650,651,653), Itariri (535), Embu das Artes (524,535) and São Paulo (549,551).Due to these values we could not assign a clear barcode gap between B. tenagophila and B. occidentalis and B. t. guaibensis (Fig. 2 b, d, f).However, the Bayesian tree inferred from COI data (Fig. 3), as well as the ABGD and both bPTP and GMYC analyses recovered these close related taxa as distinct groups.We suggest that in geographical areas where B. tenagophila species complex have the same geographical distribution.
The application of the 3-4% cutoff value for maximum intraspecific divergence may be appropriate for our dataset as 36% of the intraspecific comparisons reached this value (Table 4).The highest values for intraspecific divergence (>3%) do not appear to be a consequence of geographic distance given that the greatest divergence in B. tenagophila was between closely proximal localities in São Paulo state (Fig. 3).Biomphalaria glabrata and B. tenagophila, are differentiated by the renal ridge, which is present in the former and absent in the latter.Paraense and Deslandes (1959) described a false ridge that runs obliquely to the renal tube and is attached to the pneumostome, in specimens of B. tenagophila from Macaé, RJ.The presence of this false ridge in B. tenagophila may lead to incorrectly identify this species, particularly in juvenile specimens or specimens that have not been properly fixed.Five specimens  of B. tenagophila in our study (three from São Lourenço da Serra and two from São Paulo) had a membrane on the renal tubes similar to that described by Paraense and Deslandes.The genetic distance of 9% between B. glabrata and B. tenagophila observed with both genetic distance and tree-based approaches show that DNA barcoding is an important tool for identifying these closely similar taxa.
The ABGD analysis partitioned B. glabrata into five distinct groups, while the GMYC analysis yielded a more cohesive group.Despite the pronounced COI divergence within B. glabrata, in all the specimens analyzed here the renal ridge has been considered a robust and consistent taxonomical character, suggesting that morphology is more effective than DNA barcode in this case.However, bBTP analysis and phylogenetic reconstruction supported B. glabrata as single and well supported MOTU, a result congruent with the morphological identification.
Another group of morphologically similar and frequently misidentified congeners includes B. intermedia and B. straminea; the latter a natural intermediate host of S. mansoni.Of the seventeen diagnostic characters common to B. straminea and B. intermedia, the degree of corrugation in the dorsal wall of the vagina has been used to these taxa as a species complex (Paraense 1975) .The vaginal corrugation, which is markedly wavy in B. straminea appears as swollen in B. intermedia.The large genetic divergence between B. straminea and B. intermedia, which was 9% greater than the intraspecific values in both species, indicates that these two species can be identified by DNA barcode.Note, however, that our study only included two of the three species in the B. straminea complex, as B. kuhniana does not occur in São Paulo state (Paraense 1988, Teodoro et al. 2010).In addition, we were unable to collect many specimens of B. intermedia owing to its rareness in São Paulo state.
Our findings show Biomphalaria species delimitation by phylogenetic approaches and bPTP yielded the same groups identified by traditional taxonomy.The use of DNA barcode to identify species in conjunction with Biomphalaria surveys requires the application of both evolutionary and bioinformatics criteria, making it a timeconsuming approach that is dependent on specialist knowledge.Morphological identification also requires specialist knowledge.However, as shown in this study DNA barcoding can identify subtle (genetic) differences between intraspecific populations that are not detectable by traditional morphological study.
Furthermore, morphological identification of Biomphalaria species depends on subjective interpretation of anatomical variations, as these are measured in terms of relative rather than absolute sizes.We therefore agree with Hebert and Gregory (2005, p. 853), who stated that by reversing the logic of standard taxonomic approaches that "operate in an a priori fashion-seeking…morphological discontinuities", DNA barcoding may, as "a posteriori approach", direct the study of morphological variation in genetically divergent groups of Biomphalaria.

Figure 2 .
Figure 2. A histogram showing pairwise Kimura 2-parameter intraspecific and interspecific distances for 104 Biomphalaria cytochrome oxidase I sequences B-H pairwise distances between each species and the other taxa analyzed.

Figure 3 .
Figure 3. Bayesian phylogram.Support values for individual branches are given as Bayesian credibility/ ML bootstrap/NJ bootstrap and are depicted above each node.The different shades of gray identify morphological species.The red, green and blue bars indicate species delimitations based on the distance-based (ABGD) and tree-based (bPTP and GMYC) models, respectively.

Table 1 .
Collection localities, sample information, and GenBank accession numbers for COI sequences used in this study.

Sites/ Country Map locality Municipality Latitude (S) Longitude (W) COI sequence GenBank accession number
* M-Line refers to a laboratory strain of Biomphalaria glabrata derived from a Puerto Rico pigmented snail and an albino Brazilian snail

Table 2 .
Morphological characters used to identify 118 Biomphalaria specimens from the state of São Paulo.