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