New records of the endemic Sicilian land snail species Marmorana (Murella) muralis (O. F. Müller, 1774) from the north of Tunisia (Pulmonata, Gastropoda)

Abstract Marmorana (Murella) muralis is known as an endemic species of Sicily Island, which is introduced in many European countries. Here, M. (M.) muralis is recorded from the north of Tunisia. In order to confirm the identification of samples collected from several localities, shell morphology, details of genital organs and two mitochondrial markers (COI and 16S) were investigated. The results of the molecular study, as well as the morphological and anatomical studies confirm the identification of all Tunisian samples as M. (M.) muralis. The analysis of mitochondrial markers shows a low divergence between Sicilian and Tunisian samples suggesting a recent introduction of M. (M.) muralis to the North of Tunisia. The comparison of morphological characters of M. (M.) muralis with shell characters of Murella nicollei described by Pallary (1926) confirms that the latter should be considered as synonym of M. (M.) muralis.


Introduction
Land snails compose a group of invertebrates which are characterized by low mobility and dispersal capacity. The evolution of morphological characters within land snail species is widely influenced by the environmental and ecological conditions. Marmorana (Murella) muralis is an endemic helicid species from Sicily Island, which is characterized by an extremely high variability of shell morphology as well as molecular characters (Fiorentino et al. 2013). It was demonstrated that paleogeographical factors and environmental changes affected the shell morphology of M. (M.) muralis in Sicily (Fiorentino et al. 2013). This species was introduced by humans to many other European areas such as Tuscany in Italy, Sardinia, the Baleares, Portugal and Bouchesdu-Rhône in France. Tunisia is a quite well sampled area as evidenced by Letourneux and Bourguignat (1887), who documented land snails from a plethora of localities. Interestingly, they never recorded the presence of Marmorana (Murella) Pfeiffer, 1877 in the area. It was Pallary in 1926, who was the first to describe a Murella, Murella nicollei, from Tabarka in northwest Tunisia.
Recent sampling efforts by the senior author revealed the presence of a Marmorana (Murella) taxon at several localities in the north of Tunisia. The present study aims to 1) identify the samples collected from Tunisia based on morphological and molecular characters, 2) determine the possible origin of each Tunisian population known and 3) clarify the status of Murella nicollei Pallary, 1926.

Materials and methods
Living specimens were collected by hand at several localities in Tunisia during two periods: spring 2014, and winter 2015/2016. Geographic coordinates were recorded using a GPS device. For subsequent molecular analyses, specimens were preserved and stored in 80% ethanol until dissection and DNA extraction.

Morphological and anatomical studies
First assessments of the shell morphological characters were done by using simple magnifying glasses. Preserved animals were dissected under a LEICA M212 stereo microscope using thin tweezers. The genital organs of the specimens were removed from the body, and the outer morphology of the complete hermaphroditic genital organ (situs) and further morphological details were investigated. After that, shells, genital situs, and details of the genital organs were photographed with a LEICA DFC 425 camera combined with a LEICA M205 C stereo microscope. The multifocal images were processed by using Imagic IMS software (Imagic, Switzerland).

DNA extraction, PCR amplification and sequencing
Total genomic DNA was extracted from foot muscle tissue of each specimen using a standard phenol chloroform method (Estoup et al. 1996). Two mitochondrial gene fragments were chosen for analyses in the present study: cytochrome c oxidase subunit I (COI) of 711 base pairs (bp) length and the gene of the 16S ribosomal RNA subunit (16S rRNA) for an approximately 470-480 bp fragment. Polymerase chain reactions (PCR) were performed in a reaction mixture containing 15 ng of DNA template, 1X reaction buffer (1.5 mM), 0.1 mM of each primer pair, 0.2 mM dNTPs, Taq polymerase (1.25 U) and adjusted till a total volume of 25 μl with DNAase free water/sterilized water (UNIMED) (H 2 O). PCR reactions were run under the following conditions: 3 min at 95 °C, followed by 35 cycles of 1 min at 95 °C, 1min at 40 °C and 1 min at 72 °C and finally, 5 min at 72 °C for COI. For 16S the amplification conditions were: 3 min at 95 °C, followed by 35 cycles of 1 min at 95 °C, 1 min at 50 °C and 1 min at 72 °C, and finally, 5 min at 72 °C. PCR products were sequenced using automated and standardised ABI 3730 XL sequencing run with a read length up to 1.100 bp (PHRED20 quality) and using the same primers as for the PCR (Table 2). 5'-CGCCTGTTTATCAAAAACAT-3' 5'-CCGGTCTGAACTCTGATCAT-3' Simon et al. 1994

Sequence analyses
Forward and reverse sequences were assembled, checked for ambiguities and aligned using the default settings of the ClustalW multiple alignment algorithm as implemented in Bioedit V 7.2.5 (Hall 1999) and trimmed for 655 bp and 414 bp respectively for COI and 16S. Obtained sequences were deposited in GenBank under the accession numbers MG780362-MG780370 and MG774439-MG774448 (Table 1). Aligned Tunisian and Sicilian sequences were analysed using DnaSP v5.10.01 (Librado and Rozas 2009) to estimate the number of informative sites and nucleotide diversity for each marker. The K2P values were estimated using Mega v.6 (Tamura et al. 2013). The relationships of inferred haplotypes of Tunisian and Italian M. (M.) muralis were estimated using the Minimum Spanning Network (MSN) method (Bandelt et al. 1999) implemented PopART v1.7 (Leigh et al. 2015). Because of lack of sequences available on GenBank, we produced the haplotype network separately for COI and 16S markers.

Phylogenetic analysis
Concatenated sequences of the two mitochondrial markers were analysed by Bayesian inference of phylogeny. The sequence data was initially partitioned into four partitions: three partitions corresponding to the codon positions of COI and one partition for16S. Based on the Akaike Information Criterion (AIC), the substitution models F81, K81uf+G, TrN+I and HKY+G were chosen as best models, respectively, for the first, second and third codon positions of COI and for 16S by PartitionFinder v 1.1.1 (Lanfear et al. 2012). For the Bayesian Inference, we used Mr Bayes v3.2.2 (Ronquist and Huelsenbeck 2003) using the partition scheme and substitutions models suggested by PartitionFinder. Four independent runs were conducted for 10 6 generations, sampling every 1000. The first 25% trees were discarded as default burn-in and a majority-rule consensus tree was calculated from the remaining trees. Convergence between runs was assessed by comparing the traces using Tracer v1.5 (Rambaut and Drummond 2007). The topology obtained, and the posterior probabilities for each node were displayed with Figtree V1.4.0 (Rambaut 2012).
Male genital anatomy. Penis club-shaped, thick; epiphallus as long as penis; retractor muscle inserting into the distal part of the epiphallus; flagellum twice the length of epiphallus; penial papilla elongated, with a slit-like pore on one side.
Female genital anatomy. Dart sac simple, well developed, two glandulae mucosae, non-ramified, inserting into the middle part of the vagina near the base of the dart sac; bursa copulatrix and diverticulum inserting into the proximal part of the vagina (Fig. 2).

Haplotype network and genetic diversity
Among nine Tunisian and 58 Italian partial COI sequences of M. (M.) muralis (Fig.  3), 46 distinct haplotypes were found, suggesting an extremely high haplotype diversity (Hd = 0.9815) (Fig. 4A). With 45 haplotypes detected, Italian sequences are   highly diverse (Hd = 0.9903), while only 3 haplotypes were found in Tunisia (Hd = 0.5596). The haplotype network therefore suggests a relatively low genetic variability of COI sequences from Tunisian specimens compared to sequences from Italian specimens. Tunisian and Italian specimens share two haplotypes: the first is represented by the COI sequences of the samples collected from Manzel Jemil, Manzel Abderrahmen and the sequence from Selinunte. The second haplotype is represented by the sequence of the sample collected at Haouaria and the sequences H1, H4, and H50 from Erice (Fiorentino et al. 2013). The sequences of the samples collected at Kelibiya represent a unique haplotype that is neither shared with the other specimens from Tunisia nor with any of the specimens from Italy. Within the Tunisian sequences, the highest K2P value (0.078) was recorded between the haplotype of the sequence from Haouaria and the sequences from Kelibiya however; the lowest value (0.01) was reported between the sequences from Kelibiya and those from Manzel Jemil and Manzel Abderrahmen. Between Tunisian and Italian populations, the highest K2P value (0.081) was registered between the sequences from Kelibiya and the sequences from Erice (H5, H10) and Monte Monaco (H22, H28, H30, H38, H52). The lowest was recorded between the sequence from Haouaria and the sequences H1, H4 and H50 (Fiorentino et al. 2013) on the one hand and the sequences from Manzel Jemil-Manzel Abderrahmen and the sequences from Selinunte on the other hand (Fiorentino et al. 2010). The nucleotide divergence within Tunisian population reached a value of 0.0178 but the divergence was slightly higher (0.0316) within Italian population. Moreover, the divergence between Tunisian and Italian populations reached a value of 0.0353.
The analysis of ten Tunisian and six Italian 16S partial fragments shows five haplotypes suggesting a low haplotype diversity (0.450) (Fig. 4B). Sequences of Italian specimens represent four haplotypes (Hd = 0.80), while sequences from Tunisian specimens represent only two haplotypes (Hd = 0.20). Tunisian and Italian samples share one haplotype, which was represented by sequences of Tunisian specimens from Kelibiya, Manzel Jemil, Manzel Abderrahmen and sequences of Italian specimens from Lazio, Selinunte and Marsala. The sequences of specimens from Haouaria, Fiumedi- nisi, Caltabellotta and Joppolo represented four different haplotypes. The maximum value of K2P distance (0.02), within Tunisian sequences, is recorded between the sequence from Haouaria and the rest. While the maximum value recorded between Tunisian and Italian populations is 0.044 between the sequence from Haouaria and the sequences from Fiumedinisi and Caltabellotta. The nucleotide divergence of the 16S partial fragment is remarkably low within Tunisian population (0.00409), as well as, between Tunisian and Italian populations (0.00828).
The analysis of nine Tunisian and six Italian concatenated sequences (COI, 16S) recovered seven different haplotypes among them (Fig. 4. C). One haplotype is shared by Tunisian and Italian populations. The rest is divided into two Tunisian and four Italian haplotypes.

Phylogeny
The topology, obtained by Bayesian inference based on the concatenated COI and 16S data set was rooted with Macularia sylvatica and Macularia niciensis as outgroups (Fig. 5)

Morphology and anatomy
Marmorana (Murella) muralis is known as an endemic species of Sicily but it was introduced to several localities in southern Europe (Fiorentino et al. 2008). The morphological and anatomical characters of the Tunisian samples show the same morphological and genital anatomical traits presented by Fiorentino et al. (2010). Thus, these specimens are here considered to represent M. (M.) muralis. In Tunisia, this taxon was first recorded by Pallary (1926: 49, pl. VIII, fig. 9) under the name Murella nicollei from Tabarka. The photo of Murella nicollei (Fig. 1C) Tabarka, or 2) the fragmentation and urbanisation of its habitat by human activities, which easily could reduce the population. Being an alien species to Tabarka it is quite possible that it could not well disperse in the area. As a result, the population is gone extinct. However, despite its extinction in Tabarka, it does well in the other Tunisian localities recorded here. Fiorentino et al. (2013) demonstrated that the shell morphology is highly affected by environmental changes in Sicily Island; the Tunisian populations seem not yet to be influenced by the new environment so far. The absence of any environmental effects on the shell supports the hypothesis that the species was quite recently introduced to the country.

Network haplotype and genetic diversity
The nucleotide divergence of the COI sequences reaches a maximum value of 0.0316 (3.16%) between Tunisian and Italian populations. This value does not exceed the threshold of intraspecific divergence of land snails (4%) as suggested by Davison (2009), and is comparable to the threshold of 3% suggested by Hebert et al. (2003) to characterize animal species in general. Furthermore, this value is smaller than the intraspecific divergence of the Tunisian Xerocrassa latastei reported by Ezzine et al. (2017). The comparison of the nucleotide divergence, the haplotype diversity, and the K2P value between Tunisian and Italian COI sequences shows a high diversity of this marker. The divergence between Tunisian and Italian populations might be the result of the isolation caused by the Mediterranean Sea, which can be considered a geographical barrier causing the restriction of passive gene flow between the two populations.
The analysis of the results obtained by the 16S sequences shows low values of nucleotide divergence, haplotype diversity, and K2P distance between Italian and Tunisian sequences, suggesting a weak diversity of this marker. The comparison of the parameters of the COI and 16S and the haplotype network show that COI is more polymorphic than 16S. COI seems to be suitable to estimate the divergence not only on species but also on population level. The haplotype network of the concatenated data confirms the results obtained by COI and 16S separately and shows that Italian populations are more diversified than the Tunisian ones. This supports the hypothesis of a recent introduction to Tunisia.
The haplotype network of COI sequences shows that the haplotype from Manzel Jemil and Manzel Abderrahmen is similar to the haplotypes from Selinunte and Joppolo, the haplotype from Haouaria is similar to the sequences from Erice, which can be interpreted as a hint to the origin of these particular populations. Interestingly, the haplotype from Kelibiya is unique and not shared with Italian populations. The divergence of the haplotype of Kelibiya may have two reasons: 1) these snails have been introduced from a genetically unknown population on Sicily, or 2) or it could be the result of the isolation of the population inside the castle. In fact, we visited Kelibiya several times, the population seems to be isolated but well adapted to the environment within the castle. The species does not live outside the castle. Geographical isolation is widely accepted to represent the main cause of genetic divergence within a species (Graybeal 1995;Baum and Shaw 1995;Olmstead 1995). However, this is a process that requires many generations and might lead to changes in shell morphology as seen in Sicily. This is not the case here, so we assume that the first hypothesis has a higher probability.

Phylogeny
The analysis of the topology obtained by the Bayesian Inference method shows that Tunisian specimens form one well supported clade (PP = 1) together with the Italian samples of M. (M.) muralis (Fig. 5) and thus proves that the Tunisian samples are conspecific with this species. The topology obtained did not divide the samples used into separate Sicilian and Tunisian clades, and the presence of a Tunisian or Sicilian ancestral clade could not be shown. Additionally, the shell morphology seems not to be affected by the environmental difference between Sicily and Tunisia, as might have been expected in case of a longer presence of the species in Tunisia.
To better understand the population dynamics of this species, more studies including more samples from Tunisia and from Italy will be needed.

Conclusions
Based on morphological, anatomical, and mitochondrial markers, the present study confirms that the recently collected Tunisian samples of a Marmorana species belong to M. (M.) muralis. The absence of this species in the collection of Letourneux and Bourguignat (1887) leads to the hypothesis that the species may have recently been introduced to Tunisia, i.e. earliest after 1887. The first record for the species comes from Tabarka (Pallary 1926), but the species has gone extinct there. The recent populations from Tunisia share some Sicilian haplotypes indicating an origin from Selinunte, Erice and other well-known populations on Sicily; the population from Kelibiya is more isolated and does not relate to any genetically known population on Sicily. The haplotype networks of the COI, 16S and concatenated fragments prove that Italian populations are more diversified than the Tunisian. The shell morphology of the Tunisian populations is rather homogenous. We therefore conclude that the present distribution pattern is result of a recent anthropogenic introduction of the species in the north of Tunisia, which occurred sometime in the last 90 years. The species has to be considered a neozoon for the Tunisian malacofauna. It has to be emphasized that the development of the hitherto known four populations and their future dispersal in the country need to be observed. The impact of this alien species on the endemic land snail fauna of Tunisia needs serious future monitoring.