DNA barcoding using skin exuviates can improve identification and biodiversity studies of snakes

Abstract Snakes represent a taxonomically underdeveloped group of animals in India with a lack of experts and incomplete taxonomic descriptions being the main deterrents to advances in this area. Molecular taxonomic approaches using DNA barcoding could aid in snake identification as well as studies of biodiversity. Here a non-invasive sampling method using DNA barcoding is tested using skin exuviates. Taxonomically authenticated samples were collected and tested for validation and comparisons to unknown snake exuviate samples. This approach was also used to construct the first comprehensive study targeting the snake species from Maharashtra state in India. A total of 92 skin exuviate samples were collected and tested for this study. Of these, 81 samples were successfully DNA barcoded and compared with unknown samples for assignment of taxonomic identity. Good quality DNA was obtained irrespective of age and quality of the exuviate material, and all unknown samples were successfully identified. A total of 23 species of snakes were identified, six of which were in the list of Endangered species (Red Data Book). Intra- and inter-specific distance values were also calculated, and these were sufficient to allow discrimination among species and between species without ambiguity in most cases. Two samples were suspected to represent cryptic species based on deep K2P divergence values (>3%), and one sample could be identified to the genus level only. Eleven samples failed to amplify COI sequences, suggesting the need for alternative PCR primer pairs. This study clearly documents how snake skin exuviates can be used for DNA barcoding, estimates of diversity and population genetic structuring in a noninvasive manner.


Introduction
The snakes are known to be one of the most successful group of predatory vertebrates. They occupy a wide range of habitats and almost 3500 species have been described (Burbrink & Crother, 2011;Vitt & Caldwell, 2009). However, limited numbers of useful morphological characters, coupled with inadequate sampling and availability of qualified experts have been major curbs impeding advances in understanding snake diversity (Cox & Davis, 2013;Gould et al., 1977;Nee et al., 1994). In addition, morphological plasticity in shape or colour patterns at the local or geographical scale is known to possibly contribute to misidentification (Cox & Davis, 2013). Despite these and other known limitations inherent in reliance on morphological approaches to characterise diversity of such species, possible DNA-based solutions to address these taxonomic and diversity issues have thus far attracted little attention.
Snakes are also a challenging group for biological sampling due to their deadly poisonous bite and aggressive attacking behaviour. The small body size of these species also makes it difficult to draw blood from snakes. If they are injured during handling, snakes often do not survive because they fall prey to predatory ants or sondary infections. As a result, snakes are often killed for acquisition of samples (Eguchi & Eguchi, 2000).
In India, however, all snake species are protected under the Indian Wildlife Protection Act, 1972, and cannot be touched or handled for this type of sampling. Although dead, degraded specimens or trace biological are sometimes used, morphologically based taxonomic identification of this material is challenging. Alternatively, molecular identification techniques have proved to be an effective tool for species identification in such cases (Dubey et al., 2009;Guha, 2005;Hebert et al., 2003;Jerome et al., 2003;Purcell et al., 2004;Teletchea et al., 2005;Wong et al., 2003). Among the available molecular techniques, DNA barcoding is considered to be among the most advanced, reliable and inexpensive method for biological species identification. The barcoding concept typically uses the 648 bp sequence of the Cytochrome Oxidase I (COI) gene (Barrett & Hebert, 2005;Hebert et al., 2003Hebert et al., , 2004aHajibabaei et al., 2006;Smith et al., 2006;Ward et al., 2005). In addition to identification, the intraand inter-species distance values obtained from barcode sequence data can advance our understanding of species evolutionary relationships. Data from DNA studies is easily accumulated and can be developed for populations of any species. Estimates of diversity made at global and/or regional scales can also be critical for decisions on conservation priorities.
Being a highly Vulnerable group for Extinction due to rapid urbanization, population pressure, habitat alteration, poaching to some extent and illegal international trading, our knowledge of snake diversity needs to be updated in this regard. Biological realities specific to snakes make it imperative to have some reliable method to study their diversity and make correct species identification without killing, physical handling, or otherwise disturbing their resting burrows.
In the present study, we have demonstrated that a non-invasive method for biological sampling to obtain DNA from snake skin exuviates can be used to identify species using the COI gene. This method also provides for applications in biodiversity and population studies in the future. The results of this study using DNA barcoding from snake skin exuviates can also be used to complement and corroborate other efforts to identify snake species in India.

Ethical statement
This study does not include any experimentation on live animals. Also, we have not handled or sampled any biological material directly from live animals. Exuviate samples used in this study are considered as animal waste. Therefore no permissions are required from any ethical committees on animal experimentation or conservation authorities in India.

Taxon sampling, DNA isolation, PCR reaction and sequencing
Taxonomically authenticated snake exuviate samples were procured from Aurangabad Municipal Corporation Zoo, Aurangabad (n ¼ 34) and Katraj Snake Park, Pune (n ¼ 21), India. Also, 37 unknown exuviate samples were obtained from public collections (Appendix 1). The exuviate samples used range from recent to very old and partially degraded and difficult to handle due to excessive weathering. Samples were repeatedly washed with saline solution (8%) followed by ethanol (97-100%) to remove traces of any foreign material. Labeled samples and sub samples taken for preservation were stored in absolute alcohol at À20 C until being processed in the laboratory.
DNA was purified from a small section of the snake skin exuviates ($2.0 cm 2 ) using the Promega Wizard Õ nucleic acid purification kit (Cat no. A1620, Promega Corporation, Madison, WI) following manufacturer's instructions. Isolated DNA was quantified using a NanoDrop1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). DNA concentration was usually found to be in the range of 40-60 ng/ml with 260/280 ratios in the range of 1.75-1.85. Quality was also checked on 0.8% agarose gel electrophoresis, diluted (if necessary) to final concentration of 50 ng/ml and stored at À20 C for further use.
PCR amplification and DNA sequencing PCR amplifications were conducted in 25 mL volumes having 1 Â PCR buffer, 200 mM each dNTP, 2.5 mM MgCl 2 , 0.5 mM each primer, 1.00 unit of Taq polymerase (Gold Taq, Invitrogen), and 50 ng tamplet DNA. Primers used for amplification were C_VF1LFt (VF1_t1:VF1d_t1:LepF1_t1:VF1i_t1), C_VR1LRt1 (VR1_t1:VR1d_t1:LepR1_t1:VR1i_t) (Ivanova et al., 2007) having first step: 5 cycles of 94 C for 30 s, 50 C for 40 s, 72 C for 1 min, second step: 35 cycles of 94 C for 30 s, 53 C for 40 s, 72 C for 1 min, 72 C for 10 min. Another primer pair used for unsuccessful samples was RepCOI-F-5 0 -TNTTMTCAACNA ACCACAAAGA-3 0 and RepCOI-R: 5 0 -ACTTCTGGRTGKCCA AARAATCA-3 0 (Nagy et al., 2012), PCR conditions includes initial denaturation step of 3 min at 94 C; 40 cycles of 94 C for 40 s, 49 C for 30 s and 72 C for 60 s followed by final extension step of 7 min at 72 C. Fragments of the mitochondrial genes cytochrome C oxidase subunit I (COI; 652 bp) were amplified and cycle sequenced. Details on the samples applied are available and can be found in the dataset at: dx.doi.org/10.5883/DS-SDBI. Amplified products were visualized on 1.5% agarose gels. Both strands of the gene segments were amplified using a Verity thermal cycler and sequenced on an ABI 3130 automated sequencer in 10 mL volumes using ABI Prism BigDye Terminator chemistry (ABI Biosystems, Grand Island, NY).
Sequences obtained were (in some cases) corrected manually for misreads and merged into single sequence files using BioEdit Version 5.0.1 (Hall, 1999). Data obtained was analyzed using BOLD online tools (www.boldsystems.org). The genetic distances among and within species were calculated based on the Kimura two-parameter (K2P) substitution model (Kimura, 1980). A neighbor-joining (NJ) tree of K2P distances was created to provide a graphic representation of the patterning of divergence between species with the MEGA v5.0 software package (Tamura et al., 2011).

Results
DNA was obtained from all 92 snake exuviate samples (Table 1; Figure 1). It appears that either old or degraded samples and relatively fresh exuviate samples did not differ significantly in terms of quality and quantity of DNA. From these, we obtained a total of 81 high quality mitochondrial COI sequences ($652 bp; BOLD project code ISDB; www.boldsystem.org; DOI dx.doi.org/ 10.5883/DS-SDBI). The contig assemblies generated did not exhibit any evidence of insertions, deletions, or stop codons. Irrespective of DNA quality, eleven samples failed to generate DNA barcode. Five of these were from the set of authenticated specimens and six were from the set of unknown skin exuviates (Table 1). This suggests the need for new and/or additional clade specific primers for these snake species.
For the taxonomic identifications, at least one taxonomically authenticated sample was available for all taxa under this study, except for Amphiesma stolatum. These included 55 taxonomically authenticated snake samples from two Zoos. Also, 35 unknown snake exuviate samples were classified to the species level without ambiguity. We confirmed that these specimens represent a total of twenty three snake species belonging to nineteen genera, five families and one order (DOI dx.doi.org/10.5883/DS-SDBI). Two samples were suspected to represent possible cryptic species based on deep K2P divergence (43%), and one sample could be identified to the genus level only. Fifteen species were represented by more than one specimen (average, 4 specimens/species) in this dataset, while eight species were represented by single individuals only (Table 1). Our sampling includes six snake species that are listed as having threatened status in the Red Data Book of endangered and threatened species.

Distance analysis
All species were discriminated by the DNA barcode sequences forming a cohesive group in NJ tree is presented in Figure 2. The average K2P distance of specimens within species was 0.67% compared with 2.34% within genera and 19.78 within family ( Table 2). The overall variation is in increasing order with higher taxonomic rank, supporting a marked change of genetic divergence at the taxonomic boundaries (Table 2).
No overlap was found between the pairwise intra-and interspecific distance values. The intra-specific distances ranged from 0.0006 to 0.0349 (Table 3), while the inter-specific values ranged from 0.00 to 0.46 (Table 4). The inter-specific values are more than an order of magnitude higher than the intra-specific distance.
It is also appears that the intra-and inter-specific distance are not dependent on the number of specimens analysed.

Sequence divergence
The overall mean GC content was 45.71% whereas the mean GC content at codon positions 1-3 was 50.79%, 42.32% and 43.88% respectively (Supplementary Table S1).
The nearest-neighbor distance analysis (NND), namely the minimum genetic distance between a species and its closest neighbor species, was carried out to analyze the distribution of genetic divergence. The analysis showed only two cases with K2P genetic divergence lower than 2% (Figure 3), but enough to separate these species. The NND analysis also showed that 98% of con-specific comparisons were lower than 2% genetic distance ( Figure 4).

Intra and inter species distances
In this study. to validate the effectiveness of COI barcode data for the nearest related species of snakes we have analyzed two subspecies belonging to the genus Coelognathus helena ( Figure 5). Both of these species, Coelognathus helena helena and Coelognathus helena monticollaris, formed separate clades in the NJ tree ( Figure 5). The mean intra species difference here was 0.51% and the intra genus difference was 2.28%. This is sufficiently high for discriminating these subspecies.
We have also analyzed the sequence divergences for the Argyrogena fasciolata specimens collected from distant localities ($300 kms away). In cases where the samples originated from the same region they were clustered together and formed a cohesive group ( Figure 6) with intra species K2P distances of 1.07% (0.0107 log distance) (Table 3).

Blast analysis
A few COI sequences from the taxa included in this study were available in the NCBI GeneBank (Amphiesma stolatum, Coelognathus helena monticollaris, Cyclophiops major, Daboia russellii, Echis carinatus, Naja naja, Ophiophagus hannah and Ramphotyphlops braminus) or BOLD databases (Naja naja, Ophiophagus Hannah and Ramphotyphlops braminus). Most of the available records are not from India (Supplementary Table  S2).

Suspected cryptic diversity
Zoo managers authenticated the species from the exuviate samples for almost all from the original specimens used in this study. Two samples were authenticated as Gongylophis conicus (NS6) and Ptyas mucosa (ZSS69). These showed deep intraspecific (K2P) divergence (43) when compared to their nearest species, and these two samples were recorded as possible cryptic species (Figure 2).

Discussion
The ability to make correct species identifications is a critical component of any biodiversity estimation as well as for enforcement of conservation protocols. Morphological plasticity and regional scale variations can confound the use of traditional taxonomical procedures for identification of species and requires enormous numbers of specimens. The enormous numbers of specimens needed to use this methodology has also led to the opinion that such animal collections can be responsible for lowering the biodiversity of some vulnerable group of species  (Akbarsha, 2007;Rosse, 1995;Sathyanarayana, 2009Sathyanarayana, , 2013. Globally, there is also great pressure to minimize the handling and use of animals in experimentation. In addition, compliance with ethical procedures in animal experimentation is becoming more and more time consuming. The use of DNA barcoding as an alternative identification and study method signifies a shift from the nearly exclusive dependence on morphological characters to identification of species. This study demonstrates the usefulness of DNA barcoding of snake species identification in a non-invasive way because we have successfully used snake exuviates for obtaining high quality DNA and for amplification of COI gene sequences to obtain the full length (652 bp) DNA barcodes needed for the identification of species.
Exuviate samples work well with squamates in concert with easy-to-use protocols for obtaining good quality DNA (Eguchi &   Eguchi, 2000). One of the main goals of our study was to use exuviate samples to establish a DNA barcode database for snake species from Maharashtra state and to test the utility of these COI sequences for identifying species and assigning them to major units (corresponding to genera and families). In some case the samples are not reliably amplified, suggesting that multiplexing with additional primers or designing clade-specific primers will be necessary for future expansion of the databases. Nagy et al. (2012) made similar observations during his study on reptiles of Madagascar.
Generation of a valid library of sequences is prerequisite for barcode studies (Ekrem et al., 2007;Hickerson et al., 2006;Meyer & Paulay, 2005;Prendini, 2005). For this we obtained morphologically validated snake exuviate samples from zoos in Maharashtra state of India to begin compiling a standard barcode library of snake species. Although some COI barcode data is available through the NCBI as well as BOLD databases as accessed on 17 January 2014 (Supplementary Table S2), the inadequacy of reference snake sequences in these databases has been noted by others (Dubey et al., 2011). The difficulty of procuring enough specimens due to restrictions on handling and       DOI: 10.3109/19401736.2014.905830 collecting of live specimens has no doubt contributed to this paucity of material. This further emphasizes the value of the approach we used which requires only the exuviate material which would otherwise be discarded. All species analyzed in this study processed distinct genetic characteristics. No overlap was found between pairwise intra-and inter-specific distance values, and this also validates the usefulness of a barcoding approach for studies of these taxa (Baker & Palumbi, 1994;Clarke et al., 2006;Fong et al., 2007;Hebert et al., 2004b;Mikkelsen et al., 2007;Ntie et al., 2010). The NJ trees also showed cohesive groupings of nearly al taxa without ambiguity.
Two species studied here (Gongylophis conicus and Ptyas mucosa) are candidates for designation as cryptic species based on the magnitude of the deep K2P divergence values found (Hebert et al., 2004b). Previously, Gongylophis conicus was described as a single species from India by Murray (1884) followed by Boulenger (1890). In case of the genus Ptyas, three species have been described from India: (Ptyas mucosa (Linnaeus, 1758), Ptyas korros (Schlegel, 1837) and Ptyas nigromarginata (Blyth, 1854), two of which have geographic distributions restricted to North Eastern regions of India, Bangaladesh, Myanmar and China (Farid & Parvin, 2001;Boulenger, 1890). Only P. mucosa is reported from the rest of the Indian regions. The suggestion of the possibility of cryptic Pytas species in our study may indicate need to verify the range of distribution of these species, as previously suggested by Farid & Parvin (2001).
DNA barcodes of the different type derived also have potential applications for conservation (Ardura et al., 2010;Francis et al. 2010) and for deterrents of illegal trading (Eaton et al., 2010;Welton et al., 2012). Previously in India, data on snake diversity and taxonomic identification has been impeded for several reasons (Gower & Winkler, 2007). This has resulted in species being included in the Red Data Book for Threatened species even though their diversity status has not been adequately evaluated (Sharma, 2003). The ability to reliably identify snake species using material from all life stages of snake species should be a high priority to set up a properly functioning biodiversity monitoring system. Five of the snake species studied here are Critically Endangered, two are Endangered and three are listed as Threatened (from the appendices of the Convention on the International Trade of Endangered Species) and their commerce needs to be internationally monitored. India also ranks high for illegal trading of snake skin (Fitzgerald, 1989), and such activity would be much easier to detect using Barcode methodology.
In our study it is shown that the use of skin exuviates provides quality DNA that is reliable and easy to obtain (Eguchi & Eguchi, 2000) for the DNA barcoding method. This allows for noninvasive sampling and has significant benefits for animal welfare. The COI database we developed provides barcode data for 81 snakes, and it provides a solid basis for future biodiversity and population studies using molecular identification methods. In addition, this data can be used for detection and regulation of unauthorized trading of live specimens, skins or body parts or products from snakes.

Conclusions
In the absence of whole animals, snake exuviates are useful for species identifications based on the use of COI gene sequences for biodiversity assessments and for the possible detection of illegal activities. This study also contributed new barcode records for 20 snake species to the BOLD database as well as 15 to the NCBI Genebank from India. The suspected cryptic diversity for two snake species identified here may also promote re-evaluation of the status of these taxa. This could be accomplished by in-depth investigations using additional molecular markers as well as morphological characters on new collections of material.

Declaration of interest
Authors of this manuscript declare that they do not have any conflict of interests regarding this work. This study was carried out using funding and facilities from the Paul Centre for DNA Barcoding