Molecular Approach to Identifying Three Closely Related Slug Species of the genus Deroceras (Gastropoda: Eupulmonata: Agriolimacidae)

© 2020 Academia Sinica, Taiwan Open Access Molecular Approach to Identifying Three Closely Related Slug Species of the genus Deroceras (Gastropoda: Eupulmonata: Agriolimacidae) Kamila S. Zając1,* and Daniel Stec2 1Institute of Environmental Sciences, Faculty of Biology, Jagiellonian University, Gronostajowa 7, 30-387 Kraków, Poland. *Correspondence: [email protected] (Zając) 2Institute of Zoology and Biomedical Research, Faculty of Biology, Jagiellonian University, Gronostajowa 9, 30-387 Kraków, Poland. E-mail: [email protected] (Stec) Received 27 July 2020 / Accepted 6 September 2020 / Published 12 November 2020 Communicated by Benny K.K. Chan Some species of slugs belonging to the genus Deroceras are invasive and cause severe agricultural damage. Despite extensive knowledge about their invasiveness, data on the molecular differentiation of these morphologically similar species are lacking. Here we present a molecular approach to identifying three closely related species of the genus Deroceras—D. agreste (L., 1758), D. reticulatum (O. F. Müller, 1774) and D. turcicum (Simroth, 1894) (Gastropoda: Eupulmonata: Agriolimacidae)—based on sequences of multiple molecular markers: cytochrome c oxidase subunit I (COI), cytochrome b (cyt-b), internal transcribed spacer 2 (ITS-2) and 28S ribosomal RNA (28S rRNA). We also provide detailed photomicrographs of the penis and penial gland of the three species, as it is the latter that holds the most important phenotypic characters for distinguishing between these taxa. Since identification of the studied species based solely on morphology is considered challenging, contributing a means of molecular differentiation will aid further ecological and biodiversity surveys of these important pests. Key words: Taxonomy, Deroceras, Barcoding, Gastropods, Slugs. BACKGROUND Species identification is central to estimating biodiversity. Traditional taxonomy is based mainly on the morphology of the investigated organism. Since the recognition of cryptic diversity (i.e., complexes of morphologically similar species), which is a consequence of applying molecular techniques to the problem of species identification, morphology-based delineation became insufficient to describe biological diversity (Dayrat 2005). The modern methods for estimating species richness such as DNA-barcoding, metabarcoding and DNA-based taxonomy are now the standard and have been used in meiofauna, ants, decapod crustaceans, lizards and millipedes (Fontaneto et al. 2015; Troncoso-Palacios et al. 2018; Hosoishi and Ogata 2019; Shih et al. 2019; Nguyen et al. 2019). However, all of these methods suffer from limited correspondence between the libraries of DNA sequences and the correctly identified species. Among terrestrial gastropods, many species of slug are especially challenging to identify, as some taxa can only be distinguished by minute morphological differences in their reproductive organs, and the anatomical sections, necessary for proper identification, are often problematic for non-specialists. Deroceras Rafinesque, 1820 is a type genus within the family Agriolimacidae and comprises more than 120 species (Wiktor 2000). Some of the Deroceras species are invasive and drastically affect biodiversity and agriculture by reducing native species richness and crop infestation (Barker 2002). Among Deroceras Citation: Zając KS, Stec D. 2020. Molecular approach to identifying three closely related slug species of the genus Deroceras (Gastropoda: Eupulmonata: Agriolimacidae). Zool Stud 59:55. doi:10.6620/ZS.2020.59-55. Zoological Studies 59:55 (2020) doi:10.6620/ZS.2020.59-55 1 © 2020 Academia Sinica, Taiwan species, D. reticulatum (O. F. Müller, 1774) is a major agricultural pest and is known to cause considerable damage to crops of winter rape, winter wheat and many others (Kozłowski 2010). The majority of studies are concerned with the invasiveness and significance of D. reticulatum in agriculture (Ferguson et al. 1988; Birkett et al. 2004; Berman et al. 2011), whereas the molecular data supporting its pattern of distribution are lacking. Additionally, sequences of two closely related and morphologically very similar species, Deroceras turcicum (Simroth, 1894) (locus typicus: Lake Ohrid in Macedonia) and Deroceras agreste (L., 1758) (terra typica: Sweden) are unknown and uncertain, respectively; what is more, these species have overlapping distribution ranges. The COI sequences of D. agreste are known from a single study about gastropods of Britain and Ireland (Rowson et al. 2014). In that study, the authors found individuals with the genital morphology characteristic of D. reticulatum but mtDNA sequences characteristic of D. agreste. These results were explained to be an example of a one-way introgression of D. reticulatum genes into a D. agreste population (Rowson et al. 2014). Notably, the population of D. agreste used in the study was not from the species’ terra typica, thus the validity of the proposed conclusions remains questionable. Deroceras reticulatum (O. F. Müller, 1774) is a widely distributed synanthropic species occurring in Europe (Wiktor 2004). The origin of this species is unknown, but it likely comes from Central Europe (Kozłowski 2010). It was described from Frideriksdal in Denmark (locus typicus) in 1774. Recently, it was introduced to North and South America, Tasmania, New Zealand and Central Asia (e.g., Crowley and Pain 1977; Hausdorf 2002; Tulli et al. 2009; Welter-Schultes 2012). Notably, the majority of these records have not been verified genetically. For instance, Hausdorf (2002) studied terrestrial gastropods (both slugs and snails) in Colombia and determined slug specimens anatomically. Crowley and Pain (1977) published a study about molluscs of Saint Helena in the last century, identifying them based on classical taxonomic methods. Tulli et al. (2009) conducted research on predation on D. reticulatum by the carabid Scarites athracinus and collected slugs in Argentina. This study did not include any information about slug identification, thus it seems probable that species identification was based on morphological characters and/or external appearance. Finally, a book by Welther-Schultes (2012) summarizes data on the occurrence and distribution of particular species in Europe. Deroceras reticulatum is regarded by the DAISIE organization (Delivering Alien Invasive Species Inventories for Europe) as one of the most invasive species in Europe. It mostly occurs in cultivated areas, open habitats, meadows, gardens, and cemeteries—all in places where human activity is observed (Wiktor 2004; Welter-Schultes 2012). Deroceras reticulatum has a creamy or light coffee cream body pigmentation; it is sometimes blackish and spotted, in which case the dark spots behind the mantle form a reticulate pattern. However, Wiktor (1989) noted considerable variability in the distribution and pigmentation of spots and stated that the most common D. reticulatum individuals are spotted, even though juveniles may be uniform in color. All this makes species identification based on the external appearance extremely difficult. The second species examined in our study, D. agreste, was originally described from Sweden. This species likely occurs throughout Europe, but its precise distribution is unknown (Wiktor 2000; Welter-Schultes 2012). There are currently no data relating to where D. agreste has been introduced and where it is indigenous. These slugs prefer open habitats, especially gardens, meadows, and less often farmland. It can be described as synanthropic and can be found together with, and in similar habitats to, D. reticulatum and D. sturanyi (Wiktor 2000). Its body is yellowish white to light greyish yellow, transparent, without spots (WelterSchultes 2012). Our third species, D. turcicum, occurs in the Balkan region (Wiktor 1989). It was originally described from Lake Ohrid in Macedonia in 1894. Deroceras turcicum has expanded its distribution, and is currently present in Austria, Hungary, Slovenia, Croatia, Bosnia and Herzegovina, Serbia, Montenegro and probably Albania, Greece, and Turkey (Wiktor 2000). Reise and Hutchinson (2001) reported D. turcicum in the Czech Republic and Slovakia, whereas Reise et al. (2005) recorded new localities of this species in Poland. Deroceras turcicum inhabits mostly mixed or beech woodland, parks, and gardens. Its body is dirty-cream or light cream, with dark spots covering the mantle; overall, it is extremely similar in appearance to D. reticulatum (Welter-Schultes 2012). Most reports about the occurrence of these slugs are based on observations and/or species identification based on traditional taxonomic methods such as anatomical sections. In the majority of reports, species identification was not supported by any molecular analysis, e.g., by comparisons of COI sequences (Szybiak 2004; Reise et al. 2005; Sulikowska-Drozd 2007; Hutchinson and Reise 2015). Thus, cited records should be treated with caution, as the species under consideration are morphologically very similar. In this study, we molecularly characterise and compare three closely related Deroceras species (D. agreste, D. reticulatum, D. turcicum), providing for page 2 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan each of them DNA sequences of four molecular markers (cytochrome c oxidase subunit I (COI), cytochrome b – (cyt-b), internal transcribed spacer 2 (ITS-2), 28S ribosomal RNA (28S rRNA)). We also provide comparative genetic analyses of p-distances to show differences between the studied species. Moreover, we visualise phylogenetic relationships between Deroceras species based on available COI sequences for the named species deposited in the GenBank database as well as the sequences obtained in this study. Finally, we also provide a set of photomicrographs that detail the morphology of the penis and penial gland in the studied species. Our contribution will be helpful to other researchers using modern approaches to identify the difficult to distinguish Deroceras species. Our findings may also have significance in biological conservation, as these slugs have been shown to have considerable economic importance in agriculture. MATERIALS AND METHODS Material collection and species identification All specimens of Deroceras agreste, D. turcicum and D. reticulatum were obtained from a museum collection deposited in the Museum of Natural History, Wrocław University, Poland (Table 1). Gastropods were preserved in 70–75% ethyl alcohol (Wiktor and Jurkowska 2007). All specimens were identified and labelled by Professor Andrzej Wiktor (Wrocław University, Poland) and further verified by us based on anatomical characteristics under Nikon SMZ1500 stereomicroscope (Table 2). From each slug, a small piece of tissue, preserved in 99.8% ethanol, was taken for DNA extraction. The photos of penial glands were taken using a digital camera (Nikon DSFil-U2) associated with the stereomicroscope and later assembled into figures in Corel Photo-Paint X6, ver. 16.4.1.1281. Schematic drawings of the genitalia of each species were made using InkScape software based on Wiktor’s (2000) monograph. DNA extraction, amplification and sequencing DNA was extracted from the collected tissue using Sherlock AX kit (A&A Biotechnology, Poland) according to the manufacturer’s protocol. This kit is intended for material with trace amounts of DNA. To obtain DNA sequences of the COI, ITS-2, 28S rRNA, cyt-b genes, PCR reactions were run using the primers listed in table 3. For COI, ITS-2 and 28S rRNA specific primers for Deroceras species were designed de novo based on sequences of Deroceras (D. laeve, D. reticulatum, D. invadens), Arion (Arion sp., A. silvaticus, A. vulgaris, A. distinctus, A. rufus) and Limax (L. maximus, L. flavus, L. cinereoniger) species available in GenBank (for exact accession numbers, see Table S1). A PCR cocktail and profiles/conditions for specific markers are given in the supplementary materials (Tables S2 and S3). To check the DNA quality, 3 μl sample of PCR product was run on a 1.5% agarose gel for 30 min at 100 V. PCR products were cleaned using NucleoSpin Gel and PCR Clean-up (MachereyNagel, Germany). A sequencing reaction was performed using a 10-μl reaction mixture consisting of 2 μl of PCR Table 1. Collection data on Deroceras specimens used in this study (N: number of sequenced individuals from each locality) Species N Locality Collector Collection date Fig. 1. GenBank Accession numbers COI cyt-b ITS-2 28S rRNA D. reticulatum 3 Aragón, Pirineos Mts, Diazas near Torla (N of Huesca), Spain A. Wiktor 07.07.1987 A MN934414 MN931235 MT361823 MN930516 3Leon province, Villaturel and Villasabariego (NE of Leon), Spain A. Wiktor 19.06.1987 B, C D. agreste 2 Västergötland, Sweden H. W. Walden 30.09.1964 D MN934413 MN931234 MT361821 MN930515 2 Bohuslän, Sweden H. W. Walden 10.10.1965 E 2 Cantabria, Spain A. Wiktor 12.11.1988 F D. turcicum 2 near Maçka in Trabzon, Turkey A. Riedel 07.12.1985 G MN934415 MN931236 MT361822 MN930517 2Zafanos, Trabzon, Turkey A. Riedel 01.12.1985 H 2 near Maçka in Trabzon, Turkey A. Wiktor 29.04.2001 I page 3 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan Fig. 1. Genitalia of D. reticulatum (A: Aragón, Pirineos Mts, Diazas near Torla (N of Huesca), Spain; B–C: Leon province, Villaturel and Villasabariego (NE of Leon), Spain), D. agreste (D: Västergötland, Sweden; E: Bohuslän, Sweden; F: Cantabria, Spain) and D. turcicum (G: near Maçka in Trabzon, Turkey; H: Zafanos, Trabzon, Turkey; I: near Maçka in Trabzon, Turkey) detailing penis and penial gland morphology. Penial gland and constriction in the penis are indicated by flat and indented arrowheads respectively. Scale bars in mm. Table 2. Morphological differences in reproductive system between D. reticulatum, D. agreste and D. turcicum according to Welter-Schultes (2012) D. reticulatum D. agreste D. turcicum penis fleshy and with a silky sheen, in the shape of an irregular sac, divided into 2 parts by a deep lateral constriction oval in juveniles, with lateral constriction in adults strong swelling on the anterior part of the penis penial gland variable shape; usually a few branches or a single long branch finger-shaped, never branched, smooth, without glandular papillae variable or reduced stimulator large, conical and narrow cone-shaped, with narrow base conical rectal caecum large caecum on rectum well-developed, clearly longer than wider well-developed vas deferens opens into the penis wall facing the external body side opens laterally close to penial gland opens at the base of the penial gland page 4 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan product, 0.15 μl of primer, 1 μl of sequencing buffer (BrilliantDye Terminator Sequencing Kit, Nimagen, The Netherlands), 5.85 μl of ddH2O and 1 μl of Terminator (BrilliantDye Terminator Sequencing Kit, Nimagen, The Netherlands). The sequencing program consisted of four steps: 1 min initial denaturation at 96°C, followed by 10 s denaturation at 96°C, 5 s annealing at 55°C, 4 min elongation at 60°C for 25 cycles. Sequencing products were cleaned using ExTerminator (A&A Biotechnology, Poland) and sequenced in both directions in Genomed company (Warsaw, Poland). Obtained sequences of single haplotypes per marker and species were deposited in GenBank. Data analysis All obtained COI sequences were blasted by NCBI BLAST (Altschul et al. 1990) to verify species identification and find similar homologous sequences. In order to visualize phylogenetic relationships between the studied species, selected COI sequences (the longest and with published status) of six named Deroceras species—D. golcheri (Accession numbers: JN248291293), D. laeve (MG421043, MG423214, KX959499501, HM584699), D. panormitanum (JN248304313), D. reticulatum (MG421618, MG421373, MG421157, MG421125, MG421099, MF545181, MF545161, MF545125, MF545107, FJ481179), D. agreste (KF894312, KF894346-247, KF894375) and D. invadens (JN248295-297, JN248301-303, JN248314315, KX959490)—were downloaded from GenBank (Eskelson et al. 2011; Reise et al. 2011; Rowson et al. 2014; Araiza-Gómez et al. 2017) and aligned with sequences obtained in this study for D. reticulatum, D. agreste and D. turcicum in BioEdit 5.0.0 (Hall 1999) with the multiple alignment function of ClustalW (Thompson et al. 1994). Aligned sequences were trimmed and translated into protein sequences in MEGA v. 7 to check against pseudogenes. Uncorrected pairwise distances between species were calculated in MEGA7 (Kumar et al. 2016). Uncorrected genetic distances for COI were calculated based on all sequences obtained for D. agreste, D. reticulatum and D. turcicum together with those used for phylogenetic tree calculation (alignment length: 402 bp). For this molecular marker, genetic distances were calculated in two ways, first between all sequences in the dataset (see Table S4), and second between species defined as groups (Table 4). In the case of cyt-b, ITS-2 and 28S rRNA, only the newly obtained sequences for D. agreste, D. reticulatum and D. turcicum were used to calculate genetic distances (lengths of alignments: 299 bp, 280 bp and 352 bp, respectively) (Table 5). The concatenated dataset comprises 79 sequences in total (COI sequences of six Deroceras species from GenBank (42 sequences), COI, cyt-b, ITS-2 and 28S rRNA sequences for studied species obtained in this study (36 sequences) and one sequence of Limax maximus (Accession Number: KM612139) as the outgroup). The obtained alignment of the concatenated sequences (1619 bp) were divided into eight data blocks constituting the ITS-2, 28S rRNA as well as COI and cytochrome b that were separated into three codon positions using PartitionFinder v. 2.1.1 under the Bayesian Information Criterion (BIC) (Lanfear et al. 2016). PartitionFinder retained five out of eight predefined partitions, the best fit of which were (1) F81 for 28S rRNA, (2) TIM+G for the 2nd codon position in COI and 3rd codon position in cytochrome b, (3) TRN+G for the 3rd codon position in COI and 1st codon position in cytochrome b, (4) TRN+I for the 1st codon position in COI and 2nd codon position in cytochrome b, (5) JC for ITS-2. For the maximum-likelihood (ML) analysis in RAxML, models GTR, GTR+I, GTR+G and GTR+I+G were also tested (Stamatakis 2014). The best fit model for all partitions in this analysis was GTR+G. ML topology was constructed using RAxML v. 8.0.19 (Stamatakis 2014). The tree branches were supported by bootstrap analysis with 1,000 replicates. Bootstrap support values ≥ 70% were regarded as significant statistical support. Bayesian inference (BI) marginal posterior probabilities were calculated in MrBayes v. 3.2 (Huelsenbeck and Ronquist 2001, Table 3. Primers used for amplification of the four DNA fragments sequenced in the study Fragment Primer name Primer sequence (5’-3’) References COI COI_Der_Ff TATATATAATTTTTGGGGTTTGATGTGG This study COI_Der_Rr CAAAAAGATGTTGATATAAAATAGG CytB 151-F TGTGGRGCNACYGTWATYACTAA Merritt et al. 1998 270-R AANAGGAARTAYCAYTCNGGYTG ITS-2 ITS2_Der_Ff GTCGGCTAGTCWAAAGCAATCG This study ITS2_Der_Rr CCGCTTCACTCGCCGTTACT 28S rRNA 28S_Der_Ff GCTAAATACTTGCACGAGTCCG This study 28S_Der_Rr ACGGTTGCCCAGTCTCTCC page 5 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan Huelsenbeck et al. 2001) with one cold and three heated Markov chains for 10 million generations and trees were sampled every 1000 generations. In the BI consensus tree, clades recovered with posterior probability (PP) between 0.95 and 1 were considered well supported, those with PP between 0.90 and 0.94 were considered moderately supported and those with lower PP were considered unsupported. Obtained trees were visualized in FigTree v.1.4.3, available at: http://tree.bio.ed.ac.uk/ software/figtree. RESULTS All specimens previously identified by prof. Andrzej Wiktor were confirmed as Deroceras reticulatum, D. agreste and D. turcicum based on morphology of their genitalia, especially penis and penial gland appearance (Figs. 1A–I and 2A–C). Deroceras reticulatum has a pronounced laterally constricted penis with both parts inscribed in an oval and with a large penial gland on the posterior end of the penis (Fig. 1A–C). The penial gland can vary in shape (Fig. 1A–C). It can be represented by a single process with both lateral sides covered by easily identifiable papillae-like structures (Fig. 1A) or made up of a few branches covered by glandular papillae and merged into a common short trunk (Fig. 1B). However, these papillae-like structures or glandular papillae on the lateral sides of the penial gland can sometimes be only faintly visible (Fig. 1C). Deroceras agreste is characterized by sack-shaped penis, which can be laterally narrowed in older specimens. The penial gland does not show variability in its shape, it is visible as a small conical nodule or unbranched finger with a smooth surface without any papillae-like structures on the lateral sides (Fig. 1D–F). Deroceras turcicum has a laterally constricted penis with the anterior part heavily bloated and a narrower posterior part. The penial gland can be a complex, multifurcated structure or may be significantly reduced (Fig. 1G–I). More specifically, penial glands can consist of multiple knots, with the lateral sides covered by convexity/constrictions that can resemble papilla-like structures (but usually bigger than in D. reticulatum) (Fig. 1G–H), or it can be reduced and comprise only a few, usually short nodular processes set directly on the penis or on a short common stem (Fig. 1I). Good quality sequences for all four molecular markers were obtained from all analyzed specimens. Only single haplotypes were recovered for each marker in each species. For the three analyzed species, the Table 4. Genetic distances (%) between Deroceras species calculated from COI sequences 1234567 1D. agreste 2D. golcheri 17.1 3D. invadens 14.1 9.8 4D. laeve 14.4 10.3 7.8 5D. panormitanum 15.7 6.0 7.4 8.5 6D. reticulatum 5.6 15.1 14.7 13.3 13.6 7D. turcicum 9.8 16.6 14.2 14.4 15.2 8.6 Table 5. Genetic distances (%) among three Deroceras species calculated from cytochrome b, ITS-2 and 28S rRNA sequences 123 cytochrome b1D. agreste - 2D. reticulatum 6.8 - 3D. turcicum 10.8 12.1 - ITS-2 1 D. agreste - 2D. reticulatum 7.6 - 3D. turcicum 8.0 3.5 - 28S rRNA 1 D. agreste - 2D. reticulatum 0.6 - 3D. turcicum 1.2 1.7 - page 6 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan following sequences were obtained: • D. agreste: the COI sequence (GenBank: MN934413), 512 bp long, the cytochrome b sequence (GenBank: MN931234), 390 bp long, the ITS-2 sequence (GenBank: MT361821), 280 bp long, the 28S rRNA sequence (GenBank: MN930515), 400 bp long. • D. reticulatum: the COI sequence (GenBank: MN934414), 651 bp long, the cytochrome b sequence (GenBank: MN931235), 299 bp long, the ITS-2 sequence (GenBank: MT361823), 352 bp long, the 28S rRNA sequence (GenBank: MN930516), 352 bp long. • D. turcicum: the COI sequence (GenBank: MN934415), 420 bp long, the cytochrome b sequence (GenBank: MN931236), 299 bp long, the ITS-2 sequence (GenBank: MT361822), 350 bp long, the 28S rRNA sequence (GenBank: MN930517), 409 bp long. Regarding the D. reticulatum COI sequence, the verification in BLAST resulted in fitting to D. reticulatum (KY765589 (Ahn et al. 2017); Query cover: 100%; E-value: 0.00; Perc. ident: 98.92%). The best hits for the D. agreste sequence were D. reticulatum and D. agreste (LS974196, (Kropf, unpublished) and KF894312 (Rowson et al. 2014), respectively; for both species: Query cover: 100%; E-value: 0.00; Perc. ident: 99.41%). For the D. turcicum sequence, the best hit was D. reticulatum (MG421618 (Dewaard, unpublished); Query cover: 99%; E-value: 1e-162; Perc. ident: 91.87%). Genetic distances calculated for COI sequences between Deroceras species ranged from 5.6% (between D. agreste and D. reticulatum) to 17.1% (between D. agreste and D. golcheri) (12.01% on average) (Table 4). Regarding other molecular markers, genetic distances for cytochrome b range from 6.8% (between D. agreste and D. reticulatum) to 12.1% (between D. reticulatum and D. turcicum); for ITS-2, they range from 3.5% (between D. reticulatum and D. turcicum) to 8.0% (between D. agreste and D. turcicum); and for 28S rRNA, genetic distances range from 0.6% (between D. agreste and D. reticulatum) to 1.7% (between D. reticulatum and D. turcicum) (Table 5). Maximum-Likelihood (ML) and Bayesian (BI) phylogenetic trees resulted in similar tree topologies (Fig. 3). Specimens belonging to particular species constituted seven separate clades that, in most cases, were strongly supported in BI analysis and were moderately or not supported by ML analysis (Fig. 3). The three species that we focused on in this study (D. reticulatum, D. agreste and D. turcicum) have been found to be distinct but closely related to one another and together form a clade that is sister to the other analyzed species (Fig. 3). However, due to polytomy, specific relationships between these three species could not be inferred. Individuals identified based on anatomical studies conducted by Wiktor and us (GenBank: MN934413, this study) and by Rowson et al. (2014) (GenBank: KF894312, KF894346-247, KF894375) as D. agreste clustered together, thus tentatively providing an additional and independent confirmation of species affiliations (Fig. 3). The remaining four Deroceras species formed another big clade, sister to (D. reticulatum + D. agreste + D. turcicum). Here, similarly, the specific relationships between these species could not be fully recovered due to the several weakly supported nodes that indicate polytomy (Fig. 3). DISCUSSION In this study, we showed that the three Fig. 2. Schematic drawings of genitalia of D. reticulatum (A), D. agreste (B) and D. turcicum (C) after Wiktor (2000). Bc, bursa copulatrix; Pl, penial gland; P, penis; At, atrium. Scale bars = 1 mm. (B) (A) (C) At At At PP P Pl Pl Pl Bc Bc Bc page 7 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan morphologically similar species, Deroceras reticulatum, D. agreste, D. turcicum, are genetically easily differentiable. The results from the BLAST search clearly underline the problems with identifying these taxa, namely that D. agreste and D. turcicum are often wrongly designated as the more common D. reticulatum. The DNA sequences presented in our study are valuable as they came from specimens deposited in a museum collection and were identified by Professor Andrzej Wiktor, a specialist in slug taxonomy, and verified in this study. Additionally, sequences for D. agreste were obtained from the specimens collected from terra typica in Sweden. Our phylogenetic results, presented above, further confirmed the molecular distinctives and monophyly of each studied species. As the number of traits associated with the external appearance of slugs is usually extremely limited, the morphology of the reproductive system is considered crucial for precisely identifying many slug species (e.g., Reise and Hutchinson 2001). In accordance with Wiktor (2000), we identified penis and penial gland appearance to be the most useful character for at least a provisional, and not so complicated identification (even for nonspecialists), due to the ease by which these structures can be located. In concordance with previous studies, we also confirmed a considerable morphological variability in the penial glad in D. reticulatum (Wiktor 2000) and D. turcicum (Wiktor 2000; Reise and Hutchinson 2001), as well as almost no variation in D. agreste (Wiktor 2000). However, in dubious cases, the application of genetic data and cross-verification with the specific DNA sequences provided in our study will be of great help to other researchers. Moreover, this may have increased significance when only juveniles are present in a given sampling area, since anatomical studies of their genitalia will not be possible. The mitochondrial COI and cytochrome b were the most variable molecular markers and thus most Fig. 3. A phylogenetic tree constructed from concatenated sequences of COI, cyt-b, ITS-2 and 28S rRNA of Deroceras species. Only topology of the Bayesian tree is shown. Bolded and underlined GenBank accession numbers indicate sequences obtained in this study. Numbers at the nodes indicate Bayesian posterior probability and bootstrap values separated by the “/” symbol. Values less than 0.90 for BI and less than 70 for ML are not shown. The scale bar represents 0.02 substitutions per nucleotide position. KM612139Limax maximus 0.02 FJ481179 MF545107 MF545125 MF545161 MF545181 MG421099 MG421125 MG421157 MG421373 MG421618 MN934414+ MN931235 MT361823+ MN930516 0.78 / 86 1.00 / 84 KF894247 KF894312 KF894346 MN934413+ MN931234 MT361821+ MN930515 KF894375 0.96 0.94 / 89 D. reticulatum D. agreste MN934415+ MN931236 MT361822+ MN930517 1.00 / 100 0.90 JN248295 JN248296 JN248297 JN248301 JN248302 JN248303 1.00 / 100 JN248314 KX959490 JN248315 D. invadens 1.00 / 98 JN248304 JN248305 JN248311 JN248307 JN248308 JN248309 JN248310 JN248306 JN248312 JN248313 0.91 / 77 1.00 / 97 JN248291 JN248292 JN248293 1.00 / 100 HM584699 KX959499 KX959500 KX959501 MG423214 MG421043 1.00 D. turcicum D. golcheri D. panormitanum D. laeve page 8 of 11Zoological Studies 59:55 (2020) © 2020 Academia Sinica, Taiwan informative in terms of species distinctions. The first of these is already known to be the most useful molecular marker for species identification (Fontaneto et al. 2015); however, as we show, cytochrome b is also variable enough to distinguish Deroceras species. Moreover, cytochrome b was successfully used for identifying species and deliminating members of other gastropod families (e.g., Merritt et al. 1998; Groenenberg et al. 2016). Notably, amplification and/or sequencing of these two DNA fragments are often problematic because of their length, so obtaining good quality sequences is challenging but possible in most cases with specific primers designed for short DNA fragments. Moreover, specimens deposited in museum collections are usually old and thus their DNA is usually degraded and highly fragmented. Work with such material often requires specific kits for DNA extraction that are intended for material with trace DNA content. Therefore, by presenting reliable sequences for cyt-b and COI, we increased the chances for correct species identifications in future research. Regarding nuclear markers, ITS-2 and 28S rRNA, it has already been shown that they are more conservative than the mentioned mitochondrial cyt-b and COI (Wade et al. 2006; Zając et al. 2020). This makes them especially suitable for resolving deeper phylogenetic nodes, especially those between genera, and of higher taxonomic levels (e.g., Wade et al. 2006; Neiber and Hausdorf 2015). Similar to other studies, our results show that ITS-2 is more variable than 28S rRNA. As such, ITS-2 can be used as additional support for identification by mitochondrial markers as well as in phylogeographic research (Schniebs et al. 2013; Zhou et al. 2017). Our phylogenetic analysis revealed that the tree topology for Deroceras species is similar to that found in the recent study conducted by Rawson et al. (2014). The addition of new sequences of D. turcicum enabled us to show a close relationship with the morphologically similar D. reticulatum and D. agreste. Even though the phylogenetic position of the genus Deroceras is known (Bouchet and Rocroi 2005), the specific and detailed relationships between its species have only been examined superficially (e.g., Koene and Schulenburg 2005). As mentioned above, molecular data regarding these slugs are extremely limited. Currently, the existing data enabled the testing of relationships among only seven named Deroceras taxa, from a genus comprising more than 120 species (Wiktor 2000). Importantly, many of them (including specimens of the type series) are deposited in museum collections, thereby providing an opportunity to obtain tissue for genetic studies from accurately identified specimens. This might be important, as morphological species determination is challenging and it may be difficult to obtain fresh material from the field because of the species’ status (e.g., endangered, endemic) or its rare occurrence. Thus, museum collections constitute a valuable source of reliable barcodes which might help enhance species identification and facilitate biodiversity estimations (Neiber and Hausdorf 2016; Neiber et al. 2017). CONCLUSIONS Our study presents a molecular approach to identify three closely related species of Deroceras slugs by providing, for each of them, definitive DNA sequences of four molecular markers. We designed de novo specific primers for three molecular markers (COI, ITS-2, 28S rRNA), which successfully amplified DNA fragments, and which can be used on the other Deroceras species held within museum collections. The DNA sequences of D. turcicum are presented in our study for the first time, providing new and compelling evidence of its close relationship with the two morphologically similar taxa (D. reticulatum and D. agreste). The sequences provided in this study are reliable barcodes which may aid in correct species identification, even by non-specialists. This study underlines the important role of natural history museums and identifies their collections as priceless resources of global biodiversity. Finally, our results have application in studies concerning the distribution, conservation biology and introgression of Deroceras species. Acknowledgments: We thank Jolanta Jurkowska and the Museum of Natural History, Wrocław University (Poland) for making the collection of Prof. Andrzej Wiktor available for us. We are also grateful to Brian Blagden (Scottish Environment Protection Agency, UK), Krzysztof Miler (Institute of Systematics and Evolution of Animals, Polish Academy of Sciences) and Dorota Lachowska-Cierlik (Institute of Zoology and Biomedical Research, Jagiellonian University) and the anonymous reviewer for their critical reading of the manuscript and comments, which improved our paper. This work was supported by a grant from the Jagiellonian University (DS/MND/WB/INOŚ/12/2018). Authors’ contributions: KSZ designed the study, analyzed the data and wrote the manuscript. DS designed the primers and contributed to the data analysis and manuscript writing. Competing interests: KSZ and DS declare that they have no conflicts of interest. page 9 of 11Zoological Studies 59:55 (2020)