Application of a DNA barcode using the 16S rRNA gene to diagnose pest Arion species in the USA

Arion Férussac, 1819, is a large, predominantly western Palearctic genus of slugs with a number of species that are of considerable agricultural and forestry concern. Many Arion species have been introduced throughout the world through commerce, and have established themselves as pests in predominantly synanthropic and other disturbed environments in the USA and Canada (Robinson, 1999). Since 1999, nationwide surveys under the umbrella of the USDA North American Slug Project have been conducted to document the current distributions of non-native slug species. Arion species are as yet found to be absent from Puerto Rico, and surveys in Alaska and the Hawaiian Islands are incomplete. In view of the potential threat of new introductions of non-established Arion species, all specimens intercepted from incoming shipments from abroad must be identified by quarantine authorities in order to determine whether or not the associated commodity can be allowed to enter the country. Arion slugs are frequently intercepted at the US ports on commercial imports from Europe, particularly associated with shipments of mushrooms, cut flowers, flower bulbs, horticultural plants and vegetables. However, most intercepted slugs are immature, lacking the diagnostic genitalia that would facilitate their identification. Without internal examination, several species are difficult to identify, as the plethora of synonyms in the malacological literature would suggest. The adults of some species resemble the immatures of others, and the external appearance of a species may also differ according to its geographic origin. The existence of cryptic species complexes that can only be identified by molecular means further complicates the situation. For example, the Arion subfuscus/fuscus complex has been only recently characterized using genetic methods (Pinceel et al., 2004; Pinceel, Jordaens & Backeljau, 2005a). In addition, recent analysis of the subgenus Carinarion indicates that it may include cryptic species or it may represent a species complex (Geenen, Jordaens & Backeljau, 2006). Similarly, within Arion s. s. is what we are calling the ater-rufus-vulgaris complex (ARVC), comprising Arion ater (Linné, 1758), Arion rufus (Linné, 1758) and Arion vulgaris Moquin-Tandon, 1855. The latter species is a serious invasive agricultural and forestry pest throughout most of Europe and is of particular concern to the US quarantine authorities as it could be devastating if introduced into North America. Besides being part of the ARVC, Arion vulgaris by itself has a confused taxonomic status as it has long been known in the European malacological literature as Arion lusitanicus Mabille, 1868 (Falkner, Ripken & Falkner, 2002; Anderson, 2005; Quinteiro et al., 2005; Bank, Falkner & von Proschwitz, 2007). Some workers indicate that Mabille's taxon should be restricted to an Arion (Mesarion) species that is endemic to the Iberian peninsula, and that the widespread synanthropic pest should be called A. vulgaris as this name is the first that unequivocally applies to it (Anderson, 2005; Quinteiro et al., 2005). For the purposes of this research note, we follow the latter arrangement with the caveat that the taxonomic status of the slug remains unresolved. To further complicate the situation, putative hybrids within ARVC species have recently been documented (Schander & von Proschwitz, 2003; Hagnell, von Proschwitz & Schander, 2006). There is increasing interest among taxonomists in using molecular methods in species identification. The application of DNA technology can be used to circumvent the problems associated with identifying a juvenile, damaged specimen or cryptic species (Avise, 2004; Jordaens et al., 2006). An increasingly used technology for taxonomic purposes of a variety of taxa is DNA bar-coding (Hebert et al., 2003). In this methodology, a gene of interest is sequenced from multiple representatives of each described species and the DNA sequences are used to generate a reference database. When a query specimen needs identification, it is sequenced and compared with the reference database using calculations of genetic distances, haplotypic data or topology visualized on phylogenetic trees (De Salle, Egan & Siddall 2005; Savolainen et al., 2005). Unlike phylogenetic studies that test taxonomic classifications and species limits, molecular diagnostic studies characterize the genetic variation in recognized taxa to identify molecular markers that can function as surrogates for morphological characters. Previous studies have demonstrated that the 16S rRNA gene is variable within cryptic Arion species (Pinceel et al., 2004, 2005a, b, c) and useful for identification of snails (Helicidae) and amphibians (Steinke et al., 2005; Vences et al. 2005). In order to assess the utility of the 16S gene for Arion diagnostics, we generated a reference database with representatives of 23 species (with 18 species represented by >1 specimen; Table 1). Other Arion species were not included because specimens or sequences were not available. Taxa based on nomenclature and genetic clusters reported in the NJ tree. Sample sizes (n), bootstrap percentage support (btsp) and interior branch test (IBT) values, the maximum number of differences between two specimens of the same taxon (Max D.) and number of gaps (i.e. total number of nucleotide sites within all INDELs) for this comparison, the minimum number of differences (and percentage) separating two specimens from different taxa, and the documented presence in the U.S. mainland are reported. Distance values for within group comparisons are not applicable (N/A) for taxa represented by a single specimen. aThis genetic type is represented by a single specimen with an A. rufus genital morphology (the A. ater and A. vulgaris morphologies not yet detected in USA). Taxa based on nomenclature and genetic clusters reported in the NJ tree. Sample sizes (n), bootstrap percentage support (btsp) and interior branch test (IBT) values, the maximum number of differences between two specimens of the same taxon (Max D.) and number of gaps (i.e. total number of nucleotide sites within all INDELs) for this comparison, the minimum number of differences (and percentage) separating two specimens from different taxa, and the documented presence in the U.S. mainland are reported. Distance values for within group comparisons are not applicable (N/A) for taxa represented by a single specimen. aThis genetic type is represented by a single specimen with an A. rufus genital morphology (the A. ater and A. vulgaris morphologies not yet detected in USA). For the purposes of this study, we treat Arion (Mesarion) subfuscus (Draparnaud, 1801) and Arion (M.) fuscus (Müller, 1774) as separate species and recognize the population designations ‘S1-5’ for the five A. subfuscus clades (Pinceel et al., 2005a). We recognize the current taxonomic classification of Arion (Carinarion) fasciatus (Nilsson, 1823), Arion (C.) silvaticus Lohmander, 1937 and Arion (C.) circumscriptus Johnston, 1828 based on morphospecies descriptions (Geenen et al., 2006). However, we exclude data from specimens that do not have a species-level identification assigned (i.e. referred to with ‘cf.’ in GenBank). Although these three names might represent subspecies or populations, a 16S gene tree does support three clades that correlate with morphospecies descriptions and a complete revision of Carinarion taxonomy is required (Geenen et al., 2006). We treat Arion anguloi Martín & Gomez, 1988 as a synonym of A. (Mesarion) urbiae De Winter, 1986 based on morphological and genetic evidence (Backeljau et al., 1994). Details of morphological and genetic variation within the species Arion (Mesarion) transsylvanus Simroth, 1885 are currently being studied (Pinceel, 2005). Our reference database includes only 16S sequences for Arion specimens that have been identified to species by experts. In total, 269 sequences are incorporated into our reference database, including: 150 previously published sequences (Pinceel et al., 2004; 2005a, b, c; Geenen et al., 2006), 53 available sequences from GenBank (AY860798-807, AY968073-76, AY947346-83), and 66 new sequences submitted to GenBank (Accession numbers: EU541904-EU541969; specimen vouchers at Academy of Natural Sciences, Philadelphia) from collections made in Europe and the continental USA and identified using adult morphology (Kerney & Cameron, 1979; Likharev & Wiktor, 1980; Castillejo, 1998). DNA was extracted from specimens using either the hot sodium hydroxide and Tris (HotSHOT) extraction method (Truett et al., 2000) or the DNeasy tissue extraction kit of Qiagen (Valenia, CA, USA). According to the HotSHOT protocol, the tissue was placed into a tube with 75 µl of alkaline lysis solution (25 mM NaOH, 0.2 mM EDTA), heated to 95°C for 30 min, cooled to 4°C and then mixed with 75 µl of neutralizing reagent (40 mM Tris–HCl). In the DNeasy extraction protocol, foot tissue was incubated at 55°C for 30–60 min in ATL buffer containing proteinase K (Qiagen). 16S rDNA sequences were generated using variations of primers LR-N-13398 and LR-J-12887: 16S-1 (5′-CGC CTG TTT AAC AAA AAC AT) and 16S-2 (5′-ACG CCG GTT TGA ACT CAG ATC), respectively (Simon et al., 1994). Amplification was performed using Qiagen Taq (with final reaction concentrations of 2 mM of MgCl2, 0.2 mM of each dNTP and 0.2 µM of each primer) and a program file of an initial step of 3 min at 94°C, 32 cycles of 94°C, 20 s/50°C, 20 s/72°C 30 s, and a final step of 5 min at 72°C. PCR product was purified using the QiaQuick purification kit (Qiagen). Sequencing reactions were performed using LARK Technologies Inc. (Houston, TX, USA) or DTCS Quick Start Kit (Beckman-Coulter) on a Beckman-Coulter CEQ8000 sequencer. Sequences (c. 410 bp) were edited with the program Sequencher (v4.6; Gene Codes). Alignments of edited sequences were performed using default ClustalW settings as implemented in MEGA3.1: gap opening and extension penalties of 15 and 6.66, respectively, and DNA weight matrix IUB, transition weight of 0.5 (Kumar, Tamura & Nei, 2004). Data were analysed in MEGA using uncorrected distance estimates (i.e. p-distance) with pair-wise deletion of gaps and Neighbour-Joining (NJ) tree construction. Uncorrected distances were used to generate a conservative estimate of divergence, and the pair-wise deletion option was selected to maximize information between a query and reference sequence. Bootstrap percentage support (btsp) and Interior Branch Test (IBT) values were calculated using 1,000 replicates. Monophyletic groups in the NJ tree are referred to as ‘clusters’ in the current study rather than as ‘clades’ to make a distinction between a DNA barcode NJ tree and evolutionary relationships on a cladogram. Although in an evolutionary context a gap is defined as an independent insertion/deletion (INDEL) event, we count each nucleotide site within a single INDEL event as a gap to emphasize differences in sequences for identification purposes. The maximum variation within clusters (taxa) using substitutions and gaps and the minimum distance between specimens of differing clusters (taxa), excluding gaps, are reported in Table 1. Excluding species represented by a single specimen (which cannot cluster with conspecifics), our NJ tree (not shown) produced topological groups for most of the predefined taxonomic units. Most of these branches had >95% btsp and IBT values (Table 1). A. (M.) subfuscus populations (S1-5) do not form a monophyletic cluster, but rather form two distinct sub-clusters, viz. an S1 + S4 + S5 cluster (86% btsp; 94%, IBT) and an S2 + S3 cluster (89% btsp; 94% IBT). The ARVC specimens form a genetic cluster (89% btsp; 99% IBT) that is comprised of four sub-clusters: designated as A–D (tree not shown). The B sub-cluster is genetically diverse and comprises three main genetic types: B1 (n = 6, 73% bootstrap), B2 (n = 3, 99% bootstrap) and B3 (n = 1). On the basis of the current methods of species identification using characters of the reproductive system, the four sub-clusters do not clearly correspond to the three putative morphological species. The lack of congruence between morphology and genetic type could be the result of incorrect taxonomy, hybridization and/or incomplete sorting of shared ancestral polymorphism. The three morphospecies are represented in the following sub-clusters: A. ater (sub-cluster A), A. rufus (sub-clusters A, B and D) and A. vulgaris (sub-clusters A and C). Our reference database was used to diagnose 57 Arion specimens collected from disturbed habitats in the continental USA and 41 specimens intercepted at ports of entry (data available from N. Barr). This database was developed to determine if the variation captured in using our 16S barcode is sufficient to make identifications, which is distinct from assessing the correctness (i.e. false positives and negatives) of our identifications. Only 28% of the intercepted specimens and 62% of the US-collected specimens were suitable for morphological confirmation (based on maturity and specimen preservation). The morphological identifications for these specimens did not contradict the genetic results. Identifications were made using the 269-sequence database and a modified consensus database in which variable sites (and gaps) within a predefined taxon were converted into missing data (Ns) to create a single representative sequence per taxon. Because high variation within ARVC-B and A. (M.) flagellus Collinge, 1893 (Table 1) generated consensus sequences that contained over 10% missing data, additional consensus sequences were generated using sub-clusters of these groups for use in identification (i.e. consensus ARVC-B forms B1, B2 and B3 and A. flagellus forms A and B). In addition, databases were tested using 0% (i.e. perfect match), 1% and 2% cut-off thresholds, including and excluding gaps (Table 2). Based on variation among taxa, a threshold >2% is considered too high between a query and a reference sequence (Table 1). Identification using the topology of an NJ tree is often not straightforward, and application of a threshold value for distance estimates can be misleading if selected arbitrarily (Meier et al., 2006). However, from a quarantine perspective, evaluation of a genetic distance with set rules for taxonomic assignment (such as thresholds that are based on the available data) is a practical approach to identification that allows reproducibility and consistency in a test. Percentage of taxonomic identifications assigned using various barcode criteria (i.e. inclusion of gaps as differences and percentage cut-off values) and databases. Percentage of taxonomic identifications assigned using various barcode criteria (i.e. inclusion of gaps as differences and percentage cut-off values) and databases. For comparison of threshold values calculated excluding gaps, percentage divergences for the 269-sequence reference database were calculated using the following formula: (no. of mismatches + no. of gaps)(no. of bp aligned − no. of Ns in query)−1. Similar calculations are not applicable to the consensus database because it already incorporates gaps as ‘Ns’ based on intraspecific variation. We found that INDEL events greater than a single-nucleotide long were rare (<1%) for our intraspecific comparisons. However, if a single large gap (e.g. six nucleotides) existed between sequences in a query search using the above formula, this method would indicate each site in the INDEL as a difference (e.g. six differences). Gaps in DNA alignments are more commonly found in ribosomal genes than in protein-encoding genes (e.g. COI) because of differences in molecular evolution and technical difficulties associated with sequencing stretches of nucleotide repeats. Therefore, in addition to accommodating for incomplete sampling of intraspecific variation, a low threshold can be a pragmatic method to compensate for sequencing and editing errors. It should be noted, however, that gaps can be diagnostically informative characters. Although a COI barcode was initially considered for this study, preliminary work revealed that 16S is a more reliable marker based on PCR success rates and that both genes have comparable levels of variation (N. Barr, unpubl.). All of the US specimens were identifiable using a 1% threshold (43 identical matches with a sequence in the 269-sequence database and 14 differed by ≤3 substitutions/gaps). This includes the first confirmed identification of A. (M.) fuscus and type S2 of A. (M.) subfuscus in the USA (Morphological analyses of our A. fuscus specimens from Canada, U.S. and Europe confirm the gonad characterization indicated by Pinceel et al., 2004.) The ARVC sub-clusters such as A, B and D are represented in the USA, based on our query and reference samples. The consensus and 269-sequence reference databases produced nearly identical results (Table 2). Of the intercepted specimens, three failed to provide sufficient DNA for analysis (7% failure rate) and one specimen sequenced reliably only in one direction. Sequence data were generated for 38 of the 41 slugs (Table 2). Sixteen of the interceptions were found to be identical matches to reference sequences, nine had ≤4 substitutions (no gaps), eight had ≤4 deletions (no mismatches), two were similar to Arion (Kobeltia) intermedius Normand, 1852 (five mismatches, no gaps), and three were different from species in the database (>40 substitutions, >10% divergence). However, all three divergent sequences cluster with Arion specimens in a tree when other gastropod groups were included (N. Barr, unpubl.). One of these divergent sequences was most similar to A. flagellus (47 substitutions; >80% btsp, 99% IBT). Based on the reference database A. flagellus contains considerable variation, so the query (intercepted from Spain) could be an undocumented form of a diverse species. According to the reference specimens, diversity within A. flagellus correlates with geographic separation: the specimen from Spain (AY947361) differed from the UK specimens (Table 1) by 38 substitutions. Therefore, it is also possible that these slugs represent several separate lineages. Comparison of the consensus sequence for this taxon and the 269-sequence reference databases (excluding gaps) produced similar results (Table 2). Although it did not affect identification rates, the consensus database inflated the number of mismatches (by one) for several of the A. (M.) subfuscus (S1) queries because gaps treated as Ns affect the subsequent alignment with the query. In terms of database performance, only analyses including gaps and a strict 0% cut-off produced notably lower identification rates. Although highly divergent sequences (≥10%) possess both substitution and gap differences, our data do not suggest any correlation between the presence of substitutions and gaps when the overall sequence divergence is ≤5% (contingency analysis of interceptions: X2 = 1.93, df = Therefore, of gaps in the analysis should have a when a 1% threshold value or consensus database approach is On the basis of and our genetic identifications, several Arion pests genetic are not in the US It should be noted, however, that of evidence is not evidence of work this diagnostic method is to the of these pests within the USA, as as the utility of this method for of these these data already provide useful information for purposes. For example, A. (M.) flagellus has been reported in the US mainland and should be as Similarly, only two and of A. (M.) subfuscus have been detected in the USA, that introductions of other could to greater pest considerable genetic variation within the USA for the ARVC, sub-cluster has not been reported. Since this genetic form correlates with the A. vulgaris this method can be used to provide data in the of this We specimens or identifications to this von Proschwitz, A. Schander, and We also Backeljau, Jordaens and two for on of the and for in specimens for