Stand out from the Crowd: Small-Scale Genetic Structuring in the Endemic Sicilian Pond Turtle

The geographical pattern of genetic diversity was investigated in the endemic Sicilian pond turtle Emys trinacris across its entire distribution range, using 16 microsatellite loci. Overall, 245 specimens of E. trinacris were studied, showing high polymorphic microsatellite loci, with allele numbers ranging from 7 to 30. STRUCTURE and GENELAND analyses showed a noteworthy, geographically based structuring of the studied populations in five well-characterized clusters, supported by a moderate degree of genetic diversity (FST values between 0.075 and 0.160). Possible explanations for the genetic fragmentation observed are provided, where both natural and human-mediated habitat fragmentation of the Sicilian wetlands played a major role in this process. Finally, some conservation and management suggestions aimed at preventing the loss of genetic variability of the species are briefly reported, stressing the importance of considering the five detected clusters as independent Management Units.

The taxonomic status of the Sicilian pond turtle is currently debated. When it was described, species status was inferred by nuclear-genomic fingerprinting [2] and later confirmed using nuclear microsatellite markers [5,7]. However, Speybroeck et al. [11] proposed treating the taxon as a subspecies of E. orbicularis, reflecting the weak genetic divergence of E. trinacris and E. orbicularis compared to other European reptiles. On the other hand, Speybroeck et al. [11] conceded that further distinct species could be currently subsumed under E. orbicularis. Among turtles, several taxa with less genetic divergence are currently recognized as full species (e.g., Graptemys [12][13][14]). To avoid unnecessary nomenclatural changes, we here continue to treat E. trinacris as a distinct species until new genetic evidence becomes available.
Regardless of the taxonomic status, the Sicilian endemic pond turtle is listed, as well as Emys orbicularis, as "EN" (i.e., "Endangered") in the Italian IUCN Red List of Threatened Species [15].
Moreover, it is included in both Appendix II of the EU Council Directive 92/43/EEC ("Habitats Directive") and in Appendix II of the "Bern Convention on the Conservation of European Wildlife and Natural Habitats".
Based on a mitochondrial marker (cytochrome b gene) and 15 microsatellite loci, a strong genetic structuring of the Sicilian pond turtle was found throughout the island by Vamberger et al. [7]. In agreement with these results, attention should be paid when individuals are translocated. Such an approach is indicated as pivotal to avoid the loss of genetic variability through a homogenization process. Even though the sampling was adequate in the study by Vamberger et al. [7], some locations were not sampled or represented only by few individuals. Based on 16 nuclear microsatellite loci and a significantly increased sampling, we here reassess the genetic structuring of E. trinacris populations in Sicily. We provide a sound population genetic database covering the whole distribution of E. trinacris in Sicily for management and conservation actions for this threatened pond turtle.

Sampling, DNA Extraction and Selected Loci Amplification
Sicilian pond turtles were collected in permanent water bodies throughout Sicily ( Figure S1; Table S1; no precise geographical coordinates are provided since the species is subject to poaching). Turtles were captured using baited hoop traps, immersed in water and equipped with a flotation system, so that turtles and other air-breathing animals that were eventually caught could breathe at the surface [28]. The traps were left fishing for about 48 h per session, and then checked for the possible presence of captured turtles.
Tissue, blood or saliva samples were collected from each turtle. DNA extraction, PCRs and genotyping followed the protocols described in Vamberger et al. [7].
The fragment lengths were determined on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster, CA, USA) using the GeneScan-600 LIZ Size Standard (Applied Biosystems, Foster, CA, USA) and the software PEAK SCANNER 1.0 (Life Technologies, Carlsbad, CA, USA).

Genetic Cluster Analysis and Diversity Indices
Microsatellite genotypes and sample spatial location data were analysed for all loci and samples using two Bayesian clustering software STRUCTURE 2.3.3 [29,30] and GENELAND 4.0.9 [31]. In STRUCTURE, the admixture model and correlated allele frequencies were used, setting an arbitrary upper bond value of K = 10; the most likely number of clusters (K) was determined using the ∆K method [32] with the software STRUCTURE HARVESTER [33] through its online interface (http://taylor0.biology.ucla.edu/structureHarvester/). Calculations for each K were repeated 10 times using a MCMC chain of 750,000 generations for each run, including a burn-in of 250,000 generations. Population structuring and individual admixture were visualized using the software DISTRUCT 1.1 [34]. Following Randi [35], individuals with proportions of cluster membership below 80% were treated as having mixed ancestries. Afterwards, data of admixed individuals were removed to avoid introducing noise from alleles of other clusters.
In GENELAND, the geographic information was used to detect spatial delineation of genetic discontinuities. The number of clusters (K) was inferred using 500,000 MCMC iterations with 100 burn-in generations.
In addition, microsatellite data of the inferred clusters were analysed by calculating population genetic diversity indices using CONVERT 1.31 [36] and GENETIX 4.05.2 [37]. Hardy-Weinberg equilibrium (HWE) for all loci and all clusters was tested in GENEPOP ON THE WEB 4.7 [38,39] and the pairwise F ST values were computed through the software ARLEQUIN 3.5.1.2 [40]. Tests for recent bottlenecks were performed with the program BOTTLENECK 1.2.02 [41] using the strict one-step stepwise mutation model (SMM; Ohta and Kimura [42]) as suggested by Piry et al. [41] for microsatellite loci. In order to compare our results, a Discriminant Analysis of Principal Components (DAPCs) was computed using the same dataset, with the package ADEGENET [43,44] for CRAN R 3.6.3.

Results
Altogether, 245 samples of the Sicilian pond turtle from 21 different sampling sites were studied ( Figure S1; Table S1), including 62 new samples and 183 samples from Vamberger et al. [7], thus covering the whole distribution of the species in Sicily. All studied microsatellite loci were highly polymorphic, with allele numbers ranging from 7 to 30, and a total allele number of 293.
The highest support in both STRUCTURE and GENELAND analyses was found for K = 5, suggesting the split of the Sicilian pond turtles into five different clusters: the Nebrodi mountains (cluster A), the Sicilian south-western area (cluster B); south-eastern Sicily (cluster C); the Sicani mountains (province of Palermo; cluster D); north-western Sicily (cluster E) (Figure 1 and Figure S2). DACP analysis (Figure 2) supported the differentiation of the five clusters recognized by STRUCTURE and GENELAND ( Figure S2). In addition, microsatellite data of the inferred clusters were analysed by calculating population genetic diversity indices using CONVERT 1.31 [36] and GENETIX 4.05.2 [37]. Hardy-Weinberg equilibrium (HWE) for all loci and all clusters was tested in GENEPOP ON THE WEB 4.7 [38,39] and the pairwise FST values were computed through the software ARLEQUIN 3.5.1.2 [40]. Tests for recent bottlenecks were performed with the program BOTTLENECK 1.2.02 [41] using the strict one-step stepwise mutation model (SMM; Ohta and Kimura [42]) as suggested by Piry et al. [41] for microsatellite loci. In order to compare our results, a Discriminant Analysis of Principal Components (DAPCs) was computed using the same dataset, with the package ADEGENET [43,44] for CRAN R 3.6.3.

Results
Altogether, 245 samples of the Sicilian pond turtle from 21 different sampling sites were studied ( Figure S1; Table S1), including 62 new samples and 183 samples from Vamberger et al. [7], thus covering the whole distribution of the species in Sicily. All studied microsatellite loci were highly polymorphic, with allele numbers ranging from 7 to 30, and a total allele number of 293.
The highest support in both STRUCTURE and GENELAND analyses was found for K = 5, suggesting the split of the Sicilian pond turtles into five different clusters: the Nebrodi mountains (cluster A), the Sicilian south-western area (cluster B); south-eastern Sicily (cluster C); the Sicani mountains (province of Palermo; cluster D); north-western Sicily (cluster E) (Figure 1 and Figure S2). DACP analysis (Figure 2) supported the differentiation of the five clusters recognized by STRUCTURE and GENELAND ( Figure S2).  (Table  S1). Within each cluster, an individual turtle is represented by a vertical segment that reflects its ancestry. Mixed ancestries are indicated by differently coloured sectors corresponding to inferred genetic percentages of the respective cluster. Colours of pooled sampling sites in the map correspond to STRUCTURE clusters; slices represent turtles with mixed ancestries or conflicting cluster assignment (percentages).  (Table S1). Within each cluster, an individual turtle is represented by a vertical segment that reflects its ancestry. Mixed ancestries are indicated by differently coloured sectors corresponding to inferred genetic percentages of the respective cluster. Colours of pooled sampling sites in the map correspond to STRUCTURE clusters; slices represent turtles with mixed ancestries or conflicting cluster assignment (percentages). Observed heterozygosity (Ho) values (range 0.54-0.66) were lower than the expected ones (He, range 0.60-0.76) ( Table 1) In addition, all clusters showed evidence of recent bottlenecking (see Table  1). Table 1. Genetic diversity indices and BOTTLENECK analysis of STRUCTURE clusters of Emys trinacris. Colours refer to clusters (see Figure 1). The category admixed individuals comprises turtles with admixed genotypes irrespective of their collection sites. Moreover, all identified clusters of E. trinacris showed a significant deviation from HWE at all loci and within all clusters defined by STRUCTURE and GENELAND (Table 2). Pairwise FST values between the five a priori defined cluster, based on the STRUCTURE results, ranged from 0.075 to 0.160 (Table 3); all FST values were significant with a p-value of 0.05. FST values indicated moderate genetic differentiation among sampled Sicilian pond turtles, where the highest value was observed between the comparison of the clusters "D" and "C". Observed heterozygosity (Ho) values (range 0.54-0.66) were lower than the expected ones (He, range 0.60-0.76) ( Table 1) In addition, all clusters showed evidence of recent bottlenecking (see Table 1). Table 1. Genetic diversity indices and BOTTLENECK analysis of STRUCTURE clusters of Emys trinacris. Colours refer to clusters (see Figure 1). The category admixed individuals comprises turtles with admixed genotypes irrespective of their collection sites. Moreover, all identified clusters of E. trinacris showed a significant deviation from HWE at all loci and within all clusters defined by STRUCTURE and GENELAND (Table 2). Pairwise F ST values between the five a priori defined cluster, based on the STRUCTURE results, ranged from 0.075 to 0.160 (Table 3); all F ST values were significant with a p-value of 0.05. F ST values indicated moderate genetic differentiation among sampled Sicilian pond turtles, where the highest value was observed between the comparison of the clusters "D" and "C".

Discussion
The results obtained in this study suggest the presence of a clear structuring of Emys trinacris populations, with a moderate degree of genetic differentiation among the five identified clusters ( Table 3). The clusters are supported by moderate genetic diversity indices among them (see Table 1); a similar moderate genetic structuring can be also observed in E. orbicularis [45]. Surprisingly however, the genetic structuring only partially agrees with that revealed by Vamberger et al. [7]. While Vamberger et al. [7] found eight well differentiated clusters within E. trinacris, our increased sampling covering some new sampling localities concluded that only five clusters exist. Nevertheless, three of them were in full agreement with the ones from Vamberger et al. [7], i.e., clusters A, C, D. Both studies highlighted a moderate genetic structuring of E. trinacris in Sicily, with a slightly different grouping of the studied samples, possibly due to an increased sampling scheme and the implementation of a further microsatellite marker, suggesting a long-lasting persistence of turtles in the respective areas of Sicily. In contrast, the subspecies of E. orbicularis do not show such a pronounced substructuring (e.g., [7,8]). Stöck et al. [46] found low levels of genetic structuring in the studied Sicilian herpetofauna using two mitochondrial markers, the cytochrome b and 16S rDNA genes, and a nuclear marker, the tropomyosin intron. However, it is well known that the current genetic structuring of many other freshwater taxa in Sicily can be ascribed to paleogeographic and paleoclimatic events [47][48][49]. The pronounced genetic structure in the Sicilian pond turtle could thus be due to the role of Sicily as a "refuge area" during the Plio-Pleistocene climatic fluctuations [50]. The local physiography can maintain genetic diversity (see Vamberger et al. [7]), corresponding to a "refugia within refugia" pattern [51].
In addition, increased habitat fragmentation and the consequent reduction in population size and gene flow levels [52] could also contribute to the current structuring, which shows evidence of recent bottleneck experienced by all the local clusters of the species (Table 1). The paucity of water bodies and habitats suitable for the presence of the species in Sicily prevents the populations from establishing sound metapopulation dynamics, thus precluding significant inter-population gene flow. In some areas, Emys trinacris is locally abundant but these regions are separated from other inhabited areas by large gaps where the species is completely absent (e.g., in the Madonie and Peloritani mountain chains).
In recent centuries, the wetlands in Sicily have been significantly modified due to remediation measures to reclaim agricultural land and to fight against malaria ("Provveditorato alle opere pubbliche per la Sicilia" [53]; see also Naselli-Flores and Marrone [23]). During those years, wetlands across Sicily were altered completely by local authorities by re-shaping the surrounding landscape. Furthermore, according to ISPRA ("Istituto Superiore per la Protezione e la Ricerca Ambientale" [54]), in the last century Sicily has undergone a severe degradation of the wetlands, where 34.5% of all the Sicilian wetlands are in a state of complete degradation and only 5% remain untouched (cf. Table 10 in [54]).
In conclusion, to preserve the genetic variability of E. trinacris, protection actions by local management entities are required. According to our improved genetic database, which can now serve in conservation actions, we propose that the five population clusters highlighted by present study should be considered as independent Management Units for Sicily. As already proposed by Vamberger et al. [7], it is necessary to avoid translocating turtles, even when events or activities regarding the restocking or release of treated or seized individuals are carried out. It is thus of outstanding importance to genetically characterize all individuals before any translocation actions can be approved.
Supplementary Materials: The following materials are available online at http://www.mdpi.com/1424-2818/12/9/ 343/s1, Figure S1: Pooled sampling sites of Emys trinacris. Consecutive numbers of pooled sites refer to Table S1, Figure S2: Genotypic structuring of 245 Sicilian pond turtles (Emys trinacris) from 21 sampling sites for K = 5 using 16 microsatellite loci. Samples are arranged according to the numbers of the pooled sampling sites (see Table S1). Colours of pooled sampling sites in the map correspond to GENELAND clusters; slices represent turtles with mixed ancestries or conflicting cluster assignment (percentages). Table S1: Origin of the studied Emys trinacris samples. Sample sites values refers to number showed in Figure S1.