To refer to this article use this url:

Contribution to Zoology, 75 (3/4) (2006)

Historical biogeography of Western Palaearctic pelobatid and pelodytid frogs: a molecular phylogenetic perspective

M. Veith1, L. Fromhage2, J. Kosuch3, M. Vences4

1.  Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam Kruislaan 318, 1098 SM Amsterdam, The Netherlands

2.  Biozentrum Grindel und Zoologisches Museum, Universität Hamburg, Martin-Luther-King Platz 3, 20146 Hamburg, Germany

3.  Universität Trier, Fachbereich VI Biogeographie, D-54286 Trier;

4.  Institute for Biodiversity and Ecosystem Dynamics, Zoological Museum, University of Amsterdam, Mauritskade 61, 1092 AD Amsterdam, The Netherlands

      current address: Zoological Institute, Technical University Braunschweig, Spielmannstr. 8, 38106 Braunschweig, Germany

Keywords: Phylogeny, paleobiogeography, speciation, regression-based molecular clock, Bayesian molecular dating, Messinian salinity crisis, 16S, 12S.


Spadefoot toads (Pelobates) and Parsley frogs (Pelodytes) are an enigmatic group of Western Palaearctic anurans. In the genus Pelobates, a fossorial lifestyle has enforced a conserved bauplan that masks their intraspecific evolutionary history. We used partial sequences of the mitochondrial 16S and 12S rRNA genes to infer a paleobiogeographic scenario of speciation events in these two anuran genera. Based on two alternative, mutually exclusive calibrations of the Iberian-African split within Pelobates (Pb. cultripes and Pb. varaldii), the disjunction of the Betic Cordillera ca. 14-16 million years ago (mya), and the end of the Messinian Salinity crisis 5.33 mya, we inferred alternative scenarios for species evolution within both genera applying regression-based dating and Bayesian molecular dating. Pelobates and Pelodytes are both monophyletic genera. Interspecific relationships among spadefoot toads are poorly resolved, and only an Iberian-African Pb. cultripes/Pb. varaldii clade consistently emerges from our analyses. An evolutionary scenario based on the Messinian divergence of African and Iberian Pelobates lineages becomes plausible in the light of geological and paleontological data. Consequently, Pelobates species are likely to have originated from the Miocene. Speciation around the Oligocene/Miocene boundary is inferred for the Iberian-Caucasian Pelodytes, and a Messinian divergence has to be invoked to explain intraspecific diversification of Iberian parsley frogs. There is indication that the different Pb. syriacus lineages may not form a monophylum.


Frog diversity is very unevenly distributed among major lineages. According to current species counts ( as of November 2005), the monophyletic frog crown group (Neobatrachia) contains more than 5000 species, whereas the basal lineages, subsumed in the Archaeobatrachia, include only 200 species. Whether archaeobatrachians are paraphyletic (e.g., Ford and Cannatella 1993) or monophyletic (e.g., Feller and Hedges 1998) has long been disputed, but recent studies using large datasets of nuclear DNA sequences point to a paraphyly of archaeobatrachians (Hoegg et al. 2004; San Mauro et al. 2005; Roelants and Bossuyt 2005). Biogeographically, several families of the Archaeobatrachia show an exclusively Holarctic distribution (such as Pelobatidae, Ascaphidae), with some being restricted to the Palearctic realm (Discoglossidae, Bombinatoridae, Pelodytidae, Megophryidae). Three of these families, the Megophryidae, Pelobatidae and Pelodytidae were formerly subsumed in a single family named Pelobatidae (e.g., Gilsen 1936). However, molecular data showed non-monophyly of the former Pelobatinae, resulting in the assignment of familial rank to each of these groups (e.g. Ford and Cannatella 1993; Hoegg et al. 2004), and to the Scaphiopodidae, containing the Nearctic spadefoot toads, Scaphiopus and Spea (e.g., Roelants and Bossuyt 2005). Hence, the Pelobatidae consist of a single genus, the Western Palearctic spadefoot toads ( Pelobates Wagler), similar to the Pelodytidae which contain a single genus, the Western Palaearctic parsley frogs ( Pelodytes ).

Within Pelobates , four species are currently recognised: Pb. fuscus (Laurenti, 1768) (with two subspecies, the widespread Pb. f. fuscus and the northern Italian endemic Pb. f. insubricus Cornalia, 1873), Pb. varaldii Pasteur and Bons 1959 (Morocco, Africa), Pb. cultripes Cuvier 1929 (mainly the Iberian Peninsula) and Pb. syriacus Boettger, 1889 (with four recognised subspecies, Pb. s. syriacus , Pb. s. boettgeri Mertens 1923, Pb. s. balcanicus Karaman 1928 and Pb. s. transcaucasicus Delwig 1928). Spadefoot toads are the most fossorially adapted European amphibians. They prefer sandy areas where they burrow themselves during the daytime, using their large and sharp metatarsal tubercles.

Systematic relationships among Pelobates species have recently been studied by García-Paris et al. (2003). They showed that a Pb. cultripes-Pb. varaldii clade is sister to a Pb. fuscus - Pb. syriacus clade. However, their sampling did not allow for detailed assessment of intraspecific phylogeny, and their approach did not include a molecular dating of speciation events in Pelobates or Pelodytes .

The intrageneric taxonomy of Pelobates has changed several times, mainly with regard to the eastern Pb. syriacus group. The validity of the Pb. syriacus subspecies has been repeatedly questioned (e.g. Eiselt and Schmidtler 1973, Terentiev and Chernov 1965, Kuzmin 1999), and Eiselt (1988) even considered Pb. syriacus a monotypic species. Morphological evolutionary stasis, probably linked to their fossorial lifestyle, may account for the lack of diagnostic morphological characters among Pelobates taxa (e.g. between Pb. varaldii and Pb. cultripes ; Busack et al. 1985). However, based on a morphological analysis, Ugurtas et al. (2002) recently emphasised that Anatolian and European Pb. syriacus comprise at least three taxa: Pb. s. syriacus , Pb. s. balcanicus and a third lineage that they collected in Serbia. In addition, Borkin et al. (2001, 2003) discovered two well-differentiated Pb. f. fuscus lineages with non-overlapping genome sizes (measured by the total amount of nuclear DNA per cell). This finding was supported by allozyme data, so they regarded them as differentiated at the species level (provisionally named the “western” and the “eastern” form of P. f. fuscus ). Again, they attributed this cryptic species diversity to the morphological stasis of spadefoot toads.

Parsley frogs live in all kinds of stagnant and slow-moving waters of the Iberian Peninsula and the Caucasus mountains (Gasc et al. 1997; in contrast to the Iberian species the Caucasian parsley frog prefers habitats with semi-current water). Their distribution is disjunct. While Pd. caucasicus (Boulenger, 1896) is restricted to the Caucasus, two species, Pd. ibericus Sánchez-Herráiz, Barbadillo, Machordom and Sanchiz, 2000 and Pd. punctatus (Daudin, 1803) are originally Iberian endemics, with the latter having dispersed postglacially into large parts of France (Gasc et al. 1997). This Caucasus-Iberian disjunction is mirrored by the almost identical distribution of Mertensiella caucasica and Chioglossa lusitanica , two sister taxa within the Salamandridae (Titus and Larson 1995, Veith et al. 1998).

In the present paper we use molecular sequence data of mitochondrial genes to analyse the intrageneric phylogenetic relationships of Pelobates and Pelodytes . Based on two different calibrations of a molecular clock, we infer paleogeographic scenarios of their evolution in the circum-Mediterranean region.

Materials and methods

We sequenced 34 specimens from eight taxa of Pelobates and Pelodytes (see appendix). For hierarchical outgroup rooting we included Scaphiopus (from GeneBank), several archaeobatrachian representatives (Pipidae, Megophryidae, Discoglossidae, Leiopelmatidae and Ascaphidae; partially from GenBank), two Neobatrachian species of the Ranidae (GenBank) and Bufonidae (GenBank), and the newt Triturus vulgaris (GenBank).

DNA extraction, sequencing and sequence alignment

DNA was extracted from ethanol-preserved tissue samples using the QiAmp tissue extraction kit (Qiagen). We amplified via polymerase chain reaction (PCR) fragments of two mitochondrial genes coding for fractions of the 16S rRNA and 12S rRNA (Tab. 1).

16S: 16SA (light chain; 5’ - CGC CTG TTT ATC AAA AAC AT - 3’) and 16SB (heavy chain; 5’ - CCC GTC TGA ACT CAG ATC ACG T - 3’) of Palumbi et al. (1991) amplified a ca. 580 bp section of the mitochondrial 16S rRNA gene homologous to positions 3976-4554 of the Xenopus laevis mitochondrial genome (Roe et al. 1985).

12S: 12SA-L (light chain: 5’ - AAA CTG GGA TTA GAT ACC CCA CTA T - 3’) and 12SB-H (heavy chain: 5’ - GAG GGT GAC GGG CGG TGT GT - 3’) of Goebel et al. (1999) amplified a ca. 490 bp section of the mitochondrial 12S rRNA gene homologous to positions 2510-2997 of the Xenopus laevis mitochondrial genome (Roe et al. 1985).

Polymerase chain reaction cycling procedures were as follows: 16S rRNA gene: denaturation for 45 s at 94°C, primer annealing for 45 s at 55°C, extension for 60 s at 72°C; after 35 cycles final step at 72°C for 10 min.; 12S rRNA gene: denaturation for 45 s at 92°C (initial denaturation for 120 s at 94°C), primer annealing for 60 s at 50°C, extension for 90 s at 72°C; after 35 cycles final step at 72°C for 10 min.

PCR products were purified using the ”High Pure PCR Product Purification Kit” (Roche diagnostics). We sequenced single-stranded fragments using the “Big dye terminator cycle sequencing kit” of Applied Biosystems Inc. on an ABI 377 automatic sequencer using standard protocols.

Sequences were aligned automatically using the Clustal X software for Windows, version 1.81 (Thompson et al., 1997). Alignment penalties were arbitrarily set to 12 for gap opening and to 5 for gap extension. Redundant haplotypes were excluded from analyses.

Phylogenetic analyses

We tested for partition (=gene) combinability (Huelsenbeck et al., 1996) using the maximum parsimony procedure of Farris et al. (1994; incongruence of length differential test, ILD) as implemented in Paup* (Swofford, 2001; 100 replicates, heuristic search using the tree bisection-reconnection (TBR) branch-swapping algorithm). Recently, the ILD test has been criticised (e.g., Yoder et al., 2001; Barker and Lutzoni, 2002). We therefore also analysed each gene separately and compared the resulting topologies with the combined analysis.

To avoid topology formation by variation in base composition among taxa (Steel et al., 1993) we applied a χ2-test as implemented in PAUP*. Unequal base frequencies among taxa would necessitate a LogDet transformation of sequence differences that allows consistent recovery of the correct tree when sequences evolve under simple asymmetric models that can vary between lineages (Lockhart et al., 1994). Equal base distribution among taxa would allow application of a specific substitution model, as was recently recommended by Sullivan and Swofford (2001).

Additionally, we applied a hierarchical likelihood ratio test for the goodness-of-fit of nested substitution models (ingroup taxa only; for details see Huelsenbeck and Crandall, 1997) using the program Modeltest version 3.02 of Posada and Crandall (1998). To test for the possibility that some types of nucleotide substitutions have become saturated, we plotted uncorrected p distances versus molecular distances derived from the specific substitution model (transitions and transversions, separately for each gene).

We defined the newt Triturus vulgaris as the outgroup and subjected our alignment to four different methods of phylogenetic reconstruction: (i) neighbor joining; (NJ; Saitou and Nei, 1987) using the specific substitution model; (ii) maximum parsimony with gaps treated as fifth character state (within ingroups, gaps occurred only in the 16S rDNA and mostly consisted of single indels); 10 random additions of haplotypes; transitions and transversions were given equal weight; heuristic search with the TBR branch swapping algorithm; (iii) maximum likelihood (ML; Felsenstein, 1981) based on the specific substitution model, (iv) Bayesian inference (Rannala and Yang, 1996, Huelsenbeck et al., 2001), which is based upon the notion of posterior probabilities of a phylogenetic tree. With the exception of the Bayesian approach (MrBayes; Huelsenbeck and Ronquist, 2001), all analyses were done with Paup* (Swofford, 2001).

Robustness of NJ and MP tree topologies was tested by bootstrap analyses (Felsenstein, 1985), with 2000 replicates each (Hedges, 1992). Due to computational constraints, we used Quartet Puzzling (Strimmer and von Haeseler, 1996) with 100,000 permutations to infer reliability values for ML tree topologies. To gain confidence in Puzzle support values we compared them to ML bootstrap values derived from 100 replicates. We applied the Bayesian method using the general time reversible model of nucleotide substitution (GTR; Rodríguez et al., 1990) with the proportion of invariable sites estimated from the data. We ran four simultaneous Metropolis-coupled Monte Carlo Markov chains for 500,000 generations. We repeated Bayesian analysis four times and plotted likelihood values against generation number to ensure that the likelihood converged at the same value (Leaché and Reeder, 2002) and to determine the necessary burn-in. All runs with the Bayesian approach converged at a likelihood of ca. -lnL =4660, indicating that this marks a global optimum. We sampled a tree every 100 generations and calculated Bayesian posterior probabilities for 4000 trees by omitting the first 1000 trees (burn-in).

For simplicity, we use the term “support values” when referring to bootstrap values for NJ and MP, puzzle support values for ML and Bayesian posterior probabilities. We regard bootstrap proportions and Puzzle support values of ≥70% as indicators for sufficiently resolved topologies since they usually correspond to a probability of ≥95% that the corresponding clade is valid (Hillis and Bull 1993, but see constraints of this interpretation discussed therein). We accept Bayesian posterior probabilities ≥95% as indicators of sufficient node support, since they correspond to a ≤5% probability of committing a type I error (i.e. of drawing an erroneous conclusion).

Molecular clock calibration

We tested for rate constancy between Pelobates and Pelodytes using a likelihood ratio test (Tree Puzzle; Schmidt et al., 2002) and defining Discoglossus montalentii as the single outgroup (TrN substitution model with eight gamma rate categories). This test compares the log-likelihood of the most likely tree with and without a molecular clock enforced. In addition we applied the Tajima χ2-test (Tajima, 1993) for substitution rate heterogeneity among all possible 153 pairs of ingroup haplotypes (including Scaphiopus ), again using D. montalentii as the outgroup.

We used the specific substitution model to calculate pairwise corrected molecular distances among haplotypes. We calibrated two molecular clocks. Calibration I is based on the hypothesis that the population of the ancestor of Ibero-African Pelobates was separated by the formation of the Strait of Gibraltar (Busack, 1986; García-Paris et al., 2003) into the Iberian Pelobates cultripes and the African P. varaldii . This event was precisely dated to 5.33 mya, when the re-flooding of the Mediterranean basin marked the end of the Messinian Salinity Crisis (Krijgsman et al., 1999). Calibration II is based on the Betic Crisis. This vicariance event has often been invoked to explain divergent evolution of African and Iberian lineages (García-París and Jockusch 1999, Fromhage et al. 2004, Veith et al. 2004). During the Betic Crisis 16-14 mya, marine transgressions through the Betic Strait (corresponding to today’s Guadalquivir basin in southern Iberia) separated the Betic-Rif mountain belt, which developed during, and partly in response to, Late Mesozoic and Tertiary convergence between Africa and Iberia (Lonergan and White, 1997), from the Iberian mainland. Promoted by the further orogenic uplift of the Alboran basin between Iberia and Africa, the southernmost part of the insular Betic region connected to the African continent and formed today’s Rif Mountains (De Jong, 1998).

We applied the bootstrap method to compute variances for molecular distances between two sequences. This requires no assumption about the underlying distribution of evolutionary distances except that each site evolves independently, an assumption that is usually met when the number of sites is >100 (Nei and Kumar, 2000). We used the software MEGA, version 2.1 (Kumar et al., 2001), to estimate times of divergence of splits. 95% confidence limits of estimates of divergence time were calculated via 1000 bootstrap replicates.


Fig. 1. Maximum parsimony tree of Western Palearctic Pelobatidae and Pelodytidae based on 809 bp of the 16S and 12S rRNA genes (TrN substitution model); support values >50% are given for MP (upper), NJ (upper half), ML using quartet puzzling (lower half) and Bayesian inference (lower).

Since estimation of rate constancy using the likelihood approach as implemented in Tree Puzzle may strongly depend on the choice of outgroup (in our case some possible outgroups resulted in a rejection of the assumption of rate constancy) we also applied Bayesian molecular dating (Thorne et al., 1998) using the multidivtime package of Thorne and Kishino (2002). This approach accounts for rate heterogeneity among lineages. We consistently used default settings of multidivtime as recommended by Rutschmann (2004). Again we defined two calibrations: the vicariance of Europe from Africa following the end of the Messinian salinity crisis 5.33±0.02 million years ago (Krijgsman et al. 1999), and the Betic cordillera disjunction between 16 and 14 mya (Lonnergan & White, 1997). To avoid biased estimates of divergence times due to overrepresentation of some lineages (Yoder & Yang 2004) we only included one sequence per species, except for P. fuscus , where we included a representative of each genome size form.


Alignment statistics and model selection

According to the ILD test, data partitions contained no incongruent phylogenetic signals, so we combined 12S and 16S fragments into a single alignment. It consistently contained 809 bp for all specimens, with 333 variable and 245 parsimony informative sites (143 and 98, respectively, when regarding ingroup taxa only). Empirical base frequencies were πA = 0.356, πC = 0.233, πT = 0.249 and πG = 0.162. The strong anti-G bias indicates that no nuclear pseudocopies of the genes have been analysed (Zhang and Hewitt, 1996). Bases were homogeneously distributed among ingroup haplotypes (χ2-test: p = 1.00). The substitution model that fitted our alignment best was the Tamura-Nei model (Tamura and Nei, 1983) with a proportion of invariable sites of I = 0.3674, a gamma distribution shape parameter α = 0.6172 and substitution rate= 3.0767 and rate= 6.1846 (-lnLTrN+I+G = 4690.2271).

Clade formation

Only a few topologies were sufficiently supported by all tree-building approaches (Fig. 1): (i) monophyly of each of Pelobates , Pelodytes , Leptobrachium and Leptolalax ; (ii) monophyly of Pb. varaldii and Pb. cultripes ; (iii) monophyly of Pd. ibericus and Pd. punctatus .

Basal split among families are poorly resolved. Scaphiopus stands together with Pelodytes , although with low support values. Within Pelobates , two well-differentiated lineages of Pb. syriacus (ssp. transcaucasicus and balcanicus ) emerged, with no support for monophyly of the species.

Molecular clock calibration and application

Plots for uncorrected distances against corrected TrN+I+G distances (Fig. 2) showed signs of saturation for neither transitions (ti’s) nor transversions (tv’s). The log-likelihood ratio test with D. montalentii as outgroup indicated rate constancy among all taxa of Pelobates and Pelodytes . The log-likelihood of the clocklike tree (-logL=2332.95) was not significantly smaller than the log-likelihood of the tree with no clock enforced (-logL=2342.67; the simpler (clocklike) tree is therefore not rejected on a significance level of 5%). In addition, no pair of ingroup haplotypes showed a significant deviation from substitution rate constancy when applying the Tajima test.

Calibrating a molecular clock based on molecular distances and using the end of the Messinian Salinity Crisis (calibration I) for the vicariance of Pb. cultripes and Pb. varaldii results in mean separation times among basal Pelobates lineages of 12.12 my (Tab. 1). Separation of the Eastern and Western genome size lineages of Pb. f. fuscus are estimated to the Pleistocene ca. 1.69 mya. However, contemporary separation is within their lower 95% CI’s. The split between Iberian Pelodytes is dated exactly to the end of the Messinian Salinity Crisis 5.33 mya, while separation of Iberian and Caucasian Pelodytes lineages is much older (22.3 mya). Divergence time estimates based on Bayesian methods (Tab. 2) are not in line with those based on molecular distances (Tab. 1). The split between the eastern and western genome size lineages of Pb. fuscus is estimated to 2.4 mya, while the basal split among Pelobates is dated to 6.8 mya. Estimates among Pelodytes lineages are much younger under a Bayesian framework, with the intra-Iberian split estimated to only 1.15 mya and the Iberian-Caucasian separation being estimated to the late Pliocene 2.87 mya.


Fig. 2. Saturation plots of uncorrected distances against TrN+I+G molecular distances for different genes and portions of genes; deviation from the bisector of an angle indicates the degree of saturation.


Phylogeny of Pelobatidae and Pelodytidae

The phylogenetic hypotheses (Fig. 1) do not provide an adequate resolution of the basal splits among archaeobatrachian representatives. In contrast to García-Paris et al. (2003) but in agreement with the nuclear DNA sequence analysis of Hoegg et al. (2004) and Roelants and Bossuyt (2005), our data do not support a Mesobatrachia clade that includes the Pipidae, since our representative, Pipa parva , clusters outside a (Megophryidae, Pelobatidae, Pelodytidae, Scaphiopodidae) clade. While the monophyly of each included family receives moderate to high support, and relationships within Pelodytes are well supported, the topology within Pelobates was not consistent among tree building approaches, resulting in only weak support for several of the clades shown in Fig. 1.

Reliability of the molecular clocks

In our dataset, a clocklike behaviour of the DNA sequence evolution was not rejected. This is a privileged situation in amphibians, where large differences in substitution rates among clades are common, especially in mitochondrial genes (Hoegg et al. 2004). This pattern also allows for the application of a regression-based molecular clock approach, in which calibration times are plotted against pairwise distances among taxa, and the age of unknown splits is deduced from the regression slope (e.g., Hillis et al. 1996). It is remarkable that this method, in the present dataset, produces results largely inconsistent with those of the Bayesian approach that is able to handle differences in substitution rates among lineages (Thorne et al. 1998). It needs further exploration how good this method is able to handle datasets of low to very low divergences among taxa, as in pelobatid and pelodytid frogs. Similar cautions may apply to likelihood methods of phylogeny reconstruction; Fromhage et al. (2004) were unable to recover plausible phylogenetic topologies for closely related discoglossid frogs using ML methods. Overparametrization and, consequently, usage of too complex nucleotide substitution models may explain the failure of likelihood-based methods in these examples, which is why we base the following discussion on the regression-based molecular clock estimates using pairwise distances only.

Which paleogeographic scenario can best explain Pelobates and Pelodytes evolution?

Estimated times of divergence of our calibrations differ by a factor of ca. 2.6. In a simplistic approach to molecular clock calculations, and considering the uncorrected pairwise divergences between Pelobates cultripes and Pb. varaldii of 0.9%, calibration I would result in a rate of pairwise sequence divergence of approximately 0.17%/million years (my), whereas calibration II would result in a rate of 0.06%/my. The latter rate is much lower than rates usually assumed for mitochondrial rRNA genes in vertebrates (e.g., Veith et al. 1998: 0.7%/my in salamanders; Vences et al. 2001: between 0.100 and 0.500%/my in fish and mammals). This would argue in favour of calibration I, however, from biogeographical and geological points of view, there is support for and contradiction to both scenarios.

Assuming a Pliocene divergence of Pb. varaldii and Pb. cultripes 5.33 mya, calibration I in the regression approach dates the major splits among Pelobates to the Miocene. A separation of an Ibero-African lineage ( Pb. cultripes and Pb. varaldii ) from an E European-W Asian lineage ( Pb. fuscus and Pb. syriacus ) can be explained by the middle Miocene disjunction of land masses that formerly separated the Tethys and Paratethys ocean (biogeographic characteristic no. 4 of Oosterbroek and Arntzen 1992). However, assuming this we have to invoke the upper 95% confidence limits of our time estimates. Separation within the E European-W Asian lineage (split 2 in Tabs. 1 and 2) fits the onset of the first glaciation cycles at ca. 3.4-2.5 mya (Wilson et al. 2000) only if we invoke the respective lower 95% confidence limit (see also Veith et al. 2003a,b for the effect of Pleistocene glaciations on amphibian speciation in the Western Palearctic). Problems occur with the explanation of the intra-Iberian speciation of Pelodytes (split 5) since the end of the Messinian Salinity Crisis was not linked to any known Intra-Iberian disjunction. The current distribution of Pd. punctatus and Pd. ibericus instead corresponds to a pattern found in other pairs of vicariant amphibian taxa (e.g. in Alytes , Discoglossus , and Salamandra ; Arntzen and Garcia-Paris 1997, García-Paris et al. 1998, García-Paris and Jockusch 1999), which are often explained by the disjunction of the Betic Cordillera from the Iberian mainland through the formation of the Betic Sea Strait ca. 14 mya. However, a second reopening of the Betic Strait ca. 10-6 mya (López Martínez 1989 in Martínez-Solano et al. 2004) falls well into the 95% CI of split 5. In addition, it needs to be taken into account that data currently being assembled by other research groups support the presence of a further Pelodytes lineage, probably an undescribed species, in the Iberian Peninsula and the available distribution data for Pd. ibericus probably subsume two lineages of yet unknown relationships to each other. The time estimate for the Iberian-Caucasian vicariance within Pelodytes , according to our data, was estimated to the Oligocene/Miocene boundary. It matches well the time period of 15-20 mya when Mertensiella caucasica and Chioglossa lusitanica were assumed to have diverged (Veith et al. 1998; but see Steinfartz et al. 2000). Although also this estimate was deduced from molecular data (16S rDNA), it is plausible to attribute the evolution of a pair of vicariant amphibian species with such a strange but identical disjunct distribution to the same vicariance event.


Table 1: Regression-based time estimates of divergence among lineages of Pelobates and Pelodytes, based on two molucelar clock calibrations for the split Pelobates cultripes - Pb. varaldii and applying molecular distances; standards errors of TrN distances were calculated via 1000 bootstrap replicates; mean±1.96 S.E. resulted in lower and upper 95% confidence limits.

Calibration II may easily explain split 5. Separation of the former Betic region from the Iberian mainland 14 mya ago may have left behind ancestors of both Pb. varaldii and Pd. ibericus on the Betic Cordillera. The southern part, today’s Rif mountains, that carried Pb. varaldii, drifted southwards and subsequently connected to the African continent, while the north-eastern part, harbouring Pd. ibericus , reconnected to the Iberian mainland after the Betic Sea Strait closed again. However, if this scenario holds, two extinction events on both parts of the formerly isolated Betic mountain chain have to be assumed: extinction of the ancestral Pb. varaldii on its north-eastern part, and extinction of Pd. ibericus on its southern part. Evolution of major Pelobates lineages may be explained by fragmentation of today’s E European/Asian Minor landmasses in the course of the restoration of marine conditions between the Tethys and Paratethys during the early Oligocene (see Oosterbroek and Arntzen 1992[b). Split 4 can easily be explained by any appropriate event during the last ten million years (see its broad CI ranging from 0 to almost 10 my), including the onsets of Pliocene and Pleistocene glaciation cycles; in a strict sense this makes it almost entirely non-informative. Problems arise with the Iberian-Caucasian separation within Pelodytes more than 60 mya. Even its lower 95% CI by far predates the above mentioned time of divergence of 15-20 mya between the Iberian-Caucasian sister taxa Mertensiella caucasica and Chioglossa lusitanica and is not supported by paleontological data (Sanchiz 1998, Rage & Rocek 2000) according to which Pelodytes are not known from periods earlier than the Eocene. An extinct species of Pelodytes ( P. arevacus ) is known from Miocene deposits of Spain and belongs to a clade with P. punctatus and P. ibericus (Sanchiz 1998; Sanchiz et al. 2002) or may even be conspecific with P. punctatus (Rage & Rocek, 2000). One may also assume a pre-Miocene extinction of intermediate linking Pelodytes populations between the Caucasus and Iberia. Alternatively, we could invoke an Early Paleocene disjunction of landmasses in the Tethys Ocean (Oosterbroek and Arntzen 1992) to explain the evolution of eastern and western Pelodytes lineages.


Table 2: Bayesian molecular dating of times of divergence among lineages of Pelobates and Pelodytes, based on two molecular clock calibrations for the split Pelobates cultripes - Pb. varaldii. Much older estimates of divergence times result if calibration II is invoked. If Pb. varaldii and Pb. cultripes are of early Miocene origin (14 mya), most other lineages should have formed in the middle and late Miocene ca. 20 and 35 mya (distance based estimates). The split between eastern and western Pb. fuscus lineages would be of pre-Pleistocene origin (4.4 mya). Caucasian and Iberian Pelodytes would have diverged ca. 60 mya. Again, Bayesian time estimates largely disagree with distance based ones.

Taxonomic implications

Based on an evolutionary species concept several taxonomic implications emerge from our analyses:

(1) The monophyly of Pb. syriacus is not unambiguously supported. Non-monophyly would be a possible indication for the existence of two species. The degree of differentiation of Pb. syriacus clade 1 ( Pb. s. transcaucasicus ) and Pb. syriacus clade 2 ( Pb. s. balcanicus ) would justify assignment of species rank to both of them. This is corroborated by the sympatric occurrence without hybridisation of other Pelobates taxa that are differentiated at a similar level: Pb. cultripes and Pb. f. fuscus in France (Lescure 1984), and Pb. f. fuscus and Pb. syriacus in different areas (Eiselt 1988, Eggert et al. 2006).

(2) The eastern and western genome size types of Pb. f. fuscus are not differentiated on a level typical for other Pelobates species. Consequently, our molecular data do not support species rank assignment to them as was suggested by Borkin et al. (2001, 2003). A final conclusion must await data on the importance of the karyotypical characters for preventing viable and fertile hybridization between these forms. Final taxonomic conclusions in this complex, however, need to await a comprehensive phylogeographic sampling of Pd. syriacus, including all four currently accepted subspecies.


We are grateful to T. Amann, F. Andreone, W. Bischoff, L. Borkin, M. Franzen, F, Glaw, S. Litvinchuk, U. Joger, A. Kiefer, K.-D. Kühnel, R. Malkmus, B. Samraoui and F.J. Schmidtler for providing tissue samples.


Arntzen JW, García-París M. 1997. Phylogeny and biogeography of midwife toads ( Alytes , Discoglossidae): a rebuttal. Contrib. Zool. 66: 263-268.

Barker FK, Lutzoni FM. 2002. The utility of the incongruence length difference test. Syst. Biol. 51: 625-637.

Borkin LJ, Litvinchuk SN, Rosanov JM, Milto KD. 2001. Cryptic speciation in Pelobates fuscus (Anura, Pelobatidae): evidence from DNA flow cytometry. Amph.-Rept. 22: 387-396.

Borkin LJ, Litvinchuk SN, Rosanov JM, Khalturin MD, Lada GA, Borissovsky AG, Faizulin AI, Kotserzhinskaya IM, Novitsky RV, Ruchin AB. 2003. New data on the distribution of two cryptic forms of the common spadefoot toad ( Pelobates fuscus) in eastern Europe. Russian J. Herpetol. 10 (2): 115-122.

Busack SD. 1986. Biogeographic analysis of the herpetofauna separated of the Strait of Gibraltar. Nat. Geographic Res. 2: 17-36.

Busack SD, Maxson LR, Wilson MA. 1985. Pelobates varaldii (Anura: Pelobatidae): a morphologically conservative species. Copeia 1985:107-112.

De Jong H. 1998. In search of historical biogeographic patterns in the western Mediterranean terrestrial fauna. Biol. J. Linn. Soc. 65: 99-164.

Eggert C, Cogalniceanu D, Veith M, Dzukic G, Taberlet P. 2006 in press. The declining Spadefoot toad, Pelobates fuscus (Pelobatidae): Recent environmental changes as a major influence on current population structure and status. Conserv. Genet.

Eiselt, J. 1988. Krötenfrösche ( Pelobates gen., Amphibia, Salientia) in Türkisch-Thrakien und Griechenland. Ann. Naturhist. Mus. Wien 90B:51-59.

Eiselt J, Schmidtler JF. 1973. Froschlurche aus dem Iran unter besonderer Berücksichtigung außeriranischer Populationen. Ann. Naturhist. Mus. Wien 77: 181-243.

Farris JS, Källersjö M, Kluge AG, Bult C. 1994. Testing significance of incongruence. Cladistics 10: 315-319.

Feller AE, Hedges SB. 1998. Molecular evidence for the early history of living amphibians. Mol. Phyl. Evol. 9: 509-516.

Felsenstein J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17: 368-376.

Felsenstein J. 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39: 783-791.

Ford LS, Cannatella DC. 1993. The major clades of frogs. Herpetol. Monographs 7: 94-117.

Fromhage L, Vences M, Veith M. 2004.Testing alternative vicariance scenarios in Western Mediterranean discoglosside frogs. Mol. Phyl. Evol. 31: 308-322.

García-París M, Alcobendas M, Alberch P. 1998. Influence of the Guadalquivir river basin of mitochondrial DNA evolution of Salamandra salamandra (Caudata: Salamandridae) from Southern Spain. Copeia 1:173-176

García-París M, Jockusch EL. 1999. A mitochondrial DNA perspective on the evolution of Iberian Discoglossus (Amphibia: Anura). J. Zool. 248(2): 209-218.

García-París M, Buchholz DR, Parra-Olea G. 2003. Phylogenetic relationships of Pelobatoidea re-examined using mtDNA. Mol. Phyl. Evol. 28: 12-23.

Goebel AM, Donelly JM, Atz ME. 1999. PCR primers and amplification methods for 12S ribosomal DNA, the control region, cytochrome oxidase I, and cytochrome b in bufonids and other frogs and an overview of PCR primers which have amplified DNA in amphibians successfully. Mol. Phyl. Evol. 11: 163-199.

Gasc J-P, Cabela A, Crnobrnja-Isailovic J, Dolmen D, Grossenbacher K, Haffner P, Lescure J, Martens H, Martínez-Rica JP, Maurin H, Oliveira ME, Sofianidou TS, Veith M, Zuiderwijk A. 1998. Atlas of amphibians and reptiles in Europe.Paris: Societas Europaea Herpetologica and Muséum national d’Histoire naturelle.

Hedges SB. 1992. The number of replications needed for accurate estimation of the bootstrap P value in phylogenetic studies. Mol. Biol. Evol. 9: 366-369.

Hillis DM, Bull JJ. 1993. An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis. Syst. Biol. 42(2): 182-192.

Hillis DM, Mable BK, Moritz C. 1996. Applications of molecular systematics. In: Hillis DM, Moritz C, Mable BK, eds. Molecular Systematics, 2nd ed., Sinauer Associates, Sunderland, MA, 515-543.

Hoegg S, Vences M, Brinkmann H, Meyer A. 2004. Phylogeny and comparative substitution rates of frogs inferred from three nuclear genes. Mol. Biol. Evol. 21: 1188-1200.

Huelsenbeck JP, Crandall KA. 1997. Phylogeny estimation and hypothesis testing using maximum likelihood. Ann. Rev. Ecol. Syst. 28: 437-466.

Huelsenbeck JP, Bull JJ, Cunningham CW. 1996. Combining data in phylogenetic analysis. TREE 11: 152-158.

Huelsenbeck JP, Ronquist FR. 2001. MrBayes: Bayesian inference of phylogenetic trees. Bioinformatics 17: 754-755.

Huelsenbeck JP, Ronquist F, Nielsen R, Bollback JP. 2001. Bayesian inference of phylogeny and its impact on evolutionary biology. Science 294: 2310-2314.

Leaché AD, Reeder TW. 2002. Molecular systematics of the Eastern fence lizards (Sceloporus undulates): a comparison of parsimony, likelihood, and Bayesian approaches. Syst. Biol. 51: 44-68.

Lescure J. 1984. Répartition des Pelobates en France au XIXe et XXe siècle. Bull. Soc. Herpetol. France 29:45-59.

Lockhart PJ, Steel MA, Hendy MD, Penny D. 1994. Recovering evolutionary trees under a more realistic model of sequence evolution. Mol. Biol. Evol. 11: 605-612.

Lonergan L, White N. 1997. Origin of the Betic-Rif mountain belt. Tectonic 16: 504-522.

Krijgsman W, Hilgen FJ, Raffi I, Sierro FJ, Wilson DS. 1999. Chronology, causes and progression of the Messinian salinity crisis. Nature 400: 652-655.

Kumar S, Tamura K, Jakobsen IB, Nei M. 2001. MEGA2: Molecular Evolutionary Genetics Analysis software. Tempe, Arizona, USA: Arizona State University.

Kuzmin SL. 1999. The amphibians of the former Soviet Union. Sofia/Moscow: Pensoft.

Nei M, Kumar S. 2000. Molecular evolution and phylogenetics. Oxford: Oxford Univ. Press.

Palumbi S, Martin A, Romano S, McMillan WO, Stice L, Grabowski G. 1991. The simple fool’s guide to PCR, version 2. Honolulu, Hawai.

Posada D, Crandall KA. 1998. MODELTEST: testing the model of DNA substitution. Bioinformatics 14: 817-818.

Rage J-C, Rocek Z. 2000. Tertiary anura of Africa, Asia, Europe and North America. In: Heatwhole H, ed. Amphibian Biology, volume 4: Palaeontology. Chipping Norton: Surrey Beatty & Sons.

Rannala B, Yang Z. 1996. Probability distribution of molecular evolutionary trees: a new method for phylogenetic inference. J. Mol. Evol. 43: 304-311.

Rice WR. 1989. Analyzing tables of statistical tests. Evolution 43: 223-225.

Rodríguez F, Oliver JF, Marín A, Medina, JR. 1990. The general stochastic model of nucleotide substitution. J. Theoret. Biol. 142: 485-501.

Roe BA, Din-Pow M, Wilson RK, Wong JF. 1985. The complete nucleotide sequence of the Xenopus leavis mitochondrial genome. J. Biol. Chem. 260: 9759-9774.

Roelants K, Bossuyt F. 2005. Archaeobatrachian paraphyly and Pangaean diversification of crown-group frogs. Syst. Biol. 54: 111-126.

Rutschmann F. 2004. Bayesian molecular dating using PAML/multidivtime – A step-by-step manual.

Saitou N, Nei M. 1987. The neighbour-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4: 406-425.

Sanchez-Herraiz MJ, Barbadillo LJ, Machordom A, Sanchiz B. 2000. A new species of pelodytid frog from the Iberian Peninsula. Herpetologica 56: 105-118

Sanchiz B. 1998. Encyclopedia of Palaeoherpetology, Part 4. Salientia. München: Pfeil.

Sanchiz B, Tejedo M, Sánchez-Herráiz MJ. 2002. Osteological differentiation among Iberian Pelodytes (Anura, Pelodytidae). Graellsia 58: 35-68.

San Mauro D, Vences M, Alcobendas M, Zardoya R, Meyer A. 2005. Initial diversification of living amphibians predated the breakup of Pangaea. Am. Nat. 165: 590-599.

Schmidt HA, Strimmer K, Vingron M, Haeseler A von. 2002. TREE-PUZZLE, version 5.0. München.

Steel MA, Lockhart PJ, Penny D. 1993. Confidence in evolutionary trees from biological sequence data. Nature 364: 440-442.

Steinfartz S, Veith M, Tautz D. 2000. Mitochondrial sequence analysis of Salamandra taxa suggests old splits of major lineages and postglacial recolonizations of Central Europe from distinct source populations of Salamandra salamandra . Mol. Ecol. 9: 397-410.

Strimmer K, Haeseler A von. 1996. Quartet puzzling: a quartet maximum likelihood method for reconstructing tree topologies. Mol. Biol. Evol. 13: 964-969.

Sullivan J, Swofford DL. 2001. Should we use model-based methods for phylogenetic inference when we know that assumptions about among-site rate variation and nucleotide substitution pattern are violated? Syst. Biol. 50: 723-729.

Swofford DL. 2001. Paup*: Phylogenetic Analysis Using Parsimony (and Other Methods), version 4.06b. Sunderland/Massachusetts: Sinauer Associates.

Tajima F. 1993. Simple methods for testing the molecular evolutionary clock hypothesis. Genetics 135: 599-607.

Terentjev P, Chernov SA. 1965. Key to amphibians and reptiles. Jerusalem: Israel Program for Scientific Translations.

Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG. 1997. The Clustal X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 24: 4876-4882.

Thorne JL, Kishino H. 2002. Divergence time and evolutionary rate estimation with multilocus data. Syst. Biol. 51: 689702.

Thorne JL, Kishino H, Painter IS. 1998. Estimating the rate of evolution of the rate of molecular evolution. Mol. Biol. Evol. 15: 1647-1657.

Titus TA, Larson A. 1995. A molecular phylogenetic perspective on the evolutionary radiation of the salamander family Salamandridae. Syst. Biol. 44: 125-151.

Ugurtas IH, Ljubisavljevic K, Sidorovska V, Kalezic ML, Dzukic G. 2002. Morphological differentiation of eastern spadefoot toad ( Pelobates syriacus) populations. Israel J. Zool. 48: 13-32.

Veith M, Steinfartz S, Zardoya R, Seitz A., Meyer A. 1998. A molecular phylogeny of ‘true’ salamanders (family Salamandridae) and the evolution of terrestriality of reproductive modes. J. Zool. Syst. Evolut. Res. 36: 7-16.

Veith M, Kosuch J, Vences M. 2003. Climatic oscillations triggered post-Messinian speciation of Western Palearctic brown frogs (Amphibia, Ranidae). Mol. Phyl. Evol. 26: 310-327.

Veith M, Schmidtler JF, Kosuch J, Baran I, Seitz A. 2003. Paleoclimatic changes explain Anatolian mountain frog evolution: a test for alternating vicariance and dispersal events. Mol. Ecol. 12: 185-199.

Veith M, Mayer C, Samraoui B, Donaire Barroso D, Bogaerts S. 2004. From Europe to Africa and vice versa: evidence for multiple intercontinental dispersal in ribbed salamanders (genus Pleurodeles). J. Biogeography 31: 159-171.

Wilson RCL, Drury SA, Chapman JL. 2000. The Great Ice Age. London and New York: The Open University.

Yoder AD, Irwin JA, Payseur BA. 2001. Failure of the ILD to determined data combinability for slow loris phylogeny. Syst.Biol. 50: 408-424.

Yoder A, Yang Z. 2004. Divergence data for Malagasy lemurs estimated from multiple gene loci: genealogical and evolutionary context. Syst. Biol. 13: 757-773.

Zhang D-X, Hewitt GM. 1996. Nuclear integrations. Challenges for mitochondrial DNA markers. TREE 11: 247-251.