Phylogenetic Relationships among Populations of the Vineyard Snail Cernuella virgata ( Da Costa , 1778 )

Cernuella virgata (Da Costa, 1778) (Mollusca: Hygromiidae), commonly known as the “vineyard snail,” is endemic species inMediterranean andWestern Europe including the British Isles, but in the Eastern USA and Australia it represents an introduced invasive species. The present work examines the genetic variability and phylogenetic relationships among the four populations of this land snail sampled along the east Adriatic region of Croatia usingmitochondrialmarkers (partial 16S rDNA andCOI gene) in addition to traditional methods of shell’s shape analysis. All the three molecular-phylogenetic approaches (median joining haplotype network analysis and Bayesian analysis, as well as maximum likelihood analysis) revealed two-three major subnetworks for both 16S rDNA and COI, with a clear distinction between south Adriatic haplotypes (Pisak) and north Adriatic haplotypes (Krk and Cres). The population from Karlobag was comprised of both north and south haplotypes, thus representing a putative contact zone between these two groups. The morphometric analysis showed that individuals from Cres island population were statistically significantly wider and higher than individuals from Pisak population. Analysis of the SW/SH ratio and the relationship between shell width and shell height showed no differences in shell growth between the two examined populations, indicating equal shell growth and shape, which gives the possibility that differences in size of individuals between those two populations could be influenced by biotic (physiological) or abiotic (environmental) factors.This study represents the first analysis of genetic variability and relatedness among native populations of C. virgata.


Introduction
Cernuella virgata (Da Costa, 1778) (Mollusca: Hygromiidae) is a terrestrial mollusks gastropod, commonly known as the "vineyard snail." Terrestrial gastropods represent a highly diverse group of mollusks.Solem [1] estimated about 24.000 species, of which about 20.500 are pulmonate stylommatophorans, and about 3.650 belong to other groups of terrestrial gastropods.Recent studies report high diversity of this group of animals and emphasize the need for more phylogenetic studies.The molecular techniques have proven useful in cases where morphological features failed to provide unequivocal phylogenetic signal by providing an independent set of characters [2].Among many different molecular markers, the cytochrome oxidase subunit I (COI gene) has been recently recommended as particularly useful for DNA barcoding of animals [3].
C. virgata is frequent in coastal regions.This species has the potential to be widely distributed because of the ability to endure long periods of warmth, dryness, fasting, and light [4].The shells of C. virgata ranging from 6 to 19 mm in height and from 8 to 25 mm in width [5] can be streaked with a variable number of pale to darker brown markings and can be extremely variable in colour pattern [6,7].
C. virgata is known to be native and endemic species to Mediterranean and Western Europe including the British Isles, but the occurrence of this species was reported in the Eastern USA [8] and Australia [6,[9][10][11][12], where it represents an introduced invasive species.In the USA, C. virgata is considered of priority quarantine importance [8], while in Australia the species has the status of a serious pest species in agriculture [10].The snails are principally regarded as pests because they foul grain harvests when aestivating on the heads and stalks of crops [12].According to the available   literature, there are no data about the genetic variability of the C. virgata populations.In this study, we performed the first combined analysis of molecular and morphological data for four C. virgata populations sampled along Croatian coast and islands-part of the native area of distribution of this snail in the Mediterranean region.

Material and Methods
2.1.Morphological Analysis.Adult individuals of Cernuella virgata have been sampled at four locations along the Croatian coastal part and islands (Table 1, Figure 1) during 2012.The two largest populations (Cres and Pisak) were subjected to the morphological analysis.Eighteen specimens of different sizes were collected from Cres island, and fourteen specimens were collected from Pisak.Two shell morphological characteristics including shell width (SW) and shell height (SH) were measured for each individual to the nearest 0.1 mm using Vernier callipers.Shell width (SW) was measured as the maximum width axis (A-B), and shell height was measured along an axis passing through the apex to the bottom of the shell (C-D) (Figure 2).Relative shell heights (SW/SH ratio) were calculated as indices of "globularity" and used to compare shell growth of both populations.A -test for independent samples was used to determine differences in two shell axes and SW/SH ratio between Pisak and Cres populations.Prior to this analysis, data were tested for homogeneity of variance using Levenès test.A linear regression was used to estimate the relationship between two shell axes. values <0.05 were considered statistically significant.
2.2.DNA Isolation, PCR Amplification, and Sequencing of DNA.Two mitochondrial genes (16S rRNA and COI) were PCR amplified with the primers described [2].The specific PCR products were purified using "QIAEX II Gel Extraction Kit 150" (Quiagen, Hilden, Germany).Samples were subjected to cycle sequencing using the ABI PRISM 3100-Avant Genetic Analyzer.

Sequence Analysis.
DNA sequences were assembled and prealigned using BioEdit ver.7.0.5.3 [13].DNA sequences were aligned in ClustalW [14] and implemented in MEGA5 [15], and the alignment was refined manually.The sequences were deposited in GenBank (Table 2).In order to avoid multiple submissions of identical sequences, we sent only one sequence of each haplotype.The number of haplotypes and the haplotype and nucleotide diversity within populations (Table 3) were determined using DnaSPv5 [16].We included into our analysis two C. virgata sequences of 16S rDNA and  a single COI sequence that were available in the GenBank [2,17], as well as several sequences that we used as outgroups [2,17,18].Phylogenetic analyses were carried out using maximum likelihood (ML) and Bayesian-based inference (BI) methods.The most suitable model of nucleotide evolution was determined by the Akaike Information Criterion, AIC, [19] as implemented in JModelTest [20].Bayesian analysis was performed with MrBayes 3.1.[21] with 4 chains of 1,000,000 generations, trees were sampled every 100 generations and burnin value set to 25% of the sampled trees.Maximum likelihood analysis using starting trees obtained by neighbor joining and TBR branch swapping with model parameters was performed using PAUP * 4.0b10 [22].The number of bootstrap replicates was set to 1000.Phylogenetic trees were displayed in FigTree v1.3.1.
Alternative procedures to phylogenetic tree reconstruction were also conducted to depict the relationships among the C. virgata haplotypes.We used the network-based approach which seems better suited for representing relationships among closely related sequences [23].The median joining (MJ) algorithm implemented in Network v4.6.1.0.software [24] was used with for constructing networks (weight = 10 and  = 10).Median networks that contained all possible equally short trees were simplified by running the maximum parsimony (MP) calculation option to eliminate superfluous nodes and links.
In order to infer the demographic history of populations, that is, to detect past population expansion, we used Fu's   statistics [25], which are based on the distribution of haplotypes.  statistics were estimated for overall Croatian population and for each of the four populations for the two genes analyzed using DnaSPv5 [16], and significances were assessed using 1000 coalescent simulated resampling. ) ranged from 6.6 mm to 14.0 mm, with a mean (±SD) of 11.3 ± 1.9 mm, whereas Pisak population ( = 14) ranged from 6.2 mm to 8.5 mm, with a mean (±SD) of 7.8 ± 0.7 mm.A -test for independent samples showed that the differences in both shell axes between Cres and Pisak populations were significant ( > 0.05).

Results and Discussion
The estimated SW/SH ratio of analyzed shells from Cres island ranged between 0.71 and 0.93 with a mean of 0.78±0.05,whereas SW/SH ratio of Pisak population ranged between 0.69 and 0.82 with a mean of 0.76 ± 0.04.There were, thus, no statistically significant differences in the estimated SW/SH ratio of both populations ( = 1.15;  = 0.259), indicating the same shell growth pattern for Cres and Pisak populations; that is, as the shell width increased the shell height increased equally in both populations (Figure 4).According to this, a high positive correlation was found between shell width and shell height of Cernuella virgata from Cres ( = 0.9337) and Pisak ( = 0.8530).These results indicate the possibility that differences in size of individuals from those two populations could be influenced by biotic (physiological) or abiotic (environmental) factors.By comparing the size of shells between  the Australian C. virgata populations, Swofford [22] concluded that it was inversely related to population density.

Genetic Variability of C. virgata Populations along the East
Adriatic Coast.The genetic variability of C. virgata based on four populations and 39 specimens collected along the east Adriatic coast and islands of Croatia was analyzed using two mitochondrial molecular markers: 16S rDNA and cytochrome oxidase I (COI).The genetic diversity information in the mtDNA is given in Table 3.The 16S rRNA gene sequences were 407 bp long, and the COI gene sequences were 641 bp long.Comparison of sequence alignments showed COI as more variable than 16S rRNA.The nucleotide diversity () was 2.9 times higher for COI than for 16S rRNA genes, and we identified nine COI haplotypes as opposed to six 16S rDNA haplotypes.Among the four populations, the highest nucleotide diversity () and haplotype diversity (  ) for both genes had population from Karlobag.Despite numerous base substitutions, there was no change in the amino acid composition of the COI protein in all the 39 specimens.
Demographic parameters used to test for demographic events were estimated for each gene, four populations, and overall Croatian haplotypes.Fu's test was significant only for the population from Cres for both genes:   = 5.809,  < 0.05 for 16S and   = 9.778,  < 0.01 for COI, Table 3. Fu's tests were performed at a higher geographic scale; that is, all four C. virgata populations also rejected the null hypothesis of constant size (  = 5.964,  < 0.01 for 16S;   = 18.723,  < 0.001 for COI, Table 2).Since the result which indicates an expansion of the overall Croatian C. virgata populations is based on results for Cres population, more populations should be analyzed before the final conclusion.One should bear in mind that transportation by human activities could also influence the stability of the C. virgata populations.

Phylogenetic and Network Relationships among Croatian C. virgata Haplotypes.
To describe the fine-scale spatial genetic structure of Croatian populations of C. virgata, we performed phylogenetic analysis of the newly obtained Croatian 16S rRNA and COI haplotypes using three different approaches: median joining network analysis (MJ), Bayesian inference (BI), and maximum likelihood (ML) analyses.
MJ network revealed three 16S rDNA and two COI major subnetworks (Figure 5), with a clear distinction between south Adriatic haplotypes (Pisak) and north Adriatic haplotypes (Krk and Cres).The topology of the 16S rDNA tree (Figure 6) and COI tree (Figure 7) confirmed the haplotype  2.
grouping obtained by MJ.South Adriatic clade comprised of individuals from Pisak and two individuals from Karlobag was better supported than the north Adriatic clade containing sequences from islands of Krk, Cres, and Karlobag.Both the northernmost population (Krk) and the southernmost population (Pisak) had low haplotype diversity and were quite uniform.On the contrary, the populations of Cres and Karlobag were comprised of two and three distinct haplotypes, respectively, and possessed higher haplotype and nucleotide diversity.Based on these results, the population from Karlobag could represent a putative contact zone between these two groups.The lack of star-like structure indicated that any Croatian haplotype is presumed to be ancestral.
Our COI tree shows that the two snail species, Suboestophora hispanica (Gude, 1910) and Oestophora tarnieri (Morelet, 1854), which we used as outgroups, were more closely related to C. virgata than two species from the same genus: Cernuella cisalpina (Rossmässler, 1837) and Cernuella neglecta (Draparnaud, 1805).These results are congruent with the results of Steinke et al. [2], who made the first comprehensive molecular analysis of the Helicidae s.l.phylogeny based on two mitochondrial (COI, 16S rDNA) and two nuclear (ITS-1, 18S rDNA) fragments.Their results showed that the members of the genus Cernuella appeared in several clades within the Helicellinae, indicating that this genus might be polyphyletic.
Based on the morphological and molecular-phylogenetic analyses of the four populations of C. virgata, sampled along Croatian coast and islands, we found differences between the north Adriatic and the south Adriatic populations.Such a result is concordant with the results of phylogeographic analyses of two other land snails distributed along Croatian coast and island: Eobania vermiculata (O.F. Müller, 1774) [26,27]   and Cornu aspersum (Müller, 1774) (Puizina et al., unpublished).We assume that such a spatial distribution of the snail haplotypes could be the consequences of climate changes and glaciations in the Quaternary [28][29][30], but additional analyses of more populations of all the three snail species are necessary to confirm this hypothesis.

Figure 1 :
Figure 1: A map of Croatia showing the locations of the sampling sites: Krk island, Cres island, Karlobag, and Pisak.

Figure 2 :
Figure 2: The shells of Cernuella virgata from (a) and (b): Cres island and (c) and (d): Pisak; measurements made on shells: (A-B) the maximum shell width axis and (C-D) the maximum shell height axis.Scale bar = 10 mm.

3. 1 .
Morphological Analyses.The morphometric data obtained for individuals from both populations of C. virgata are summarised in Table 1 and illustrated in Figure 3. Shell width of individuals from Cres population ( = 18) ranged from 8.5 mm to 18.2 mm, with a mean (±SD) of 14.6 ± 2.5 mm, whereas individuals from Pisak population ( = 14) ranged from 9.0 mm to 11.2 mm, with a mean (±SD) of 10.2±6.6 mm.Shell height of Cres population ( = 18

Figure 5 :
Figure 5: Median joining networks for the 16S rDNA haplotypes and COI haplotypes.Each circle represents a haplotype, and circle size is shown proportional to haplotype frequency.Colors indicate the geographic origin of haplotypes.Median vectors (small red dots) are produced by the network program, representing missing or not sampled haplotypes.Branch lengths are approximately equal to inferred mutational steps.Haplotype codes according to those are represented in Table2.

Figure 6 :
Figure 6: Phylogenetic trees resulting from a Bayesian analysis of 16S rDNA sequences.The numbers above the branches depict Bayesian posterior probabilities, and the numbers below the branch indicate bootstrap support values from maximum likelihood analysis (in the case of nodes not supported by all methods, the respective missing support values are indicated by "n.a.").Bar indicates substitutions/site.

Figure 7 :
Figure 7: Phylogenetic trees resulting from a Bayesian analysis of COI sequences.The numbers above the branches depict Bayesian posterior probabilities, and the numbers below the branch indicate bootstrap support values from maximum likelihood analysis (in the case of nodes not supported by all methods, the respective missing support values are indicated by "n.a.").Bar indicates substitutions/site.

Table
: Morphometric data and SW/SH ratio obtained for Cernuella virgata from Cres island and Pisak.Values are presented as means ± standard deviations (minimum-maximum).

Table 2 :
16S rRNA and COI haplotypes found in Cernuella virgata populations analyzed with their corresponding Genbank accession numbers.N: number of individuals.

Table 3 :
Genetic diversity estimates and results of demographic tests for 16S rDNA and COI of specimens of Cernuella virgata sampled in east Adriatic coast and islands.