Multilocus phylogeography of Italian Moorish geckos adds insights into the evolutionary history of European populations

Geckos of clade III of the Tarentola mauritanica species complex are widespread throughout southern Europe and northern Africa. We investigated the genetic variability of the Italian populations by performing a widespread sampling throughout the mainland and the two main islands of Sicily and Sardinia. We analysed 199 newly generated sequences of the mitochondrial 16S rRNA gene and 269 nuclear genotypes inferred from nine microsatellite loci from 307 individuals. We found 13 new mitochondrial haplotypes in Italy, whereas previous findings reported a single haplotype widespread throughout the country and in the rest of Europe, which currently make Italy the centre of genetic diversity of this taxon. There was no evidence of mitochondrial DNA structuring with geographic correlation. At the population genetic level, our multilocus approach based on nuclear markers returned a shallow genetic structure. Nonetheless, we disclosed the presence of at least four distinct genetic clusters (namely the Adriatic, two Tyrrhenian and the Calabrian clusters). Our findings do not support the two hypotheses proposed to explain the low level of mitochondrial polymorphism in this taxon, namely the genetic hitch‐hiking due to selective sweep and the historical human‐mediated colonization hypotheses. Based on the fossil record, the presence in Italy of this taxon since the Pleistocene Epoch is plausible. Given the discordance in genetic structure between mitochondrial and nuclear DNA, the exact role of the Italian Peninsula in shaping the observed patterns of genetic diversity during the Pleistocenic climatic oscillations needs further investigation.

The six clades composing the T. mauritanica complex, including the species T. angustimentalis from the Canary Islands, had been identified based on strong reciprocal divergence at the mitochondrial DNA (mtDNA) (Harris et al., 2004a(Harris et al., , 2004b)).Subsequently, the integration of nuclear data suggested that these clades represent unconfirmed candidate species, mostly distributed in north-western Africa (Rato et al., 2010(Rato et al., , 2012(Rato et al., , 2016)).Accordingly, this region has been identified as the centre of origin of the complex, whose diversification from the rest of the genus Tarentola was dated in the mid-late Miocene (ca.8.69 Ma), a period coinciding with the rise of the current Atlas Mountains.These six clades subsequently diversified during the paleogeographic events linked to the Messinian Salinity Crisis (ca.5.96-5.33Ma), when the level of the Mediterranean Sea decreased sharply and rapidly due to the interruption of the connection with the Atlantic Ocean (Rato et al., 2012).The reopening of the Strait of Gibraltar and the climatic and sea level oscillations during the Pleistocene further influenced the evolutionary history of the complex (Rato et al., 2015).As a result of this, the phylogenetic and phylogeographic patterns that we observe today in these clades are likely the result of a complex combination of allopatric and adaptive speciation events (Rato et al., 2015).
Among the six lineages of the T. mauritanica species complex, clade III received much scientific attention (see Rato et al., 2010Rato et al., , 2012Rato et al., , 2016)).This lineage, indeed, has the widest distribution, inhabiting Morocco, Algeria, Tunisia, the Iberian Peninsula, southern France, Italy, the Balkan Peninsula and Greece.Nonetheless, it is characterized by very low levels of diversity at the mtDNA, particularly in Europe, where only one haplotype (H2) is currently known from the entire distributional range (Rato et al., 2010).It was initially proposed that the European populations resulted from an anthropogenic introduction from northern Africa in historical times (Harris et al., 2004a(Harris et al., , 2004b)).However, the discovery of clades II and VI in the Iberian Peninsula, along with the detection of a higher level of diversity in the nuclear DNA (nuDNA), led to the hypothesis that a mechanism of genetic hitch-hiking of the mtDNA due to selective sweep took place (Perera & Harris, 2008;Rato et al., 2010Rato et al., , 2012Rato et al., , 2013)).In this process, a neutral allele can strongly increase its frequency in a population, theoretically reaching complete fixation, due to the linkage with another allele under positive selection (Gillespie, 2000;Losos, 2014).
Genetic data indicate that within the Italian Peninsula populations of T. mauritanica belong to clade III (Rato et al., 2012;Sindaco et al., 2006).Individuals of this lineage are found in areas with Mediterranean climate and especially in coastal regions, and most occurrences are reported between 0 and 400 m a.s.l.However, in recent years a growing number of records of T. mauritanica outside its native range (i.e., northern Italy) have been reported (Sindaco et al., 2006).Populations of clade III also occur in the islands of Sardinia and Sicily, whereas clades belonging to the T. fascicularis/deserti complex are known from Lampedusa (clade VIII) and the close Conigli islands (clades VIII and IX), which are located between Sicily and Tunisia (Harris et al., 2009).
Like the other southern European peninsulas, Italy played a major role as a climatic refugium during the alternating Pleistocene glacial events (Hewitt, 2000).Multiple evidence from amphibians and reptiles have shown that when retreating within the Italian Peninsula during glacial periods, populations remained geographically isolated in smaller refugia.This process, commonly referred to as the 'Refugia within refugia' paradigm (Gomez & Lunt, 2007), shaped both the inter-and intra-specific diversity of the amphibian and reptile species inhabiting Italy (e.g., Canestrelli et al., 2014;Kindler et al., 2013;Rato et al., 2009;Salvi et al., 2013;Senczuk et al., 2017).
In the present study, we increase the sampling of T. mauritanica from Italy to investigate the population genetic structure at a finer scale.We use a multilocus approach to unveil the phylogeographic and diversity patterns of the Italian populations and contribute shedding light on the evolutionary history of the entire clade.

| Sampling
We sampled 306 individuals from 51 different localities (Figure 1; Tables S1, S2) within the Italian native range of the species between 2012 and 2015.Animals were spotted during nocturnal and diurnal opportunistic searches and caught by hand and/or noosing.After capture, the extremity of the tail was clipped and preserved in 96% ethanol.No animals were sacrificed to accomplish this study.
One museum specimen from the Natural History Museum of Genova (Italy) was included in this study (ABZC01273).This individual was originally collected from Lampedusa Island (locality 39) in 2001.The final data set includes 307 individuals (Tables S1, S2).

| DNA extraction, amplification and sequencing
Total genomic DNA was extracted following a standard proteinase K and saline solution extraction protocol (Bruford et al., 1992).Polymerase Chain Reaction (PCR) amplifications were carried out for a fragment of ca.450 bp of the mitochondrial 16S rRNA gene (16S) using available primers (Palumbi et al., 1991).We amplified this molecular marker since it has been used in several previous phylogeographic studies of the genus, and the reference database available on GenBank is particularly large.We amplified to five individuals for each locality.We tried to amplify all the individuals of each locality only when more than one haplotype was detected within the first five sequenced samples.Exceptions to this procedure were localities 7, 8 and 22, 23 as they are only a few kilometres apart and the five individuals were counted over the two pairs of areas (Table S1).We also amplified nine microsatellite loci of nuclear DNA (Arranz et al., 2013) for F I G U R E 1 Haplotype phylogenetic network inferred with TCS on the 319 sequences of the 16S gene (comprising both sequences newly generated in this study and GenBank accession numbers) and geographic distribution of clade III of the Tarentola mauritanica species complex.Each circle in the map represents a locality and each colour a specific haplotype.Numbers in the map refer to the sampled localities and follow the numeration in Table S1.Squares represent accession numbers from Italy retrieved from GenBank.The proportions of different colours within squares and circles represent the actual haplotype frequencies within each locality.Madeira's archipelago is represented in the square on the low-left part of the figure.Sample ABZC01273 from Lampedusa Island (locality 39) belonging to Clade VIII of the T. fascicularis/deserti species complex (Figure S1) is shown in yellow.In the haplotype phylogenetic network, dots represent unsampled or extinct haplotypes and segments between dots are estimated mutational steps between haplotypes.a subset of 28 localities which had at least five sampled individuals, in order to build a sufficiently informative dataset to perform genetic structure analyses (Tables S1,  S2; primers and amplified loci are listed in Table S3; see Text S1 for details on laboratory work and sample amplification).16S chromatograms were checked by eye in BioEdit 7.2.6 (Hall, 1999).All sequences amplified for this study have been submitted to GenBank (accession numbers OR940307-OR940505; Table S2).Microsatellite alleles were determined using GENEMAPPER 6.0 (Applied Biosystems) and checked manually.

| Phylogenetic analysis
We performed a phylogenetic analysis to assign our samples to the different clades of the T. mauritanica and T. fascicularis/deserti species complexes (Rato et al., 2012).The 16S sequences generated in this study were complemented with 16S accession numbers retrieved from GenBank.Sequences of Tarentola boehmei Joger, 1984, Tarentola boettgeri Steindachner, 1891 and T. deserti were used as outgroups (Tables S2, S4; File S1).Sequences were aligned with the MUSCLE algorithm (Edgar, 2004a(Edgar, , 2004b) ) and collapsed into unique haplotypes using the online web tool 'DnaCollapser and converter' available on the FaBox 1.5 online platform (Villesen, 2007).We performed phylogenetic analyses under a Bayesian Inference model using MrBayes 3.2.6 (Ronquist et al., 2012) and visualized the 50% majority rule consensus tree in FigTree 1.4.3 (Rambaut, 2009) (see Text S1 for details on the analysis).

| Haplotypes phylogenetic network
We estimated a phylogenetic network of mitochondrial haplotypes to characterize the intra-specific relationships among sampled localities using the newly generated 16S sequences and other sequences of the same marker retrieved from GenBank for clade III (Tables S2, S4; File S2), following the results of the previous phylogenetic analysis.Sequences were aligned with the ClustalW algorithm implemented in BioEdit (Thompson et al., 1994).The phylogenetic network was computed with TCS 1.21 (Clement et al., 2000), which implements a parsimony criterion as defined in Templeton et al. (1992).We set a 95% cut-off probability and treatment of gaps as the fifth state to estimate the number of mutational steps by which pairs of haplotypes differ.The final network was obtained with the online tool tcsBU (dos Santos et al., 2015), which improves the visualization of TCS outputs and makes the interpretation easier.We tested the nine microsatellite loci for deviations from Hardy-Weinberg and linkage disequilibrium using the complete data set of all the genotyped individuals (n = 269) from 28 localities (Table S1; File S3).We used the software GENEPOP 4.7 (Raymond & Rousset, 1995;Rousset, 2008) with a probability test of 1000 generations of the Markov chain Monte Carlo algorithm and 95% confidence interval.The software MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004) was used to statistically test for the presence of genotyping errors with 1000 generations of the MCMC algorithm and 95% of confidence interval.
The following population genetic analyses were carried out on a reduced dataset with higher statistical power, including individuals with <10% of missing data and using localities with at least 10 sampled individuals.The final reduced data set comprised eight localities, spanning most of the geographic extension of peninsular Italy and 93 individuals (Table S1; File S4).

| Genetic structure
The genetic structure of the Italian population of T. mauritanica clade III was first explored with a principal component analysis (PCA) on the reduced microsatellite data set of 93 individuals and eight localities (Table S1), using the packages adegenet 2.1.1 and ade4 1.7-22 implemented in R 4.2.3 (Dray & Dufour, 2007;Jombart, 2008;R Core Team, 2023).
We investigated the genetic structure of the sampled populations with the Bayesian Clustering algorithm implemented in Tess 2.3.1 (Chen et al., 2007).This software accounts for both geographic and genetic data of the samples and can outperform similar software when the genetic structure is not particularly deep (Chen et al., 2007).We tested genetic clusters (K) from 2 to 15, running 50,000 sweeps (20,000 discarded as burn-in), 100 replicates for each K value, conditional autoregressive (CAR) admixture model, spatial interaction parameter set to the default value (0.6) with the option of being updated.The K value that best fitted the data set was chosen based on the deviance information criterion (DIC) associated with each run.The K values were plotted against the DIC values averaged over the 100 replicates for each K.When the average DIC decreased its slope towards a plateau, the associated K value was considered as the most suitable value.For the selected K value, we averaged the membership coefficient matrices of the 20% of the runs with the lowest DIC values with CLUMPP 1.1.2(Jakobsson & Rosenberg, 2007).The final membership coefficient matrix was displayed with a bar plot generated in DISTRUCT 1.1 (Rosenberg, 2004).We finally generated a posterior predictive map with Tess representing the genetic clusters spatially interpolated on the geographic map of Italy that was visualized with TESS Ad-Mixer 1.0 (Mitchell et al., 2013).

| Genetic diversity
The following indexes of genetic diversity were inferred on the reduced microsatellite dataset (eight localities; 93 samples) at both intra-populational and inter-populational levels: observed and expected heterozygosity, average number of alleles, average allelic range, number of polymorphic loci and genetic differentiation (F ST ).All the calculations were performed with Arlequin 3.5.2.2 (Excoffier & Lischer, 2010) with 16,000 random permutations of the significance test at 95% of confidence interval.To better visualize F ST values, we performed a multidimensional scaling plot using the function 'cmdscale' implemented in R.

| Mitochondrial phylogenetic tree and haplotype network
The final data set of the 16S gene used in the phylogenetic analysis included 469 sequences divided between 199 sequences generated in this study (Table S2) and 270 GenBank accession numbers (Table S4; File S1).After excluding hypervariable and misaligned portions (Text S1), the final alignment comprised 482 characters.The best model of sequence evolution under all Information Criteria (AIC, AICc, BIC and DT) was GTR + G.We recovered eight clades in total, five belonging to the T. mauritanica species complex (clades I to VI) and three to the T. fascicularis/deserti species complex (clades VII, VIII, IX) (Figure S1).
Inter-clade relationships were overall supported.Clades VIII and IX of the T. fascicularis/deserti species complex were fully supported (PP = 1.00) and belonged to a clade which also contains clade VII (PP = 1.00).The T. fascicularis/deserti species complex was recovered in sister position to the T. mauritanica species complex (PP = 0.91).All but clade IV of the T. mauritanica species complex were fully supported (PP = 1.00).Clades I and IV were recovered as sister (PP = 0.93).Clade III was sister to clades I and IV (PP = 0.90), and the group composed of these three clades was in sister position to clade II (PP = 1.00).Finally, clade VI was recovered as sister (PP = 1.00) to the group composed of all previous clades of the T. mauritanica species complex.All samples of Italian Moorish gecko newly generated in this study belonged to clade III of the T. mauritanica species complex except the individual ABZC01273 from Lampedusa Island, which belonged to clade VIII of the T. fascicularis/deserti species complex, known to occur in Libya, Lampedusa and Conigli islands (Figure S1).This sample was consequently excluded from the following analyses.
The haplotype phylogenetic network of 16S included all individuals analysed in this study of the clade III of T. mauritanica (n = 197) and 122 GenBank accession numbers of the same clade (Tables S2, S4; File S2).Sequences of one sampled individual (ABZC01654) from locality 28 and three accession numbers (FJ156045, FJ156046 and KY762013) were not included in the haplotype phylogenetic network analyses because they were shorter than the rest of the alignment.However, based on the phylogenetic analysis, they belonged to clade III (Figure S1).
The analysis recovered 18 haplotypes (Figure 1; Tables S1, S2 and S4).The haplotype phylogenetic network was characterized by a moderate star-like shape (Figure 1).Some haplotypes were closely related to the central haplotype (H2) either via a single mutational step or through short branches, whereas others showed some articulated branching among each other and relative to the central H2.The common haplotype H2 had high frequency (91.8% relative to the analysed individuals) and occurred across the entire distributional range of clade III on both European and African sides (Figure 1).Outside Italy, this was the only haplotype known from the entire European range.Within Italy, this haplotype was present in all sampled localities except at Timpazzi (Messina, Sicily; locality 30) and Torre Stelle (Cagliari, Sardinia; locality 34), which, anyway, were composed of only two and one individuals, respectively.We found four haplotypes (H15, H16, H17 and H18) exclusive from northern Africa.All of them were present in Morocco, and only H16 occurred in Tunisia as well.13 haplotypes are currently known only from Italy (H1, H3-H14), which were never found in previous studies and showed relatively low frequencies, never above 1.5% relative to all the included sequences.They are often represented by a single individual, and they are normally restricted to a single geographic area (exceptions to this are the haplotypes H5 and H11).The areas of Sicily (localities 28-32: haplotypes H2, H11, H12 and H13), Liguria (locality 1: haplotypes H1, H2, H3 and H4) and Sardinia (localities 34-38 and from 40 to 51; haplotypes H2, H5, H11, H14) harboured higher mtDNA diversity.

| Microsatellite genetic structure and diversity
In the whole nuclear data set (n = 269; File S3), one single locus (MT13) resulted in Hardy-Weinberg disequilibrium (p-value: .0001).Regardless, we decided to keep this locus in the following population genetic analyses as eight loci do not represent a large number of markers and may not provide enough information to recover a reliable genetic structure, especially when considering that the loci were not particularly informative as shown by the average number of alleles (see Table 1).We did not detect linkage disequilibrium for any pair of loci, nor potential genotyping errors.

| Genetic structure
The PCA recovered an eastern Adriatic coast group, and a western Tyrrhenian coast group (Figure 2).Localities 6 (Marche), 18 and 20 (Apulia) were part of the same Adriatic group, while localities 1 (Liguria), 2 (Tuscany), 8 and 9 (Latium) belonged to the western Tyrrhenian coast group.The Calabrian locality 25 was positioned between the two groups showing a possible genetic affinity with both.In the Axis 1 × Axis 3 plot (Figure 2b), the western Tyrrhenian coast group was more expanded than in Axis 1 × Axis 2 plot (Figure 2a), showing some affinity between localities 8 and 9, while localities 1 and 2 grouped more singularly, slightly separated from each other and the rest of the localities.
In Tess analysis, according to the plot of K values against DIC values, the most probable number of clusters was five (Figure 3b).We can interpret both member coefficients plot (Figure 3c) and the posterior predictive map (Figure 3a) as showing four actual different genetic clusters, three of which are on the Tyrrhenian Sea coast and one in the Adriatic Sea coast.This subdivision mostly mirrored PCA results (Figure 2).Localities 1 and 2 formed a northern Tyrrhenian cluster (red), with the first (light purple, Figure 3c) showing strong admix with the second.
F I G U R E 2 Principal component analysis (PCA) performed with the R packages adegenet and ade4 on 93 individuals from eight Italian localities that were genotyped for nine microsatellite loci (File S4).Note that coloured circles highlight the five genetic clusters that were found in the genetic structure analyses performed with Tess (Figure 3c).Both Axis 1 × Axis 2 (a) and Axis 1 × Axis 3 (b) plots are shown.The first axis accounts for 14.24% of the total variation, the second for 10.12% and the third for 8.68%.Axis 1 divides the sampled localities into an eastern Adriatic coast group (on the left) and a western Tyrrhenian coast group (on the right).
T A B L E 1 Standard genetic diversity indexes calculated with Arlequin on the dataset of eight Italian localities, 93 individuals and nine microsatellite loci used for genetic structure analyses of the clade III of the Tarentola mauritanica species complex.Localities numeration follows Table S1.The second Tyrrhenian cluster (green) was positioned in the central portion of the coast with localities 8 and 9.The last Tyrrhenian cluster (light blue) was located at the extreme south (locality 25) with some admixture and affinity to the green central Tyrrhenian cluster.Finally, our results showed the presence of a geographically extended Adriatic cluster (yellow) composed of three investigated localities ranging from the northern-central part of the Adriatic coast (locality 6), passing through locality 18 until the extreme south in Apulia (locality 20).Average member coefficient values for each dominant genetic cluster were 0.65, 0.88, 0.89, 0.92 and 0.73 for light purple, red, yellow, green and light blue clusters, respectively.In order to provide further evidence of the analysis results, we report alternative scenarios for K = 3 and K = 8 in Supporting Information (Figures S2, S3), both recovering the same four-group structure of the scenario with K = 5.

| Genetic diversity
In the nuclear-reduced data set, the average number of alleles per locus were not particularly high, ranging from 2.11 to 2.56, with an overall average of 2.28 for the eight localities (Table 1).The average number of polymorphic loci was 8.25 among the eight localities and the average allelic range 9.16.Observed and expected heterozygosity were moderate and, in many cases, similar to each other within the same locality, ranging from 0.23 to 0.51 and from 0.31 to 0.45, respectively.F ST values calculated between the same localities were all statistically significant and overall high, with only six comparisons lower than 0.20, confirming a relatively high differentiation between localities (Table S5).The highest F ST value was between localities 6-9 (0.42), 1-6 (0.38) and 1-20 (0.37), whereas the lowest values were between localities 6-18 (0.05), 18-20 (0.13) and 18-25 (0.11).When considering the genetic structure results, in a few cases localities belonging to different genetic clusters had more similar F ST than localities from the same group.For instance, localities 20 and 2, which belonged to different genetic clusters, had F ST values of 0.17, whereas the same locality 20 showed a slightly higher F ST (0.18) with locality 6 from the same Adriatic cluster.This pattern is better summarized in the multidimensional scaling plot computed on the same F ST values (Figure 4).Although the genetic clusters highlighted in Tess analysis were recovered, in some cases the distance between localities did not totally reflect that genetic structure, which could provide additional information on the reciprocal affinities among sampled localities.This was the case for localities 6 and 18, which showed less affinity with locality 20, belonging to the same Adriatic cluster, than locality 25 from Calabria belonging to the southern Tyrrhenian cluster.Also, locality 1 resulted quite far from locality 2, despite the affinity shown in the member coefficient plot (Figure 3c), the latter being closer to the central (green coloured) Tyrrhenian cluster.previous findings (Harris et al., 2009;Rato et al., 2012), although none of the past studies included such a dense and extended sampling from Italy.There is another clade in Italy, known only from Conigli Island (clade IX of the T. fascicularis/deserti species complex; Joger & Bshaenia, 2010;Rato et al., 2012), but no samples from this islet were included in the present study.

| DISCUSSION
We found 13 new mitochondrial haplotypes within clade III, which, to date, are known only from Italy and make the Italian Peninsula the centre of mitochondrial haplotype diversity of this taxon, followed by Morocco (five total haplotypes) and Tunisia (two haplotypes).Haplotype H2, which has been so far the only known European haplotype (Rato et al., 2010), has by far the highest frequency and the largest distribution in Italy (Figure 1; Tables S1,  S2 and S4).
The widespread presence of the H2 haplotype in almost all populations indicates the absence of a genetic structure with geographic correlation at the 16S level.The interpretation of the geographic distribution of haplotypes within the Italian territory provides contrasting signals: in some cases, closely related haplotypes are distributed in contiguous areas (e.g., haplotypes H8, H9, H12 and H13) that span over a long geographical continuum between locality 10 (Latium), in central Italy, and the several localities that were sampled in Sicily around Messina (localities 28-32) (Figure 1).In other cases, the same haplotype can occur at hundreds of kilometres of distance (e.g., H5, H11).This might be explained by the synanthropic nature of T. mauritanica, which can facilitate its translocation to different areas.
Our study suggests that a denser sampling from each locality is needed to uncover the genetic variability of clade III.Most of the genetic data available on GenBank originates from the sequencing of one or two sampled individuals at each site.Given the high frequency of the haplotype H2 within localities and its ubiquitous presence, the sampling of multiple individuals from the same locality is suggested to identify new haplotypes.
In contrast to mtDNA, microsatellite data revealed a different scenario and identified a higher genetic diversity at the nuclear level, as highlighted in previous studies on clade III (Rato et al., 2010(Rato et al., , 2013)).PCA, F ST and Tess analyses were overall reciprocally concordant in suggesting the presence of a genetic structure with geographic correlation.However, the few dissimilarities among these analyses might indicate that the detected structure is not particularly deep (Figures 2-4).This might be because the analysed loci are not particularly informative, as exemplified by the average number of alleles per locus (Table 1).A widespread cluster was recovered throughout the Adriatic Sea coast, on the eastern part of the Italian Peninsula (yellow cluster in Figure 3a) opposed to the higher differentiation characterizing the Tyrrhenian Sea coast.At least three genetic clusters can be recognized here: a southern (light blue) cluster, a central (green) cluster and a northern (red) cluster.The latter was the least consistent among the three analyses.In both PCA and F ST analyses, individuals from localities 1 and 2 showed much less reciprocal affinity than in the member coefficient bar plot (Figure 3c).The southern Tyrrhenian cluster (light blue) showed in all analyses a clear affinity with both the Adriatic cluster (yellow) and the rest of the Tyrrhenian clusters (green and red).Similar to mtDNA analyses, the genetic structure uncovered with microsatellites highlighted possible signs of occasional translocations that can be interpreted from the admixture of individuals among different genetic clusters as, for instance, individual ABZC02568 belonging to locality 9 (Latina-Borgo Montello, Latina) clearly admixed with the red northern Tyrrhenian cluster (Figure 3c).
The overall genetic structure suggests a predominant role of the Apennine Mountain chain as a potential barrier to the species gene flow.This agrees with the Ecological Niche Model study from Rato et al. (2015), showing the low suitability of montane environments for clade III in the Italian Peninsula.The affinity of the southern Tyrrhenian cluster with both the Adriatic and the other Tyrrhenian groups is probably explained by the intermediate position of Calabria in terms of connectivity between the two coasts.The importance of the Apennines as a geographical barrier to gene flow triggering allopatric differentiation, especially between Adriatic and Tyrrhenian coasts, has also been reported for other species of Italian reptiles, although this mountain range did not generally act as a barrier across its entire latitudinal extension, as suggested by our results.Among these, we report the examples of Hierophis viridiflavus (Rato et al., 2009), Natrix natrix (Kindler et al., 2013) and Podarcis muralis (Salvi et al., 2013), in which multiple glacial refugia were hypothesized, some of which in the central and northern parts of the Peninsula.Specifically, the genetic pattern we found, characterized by a unique Adriatic cluster and multiple Tyrrhenian groups, has been previously detected by Senczuk et al. (2017) for Podarcis siculus, who defined the concept of 'western richness and eastern purity'.The higher climatic stability of the Tyrrhenian coast during the Quaternary climatic fluctuations favoured the persistence and differentiation of isolated populations of P. siculus in several coastal refugia, as opposed to the higher genetic homogeneity of the Adriatic coast.Following these considerations, the Italian Moorish gecko may have followed a similar scenario of multiple refugia within Italy.Although the probable presence in Italy of this gecko during the Pleistocene and the consequent influence that the climatic oscillations of that epoch may have had in shaping the phylogeographic patterns of the species (see below), any hypothesis on the location and number of glacial refugia of T. mauritanica within the Italian Peninsula would be only tentative.Moreover, the absence of an analogous genetic structure with geographic correlation at the 16S currently prevents us from hypothesizing any plausible scenario, although other Mediterranean reptile groups show a mismatch between genetic structure and geographic distribution due to complex dynamics of range shifts during the Pleistocene climatic oscillations that caused repeated isolation of populations followed by expansions resulting in secondary contacts (Rato et al., 2021;Velo-Anton et al., 2018).

| Implications on the evolutionary history of the European population of clade III
At the time of this work, two different hypotheses had been proposed to explain the origin of the European population of clade III of the T. mauritanica species complex and the presence of a unique haplotype (H2) in the overall European range.The first hypothesis suggested an anthropogenic introduction in historical times (Harris et al., 2004a(Harris et al., , 2004b)).However, the following discovery of two additional clades of the species complex in the Iberian Peninsula and the detection of a much higher genetic diversity at the nuclear level led to the hypothesis of genetic hitch-hiking of mtDNA caused by selective sweep (Gillespie, 2000;Rato et al., 2010Rato et al., , 2012Rato et al., , 2013)).
According to the fossil record, remains assigned to the genus Tarentola and T. mauritanica, or morphologically closely related taxa, are known from Spain since the Early Pleistocene (Gelasian, ca.2.58-1.80Ma and Calabrian, ca.1.80-0.78Ma, respectively), attesting the ancient occurrence of this taxon in Europe (see Villa & Delfino, 2019 and references therein).Younger Pleistocenic remains of T. mauritanica are known from France and England, while in Italy remains of this taxon from Sicily (San Vito lo Capo) dating to the Late Pleistocene-Holocene pose the possibility that the species was present in southern Italy since at least the Upper Pleistocene (Villa & Delfino, 2019).In light of this, our molecular results can provide some insights into the two alternative hypotheses explained above (Harris et al., 2004a(Harris et al., , 2004b;;Rato et al., 2010Rato et al., , 2012Rato et al., , 2013)).
Following the theory of genetic hitch-hiking from selective sweep, an allele can increase its frequency within a population, potentially reaching complete fixation.Although our analyses were not aimed at specifically testing this hypothesis, which should be better addressed with different methodological approaches, the discovery of 13 new haplotypes is in disagreement with the occurrence of this process.However, it is true that until the process reaches total completion, other haplotypes can still occur (partial sweep) while gradually decreasing their frequency (Gillespie, 2000;Losos, 2014).
On the other hand, the discovery of new mitochondrial haplotypes, the detection of a genetic structure at the microsatellite level and the presence of Tarentola sp. and T. mauritanica remains in the fossil record of southern Italy and other European countries (Villa & Delfino, 2019) do not support an anthropogenic introduction in historical times (human-mediated colonization hypothesis).On one side, the processes involved in determining the observed genetic patterns are expected to require much more time to take place (Koskinen et al., 2002;Richard & Thorpe, 2001).In addition, the shallow but significant population structuring detected along the peninsula with nuclear markers supports a scenario of diversification taking place in situ during the last glacial cycles.
We analysed the 16S gene due to the large reference database available for the T. mauritanica species complex to identify the clades to which our samples belonged and compare our findings with the previous studies on this group (Rato et al., 2010(Rato et al., , 2012(Rato et al., , 2013)).However, the low variability of this marker within the Italian population suggests that the same is unsuitable to investigate phylogeographic patterns.To test the validity of our phylogeographic hypothesis, we suggest that future studies should consider other mitochondrial loci ideally more suitable for this scope, as the cytochrome b gene, which has been shown by Rato et al. (2013) to be among the most variable markers in the mitogenome of this taxon.Alternatively, cost-effective genomic techniques such as RADseq targeting hundreds of loci (Andrews et al., 2016) could uncover more deeply the genetic structure of the Moorish gecko populations and its geographic correlation.Similarly, the investigation of the phylogeography of clade III in other southern European peninsulas, such as Iberia and the Balkans (Hewitt, 2000), might help clarifying the origin of the European populations of this clade.Finally, we suggest increasing the sampling size at each new Italian locality, and especially in the diverse Tyrrhenian coast and Calabria for its focal position in the connection between the eastern and western coast sides of the peninsula.

AUTHOR CONTRIBUTIONS
AB and AC conceived the ideas and designed methodology; AB, FB, DPR, WC, CL and CR performed research; FB and AB analysed data; FB and AB led the writing of the manuscript.All authors contributed critically to the drafts and gave final approval for the publication.

F
Genetic structure analyses performed with Tess on the data set composed of 93 individuals, eight Italian localities and nine microsatellite loci.(a) Posterior predictive map representation of the genetic clusters interpolated on the geographic map of Italy for the K = 5 scenario visualized with TESS Ad-Mixer.The distribution of the different cluster colours is especially reliable in the areas surrounding the sampled localities.(b) Number of K clusters plotted against the DIC values averaged over the 100 runs for each K value.(c) Member coefficients bar plot obtained with CLUMPP and DISTRUCT for the K = 5 scenario showing the subdivision of the individuals of the eight localities with colours matching the posterior predictive map.

4. 1 |
Phylogenetic history of the Italian populations of T. mauritanica The recovered genetic and geographic subdivision of the Italian populations between clades III of T. mauritanica and VIII of T. fascicularis/deserti complexes agrees with F I G U R E 4 Multidimensional scaling plot performed with the function 'cmdscale' in R on the F ST values on the dataset composed of 93 individuals, eight Italian localities and nine microsatellite loci computed with Arlequin.Colours are the same as the ones used for the genetic clusters inferred with Tess for the scenario K = 5 (Figure 3c).