The European pond turtle (Emys orbicularis) is an endangered Old World species. The phylogeographic history of E. orbicularis has been extensively studied throughout its range. While local genetic variation holds valuable information for conservation efforts, specific regional details have yet to be thoroughly examined everywhere. The Valencian region on the Iberian Peninsula is one such area. Here, different evolutionary lineages of E. orbicularis meet, forming a part of the western Mediterranean hybrid zone. In this study, we comprehensively sampled all Valencian localities where E. orbicularis occurs. Using mitochondrial DNA and nuclear microsatellites, we determined the genetic structure of E. orbicularis populations. Our data reveals that the mitochondrial haplotype originating from the glacial refugium in the Apennine Peninsula is primarily restricted to the northern part of the Valencian region. Additionally, human-mediated dispersal processes may have contributed to the complex relationships between the evolutionary lineages of E. orbicularis in the natural hybrid zone. On a finer scale, E. orbicularis in the Valencian region can be classified into five distinct genetic populations. Morphometric analyses revealed sexual dimorphism in the plastron shape and slight variation among genetic populations. Notably, female, but not male, plastron shape correlates with individual heterozygosity. Given the ongoing threat posed by exotic turtles, we propose implementing targeted management strategies to mitigate the presence of alien turtles. These strategies should be directed towards localities that represent each genetic population within the region, reflecting the extensive population structure observed in our study.
Introduction
Turtles are among the most endangered groups of vertebrates (Baillie et al. 2004). Among the most prominent threats to turtles are loss of habitats, alien species, hunting, climate change, and direct and indirect human-mediated loss of genetic diversity due to displacement or hybridisations (Stanford et al. 2020). Fortunately, this fact is now well recognised (Velo-Antón et al. 2021), and different conservation groups worldwide and national agencies are taking action to protect turtle habitats, mitigate population numbers decrease, and preserve the natural diversity of the turtle lineages. However, this effort is, in many places, complicated by the fact that alien turtles are still being introduced into the natural habitats of many endangered species. Hence, for conservation authorities, it might be important to know and understand not only the distribution and number of turtle individuals but also their ancestry, genetic structure, and morphological variability to decide on the conservation priorities within the region of an authority's scope (Velo-Antón et al. 2021).
Emys orbicularis (Linnaeus, 1758) is a semi-aquatic species of Emydid turtle (family Emydidae). The current distribution of E. orbicularis covers Maghreb, Portugal, the Northern Mediterranean coast, Parts of western and central Europe, the Baltic region, the Balkan Peninsula, the Black Sea region, Asia Minor, the Caspian Sea region, parts of Russia, and Central Asia (Fritz 2003). Due to the variation in climatic conditions and subsequent dynamics of the distribution of E. orbicularis, multiple evolutionary lineages evolved (Lenk et al. 1999). Some of those lineages are currently in secondary contact and naturally hybridise (Vamberger et al. 2015, Raemy et al. 2017). One such zone of secondary contact lies in the northeastern Iberian Peninsula, where three evolutionary lineages of E. orbicularis meet each other.
It is believed that three evolutionary lineages of E. orbicularis that currently occur on the Iberian Peninsula were at some point isolated in glacial refugia each at the three major Mediterranean peninsulas: the Balkan peninsula, the Apennine Peninsula, and the Iberian Peninsula (Pedall et al. 2011). Emys o. orbicularis (mt Halotype II sensu (Fritz et al. 2007)) currently spreads from the central Balkan peninsula and also occurs in the Pannonia, an isolated region in the Elbe valley in Northeast Germany, Western France but also in the North of the Iberian Peninsula throughout the foothills of Pyrenees in the valley of the River Ebro and the Catalan Coastal region. Emys o. galloitalica (mt Haplotype V) currently spreads from the western Apennine Peninsula on the Mediterranean coast of France, occurs on Corsica and Sardinia, and reaches the Mediterranean coast of the Northeastern Iberian Peninsula. Emys o. occidentalis (haplotype VI) is the lineage that occurs throughout the Iberian Peninsula and a small region in Morocco (Stuckas et al. 2014, Velo-Antón et al. 2015).
Previous research showed mitochondrial introgression, but low levels of nuclear admixture between different evolutionary lineages of E. orbicularis occur in the Iberian contact zone (Velo-Antón et al. 2015). A more detailed study (Pöschel et al. 2018) confirmed the patterns observed in the contact zone, the mitochondrial introgression and weak nuclear admixture between the three evolutionary lineages. Moreover, the results of Pöschel et al. (2018) showed that there are two regions where the hybridisation patterns are the most pronounced: the Atlantic transect from west to east of the Iberian peninsula, including the Ebro valley (between Haplotype II-Haplotype VI) and the Mediterranean coast of Iberian Peninsula (between Haplotype V-Haplotype VI). The hybrid zone along the Mediterranean coast shows comparatively more pronounced introgression of allochthonous haplotype to Iberian populations, which is very recent (Pöschel et al. 2018).
The Valencian region is a part of the Mediterranean hybrid zone, and it is considered to have an admixed origin, containing individuals of Haplotypes V and VI overall. However, in all previous studies, only four (two localities described as Valencia: Valencia: Tancat de la Pipa; and Valencia: Castelló: Plana Alta; Torreblanca: Marjales de Torreblanca (Pedall et al. 2011, Pöschel et al. 2018); and two localities described as Valencia: Sagunto, Burriana and Moro (Velo-Antón et al. 2015)) out of 11 available localities were analysed. In the Valencian region, characteristic habitats of E. orbicularis are coastal marshes that span the entire coast from south to north of the Valencia region (Bataller et al. 2008, Sancho & Ramia 2008). The system of coastal marshes was originally divided only by multiple rivers that flow from inland to the Mediterranean Sea and several mountain ranges that reach the coast. However, with the rise of human settlement and civilisation, anthropogenic barriers such as settlements and infrastructure may also act as obstacles for migrating turtles. Consequently, these barriers can impede the gene flow among turtle populations. There are also several inland localities isolated from the system of coastal marshes; however, during wetter seasons of the year, these may also be connected by the inundation system with the coastal marshes. Traditionally, Valencian E. orbicularis has been named E. o. fritzjuergenobsti, but recognition of this subspecies has been abandoned based on molecular data (Stuckas et al. 2014). The turtles in the Valencian region are reportedly smaller than in the rest of the Iberian Peninsula (Bataller et al. 2008). In the region, competitive alien turtle species are abundant, Trachemys scripta (Sancho & Lacomba 2016), and extensive measures have been taken to protect the localities of E. orbicularis. For example, extensive head-starting and other population-aid measures are taken, and the possibility of the repatriation of native E. orbicularis lineages is considered (Ayres et al. 2013). Nevertheless, until now, the detailed genetic variability of E. orbicularis in Valencia is unknown.
In this study, we have sampled all Valencian localities with the occurrence of E. orbicularis, and using mitochondrial DNA (mtDNA) and nuclear microsatellites, we determined the genetic structure of E. orbicularis populations in the region. This study primarily aims to determine if haplotype V, originating from the Apennine Peninsula, extends south of the River Ebro as part of the Mediterranean introgression zone. Additionally, the study seeks to identify any localities in the Valencia region unaffected by this hybridisation. Subsequently, we investigated if and how the population genetic structure relates to morphological variability of the shape of turtle plastron. Finally, we ask: are there any localities that would need more conservation attention than others?
Material and Methods
Sampling
Individuals of E. orbicularis (n = 204; females = 86, males = 77, juveniles = 41) were captured during regular monitoring of the populations during spring and summer in 2016 using baited and unbaited hoop traps (García-Díaz et al. 2017) (see sampling sites in Fig. 1). Small amount of blood was sampled using ultra-fine syringes from the subcarapacial vein while fixing head of the turtles retracted inside the carapace. Blood samples were immediately stored in 96% molecular grade ethanol and stored at –20 °C the same day.
Further, the standardised photographs of turtle plastron were taken using a Sony Nex digital camera in the following way: the camera was fixed on a tripod at 98 cm above the scene; the focal length was calculated using the crop factor of the camera to 78 mm to correspond to 52 mm focal length in 35 mm film; colour scale and the ruler were placed in the fix position on the scene, artificial lighting was used during the photographing. To fix the turtle position, the adult turtles were placed on the silver tape roll so they could not touch the ground with their feet or head to prevent them from successfully righting. Juveniles were placed in plastic pin boxes for the same reason. In this way, no force is imposed on the animals, and images can be standardised. The local government approved all the animal handling.
DNA extraction, PCR and sequencing
DNA was extracted using the Tissue Genomic DNA Mini Kit (GT300, Geneaid) according to the manufacturer's instructions. To determine the mitochondrial haplotype of the turtles, we used the most reported mitochondrial marker for the species, cytochrome b (cytb). To amplify the cytb, we used forward cytbG (Spinks et al. 2004) and mt-E-Rev2 primers (Stuckas et al. 2014). PCR was performed in a final volume of 25 µl. PCR conditions were the following: 420 s at 94 °C for initial incubation; 40 cycles of 40 s at 94 °C, 30 s at 50 °C, 60 s at 72 °C, and a final extension for 7 min at 72 °C. This procedure was followed by PCR product purification (Šanda et al. 2008). Sequencing was carried out by Macrogen Service Centre Europe (Amsterdam, Netherlands) using amplification primers.
To study the population structure, we used primers of eight previously characterised unlinked and polymorphic microsatellite loci (Pedall et al. 2009). After pilot experiments, the fluorescently labelled forward primers were divided into two multiplexes according to optimal annealing conditions and putative lengths of alleles (Table S1). Fragmentation analyses using PCR product, formamide, and size standard (Gene ScanTM 500 LIZ Size Standard; Applied Biosystems) were run on an ABI Prism 3100 Avant Genetic Analyzer (Applied Biosystems) in the institutional Laboratory of DNA sequencing at the Faculty of Science, Charles University.
Mitochondrial DNA analyses
From the GenBank database (Benson et al. 2012), we downloaded the sequences of E. orbicularis to compare with sequences of cytb from our samples and membership to haplotypes. All sequences were manually checked for quality and homologous regions aligned using the ClustalW algorithm implemented in Geneious 9.0.5. (Kearse et al. 2012). Haplotypes of the mitochondrial cytb were determined using DNAsp 5.10 (Librado & Rozas 2009), and the haplotype network was calculated in TCS 1.21 (Clement et al. 2000) using 95% connection limit and gap as the fifth character. The network was visualised using tcsBU (Múrias dos Santos et al. 2016). The ratios of presence of mtDNA haplotypes were interpolated into a map using R using package IDW (Pebesma 2004) as described in Brejcha et al. (2021).
Fig. 1.
Sampling of the mtDNA analyses and the haplotype structure of E. orbicularis in the Valencian region using mtDNA cytb. a) Emys orbicularis. b) Sampling sites and the interpolated ratio of haplotype V in the Valencian region. Green marks regions where haplotype V is dominant and blue where haplotype VI is dominant. 1 – coastal marsh Les Deveses in municipality Pego-Oliva, 2 – continuation of coastal marsh in municipality Pego-Oliva protruding further into Gallinero valley, 3 – Ullal de l'Estany del Duc in municipality Gandia, part of La Safor marsh, 6 – coastal marsh La Safor in municipality Xeresa, 7 – Marjal dels Moros, 8 – Almenara, 9 – Nules, 10 – Borriana, 12 – Castellón, 14 – La Pobla Torrnesa, 15 – Prat de Cabanes, 16 – Barranc de Cabanes, 17 – Vilanova d'Alcolea, 18 – marsh at Peniscola and mouth of the River Polpis in municipality Benicarlo. c) A statistical parsimony haplotype network (connection limit 95%) for 1,000 bp partial cytb sequence of studied E. orbicularis calculated using the TCS algorithm. In this network, branches represent mutations, circles represent haplotypes. Green fill marks haplotype V, blue fill marks haplotype VI, orange outline marks samples collected in this study. Haplotypes not collected in this study were downloaded from the NCBI Genbank database.
Population genetic analyses
The size of microsatellite alleles was determined using the microsatellite plugin in Geneious Prime 2020 ( https://www.geneious.com/prime), which also performed data binning. A Micro-Checker was run to check for genotyping errors like stutter bands, large allele dropouts, and null alleles (Van Oosterhout et al. 2004). To assess the genetic structure in the data, Structure 2.2. software was used (Pritchard et al. 2000) with a Model of Correlated frequencies using 1,000,000 MCMC steps with a burn-in period of 20% of total MCMC steps (200,000) and five iterations for each K. The number of tested populations was set within an interval of 1-18, corresponding to 18 localities. The result of three separate runs of Structure was evaluated in Clumpak (Kopelman et al. 2015), and the most supported K was determined using the Evanno method (Evanno et al. 2005) and median value of posterior probability Prob(K). The results were visualised using Distruct 1.1 (Rosenberg 2004). For the identified number of populations, the inbreeding coefficient was counted in Genepop (Rousset 2008) and FST between populations to assess the level of inbreeding of subpopulations relative to the total population. The number of alleles (Na) and the number of private alleles (AP) were estimated in GeneAlex (Peakall & Smouse 2006). Expected (He) and observed heterozygosity (Ho) for each population were computed in Genetix (Belkhir et al. 1999). To quantify heterozygosity for each individual, we calculated the inbreeding index using the Fbar function from the package LEA (Frichot & François 2015).
To confirm the results of Structure further, we have calculated admixture coefficients using sparse Non-Negative Matrix Factorization (snmf) implemented in R package LEA(Frichot & François 2015). We have also calculated cluster membership using Discriminant analysis of principal components (DAPC) from an R package Adegenet (Jombart 2008). The snmf procedure was performed for K 1-10 over 1,000 iterations. The optimal number of K was assessed using the cross-entropy criterion implemented in the LEA R package (Frichot & François 2015), and the run with the lowest entropy value for a particular K was selected. Further, the Principal component analysis (PCA) was performed, and subsequently, the function find. cluster was run for K1-18, which helped to identify the number of populations present in the data based on the results of PCA. After identifying the optimal number of clusters in the dataset, DAPC (Jombart 2008) was performed using individual assignment into the population. DAPC aims to assess population structure by minimising the distance between individuals a priori assigned into populations and maximising the distance between populations. The number of PC components for DAPC was assessed using function optim.a.scores in Adegenet. To evaluate whether the clustering resulting from Structure, snmf, and DAPC was congruent, we plotted group memberships for each individual from each type of analysis and compared the resulting barplots visually.
The geographic distribution of genotypes over the sampled area was assessed using the R package TESS3 (Caye et al. 2016). Similarly to PCA and DAPC, TESS3 estimation is model-free – it estimates the Q coefficient of individual ancestry and G matrix, which represents the ancestral genotypic frequencies, and it further establishes that each individual genotype is sampled from K pools of ancestral genotypes and that the sampling probabilities correspond to their admixture coefficients (Caye et al. 2016). During the spatial estimation of genotype distribution, spatial constraints ensure that individuals who are closer geographically share the ancestral genotypes more likely than those far from each other. TESS3 does not assume linkage disequilibrium or Hardy-Weinberg equilibrium and thus is prone to deviation from those assumptions introduced, for example, by population structure (Caye et al. 2016). The Q matrix for K 2-5 was plotted using the individual geographic coordinates in the sampling area. Similarly, we have used q-matrix resulting from Structure and interpolated it to geographic space using the ggtess3Q function from package TESS3. Four divisions into the population (K 2-5) were used because each K provides valuable information about the clustering hierarchy. The percentage of individuals from each of the five clusters resulting from DAPC for each population was calculated. The results were then interpolated into geographical space using the idw function from package gstat (Pebesma 2004). We have visually compared geographical projections from all three approaches: Structure, TESS3, and DAPC.
The Mantel test (Mantel 1967) was performed in GeneAlex over 999 iterations to test the presence of isolation by distance between samples. The Mantel test is commonly used to test for a relationship between geographic and genetic distances between individuals when a positive correlation indicates increasing genetic distance over increasing geographic distance. AMOVA analyses were used in GeneAlex to explore the variance in genetic variability among samples and populations (Peakall & Smouse 2006).
Geometric morphometric analyses
Data preparation for geometric morphometric analyses followed standard procedures. We defined 37 landmarks, including 18 semi-landmarks on each specimen (see Fig. 3a for landmark positions). The General Procrustes Analysis (GPA) was separately run on male, female, and juvenile shape configurations. Procrustes fit (GPA) was done using the gpagen function in the Geomorph package in R (Adams & Otárola-Castillo 2013). Procrustes-aligned configurations were projected into a tangent space. Semi-landmark positions were optimised to minimise the bending energy between corresponding points. Subsequently, we symmetrize aligned all plastron configurations so that the left and right sides were reflected, and the original and mirrored objects averaged. This procedure was performed using the symmetrize function in the Morpho package (Schlager 2017). Symmetrization was performed after semi-landmarks were moved along their tangent directions.
To test if the plastral shape variation associates with haplotype membership, symmetrised Procrustes residuals were regressed on this independent variable (ANOVA using residual randomisation) while controlling for centroid size and sex (males, females, and juveniles) as covariates. This procedure was accomplished using the procD.lm function from the Geomorph R package, with a significance test based on 9,999 permutations. Predicted shape variation was visualised by a thin-plate spline interpolation function using the plotReftoTarget function available in the Geomorph package as a deviation from the average facial shape.
The symmetrised Procrustes residuals were also regressed on microsatellite based DAPC group membership. Further, we tested for differences in plastral shape among turtles with different microsatellite DAPC group memberships using group-PCA implemented in group-PCA function, with the significance of group DAPC separation based on 10,000 permutations. We have chosen to use DAPC membership because the assignment to these groups is categorical compared to the other two methods of population assignment (Structure, snmf), while the overall patterns of population membership were similar for all three methods. Because previous tests revealed significant differences in plastron shape between males, females, and juveniles, we have analysed these groups separately.
To see whether genetic variability relates to plastron shape in turtles, we tested whether genetic heterozygosity leads to a more robust phenotype (Lerner 1954). We calculated plastron averageness to see whether the measure of heterozygosity calculated for each individual (see above) correlates with plastron shape (Kleisner et al. 2017). Based on symmetrised Procrustes residuals of plastron shape configurations, averageness was calculated as the Procrustes distance from the group mean. The greater the distance, the less average the object.
Results
The final dataset contained 140 individuals: 56 females, 53 males, and 31 juveniles. Some individuals either failed to provide a sufficient quantity of DNA, failed to amplify, or were missing photographic documentation. Please see the supplementary data files for further details ( https://doi.org/10.17605/OSF.IO/JA8MF).
MtDNA analyses
Sequences of cytb produced by Sanger sequencing displayed good quality and were approximately 1,000 bp long, depending on the quality of extracted DNA. We detected six haplotypes in our samples, of which three were not described before (see haplotype network in Fig. 1c). One haplotype, Va, belonged to the Apennine lineage. Five haplotypes belonged to the Iberian haplotype VI lineage. Two haplotypes, haplotype VIa and haplotype VIj, were previously described in the region (Velo-Antón et al. 2015). Furthermore, we demonstrated the presence of three undescribed haplotypes, each characterised by a single nucleotide mutation in addition to mutations characteristic of haplotype VIa. One was found in a single individual at Hort de Miralles, part of the Burriana locality. The second novel, haplotype VI, was found in two individuals from Marjal dels Moros. The third unrecognised haplotype was found in five individuals from southern populations (one from the locality coastal marsh Les Deveses in the municipality Pego-Oliva, two from Ullal de l'Estany del Duc in the municipality Gandia, and two from coastal marsh La Safor in the municipality Xeresa). Individuals bearing haplotype V were restricted to the northern part of the region, including two inland localities (see map in Fig. 1). The locality with the highest proportion of individuals bearing haplotype V was Vilanova d'Alcolea where there were no individuals bearing haplotype VI (out of ten captured). A high proportion of individuals bearing haplotype V (21 out of 25) was found at the coastal marshland Prat de Cabanes. Ten individuals (out of 30) corresponded to haplotype V at an inland locality, Barranc de Cabanes. Unfortunately, in the two northernmost populations, we captured only eight individuals during the whole period of sampling. Of these, five individuals corresponded to haplotype VI at marshland Peniscola, and a single individual from Desembocadura del barranco de Polpis in the municipality Benicarlo also corresponded to haplotype VI.
Fig. 2.
Results of nuclear microsatellite data analyses. a) results of admixture analyses conducted using methods implemented in Structure for K2-K5 (top to bottom). The individuals in barplots are grouped by locality from south to north. Interpolation of the ancestry coefficients into the geographic space using methods from the TESS3. b) results of admixture analyses conducted using methods implemented in LEA for K2-K5 (top to bottom). The individuals in barplots are grouped by locality from south to north. Interpolation of the ancestry coefficients into the geographic space using methods from TESS3. c) results of DAPC clustering analyses implemented in the R package Adegenet. The ordination plot shows the distances of individuals of the five clusters found. The barplot shows the cluster membership of each individual. The individuals are grouped by locality from south to north. Interpolation of the ratio of individuals belonging to each cluster for each sampled locality into the geographic was done using functions implemented in the R package gstat. The number below each barplot denotes geographic localities. 1 – coastal marsh Les Deveses in municipality Pego-Oliva, 2 – continuation of coastal marsh in municipality Pego-Oliva protruding further into Gallinero valley, 3 – Ullal de l'Estany del Duc in municipality Gandia, part of La Safor marsh, 6 – coastal marsh La Safor in municipality Xeresa, 7 – Marjal dels Moros, 8 – Almenara, 9 – Nules, 10 – Borriana, 12 – Castellón city parks, 13 – Castellón Grao 14 – La Pobla Torrnesa, 15 – Prat de Cabanes, 16 – Barranc de Cabanes, 17 – Vilanova d'Alcolea, 18 – marsh at Peniscola and mouth of the River Polpis in municipality Benicarlo.
Table 1.
Parameters of genetic variability for the five genetic populations. n = number of individuals, Na = number of alleles, AP = number of private alleles, He = expected heterozygosity, Ho = observed heterozygosity, FIS = inbreeding coefficient.
Microsatellites
Micro-Checker did not detect a genotyping error, which may result from stuttering or large allele dropout among all 140 genotyped individuals split intofivepopulationsidentifiedwithTESS3.Onelocus (GmuD55) displayed an excess of homozygotes in all five populations, probably because of missing data at the locus. All microsatellite loci were polymorphic, and the number of alleles at the loci varied from 3 to 17 when evaluating all samples (Table S2).
Fig. 3.
Results of geometric morphometric analyses. a) Landmarks and semi-landmarks. b) PCA with marked sex factor level. c) Interpolation of sexual shape differences. d) Correlation between plastron shape averageness and genetic heterozygosity for males, juveniles and females.
From the number of 1-18 tested populations with Structure, a division into two populations (K = 2) was most supported by the Evanno method (Evanno et al. 2005) and into five populations (K = 5) by the median value of posterior probability Prob(K). When divided into two populations, a population from the locality of Cabanes Interior split from the rest and partially also the southernmost population from Pego (Fig. 2). When divided into five populations, both populations from Pego and Cabanes Interior split from the rest (Fig. 2). The three remaining populations show a geographic gradient from North to South, but the level of admixture between those three populations was substantial. DAPC analyses in Adegenet support a similar division into populations as Structure did (Fig. 2), with the most supported number of clusters identified as five.
Spatial analyses from TESS3 display the following hierarchical division into populations: 1) Cabanes Interior, 2) Pego, 3) Cabanes Coast and Vilanova, 4) Peniscola, Pobla Tornesa, Gandia, Castellon, Nules, 5) Safor, Gandia, Almenara, Marjal de Moros. This division largely agrees with the non-spatial analyses provided by Structure, LEA and DAPC (Fig. 2).
Genetic variability parameters are reported for five populations described above (Table 1). The observed heterozygosity was lower than expected in all five populations. However, the results were not significant in any tested populations. The inbreeding coefficient was highest for the population grouped from localities of Nules, Borriana, Grao Castellon, Pobla cornea and Peniscola (0.175), followed by the population from Gandia, Safor, Marjal de Moros and Almenara (0.135) and population from Cabanes coast and Vilanova (0.124). On the contrary, the two most isolated populations from Cabanes Interior and Pego displayed the lowest level of inbreeding.
Fig. 4.
Differences in plastron shape among genetic clusters determined by DAPC. The ordination plot shows the first two PCA components describing the variability in the shape of the plastron. The table shows significance values for the test of differences in plastron shape for individual clusters.
The analyses in AMOVA showed that most genetic variability was found among individuals (87%), and the amount of variability that can be explained by population structure was 10%. Only 3% of the variability was explained by within-individual variability (Fig. S1). These results suggest a high level of gene flow between localities or mixing of introduced and local gene pools.
The Mantel test did not reveal a significant correlation between genetic differentiation and geographic distance between samples (P = 0.38; Fig. S2). This finding agrees with the results of AMOVA in that the level of variability differentiation among samples is not explained by geographic distance between localities. This result can be explained by the introduced population in Cabanes Interior, which resulted in bigger genetic differences between individuals from the same locality, with large differences in genetic variation over short geographic distances. Overall, the pattern suggests isolation of the population from Cabanes Interior and Pego (which may be separated from the nearby population by the River Serpis, see Discussion) but an otherwise substantial capacity for migration of Emys turtles in the rest of the sampled area.
Geometric morphometric analyses
The results of analyses of plastron shape revealed some interesting population differences. When testing for differences among turtles with different mitochondrial haplotypes, the results showed that size (P = 0.041, R2 = 0.015) and sex (P = 0.008, R2 = 0.041) (see average plastron shape of females, males and juveniles in Fig. 3c), but not haplotype (P = 0.550), significantly correlate with plastral shape. The plastron shape differences between juveniles, males and females are depicted in Fig. 3b. When testing for differences among turtles with different microsatellite-based DAPC cluster membership, the results showed that size (P < 0.001, R2 = 0.23) and sex (P < 0.001, R2 = 0.076), and DAPC membership (P = 0.022, R2 = 0.012) significantly correlate with plastral shape. Separate group PCA further for both sexes and juveniles revealed that among males, there is a significant difference in plastron shape between cluster 1 and cluster 4; among juveniles, there are significant differences in plastron shape between cluster 2 and 5, and cluster 3 and 5; and among females, there are significant differences between cluster 1 and cluster 3, cluster 2 and cluster 5, cluster 3 and cluster 5, and cluster 4 and cluster 5 (Fig. 4). See geographical interpretation of DAPC clusters in Fig. 2c. The plastron shape changes along the PCA axes and the group PCA results for individual sexes are summarised in Fig. 4.
Fig. 5.
Overview of the results in the context of the geography of the Valencian region. The thick line along the coast marks the occurrence of the two haplotypes (V, VI), the green marks haplotype V, blue haplotype VI, and grey marks regions where the occurrence of haplotypes is ambiguous based on a comparison of our results with the literature. More effort should be made to sample for the presence of haplotype V at these localities. The colours of points at sampled localities correspond to the colour coding of the population genetics results. Locality Vilanova d'Alcolea is marked by icon of roman settlement and question mark.
Overall, heterozygosity (P = 0.171) showed no significant relationship with the plastral shape averageness. However, when both sexes and juveniles were treated separately, the plastron shape of females, but not males and juveniles, showed a significant positive relationship with the individual heterozygosity (P = 0.044, resp. P = 0.893, P = 0.143) (Fig. 3d).
Discussion
Patterns of distribution of individual evolutionary lineages of E. orbicularis in the western Mediterranean are intriguing, especially south to the Pyrenees (Pedall et al. 2011, Pöschel et al. 2018). In Catalonia, haplotype II, which spreads far from the Balkan peninsula, meets with haplotype V, which spreads from the Apennine peninsula, and haplotype VI, which is the Iberian Peninsula native lineage of E. orbicularis since the last ice age. Further south along the Mediterranean coast, only Apennine haplotype V continue to co-occur at localities with Iberian haplotype VI and form a zone of introgression (Pöschel et al. 2018). In the Valencia region, admixed populations have been recorded before. However, the previous sampling was scarce and problematic. For example, one locality out of three sampled for studies before, Tancat de la Pipa (Fritz et al. 2007, Pöschel et al. 2018), composed of captive-bred individuals reintroduced during management efforts. Another example of previous poor sampling of the region is that the population studied before, ‘Burriana and Moro’, encompasses multiple localities more than 30 km apart (Velo-Antón et al. 2015). Hence, detailed knowledge of the genetic structure of E. orbicularis in the Valencia region was lacking until now.
Our data show that mitochondrial haplotype V originates from glacial refugium on the Apennine peninsula and is limited to the northern part of the Valencia region (Figs. 1, 5). However, its primary distribution appears to be inland locality, which is not the typical habitat for E. orbicularis in this area (Bataller et al. 2008). Our belief stems from the fact that this locality exclusively harbours individuals with haplotype V, a phenomenon not observed elsewhere, but despite this unique feature, the locality does not appear isolated, as evidenced by microsatellite markers indicating gene flow. Analysing the nuclear microsatellite markers, we have observed that the E. orbicularis population in the Valencia region can be categorised into two major clusters (Fig. 2). Notably, one of these clusters is associated with another single inland locality. Further differentiation in the clustering hierarchy provides evidence for subpopulations within the Valencia region (Fig. 5). None of the populations suffer from a heterozygosity deficiency or genetic diversity, and none seem to be completely isolated from the rest. Compared to most other populations, the populations that exhibit the most pronounced divergence in nuclear microsatellite markers also demonstrate distinct variations in the plastral shape of females and juveniles (Fig. 4). However, this divergence in plastral shape is not observed in males within these populations. The fine population differentiation exhibits a discernible north-to-south pattern, potentially suggesting the influence of haplotype V expansion southwards along the Mediterranean Sea coast. Nevertheless, our findings also shed light on the complexity of potential lineage admixture in the Valencia region's E. orbicularis, possibly entwined with human-mediated dispersal processes contributing to the spatial intricacies observed.
Haplotype V occurs only in the northern part of the Valencian region, while southern populations are composed solely of individuals of haplotype VI (Fig. 5). We have found the southernmost limit of haplotype V is in the coastal marsh Prat de Cabanes and a small inland stream in Barranc de Cabanes in the municipality of Cabanes. However, a previous study found individuals of haplotype V further south in the Burriana and Moro population (Velo-Antón et al. 2015). Even though we sampled localities from Marjal dels Moros, Els Estanys d'Almenara, L'Estany de Nules, until El Clot, and El Grau de Burriana, it may be possible that haplotype V is relatively rare at these localities and escaped our sampling efforts. Interestingly, the highest density of turtles with haplotype V occurs in the small stream in Barranc de Traver near the village Vilanova d'Alcolea (marked by question mark in Fig. 5), where no turtles possess Iberian haplotype VI. The locality of Prat de Cabanes is also rich in the presence of haplotype V. Haplotype V-rich localities are sometimes connected, but only during heavy rains when the periodic stream from Vilanova d'Alcolea enters the River San Miguel near the village of Les Coves de Vinromá. During most of the year, however, the putative migration route is dry, rocky, and probably also quite steep for turtles to overcome, especially in the direction from Prat de Cabanes. Thus, the migration of turtles is more probable only in one direction, from inland towards the coast. This conclusion raises the question of which locality is the source of haplotype V in the region. Local conservation authorities believe that the population in Vilanova could result from a recent introduction, possibly originating from the Cabanes marshland. However, if this were the case, the Vilanova d'Alcolea population should exhibit similar nuclear characteristics as the population from the Barranc de Cabanes site (cluster 5), which we thought to be a product of a recent introduction. Yet, the turtles at the Vilanova d'Alcolea display microsatellite characteristics that are more typical of a long-established population, and intriguingly, they might even engage in occasional interactions with the Prat de Cabanes area.
Another possibility worth mentioning is that turtles at locality Vilanova d'Alcolea were introduced in the past, and from there, turtles of haplotype V spread to at least surrounding localities. Interestingly, the locality corresponds with a recently discovered Roman settlement, the Ildum (Arasa 2008). Vilanova d'Alcolea is located on the route of the Roman road Via Augusta, which connected Hispanic colonies with Rome and, more recently, Camino Real, which was part of the road system connecting major cities from the 16th to 19th Centuries. The Ildum was a system of buildings that provided services to travellers passing through the adjacent road. Historically, turtles have been attributed with different importance (Harris 1989), and Roman and later Medieval times were no exception. Following the Greek tradition, turtles symbolised a good home and fertility (el-Kady 2011). Even though in Roman times turtle meat was only consumed exceptionally, it may still serve as food for lower social classes (Toynbee 1973). Moreover, the Roman period, along with the Neolithic, the era of the Crusades, and Colombian, is one of the four main eras of human-mediated historic biological exchange, and Romans are well known to be a source of multiple introduced species (Witcher 2013, Detry et al. 2018). Roman settlement has been implied in introducing E. orbicularis to the Balearic Islands (Braitmayer et al. 1998, Valenzuela et al. 2016). However, we mention this hypothesis here mainly to give an incentive to investigate the whole situation, and more research is needed to determine the origin of the population in Vilanova d'Alcolea and its relationship to the population in Prat de Cabanes.
On a finer scale, individuals of E. orbicularis in the Valencia region can be grouped into five populations according to clustering based on nuclear microsatellite markers (Figs. 2, 5). 1) The most southern part of the regional distribution of E. orbicularis covers coastal marshland in municipalities Pego and Oliva (Fig. 2c – green, DAPC cluster 4). This population is geographically well delimited from the rest of the populations by the River Serpis to the south, and unlike any other population, it fell within the range of the Baetic mountain system. 2) The second southernmost population has the broadest distribution in the region because it occupies the central part of the Valencia region. The population is delimited geographically from coastal marshland La Safor, i.e. northern to the River Serpis, until ancient channel El Clot de la Mare de Déu in municipality Burriana, southern to the River Mijares (Fig. 2c – blue, DAPC cluster 3). 3) In between the northern and central populations, there is population characteristic for channels and pools around Castellón de la Plana, coastal marsh in municipalities Cabanes and Torreblanca, and two inland localities at municipalities la Pobla Tornesa (in small stream Barranco de la Pobla) and Vilanova d'Alcolea. From a physical geographic point of view, this population spread north to the River Mijares and south to the River de les Coves and the Serra d'Or. Notably, this population is divided by a small (approx. 7 km long) mountain range of Serra d'Orpesa (Fig. 2c – yellow, DAPC cluster 1). 4) Inland population at Barranc de Cabanes near the village of Cabanes (Fig. 2c – red, DAPC cluster 5). 5) The most northern part of the regional distribution is formed by a characteristic population covering coastal marshes in municipalities Peniscola, Benicarlo and Vinaros. Geographically, this population is delimited by the Serra d'Irta mountains and the River San Miguel (Fig. 2c – grey, DAPC cluster 2). This structure only partially corresponds to biogeographic provinces delimited by the similarity of freshwater fish assemblages (see Fig. 3 in Filipe et al. (2009)). For example, population 1 (DAPC cluster 4) lies within an isolated region of the north-east province, population 3 (DAPC cluster 1) lies within the province Túria-Mijares, or the border between population 5 (DAPC cluster 2) and population 3, corresponding to the approximate border between the two provinces. While the borders of the provinces generally coincide with shifts in population characteristics, there are instances where the overall pattern of population structure of E. orbicularis deviates from a strict correspondence with fish-based biogeography. For example, a particular population may be encompassed by multiple provinces, or two distinct populations might fall within the boundaries of a single province. Similarly, the correspondence of the population structure of E. orbicularis in the Valencia region is only in partial agreement with the phytogeographic structure of the region (Rivas-Martínez et al. 2017). Thus, even though the population structure of E. orbicularis demonstrates certain similarities with the established biogeographical divisions of the region, factors beyond those used in the current biogeographical delineation have likely influenced the diversity of E. orbicularis in this area. Exploring these unknown factors in future studies may also shed light on shared biogeographical structuring elements for other species.
We have found overall sexual differences in plastron shape (Fig. 3b, c). From the juvenile shape, males grow their plastron by widening the humeral region, while females grow their plastron by widening the femoral region. As reproductive strategies are known to be influenced by shell shape in E. orbicularis (Zuffi et al. 2007), the relatively wider femoral region in females may relate to the general tendency in female turtles to maximise the space to lay eggs (Bonnet et al. 2001). The tendency for a wider humeral region in males could be explained in terms of exploratory and mating behaviour because it may enable males to move their forelimbs more freely. Males of E. orbicularis tend to move more than females (Lebboroni & Chelazzi 1991) and use their forelimbs to mount females during copulation (Liu et al. 2013). Thus, the plastron shape of E. orbicularis reflects different life-history characteristics of the sexes.
Further, we have found only subtle but significant population-specific differences in plastron shape (Fig. 4). Previously, differences in plastron shape among populations of E. orbicularis have been linked to habitat type (Zuffi et al. 2007). Given that the shape of the turtle shell is developmentally plastic (Zuffi et al. 2017), the notable difference observed between population 4 (DAPC cluster 5) and most other populations may result from differences in conditions turtles face during growth. Population 4 is found at a single locality at channel Barranc de Cabanes, which differs in its characteristics from coastal marshes, e.g. by water regime. Interestingly, these differences were notable only for females and juveniles, not males.
Notably, female plastron shape averageness, but not the male plastron shape averageness, correlates with individual heterozygosity (Fig. 3d). Genetic variability intrinsically relates to morphological variability. Specifically, a hypothesis exists that increased genetic heterozygosity should lead to a more robust phenotype (Lerner 1954). This hypothesis means, for example, that the individuals that are descendants of admixture between two unrelated parents should have less extreme, average, phenotype. Even though such a hypothesis seems plausible, it is difficult to draw any conclusions because the support of our results is not strong, and the literature on the topic is scarce, with ambiguous results (Vøllestad et al. 1999). The hypothesis has been challenged, e.g. by observations of the three-spined stickleback, which showed neutral heterozygosity does not influence averageness of the shape of the body but that heterozygosity in loci related to growth does (Morris et al. 2019). In E. orbicularis, neutral genetic variability is a good predictor of anomalous scutes on the carapace (Velo-Antón et al. 2011). However, if neutral genetic variability relates to the averageness of phenotype in turtles, it would be interesting to test whether these two traits are causally linked. Thus, further studies are needed to clarify the relationship between heterozygosity and phenotype.
From a conservation biology perspective, it is noteworthy that our sampling reflects the population status (and behavioural characteristics of turtles) at sampled localities. Our efforts resulted in the capture of tens of individuals at some localities with similar effort. In contrast, we struggled to obtain even a single individual at other locations, requiring substantial effort to collect samples. This difference was most apparent at localities in the northernmost part of the region, such as the municipalities of Peniscola, Benicarlo, and Vinaros. Paradoxically, the localities where individuals were difficult to find hold significant importance for further research regarding our understanding of the connection of Valencian populations to more northern regions. From a management perspective, it would be worthwhile to closely follow the population structure of E. orbicularis in the region (Fig. 5), at least until the migration of turtles within the region is better understood. This applies especially in the South of the distribution, where it is a less dynamic region, i.e. without the influence of the Apennine haplotype and with private Iberian sub-haplotypes. To better understand the population dynamics of E. orbicularis in the region, it would be beneficial for the local authorities to support more behavioural studies and continue in a mark-recapture program to see how the turtles migrate through the Valencian region. Addressing the introduction of E. orbicularis into new localities poses a complex challenge, often influenced by bioethical considerations and local politics, transcending the purview of biologists. Nonetheless, a standard general recommendation is to plan all introduction initiatives meticulously, conducting a comprehensive cost-benefit analysis encompassing multidisciplinary collaboration among academics, conservation authorities, and field experts. Considering the ongoing threat to populations posed by exotic turtles (Sancho & Lacomba 2016), we propose that, due to the extensive population structure within the region, targeted management strategies to mitigate the presence of alien turtles should be directed towards localities representative of each genetic population in the region. Therefore, in this context, the recommendation advocates for the continuation and persistence of the already established intensive efforts within the region.
In conclusion, our study sheds light on the distribution patterns and genetic structure of E. orbicularis in the Valencia region. We have uncovered intriguing patterns of co-occurrence of different haplotypes, notably haplotype V from the Apennine peninsula and haplotype VI, the native lineage of the Iberian Peninsula. Using nuclear microsatellite markers, we identified distinct populations within the Valencia region, revealing subpopulations and differentiation patterns. Our results raise the question of human influence on the distribution of genetic lineages of E. orbicularis. Additionally, our findings highlight sexual differences in shell shape, which may be related to heterozygosity and raise intriguing questions about the interplay between genetics and morphology in this species. Further research that includes comparative analyses with populations from the Ebro region and other areas will be valuable. Integrating advanced genetic techniques and broader sampling efforts will contribute to a more comprehensive understanding of the evolutionary history and population dynamics of E. orbicularis in the Valencia region and beyond. These future investigations hold the potential to uncover additional insights into the species' origin, migration patterns, and adaptive traits, enhancing conservation efforts and informing management strategies for this unique and ecologically important turtle species.
Acknowledgements
J. Brejcha and K. Eliášová dedicate this study to their sons Jindra and Vilém. The work of J. Moravec was financially supported by the Ministry of Culture of the Czech Republic (DKRVO 2019-2023/6.VII.e, National Museum of the Czech Republic, Prague, 00023272). The study was part of the project ‘Demonstration strategy and techniques for the eradication of invasive freshwater turtles’ (LIFE+Trachemys LIFE09 NAT/ES/000529). The work of J. Moravec was financially supported by the Ministry of Culture of the Czech Republic (DKRVO 2024-2028/6.I.a, National Museum of the Czech Republic, 00023272).
This is an open access article under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits use, distribution and reproduction in any medium provided the original work is properly cited.
Author Contributions
J. Brejcha, J. Moravec – conceived the study. J. Brejcha, K. Eliášová – designed the sampling, collected the data, did laboratory work, performed data analyses, wrote the manuscript. J.V. Bataller – designed and managed the fieldwork, collected the data and edited the manuscript. K. Kleisner – performed geometric morphometric analyses, contributed to drafting the manuscript. J.I.L. Andueza, V.S. Alcayde – assisted with the fieldwork design, edited the manuscript. L. Čtrnáctová, P. Nevečeřalová – contributed to analyses. J. Moravec – provided financial support to the project, contributed to sampling design, wrote the manuscript.
Data Availability Statement
The data supporting the findings of this study are available at the Open Science Framework (OSF) repository: https://doi.org/10.17605/OSF.IO/JA8MF.
Literature
Appendices
Supplementary online material
Table S1. Fluorescently labelled forward primers in two multiplexes according to optimal annealing conditions and putative lengths of alleles.
Table S2. Number of alleles at the microsatellite loci when evaluating all samples analysed.
Fig. S1. AMOVA analysis of the relative distribution of genetic variability. Most of the genetic variability was found among individuals (87%). Population structure explained 10% of the total variability. Only 3% of the variability was explained by within-individual variability.
Fig. S2. Correlation between genetic differentiation (GGD) and geographic distance (LinGD) between samples.
( https://www.ivb.cz/wp-content/uploads/JVB-vol.-73-2024-Brejcha-et-al.-Tables-S1-S2-Figs.-S1-S2.pdf)
