Sampling, laboratory procedures and molecular data

The sampling took place from March 2013 to June 2015, during which we collected 277 samples from 115 localities covering most of the distribution area of Podarcis siculus in the Italian Peninsula and Sicily. The tissues were obtained by inducing autotomy after light pressure and were then stored in 96% pure ethanol. All lizards were released at the capture location. The sampling was planned to complement already available data in order to achieve a better geographic coverage of the species distribution [30]. The geographic coordinates were recorded and high-resolution photos were taken for each individual. The geographic references and sample size of all sampled populations are given in Additional file 1: Table S1 and shown in Fig. 1.

Map of the study area showing the principal geographic features mentioned in the main text. The area in red corresponds to the geographic distribution of Podarcis siculus while the 115 new sampled locations are shown in black dots

Genomic DNA was extracted from all the tissue samples by means of the universal extraction protocol of Salah et al. [31] with incubation at 56 °C with proteinase K and DNA precipitation with isopropanol. The total genome was visually checked in electrophoretic runs. We used standard PCR protocols to amplify the mitochondrial fragment of cytochrome b (cytb) and two nuclear fragments, the Melanocortin receptor 1 (mc1r) and the β-fibrinogen intron 7 (β-fibint7). The primers used for the amplifications and their respective references are reported in Table 1. As the presence of a pseudogene for cytb has been reported in this species [32], we modified the primers to be strongly specific to amplify only the mitochondrial cytb sequences.

Pairs of primers used in this study with relative references

The resulting PCR products were purified with a Sure Clean (Bioline) purification kit. Sequencing reactions were run under Big-Dye TM Terminator cycling conditions by the commercial company Macrogen (www.macrogen.com). The electropherograms were then checked using FinchTV and minor changes were made by eye. IUPAC ambiguity codes were used to represent polymorphism of heterozygote individuals in nuclear sequences. For the nuclear sequences, we solved the gametic phases using the coalescent-based Bayesian algorithm in PHASE v.2.1 [33, 34] as implemented in DnaSPv.5.1 [35]. Since the β-fibint7 showed the presence of INDELs polymorphisms, we first used the method described by Flot et al. [36] to determine the phase for sequences that were heterozygous for insertion or deletions (12 individuals). We then implemented the known phases in the coalescent-based Bayesian reconstruction. Three independent runs were conducted with burn-in at 1000 and 10,000 iterations, and thinning at each 100 steps. Only sequences with posterior probability >0.75 were included in the analysis. A final consensus alignment was computed for each marker with MEGA 5.0 [37]. When the final alignments were obtained, the number of haplotypes (H), nucleotides (π) and haplotype (h) diversity were estimated using DnaSPv.5.1 for the overall mtDNA and nuDNA datasets and for each mtDNA clade revealed by the phylogenetic analysis. Since the statistical power of tests for recombination is generally low, two different methods were used to assess for nuclear recombination. The four-gamete test was performed to estimate the minimum number of recombination events obtaining confidence intervals at 95% by the coalescent algorithm implemented in DnaSP v.5.1.

Moreover, we also test the occurrence of recombination events through the Pairwise Homoplasy Index (phi) test implemented in the program splitstree v. 4 [38].