Concordant morphological and molecular clines in a contact zone of the Common and Spined toad (Bufo bufo and B. spinosus) in the northwest of France
© The Author(s). 2016
Received: 14 July 2016
Accepted: 28 November 2016
Published: 19 December 2016
Hybrid zones are regions where individuals of two species meet and produce hybrid progeny, and are often regarded as natural laboratories to understand the process of species formation. Two microevolutionary processes can take place in hybrid zones, with opposing effects on population differentiation. Hybridization tends to produce genetic homogenization, reducing species differences, whereas the presence of mechanisms of reproductive isolation result in barriers to gene flow, maintaining or increasing differences between taxa.
Here we study a contact zone between two hybridizing toad species, Bufo bufo and B. spinosus, through a combination of molecular (12 polymorphic microsatellites, four nuclear and two mitochondrial SNP markers) and morphological data in a transect in the northwest of France. The results show largely concordant clines across markers, defining a narrow hybrid zone of ca. 30 km wide. Most hybrids in the centre of the contact zone are classified as F2 or backcrossed individuals, with no individuals assigned to the F1 hybrid class.
We discuss the implications of these results for our understanding of the evolutionary history of these species. We anticipate that the toad contact zone here described will become an important asset in the study of hybrid zone dynamics and evolutionary biology because of its easy access and the abundance of the species involved.
Speciation is usually accompanied by the gradual accumulation of genetic incompatibilities between diverging gene pools through time . Thus, there are temporal windows during the speciation process where reproductive isolation is sufficient for species to remain differentiated while they are still capable of producing hybrids. This process of natural hybridization allows the incorporation of alleles of one species into the genetic background of another (introgression). Often, alleles introgressing most easily across species boundaries are those associated with an adaptive advantage, whereas those subject to divergent selection or related to reproductive isolation will show little or no introgression. Studying patterns of differential introgression across markers in hybrid zones thus provides a means to identify genes associated to local adaptation and reproductive isolation, which is key to understand how species boundaries persist in the face of gene flow . When barriers to hybridization are incomplete, as in many secondary contact zones, a wide range of possible outcomes can be expected from an evolutionary perspective, from genetic homogenization to the maintenance of independent evolutionary trajectories following the evolution of full reproductive isolation. But it is the grey area in between these extremes, where different genomic regions are shown to introgress at varying levels, that provides the most information about the speciation process. This is the most common situation in nature and the reason why hybrid zones are regarded as natural laboratories for the study of speciation .
A powerful approach to characterize the outcomes of post-divergence hybridization and understand the underlying evolutionary processes is provided by cline theory. Geographically, hybrid zones can be defined as clines in characters that show different frequencies in the parental populations. Therefore, it is possible to study hybrid zones by sampling along a transect perpendicular to the contact zone and examining patterns of variation across multiple independent markers [4–6]. Several theoretical models have been proposed to explain the structure and temporal maintenance of hybrid zones, differing in the relative roles of extrinsic (environmental) and intrinsic (organismal) factors in maintaining some degree of reproductive isolation between hybridizing species. One of the best known examples are tension zones, where selection is environment independent and phenotypic and genetic characters usually show concordant, steep clines that are often narrow relative to the geographic distributions of the parental populations and the dispersal capacity of the organisms under study. Tension zones tend to remain stable through time because of incoming genetic flow from the parental populations at both ends of the contact zone and the elimination of hybrid individuals at the centre by natural selection [5, 7–10]. However, when tension zones form in areas of low population density, their position can shift through time if environmental conditions trigger concomitant changes in density .
Since the publication of some, now classical synthesis papers on the role of Quaternary Ice Ages in shaping current patterns of genetic diversity and structure across taxa [11, 12], there has been much progress in documenting hybrid zones in Europe from a comparative, cross-taxon perspective. This has revealed patterns of concordance in the geographic location of contact zones, which has often been interpreted in terms of a parallel evolutionary history (linked to inferred routes of post-glacial dispersal) on a shared geographic background, especially in taxa with similar dispersal abilities. One of these well-known “suture zones” runs from NW to SE France, with a pair of pond-breeding amphibians (Triturus cristatus/T. marmoratus ) as a primary example. More recently, a study on the evolutionary history of another deeply divergent amphibian species pair (the Common toad Bufo bufo and the Spined toad B. spinosus) has revealed significant geographic concordance with the Triturus hybrid zone .
Recent molecular studies on the systematics of the Bufo bufo species group have resulted in the recognition of four species, two of which were formerly regarded as subspecies of B. bufo [14–16]. One of these, B. spinosus, present in Iberia and North Africa, contacts B. bufo in France, forming an extensive contact zone along a NW-SE line, approximately from Caen to Lyon . These species split in the Late Miocene at around 9.2 MYA , but are morphologically similar. They can, however, be differentiated in their contact zone and recent studies have shown overall agreement between morphological characters, mitochondrial and nuclear DNA markers in delineating species boundaries in France . Occasional discordance is mostly restricted to slowly evolving nuclear markers and interpreted as resulting from incomplete lineage sorting, although hybridization-mediated gene flow cannot be completely ruled out. In fact, we recently isolated and optimized a suite of 12 microsatellite loci cross-amplifying in both species and, in a preliminary assessment of patterns of gene flow across species, identified a hybrid population in a section of the contact zone in northwestern France . We therefore decided to perform a more detailed analysis on a fine-scale transect across the contact zone using morphological and molecular traits (12 polymorphic microsatellites, four nuclear and two mitochondrial SNP markers) in a cline analysis framework. We here describe results on a transect across the contact zone between B. bufo and B. spinosus in northwestern France. The specific questions that we sought to answer are: i) can clines be observed for morphological, nuclear and mitochondrial markers?, ii) do clines of the different traits coincide? and iii) can the position of the contact zone be linked to features of the landscape? On a broader perspective, our goals include testing: iv) whether this contact zone is a hybrid zone, and v) whether the zone is maintained by selection, with a significant degree of reproductive isolation.
Location, position in the transect and sample size for 17 Bufo populations included in the study. Toads studied are adults except for populations 1–3 and 5 where tadpoles were sampled
Altitude (m.a.s.l.) from Google Earth
Straight line distance from Beaulieu (km)
Morphology (males, females)
Le Poitrineau, Saumur
Étang du Perré
Les Tesnières, Durtal
Carrefour du Poteau, Jublains
27 (13, 14)
Les Fontaines, Pré-en-Pail
Foret de Multonne, Mont des Avaloirs
20 (18, 2)
18 (15, 3)
La Couvendière, Mortagne au Perche
22 (20, 2)
La Rosière, Foret du Perche et de la Trappe
20 (18, 2)
52 (26, 26)
Chateau des Bois Francs
2 (1, 1)
Les Quatre Vouges
15 (12, 3)
2 (2, 0)
Mare du Bois, La Houssaye
23 (22, 1)
Les Puits du Sarrasin, Erloy
25 (15, 10)
Les Ajoncs, Audresselles
22 (11, 11)
266 (190, 76)
Missing data (%)
Diagnostic SNPs were developed on 16S rRNA (16S) and Cytochrome B (cytb) mitochondrial sequences from , yielding cytb_946 allele specific primers GGGGAGTGACTAGWGGATTTGC(G/A) and common primer ACCTCCTGGGCGACCCAGATAA and 16S_307 allele specific primers GGGTGACCACGGAGCACA(G/A) and common primer GCTTAGAAATTTTCTTTCAGCATGGAGGTT for KASP genotyping of haplotypes. For nuclear markers sequences of  were used to derive the following species specific SNP markers: brain-derived neurotrophic factor (BDNF) with the allele specific primers: ACTTTCTTTTTGCTATCCATGGTGAG(T/A) and common primer: TACTGGAACTCGCAATGCCGAACTA, for proopiomelanocortin (POMC) the allele specific primers: GTCATCGGAAATCTACCCAACAGA(T/G) and common primer: TCAAACTCTACAGACAACTCTCTTCTGTA, and for 60S ribosomal protein L3 (RPL3) the allele specific primers TATTCCATGACCTCATCACAGAACA(A/G) and common primer: GCAGCCMATGACTTGTGC. For the recombination activating gene 1 (RAG1) the marker sequences were from  with the allele specific primers: TAATCACCACAATAAATGGCCCAGA(A/G) and common primer: GGTTTAAAGTCCAAAGACCTAGAGGAATA.
KASP PCRs for both nuclear and mitochondrial SNPs were carried out in 1536 wells plates with circa 1 ng of template DNA and 1 μl of reaction mix containing KASP v4.0 reaction mix and primers following the manufacturer’s directions. The conditions of the KASP PCR were as follows: 15 min denaturation at 94 °C followed by 10 cycles of 20s at 94 °C denaturation and 60s at 61 °C of annealing and extension and lowering the annealing and extension temperature by 0.6 °C each cycle, followed by 26 cycles of 20s at 94 °C denaturation and 60s at 55 °C of annealing and extension in a hydrocycler. After 36 cycles of this touch down protocol fluorescence was measured on a Pherastar plate reader. PCR was continued and after 39, 42 and 45 cycles fluorescence was measured again to follow the trajectory of all samples. Genotypes were called automatically by the module Kraken™ of LIMS controlling the LGC genomics SNP genotyping line, visually inspected and if necessary manually corrected.
Twelve microsatellite loci were studied in 300 individuals in four multiplex reactions, following . The number of alleles (Na), the observed (Ho) and expected heterozygosity (He) were calculated with GenePop 4.2  per locus and population. The same program was used to test for deviations from Hardy-Weinberg equilibrium and for linkage disequilibrium, under a sequential Bonferroni correction for multiple tests . For the microsatellite data the presence of null alleles was investigated with Micro-Checker 2.2.3 .
The software Structure  was used to describe population structure and identify genetically admixed individuals for SNP and microsatellite data separately. Runs consisted of 500,000 (burn-in) and 1,000,000 generations under the admixture model with correlated allele frequencies, for K values of 1–14, with five replicates per value of K. The optimal clustering scheme was selected through application of Evanno’s criterion  as implemented in Structure Harvester . Results were processed with the pipeline CLUMPAK . Individual genotypes were also assessed with NewHybrids software under default settings  to obtain the posterior probability that they fall into one of the following categories: B. spinosus, B. bufo, F1-hybrids, F2-hybrids and backcrosses in either direction. Structure and NewHybrids provide different but complementary information about genetic admixture and for visual inspection we constructed bivariate Structure * NewHybrid plots (see ).
The position, width and shape of geographical clines over the transect were determined with the R package HZAR  from population data for i) morphological characters (Ps, Parotoids and Metatarsal tubercle, males and females combined), ii) mitochondrial DNA (Fs, frequency of ‘spinosus’-haplotypes), iii) the loading on the axes of a principal component analysis derived from the panel of microsatellites (PCA1, etc.) and iv) for the nuclear genes BDNF, POMC, RAG1 and RPL3 (Fs, the frequency of B. spinosus specific SNPs). As reference cline models we choose the Structure Q-score of belonging to either parental class, and the NewHybrid posterior probabilities of belonging to a hybrid class versus either of the parental species, derived from the microsatellite data. Clines are considered significantly displaced if the two log-likelihood unit support limits of the cline centre do not overlap with the Structure Q-score (Qb = 1 − Qs). Individual microsatellite loci conveyed little information on species identity (results not shown) and GenoDive  was used to summarize the molecular allelic information in Principal Component vectors over all loci. The other statistical procedures were done with SPSS 21 .
The number of inferred gene pools (K) in Structure was K = 2 (see Additional file 2: Appendix II), with low, intermediate and high Qb-scores at the southeastern, central and northwestern parts of the transect, respectively (Table 2A). In NewHybrids high values for the posterior probability of belonging to a hybrid class were observed in the central part and low scores at either side of the transect. The probability of being classified as an F2-hybrid was several orders of magnitude larger than for the other hybrid classes discerned, suggesting a preponderance of hybrids at deeper levels of introgressive hybridization. The results for the SNP data show the same overall pattern as that for the microsatellites, with no significant results for HW- and LD-tests (Table 2B).
Parameter estimates for the geographical clines shown in Fig. 4. The results for Structure and NewHybrids on the panel of microsatellites are taken as the reference. Significant differences with the 2log-likelihood unit support limits falling outside the range observed for the reference clines are shown in boldface type
Microsatellite reference clines
NewHybrids, posterior probability
- B. spinosus side of transect
- B. bufo side of transect
PCA on microsatellites, loading first axis
mtDNA, haplotype frequency
Parotoids, classification probability
255.1 (107.7–> 400)
Metatarsal tubercle, idem
298.7 (164.4–> 400)
BDNF, allele frequency
Secondary hybrid zones form when taxa that have evolved in allopatry meet and exchange genes. The outcome of the species interaction depends on the evolution of mechanisms of reproductive isolation, both pre- and post-zygotic, in which isolation usually increases with time since divergence, with examples all along the continuum between genetic merger (as in ‘hybrid species’) and full reproductive isolation. For amphibian examples in which the extent of natural hybridization in contact zones correlates with divergence times see refs. [28, 32, 33].
In the hybrid zone studied here, most hybrid individuals were classified by NewHybrids as F2s or backcrosses. While this could represent differences in fertility between hybrid classes, it is also likely that later generation hybrids and backcrosses are being pooled together in the F2 class, since the confidence of assignment to different hybrid classes depends on the number of generations elapsed after secondary contact  and the contact is probably old (see below). The lack of full reproductive isolation may seem surprising given that our study system is formed by two non-sister toad species that diverged in the Miocene [14, 35]. Thus, some degree of reproductive isolation due to the evolution of genetic incompatibilities was expected, but long lasting genetic compatibility is not uncommon in amphibians (e.g. ). Initial assessments of patterns of allele sharing across species in this system were inconclusive as to whether gene flow or incomplete lineage sorting best explained observed discordance between mtDNA and nDNA data [14, 16, 17, 19]. This, and the characterization of a hybrid population based on morphological and molecular genetic data , raised questions about the extent of reproductive isolation in this species pair. However, our comprehensive transect sampling has revealed sharp, spatially concordant and narrow genetic clines (Fig. 4), suggestive of barriers to hybridization. Tentatively, this hybrid zone can be located along the ‘Collines de Normandie’, with B. spinosus at higher altitudes than B. bufo, and a centre close to locality 10 (Beaulieu) in our transect. Nevertheless, it is unclear the extent to which selection falls on the parentals and/or the hybrids and whether environmentally mediated selection plays a role in stabilizing the hybrid zone. The taxonomic status of B. spinosus as a species different from B. bufo is recent [14, 16, 18, 19, 36–39] and not entirely without debate [17, 35]. The results showing limited gene flow across a narrow hybrid zone (despite no barriers to dispersal) bring evidence that reproductive isolation is probably involved and further argues for species status of the Spined toad, B. spinosus.
In general, allopatric populations of the two parental species have lower genetic diversity than populations near the centre of the contact zone. The inspection of allele frequencies shows a clear trend in some loci, with some alleles presenting high frequencies on opposite sides of the contact zone and intermediate frequencies co-occurring in the centre of the contact zone. Genetic homogenization due to hybridization is spatially restricted. This is also apparent upon inspection of Structure and NewHybrids analyses, which show increasing levels of admixture close to the centre of the contact zone. Thus, while we have found evidence for ongoing hybridization between the two Bufo species, evidence for hybrids is restricted to a narrow area about 30 km wide. Most hybrids were classified as F2 individuals by NewHybrids, which is often the case in relatively old contact zones where many generations of backcrossing complicate precise assignment of individuals to specific hybrid classes [8, 18, 33, 34]. While the original position of this hybrid zone is unknown, the default scenario is that it formed at or close to its current location within the last glacial cycle.
Phylogeographic analyses on the basis of two different climate reconstructions for the Last Glacial Maximum (LGM), MIROC and CCSM, lead to drastically different scenarios . Under MIROC, climate conditions are favourable over most of Europe south of 50°N, therewith suggesting that the B. bufo - B. spinosus contact pre-dates the LGM. Conversely, under CCSM favourable areas are situated at more southern latitudes, with small fragments in coastal areas (both at the Atlantic and the Mediterranean) and large continuous areas restricted to southern Europe, suggesting that the Bufo hybrid zone formed more recently. We propose a scenario in which B. spinosus reached the Normandy region around 8000 YBP, because the species is absent from the Channel island of Guernsey that became isolated from the continent shortly before that date, but made it to Jersey . Jersey was the last of the Channel Islands to be separated from the European mainland, with the final inundation of the English Channel in roughly that period . In parallel and over the same period B. bufo, with a glacial refugium in the Balkan peninsula, reached the landmass that is currently England, Scotland and Wales, but it didn’t reach the Isle of Man or Ireland. We therefore assert that the species contact is thousands of years old, but further studies are required to reconstruct the species’ postglacial dynamics and to track the spatial and temporal stability of the hybrid zone.
The concordance of mtDNA and nuclear clines is expected when there is selection against F1 hybrids, but sex-biased dispersal or asymmetric genetic incompatibilities among the sexes can produce discordant clines. While data for B. spinosus are scarce, a study on British and Swedish populations of B. bufo found that 93% of females and 96% of males that survived between years returned to the same breeding ponds , together with direct observations  suggesting the absence of sex biased dispersal. Furthermore, a cline as narrow as observed cannot result from neutral processes exclusively. In the absence of selection, the width (w) of the cline can be predicted from a diffusion model as a function of dispersal distance (d) and the length of time since contact (t), as w = 2.51 d √t . The maximum dispersal distance recorded for B. bufo/B. spinosus averages at 1.32 km in eight studies (range 0.12–3.62 km) [42, 44, 45] Generation time is 2 years for male and 3 years for female Bufo . At a low average dispersal of 1 km per generation, cline widths would exceed the 30.3 km width of the Structure Q reference cline (Table 3) in ca. 400 years and at a high average dispersal of 5 km per generation this would last only 17 years. Arguably, the hybrid zone is much older than this and cline width appears to have been kept in check over several thousands of years, presumably through selection against hybrids, especially F1’s, and/or the parentals, given their absence and low frequency at the centre of the zone, respectively.
Despite the overall congruence of genetic clines, with similar patterns between nuclear and mitochondrial markers, some clines show smoother transitions between pure populations at both ends of the transect. One remarkable discrepancy is between the sharp species transition documented by the Structure and NewHybrid analyses of the microsatellite data (Fig. 4A) versus the smooth transition revealed by the first PCA-axis obtained from same data (Fig. 4B). Given the general evidence for a narrow, ca. 30 km wide hybrid zone, this result suggests that PCA does not perform well in extracting the species-specific character states. Also the morphological clines are smooth (Fig. 4DE), which can be due to a variety of factors, all of which are worth further exploration. Possibly, some species diagnostic morphological character states are yet to be discovered. On the other hand, morphological features distinguishing the two toad species are probably the result of additive effects from many genes and complex interactions with the environment, so perhaps a looser connection than with SNP or microsatellite loci should be expected. Finally, sampling variance or micro-environmental selection could also have a role in smoothing the morphological cline. Whatever the case, the presumably extensive contact zone offers ample opportunities to test the consistency of previously identified morphological and genetic characters in diagnosing species across different sections of their contact zone.
We have characterized the common toad system as a hybrid zone, with narrow and coincident clines revealing barriers to gene flow across species, which opens new and exciting lines of research to further pursue with the help of genomic tools. While the relative roles of intrinsic and extrinsic factors in shaping genetic and morphological clines are still open to further scrutiny, the geographical extent of this contact zone, encompassing France from the northwest to the southeast, allows the comparative study of multiple transects, with great opportunity for the identification of consistently versus locally-introgressed genomic regions and thus for discerning the relative roles of indeterminate drift and selection in the origin and maintenance of species boundaries. Thus, we assert that the hybrid zone between B. bufo and B. spinosus will become a prime study system in speciation research.
We thank M. Dijkstra, P. Steenbergen and J. Ziermann for help in the laboratory and A. Sánchez Vialas and A. Zuiderwijk for help with field work. We thank two anonymous reviewers for constructive comments on an earlier version of the manuscript.
The research was supported by grants CGL2008-04271-C02-01/BOS and CGL2011-28300 from the ‘Ministerio de Ciencia e Innovación’, the ‘Ministerio de Economía y Competitividad’ and FEDER and by grant PPII10-0097-4200 from the ‘Junta de Comunidades de Castilla la Mancha, to IMS. JGR was supported with a JAE pre-PhD fellowship by the ‘Consejo Superior de Investigaciones Científicas of Spain’ and the European Social Fund. IMS was supported by the project ‘Biodiversity, Ecology and Global Change’, co-financed by the North Portugal Regional Operational Programme 2007/2013 (ON.2–O Novo Norte), under the National Strategic Reference Framework, through the European Regional Development Fund, by a Naturalis ‘Temminck-fellowship’ and currently receives funding from the Spanish Severo Ochoa Program (SEV-2012-0262).
Availability of data and materials
The molecular data generated and analyzed during this study are included in this published article as supplementary information files; the morphometric data are available from the corresponding author upon request.
Study design was by JWA and IMS, sampling was by JWA, RB and IMS, analyses were by JWA, TT, JGR and IMS, laboratory work was by TT and JGR for the microsatellite data and by RB, KV and OS for the SNP data. Writing and manuscript production were by JWA and IMS with input by TT and KV. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Open AccessThis 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.
- Seehausen O, Butlin RK, Keller I, Wagner CE, Boughman JW, Hohenlohe PA, Peichel CL, Saetre GP, Bank C, Brännström Å, Brelsford A. Genomics and the origin of species. Nat Rev Genet. 2014;15:176–92.View ArticlePubMedGoogle Scholar
- Harrison RG, Larson EL. Hybridization, introgression, and the nature of species boundaries. J Hered. 2014;105:795–809.View ArticlePubMedGoogle Scholar
- Hewitt GM. Hybrid zones: natural laboratories for evolutionary studies. Trends Ecol Evol. 1988;3:158–67.View ArticlePubMedGoogle Scholar
- Endler JA. Geographic variation, speciation and clines. New Jersey: Princeton University Press; 1977.Google Scholar
- Barton NH, Hewitt GM. Analysis of hybrid zones. Annu Rev Ecol Evol Syst. 1985;16:113–48.View ArticleGoogle Scholar
- Hewitt GM. The subdivision of species by hybrid zones. In: Otte D, Endler J, editors. Speciation and its consequences. Massachusetts: Sinauer Associates, Sunderland; 1989.Google Scholar
- Barton NH, Hewitt GM. Adaptation, speciation and hybrid zones. Nature. 1989;341:497–503.View ArticlePubMedGoogle Scholar
- Szymura JM, Barton NH. The genetic structure of the hybrid zone between the fire-bellied toads Bombina bombina and B. variegata: comparisons between transects and between loci. Evolution. 1991;45:237–65.View ArticleGoogle Scholar
- Harrison RG. Hybrid Zones and the Evolutionary Process. Oxford: Oxford University Press; 1993.Google Scholar
- Rice AM, Rudh A, Ellegren H, Qvarnström A. A guide to the genomics of ecological speciation in natural animal populations. Ecol Lett. 2010;14:9–18.View ArticlePubMedGoogle Scholar
- Hewitt GM. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996;58:247–76.View ArticleGoogle Scholar
- Hewitt GM. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.View ArticlePubMedGoogle Scholar
- Arntzen JW, Wallis GP. Restricted gene flow in a moving hybrid zone of the newts Triturus cristatus and T. marmoratus in western France. Evolution. 1991;45:805–26.View ArticleGoogle Scholar
- Recuero E, Canestrelli D, Vörös J, Szabó K, Poyarkov NA, Arntzen JW, Crnobrnja-Isailovic J, Kidov AA, Cogălniceanu D, Caputo FP, Nascetti G, Martínez-Solano I. Multilocus species tree analyses resolve the radiation of the widespread Bufo bufo species group (Anura, Bufonidae). Mol Phylogenet Evol. 2012;62:71–86.View ArticlePubMedGoogle Scholar
- Litvinchuk S, Borkin L, Skorinov D, Rosanov J. A new species of common toads from the Talysh mountains, south-eastern Caucasus: genome size, allozyme, and morphological evidences. Russ J Herpetol. 2008;15:19–43.Google Scholar
- Arntzen JW, McAtear J, Recuero E, Ziermann JM, Ohler A, van Alphen J, Martínez-Solano I. Morphological and genetic differentiation of Bufo toads: two cryptic species in Western Europe (Anura, Bufonidae). Contrib Zool. 2013;82:147–69.Google Scholar
- Arntzen JW, Recuero E, Canestrelli D, Martínez-Solano I. How complex is the Bufo bufo species group? Mol Phylogenet Evol. 2013;69:1203–8.View ArticlePubMedGoogle Scholar
- Trujillo T, Gutiérrez-Rodríguez J, Arntzen JW, Martínez-Solano I. Morphological and molecular data to describe a hybrid population of the Common toad (Bufo bufo) and the Spined toad (Bufo spinosus) in western France. Contrib Zool. 2017;86:1.Google Scholar
- Arntzen JW, Wilkinson JW, Butot R, Martínez-Solano I. A new vertebrate species native to the British Isles: Bufo spinosus (Daudin, 1803) in Jersey. Herpetol J. 2014;24:209–16.Google Scholar
- Rousset F. Genepop’007: a complete reimplementation of the Genepop software for Windows and Linux. Mol Ecol Resour. 2008;8:103–6.View ArticlePubMedGoogle Scholar
- Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.View ArticleGoogle Scholar
- van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.View ArticleGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.PubMedPubMed CentralGoogle Scholar
- Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.View ArticlePubMedGoogle Scholar
- Earl DA, VonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.View ArticleGoogle Scholar
- Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. CLUMPAK: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015;15:1179–91.View ArticlePubMedPubMed CentralGoogle Scholar
- Anderson EC, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002;160:1217–29.PubMedPubMed CentralGoogle Scholar
- Arntzen JW, Wielstra B, Wallis GP. The modality of nine Triturus newt hybrid zones assessed with nuclear, mitochondrial and morphological data. Biol J Linn Soc. 2014;113:604–22.View ArticleGoogle Scholar
- Derryberry EP, Derryberry GE, Maley JM, Brumfield RT. HZAR: hybrid zone analysis using an R software package. Mol Ecol Resour. 2014;14:652–63.View ArticlePubMedGoogle Scholar
- Meirmans PG, van Tienderen PH. GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Mol Ecol Notes. 2004;4:792–4.View ArticleGoogle Scholar
- SPSS. Statistical Package for the Social Sciences. Chicago: SPSS Inc.; 2013.Google Scholar
- Hickerson MJ, Meyer CP, Moritz C. DNA barcoding will often fail to discover new animal species over broad parameter space. Syst Biol. 2006;55:729–39.View ArticlePubMedGoogle Scholar
- Dufresnes C, Brelsford A, Crnobrnja-Isailović J, Tzankov N, Lymberakis P, Perrin N. Timeframe of speciation inferred from secondary contact zones in the European tree frog radiation (Hyla arborea group). BMC Evol Biol. 2015;15:155.View ArticlePubMedPubMed CentralGoogle Scholar
- Pereira RJ, Wake DB. Genetic leakage after adaptive and nonadaptive divergence in the Ensatina eschscholtzii ring species. Evolution. 2009;63:2288–301.View ArticlePubMedGoogle Scholar
- Garcia-Porta J, Litvinchuk SN, Crochet P-A, Romano A, Geniez PH, Lo-Valvo M, Lymberakis P, Carranza S. Molecular phylogenetics and historical biogeography of the west-palearctic common toads (Bufo bufo species complex). Mol Phylogenet Evol. 2012;63:113–30.View ArticlePubMedGoogle Scholar
- AmphibiaWeb. University of California, Berkeley, CA, USA. Electronic database accessible at http://amphibiaweb.org. Accessed 13 Oct 2016.
- Carretero MA, Martínez-Solano I, Ayllón E, Llorente G. Lista patrón de los anfibios y reptiles de España (Actualizada a diciembre de 2014). Barcelona: Asociación Herpetológica Española; 2011.Google Scholar
- Frost D. Amphibian Species of the World: an Online Reference. Version 6.0. American Museum of Natural History, New York, USA. Electronic database accessible at http://research.amnh.org/herpetology/amphibia/index.html. Accessed 13 Oct 2016.
- Glandt D. Die Amphibien und Reptilien Europas. Wiebelsheim: Quelle & Meyer; 2015.Google Scholar
- Johnston DE. The Channel Islands. An Archaeological Guide. London and Chichester: Phillimore & Co.; 1981.Google Scholar
- Reading CJ, Loman J, Madsen T. Breeding pond fidelity in the common toad, Bufo bufo. J Zool. 1991;225:201–11.View ArticleGoogle Scholar
- Smith MA, Green DM. Dispersal and the metapopulation paradigm in amphibian ecology and conservation: are all amphibian populations metapopulations? Ecography. 2005;28:110–28.View ArticleGoogle Scholar
- Barton NH, Gale KS. Genetic Analysis of Hybrid Zones. Pp. 13–45 in Harrison RG. Hybrid Zones and the Evolutionary Process. Oxford: Oxford University Press; 1993.Google Scholar
- Daversa DR, Muths E, Bosch J. Terrestrial movement patterns of the Common Toad (Bufo bufo) in Central Spain reveal habitat of conservation importance. J Herpetol. 2012;46:658–64.View ArticleGoogle Scholar
- Trochet A, Moulherat S, Calvez O, Stevens VM, Clobert J, Schmeller DS. A database of life-history traits of European amphibians. Biodivers Data J. 2014;2:e4123.View ArticleGoogle Scholar
- Hemelaar A. Age, growth and other population characteristics of Bufo bufo from different latitudes and altitudes. J Herpetol. 1988;22:369–88.View ArticleGoogle Scholar