Open Access

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

  • Jan W. Arntzen1,
  • Tania Trujillo2,
  • Roland Butôt1,
  • Klaas Vrieling3,
  • Onno Schaap3,
  • Jorge Gutiérrez-Rodríguez2 and
  • Iñigo Martínez-Solano1, 4, 5, 6Email author
Frontiers in Zoology201613:52

DOI: 10.1186/s12983-016-0184-7

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.


Cline analysis Hybridization Microsatellites Morphology mtDNA SNPs


Speciation is usually accompanied by the gradual accumulation of genetic incompatibilities between diverging gene pools through time [1]. 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 [2]. 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 [3].

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 [46]. 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, 710]. 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 [5].

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 [13]) 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 [14].

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 [1416]. 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 [17]. These species split in the Late Miocene at around 9.2 MYA [14], 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 [16]. 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 [18]. 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.


We gathered morphological and molecular genetic data for 344 Bufo toads in 17 populations in a transect across the contact zone of B. spinosus and B. bufo in northwestern France (Table 1, Fig. 1) adjacent to a locality where we documented a B. spinosus – B. bufo hybrid population [18]. Morphological characters studied were those that help to distinguish B. spinosus from B. bufo [16], namely Snout-Urostyle length (SUl), anterior and posterior Parotoid distance (Pda, Pdp), length and width of the inner Metatarsal Tubercle (MTl, Mtw) and Paratoid angle (Pa) and the derived measures Paratoid divergence (Pd = Pda/Pdp), Metatarsal Tubercle size (MTsize = MTl/SUl) and shape (MTshape = MTw/MTl). For illustrations of these characters see [16]. Adult males (N = 190) and adult females (N = 76) were analyzed separately. For each individual we calculated the probability to belong to B. spinosus (Ps), with the logistic regression formulae derived for this part of France (see [16]).
Table 1

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)

Sample size



Northern latitude

Eastern longitude

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)




Montmean, Bazoches-sur-Hoëne





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)




le Cottin





18 (17,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 (%)





Fig. 1

Location of 17 Bufo spinosus and B. bufo populations over a transect in northwestern France. The pie diagrams at the bottom show the observed mtDNA haplotype frequencies, with haplotypes typical for the species in red (B. spinosus) and blue (B. bufo). The B. bufoB. spinosus hybrid population reported by [18] is shown by an X

Diagnostic SNPs were developed on 16S rRNA (16S) and Cytochrome B (cytb) mitochondrial sequences from [14], 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 [14] 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 [19] 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 [18]. The number of alleles (Na), the observed (Ho) and expected heterozygosity (He) were calculated with GenePop 4.2 [20] 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 [21]. For the microsatellite data the presence of null alleles was investigated with Micro-Checker 2.2.3 [22].

The software Structure [23] 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 [24] as implemented in Structure Harvester [25]. Results were processed with the pipeline CLUMPAK [26]. Individual genotypes were also assessed with NewHybrids software under default settings [27] 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 [28]).

The position, width and shape of geographical clines over the transect were determined with the R package HZAR [29] 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 [30] was used to summarize the molecular allelic information in Principal Component vectors over all loci. The other statistical procedures were done with SPSS 21 [31].


At the southern side of the transect adult toads show the morphological character states typical for B. spinosus (metatarsus tubercle long and narrow, parotoids divergent) and at the northern side of the transect they show states typical for B. bufo (metatarsus tubercle small and round, parotoids parallel). Among males, morphologically intermediate individuals are most frequent in populations 9, 10, 13 and 15 (metatarsus tubercle, Fig. 2a) or in populations 12 and 15 (parotoids, Fig. 2b).
Fig. 2

Species specific morphometric character states across a Bufo spinosus to B. bufo transect in the northwest of France. Data points shown are average values for males in 13 populations (see Table 1). a size and shape of the metatarsus tubercle and b positioning of the parotoids

For the microsatellite data, the observed number of alleles (Na), expected (He) and observed (Ho) heterozygosities and allele size range per population are summarized in Additional file 1: Appendix I. In general, genetic diversity is higher in the centre of the transect and lower at the outer sides of the transect (Table 2A). Significant instances of linkage disequilibrium (LD) were only detected in the Bspi 3.11 × Bspi 3.19 locus combination in population 6 and in the Bspi 4.28 × Bspi 4.29 combination in population 12. Across populations LD was significant for Bs4.28 * Bs4.29 and at the centre of the transect (populations 7–12) LD was significant for Bspi4.25 * Bspi4.28 and Bpis4.27 * Bspi4.28. Micro-Checker inferred the likely presence of null alleles 12 times in nine populations. Accordingly, because of the likely presence of null alleles in loci Bspi 4.24 and 4.25, analyses with Structure and NewHybrids were repeated with these markers excluded. The results were very similar (not shown) and only results based on the full dataset are reported. Significant deviations from Hardy-Weinberg equilibrium (H1: heterozygote deficit) were detected eight times in seven populations, involving different loci and with no correspondence with potential null alleles (Table 2A). The percentage of variation explained along the principal component axes was 14.3, 4.7 and 4.1% for the first, second and third axis, respectively. Only the first axis was retained, since the others did not discriminate between B. bufo and B. spinosus.
Table 2

Molecular genetic characteristics across 17 populations in a Bufo spinosus to B. bufo transect in the northwest of France for 12 microsatellite loci (A) and four SNP loci (B). Within sections panels represent: upper panel–variability measures, middle panel–deviations from neutral expectations, and lower panel–proportion of populations falling into particular Structure and NewHybrids categories. Green and yellow shadings indicate values falling into the lower and upper quartiles, respectively. For locality information see Fig. 1 and Table 1. For data per locus see Additional file 1: Appendix I. The unprocessed data are in Additional file 3: Appendix III and Additional file 4: Appendix IV

# Microsatellite loci are 1 - Bspi3.11, 2 - Bspi3.19, 3 - Bspi4.16, 4 - Bspi4.24, 5 - Bspi 4.25, 6 - Bspi4.28 and 7 - Bspi4.29

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).

Bivariate Structure * NewHybrid plots show typical unimodal distributions for the microsatellite and the nuclear SNP data (Fig. 3). The bell-shaped curves suggest that the analytical results are strongly correlated between methods. Following the Structure axis, populations 1–17 are in geographical order, with few exceptions. Low hybrid values (pp < 0.2) are shown for populations 1–6 and for populations 15–17 (microsatellites) or 16 and 17 (SNPs). These populations show extreme Structure Q-scores also and represent pure or nearly pure B. spinosus and B. bufo, respectively. High hybridity values (pp > 0.8) are shown by populations 10 and 12, with intermediate values for the remaining populations. The average Qb-score for population 10 (Beaulieu) is 0.66 (standard deviation, SD = 0.150) suggesting that it is situated just off the centre of the contact, in the direction of B. bufo (see also the geographical cline analysis below). Populations 9, 11 and 12 adjacent to Beaulieu contain genetically mixed individuals with pp- and Q-scores as in Beaulieu, as well as individuals with lower pp- and intermediate Q-scores for which the status as hybrid or a parental is equivocal. The microsatellite plot (Fig. 3a) has ‘high shoulders’ (i.e. it is platykurtic), where high NewHybrids scores are coupled with extreme Structure scores. Thirty-eight out of 40 (95%) of the individuals from Beaulieu have NewHybrids scores pp > 0.975 and Structure scores 0.1 < Q < 0.9, suggesting that they have an intermediate genetic composition.
Fig. 3

Genotype profile over 12 microsatellite markers (a) and four nDNA SNPs (b) in Bufo spinosus and B. bufo over a transect in northwestern France. The classifications are with Structure (Q-score, horizontal axis) and with NewHybrids (pp, posterior probability of falling in the pooled hybrid class, vertical axis). Data points represent individuals (solid round symbols) or population means (open round symbols). For locality information see Table 1

A sharp cline is fitted to the Structure Q-scores for the microsatellites, accompanied by wider and less pronounced clines for the NewHybrid transitions from B. spinosus to the hybrid class (Fig. 4a, left hand cline) and from the hybrid class to B. bufo (Fig. 4a, right hand cline). In the central part of the transect, we tentatively link the position of the clines to the altitudes of the populations, that are high (≥285 m a.s.l.) for populations 5 and 6 at the B. spinosus side, low (≤175 m a.s.l.) for populations 13–15 at the B. bufo side of the transect and intermediate for the central populations 7–12 (201–277 m a.s.l.). For five out of eight characters the position of the cline is similar to that of the reference and for three out of eight also the width is similar (Table 3). Significantly displaced and significantly wider clines than the reference are observed for the two morphological characters and for the nuclear POMC locus. In all three cases the displacement is towards the B. bufo side of the transect (Fig. 4, Table 3). Clines wider than the reference are also found for the microsatellite information expressed in the first PCA and for the locus BDNF.
Fig. 4

Geographical cline analysis for Bufo spinosus (left, with Qs-values close to unity at vertical axis) and B. bufo (right, with Qs-values close to zero at vertical axis) in a transect in northwestern France. Study populations are as in Table 1 and small solid dots are population averages. Clines were fitted to a) Structure Q-scores (with blue shading) and NewHybrids posterior probabilities for the species (grey shadings with the transition towards B. spinosus left and towards B. bufo right). Note that the graph to the right shows the central part of the transect in detail. Furthermore, b) the loading on the first PCA axis, based on the panel of 12 microsatellites, c) the B. spinosus frequency of the mtDNA marker, d, e) morphological identification probabilities (details see text) and f–i) the B. spinosus frequency of four SNP markers. The 95% credible cline regions are highlighted by grey shadings. The vertical blue line gives the position of population 10, at Beaulieu. Locality altitudes are plotted in red

Table 3

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



Character, unit




δ L

δ R

τ L

τ R

p min

p max

Microsatellite reference clines

 Structure, Q


−1.0 (−1.0–0.7)

30.3 (21.5–44.4)




 NewHybrids, posterior probability

   - B. spinosus side of transect


−1.0 (−1.0–−0.2)

0.6 (0.0–1.3)







   - B. bufo side of transect


10.8 (8.4–20.0)

12.2 (6.0–29.9)







 Thresholds–min, max


−1.0, 20.0

0.0, 44.4


Tested clines

 PCA on microsatellites, loading first axis


4.6 (−1.0−15.7)

104.9 (64.4–179.7)




 mtDNA, haplotype frequency


−1.0 (−1.0–2.4)

43.5 (31.0–64.0)




 Parotoids, classification probability


65.5 (27.3–122.8)

255.1 (107.7–> 400)




 Metatarsal tubercle, idem


96.6 (54.2–145.3)

298.7 (164.4–> 400)




 BDNF, allele frequency


3.8 (−0.9–18.6)

58.9 (46.9–115.8)







 POMC, idem


22.3 (9.2–41.5)

90.7 (52.6–176.7)




 RAG1, idem


4.0 (−0.7–9.1)

35.7 (23.4–56.2)







 RPL3, idem


−0.9 (−1.0–3.9)

26.8 (23.0–49.1)







For the chosen models p min and p max are fixed to 0 and 1 (typ models), or their empirical values (fix models), or p min and p max are fitted (opt model). Tail fitting encompassed right (R), left (L), none (N) or both fitted (B). Cline position is the centre of the contact zone and ca. equals the distance to Beaulieu (in km, population 10, see Table 1). Cline width (w) is calculated as 1/maximum slope. Two log-likelihood unit support limits are presented in parentheses. δ and τ are the shape parameters for the left and right tails, and p min and p max are the character states at either end of the transect


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 [34] 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. [28]). 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 [18], 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, 3639] 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 [35]. 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 [19]. 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 [40]. 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 [41], together with direct observations [42] 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 dt [43]. 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 [46]. 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.

Authors’ contributions

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.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, 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 ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

Naturalis Biodiversity Center
Museo Nacional de Ciencias Naturales, CSIC
Plant Cluster, Institute of Biology, Sylvius Laboratory, Leiden University
Instituto de Investigación en Recursos Cinegéticos (IREC-UCLM-CSIC-JCCM)
CIBIO-InBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, Campus Agrário de Vairão, Universidade do Porto
Present address: Ecology, Evolution and Development Group, Department of Wetland Ecology, Doñana Biological Station, CSIC


  1. 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
  2. Harrison RG, Larson EL. Hybridization, introgression, and the nature of species boundaries. J Hered. 2014;105:795–809.View ArticlePubMedGoogle Scholar
  3. Hewitt GM. Hybrid zones: natural laboratories for evolutionary studies. Trends Ecol Evol. 1988;3:158–67.View ArticlePubMedGoogle Scholar
  4. Endler JA. Geographic variation, speciation and clines. New Jersey: Princeton University Press; 1977.Google Scholar
  5. Barton NH, Hewitt GM. Analysis of hybrid zones. Annu Rev Ecol Evol Syst. 1985;16:113–48.View ArticleGoogle Scholar
  6. 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
  7. Barton NH, Hewitt GM. Adaptation, speciation and hybrid zones. Nature. 1989;341:497–503.View ArticlePubMedGoogle Scholar
  8. 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
  9. Harrison RG. Hybrid Zones and the Evolutionary Process. Oxford: Oxford University Press; 1993.Google Scholar
  10. 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
  11. 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
  12. Hewitt GM. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.View ArticlePubMedGoogle Scholar
  13. 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
  14. 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
  15. 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
  16. 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
  17. 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
  18. 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
  19. 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
  20. 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
  21. Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.View ArticleGoogle Scholar
  22. 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
  23. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.PubMedPubMed CentralGoogle Scholar
  24. 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
  25. 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
  26. 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
  27. Anderson EC, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002;160:1217–29.PubMedPubMed CentralGoogle Scholar
  28. 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
  29. 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
  30. 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
  31. SPSS. Statistical Package for the Social Sciences. Chicago: SPSS Inc.; 2013.Google Scholar
  32. 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
  33. 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
  34. 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
  35. 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
  36. AmphibiaWeb. University of California, Berkeley, CA, USA. Electronic database accessible at Accessed 13 Oct 2016.
  37. 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
  38. 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 Accessed 13 Oct 2016.
  39. Glandt D. Die Amphibien und Reptilien Europas. Wiebelsheim: Quelle & Meyer; 2015.Google Scholar
  40. Johnston DE. The Channel Islands. An Archaeological Guide. London and Chichester: Phillimore & Co.; 1981.Google Scholar
  41. Reading CJ, Loman J, Madsen T. Breeding pond fidelity in the common toad, Bufo bufo. J Zool. 1991;225:201–11.View ArticleGoogle Scholar
  42. 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
  43. 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
  44. 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
  45. 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
  46. 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


© The Author(s). 2016


By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.