Mitochondrial phylogeography of the Iberian endemic frog Rana iberica, with implications for its conservation.

Genetic characterization of species using phylogeographic approaches represents a basic reference to understand their evolutionary history as well as to identify conservation priorities to protect areas of particular interest regarding evolutionary potential. Even in well-studied regions such information is lacking for the majority of species, including many endemic species with reduced distribution ranges. We investigate the phylogeographic pattern of the Iberian frog Rana iberica, an endemic amphibian restricted to Central and North-Western Iberian Peninsula. Using mitochondrial sequences, we reconstruct the phylogeographic history of the species to test the effect of Quaternary climate changes on the evolutionary diversification of lineages, that is, the differentiation of mitochondrial lineages and the formation of genetic diversity melting pots, and integrate phylogeographic evidence for future conservation planning. Our results indicate the existence of 3 main mitochondrial lineages differentiated during the Upper Pleistocene. Both historical demographic analyses and climatic niche modeling show a strong effect of glacial climate changes, suggesting recurrent range contractions and expansions. Under such circumstances, differentiation took place most likely by isolation in allopatric interglacial refugia. Secondary lineage admixture in northern Portugal generated a broad mixed zone with highest nucleotide diversity. Given its particular evolutionary potential, its reduced distribution and eventual threats under current climate change scenario, conservation priorities should focus on the isolated lineage from Sierra de Guadalupe.

Conservation genetics takes advantage of the non-uniform geographic distribution of genetic diversity to select and design the best-fitting areas to conserve a species while preserving its genetic signal. Two of the most used criteria to define these conservation areas are the existence of genetically distinct, but geographically isolated populations (see definitions by Moritz 2002) and the recognition of hotspots of intraspecific genetic diversity, that is geographic areas where genetic diversity is higher than in the remaining areas (Canestrelli et al. 2014). While the first criterion has been extensively dealt with (Lesica and Allendorf 1995; Moritz 2002; Vandergast et al. 2008; Thomassen et al. 2011), the second is often treated lightly, forgetting about the evolutionary potential of these genetically complex areas. The nature and origin of intraspecific genetically diverse areas are subject to intense study in the fields of ecology and evolutionary biology (centers of origin; climatic sanctuaries; secondary contact and admixture, etc.), but often overlooked for conservation purposes, since delimitation of protected areas is often focused on current patterns of organismal and ecosystem-level biodiversity, typically ignoring the evolutionary processes involved in the origin of that biodiversity (Vandergast et al. 2008). Neglecting the nature of genetically diverse areas makes difficult to decide on the future conservation value of the chosen area but might also alter the evolutionary potential of a given species.

Some generalities have been observed when comparing phylogeographic structure among species within large glacial refugia, such as the Mediterranean peninsulas (Gómez and Lund 2007), although patterns of genetic diversity are often varied and reflect different, specific responses to the same events (Nieto Feliner 2011). One is that the majority of species exhibit a previously unanticipated level of fragmentation, indicating the occurrence of multiple refugia within the peninsulas. This is in accordance with the hypothesis that refugia are characterized by stable ecological conditions that have allowed the survival of large populations and, in consequence, the persistence of a relatively large part of their genetic legacy (Abellán and Svenning 2014). In several cases it has been possible even to identify pre-Pleistocene lineages, usually restricted to particular geographical areas that can be defined as sanctuaries of intraspecific diversity (Recuero and García-París 2011). In others, differentiation took place by allopatric fragmentation during the repeated Pleistocene climate changes (Knowles 2001; Canestrelli et al. 2014). Every genetic lineage may represent an evolutionary unit with potentially unique genetic characteristics and adaptive traits within a species, and therefore should be taken into consideration in any specific conservation plan (Hampe and Petit 2005). Another assumed generalization is that population expansions from multiple refugia generate areas of low genetic diversity but, in some cases, may produce admixture zones between divergent lineages. In the later case, high levels of genetic diversity can be reached in what are called melting pots of diversity (Petit et al. 2003; Dufresnes et al. 2016), in which genetic diversity is enriched by the immigration of individuals from different genetic lineages. These highly genetically diverse populations should also be considered for their conservation relevance, considering that they should present wider alternatives for coping with both natural and anthropogenic environmental changes (Sgrò et al. 2011). Endemic species, with reduced distribution areas and often specific ecological constrains, are particularly sensitive to current levels of anthropogenic habitat alterations and in many cases subject to different levels of protection. Phylogeographic studies for such taxa should be one of the first steps to efficiently plan conservation strategies. However, for most species genetic diversity characterization is incomplete if not totally missing.

This is the case of Rana iberica (Boulenger 1879), an endemic species restricted to the northwestern Iberian Peninsula and categorized as Near Threatened by the IUCN for its currently decreasing population trend (Tejedo et al. 2009). This species preferentially inhabits swift streams in humid, mountainous areas, where it often co-occurs with other endemics such as Chioglossa lusitanica (Bocage 1864), Lissotriton boscai (Lataste, 1879), and Lacerta schreiberi (Bedriaga 1878; Teixeira 2008). In addition, R. iberica also has population isolates in its southern distribution range, possibly raising relevant conservation issues. In spite of the potential phylogeographical interest of this species, only a few studies concerning the population genetics of R. iberica are available (Arano et al. 1993; Martínez-Solano et al. 2005), none covering the whole distribution range of the species and characterizing phylogeographical patterns.

In this study, we investigate the mitochondrial DNA (mtDNA) phylogeographic patterns in a large sample of R. iberica comprising its complete geographical distribution. Specifically, we test the effect of Quaternary climate changes on the evolutionary diversification of lineages, that is, the differentiation of mitochondrial lineages and the formation of genetic diversity melting pots by secondary admixture, and integrate the phylogeographic evidences here presented for future conservation planning and immediate actions for this species.

Materials and Methods

Sampling

We analyzed 162 specimens of R. iberica from 57 localities across northwestern Iberia (Figure 1 and Table 1), covering all the species range across Portugal, Northern Spain, and all known population isolates of the species. Tissue samples consisted of tail ends of larvae or toe clips of adult specimens. All individuals were released in the wild after tissue collection. Sample sizes per locality are indicated in Table 1.

Figure 1.

Distribution map of Rana iberica in the Iberian Peninsula (shaded), with indication of the geographical origin of samples analyzed in this study and proportion of lineages by locality: lineage A in gray, lineage B in white (Sistema Central sublineage gridded), lineage C in black. Population codes as in Table 1.

Table 1.

Sampling localities, population codes, sample size (n), haplotypes, and GeneBank accession numbers for the Cox1 sequences obtained in this study

Locality	Code	n	haplotypes	GeneBank Accession Nos
Portugal: Viana do Castelo	VIA	3	H71	MG813567–MG813569
Portugal: Ponte de Lima	LIM	4	H19, H43	MG813570–MG813573
Portugal: Paredes de Coura	COU	3	H11, H44, H71	MG813574–MG813576
Portugal: Castro Laboreiro	LAB	2	H18, H70	MG813577–MG813578
Portugal: Gerês	GER	6	H11, H15, H46, H54, H71	MG813579–MG813584
Portugal: Boticas	BOT	2	H23, H41	MG813586–MG813587
Portugal: Montesinho	MON	6	H16, H23, H25, H42, H46	MG813588–MG813593
Portugal: Vinhais	VNH	3	H23, H27	MG813594–MG813596
Portugal: Santa Comba de Rossas	NOG	2	H23, H46	MG813597–MG813598
Portugal: Santo Tirso	STS	2	H13, H28	MG813599–MG813600
Portugal: Braga	TIB	1	H21	MG813601
Portugal: Vila da Feira	FEI	1	H68	MG813602
Portugal: Lousada	LSD	2	H14, H40	MG813603–MG813604
Portugal: S. Pedro do Sul	SPS	2	H58	MG813605–MG813606
Portugal: Ribadouro	RIB	2	H23, H61	MG813607–MG813608
Portugal: Fafe	FAF	1	H23	MG813609
Portugal: Montemuro	MTM	6	H15, H47, H48, H55, H59, H60	MG813610–MG813615
Portugal: Amarante	AMA	3	H23, H38	MG813616–MG813618
Portugal: Ribeira de Pena	RPN	1	H69	MG813619
Portugal: Valongo	VAL	4	H23, H24, H28, H62	MG813620–MG813623
Portugal: Penacova	PNC	1	H45	MG813585
Portugal: Sever do Vouga	SVG	2	H66, H67	MG813624–MG813625
Portugal: Vila Pouca de Aguiar	VPA	3	H23, H38, H71	MG813626–MG813628
Portugal: Meda	MED	1	H57	MG813629
Portugal: Cernache do Bonjardim	BNJ	3	H55, H64	MG813630–MG813632
Portugal: Lousã	LOU	5	H31, H33, H63	MG813633–MG813637
Portugal: Arganil	ARG	2	H33	MG813638–MG813639
Portugal: S. Pedro de Moel	SPM	7	H55, H65	MG813640–MG813646
Portugal: Pomar, Alvito da Beira	MRD	2	H38	MG813648–MG813649
Portugal: Penalva do Castelo	PVC	1	H56	MG813647
Portugal: Castelo Branco	CBR	2	H34, H36	MG813650–MG813651
Portugal: Fundão	FUN	1	H32	MG813652
Portugal: Manteigas	MTG	1	H33	MG813653
Portugal: Guarda	GRD	1	H35	MG813654
Portugal: Malcata	MAL	5	H33, H49	MG813655–MG813659
Portugal: Vila de Rei	REI	2	H39, H55	MG813660–MG813661
Portugal: S. Mamede	SMA	8	H37	MG813662–MG813669
Spain: Pontevedra	PON	6	H4, H5, H6, H11	MG813670–MG813675
Spain: Santiago de Compostela	COM	4	H7, H11	MG813676–MG813679
Spain: A Coruña	COR	3	H1, H11	MG813680–MG813682
Spain: Lugo	LUG	5	H9, H11, H22	MG813683–MG813687
Spain: Donís	DON	1	H26	MG813688
Spain: Ponferrada	PNF	4	H11, H23	MG813689–MG813692
Spain: Ourense	OUR	2	H23	MG813693–MG813694
Spain: Monforte de Lemos	MNF	2	H8, H12	MG813695–MG813696
Spain: Toreno	TOR	2	H23	MG813697–MG813698
Spain: Navia	NAV	3	H2, H22	MG813699–MG813701
Spain: Oviedo	OVI	3	H3, H11	MG813702–MG813704
Spain: Peña de Francia	FRA	3	H50	MG813705–MG813707
Spain: Sierra de Guadalupe, Navezuelas	GUA	7	H29, H30	MG813708–MG813714
Spain: Álava	ALA	3	H10	MG813715–MG813717
Spain: Madrid	MAD	3	H53	MG813718–MG813720
Spain: El Tiemblo, Ávila	AVI	2	H53	MG813721–MG813722
Spain: Mengamuñoz, Ávila	AVL	1	H53	MG813723
Spain: Plasencia	PLA	2	H51, H52	MG813724–MG813725
Spain: Puebla de Sanabria	SAN	2	H17, H23	MG813726–MG813727
Spain: León	LEO	1	H20	MG813728

DNA amplification and sequencing

DNA was extracted from alcohol-preserved tissues following the standard proteinase K protocol (Sambrook and Russell 2001). Amplification of a 1,294-base pair fragment of cytochrome oxidase I (Cox1) was performed via PCR using the light chain primer P2F4 designed specifically for this work (5′-GGCCACTTTACCCGTGATATT-3′) and the heavy chain primer P3R (5′-ATRATATTTATTATYTGAGAAGC-3′) (San Mauro et al. 2004). Amplification conditions were as in Recuero et al. (2010), with annealing at 48°C. Sequences were determined in the ABI PRISM 3700 Genetic Analyzer. Sequences were aligned manually with BioEdit 7.0.9 (Hall 1999). The nucleotide sequences data reported here are registered in the GenBank Nucleotide Sequence Database under the accession numbers MG813567-MG813728.

Data analysis

We used PartitionFinder v. 1.0.1 (Lanfear et al. 2012) to select the best-fit model of nucleotide substitution considering 3 possible partitions for each codon position, using a greedy algorithm and the Bayesian information criterion. Accordingly, 3 independent partitions were defined for the first (HKY), second (TrN+G), and third (K80+I) positions.

Phylogenetic relationships were inferred using MrBayes v3.2 (Huelsenbeck and Ronquist 2001; Ronquist et al. 2012). Analyses were run for 10 million generations, sampling every 1,000 generations. Convergence was assessed through examination of the standard deviation of split frequencies, which was well below recommended thresholds. A consensus phylogram was computed after discarding trees reconstructed during the default burn-in period. This analysis was run in the Cipres Science Gateway (Miller et al. 2010). Maximum-likelihood (ML) analyses were performed with Garli v2.01 (Zwickl 2006) using the heuristic search algorithm with model parameters estimated with PartitonFinder. We used non-parametric bootstrapping (100 pseudoreplicates) to assess the stability of internal branches. We also explored our mtDNA data by estimating relationships between haplotypes using the program TCS version 1.21 (Clement et al. 2000). We used BEAST 1.8.4 (Drummond and Rambaut 2007) to estimate phylogenetic relationships, divergence times across clades, and demographic history of the different lineages, under a Bayesian Skyline tree prior (Drummond et al. 2005). A lognormal relaxed clock model was implemented preliminarily to test the fit of the data to a molecular clock model. As both ucld.mean and coefficient of variation priors tend to zero we cannot reject a strict molecular clock. Further analyses were performed under a strict molecular clock for the whole fragment. Lacking internal calibration points, time estimates were based on an uniform prior for the substitution rate, bounded at 0.012 and 0.014 substitutions/site/My, according to previous molecular studies in other European species of Rana (Canestrelli et al. 2014). Nucleotide substitution models were unlinked according to the PartitionFinder results. Analyses were run for 150 million generations, sampling every 15,000 generations. This analysis was run in the Cipres Science Gateway (Miller et al. 2010). We used Tracer v. 1.6 to check that effective sample sizes (ESSs) of parameters were >200 and convergence of parameter estimates across runs, and to build the final Bayesian skyline plots (BSP). A maximum clade credibility tree (MCCT) was reconstructed with TreeAnnotator.

In addition, pairwise mismatch distribution analyses (Rogers and Harpending 1992), Tajima’s D (Tajima 1989), and Fu’s Fs (Fu 1997) statistics were carried out in Arlequin V. 3.11 (Excoffier et al. 2005), to test for population expansion in the total data and within mitochondrial lineages. Mismatch distributions describe the frequency of pairwise substitutional differences among individuals. These distributions can be informative about the demographic history, as in populations that undergone a recent bottleneck and rapid expansions a unimodal mismatch distribution is expected (Rogers and Harpending 1992). Tajima’s D statistic tests for a departure from neutrality as measured by the difference between the number of segregating sites (θ) and the average number of pairwise nucleotide differences (π). Population expansions can cause significant negative departures of Tajima’s D from zero (Tajima 1989). Fu’s Fs statistic provides a test for population growth as it identifies an excess of rare alleles in an expanding population when compared with the number expected in a stationary population (Fu 1997). We carried out the analysis with 5,000 simulations.

The nucleotide diversity index (π) and the haplotype diversity (Hd) were used to estimate the genetic diversity for the overall data, within distinctive phyloclades, geographical subsets, and populations with enough sample size (n > 4) using DnaSP v. 5.10 (Rozas et al. 2003). Population structure was inferred by the analysis of molecular variance (AMOVA) (Excoffier et al. 1992) provided in the program package Arlequin V. 3.11 (Excoffier et al. 2005).

We also used BEAST v.1.8.4 to perform continuous diffusion analysis. Geographical coordinates for each sequence were included. A jitter module (window size: 0.01) generating random coordinates was included for samples from the same population. Analyses were run for 20,000 million generations under a strict molecular and a Yule tree prior. We used a Cauchy model to infer the geographical diffusion process through time across branches of the inferred gene tree (Lemey et al. 2010). Convergence was assessed using Tracer. A MCCT was reconstructed with TreeAnnotator and used to generate the reconstruction of the diffusion process using the modules “Continuous Tree” and “Time Slicer” in Spread 1.0.7 (Bielejec et al. 2011).

Contemporary and past niche modeling was performed with Maxent v.3.3.3k (Phillips et al. 2006; Phillips and Dudík 2008). We gathered a total of 510 localities from the whole distribution area, including the complete sampling plus personal observations and some of the GBIF records (many records were duplicated or based in obvious misidentifications; doi:10.15468/dl.qp7txp). Bioclimatic layers (30 arcsec), covering the range of the species, were downloaded from the WorldClim database (http://www.worldclim.org) and analyzed with ENMTools 1.4.3 (Warren et al. 2010) to check for correlations among them using the Pearson correlation coefficient. Past climate simulations followed the Community Climate System Model (http://www2.cesm.ucar.edu/). Taking into consideration the biological relevance of the different bioclimatic variables and the observed correlations (pairs of variables with Pearson coefficient value ≥0.8 were considered highly correlated), we considered Bio4 (Temperature Seasonality), Bio10 (Mean Temperature of Warmest Quarter), Bio12 (Annual Precipitation), Bio15 (Precipitation Seasonality), and Bio17 (Precipitation of Driest Quarter) to generate the model, using 70% of the localities as training data and 30% to test the model. The model was evaluated using “area under the curve” (AUC) values obtained for the test data (Swets 1988; Elith 2000).

Results

mtDNA structure and diversity

A total of 67 Cox1 variable sites and 71 haplotypes were detected among the 162 analyzed sequences. Three main matrilineal lineages were recovered by phylogenetic approximations (labeled A–C, Figure 2), with marked phylogeographic structure. With the exception of a single individual from the Serra de Montemuro (MTM), all samples in lineage A were found north of the Douro River, being the only one present in all populations from northern Spain, including those from the Basque Country. Lineage B was found across most of the range of R. iberica, mostly in the samples originated from Central Portugal and the Spanish Sistema Central Mountains. However, this lineage is also widely present in populations from Portugal north of the Douro River, in sympatry with lineage A. Lineage C is exclusive from the populations of Sierra de Guadalupe (GUA) in the southern limits of the species distribution in Spain (Figures 1 and 2).

Figure 2.

Mitochondrial diversity in R. iberica. Upper: Maximum clade credibility tree recovered by BEAST, showing time estimates for supported clades. Node bars represent 95% HPD intervals for node ages. Also shown are MrBayes/BEAST/ML posterior probabilities and bootstrap values for main nodes (only values >0.9 in at least one of the analyses or bootstrap values >50 are shown). The scale is in thousands of years. Lower: Maximum parsimony network for Cox1 haplotypes observed in R. iberica showing the 3 haplogroups recovered. Each circle represents a different haplotype (codes as in Table 1) with size proportional to its relative frequency. Small dots represent hypothetical unsampled haplotypes.

All 3 groups are supported as monophyletic, but the sister relationship between lineages A and C is only supported by MrBayes and ML, but not by BEAST analyses (Figure 2). Additionally, haplotypes of lineage B from the Sistema Central Mountains form a shallow but well-supported clade differentiated from the Portuguese populations. This clade is geographically restricted and only 1 individual from the Portuguese population of Malcata, just some tens of kilometers west of Spanish Sistema Central, falls within it. The topology of the haplotype network depicts 2 main star-type groups of haplotypes, corresponding to lineages A and B (Figure 2), with strong genetic structure within each clade. The lineage C exhibits very low genetic structure with only 2 haplotypes found. Age estimates provide wide intervals for all 3 main lineages and their most recent common ancestor (Figure 2), with differentiation likely occurring during the Middle to Upper Pleistocene.

Continuous diffusion analyses inferred an ancestral area located in the current border of north-western Portugal with Spain, with subsequent expansions to the North, followed by colonization of marginal areas in the Sistema Central Mountains, eastern Cantabric Region and Sierra de Guadalupe, and ending with a progressive expansion in Portugal into an admixture zone in northern Portugal (Figure 3).

Figure 3.

Continuous diffusion phylogeographical reconstruction in R. iberica, representing ancestral distribution and expansion of the identified mitochondrial lineages. (A) 280 Kyears before present (Kybp); (B) 200 Kybp; (C) 60 Kybp; and (D) 20 Kybp.

Historical demographic analyses

BSP results indicate demographic expansions for the whole species as well as for each of the 3 main mitochondrial lineages associated with the last glacial maximum (LGM), and a progressive stabilization in the Holocene (Figure 4). These results are generally in accordance with other historical demographic analyses: the mismatch analysis produced a bimodal distribution for the total data and unimodal distributions for all phylogroups and geographic subsets analyzed. Tajima’s D statistic detected significant deviations from neutrality in lineages A and B. With the exception of lineage C, Fu’s Fs statistics presented significant negative departures from zero in all 3 lineages and the total data set.

Figure 4.

Bayesian Skyline Plot (BSP) showing historical demographic trends in R. iberica; x-axis: time in thousands of years before present; y-axis: estimated population size [units = Net, the product of effective population size and generation length in years (log transformed)].

Mean pairwise sequence divergence among haplotypes varied from 0.080% to 1.050% and the overall nucleotide diversity (π) averaged 0.004. Lineage B exhibited the higher genetic diversity (h = 41, Hd = 0.958), followed by lineage A (h = 28, Hd = 0.887), and away down by lineage C (h = 2, Hd = 0.286). Nucleotide diversity is higher in populations from intermediate latitudes (Portugal north of the Douro River) (Figure 5), while no geographic structure was observed for Hd.

Figure 5.

Latitudinal distribution of nucleotide diversity in R. iberica along a South–North transect along the species main distribution area (Central Portugal to NW Spain).

For the AMOVA analysis, we grouped samples by lineages A–C. The AMOVA revealed significant genetic structuring across all hierarchical levels (all P < 0.0001), with 69.1% of the variation resulting of differences between clades (ΦCT = 0.691). Only 18.8% of the total variation was attributed to differences among populations within these clades (ΦSC = 0.609), while the remaining 12.1% was due to variation within populations (ΦST = 0.879).

Glacial and interglacial distribution

The models of climatic niche (Figure 6) showed high AUC test values (0.933, SD = 0.007) and significance for every binomial omission test, supporting a good performance of the model. The 5 included variables contributed unequally to the model, with 2 of them (Bio12, Bio17) accounting for 74.2% of relative contributions to the Maxent model. According to the generated models, the distribution of R. iberica has suffered a series of contractions and expansions associated to the climatic oscillations of the Pleistocene, with wider potential distribution during LGM and more restricted ones in interglacial periods, with relatively stable conditions in the North-west of the Iberian Peninsula.

Figure 6.

Climate niche models for R. iberica showing potential distribution of the species at (A) the present, (B) during the Last Glacial Maximum (LGM, 21 Ka), and (C) during the last interglacial period (LIG, up to ca. 120–140 Ka).

Discussion

Evolutionary history and phylogeographic structure

The Iberian Peninsula is one of the main biodiversity hotspots in the Western Palearctic, harboring not only a considerable number of endemic species (Médail and Quézel 1997), but also a large part of the intraspecific genetic diversity of widespread species is present in that region (Hewitt 2011). In many cases, this diversity reflects pre-Pleistocene evolutionary processes and has survived the numerous cycles of Pleistocene climatic oscillations in sanctuaries of diversity (Recuero and García-París 2011; Dufresnes et al. 2016) representing type “S” species (sensuRecuero and García-París 2011) such as the widespread species Rana temporaria Linnaeus (Vences et al. 2013, 2017).

In other cases, however, species have been severely affected by glacial and interglacial changes, presenting a comparatively reduced genetic diversity as only a small fraction of populations survived in true refugia (Recuero and García-París 2011; Martínez-Freiría et al. 2015). Rana iberica is an old lineage, sharing a Miocene common ancestor with a clade including R. temporaria, R. arvalis Nilsson, and R. pyrenaica Serra-Cobo (Yuan et al. 2016). However, the current mitochondrial phylogeographic pattern observed in R. iberica is the consequence of much recent evolutionary processes that took place only in the Middle to Upper Pleistocene, in direct association with the last 2 or 3 glacial periods. Similar patterns of reduced diversity, corresponding to type “R” species (sensuRecuero and García-París 2011) have also been observed in other species of Rana, including widespread ones such as Rana dalmatina (Vences et al. 2013) but also in other Mediterranean endemics. For example, in the Apennine Peninsula, R. italica Dubois shows 2 mitochondrial lineages with similar ages as in R. iberica (Canestrelli et al. 2008), while another Iberian endemic, R. pyrenaica, presents even higher diversity reduction with only 2 cytochrome b haplotypes differing in 1 transversion (Carranza and Arribas 2008).

Rana iberica, and probably most of Mediterranean populations of Rana, could be considered as relictual units currently surviving in interglacial refugia, as suggested also by our historical demographic analyses of R. iberica, which show that population expansions of every lineage during the last glacial period finished with, or little after, the beginning of the current postglacial era. The existence of a Pleistocene fossil record for this species in Burgos, dated at approximately 37,600 years BP (Esteban and Sanchíz 1991) and outside the current distribution of R. iberica, indicates that the species reached a wider range than its current one with subsequent extinction of peripheral populations. Of course, results presented here correspond to patterns of mitochondrial diversity and must be complemented with informative nuclear marker data to fully reconstruct the species evolutionary history, including recent demographic and dispersal events that may have not leave trace at the mitochondrial level (Zink and Barrowclough 2008).

Also, climate niche modeling for R. iberica suggests that its distribution was much extended, at least potentially, during glacial periods, spreading over a large part of the Iberian Peninsula. Potential distribution would be much reduced with interglacial conditions, confined into a range similar to that observed today. Successive events of range contractions, associated with strong demographic oscillations, provide an excellent scenario to promote the geographical isolation of more or less reduced population groups. Under such circumstances, processes of allopatric differentiation gave rise to the main lineages identified by our phylogenetic analyses. Even if more recent, this could be also the case of the Sistema Central sub-lineage. In this case, the species would have been isolated in a small refugium between Spain and Portugal, from where it expanded eastward, as suggested by the decrease of genetic diversity in that direction (Martínez-Solano et al. 2005), while westward expansion of this mitochondrial lineage would have been prevented by competitive exclusion (Waters 2011).

Alternatively, isolation by the presence of geological barriers to gene flow has been proposed in other species with overlapping distributions. The Douro River has been suggested to have acted as an important barrier for other amphibian species such as L. boscai and Alytes obstetricans (Laurenti 1768), with different phylloclades being observed northward and southward of this river (Fonseca et al. 2003; Martínez-Solano et al. 2006; Gonçalves et al. 2015). Additionally, Alexandrino et al. (2000) also indicate this river as a major barrier for C. lusitanica during its post-glacial expansion, resulting in depleted levels of genetic diversity for this species in northern Portugal and in the Spanish regions of Galicia and Asturias. However, the extensive presence of R. iberica lineage B north and south of the Douro River contradicts this hypothesis.

Lineages A and B are present in several populations in northern Portugal. This lineage combination results in the highest population nucleotide diversity across the entire R. iberica range. The most likely explanation for this is the formation of a diversity melting pot generated by admixture of the 2 lineages after range expansions, and this is supported by continuous diffusion analysis. This region is characterized by climatic niche modeling as one of the most stable areas for R. iberica, so the admixture zone may have formed early after lineage differentiation and maintained ever since without much change, allowing the accumulation and maintenance of intraspecific diversity, while its geographic limits were probably favored by competitive exclusion of mitochondrial genomes (Waters 2011).

Conservation implications

Conservation of relictual species with relatively large distribution ranges is problematic. The case of R. iberica exemplifies the situation. Despite of being an old lineage, intraspecific diversification in R. iberica is very limited and recent, suggesting that almost all units resulting from pre-Pleistocenic diversification processes along the evolutionary history of the lineage became extinct. The shallow genetic differentiation obtained for R. iberica is corroborated by previous allozyme and microsatellite studies, which detected low genetic variability among populations (DNEI = 0.000–0.020, Arano et al. 1993; FST = 0.145–0.250, Martínez-Solano et al. 2005). Evidence of a recent geographic distribution contraction, based on recent fossil records (Esteban and Sanchíz 1991), combined with a reduction in climatic suitability for amphibians and reptiles in the Iberian Peninsula (Araújo et al. 2006) makes the situation complicated. Teixeira and Arntzen (2002) predicted a reduction of favorable area for the distribution of the co-distributed and ecologically similar C. lusitanica up to 35% until the end of this century. The isolated populations of R. iberica at Serra de São Mamede (SMA), São Pedro de Moel (SPM), and Basque Country (ALA) exhibit very low genetic differentiation from the nearby populations, which suggest that they are very recent and that the species distribution is already retreating from its glacial optimum. To complicate matters, most of the mitochondrial genetic diversity of R. iberica is located in the southern part of its range, a scenario also described for several co-distributed species (Paulo et al. 2001; Alexandrino et al. 2002, 2007; Martínez-Solano et al. 2006; Teixeira et al. 2015). Lineage B exhibits by far the higher mitochondrial diversity of the species. Lineage A follows with a moderate diversity and lineage C exhibits a highly depauperate genetic diversity.

In terms of conservation efforts, it is clear that the southernmost population from the Sierra de Guadalupe deserves a special attention due to its very restricted distribution in the species rear edge, that is, marginal, isolated populations that have persisted through some of the Quaternary climate changes (Hampe and Petit 2005). In this respect, Hampe and Petit (2005) enhanced the potential critical importance of these populations as long-term stores of species’ genetic diversity and focuses of evolutionary and speciation potential. In addition, the tendency for climatic warming observed during the last decades (Hughes 2000) is expected to have a strong effect on amphibians (Pounds 2001) and to produce a poleward range shift of organisms if dispersal is unlimited (Parmesan and Yohe 2003; Araújo et al. 2006). Therefore, because of its peripheral Mediterranean location and its geographic isolation, the Sierra de Guadalupe population is prone to a higher susceptibility to climate warming and an increased risk of extinction. In this respect, the population of the Sierra de San Pedro raises an intriguing problem. Given the potential phylogeographical interest of this population isolate, at mid-distance between the population isolates of Sierra de Guadalupe (GUA) and SMA in Portugal, we realized several field trips to that area, but no larval or adult specimens were found. As so, we believe that this species does not currently occur in that mountain area, which is marked with a “?” in Figure 1. This fact can be explained either by an unlikely erroneous location of the species attributed for that area by Meijide (1985), or by an alarming evidence of local extinction of southern rear edge populations. Considering this second hypothesis, as this location constitutes the smallest population isolate with the lowest altitude amplitude from all rear edge populations, it would be the most vulnerable to local extinction events, and it highlights the risks for the rest of isolated, southern populations and particularly for the endemic mitochondrial lineage of Sierra de Guadalupe. Phylogeographic studies are a necessary tool to identify the existence of isolated lineages, with independent evolutionary trajectories from the main population core of otherwise genetically homogeneous taxa. Direct conservation measures are necessary to avoid losing these southernmost rear edge populations and unique lineages, because they might represent not only a critical fraction of the total amount of genetic diversity within the species, as it is the case for R. iberica, but also an important element of the species evolutionary potential.