Contributions to Zoology, 86 (1) – 2017
Morphological and molecular data to describe a hybrid population of the Common toad (Bufo bufo) and the Spined toad (Bufo spinosus) in western France
Tania Trujillo1, Jorge Gutiérrez-Rodríguez2, Jan W. Arntzen3, Iñigo Martínez-Solano4
Keywords: Amphibians,cross-amplification,genetic polymorphism,hybridization,microsatellites,morphology
The use of hyper-variable markers across species is often hindered by low cross-species amplification success, a reduced level of polymorphism or a high frequency of null alleles. However, optimizing sets of reliable and informative markers that can be consistently amplified and scored across taxa is key to address questions about patterns of genetic diversity and structure, hybridization and speciation. Here we present 14 newly developed microsatellite markers in the Spined toad (Bufo spinosus), assess their polymorphism in two Iberian populations and test for cross-species amplification in the closely related Common toad (Bufo bufo). We then use the 12 loci co-amplifying in both species to the study of a morphologically intermediate population (Moyaux) from the contact zone in northwest France as well as reference populations of the two species from both sides of the contact zone. Individuals from Moyaux had mtDNA haplotypes of the two species and were identified as hybrids in analyses with software NewHybrids. These results provide solid evidence for ongoing hybridization between B. bufo and B. spinosus, with no apparent restrictions to gene flow.
A widespread application of molecular markers across species is often hindered by a low cross-species amplification success, a high frequency of null alleles or a reduced level of polymorphism (Primmer et al., 2005). This is particularly true for hyper-variable markers, like microsatellites, where high mutation rates are both advantageous (because their polymorphism allows very fine-scale temporal and spatial resolution) and problematic (because mutations in flanking regions reduce cross-amplification success) (Primmer and Merilä, 2002). Thus, optimizing sets of applicable, reliable and informative markers, which are key to address questions about patterns of genetic diversity and structure, hybridization and speciation, involves selecting highly polymorphic markers that can nonetheless be consistently amplified and scored in long diverged taxa.
Recent molecular studies have resolved the phylogeny of the Bufo bufo (Linnaeus, 1758) species group, revealing the existence of four species and delimiting their respective ranges (Litvinchuk et al., 2008; Garcia Porta et al., 2012; Recuero et al., 2012; Arntzen et al., 2013a). As a result of these studies, B. spinosus Daudin, 1803 is defined as an Ibero-Maghrebian endemism, with populations in North Africa from Morocco to Tunisia, the Iberian Peninsula, the southern fringe of the British Isles (i.e., Jersey) and north of the Pyrenees across most of France, where it contacts B. bufo along a NW-SE line, roughly from Caen to Lyon (Arntzen et al., 2013b, 2014). Bufo bufo is also present in the Apennine and Balkan peninsulas and extends as far east as northern Kazakhstan and eastern Siberia (Agasyan et al., 2009).
Previous studies have used morphological characters, mtDNA and slowly evolving nuclear DNA markers to broadly delineate the contact zone between B. bufo and B. spinosus in France (Arntzen et al., 2013a, b, 2014). The results show an overall correspondence between morphology and molecules and congruence across molecular markers, with occasional discordance interpreted as resulting from incomplete lineage sorting, although hybridization-mediated gene flow cannot be ruled out.
Addressing questions about reproductive isolation versus hybridization and introgression in the two species requires characterization of molecular markers that amplify in both B. bufo and B. spinosus and the quantification of gene flow in areas of secondary contact. Bufo bufo has been the subject of a number of conservation genetics studies, especially in the UK and a number of polymorphic microsatellite loci are available (Scribner et al., 1994, 1997, 2001; Hitchings and Beebee, 1998; Brede and Beebee, 2004, 2006; Wilkinson et al., 2007). On the contrary, there are very few studies on B. spinosus and one of the reasons is low cross-amplification success of markers originally isolated from B. bufo (Brede et al., 2001; Martínez-Solano and González, 2008; Arntzen et al., 2014).
Here we describe a morphologically intermediate Bufo population from France with the help of a set of new polymorphic microsatellite markers isolated from a B. spinosus genomic library. We describe the variability of the newly developed markers in five populations of B. spinosus (two Iberian and three from western France), evaluate the cross-amplification success in two French populations of B. bufo and describe hybridization between the two species in a morphologically intermediate population in their contact zone.
Fig. 1. Sampling localities of Bufo bufo (blue) and B. spinosus (red) in France: 1) Saumur (n=12); 2) Perré (n=8); 3) Durtal (n=8); 4) Moyaux (n=22); 5) Audresselles (n=14); 6) Erloy (n=20). The population at Moyaux (open circle) is identified as a B. bufo x B. spinosus hybrid population (see results). Black dots represent populations studied in Arntzen et al. (2013b), whereas the dashed line represents the contact zone as approximated in that study.
Material and methods
A genomic library was generated from DNA of an adult female B. spinosus collected near Marbella, Málaga province, Spain, at the Sequencing Genotyping Facility, Cornell Life Sciences Core Laboratory Center, following the protocol described in Gutiérrez-Rodríguez et al. (2014) with minor modifications. Sequences containing microsatellites were scanned with Msatcommander (Faircloth, 2008) and sixty pairs of specific primers flanking regions that contained microsatellites were designed using the software Primerselect (DNASTAR, Inc.). Primers with an optimal melting temperature of 60ºC were selected to facilitate multiplexing.
PCR reactions were performed using 5 × GoTaq Flexi buffer PROMEGA®; 3 mM MgCl2 ; 0.3 mM dNTP; 0.3 μM of each primer; 0.5 U GoTaq Flexi DNA polymerase PROMEGA® and 1 μl of DNA, in a total volume of 15 μl. These initial tests included eight B. spinosus individuals from across the Iberian Peninsula and North Africa. PCR reactions consisted of initial denaturalization (95ºC; 5 minutes), 40 cycles of denaturalization (95ºC; 45 seconds), annealing (60ºC; 45 seconds) and extension (72ºC; 45 seconds) and final extension (72ºC; 10 minutes). PCR products were visualized in 2% agarose gels to verify amplification and check the negative controls for potential contamination.
Fourteen loci that consistently amplified in all individuals and were polymorphic as seen on gel images were selected for further screening. These loci were amplified in four multiplex reactions designed with Multiplex Manager 1.0 (Holleley and Geerts, 2009). Forward primers were labeled with fluorescent dyes 6-FAM, VIC, NED and PET. Multiplex reactions were performed with Type-it Microsatellite PCR® (Qiagen) kits in a total volume of 15 μl that includes 7.5 μl of Master Mix, 1.2 μl of Primer Mix with 0.3 μM of each primer and 5.3 μl of RNase-free H2 O. Multiplex reactions consisted of initial denaturalization (95ºC; 5 minutes), 30 cycles of denaturalization (95ºC; 30 seconds), annealing (during 90 seconds at different temperatures in each multiplex reaction: multiplex reactions 1, 2 and 4 had an annealing temperature of 60ºC and multiplex reaction 3 had a annealing temperature of 58ºC) and extension (72ºC; 30 seconds), and final extension (60ºC; 30 minutes).
The above markers were tested in 126 B. spinosus tadpoles collected in two different populations in central Spain (Madrid province): 94 from San Martín de la Vega (40.23º, -3.59º) and 32 from El Pardo (40.56º, -3.95º). Additionally, we sampled individuals in six localities in northern France encompassing the northwestern end of the range to test cross-amplification in B. bufo and evaluate the utility of the markers to study hybridization with B. spinosus. Our sampling includes three populations in the “spinosus” side of the contact zone: Étang du Perré (n=8, 47.53º, -0.02º), Durtal (n=8, 47.67º, -0.27º) and le Poitrineau, Saumur (n=12, 47.28º, -0.13º), two populations in the “bufo” side: Audresselles (n=14, 50.82º, 1.60º) and Autreppes (Erloy) (n=20, 49.92º, 3.83º) and a morphologically intermediate population (Moyaux, near Lisieux, n=22, 49.20º, 0.33º). Individuals from French populations were characterized with sequences of mtDNA and with diagnostic morphological characters as described in Arntzen et al. (2013b). DNA was extracted from approximately 0.25 g of tissue with Nucleospin® Tissue kits (Macherey-Nagel). Genotyping was performed on an ABI PRISM 3730 sequencer with the GeneScan 500 LIZ size standard (Applied Biosystems) and peaks were scored in GeneMapper 4.0 (Applied Biosystems).
Deviation from Hardy-Weinberg equilibrium and evidence of linkage disequilibrium across loci were tested for each population using Genepop 4.2 (Raymond and Rousset, 1995; Rousset, 2008). Using this software we also estimated the number of alleles (NA ), observed (HO ) and expected heterozygosity (HE ) for each locus and population. Significance values for multiple tests were calculated applying a sequential Bonferroni correction (Rice, 1989). The presence of null alleles was tested using Micro-Checker 2.2.3 (Van Oosterhout et al., 2004).
The software Genodive (Meirmans and Van Tienderen, 2004) was used to perform a Principal Components Analysis (PCA) on the microsatellite genotype data in order to summarize multi-allelic information in a single vector. Individual scores of Principal Component 1 were then plotted against latitude to help visualize the power of discrimination of the newly developed markers in species identification.
We conducted analyses of the French dataset with the software Structure (Pritchard et al., 2000). Runs consisted of a period of burn-in of 500,000 generations followed by 1,000,000 additional generations under an admixture model with allele frequencies correlated, for K values of 1 to 7 and 6 replicates per value of K. The optimal clustering scheme was selected through application of Evanno’s criterion (Evanno et al., 2005) as implemented in StructureHarvester (Earl and vonHoldt, 2012). Results were processed with the pipeline CLUMPAK (Kopelman et al., 2015, available at: http://clumpak.tau.ac.il/index.html).
The French dataset was also analyzed with NewHybrids software (Anderson and Thompson, 2002) to estimate the posterior probability that genetically sampled individuals fall into each of a set of six categories (parental classes, F1 , F2 and the two backcrosses with the parental species). The latter four categories were subsequently combined into a hybrid score, which is calculated as the sum of the probabilities of assignment to each category.
In 22 adult toads from Moyaux (18 males and two males and females from amplexed pairs) we studied a set of morphological characters that help to distinguish B. bufo from B. spinosus, 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 measurements Paratoid divergence (Pd=Pda/Pdp), Metatarsus Tubercle size (MTsize = MTl/SUl) and shape (MTshape = MTw/MTl). For details on the measurements, data analysis and reference populations see Arntzen et al. (2013b).
Fourteen loci consistently amplified and were polymorphic in the two populations of B. spinosus from Spain (San Martín de la Vega and El Pardo). The original sequences were deposited in GenBank under accession numbers KX271337-KX271350.
The number of observed alleles per locus across the 14 loci optimized in B. spinosus ranged from 2 to 14 (mean NA =8.1) in San Martín de la Vega and from 2 to 12 (mean NA =6.4) in El Pardo. Observed and expected heterozygosities ranged from 0.20 to 0.91 (mean HO =0.66) and 0.31 to 0.86 (mean HE =0.71) respectively in San Martín de la Vega; and from 0.13 to 0.97 (mean HO =0.66) and 0.12 to 0.89 (mean HE =0.70) in El Pardo. Micro-Checker inferred the likely presence of null alleles in loci Bspi3.19, Bspi3.20 and Bspi4.25, associated with deviations from Hardy Weinberg equilibrium due to heterozygote deficits in both populations (Bspi3.20 and 4.25) and San Martín de la Vega (Bspi3.19); therefore, caution should be taken when using these markers for population studies in B. spinosus. Four locus x locus combinations significantly deviated from linkage disequilibrium (after accounting for the Bonferroni correction, P<0.001) in Iberian populations (two in each population, with no consistent pattern across populations).
Locus Bspi4.07 only amplified in B. spinosus and locus Bspi3.20 did not consistently amplify in B. bufo. Therefore, we calculated the number of alleles (NA ), observed (HO ) and expected (HE ) heterozygosity in the remaining populations using twelve microsatellites: Bspi3.11, Bspi4.24, Bspi4.25, Bspi4.16, Bspi4.30, Bspi4.27, Bspi4.28, Bspi4.29, Bspi3.02, Bspi3.19, Bspi3.26 and Bspi4.14. Results are shown in Table 2.
Loci Bspi4.16, 3.02, 3.11 and 4.30 were monomorphic in the two B. bufo populations analyzed. Significant deviations from Hardy-Weinberg equilibrium (H1 : heterozygote deficit) were detected for loci Bspi3.19 in Durtal, Bspi4.28 in Perré and Bspi4.24 in Moyaux and Erloy. Significant linkage disequilibrium were only observed in Moyaux (four locus x locus combinations: Bspi4.24 × Bspi4.27, Bspi4.24 × Bspi4.29, Bspi4.28 × Bspi4.29 and Bspi4.28 × Bspi4.14).
All individuals analyzed from Perré (n=8), Durtal (n=5) and Saumur (n=10) had the characteristic mtDNA profile of B. spinosus, whereas all individuals from Audresselles (n=22) and Erloy (n=28) had B. bufo mtDNA. In Moyaux, six individuals had B. spinosus mtDNA and 16 individuals had B. bufo mtDNA haplotypes. We found two interspecific amplexed pairs in our sample, a male B. spinosus mating with a female B. bufo and a male B. bufo with a female B. spinosus.
The first principal component of the microsatellite genotype data explained 25.8% of variation; when individual scores were plotted against geography, a pattern emerged of three groups corresponding to 1) the populations of Saumur, Durtal and Perré; 2) Moyaux and 3) Audresselles and Erloy (Fig. 2).
Structure analyses grouped individuals in two clusters (Fig. 2), with high assignment probabilities. A first cluster, corresponding to B. spinosus, included the populations Saumur (average probability of assignment, PA : 0.99), Perré (PA : 0.97), Durtal (PA : 0.99), whereas the second cluster (B. bufo) included the populations Erloy (PA : 0.99) and Audresselles (PA : 0.99). Individuals from Moyaux were generally admixed, with a higher probability of assignment to the B. bufo cluster (PA : 0.85, range 0.54-0.98, Fig. 2).
According to the results of NewHybrids analyses (Fig. 2) most individuals in the populations of Saumur, Perré and Durtal correspond to B. spinosus (population average proportions of 1.00, 0.63 and 0.88 respectively), while the remainder of the individuals from these populations (three individuals in Perré and one individual in Durtal) was assigned to the F2 class. In the populations Audresselles and Erloy all individuals were classified as B. bufo with high probability, whereas in Moyaux, most individuals (64%) were classified as backcrosses with B. bufo, 27% were classified as F2 individuals and the remaining 9% individuals were classified as B. bufo (Fig. 2).
The morphometric characteristics of the population Moyaux were intermediate between those of B. bufo and B. spinosus, as shown by the analysis of the diagnostic characters paratoid angle / paratoid divergence and metatarsus tubercle shape and size (Fig. 3).
Fig. 3. Morphometric characteristics of adult toads from the population Moyaux (solid round symbols, 20 males, two females) relative to reference individuals from northern France, with ellipses showing one standard deviation around the mean, for B. bufo (shaded) and B. spinosus (unshaded). Reference data are taken from Arntzen et al. (2013b).
The new set of 14 polymorphic microsatellite markers is a valuable tool for population studies on B. spinosus. The species, formerly abundant throughout its Iberian range, is now undergoing a slow but sustained decline (Ortiz Santaliestra, 2014), and molecular tools can provide relevant information about patterns of genetic diversity and structure to guide conservation efforts. Furthermore, the high cross-amplification rates and observed levels of polymorphism in B. bufo make these markers valuable for fine-scale studies of the contact zone with B. spinosus as well as to complement previously developed markers in local or regional studies on B. bufo. Additional populations of the latter species should be assayed with the new markers to better characterize actual polymorphism, which was certainly underestimated in this study (for instance, markers Bspi3.02, 3.11 and 4.30 are monomorphic in the two French populations of B. bufo analyzed here, but additional polymorphism was found in a reduced sample of eastern European samples, data not shown).
Based on the new molecular data, the morphologically intermediate population of Moyaux is shown to represent a hybrid population, where mtDNA haplotypes of the two species can be found in syntopy. Co-occurrence of mtDNA haplotypes of the two species in a single site is rare, but had been previously reported (Recuero et al., 2012; Arntzen et al., 2013a). However, microsatellites also reveal a clearly admixed nuclear DNA genetic pool at Moyaux, suggestive of interspecific gene flow. This is thus the first solid evidence of hybridization between the two species (see Arntzen et al., 2013a for a discussion about alternative interpretations of allozyme data on the taxonomic status of B. bufo and B. spinosus). While more detailed transects need to be studied to better characterize the contact zone in terms of its relative width or patterns of symmetry/asymmetry of gene flow across species, our preliminary results show a consistent signal of hybridization across independent sets of markers, including species-diagnostic morphological characters, mtDNA haplotypes and microsatellite loci. Thus, at least some of the previously observed instances of nuclear DNA haplotype sharing in slowly evolving nuclear markers (POMC, RAG1) may also result from interspecific gene flow rather than from incomplete lineage sorting alone (Arntzen et al., 2013a, 2014).
The lack of F1 hybrids and the high frequency of F2 and backcrossed individuals in the hybrid population at Moyaux suggest no restrictions to gene flow (random mating), which is in accordance with the observation of amplectant mates with mtDNA from the two different species (in both directions). This contrasts with the high incidence of deviations from linkage disequilibrium in the hybrid population, which may be indicative of assortative mating or selection rather than result from random demographic effects, because hybridization tends to break down species-specific physical associations between markers in hybrid zones and create new ones (Jiggins and Mallet, 2000). In any case, the high frequency of admixed individuals and relative scarcity of representatives of the parental taxa indicates this would be a relatively old contact zone, with many generations of hybrids gradually increasing their frequency in the population and becoming the dominant parental class, although more detailed, replicated transects are required to learn more about the potential role of selection in shaping the contact zone.
Whereas the distinction of B. spinosus as a separate species from B. bufo has long been subject to confusion, close inspection of the accumulating morphological and molecular evidence for species differentiation indicates that the two species have diagnostic morphological features and are highly distinct in their mitochondrial and nuclear DNA, except, as shown here, in parts of the contact zone where some hybridization occurs (including some of the reference B. spinosus populations, which according to NewHybrids results may also include some hybrid individuals, see Fig. 2). More detailed analyses of different sections of this contact zone are required to further delineate it and investigate whether observed patterns are generalizable. Additionally, comparative data about demographic and life history traits across species may be relevant to understand the evolutionary processes maintaining species boundaries in the face of gene flow. For instance, the skeletochronological work of Hemelaar (1988) showed remarkable differences in age at maturity and growth rates in Bufo populations from different latitudes. Of the five populations compared, one from southern France represents B. spinosus as inferred from recent molecular studies, whereas the other four would correspond to B. bufo on geographic grounds. Several characteristics of the French population stand out in comparison with populations of B. bufo, for instance much faster growth rates, larger maximum body size and earlier age at maturity, suggesting genetically, rather than environmentally induced differences (Hemelaar, 1988). Also, B. spinosus appears to start breeding a few weeks earlier than B. bufo and the breeding season is also more prolonged than in B. bufo, which is usually considered an explosive breeder. These differences, in concert with potential mechanisms of species recognition, which do not seem to play a role in preventing interspecific mating in Moyaux but may be apparent elsewhere, can shape patterns of interspecific gene flow in the contact zone, ultimately determining its evolutionary fate.
The authors thank G. Rodríguez, G. Sánchez-Montes, A. Vallvé, J. van Alphen, A. Zuiderwijk and Agentes Forestales de la Comunidad de Madrid for help in collecting individuals and P. Ochoa for technical support. Two anonymous reviewers provided valuable feedback on a previous draft. This research was funded by grants CGL2008-04271-C02-01/BOS and CGL2011-28300 (Ministerio de Ciencia e Innovación), Ministerio de Economía y Competitividad and FEDER) and PPII10-0097-4200 (Junta de Comunidades de Castilla la Mancha) to IMS. J. Gutiérrez-Rodríguez was supported by the Consejo Superior de Investigaciones Científicas of Spain (CSIC) and the European Social Fund (ESF) with a JAE-pre PhD fellowship. IMS was funded by the project ‘Biodiversity, Ecology and Global Change’, co-financed by North Portugal Regional Operational Programme 2007/2013 (ON.2–O Novo Norte), under the National Strategic Reference Framework (NSRF), through the European Regional Development Fund (ERDF) and by funding from the Spanish Severo Ochoa Program (SEV-2012-0262). IMS’s stay in Leiden was supported by a Naturalis Temminck fellowship.
Received: 12 January 2016
Revised and accepted: 15 June 2016
Published online: 14 February 2017
Editor: A. Ivanović
Agasyan A, Avisi A, Tuniyev B, Crnobrnja Isailovic J, Lymberakis P, Andrén C, Cogalniceanu D, Wilkinson J, Ananjeva N, Üzüm N, Orlov N, Podloucky R, Tuniyev S, Kaya U. 2009. Bufo bufo. The IUCN Red List of Threatened Species. Version 2014.2. [http://www.iucnredlist.org]
Anderson EC, Thompson EA. 2002. A model-based method for identifying species hybrids using multilocus genetic data. Genetics 160: 1217-1229.
Arntzen JW, Recuero E, Canestrelli D, Martínez-Solano I. 2013a. How complex is the Bufo bufo species group? Molecular Phylogenetics and Evolution 69: 1203-1208.
Arntzen JW, McAtear J, Recuero E, Ziermann JM, Ohler A, van Alphen J, Martínez-Solano I. 2013b. Morphological and genetic differentiation of Bufo toads: two cryptic species in Western Europe (Anura, Bufonidae). Contributions to Zoology 82: 147-169.
Arntzen JW, Wilkinson JW, Butôt R, Martínez-Solano I. 2014. A new vertebrate species native to the British Isles: Bufo spinosus Daudin, 1803 in Jersey. Herpetological Journal 24: 209-216.
Brede E, Rowe G, Trojanowski J, Beebee TJC. 2001. Polymerase chain reaction primers for microsatellite loci in the Common Toad Bufo bufo. Molecular Ecology Notes 1: 308-310.
Brede EG, Beebee TJC. 2004. Contrasting population structures in two sympatric anurans: implications for species conservation. Heredity 92: 110-117.
Brede EG, Beebee TJC. 2006. Large variations in the ratio of effective breeding and census population sizes between two species of pond-breeding anurans. Biological Journal of the Linnean Society 89: 365-372.
Daudin FM. 1803. Histoire Naturelle, Générale et Particulière des Reptiles; Ouvrage Faisant suit à l’Histoire Naturelle Générale et Particulière, Composée par Leclerc de Buffon; et Rédigée par C.S. Sonnini, Membre de Plusieurs Sociétés Savantes. Volume 8. Paris: F. Dufart.
Earl DA, vonHoldt BM. 2012. Structure Harvester: a website and program for visualizing Structure output and implementing the Evanno method. Conservation Genetics Resources 4: 359-361.
Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters of individuals using the software Structure: a simulation study. Molecular Ecology 14: 2611-2620.
Faircloth BC. 2008. Msatcommander: detection of microsatellite repeat arrays and automated, locus-specific primer design. Molecular Ecology Resources 8: 92-94.
Garcia-Porta J, Litvinchuk SN, Crochet PA, Romano A, Geniez PH, Lo-Valvo M, Lymberakis P, Carranza S. 2012. Molecular phylogenetics and historical biogeography of the west-palearctic common toads (Bufo bufo species complex). Molecular Phylogenetics and Evolution 63: 113-130.
Gutiérrez Rodríguez J, Salvi D, Geffen E, Gafny S, Martínez-Solano I. 2014. Isolation and characterisation of novel polymorphic microsatellite loci in Iberian painted frogs (Discoglossus galganoi and D. jeanneae), with data on cross-species amplification in Discoglossus and Latonia (Alytidae). Herpetological Journal 24: 261-265.
Hemelaar A. 1988. Age, growth and other population characteristics of Bufo bufo from different latitudes and altitudes. Journal of Herpetology 22: 369-388.
Hitchings SP, Beebee TJC. 1998. Loss of genetic diversity and fitness in common toad (Bufo bufo) populations isolated by inimical habitat. Journal of Evolutionary Biology 11: 269-283.
Holleley CE, Geerts PG. 2009. Multiplex Manager 1.0: a crossplatform computer program that plans and optimizes multiplex PCR. BioTechniques 46: 511-517.
Jiggins CD, Mallet J. 2000. Bimodal hybrid zones and speciation. Trends in Ecology and Evolution 15: 250-255.
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. 2015. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Molecular Ecology Resources 15: 1179-1191.
Linnaeus C. 1758. Systema Naturae per Regna Tria Naturae, Secundum Classes, Ordines, Genera, Species, cum Characteribus, Differentiis, Synonymis, Locis. 10th Edition. Volume 1. Stockholm: L. Salvii.
Litvinchuk S, Borkin L, Skorinov D, Rosanov J. 2008. A new species of common toads from the Talysh mountains, south-eastern Caucasus: genome size, allozyme, and morphological evidences. Russian Journal of Herpetology 15: 19-43.
Martínez-Solano I, González EG. 2008. Patterns of gene flow and source-sink dynamics in high altitude populations of the common toad Bufo bufo (Anura: Bufonidae). Biological Journal of the Linnean Society 95: 824-839.
Meirmans PG, Van Tienderen PH. 2004. Genotype and Genodive: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes 4: 792-794.
Ortiz Santaliestra ME. 2014. Sapo común – Bufo spinosus. In: Enciclopedia Virtual de los Vertebrados Españoles. A. Salvador & I. Martínez-Solano (Eds.). Madrid: Museo Nacional de Ciencias Naturales. [http://www.vertebradosibericos.org/]
Primmer CR, Merilä J. 2002. A low rate of cross-species microsatellite amplification success in Ranid frogs. Conservation Genetics 3: 445-449.
Primmer CR, Painter JN, Koskinen MT, Palo JU, Merilä J. 2005. Factors affecting avian cross-species microsatellite amplification. Journal of Avian Biology 36: 348-360.
Pritchard JK, Stephens M, Donnelly P. 2000. Inference of population structure using multilocus genotype data. Genetics 155: 945-959.
Raymond M, Rousset F. 1995. Genepop (version 1.2): population genetics software for exact tests and ecumenicism. Journal of Heredity 86: 248-249.
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. 2012. Multilocus species tree analyses resolve the radiation of the widespread Bufo bufo species group (Anura, Bufonidae). Molecular Phylogenetics and Evolution 62: 71-86.
Rice WR. 1989. Analyzing tables of statistical tests. Evolution 43: 223-225.
Rousset F. 2008. Genepop’007: a complete reimplementation of the Genepop software for Windows and Linux. Molecular Ecology Resources 8: 103-106.
Scribner KT, Arntzen JW, Burke T. 1994. Comparative analysis of intra-and interpopulation genetic diversity in Bufo bufo, using allozyme, single-locus microsatellite, minisatellite, and multilocus minisatellite data. Molecular Biology and Evolution 11: 737-748.
Scribner KT, Arntzen JW, Burke T. 1997. Effective number of breeding adults in Bufo bufo estimated from age-specific variation at minisatellite loci. Molecular Ecology 6: 701-712.
Scribner KT, Arntzen JW, Cruddace N, Oldham R, Burke T. 2001. Environmental correlates of toad abundance and population genetic diversity. Biological Conservation 98: 201-210.
van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. 2004. Micro-Checker: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4: 535-538.
Wilkinson JW, Beebee TJC, Griffiths RA. 2007. Conservation genetics of an island toad: Bufo bufo in Jersey. Herpetological Journal 17: 192-198.