DNA barcodes and phylogenetic affinities of the terrestrial slugs Arion gilvus and A. ponsi (Gastropoda, Pulmonata, Arionidae)

Abstract The Iberian Peninsula is a region with a high endemicity of species of the terrestrial slug subgenus Mesarion. Many of these species have been described mainly on subtle differences in their proximal genitalia. It therefore remains to be investigated 1) whether these locally diverged taxa also represent different species under a phylogenetic species concept as has been shown for other Mesarion species outside the Iberian Peninsula, and 2) how these taxa are phylogenetically related. Here, we analysed DNA sequence data of two mitochondrial (COI and 16S) genes, and of the nuclear ITS1 region, to explore the phylogenetic affinities of two of these endemic taxa, viz. Arion gilvus Torres Mínguez, 1925 and A. ponsi Quintana Cardona, 2007. We also evaluated the use of these DNA sequence data as DNA barcodes for both species. Our results showed that ITS did not allow to differentiate among most of the Mesarion molecular operational taxonomic units (MOTUs) / morphospecies in Mesarion. Yet, the overall mean p-distance among the Mesarion MOTUs / morphospecies for both mtDNA fragments (16.7% for COI, 13% for 16S) was comparable to that between A. ponsi and its closest relative A. molinae (COI: 14.2%; 16S: 16.2%) and to that between A. gilvus and its closest relative A. urbiae (COI: 14.4%; 16S: 13.4%). Hence, with respect to mtDNA divergence, both A. ponsi and A. gilvus, behave as other Mesarion species or putative species-level MOTUs and thus are confirmed as distinct ‘species’.

between the oviduct and the pedunculus of the bursa copulatrix (unlike in A. molinae, where the pedunculus is positioned in between the epiphallus and oviduct) (figures 3-5 in Quintana Cardona 2007).
Arion gilvus (Figure 2) was described from 'Mandol' in the Spanish Province of Tarragona. However, the toponym 'Mandol' seems to be erroneous (e.g. Bech 1990) and therefore Castillejo (1990) assigned eight specimens with an A. gilvus morphology from Serra de Pandóls near Gandesa (Province of Tarragona) as topotypes [see also Castillejo and Rodríguez (1991)]. Subsequently, A. gilvus was redescribed by Garrido (1992). Afterwards, the species has also been found in the Provinces of Valencia, Teruel and Albacete [Borredà (1994), figure 15 in Castillejo (1997), figure 1 in Quinteiro et al. (2005)]. Arion gilvus reaches a length of up to 65 mm when extended. It has a yellowish to brown dorsum that gets lighter downwards at the sides and dark lateral bands that have a yellowish grey line on their upper side (Figure 1). The sole is white or evenly yellowish and the mucus is pale yellow (Torres Mínguez 1925, Bech 1990, Garrido 1992, Castillejo 1997. The epiphallus, the pedunculus of the bursa copulatrix, and the free oviduct join the atrium on a single line with the pedunculus of the bursa copulatrix in the middle, as in A. molinae, but in contrast to the latter, the epiphallus is longer than the vas deferens (Torres Mínguez 1925, Borredà 1994, Castillejo 1997, and figures 26-28 in Garrido et al. 1995).  As illustrated by Arion ponsi and A. gilvus, the alleged species-specific genital differences among the Iberian species of the A. subfuscus complex are very subtle and little is known about their intraspecific variation. Moreover, genital differences among arionid taxa do not necessarily imply reproductive isolation (Dreijers et al. 2013). Hence, if alleged species-specific phenotypic differences in arionids are to be interpreted under a phylogenetic species concept, then their correlation with reproductive isolation should be corroborated by molecular data. Molecular markers have been very effective in this respect (e.g. Pinceel et al. 2005a, b, Quinteiro et al. 2005, Geenen et al. 2006, Jordaens et al. 2010). As such, Quinteiro et al. (2005) investigated the taxonomic affinities of Iberian Mesarion species using DNA sequence data. Their analysis of the nuclear ribosomal internal transcribed spacer 1 region (ITS1) showed a polytomy of Mesarion species, yet, the analysis of the mitochondrial NADH dehydrogenase I (ND1) gene suggested a strongly bootstrap supported group of Iberian Mesarion species with a continental-Mediterranean distribution (A. paularensis, A. baeticus, A. urbiae, A. anguloi, A. wiktori, and A. gilvus), and an unsupported group of species with an Atlantic distribution (A. lusitanicus, A. nobrei, A. fuligineus, A. hispanicus and A. flagellus). In addition, the positions of three Pyrenean species (A. lizarrustii, A. iratii, A. molinae) remained unresolved. More specifically, the ND1 data placed A. gilvus as sister taxon of A. urbiae and A. anguloi. Quinteiro et al. (2005) did not study individuals from the Balearic Islands and thus probably did not include A. ponsi.
Because DNA sequence data do not only provide phylogenetic information, but can also serve as DNA barcodes for species identification (Hebert et al. 2003(Hebert et al. , 2004, we here expand on the work of Quinteiro et al. (2005) by 1) characterizing A. gilvus and A. ponsi using mitochondrial COI and 16S rDNA gene fragments, and the larger part of the nuclear ITS1 region, 2) exploring the phylogenetic affinities of A. gilvus and A. ponsi within the subgenus Mesarion, and 3) providing diagnostic COI barcodes for both species.

Material and methods
Information on the species and specimens included here is provided in Table 1. In total, we screened 45 specimens (Table 1). DNA was extracted from small parts of the foot using a NucleoSpin Tissue Kit (Macherey-Nagel, Düren) following the manufacturer's instructions. PCR reactions were done in 25 µl reaction volumes that contained 1.5 mM MgCl 2 in 1 × PCR buffer (Qiagen), 0.2 mM of each dNTP, 0.2 µM of each primer and 0.5 units of Taq polymerase (Qiagen). A fragment of the mitochondrial COI and 16S genes was amplified using primer pairs LCO1490 and HCO2198 (Folmer et al. 1994) and 16Sar and 16Sbr (Palumbi 1996), respectively. The nuclear ITS1 region (except the ± first 30 bp) was amplified using the primer pair ITS1L and 58C (Hillis and Dixon 1991). The PCR profile was an initial denaturation step of 5 min at 95 °C, followed by 35 cycles of 45 s at 95 °C, 45 s at an annealing temperature of 40 °C (COI), 42 °C (16S) or 55 °C (ITS1) and 1.5 min at 72 °C, and ending with a final extension step of 5 min at 72 °C. PCR products were purified using the GFX PCR DNA Purification Kit (GE Healthcare) following the manufacturer's instructions. Purified DNA was diluted in 15 µl of sterile water. PCR-products were bidirectionally sequenced using the ABI PRISM BigDye® Terminator v1.1 Cycle Sequencing Kit and run on a ABI3130xl Genetic Analyzer. Sequences were assembled in SeqScape v2.5 (Life Technologies) and inconsistencies were checked by eye on the chromatogram. Sequences were submitted to GenBank under accession numbers KF305196-KF305225 for COI, KF356212-KF356245 for 16S and KF385449-KF385469 for ITS1. These datasets were supplemented with DNA sequences from GenBank [including a few species of the other Arion subgenera (Table 1)]. We used those of Carinarion as outgroup.
Sequences were aligned in ClustalW (Thompson et al. 1994) with default settings and without subsequent manual adjustments. In each alignment sequences were trimmed table 1. List of specimens used in this study with specimen ID, sampling locality, GenBank accession numbers for the COI, 16S and ITS1 sequences, and collection number at the museums (if available). Neo-, para-and topotypes have been indicated. Specimen codes with an asterisk are data taken from Quinteiro et al. (2005); NA = not assessed. The specimen ID and GenBank accession numbers of newly sequenced specimens are given in bold. to equal length. The final alignments had a length of 504 bp (COI), 408 bp (16S) and 587 bp (ITS1), and of 1499 bp after concatenating the three fragments. The COI sequences were translated to amino acid sequences to check for stop codons (but none were found). The ITS1 sequences were also analysed together with those of Quinteiro et al. (2005). In this way we could extend our taxon coverage to A. hispanicus Simroth, 1886, A. fuligineus Morelet, 1845 and A. nobrei Pollonera, 1889 (Table 1). Because Quinteiro et al. (2005) used other ITS1 primers, we had to trim this dataset to a length of 378 bp. For each gene fragment, and for the concatenated dataset, we constructed Neighbour-Joining (NJ) trees (Saitou and Nei 1987) using the Kimura 2-parameter (K2P) model in MEGA v5 (Tamura et al. 2011) with complete deletion of insertions and deletions (indels). Branch support was evaluated with 1000 bootstrap replicates (Felsenstein 1985).

Arion ponsi
The four individuals of A. ponsi yielded four COI and three 16S haplotypes (Appendix, Supplementary Figures 1-2), yet two 16S haplotypes only differed by an indel of two base pairs at positions 291-292. For both genes A. molinae showed the smallest p-distance with A. ponsi (COI: mean p-distance 0.142 ± 0.014; 16S: mean p-distance 0.162 ± 0.019), but a sister species relationship with A. molinae was only well-supported by 16S. There were three ITS1 haplotypes for A. ponsi; one of these had a deletion of a poly-T stretch of six base pairs at positions 556-561; the other two differed by a deletion of a G at position 554. These three ITS1 haplotypes of A. ponsi clustered within a clade of A. subfuscus S1-5, A. lizarrustii, A. molinae, A. iratii and A. transsylvanus (Appendix, Supplementary Figure 3). The ITS1 analysis with the sequences of Quinteiro et al. (2005), placed the single remaining A. ponsi haplotype in the same clade (mean p-distance with the other taxa of this clade = 0.046 ± 0.004), but without bootstrap support (Appendix, Supplementary Figure 4). As for 16S, the concatenated tree of the three gene fragments showed a sister species relationship between A. ponsi and A. molinae (Figure 3).

Arion gilvus
The three A. gilvus specimens yielded two COI (one synonymous A-G substitution at position 366) and one 16S haplotypes. For both genes the smallest mean p-distances were observed relative to A. urbiae and A. anguloi (COI: mean p-distance = 0.145 ± 0.013; 16S: mean p-distance = 0.134 ± 0.016). The two A. gilvus ITS1 haplotypes reduced to one when considering the stretch that overlapped with the Quinteiro et al. (2005) sequences. In this stretch it differed from that of Quinteiro et al. (2005) by a deletion of a T at position 349. Separately, none of the three genes provided reliable evidence about the sister group relationships of A. gilvus (Appendix, Supplementary Figures 1-4). Yet, the concatenated tree showed a well-supported sister species relationship between A. gilvus and the A. urbiae / A. anguloi clade (mean p-distance = 0.021 ± 0.003) (Figure 3). Mean p-distances within this A. urbiae / A. anguloi clade (in which A. anguloi was paraphyletic) were P = 0.041 ± 0.006 for COI, P = 0.023 ± 0.006 for 16S, P = 0.004 ± 0.002 for ITS1 and P = 0.020 ± 0.003 for the concatenated dataset.

Discussion
The NJ-tree of the concatenated dataset confirms the major outcomes of previous phylogenetic studies, viz. 1) a strong support for the monophyly of the subgenus Carinarion (Geenen et al. 2006), 2) a clade of Arion s.s. and non-Portuguese A. lusitanicus (Quinteiro et al. 2005), 3) A. wiktori clustering with Mesarion species, in particular with A. paularensis (Quinteiro et al. 2005) instead of with Kobeltia species (Castillejo 1998), and 4) the strong differentiation within A. subfuscus s.s. that consists of, at least, five phylogenetic species (Pinceel et al. 2005a). It therefore seems that the analysis of COI, 16S and ITS1 DNA sequences yields relevant taxonomic information with respect to the characterisation of arionid species that have been described under the morphospecies concept.
Because Arion gilvus and Arion ponsi were originally described on morphological grounds they are to be interpreted as morphospecies. This phenetic morphological distinction, however, correlates well with a phenetic separation based on mtDNA distances. Indeed, the overall mean p-distance among the Mesarion MOTUs (excluding A. ponsi and A. gilvus) dealt with in this study is 16.7% for COI and 13% for 16S. As such, the mean p-distances between A. ponsi and A. molinae (COI: 14.2%; 16S: 16.2%) or between A. gilvus and A. urbiae (COI: 14.5%; 16S: 13.4%) are perfectly comparable with the mean p-distances among the other MOTUs and morphospecies in Mesarion. Hence, with respect to mtDNA divergence, both A. ponsi and A. gilvus, behave as other Mesarion species or putative species-level MOTUs.
Obviously, the strong COI differentiation among Mesarion taxa, and of A. ponsi and A. gilvus in particular, suggests that DNA barcoding may be a suitable identification tool for these animals. Yet, this may be a too simplistic conclusion, since stylommatophorans may show extremely high intraspecific mtDNA divergences of sometimes up to 27% (K2P-distances, but note the uncorrected p-distances are almost similar) (Thomaz et al. 1996, Chiba 1999). In addition, Davison et al. (2009) showed that in the Stylommatophora the mean interspecific K2P-distances (± 3%) can be substantially lower than the mean intraspecific K2P-distances (± 12%). Under these conditions, it becomes very difficult to define generally applicable thresholds that distinguish between intra-and interspecific sequence divergences. Such thresholds are normally associated with DNA barcoding gaps (Hebert et al. 2003), but Davison et al. (2009) were unable to detect DNA barcoding gaps in the taxa they studied. Nevertheless, Davison et al. (2009) suggested a pragmatic 4% threshold to separate intra-and interspecific values, but at the same time they also concluded that DNA barcoding in itself is insufficient to identify and/or detect stylommatophoran species. Unfortunately, our sample sizes were too small to explore eventual DNA barcoding gaps in Mesarion.
Because DNA barcoding on its own may be unreliable for identifying and detecting species-level taxa in stylommatophorans, it its necessary to backup this sort of data with, amongst others, phylogenetic analyses. As such, our phylogenetic trees of the DNA sequence data show that the morphospecies A. ponsi and A. gilvus, also represent phylo-genetic species, since both form well-supported clades that are "significantly" associated with well-defined, but morphologically different sister species. For A. ponsi, the sister species appears to be A. molinae, the distribution range of which is located in NE continental Spain (Castillejo 1997), i.e. north of, and facing, the Balearic Islands. Conversely, the sister taxon of A. gilvus is the "tandem" of A. urbiae and A. anguloi, two species that have been synonymized by Backeljau et al. (1994) and that jointly should be referred to as A. urbiae. Our DNA sequence data on COI, 16S and ITS1 (e.g. Figure 3), as well as those on ND1 and ITS1 of Quinteiro et al. (2005) are in line with this. As such, the distribution range of A. urbiae is situated northwest of, and probably adjacent to, that of A. gilvus. Thus, for both the species pairs A. ponsi / A. molinae and A. gilvus / A. urbiae, the distribution ranges appear at least consistent with the suggested sister group relationships.
In conclusion, the present work shows that A. ponsi and A. gilvus clearly differ from A. subfuscus or any other currently recognized arionid species. As such, former records of A. subfuscus from Menorca (e.g. Gasull and van Regteren Altena 1970, Mateo 1993, Beckmann 2007 almost certainly refer to A. ponsi. Similarly, probably all reports of A. subfuscus in the regions of Valencia and Albacete involve A. gilvus (e.g. Borredà 1994, Borredà andCollado 1996). Finally, Borredà (1994) wondered about the eventual relationship between A. subfuscus from Menorca and A. gilvus. The current data confirm unambiguously that these are two different species, with the former being A. ponsi. Yet, the overall phylogenetic relationships within Mesarion and many other A. subfuscuslike taxa remain to be resolved. In this context, one of the main questions is whether Mesarion in its present use is a monophyletic taxon. At the same time one may wonder about the relationships with the subgenus Arion s.s., with which Mesarion seems to form a well-supported clade (Figure 3).