Spotlight on islands: on the origin and diversification of an ancient lineage of the Italian wall lizard Podarcis siculus in the western Pontine Islands.

Groups of proximate continental islands may conceal more tangled phylogeographic patterns than oceanic archipelagos as a consequence of repeated sea level changes, which allow populations to experience gene flow during periods of low sea level stands and isolation by vicariant mechanisms during periods of high sea level stands. Here, we describe for the first time an ancient and diverging lineage of the Italian wall lizard Podarcis siculus from the western Pontine Islands. We used nuclear and mitochondrial DNA sequences of 156 individuals with the aim of unraveling their phylogenetic position, while microsatellite loci were used to test several a priori insular biogeographic models of migration with empirical data. Our results suggest that the western Pontine populations colonized the islands early during their Pliocene volcanic formation, while populations from the eastern Pontine Islands seem to have been introduced recently. The inter-island genetic makeup indicates an important role of historical migration, probably due to glacial land bridges connecting islands followed by a recent vicariant mechanism of isolation. Moreover, the most supported migration model predicted higher gene flow among islands which are geographically arranged in parallel. Considering the threatened status of small insular endemic populations, we suggest this new evolutionarily independent unit be given priority in conservation efforts.

Introduction

Islands epitomize a suite of simplified conditions in spatially bounded areas, which make harbored populations an ideal model to investigate several evolutionary topics[1-3]. For example, oceanic islands contributed a lot to our empirical understanding of a variety of microevolutionary processes that drive diversification and eventually speciation[4-7]. In such a context, island size, physiography, proximity among islands and dispersal capabilities are expected to shape the genetic structure of a pioneer species that colonizes an oceanic archipelago through oversea dispersal. Conversely, islands situated on continental shelves have often been reconnected with the mainland after sea level drops during glacial phases[8]. This scenario creates allopatric conditions that change through time, enabling both vicariant and dispersal mechanisms to act on genetic diversity and divergence[9,10]. In addition, groups of small geographically proximate islands may comprise an even more complex phylogeographic setting compared to single continental islands because nearby islands may experience repeated cycles of isolation (during sea level transgression) and connection (during periods of low sea level stands). Indeed, a complex geological setting dominated by considerable changes of islands size might lay the foundation for a variety of different genetic outcomes as a result of changing demographic processes and recurrent periods of fragmentation and gene flow[11-13].

With almost 5000 islands, the Mediterranean Basin provides an attractive geological history which helps to explain the insular biodiversity observed today[14]. During the Pleistocene, the basin was characterized by cyclic sea level drops, which caused the formation of transient land bridges connecting nearby islands to each other, as well as islands and the mainland. The Mediterranean insular biota contains numerous taxa formed in response to such eustatic fluctuations[15-17].

Almost half of the 23 recognized species of Podarcis lizards are endemic to Mediterranean islands and they represent an ideal model to explore the evolution of genetic diversity on islands. Here, we focus on the Italian wall lizard Podarcis siculus which, aside from inhabiting the entire Italian Peninsula, is also widespread across the greater islands of Sicily, Corsica, and Sardinia and across all the main Italian archipelagos and smaller islands (Fig. 1a). The phylogeography of this species across the Italian Peninsula is very complex[18] and a recent study[19] revealed that long term processes, giving rise to a mosaic of allopatric lineages, combined with more recent demographic dynamics, have driven the evolutionary history of this species. Here, we examine insular genetic variation of P. siculus from the Pontine Islands using mitochondrial sequences (mtDNA), nuclear sequences (nuDNA), and microsatellite loci. The Pontine Archipelago is formed by two groups of islands. The western group, which lies 32 km from the Tyrrhenian Coast, includes the oldest island of Ponza and the nearby Gavi Islet, Zannone Island, and Palmarola Island. These islands are separated from the eastern Pontine Islands by 40 km of open sea; this eastern group of islands includes Ventotene and the Santo Stefano islands, located 50 km away from the mainland (Fig. 1b). On the basis of morphological traits, different subspecies have been described within this archipelago: Podarcis s. latastei (Bedriaga, 1879) on Ponza Island, Podarcis s. patrizii (Lanza, 1952) on Zannone Island, Podarcis s. lanzai (Mertens, 1952) on Gavi Islet, and Podarcis s. palmarolae (Mertens, 1967) on Palmarola Island.

Figure 1

(a) Map of the study area showing all sampling sites on the Italian Peninsula with the geographic distribution range of Podarcis siculus highlighted in red. (b) Map of the Pontine Archipelago with relative sampling locations.

This particular insular context provides ample opportunity to investigate several factors regarding the origin and evolution of intraspecific genetic diversity within a Mediterranean archipelago. First, the volcanic origin of the islands has been thoroughly investigated and relatively precise ages for their emergence have been reported[20,21]. Presumably, they never experienced any stable contact with the mainland. Second, the surrounding shallow waters indicate that the islands have been subjected to repeated episodes of fragmentation during high sea level stands as well as reunification when sea levels were lowest, during glacial phases. Finally, the western Pontine Islands have a spatial arrangement that allows testing different hypotheses of island biogeography, specifically regarding their compatibility with the observed patterns of genetic divergence and differentiation in this Mediterranean archipelago.

In this study we aim to (i) unveil the evolutionary history of Podarcis siculus from the Pontine Islands, and (ii) assess the influence of different scenarios of inter-island genetic diversification and gene flow by comparing models of island biogeography with genetic empirical data.

Results

Phylogenetic reconstruction and time of divergence

In the concatenated mtDNA alignment of 1701 bp (nd4 = 766 bp and cytb = 935 bp;156 sequences; accession numbers are reported in Supplementary Table S4), we found 100 haplotypes with 368 polymorphic sites and an overall mean distance of 0.062 among the main clades. Uncorrected pairwise p-distances are reported in Supplementary Table S5.

The phylogenetic inference returned the same topology under the two Bayesian approaches applied, revealing the presence of eight clades (Fig. 2a and Supplementary Fig. S1), seven of which had already been described in previous works[18,19]. However, a new diverging lineage (clade P in Fig. 2a,b) was discovered to which all individuals from the western Pontine Islands belong. The phylogenetic reconstruction indicated an early split separating the S, C1, and C2 clades from all the other clades (A1, A2, T1, T2, and P). Subsequently, a second cladogenetic event involved the separation of the Pontine clade P.

Figure 2

(a) Chronogram based on mtDNA performed with beast. At each node, the posterior probabilities and the coalescence time with their relative 95% high posterior densities (HPD) are shown. Chronology of the geological epochs with the two principal volcanic episodes that formed the western Pontine Islands shown below. (b) Geographic distribution of the mtDNA clades at each sampled locality. (c) Statistical parsimony networks of the two nuclear gene fragments; circle sizes are proportional to allele frequencies. Alleles are colored based on the mtDNA clade assignment of the respective specimen.

The divergence time analysis carried out using beast showed an effective sample size over 200 for all parameters, indicating good chain mixing and adequate posterior estimates. The coefficient of variation in the average mutation rate estimates suggests a substantial heterogeneity among clades, supporting a non-clock-like model for the phylogeny (ucld.stdev parameter was 0.33 for cytb and 0.47 for nd4). According to our divergence time estimate, the Pontine clade (P) separated about 3.7 Mya (95% HPD: 5.5–2.2 Mya, Fig. 2a).

The final nuDNA alignment consisted of 252 phased sequences (600 bp) for mc1r and 164 sequences (728 bp) for β-fibint7 (accession numbers are reported in Supplementary Table S4). A test for recombination using the Pairwise Homoplasy Index (phi) did not find significant evidence for recombination in either of the two genes (p-value = 0.772 for mc1r; p-value = 0.559 for β-fibint7). The nuDNA did not show a phylogenetic pattern fully congruent with the mtDNA, but a certain level of accordance can be detected for the β-fibint7 fragments (Fig. 2c). In those specimens assigned to the new mtDNA clade P, all but one (F41) β-fibint7 allele form a monophyletic group. Moreover, all but three β-fibint7 alleles (F6, F20 and F28) were found only in specimens belonging to the same mtDNA clade. The mc1r alleles were typically shared among specimens of different mtDNA lineages. However, most alleles found in individuals of the western Pontine Islands are closely related at this locus and are mostly restricted to these islands, while they are very rare or absent elsewhere.

Among island structure and gene flow

The Podarcis siculus populations of the western Pontine Islands were characterized by 31 mtDNA haplotypes with 44 polymorphic sites. The number of haplotypes (H) and nucleotide (π) and haplotype (h) diversity are reported for each dataset and marker in Table 1.

Table 1

Genetic diversity estimates for each Podarcis siculus insular population (PO: Ponza, GA: Gavi, PA: Palmarola, ZA: Zannone, FM: Faraglioni della Madonna).

island	mtDNA					mc1r					β-fibint7
	n	p.s.	H	h + sd	π + sd	n	p.s.	H	h + sd	π + sd	n	p.s.	H	h + sd	π + sd
PO	18	11	10	0.86 ± 0.06	0.001 ± 0.0002	36	10	6	0.69 ± 0.04	0.004 ± 0.0007	22	14	4	0.72 ± 0.06	0.006 ± 0.001
FM	5	6	4	0.90 ± 0.10	0.001 ± 0.0005	8	7	4	0.78 ± 0.11	0.004 ± 0.0012	6	12	3	0.73 ± 0.15	0.009 ± 0.002
GA	20	4	5	0.44 ± 0.13	0.0003 ± 0.0001	40	9	6	0.69 ± 0.04	0.004 ± 0.0003	14	1	2	0.52 ± 0.06	0.0007 ± 0.0001
PA	10	11	8	0.95 ± 0.05	0.001 ± 0.0002	24	9	7	0.79 ± 0.04	0.005 ± 0.0004	14	3	4	0.62 ± 0.11	0.001 ± 0.0008
ZA	8	9	7	0.96 ± 0.07	0.001 ± 0.0003	14	8	5	0.72 ± 0.08	0.005 ± 0.0006	12	2	3	0.54 ± 0.14	0.0008 ± 0.0002
Tot	61	43	31	0.87 ± 0.03	0.002 ± 0.0002	122	11	12	0.79 ± 0.02	0.005 ± 0.0002	68	15	6	0.68 ± 0.04	0.0040 ± 0.0009

n; number of individuals/alleles, p.s.; number of polymorphic sites, H; number of haplotypes, h; haplotype diversity, π; nucleotide diversity, sd; standard deviation.

The mtDNA parsimony network showed the presence of two main genetic assemblages (Fig. 3a), one of which is restricted to the Zannone Island, where all specimens carried haplotypes of that assemblage. This result was also confirmed by the DAPC analysis, which suggested K = 2 as the most reliable number of genetic clusters (Supplementary Fig. S2).

Figure 3

(a) Map and bathymetry of the western Pontine Islands (PO: Ponza, GA: Gavi, PA: Palmarola, ZA: Zannone, FM: Faraglioni della Madonna). (b) Geographic location of the western Pontine Islands with the direction of the main marine currents in the Tyrrhenian Sea according to El-Geziry and Bryden[78]. (c) Statistical parsimony networks of the mtDNA (cytb + nd4) and nuDNA (mc1r and β-fibint7) fragments. Circles sizes are proportional to the haplotype frequencies and filled rectangles representing one mutation step. Haplotypes are colored by islands with different shades of green (for color code, see (c).

The nuDNA parsimony networks did not identify a concordant pattern for either mitochondrial or island repartition (Fig. 3a). The mc1r fragment showed the presence of 12 alleles and most of them were found in individuals across all island populations. Notably, of the 6 β-fibint7 alleles found, one allele (F4) was separated by 11 substitutions from the other 5 alleles. The DAPC carried out on nuclear sequences did not allow the identification of an optimal K value, returning a flat BIC curve (Supplementary Fig. S2).

All the 61 P. siculus specimens were genotyped at 11 polymorphic microsatellite loci, but the loci Pli 21 and Lv4a were excluded from further analyses because the former was monomorphic, while the latter might be influenced by positive selection. Genetic characteristics of the microsatellite loci are reported in Supplementary Table S6. Significant departures from the Hardy-Weinberg equilibrium (HWE) were observed in 2 out of 45 locus/population exact tests (Bonferroni correction applied: p < 0.001, see Supplementary Table S7). No linkage disequilibrium was detected between any pair of loci. Evidence of null alleles was found for several loci, but none was consistent across all populations, therefore we did not exclude them from further analysis. An indication of recent reduction in population size (bottlenecks) was obtained only for the Faraglioni della Madonna population by both the Wilcoxon test (P = 0.00684) and the mode-shift indicator, which revealed a distortion in the typical L-shape distribution of allele frequencies.

Both clustering analyses conducted on microsatellite data (DAPC and structure approaches) assigned K = 2 as the most likely number of genetic clusters among islands (Supplementary Figs S2 and S3). The first genetic cluster includes only the population from Gavi, while the second cluster includes all the remaining populations. Despite a further subdivision scheme, as suggested by the deltaK when K = 4, the individual assignment at the fourth cluster did not correspond to any specific island, suggesting a reliable genetic partition not higher than two (Supplementary Fig. S3). The presence of further differentiation among all P. siculus populations is suggested by FST analyses, where all comparisons yielded significant values (P < 0.005) (Bonferroni correction applied; Table 2).

Table 2

FST values (below) and respective p-values (above) for each pair of populations (PO: Ponza, GA: Gavi, PA: Palmarola, ZA: Zannone, FM: Faraglioni della Madonna).

	ZA	GA	PO	FM	PA
ZA	—	<0.00001	0.0048	0.0020	0.0010
GA	0.1144	—	<0.00001	<0.00001	<0.00001
PO	0.0466	0.0655	—	<0.00001	0.0010
FM	0.1859	0.1595	0.1081	—	0.0010
PA	0.0710	0.0807	0.0394	0.1347	—

Significant differences after Bonferroni correction (p < 0.005) are emboldened.

The analysis of current gene flow performed with BayesAss indicated no evidence of recent gene flow in any pairwise comparisons. The migration rate estimates were all below 0.07 with the exception of the migration estimate between the Zannone and Palmarola populations. This was detected above the background noise (M = 0.16). The Migrate-n analysis carried out on the microsatellite dataset under a full migration matrix model indicated the presence of asymmetric gene flow between island pairs, with higher migration rates detected between the island of Ponza and all the others (Fig. 4). The Bayes factor calculating the model probability using the log marginal likelihood of each generated migration scenario, indicated that the mtDNA migration model was the most reliable scenario for gene exchange.

Figure 4

Models of migration among western Pontine Islands tested through Bayes factor analysis (PO: Ponza, GA: Gavi, PA: Palmarola, ZA: Zannone). Inferred migration rates are only shown for the most supported model.

Discussion

Dispersal plays a key role for the biodiversity of oceanic islands[5,6]. Conversely, studies on continental islands provided more complex evolutionary histories dominated by a variety of microevolutionary processes, such as vicariance, dispersal, and gene flow upon secondary contact[11,15,22-24]. Although the Pontine Archipelago could be referred as a particular case of oceanic islands system in the biogeographic sense of the term, it is located on continental shelf where changes in island size induced by sea level fluctuations may have been drastic. Such a special context provides a compelling setting to investigate the birth and the formation of intraspecific genetic diversity within archipelagos. Our results outline a tangled evolutionary history triggered by an ancient dispersal event of P. siculus from the mainland to islands, followed by local extinctions and intra-islands vicariant processes with gene flow.

On the origin of island genetic divergence

Phylogenetic reconstruction using two mitochondrial markers confirmed the complex phylogeographic pattern across the Italian Peninsula, with seven very divergent clades whose evolutionary histories have been described in detail in previous works[18,19]. In this study, we detect the presence of a new divergent 8th P. siculus clade corresponding to populations inhabiting the western Pontine Islands of Ponza, Gavi, Palmarola, and Zannone. On the contrary, the island of Ventotene of the eastern Pontine Island group showed haplotypes belonging to the A2 clade, indicating a recent introduction of P. siculus to this island, as already reported[25].

According to our estimate of divergence time, the ancestry of the Pontine clade dates back to the early/late Pliocene transition (3.7 Mya; 95% HPD: 5.5–2.2 Mya). Although we used the same methodology, our estimates were slightly different to those obtained in a previous work[19], but comparable when the HPD for each node is taken into account. Such a discrepancy could be due to the incorporation of an additional marker (nd4) highlighting how this analysis can be influenced by the reliability of the mutation rates. However, the present reconstruction is in line with the geological evolution of the western Pontine Islands (Fig. 2a). Indeed, on the basis of potassium-argon dating, the first volcanic event that accompanied the formation of the Ponza Island occurred between 4.2 and 3.7 Mya, followed by a quiescent period up to 3.2 Mya[21]. Although the estimated dates are associated with large uncertainty values of the coalescent process, our molecular clock suggests that the western Pontine populations colonized the archipelago early after their Pliocene volcanic formation. Furthermore, global changes in ice volume and records of terrestrial and marine sedimentary sequences suggest four glaciation events during the Pliocene, with one occurring about 3.6 Mya during the early/late Pliocene transition[26]. A consequence of global cooling is the sea-level retreat, which likely relaxed the separation of Ponza Island from the mainland facilitating a possible dispersal event. A similar explanation has also proposed to explain other ancient colonization events in the Mediterranean Basin[15,27]. It must be emphasized though, that the origin of these islands only puts an upper limit to the time since colonization and the age of the respective evolutionary lineages, as there is evidence that other vertebrate species colonized these islands much later[28].

The genetic diversity of P. siculus from the Pontine Islands has been previously investigated using allozymes[29]. These authors were unable to detect any genetic differentiation of the western Pontine Island populations, neither with respect to the eastern Pontine Islands nor in comparison with the mainland populations. They subsequently suggested a possible human-mediated introduction in historical times. Our results from the nuDNA dataset clearly reject such a hypothesis, as allele frequencies of western Pontine Islands were significantly different from the other populations and several alleles of both mc1r and β-fibint7 fragments were private either for the western Pontine Islands or even for single islands. The population structure defined by β-fibint7 was more similar to the mtDNA, as all island alleles except one (F41) appear to form a monophyletic assemblage and no allele is shared with the mainland. On the contrary, most alleles of the mc1r gene occurred in specimens belonging to several mitochondrial clades. Such discrepancy is not surprising if we bear in mind that coding genes, such as mc1r, have lower mutation rates and are thus more prone to exhibit incomplete lineage sorting[30,31]. In addition, a recent study which reconstructed the phylogenetic relationship of 934 sequences of mc1r in nine Podarcis species revealed only slight differentiation, with some alleles shared even among different species[32], corroborating our interpretation of incomplete lineage sorting.

Inter-island genetic diversification

On the basis of historical sea level changes from the Early Pleistocene, the Pontine Islands have probably been subjected to repeated cycles of connection (into a single paleo-island setting) and isolation[33]. Our inter-island diversification reconstruction is in line with such a scenario only if we consider a Middle Pleistocene divergence followed by different extent of genetic exchange during subsequent periods. The absence of an earlier divergence episode is interesting considering the ancestry of these populations and the fragmentation of the islands assumed for the Middle Pleistocene periods due to sea level elevation. At the same time, it is also difficult to imagine that migration processes have been maintained for such an extended period (from 3.7 to 0.2 Mya) to prevent further divergence. A more plausible explanation is that populations diverged in allopatry, i.e., when islands were disconnected. However, some of these ancestral island populations may have become extinct during strong marine transgressions, such that their present populations originate from a more recent episode of divergence. Such an interpretation is also supported by stratigraphic and paleontological evidence, which suggest that at least two principal episodes of strong marine transgression on the Tyrrhenian Coast occurred some 0.54 and 0.35 Mya ago, which prompted a sea level rise from 55 up to 120 m a.s.l.[34]. In addition, beach deposits from the Early Pleistocene are present at an elevation of 200–250 m a.s.l. on Monte Guardia (Ponza Island), while erosional surfaces from the Middle Pleistocene are widespread at 100–120 m a.s.l[35]. Therefore, considering the low altitude of some Pontine Islands (maximum altitude: Gavi = 60 m; Zannone = 170 m; Palmarola = 180 m; Ponza = 250 m), some of them might have been totally or partially submerged. Subsequent re-colonization from nearby islands during marine regressions may then have homogenized the genetic make-up among island populations.

According to our chronogram, the first episode of divergence within the Pontine Islands involved the separation of the Zannone Island population during the Mindel–Riss interglacial period (MIS 11, Fig. 2a). This period has been described as one of the longest warm periods of the Quaternary, inducing marked footprints of global sea level rise[36]. Under such vicariant mechanisms of divergence, taking into account similar distances and bathymetric depths among islands, we should expect also a similar extent of among islands genetic diversification followed by genetic admixture as a consequence of glacial inter-island reunification. Indeed, sea depth among islands never exceeds 100 m (Fig. 3b), which is almost the same magnitude of sea level drop during the last glacial maximum[37,38].

Our multilocus dataset revealed a pattern of mito-nuclear discordance within the western Pontine Islands. The mtDNA showed the presence of two principal assemblages: one composed of individuals from Zannone Island and the other including individuals from Ponza, Gavi, and Palmarola islands. However, the nuDNA did not support strong population structure, neither when relating it to geography (different islands) nor to the mtDNA groups. Finally, the clustering analysis carried out with both DAPC and structure on the microsatellite dataset identified a K = 2 subdivision, revealing the Gavi Islet as a separate cluster. The lack of accordance between genes and genealogies represents a common issue, but can provide relevant insights into microevolutionary processes underlying the origin of small island system biodiversity[39]. In such an insular context, the causes of mito-nuclear discordance should be primarily a result of the different evolutionary trajectories of distinct markers. Indeed, the mtDNA has an inherently higher mutation rate and, due to its uniparental inheritance, is expected to reach fixation fourfold faster than the nuDNA genome because of its smaller effective population size and the associated higher susceptibility to demographic changes and genetic drift. On the contrary, the lack of congruence of the nuDNA dataset to both geographic settings and mtDNA differentiation is probably because genetic differentiation on these islands is rather recent, such that the sorting of the nuclear genealogies is not complete. Such a scenario has been repeatedly reported in many vertebrate species, including other Podarcis lizards[30,31].

Although a principal trend of mito-nuclear discordance emerged when considering the mere genealogies and the current structure of the markers, our coalescent-based analysis of gene flow allowed us to partially explain this complexity. Indeed, the most supported migration scenario using microsatellite data, was the mtDNA migration model. This model suggests reduced gene flow between the Ponza and Zannone islands, and higher gene flow between the Ponza and Palmarola islands (Fig. 3b). Such a prediction is particularly plausible because these latter two islands have a parallel geographic arrangement. Therefore, the genetic exchange could have been even more pervasive during a low sea level stand, when islands shared a larger contact area. However, we cannot rule out that the reduced gene flow between the Ponza and Zannone populations is the outcome of demographic dynamics like “density blocking”[40,41], or alternatively, is due to morphological or behavioral changes which occurred during the allopatric divergence.

Taxonomic and conservation implications

The genetic differentiation we found in the populations from the western Pontine Islands with respect to the mainland is comparable with those observed for many other Podarcis lizards species and over twice as high as between the insular endemic P. raffoneae and its sister taxon P. waglerianus[24,42,43]. On the other hand, populations from the eastern Pontine Islands seem to have colonized the Ventotene Island in recent times[25]. Considering the greater distance from the mainland and the deeper water in which the eastern Pontine Islands are located, the absence of a comparable ancient lineage is fairly surprising. Interestingly, in 1926 the herpetologist Robert Mertens described a particular phenotype endemic to Santo Stefano Island (P. siculus sanctistephani), which then became extinct and was replaced by the ordinary mainland phenotype during the first decades of the last century[44,45]. Such observations deserve special attention and may have several implications. First, it could be argued that an ancient lineage was also historically present on the eastern Pontine Islands, including Ventotene. Secondly, and perhaps more fascinating, it would be interesting to understand why the P. siculus sanctistephani became extinct. Currently, there is no trace of introgression in the current allochthonous population from the Santo Stefano Island[25,29]. In addition, a recent work focused on morphological variation of 374 individuals from the Pontine Islands and mainland, including the extinct population of Santo Stefano island, did not find intermediate phenotypes as could be expected during a complete admixture process[46]. If the local endemic P. siculus sanctistephani became extinct after the new colonizers and was indeed not hybridizing with later arriving conspecifics of the mainland type, there might have been a reproductive isolation mechanism, e.g. a behavioral or morphological change present in these island populations.

All these considerations point to a priority of the conservation of all the western Pontine Island populations, which should be referred to a different taxon. Indeed, the aforementioned example of possible population replacement underscores the vulnerability of endemic insular populations, suggesting that further introductions from the mainland could be deleterious and should be prevented in order to maintain the genetic integrity on this archipelago. In support of our conclusions, additional morphological, behavioral, and physiological studies would be needed to unravel whether the found genetic divergence is associated with ecological and/or phenotypic diversification on these islands.

Materials and Methods

Laboratory techniques and molecular data

A total of 156 specimens of P. siculus have been analyzed in this study: 63 new tail tissues were collected from six islands (Ponza, Faraglioni della Madonna, Gavi, Zannone, Palmarola, and Ventotene) while 93 individuals from 46 localities of the Italian Peninsula have been selected from a previous study[19] (see Supplementary Table S1). The tail tissues were collected in accordance with relevant guidelines in full compliance with specific permits released by the Italian Environment Ministry (Prot. 00017879/PNM-09/09/2012) and no lizards were killed.

Genomic DNA was extracted according to the universal protocol of Salah[47]. We amplified one portion of the NADH dehydrogenase subunit 4 (nd4) with flanking regions (tRNASer, tRNAHis and tRNALeu) in all 156 collected samples. Additionally, three fragments including the cytochrome b (cytb), the coding melanocortin receptor 1 (mc1r), and the β-fibrinogen intron 7 (β-fibint7) were amplified only in the 61 specimens of the Pontine Islands, as the samples from the Italian Peninsula had been previously analyzed at these loci by Senczuk et al.[19]. Details on reaction buffer, primers, and cycling conditions are reported in Supplementary Table S2. The obtained amplified products were purified with Sure Clean (Bioline) and sequenced with an automated DNA sequencer by Macrogen (www.macrogen.com).

BioEdit 7.2[48] was used to evaluate the electropherograms, encode heterozygote positions using IUPAC ambiguity codes, and compute consensus alignments. The two mitochondrial fragments (mtDNA) were concatenated. For each nuclear gene (nuDNAs), we inferred the gametic phase using the coalescent-based Bayesian algorithm in phase 2.1[49,50]. The β-fibint7 showed insertion/deletion (INDEL) polymorphisms in a region of poly-T and poly-G starting at 157 bp. Following the method described by Flot[51], this region was used to determine the phase for sequences that were heterozygous for INDELs. The known phases were then implemented in the coalescent-based Bayesian reconstruction. Three independent runs were conducted with a burn-in at 1 × 103, 1 × 104 iterations, and thinning at each 100 steps. Finally, the occurrence of recombination events was assessed for each nuclear gene through the Pairwise Homoplasy Index (phi) test implemented in splitstree v. 4[52].

The 61 samples from the western Pontine Islands were also genotyped at 11 microsatellite loci. Details on primer and PCR conditions are reported in Supplementary Table S3. Microsatellite multiplex amplifications were performed and then compared with singleplex PCR to confirm that there were no apparent allele dropouts or allele size differences. Amplification products were run on an ABI 3130xl Genetic Analyzer and allele size was scored using peak scanner 1.0 (Applied Biosystems). Amplification and genotyping were repeated twice for some 15% of the total samples (10 lizards) to verify reproducibility in microsatellite scoring. Genepop 3.4[53] was used to assess deviations from Hardy-Weinberg equilibrium and linkage disequilibrium. The statistical significance levels were corrected for multiple comparisons using Bonferroni’s method (α = 0.05). lositan[54] was used to identify candidate loci under selection. Before subsequent analyses, the presence of null alleles, large allele dropouts, and scoring errors was checked using micro-checker[55].

To infer phylogenetic relationships, the mtDNA dataset was analyzed using MrBayes v. 3.1.2[56]. The best substitution model for each partition (i.e., HKY + I + G for nd4 and GTR + I + G for cytb) was selected using JModeltest[57] under the Akaike Information Criterion (AIC). The analysis was run for 3 × 106 generations using four annealing chains, sampling every 1 × 104 and discarding the first 10% of computed trees as burn-in. Sequences of Podarics muralis retrieved from GenBank (accession numbers: HQ65293 for cytb, KF372413 for nd4) were used as an outgroup.

The software beast 1.8.1[58] was used to date the most recent common ancestors (TMRCA). In this analysis the selection of proper priors is a crucial component although it may be challenging. Because of the absence of fossil records in Podarcis lizards to calibrate the phylogeny, we used substitution rates previously estimated for both gene fragments[12,43]. We incorporated a lognormal prior distribution to the mean rate in order to constrain the 95% of the prior distribution between the upper and lower values estimated for Podarcis (μ = 0.0175, SD = 0.2 for cytb, and μ = 0.0115, SD = 0.5 for nd4). Another relevant issue is the consideration of how substitution rates can vary across branches assuming a clock-like or alternatively a relaxed-clock model. To resolve this issue, a preliminary run was performed using a lognormal relaxed-clock in order to account for branch heterogeneity by screening the ucld.stdev parameter, which provides information about whether evolutionary rates inferred from the data are clock-like. Finally, the choice of an outgroup can be difficult when the phylogeny is not fully resolved[32]. However, the use of a relaxed-clock allows estimates of the tree root even in the absence of a known outgroup[59]. beast was run for 1 × 108 generations sampling every 1 × 104 steps and the software tracer 1.6[60] was used to check parameter convergence. The final consensus tree was summarized from the stationary distribution (burn-in = 25%) in TreeAnnotator 1.8.1[58].

Furthermore, to better understand the relationship between mtDNA and nuDNA and the variation of the latter throughout the whole species distribution area, two nuDNA parsimony networks were built using the software popart[61]. For this analysis we removed 19 bp from the original β-fibint7 dataset (157–175), due to the presence of many INDELs polymorphisms in a region of poly-T and poly-G. Finally, uncorrected pairwise p-distances were calculated for each mtDNA clade obtained by the phylogenetic reconstruction using mega 7.0[62].

To infer the genetic relationship between mitochondrial and nuclear haplotypes among island populations, we used the mtDNA and nuDNA dataset to construct statistical parsimony networks under the 95% probability connection limit using the software popart.

For the mtDNA and nuDNA datasets, the number of haplotypes (H) as well as nucleotide (π) and haplotype (h) diversity were estimated using DnaSP 5.1[63]. The microsatellite genetic variability was estimated by calculating allele frequencies, expected (HE) and observed (HO) heterozygosities, number of alleles, and number of private alleles using Genetix 4.05[64,65] and Fstat 2.9.3.2[66]. All calculations were performed for the whole dataset and for each island. The software Bottleneck 1.2.02[67] was used to assess the occurrence of recent bottlenecks; the tests were performed assuming the two-phase model (TPM).

To explore the genetic structure among island populations, we performed a discriminant analysis of principal component (DAPC) implemented in the R-package ‘Adegenet’[68] using the mtDNA, nuDNA, and microsatellite datasets. For the microsatellite data we also performed a Bayesian clustering analysis implemented in structure 2.3.1[69] to infer the number of different genetic clusters (K), assuming the admixture model. Ten simulations of 1 × 105 iterations for each K (from K = 1 to K = 7) were run with a burn-in of 1 × 106. The best K value, based on the modal value ΔK[70], was detected using the online software Structure Harvester 0.6.93[71]. In addition, Arlequin 3.5[72] was used to investigate the interpopulation level of genetic differentiation by calculating F, according to Weir and Cockerham[73].

In order to assess the presence of gene flow and whether such transfer was recent or historical we employed two different approaches. The software BayesAss 3.0[74] was used to estimate the presence of gene flow within the past few generations. Several preliminary runs were performed to adjust the acceptance rate of the mixing parameters (migration rates, allele frequencies, and inbreeding coefficients). The final analysis was carried out with 5 independent replicates running for 1 × 107 iterations, sampling every 100 steps and discarding the first 1 million generations as burn-in. Convergence was assessed by comparing the accordance of the posterior mean parameter estimates across independent runs and by checking the stationary of the trace file using the software tracer. To obtain an estimate of a more ancient signal of gene flow, we used the software migrate-n[75,76]. The program implements a coalescent-based approach to estimate the mutation scaled effective population size (θ = 4Neu) and the mutation scaled long-term migration rate (M = m/u, where m is the migration rate, and u is the mutation rate). We generated a set of plausible a priori biogeographic migration models, to test microsatellite empirical data: (i) the “free gene flow model”, which implies gene exchange among all islands, (ii) gene flow among all islands except for Zannone, as indicated by our mtDNA results (“mtDNA model”), (iii) gene flow in a south-east/north-west direction as could be supposed when considering the marine currents, favoring overseas dispersal (“marine currents model”), and (iv) gene flow from the bigger island to the smaller (“big to small island model”). Bayes factors were calculated for each model and used to test the more reliable scenario following the guideline of Kass and Raftery[77]. The Faraglioni della Madonna population was excluded from this analysis because of the small sample size (n < 5). After a preliminary tuning run to estimate θ and m parameters through FST to use as starting values, we performed four independent runs for each generated model using four parallel chains with static heating (temperatures 1.0, 1.5, 3.0 100,000.0) and Brownian motion microsatellite model.