Skip to main content
  • Research article
  • Open access
  • Published:

Evolutionary history of Podarcis tiliguerta on Corsica and Sardinia

Abstract

Background

Podarcis tiliguerta is a wall lizard endemic to the Mediterranean islands of Corsica and Sardinia. Previous findings of high mtDNA and morphological diversity have led to the suggestion that it may represent a species complex. Here, we analysed mitochondrial and nuclear markers (mtDNA, 3110 bp; 6 nDNA loci, 3961 bp) in P. tiliguerta sampled from thirty-two localities across Corsica and Sardinia.

Results

We find much greater intraspecific genetic divergence than between sister species of other Mediterranean island Podarcis, i.e., between P. lilfordi and P. pityusensis. We detected three mtDNA clusters in Corsica (North, South-East and South-West) and either two or three in Sardinia (North vs. South) depending on the clustering method. Only one or two nDNA groups were identified within each main island (again, depending on the method). A Bayesian time-calibrated multispecies coalescent tree was obtained from mtDNA and provided statistical support for a Miocene origin of the species (13.87 Ma, 95% HPD: 18.30–10.77 Ma). The posterior mean divergence time for the Corsican and Sardinian lineages was 12.75 Ma ago (95% HPD: 16.94–9.04 Ma).

Conclusion

The results support the evolutionary distinctiveness of Corsican and Sardinian populations and also indicate a lack of post-divergence migration despite periods of contact being possible. Further to this, species delimitation analyses of Corsican and Sardinian lineages provided statistical support for their recognition as distinct (sister) taxa. Our results provide new insights into the biogeography of the Mediterranean biodiversity hotspot, and contribute important findings relevant to the systematics and evolution of this speciose lizard genus.

Background

Islands provide many of the world’s biodiversity hotspots. They continuously generate new species with the term “speciation machines” having been coined to describe this phenomenon [1]. These high speciation rates are facilitated by specific island characteristics that lead to a wealth of ecological, biogeographic and evolutionary processes. These include: natural fragmentation, long term isolation, high altitudes creating habitat heterogeneity, complex intra-island landscapes, and island emergence due to volcanism and eustatic sea-level changes.

The Mediterranean basin was one of the first 25 Global Biodiversity Hotspots to be named [2] and is characterized by substantial island diversity [3]. The Tyrrhenian region within the Mediterranean is particularly diverse due to the influence of complex paleogeographical and historical factors [4, 5]. It includes Sardinia and Corsica, two of the largest and highest (1834 m and 2710 m, respectively), Mediterranean Islands, which are separated by a strait of only 11 km. They contain high levels of endemism including several hundred endemic plants. Geological history has undoubtedly played a role in generating this diversity. A land bridge maintained the connection between the islands and Europe [6] until their separation, completed 9 Ma ago [7, 8]. This was subsequently lost although connections between the Corsica–Sardinia archipelago and Europe and North Africa were re-established through desiccation of the Mediterranean basin during the Messinian Salinity Crisis (MSC; 5.96–5.33 Ma) [9, 10]. Corsica and Sardinia have since been temporarily connected as a result of more minor changes in sea level during glacial periods in the Pleistocene, with the last contact being during the Last Glacial Maximum (LGM) [11].

The lizard genus Podarcis encompasses about 19 species and is widespread across the Mediterranean. Podarcis tiliguerta is endemic to Corsica and Sardinia. It exhibits a great intraspecific variability in morphological traits, especially colour pattern and melanism, the latter found within populations from small islands away from the main islands [1217]. This variation has led to the description of several microinsular subspecies [18, 19]. Morphological [20, 21] and genetic studies, based on both allozymes [22] and mitochondrial sequences [23,24,25], have evidenced the existence of two or possibly three divergent genetic clusters within the Corsica-Sardinia archipelago. The underlying reasons for this divergence remain unexplained, but its magnitude has led some authors to suggest that P. tiliguerta might represent an insular species complex [20, 24]. Similar findings have been documented for other Podarcis from much wider geographic areas, e.g., P. hispanica, from the Iberian peninsula [26] and P. erhardii, from Greece [27]. Objective coalescent-based statistical analyses of species delimitation may provide an important step towards robust assignment of new taxonomic units based on genetic groups within the putative complex.

This study aims to provide the most detailed analysis of genetic variation of P. tiliguerta to date, using both mitochondrial (mtDNA) and nuclear DNA (nDNA) markers. We use these data to reveal new insights into the patterns and causes of the genetic diversity. We examine two important aspects of the evolutionary history of P. tiliguerta, namely the number of lineages and timing of divergence between these lineages and from other Podarcis. In addition, we investigate the historical biogeography of P. tiliguerta in relation to known physical changes in the region.

Methods

Sampling

Forty-one Podarcis tiliguerta were captured under license from 15 localities in Corsica, 13 localities in Sardinia (in May 2011) and four offshore islands and islets (September 2012) by noosing, with the aim of covering the main geographical regions (Fig. 1 and Table 1). Tail-tips were removed and stored in 100% ethanol before releasing individuals at sites of capture.

Fig. 1
figure 1

Sites at which Podarcis tiliguerta were sampled in Corsica, Sardinia and adjacent islands

Table 1 Localities and P. tiliguerta specimens sampled

DNA Amplification and sequencing

A standard phenol-chloroform protocol was used for DNA extraction [28]. Six non-overlapping mtDNA fragments were amplified and sequenced for each specimen containing sequences from the following genes: i) partial 12S rRNA, ii) partial tRNAGlu and partial cytochrome b (cytb), iii) partial cytb and partial tRNAThr, iv) partial control region (CR), v) two partial subunits of the NADH dehydrogenase gene and associated tRNAs (referred to as ND1, ND2, tRNAIle, tRNAGln, and tRNAMet) and vi) partial 16S rRNA. Primers and amplification conditions are the same as those used in our previous studies of Podarcis [29,30,31] (see Additional file 1), with the exception of the light strand primer (L14143; see Additional file 1) for amplification of fragment (ii) (above) in the Sardinian specimens. The 16S rRNA primers were obtained from Carranza et al. [32].

Six partial nuclear genes were sequenced: i) melanocortin-1 receptor gene (MC1R), ii) recombination-activating gene 1 (RAG1), iii) apolipoprotein B gene APOBE28, iv) lipoprotein gene BLC9L, v) transcription factor gene KIAA2018 and vi) kinesin family member KIF24. Primers [3335] are described in Additional file 1.

Sequencing of both strands of PCR products was carried out on an automated ABI 3130 sequencer (Applied Biosystems, Foster City, CA, USA) using a BigDye® Terminator v. 3.1 Cycle sequencing kit (Applied Biosystems, Foster City, CA, USA), and edited and aligned using BioEdit v.7.0.5.2. [36].

Analyses of genetic variability and structuring

DnaSP v.5.10.01 [37] was used to obtain basic genetic diversity indices, test for neutrality and phase nuclear data. Potential recombination events were assessed using the Pairwise Homoplasy Index (phi) test [38] implemented in SplitsTree v. 4 [39].

Relationships between nuclear DNA haplotypes were explored using statistical parsimony networks, with a 95% connection limit (program: TCS v.1.21 [40]).

Analysis of molecular variance (AMOVA) (program: Arlequin v.3.1.1 [41]) was used to examine genetic structuring between primary geographic groups (Corsica vs. Sardinia and Sardinian islands).

Bayesian phylogeographical and ecological clustering (BPEC) [42, 43] was used to identify genetically distinct geographical population clusters in Corsica and Sardinia. BPEC relies on parsimony in order to reduce the number of candidate trees. Different haplotype trees and migration events are explored through the Markov chain Monte Carlo (MCMC) sampler. Migration events are assumed to occur when a haplotype (with or without a mutation from its parent haplotype) migrates to a new geographical cluster. Locations were specified as latitudes and longitudes recorded during fieldwork using a motion X-GPS HD V21.0 Build 1737R (Fullpower Technologies Inc.). For each run, a strict parsimony criterion was applied and the maximum number of migrations was 4. The MCMC chain was run for 10 million iterations with 20,000 steps saved. Separate analyses were performed on Corsica and Sardinia. A subsample of 1371 mtDNA sites were selected for this analysis due to high levels of mtDNA divergence (saturation of third codon positions and other regions were detected using an appropriate test [44, 45]). The analysed sequence comprised: 12S rRNA (359 bp), CR (465 bp) and cytb codon positions 1 and 2 (548 bp). The geographical distributions of alleles were similarly analyzed for each nuclear locus, independently.

For comparison, an alternative genetic clustering procedure implemented within the program BAPS v.5.3 [46] was used to infer population structure without reference to geographical location. It was applied to all nuclear markers together, and to mtDNA alone. The number of clusters was determined and admixture analyses were executed to infer the ancestral source of the samples (without priors on geographic location).

Species phylogeny, divergence times and species delimitation

The multispecies coalescent approach implemented in *BEAST v. 1.7.4 [47, 48]; was used to estimate divergence times and phylogenetic relationship between specific mitochondrial lineages from the main islands, Corsica and Sardinia (with surrounding islets). The method allows estimation of divergence times while taking into account ancestral polymorphism. However, it often suffers from problems associated with prior specifications [49] and achieving good mixing of the MCMC chain (particularly when multiple loci are analysed). Also, the inclusion of several poorly informative nuclear markers could have a substantial and potentially undesirable impact on the estimated divergence times [50]. For these reasons, we included only mtDNA in the analysis.

Eight partitions were used to account for heterogeneity in sequence evolution and suitable models for each partition were identified using jModeltest v.0.1.1 [51] (AICc criterion, i.e., Akaike Information Criterion corrected for small samples). We partitioned the sequence by grouping functionally similar sites from specific regions as they should have similar properties (rates in particular). These were: i) 12S rRNA, ii) CR, iii) all tRNAs, iv) 1st and 2nd codon positions of the ND1/ND2 sequences, v) 3rd codon position of ND1/ND2, vi) 1st codon position of cytb, vii) 2nd codon position of cytb, viii) 3rd codon position of cytb, When the selected model was not available in *BEAST, we used the most similar available model.

We included mtDNA from all P. tiliguerta samples together with homologous sequences from two species that comprise a sister lineage to P. tiliguerta: P. lilfordi and P. pityusensis [52]. Inclusion of these taxa allowed incorporation of a time-calibrated node in the analysis [31, 53] (see Additional file 2 for GenBank accession numbers). The age of the (P. lilfordi, P. pityusensis) node on the population was constrained using the normal distribution N (5.325, 0.0001), where 1 unit = 1 Ma. The MSC is a well-known physical event in the Mediterranean and ended with very rapidly rising sea levels after the opening of the Mediterranean-Atlantic connection. Some studies have suggested that refilling of the basin may have only taken only a few thousand years [9, 10, 54]. This would have separated the Gymnesic and Pityusic Islands. Evidence from other studies (e.g., [55]) indicates the suitability of narrow time constraints on (Gymnesic, Pityusic) nodes. Refilling has been dated at 5.33 Ma [9] and given the deep intermediate channel between the island groups and the evidence of a very rapid refilling of the basin, the narrow normal distribution around a mean of 5.325 Ma would appear suitable.

A strict clock model was used in *BEAST (see [56]), and the tree prior was specified by a Yule process. The MCMC chain was run for 300 million generations sampled every 10,000 steps. Tracer v. 1.5 [57] allowed visual inspection of marginal posteriors. TreeAnnotator (BEAST package) was used to combine and analyse the posterior tree samples at stationarity (trees were combined using the maximum sum of clade credibilities criterion).

Species delimitation was examined using three different approaches. First the Bayesian Phylogenetics and Phylogeography program (BPP, v. 3.1) [58, 59] was used to obtain posterior probabilities on specified nodes. BPP implements a reversible jump Markov chain Monte Carlo (rjMCMC) to search the parameter space that includes θpopulation divergence and gene trees. Analyses were performed with i) mtDNA, ii) nDNA and iii) all markers, to determine the impact of mtDNA polymorphism on species delimitation. Gamma distributions (G) with the following shape and different scale parameters were used to specify the population size, θs: G(1,10) τ0 G (1,10), θs G(1,10) τ0 G(1, 2000), θs G(1, 2000) τ0 G (1,10), θs G(1, 2000) τ0 G (1, 2000). The MCMC chain was run for 400,000 steps (following a burnin of 20,000 steps), sampled every 20 steps. Each analysis was run three times to confirm consistency between runs.

Two additional approaches of species delimitation analysis suitable for single locus data (only mtDNA) were also tested. The Poisson Tree Processes model (PTP) [60] uses non-ultrametric phylogenetic trees and directly incorporates the number of nucleotide substitution in the model of speciation using a Poisson distribution to describe the probability that a substitution gives rise to a speciation event. This method delimits species under the assumption that the number of substitutions between species is considerably higher than within species. The branch lengths of the input tree are supposed to be generated by two independent classes of Poisson events, one corresponding to speciation and the other to coalescence. A bPTP analysis (a recent version of the PTP that implements Bayesian support values to the delimited species) was conducted on the bPTP web server (http://species.h-its.org/ptp/). A phylogenetic tree inferred by IQ-Tree [61] was used as input. The chain was run for 150,000 MCMC generations, 10% of samples were discarded as burn-in and sampling interval was 100.

Finally GMYC [62, 63], based on time-calibrated ultrametric phylogenetic trees and using differences in branching rates to infer species was applied. The input tree was generated in BEAST v1.8.2. [48] and the specifications were the same as the ones described above for the *BEAST analysis (8 partitions, models, strict clock, time-calibration and tree prior/yule process). Independent MCMC runs of 76 million steps were sampled every 10,000 steps (10% of samples were removed as burn-in). The GMYC analysis was conducted using this consensus tree and the single-threshold method within the SPLITS package [64] on R v 3.3.2 [65].

Results

Mitochondrial DNA

A total of 3110 bp of mtDNA sequence was obtained (12S rRNA 392 bp, control region 497 bp, ND1 85 bp, ND2 510 bp, cyt b 929 bp (two fragments), tRNAs 256 bp, 16S rRNA 431 bp). Sequences have been deposited in GenBank (for accession numbers see Additional file 3). The concatenated sequence contained 646 polymorphic sites, giving rise to 37 different haplotypes, which differed by a mean pairwise difference of 216.09 (±92.03). Indices of mtDNA diversity (Table 2; analyses without 16S rRNA sequences) were higher in P. tiliguerta than corresponding values for related species [29, 31], while non-significant neutrality tests [66,67,68,69] supported neutral evolution.

Table 2 Diversity estimates for the concatenated mtDNA fragment in P. tiliguerta (without 16S rRNA sequences) compared with others species of the Western Islands Group (P. lilfordi, P. pityusensis and P. filfolensis). Parameters for Corsica and Sardinia P. tiliguerta were also estimated independently

Nuclear data

The nuclear loci analysed were: MC1R (720 bp), RAG1 (939 bp), APOBE28 (489 bp), BLC9L (627–636 bp), KIAA2018 (644–659 bp) and KIF24 (497–518 bp). Sequences have been deposited in GenBank (Additional file 3). Basic genetic diversity indices and neutrality statistical tests [66, 69] are indicated in Table 3. The phi test [38] showed no statistically significant evidence for recombination in any gene: MC1R (p = 0.164), APOBE28 (p = 0.369), KIAA2018 (p = 0.095), BLC9L (p = 0.467), KIF24 (p = 0.322), RAG1 (p = 0.193).

Table 3 Diversity estimates for six nuclear exons in P. tiliguerta (length polymorphisms are included in these analyses)

Length polymorphisms were detected in three nuclear genes, including STRs (short tandem repeats) and INDELs (insertion-deletions). We briefly describe them here. The KIAA2018 sequence contained a polymorphic microsatellite with a proline repeat (CCT) ranging from 4-9. The (CCT)7 repeat was the most frequent (see Additional file 4). Only three individuals showed an amino acid change within this repeat. These were: i) two Sardinian samples TSA8 (CTT) and TSA18 (CYT) showed a leucine substitution in the fifth repeat, and ii) the Corsica sample TCO9 showed eight repeats with a threonine substitution in the final repeat. An STR was also observed within the BLC9L sequence with a variable number of repeats for glutamine and proline. The common structure was (Gln)3-4 (Pro)5-7 (Glu)5 with the unique exception of a Corsican sample (TCO2) that showed an histidine: (Gln)4 (Pro)5 (His) (Gln)5. Three different INDELs were detected in KIF24 sequences: i) a 9 bp deletion (CATTTTGGT) in specimen TSA5 and ii) a 12 bp insertion (AAGGACTTTGGG) in TSA1–3, TSA5, TSA13 and Tf1–2, iii) 6 bp deletions (GAAAGC) in TCO6 and TCO10; and (AGAAAG) in TCO9. Note that i) and ii) occurred only in Sardinia, while iii) occurred only in Corsica.

Genetic structure of the insular populations

The mtDNA variance was partitioned as 51.01% (i.e., FST = 0.51) between Corsica and Sardinia (and offshore islands) and the remaining variance within these regions (AMOVA analysis). BAPS analysis of the mtDNA locus defined five genetic clusters (lnL = -8303.607) (Fig. 2a). The analysis distinguishes three Corsican groups: i) cluster A includes lizards from the south-western region (TCO2–4), ii) cluster B includes lizards from the south-east (TCO1, 12–14, 17, 20–21) and iii) cluster C includes lizards from the north of Corsica (TCO5–11, 18–19). Two Sardinian clusters are identified: cluster D includes specimens corresponding to localities in the south of the island (TSA12–13, 15, 18–20) and cluster E contains P. tiliguerta from Sardinian islands and islets (Molara Island, Foradada Island, Paduleddi Islet and Stramanari Islet) and from the north of the main island (TSA1–5, 10–11, 17).

Fig. 2
figure 2

Overall genetic structure inferred from a) mtDNA and b) concatenated nuclear loci across the two main geographic regions used in this study. Mixture analyses were estimated by BAPS software. In the bar plot, vertical bars represent individuals and proportions of admixture, with colours corresponding to ancestral sources

AMOVA of the nuclear loci indicated that 18.2% (i.e., FST = 0.182) of the variance occurred between Corsica and Sardinia (and offshore islands) with most of the variance being within these regions. Only two genetic clusters were defined by a BAPS analysis of concatenated nDNA (lnL = -5640.911). The first included all Corsican individuals, and the second included Sardinian and individuals from small islet populations (Fig. 2b). The TCS analysis indicated that two loci (APOBE28 and MC1R) show star-like parsimony networks, with a main haplotype shared by Corsican and Sardinian samples, and derived haplotypes belonging to both islands (Fig. 3). In contrast, RAG1 and BLC9L show compound star-like networks consisting of two central haplotypes, corresponding to Corsica and Sardinia, separated by one mutational step. Finally, the two other networks (KIF24 and KIAA2018) exhibit a clear separation between Corsican and Sardinian samples, together with a much higher differentiation between haplotypes (30 mutational steps). The KIAA2018 sequences show greatest diversity of all nuclear loci analysed, characterised by 66 different haplotypes.

Fig. 3
figure 3

Nuclear haplotype networks for Podarcis tiliguerta

Mitochondrial haplotypes from Sardinia were assigned by BPEC to three phylogeographic clusters with high posterior probabilities (>0.7) (Fig. 4a), with most likely ancestral locations being Gallura, Villasimius and Costa rei (see Fig. 1 and Table 1). Three clusters were also obtained for Corsica (Fig. 4a) with likely ancestral locations being Sant Bastiano, Corte and Vivario (see Fig. 1 and Table 1). Low posterior cluster assignment probabilities were obtained for eight Corsican haplotypes (~0.4–0.5), and one Sardinian (~0.5) haplotype. Geographical structuring in Corsica corresponded to the north, south-east and south-west of the island, while in Sardinia three geographical groupings were detected: i) north (with some offshore islands), ii) south and iii) Foradada Island (Fig. 4a). BPEC analyses of nuclear loci assigned haplotypes to four phylogeographic clusters with generally high posterior probabilities: two of them, in Corsica and two, in Sardinia (Fig. 4b). However all nuclear genes divide the main islands into northern and southern regions.

Fig. 4
figure 4

Analyses of Bayesian phylogeographical and ecological clustering (BPEC): a) using mtDNA or b) RAG1 as example of nuclear DNA. Each coloured contour plot indicates a different phylogeographical clusters. The contour plots are centred at the ‘centre’ of each population cluster, and the coloured regions show the radius of 50% concentration contours around it. White arrows show the ancestral locations with the highest posterior probability for each island

Species phylogeny, divergence times and species delimitation

Preliminary gene tree analyses were performed on all individual loci prior to the Bayesian multispecies coalescent analysis. These analyses ruled out the possibility that P. tiliguerta was polyphyletic relative to the other closely related taxa that were included (P. hipanica, P. guadarramae, P. virescens, P. liolepis, P. lilfordi, P. pityusensis, and P. filfolensis) (data not shown).

The species tree (Fig. 5), based on mtDNA, provided a posterior mean for the divergence of P. tiliguerta from (P. lilfordi, P. pityusensis) at 13.87 Ma (95% highest posterior density, HPD: 18.30–10.77 Ma), and that for the divergence of Corsica and Sardinia lineages, at 12.75 Ma (95% HPD: 16.94–9.04).

Fig. 5
figure 5

Mitochondrial populations tree chronogram estimated by *BEAST in P. tiliguerta. P. lilfordi and P. pityusensis are used as outgroups. Grey bars correspond to the 95% highest posterior density intervals for each divergence time

All BPP analyses supported the two species (Corsica, Sardinia) model (P = 1.00) irrespective of the prior specifications for θs and τ0 and whether mtDNA was included or excluded from the analyses.

The trees obtained from the bPTP analysis (PTP_ML and PTP-sh) (excluding outgroups) showed 11 groups with high posterior probabilities (≥0.91). Of the 11 clusters, 8 comprised only 1 sample. The GMYC approach detected 9 clusters within P. tiliguerta with a total of 19 entities, excluding outgroups and with posterior probabilities ≥ 0.85. The different clusters detected by these two analyses corresponded to different clades detected within Corsica and Sardinian islands, but none of these groups comprised both Corsica and Sardinia individuals, further supporting species delimitation of these two island forms. Additional information about PTP and GMYC analyses is available on request.

Discussion

The most striking finding is the high degree of genetic variability detected within the taxon Podarcis tiliguerta. This has already been suggested by two previous studies that analysed protein mobility and mtDNA [22, 24]. Other insular Mediterranean species, such as P. pityusensis, P. lilfordi or P. filfolensis show lower levels of divergence, although they inhabit smaller islands. For exemple, the mean pairwise differences for comparable mtDNA regions is more than four times higher than in P. lilfordi, and is almost 22 times greater than in P. filfolensis or P. pityusensis. In addition, intraspecific nDNA variability also appears to be greater P. tiliguerta than in other Podarcis which confirms the extreme genetic diversity observed within this species [24, 70].

The high diversity is primarily accounted for by strong geographical structuring with a deep, ancient divergence between Corsican and Sardinia lineages as well as substantial within-island divergence. Early analyses of allozymes and short mtDNA fragments indicated two genetic groupings: Corsica and Sardinia [22, 24]. Subsequent analyses, with wider sampling, revealed three mtDNA groups: two in Sardinia and one in Corsica [25, 70]. Here, we detect three geographically-structured mtDNA groups within Corsica, corresponding to the north, south-east and south-west of the island. We also detect two/three (depending on the analysis) spatial groupings of mtDNA within Sardinia: one of which corresponds to the north of the main island and associated smaller islands, while the other corresponds to the south.

As expected, the nuclear data detected fewer clusters. The BAPS analyses revealed two groups: Corsica and Sardinia, while the slightly more sensitive BPEC analyses revealed two additional clusters within each island. Bruschi et al. [20] studied P. tiliguerta using morphological characters and suggested hybridization was occurring in southern Corsican populations due to introductions of individuals from Sardinia which might contribute to these patterns. However, our nuclear and mtDNA analyses found no evidence of this. Although some nuclear alleles were shared between Corsican and Sardinian populations, the lack of any such patterns in the mtDNA suggests that incomplete lineage sorting is the most likely cause.

The ancestral species node for P. tiliguerta was dated at 12.75 Ma (Fig. 5; 95% HPD: 16.94–9.04 Ma) confirming preliminary divergence time estimates (13 Ma) [71]. P. pityusensis and P. lilfordi are related species to P. tiliguerta, including all them in the named “Western Islands group” [72]. This clearly predates the intraspecific divergence within: i) P. lilfordi which dates to the early Pleistocene [53], ii) P. filfolensis which also diverged in the Pleistocene [30] and iii) P. pityusensis which diverged in the mid-late Pleistocene [31]. The divergence time between Corsica and Sardinia far exceeds the usual limits for distinct species within the Squamata [73, 74]. Reciprocal monophylia between the islands was found for several loci, which supports our detection of species delimitation (although not for all loci, as would be expected between distinct species). This, and the general concordance with and among nuclear markers, explains why our species delimitation analysis provides such strong support for separation of P. tiliguerta into two species. In a future study we intend to provide formal (morphological and ecological) descriptions of these new taxa.

P. tiliguerta appears to have originated by allopatric speciation, from continental ancestors that colonized the Corsica-Sardinia microplate. According to Arnold et al. [52], the Lacertini spread and diversified around 16–12 Ma BP producing the ancestors of the present genera, so genus origin postdates the isolation of the Corsica-Sardinia microplate from the Iberian plate (29 Ma). This indicates that the ancestors reached it subsequently via the land bridge that it maintained with the continent during its rotation. Later, the microplate split into two islands (15–9 Ma ago [7, 8]) leading to isolation of the Sardinian and Corsican lineages, concordant with our divergence time estimate. Since separation, the heterogeneous orography of Corsica and Sardinia could have played an important role in determining intraspecific diversity.

Conclusions

In summary, we analysed mtDNA and nuclear markers in P. tiliguerta and found high genetic variability and a deep genetic structure, corresponding to deep Miocene divergence between Corsica and Sardinia. In addition, we detected substantial within-island divergence but were unable to relate it to any geological or other events in the islands past. This study will motivate and contribute evidence to a taxonomic reassessment of P. tiliguerta with respect to the Corsican and Sardinian forms.

Abbreviations

AICc:

Akaike information criterion corrected

AMOVA:

Analysis of Molecular Variance

APOBE28 :

Apolipoprotein B gene exon 28

BAPS:

Bayesian Analysis of Population Structure

BEAST:

Bayesian Evolutionary Analysis by Sampling Trees

BLC9L :

Lipoprotein gene

BPEC:

Bayesian Phylogeographical and Ecological Clustering

BPP:

Bayesian Phylogenetics and Phylogeography

CR:

Control Region

cytb :

Cytochrome b

G:

Gamma distribution

HPD:

Highest Posterior Density

INDEL:

Insertion-deletion

KIAA2018 :

Transcription factor gene

KIF24 :

Kinesin Family Member 24

LGM:

Last Glacial Maximum

MC1R :

Melanocortin-1 Receptor

MCMC:

Markov Chain Monte Carlo

MRCA:

Most Recent Common Ancestor

MSC:

Messinian Salinity Crisis

mtDNA:

mitochondrial DNA

ND :

Nicotine adenine dinucleotide gene

nDNA:

Nuclear DNA

phi:

Pairwise Homoplasy Index

RAG1 :

Recombination-Activating Gene 1

rjMCMC:

reversible jump Markov chain Monte Carlo

rRNA:

Ribosomal RNA

STR:

Short Tandem Repeat

tRNA:

transfer RNA

References

  1. Rosenzweig ML. Species diversity in space and time. Cambridge: University Press; 1995

  2. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403(6772):853–8.

    Article  CAS  PubMed  Google Scholar 

  3. Fernández-Palacios JM. Why islands. Islands Evolution. 2010;85–109.

  4. Verlaque R, Médail F, Quézel P, Babinot J-F. Endémisme végétal et paléogéographie dans le bassin méditerranéen. Geobios. 1997;30:159–66.

    Article  Google Scholar 

  5. Médail F, Quézel P. Biodiversity Hotspots in the Mediterranean Basin: Setting Global Conservation Priorities. Conserv Biol. 1999;13(6):1510–3.

    Article  Google Scholar 

  6. Meulenkamp JE, Sissingh W. Tertiary palaeogeography and tectonostratigraphic evolution of the Northern and Southern Peri-Tethys platforms and the intermediate domains of the African–Eurasian convergent plate boundary zone. Palaeogeog Palaeoclim Palaeoecol. 2003;196(1):209–28.

    Article  Google Scholar 

  7. Alvarez W. Rotation of the Corsica–Sardinia microplate. Nature. 1972;235(58):103–5.

    Google Scholar 

  8. Alvarez W. Fragmentation of the Alpine orogenic belt by microplate dispersal. Nature. 1974;248:309–14.

    Article  Google Scholar 

  9. Krijgsman W, Hilgen F, Raffi I, Sierro F, Wilson D. Chronology, causes and progression of the Messinian salinity crisis. Nature. 1999;400(6745):652–5.

    Article  CAS  Google Scholar 

  10. Duggen S, Hoernle K, Van Den Bogaard P, Rüpke L, Morgan JP. Deep roots of the Messinian salinity crisis. Nature. 2003;422(6932):602–6.

    Article  CAS  PubMed  Google Scholar 

  11. Grill A, Casula P, Lecis R, Menken S. Endemism in Sardinia. In: Phylogeography of southern European refugia. Netherlands: Springer; 2007. p. 273–96.

  12. Lanza B. The natural history of the Cerbicale islands (south-eastern Corsica) with particular reference to their herpetofauna. (Storia naturale delle Isole Cerbicale (Corsica sudorientale) con particolore riferimento alla loro erpetifauna). Natura. 1972;63(4):345–407.

    Google Scholar 

  13. Lanza B, Brizzi R. On two new Corsican microinsular subspecies of Podarcis tiliguerta (Gmelin, 1789)(Reptilia: Lacertidae). Natura. 1974;65(3-4):155–93.

    Google Scholar 

  14. Lanza B, Poggesi M. Storia naturale delle isole satelliti della Corsica. Firenze: Istituto geografico militare; 1986.

  15. Arnold EN, Ovenden D. Reptiles and amphibians of Europe, vol. 288: Princeton: University Press Princeton; 2002.

  16. Brizzi R, Lanza B. The natural history of the Macinaggio Islets (northeastern Corsica) with particular reference to their herpetofauna. Natura. 1975;66(1-2):53–72.

    Google Scholar 

  17. Lanza B. On a new wall lizard from Corsica, with notes on the Sanguinarie Islands (Reptilia: Lacertidae). Natura. 1976;67(3-4):185–202.

    Google Scholar 

  18. Lanza B, Corti C. Evolution of the knowledge on the Italian herpetofauna during the 20th century. Bollettino del Museuo Civico di Storia Naturale di Verona. 1996;20(1993):373–436.

    Google Scholar 

  19. Corti C, Cascio PL. The lizards of Italy and adjacent Areas. In: Frankfurter Beiträge zur NaturKunde, Vol. 14: Edition. Frankfurt: Chimaira; 2002.

  20. Bruschi S, Corti C, Carretero M, Harris D, Lanza B, Leviton A. Comments on the status of the Sardinian-Corsican Lacertid lizard Podarcis tiliguerta. Proc Calif Acad Sci. 2006;57(1/11):225.

    Google Scholar 

  21. Bruschi S, Corti C, M. C: Podarcis tiliguerta (Gmelin, 1789). In: Fauna d’Italia Reptilia. Edited by Corti C, Capula M, Luiselli L, Razzetti E, Sindaco R. Bologna: Edizioni Calderini de Il Sole 24 Ore Editoria Specializzata S.r.l; 2010.

  22. Capula M. Evolutionary genetics of the insular lacertid lizard Podarcis tiliguerta: genetic structure and population heterogeneity in a geographically fragmented species. Heredity (Edinb). 1996;77(5):518–29.

    Article  Google Scholar 

  23. Pinho C, Corti C, Carretero M, Harris DJ. Genetic variability within Podarcis tiliguerta: preliminary evidence from 12S rRNA gene sequences. In: Fifth International Symposium on the Lacertids of the Mediterranean Basin, Lipari, Aeolian Islands, Sicily, Italy. 2004. p. 7–11.

  24. Harris D, Pinho C, Carretero M, Corti C, Böhme W. Determination of genetic diversity within the insular lizard Podarcis tiliguerta using mtDNA sequence data, with a reassessment of the phylogeny of Podarcis. Amphibia-Reptilia. 2005;26(3):401–7.

    Article  Google Scholar 

  25. Podnar M, Mayer W. Can mitochondrial DNA draw the phylogenetic picture of Central Mediterranean island Podarcis? Herpetoza. 2005;18(1/2):73–7.

    Google Scholar 

  26. Pinho C, Ferrand N, Harris DJ. Reexamination of the Iberian and North African Podarcis (Squamata: Lacertidae) phylogeny based on increased mitochondrial DNA sequencing. Mol Phylogenet Evol. 2006;38(1):266–73.

    Article  CAS  PubMed  Google Scholar 

  27. Poulakakis N, Lymberakis P, Antoniou A, Chalkia D, Zouros E, Mylonas M, Valakos E. Molecular phylogeny and biogeography of the wall-lizard Podarcis erhardii (Squamata: Lacertidae). Mol Phylogenet Evol. 2003;28(1):38–46.

    Article  CAS  PubMed  Google Scholar 

  28. Sambrook J, Fritsch E, Maniatis T. Molecular Cloning: A Laboratory Manual. 2nd ed. Cold Spring Harbor: Cold Spring Harbor Laboratory Press; 1989.

  29. Terrasa B, Pérez-Mellado V, Brown RP, Picornell A, Castro JA, Ramon MM. Foundations for conservation of intraspecific genetic diversity revealed by analysis of phylogeographical structure in the endangered endemic lizard Podarcis lilfordi. Diversity Distrib. 2009;15(2):207–21.

    Article  Google Scholar 

  30. Rodríguez V, Brown R, Terrasa B, Pérez-Mellado V, Picornell A, Castro J, Ramon C. Genetic diversity and historical biogeography of the Maltese wall lizard, Podarcis filfolensis (Squamata: Lacertidae). Conserv Genet. 2014;15(2):295–304.

    Article  Google Scholar 

  31. Rodriguez V, Brown RP, Terrasa B, Perez-Mellado V, Castro JA, Picornell A, Ramon MM. Multilocus genetic diversity and historical biogeography of the endemic wall lizard from Ibiza and Formentera, Podarcis pityusensis (Squamata: Lacertidae). Mol Ecol. 2013;22(19):4829–41.

    Article  CAS  PubMed  Google Scholar 

  32. Carranza S, Arnold EN, Thomas RH, Mateo J, López-Jurado L. Status of the extinct giant lacertid lizard Gallotia simonyi simonyi (Reptilia: Lacertidae) assessed using mtDNA sequences from museum specimens. Herpetol J. 1999;9(2):83–6.

    Google Scholar 

  33. Buades JM, Rodriguez V, Terrasa B, Perez-Mellado V, Brown RP, Castro JA, Picornell A, Ramon MM. Variability of the mc1r Gene in Melanic and Non-Melanic Podarcis lilfordi and Podarcis pityusensis from the Balearic Archipelago. PLoS One. 2013;8(1):e53088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Mayer W, Pavlicev M. The phylogeny of the family Lacertidae (Reptilia) based on nuclear DNA sequences: Convergent adaptations to arid habitats within the subfamily Eremiainae. Mol Phylogenet Evol. 2007;44(3):1155–63.

    Article  CAS  PubMed  Google Scholar 

  35. Portik D, Wood P, Grismer J, Stanley E, Jackman T. Identification of 104 rapidly-evolving nuclear protein-coding markers for amplification across scaled reptiles using genomic resources. Conserv Genet Resour. 2012;4(1):1–10.

    Article  Google Scholar 

  36. Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.

    CAS  Google Scholar 

  37. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.

    Article  CAS  PubMed  Google Scholar 

  38. Bruen TC, Philippe H, Bryant D. A simple and robust statistical test for detecting the presence of recombination. Genetics. 2006;172(4):2665–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23(2):254–67.

    Article  CAS  PubMed  Google Scholar 

  40. Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–9.

    Article  CAS  PubMed  Google Scholar 

  41. Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol Bioinform Online. 2005;1:47–50.

    CAS  Google Scholar 

  42. Manolopoulou I, Emerson BC. Phylogeographic ancestral inference using the coalescent model on haplotype trees. J Comput Biol. 2012;19(6):745–55.

    Article  CAS  PubMed  Google Scholar 

  43. Manolopoulou I, Legarreta L, Emerson BC, Brooks S, Tavaré S. A Bayesian approach to phylogeographic clustering. Interface focus. 2011;1(6):909–21.

  44. Xia X, Lemey P. Assessing substitution saturation with DAMBE. In: The Phylogenetic Handbook: A Practical Approach to DNA and Protein Phylogeny. Edited by: Lemey P, Salemi M, Vandamme A-M. Cambridge: University Press; 2009. p. 615–30.

  45. Xia X, Xie Z, Salemi M, Chen L, Wang Y. An index of substitution saturation and its application. Mol Phylogenet Evol. 2003;26(1):1–7.

    Article  CAS  PubMed  Google Scholar 

  46. Corander J, Waldmann P, Sillanpää MJ. Bayesian analysis of genetic differentiation between populations. Genetics. 2003;163(1):367–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Heled J, Drummond AJ. Bayesian Inference of Species Trees from Multilocus Data. Mol Biol Evol. 2010;27(3):570–80.

    Article  CAS  PubMed  Google Scholar 

  48. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29(8):1969–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Heled J, Drummond AJ. Calibrated Tree Priors for Relaxed Phylogenetics and Divergence Time Estimation. Syst Biol. 2012;61(1):138–49.

    Article  PubMed  Google Scholar 

  50. Brown RP, Yang Z. Bayesian Dating of Shallow Phylogenies with a Relaxed Clock. Syst Biol. 2010;59(2):119–31.

    Article  CAS  PubMed  Google Scholar 

  51. Posada D. jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008;25(7):1253–6.

    Article  CAS  PubMed  Google Scholar 

  52. Arnold EN, Arribas O, Carranza S. Systematics if the Palaearctic and Oriental lizard tribe Lacertini (squamata: Lacertidae: Lacertinae), with descriptions of eight new genera. Zootaxa. 2007;1430:1–86.

    Article  Google Scholar 

  53. Brown RP, Terrasa B, Perez-Mellado V, Castro JA, Hoskisson PA, Picornell A, Ramon MM. Bayesian estimation of post-Messinian divergence times in Balearic Island lizards. Mol Phylogenet Evol. 2008;48(1):350–8.

    Article  CAS  PubMed  Google Scholar 

  54. Pierre C, Caruso A, Blanc-Valleron M-M, Rouchy JM, Orzsag-Sperber F. Reconstruction of the paleoenvironmental changes around the Miocene–Pliocene boundary along a West–East transect across the Mediterranean. Sediment Geol. 2006;188:319–40.

    Article  Google Scholar 

  55. Martínez-Solano I, Gonçalves H, Arntzen J, García-París M. Phylogenetic relationships and biogeography of midwife toads (Discoglossidae: Alytes). J Biogeogr. 2004;31(4):603–18.

    Article  Google Scholar 

  56. Brown RP, Yang Z. Rate variation and estimation of divergence times using strict and relaxed clocks. BMC Evol Biol. 2011;11(1):1.

    Article  Google Scholar 

  57. Rambaut A, Drummond A. Tracer 1.5. Edinburgh: University of Edinburgh. Available at: http://tree.bio.ed.ac.uk/software/tracer 2007.

  58. Yang Z, Rannala B. Bayesian species delimitation using multilocus sequence data. Proc Natl Acad Sci U S A. 2010;107(20):9264–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Rannala B, Yang Z. Improved reversible jump algorithms for Bayesian species delimitation. Genetics. 2013;194(1):245–53.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Zhang J, Kapli P, Pavlidis P, Stamatakis A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics. 2013;29(22):2869–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. Iq-tree: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74.

    Article  CAS  PubMed  Google Scholar 

  62. Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Syst Biol. 2006;55(4):595–609.

    Article  PubMed  Google Scholar 

  63. Fujisawa T, Barraclough TG. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent (GMYC) approach: a revised method and evaluation on simulated datasets. Syst Biol. 2013;62(5):707–24.

  64. Ezard T, Fujisawa T, Barraclough T. SPLITS: species’ limits by threshold statistics. R package version. 2009. http://R-Forge.R-project.org/projects/splits/.

  65. Team RC. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2015. URL http://www.R-project.org2016.

  66. Fu Y-X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  67. Fu Y-X, Li W-H. Statistical tests of neutrality of mutations. Genetics. 1993;133(3):693–709.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Fay JC, Wu C-I. Hitchhiking under positive Darwinian selection. Genetics. 2000;155(3):1405–13.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. Vasconcelos R, Harris DJ, Carretero MA, Pinho C, Corti C, Capula M, Bassu L, Spano G, Delauguerre M. Genetic diversity within Corsican and Sardinian specimens of the Tyrrhenian Wall Lizard, Podarcis tiliguerta, estimated using mtDNA sequences. In: Mainland and insular lizards a Mediterranean perspective. Florence: Firenze University Press; 2006. p. 199–207.

  71. Oliverio M, Bologna MA, Monciotti A, Annesi F, Mariottini P. Molecular phylogenetics of the Italian Podarcis lizards (Reptilia, Lacertidae). Ital J Zool. 1998;65(3):315–24.

    Article  CAS  Google Scholar 

  72. Harris DJ, Arnold EN. Relationships of Wall Lizards, Podarcis (Reptilia: Lacertidae) Based on Mitochondrial DNA Sequences. Copeia. 1999;1999(3):749–54.

    Article  Google Scholar 

  73. Hasbun CR, Gomez A, Köhler G, Lunt DH. Mitochondrial DNA phylogeography of the Mesoamerican spiny-tailed lizards (Ctenosaura quinquecarinata complex): historical biogeography, species status and conservation. Mol Ecol. 2005;14(10):3095–107.

    Article  CAS  PubMed  Google Scholar 

  74. Pinho C, Harris D, Ferrand N. Contrasting patterns of population subdivision and historical demography in three western Mediterranean lizard species inferred from mitochondrial DNA variation. Mol Ecol. 2007;16(6):1191–205.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Our acknowledgement to CERFA for the license for field work in May 2011 (permit No. 13616*01) and the Ministero dell’Ambiente e della Tutela del Territorio e del Mare, that authorized sampling of lizards in Sardinia (authorization DPN-2010-0017564). We thank M. Bassitta and C. Colom for 16S rRNA amplification and sequencing.

Funding

This work was financed by the research grant CGL2015-68139-C2-1-P of the Ministerio Español de Economía y Competitividad and European Regional Development Fund (ERDF). VR (FPI/2009) and JMB (FPI/1488/2012) were funded by fellowships (FPI/1488/2012) from the Conselleria d'Educació, Cultura i Universitats (Govern de les Illes Balears), co-financed by the European Social Fund.

Availability of data and materials

DNA sequences: GenBank Accessions numbers of the mtDNA (KM058760–KM058886) and from the nDNA (KX237247–KX237504, JX126691–93).

Authors’ contributions

CR, JAC and AP conceived and designed the experiments; VPM, MD and CC sampled; Performed the experiments and analysed the data: VR, JMB and RPB performed the experiments and analysed the data; VPM and BT contributed with reagents/materials/analysis tools; Wrote the paper: VR, JMB, RPB and CR. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval

Animal Care and Use Committee approval was not applicable for this study.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to M. M. Ramon.

Additional files

Additional file 1:

Primers sequences used in amplification and sequencing. (DOCX 17 kb)

Additional file 2:

Podarcis spp analysed as outgroup in *BEAST analyses. (DOCX 11 kb)

Additional file 3:

Sample accession numbers for mitochondrial and nuclear gene. (DOCX 24 kb)

Additional file 4:

STRs observed in KIAA2018 gene (664 bp) included a proline repeat (CCT motif). The number of repeats ranged 4 to 9 (12 bp to 27 bp). (DOCX 12 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Rodríguez, V., Buades, J.M., Brown, R.P. et al. Evolutionary history of Podarcis tiliguerta on Corsica and Sardinia. BMC Evol Biol 17, 27 (2017). https://doi.org/10.1186/s12862-016-0860-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-016-0860-4

Keywords