- Research article
- Open Access
Strong reproductive barriers in a narrow hybrid zone of West-Mediterranean green toads (Bufo viridissubgroup) with Plio-Pleistocene divergence
BMC Evolutionary Biologyvolume 10, Article number: 232 (2010)
One key question in evolutionary biology deals with the mode and rate at which reproductive isolation accumulates during allopatric speciation. Little is known about secondary contacts of recently diverged anuran species. Here we conduct a multi-locus field study to investigate a contact zone between two lineages of green toads with an estimated divergence time of 2.7 My, and report results from preliminary experimental crosses.
The Sicilian endemic Bufo siculus and the Italian mainland-origin B. balearicus form a narrow hybrid zone east of Mt. Etna. Despite bidirectional mtDNA introgression over a ca. 40 km North-South cline, no F1 hybrids could be found, and nuclear genomes display almost no admixture. Populations from each side of the contact zone showed depressed genetic diversity and very strong differentiation (FST = 0.52). Preliminary experimental crosses point to a slightly reduced fitness in F1 hybrids, a strong hybrid breakdown in backcrossed offspring (F1 x parental, with very few reaching metamorphosis) and a complete and early mortality in F2 (F1 x F1).
Genetic patterns at the contact zone are molded by drift and selection. Local effective sizes are reduced by the geography and history of the contact zone, B. balearicus populations being at the front wave of a recent expansion (late Pleistocene). Selection against hybrids likely results from intrinsic genomic causes (disruption of coadapted sets of genes in backcrosses and F2-hybrids), possibly reinforced by local adaptation (the ranges of the two taxa roughly coincide with the borders of semiarid and arid climates). The absence of F1 in the field might be due to premating isolation mechanisms. Our results, show that these lineages have evolved almost complete reproductive isolation after some 2.7 My of divergence, contrasting sharply with evidence from laboratory experiments that some anuran species may still produce viable F1 offspring after > 20 My of divergence.
One key question in evolutionary biology deals with the mode and rate at which reproductive isolation accumulates during allopatric speciation [for overview: ]. Johns and Avise  estimated the average mitochondrial DNA (mtDNA)-based genetic distance between congeneric species in amphibians to be > 7.0 My, suggesting absence of natural hybridization in taxa of that age. A few major results on intrinsic reproductive isolation in anurans come from artificial hybridization experiments. Sasa et al.  reported hybrid sterility or inviability in 46 frog species to be positively correlated with Nei's genetic distance (allozymes). Measuring albumin distances among 50 species pairs, Wilson et al.  showed that frogs could still produce viable hybrids with an average immunological distance of 7.4% (= ca. 21 My). Using Blair's  crossing experiments in Bufo, Malone & Fontenot  showed the hatching success, the number of larvae produced, and the percentage of tadpoles reaching metamorphosis to be inversely related with genetic divergence, some metamorphosing offspring being still produced with a distance of 8% (mtDNA). All of these laboratory data suggest that reproductive isolation increases gradually with phylogenetic distance, presumably driven by complex genomic processes rather than by a few speciation genes, and that very large time scales (in the order of tens of millions of years) are required to achieve hybrid infertility or inviability.
Under natural conditions, however, reproductive isolation could arise much earlier than detected in the laboratory. In frogs, as in many other taxa, "surveys of natural hybrid zones (...) in the field are needed to complement laboratory-based studies to establish the significance and strength of specific barriers in nature" . Little is known about secondary contact in allopatrically diverged lineages of anurans, where reproductive isolation may quickly arise as a result of reinforcement , in addition to genetic drift and local adaptation. Extant studies of contact zones in anurans have mostly focused on hybrid fitness [9, 10] or on mechanisms of pre- or post-mating isolation [9–15]. Such studies classically relied on allozymes [e.g. [9, 11–13]] or (more recently) nuclear and mitochondrial DNA markers [e.g. [10, 14, 15]], but often lack any molecular-based estimates of divergence times, which are at best inferred from geological information. Some phylogeographic studies include molecular-based estimates of divergence time [e.g. [16–20]], but very few have combined such estimates with multi-locus transect approaches to infer the time required to reach reproductive isolation in natural contexts [e.g. [8, 21, 22]].
The current study focuses on Palearctic green toads [Bufo viridis subgroup, ]. After range-wide phylogeographic analyses, secondary contact zones of clades were predicted [18, 23], in which possible hybridization can be examined using fast evolving molecular markers. To do this, we recently developed microsatellites for two West-Mediterranean species : B. balearicus (Boettger 1880; Peninsular Italy, north-eastern Sicily, Corsica, Sardinia, Balearic Islands), and B. siculus [; endemic to Sicily, Figure 1]. Using a Bayesian-coalescence approach (mtDNA control region and 16 S rRNA), divergence time for the two species was estimated to late Pliocene (2.7 My), with a range from the early Pliocene (4.9 My) to Pleistocene (1.1 My) . A single record of Italian mainland-origin B. balearicus in north-eastern Sicily  suggests their recent (late Pleistocene) invasion into Sicily, where they may secondarily meet the endemic B. siculus.
In this work, we combined mitochondrial and nuclear intronic sequences with multilocus microsatellite markers to examine (i) whether B. siculus and B. balearicus meet each other in north-eastern Sicily, (ii) if so, whether these two closely related species hybridize, and (iii) in such a case, what are the patterns of hybridization. In parallel, we conducted limited and preliminary experimental crosses to help interpreting field data.
Nuclear and mitochondrial DNA sequences and mitotyping
Both of the phylogenetic trees built from mitochondrial (D-loop) and nuclear (Tropomyosine intron) DNA sequences show two highly homogeneous and strongly distinct clades (Figure 2), corresponding to B. balearicus and B. siculus. The tropomyosine tree displays a clear geographic pattern: the balearicus clade includes all individuals from mainland Italy (populations 1 to 8, Figure 1, Table 1) and north-eastern Sicily, southwards to population 14 (east coast), while the siculus clade includes individuals from western and southern Sicilian populations, from population 15 (East coast) south-westwards. This pattern points to a very narrow contact zone separating populations 14 and 15, between the Mount Etna and the Ionian coast (Figure 1b).
The mtDNA clades also show a clear geographic signal, with however, some overlap. Populations from mainland Italy (pop. 1 to 8) and north-eastern Sicily (pop. 9 to 12) present only balearicus haplotypes, and populations from western and southern Sicily (pop. 16 to 24) only siculus, but haplotypes from both clades are found in populations 13 to 15, around the contact zone identified with tropomyosine.
These phylogenetic trees also provide evidence for past hybridization, as revealed by cytonuclear disequilibria (see highlighted individuals in Figure 2): one individual from pop. 14 possesses balearicus tropomyosine alleles but a siculus mtDNA-haplotype, while three individuals from pop. 15 present siculus tropomyosine alleles but balearicus mtDNA haplotypes.
These patterns of mitochondrial distribution were widely confirmed by larger-scale mitotyping (Figure 1). All populations on the Apennine Peninsula (pop. 1 to 8) and four populations (pop. 9 to 12) from the North-East of Sicily presented only B. balearicus haplotypes. All populations from western and southern Sicily (pop. 16 to 22) and the two islands off the coast of western Sicily (pop. 23, 24) presented only B. siculus haplotypes. In the three populations east of Mount Etna (pop. 13 to 15), both B. balearicus and B. siculus haplotypes were present, with a marked north-south cline (Table 2): The frequency of balearicus haplotypes declined from 93.75% in Calatabiano (pop. 13) to 68% in Giarre (pop. 14) and 50% in Gravina (pop. 15), down to 0% in Misterbianco (pop. 16).
Autosomal microsatellites and population-genetics analyses
There was no evidence for allelic dropout from any locus in any population. Null alleles at low frequencies were detected (and corrections performed) in one population each for loci C203 (pop. 14) and D105 (pop. 13), and in two populations each for loci C218 (pop. 6, 22), C223 (pop. 18, 21) and D5 (pop. 17, 18). Tests for linkage disequilibrium between loci (after sequential Bonferroni corrections) revealed four significant combinations (Bcal μ10 × C203 for pop. 9, D5 × Bcal μ10 for pop. 15, C223 × Bcal μ10 for pop. 18, D105 × C205 for pop. 23). A few B. siculus populations showed some heterozygote deficit (Additional file 1), presumably due to sampling design (substructures may arise when pooling tadpoles or adults from several nearby ponds).
Bayesian clustering assignment using STRUCTURE  largely confirmed the nuclear information from Tropomyosine. All populations from Sicily were clearly grouped into two clusters [K = 2; ] corresponding to B. balearicus and B. siculus gene pools respectively (Figure 3). All individuals from populations 9 to 14 were assigned to B. balearicus, while all individuals from populations 15 to 24 were assigned to B. siculus. Ten F1-hybrids from an experimental cross between a female balearicus (pop. 11) and a male siculus (pop. 22) were correctly assigned a 50% probability of belonging to either balearicus or siculus (pop. 25). Surprisingly, the two populations north and south of the contact zone (pop. 14 and 15) did not show any sign of hybridization or gene flow, despite harboring mtDNA from both clades. All individuals from population 14 were assigned with a 100% probability to balearicus, and all individuals from the population 15 with 100% probability to siculus. As a matter of fact, potential hybrids appear very few (altogether four individuals with assignment probabilities lower than 90% to either parental species), and largely backcrossed (assignment probabilities to the alternative parental species lower than 25%).
This pattern was confirmed by NEWHYBRIDS , which, when including all Sicilian populations, correctly assigned all experimental crosses as F1-hybrids, and identified four wild-caught individuals as possible F2-hybrids (two each from pop. 13 and 18, details in Additional file 2). When focusing on populations where hybrids occurred or were likely to do so (pop. 12 to 16 and 18), while pre-assigning pop. 9 to 11 and 17 as pure B. balearicus and B. siculus, respectively, no nuclear hybrids were detected. Finally, diagnostic alleles also suggested faint signs of past hybridizations (Additional file 3). We found B. siculus alleles in three individuals assigned as B. balearicus using STRUCTURE (one from pop. 13 and two from pop. 14), and B. balearicus alleles in six individuals assigned as B. siculus (four in pop. 15 and two in pop. 18). All of these analyses concur to suggest limited events of nuclear introgression.
In order to fine-tune our analysis of potential gene flow, we performed separate STRUCTURE analyses for B. balearicus and B. siculus populations. Substructure within B. balearicus was best explained with K = 3 (Figure 1b). Cluster b1 contained individuals from mainland Italy, cluster b2 individuals from north-eastern Sicilian populations, and cluster b3 individuals from populations close to the contact zone. Interestingly population 8 (tip of Calabria) showed admixture of mainland Italy and Sicilian genotypes. Substructure within B. siculus was similarly best explained with K = 3 (Figure 1b). Cluster s2 contained individuals from the vast majority of populations (including off coast islands, pop. 23 and 24), except for populations 21 and 22 (north-western coast of Sicily, cluster s1) and population 15 at the contact zone (cluster s3).
Hence, in both species, the populations close to the hybrid zone (showing coexistence of mtDNA haplotypes) form a cluster of their own. However, this pattern is clearly not generated by nuclear gene flow between the two species. Indeed, from pair-wise FST values, the strongest differentiation (FST = 0.52, Additional file 4) actually occurs between the two populations (pop. 14 and 15) from each side of the contact zone, as compared to an average value of 0.32 between allospecific populations (and 0.18 between conspecific populations). This unexpected result was confirmed by a principal component analysis [PCAGEN; ] aimed at extracting factors maximizing genetic differentiation among populations (Figure 4). Two factors turn out to be significant, explaining respectively 40.4 and 13.1% of the total differentiation (FST). The first one accounts for the contrast between B. siculus (left) and B. balearicus (right). The three B. balearicus clusters identified with STRUCTURE differentiate along this axis, with lowest values for b1 (mainland Italy, pop. 1 to 7), and highest values for b3 (pop. 12 to 14, close to the contact zone). The three B. siculus clusters differentiate mostly on the second factor, with lowest values for s1 (north-western coast, pop. 21 and 22) and highest values for s3 (pop. 15, close to the contact zone).
The spread of clusters correlates with geography, which translates into some isolation by distance. The relationship between genetic differentiation and geographic distance is strong and significant in both species when dropping the three populations (pop. 13 to 15) at the contact zone (r = 0.81, R2 = 66%, p = 0.0026 for B. balearicus, and r = 0.41, R2 = 17%, p = 0.035 for B. siculus), but drastically reduced when including these three populations, due to strong differentiation over short geographic distances (r = 0.23, R2 = 8.72%, p = 0.25 for B. balearicus, and r = 0.21, R2 = 4.65%, p = 0.47 for B. siculus).
This enhanced differentiation between populations close to the contact zone correlates with increased genetic drift and loss of diversity. Genetic diversity in B. balearicus populations decreases from Hs = 0.74 in mainland Italy (pop. 5 and 6) to 0.54 in Calabria and Northern Sicily, down to 0.38 at the contact zone (pop. 14). Similarly (though to a lesser extent), genetic diversity in B. siculus populations decreases from Hs = 0.75 South and West of the Mount Etna (pop. 17 and 18) to 0.62 in populations closer to the contact zone (pop. 15 and 16; Additional file 1, see also  for representative populations).
From an F1-cross B. balearicus x B. siculus obtained in spring 2007, about 80% offspring were viable and developed normally (Table 3). The remaining 20% did not hatch or produced malformed, dwarfed and/or leucistic larvae (Figure 5c-f). Most of these died at early stages or during metamorphosis (four-legged stage), and a few ones survived as never-metamorphosing "giant" tadpoles (Figure 5f). The reciprocal cross (B. siculus x B. balearicus) showed much lower survival after metamorphosis (Table 3).
We raised about 160 F1-metamorphs from the 2007 cross to a snout-vent-length of ca. 2 cm, and kept 50 of them until secondary sexual characters became visible (paler coloration and nuptial pads of males). Though sex ratio was approximately even, we noticed about 30% of dwarfed F1-males that reached only about two-thirds the size of normal individuals. Ten males and ten females were further raised until maturity (Figures 5g-h). In spring 2009, we used one F1 of each sex to produce F2-hybrids (F1 × F1) and reciprocal backcrosses with either B. siculus or B. balearicus (one new, wild-caught male and one female each, Figure 6). F2-hybrids turned out to be unviable, with all tadpoles dying a few days after hatching. While 200 out of 328 tadpoles from the backcrosses were still alive two months after spawning, they presented, a number of developmental abnormalities, including greenish individuals and a bimodal size distribution within the same cross (Figure 7), and suffered from dramatic mortality at later stages, with two individuals only surviving after metamorphosis.
Our study shows that two distinct lineages of green toads (Bufo viridis subgroup) occur parapatrically in Sicily. These lineages have diverged some 2.7 My ago , a time frame long enough to allow significant differentiation on both mitochondrial and nuclear DNA sequences (Figure 2). As our study further shows, this divergence was also sufficient to bring the speciation process close to completion.
On the one hand, the endemic B. siculus, of North-African origin [home of its sister clade B. boulengeri; ], occupies the western and southern parts of Sicily, plus two small islands off the north-western coasts (Ustica and Favignana) . On the other hand, B. balearicus occupies the north-eastern part of Sicily. Based on their geographical localization and patterns of genetic similarity with mainland Italy, we infer that these Sicilian B. balearicus populations recently originated from close-by Calabrian populations. Faunal exchange across the Strait of Messina [including amphibians ] are well documented for the Upper Pleistocene . From our genetic analyses, these two species nowadays meet at the eastern coast of Sicily, between the Mount Etna and the Ionian Sea. We cannot exclude that another contact exists along the North coast (north-west of Mount Etna), but could not find any currently occupied site in this area despite thorough examination.
Though very restricted, the documented contact zone shows signs of past hybridization, with differential introgression patterns depending on markers. Mitochondrial alleles show a clear North-South cline, where the frequency of balearicus haplotypes progressively decreases from 94% (pop. 13, Calatabiano) to 0% (pop. 16, Misterbianco) over a distance of ca 40 km. Cytonuclear disequilibrium occurred in individuals from both species, pointing to a two-way introgression. This presumably involved symmetric events of hybridization, followed by backcrossing of fertile F1-females with their paternal species. Though no F1-hybrids were detected in the field, the occurrence of rare and symmetric events of hybridization was confirmed by a few backcrosses (F2 or more) identified via STRUCTURE and NEWHYBRIDS in both B. siculus and B. balearicus populations, as well as a two-way leak of diagnostic nuclear alleles.
However, nuclear introgression was surprisingly low overall. The transition between tropomyosine alleles from the two clades was abrupt, occurring at some point between populations 14 (Giarre) and 15 (Gravina), separated by just 16 km. The same holds for autosomal microsatellites in general, since STRUCTURE assigned (with 100% probability) all individuals from pop. 14 to balearicus, and all individuals from pop. 15 to siculus (Figure 3). The sharpness of this transition is underlined by the patterns of genetic differentiation: The two populations each side of the contact zone (pop. 14 and 15), though harboring a mix of mitochondrial haplotypes from both lineages, display the highest differentiation value observed in the study area (pairwise FST = 0.52). Populations at the contact zone (clusters b3 and s3) are actually the ones most differentiated on the first PCAGEN factors (Figure 4).
Genetic drift certainly plays a role in this strong local differentiation. The B. balearicus populations at the contact zone represent the front wave of a recent expansion, as evidenced by the drastic decrease in genetic diversity from mainland Italy (Hs = 0.74) to the southernmost populations (Hs = 0.38 in population 14). Drift is certainly further amplified by the geographic localization of the contact zone: The geographical bottleneck between Mount Etna and the Ionian Sea creates a peninsular situation, largely isolating populations at the contact zone from conspecifics. This induces a strong differentiation over a small geographic scale, which somewhat blurs the overall clear pattern of isolation by distance observed in both species.
Drift might also partly account for the marked contrast between mitochondrial and nuclear introgression. Mitochondrial markers have low effective sizes (about one quarter of nuclear markers), and are therefore more prone to introgression. In small hybridizing populations, mtDNA might sometimes get fixed into foreign taxa via drift or founder effects, even with little or no gene flow at nuclear loci . This contrast should be further amplified if females display lower effective size than males, which occurs when dispersal is male biased [as expected in polygynous mating system such as found in bufonids; [32–36]]. Effective population sizes at the front wave of expansions (or in any poorly connected population) largely depend on gene flow from incoming immigrants, and thus on the mode of inheritance of markers when dispersal is sex biased [37, 38].
In addition to drift and dispersal, selection often contributes crucially to the maintenance of narrow hybrid zones . We think that the one under study is no exception, and that selection against nuclear introgression is very likely to also play a role. First, F1-hybrids from our experimental crosses showed reduced fitness (Table 3). Outcomes from reciprocal crosses were asymmetric, which often occurs in interspecific crosses (as documented e.g. in hybrids between green toads and Bufo calamita  or B. bufo [41, 42]. If the small size and altered coloration observed in males from the B. balearicus × B. siculus cross translate into lower survival or fertility in the field, then nuclear markers are indeed expected to show less introgression than mtDNA. Second, survival was drastically affected both in backcrosses and in F2-hybrids (F1 × F1), which should strongly reduce introgression. Third, the currently known ranges of both taxa in Sicily roughly coincide with the borders of semiarid (B. balearicus) and arid (B. siculus) climates . Adaptation of these two genomes to different climates may select against hybrids, and potentially stabilize the contact zone.
The sex-specific phenotypic effects in the F1 as well as the dramatic hybrid breakdown observed both in backcrosses and in the F2 generation (F1 × F1) were expected from the classical Dobzhansky-Muller model of speciation , arising from the confrontation of incompatible genes and the disruption of co-adapted sets of genes. Except for sex chromosomes in the heterogametic sex, F1-hybrids inherit complete sets of genes from both parental species and should thus suffer little from co-adaptation losses. We do not know, however, whether sex-specific differences in F1 phenotypes (dwarfed males) conform to Haldane's rule , because sex-determination mechanisms are unknown for Sicilian green toads (as for most amphibians). Results from related species suggest a XX/XY-system for B. variabilis (as "B. viridis") from Asia Minor [46, 47], but a ZZ/ZW-system for taxonomically undefined green toads from Moldavia [B. viridis/B. variabilis contact zone?, ]. Contrasting with F1, backcrosses and F2 (F1 × F1) inherit imbalanced numbers of genes from both parental species due to recombination in F1, and may thus lack crucial alleles at complementary loci. Inbreeding presumably also played a role in our F2 crosses, which may partially explain the additional mortality relative to the breakdown observed in backcrosses (Table 3). However, inbreeding is also expected to affect F1 × F1 crosses in the field, given the scarcity of hybridization events. It is also worth noting that a few backcrossed tadpoles survived metamorphosis. Among such rare individuals, females are most likely responsible for the mitochondrial introgression and signs of nuclear allele leakage observed in the field.
Important caveats obviously apply to our crossing experiments, mainly due to the practical difficulties in obtaining reproducing individuals from the field. Absence of replicates limits the power of our inferences, and we lack the exact controls for intraspecific matings (although the latter is compensated by our long experience of breeding green toads in the lab, which allowed us to provide intraspecific controls from other Western Palearctic lineages; Table 3). Results from these preliminary crosses must be clearly considered as provisional, but are worth reporting here as a support for our main conclusions gathered from field population-genetics data.
Though many hybrid zones have been documented in amphibians (see Background), few provide the data required to calibrate the speciation process, in terms of reliable time of divergence and patterns of introgression in the field. A notable exception is provided by fire- and yellow-bellied toads, which constitute one of the best-studied anuran hybrid zones. Bombina bombina and B. variegata, thought to have diverged during Upper Miocene or Lower Pliocene [3.5 Mya; [49–51]] hybridize in narrow, stable zones maintained by selection and dispersal . Strong selection against hybrids is generated both by hybrid breakdown [21, 52], and by environment-dependent selection against toads in mismatched habitats .
However, mitochondrial introgression in Bombina seems more limited than in Sicilian green toads, with mtDNA clines similar to or even steeper than those of nuclear loci [allozymes; ]. The divergence between B. siculus and B. balearicus is slightly more recent (Plio-Pleistocene, ca 2.7 Mya), but, despite higher mitochondrial DNA introgression, we found almost no admixture at supposedly neutrally evolving nuclear microsatellite loci, suggesting stronger selection to keep gene pools apart.
While anuran species that diverged > 8 Mya may exhibit partial or complete hybrid inviability in the laboratory, as recently shown for instance by a combination of experimental crosses and molecular divergence time estimates for the Fejevarya limnocharis group , there is accumulating evidence that anurans with Plio-(Pleisto)cene divergence tend to be reproductively isolated under natural conditions. Secondary contacts between Australian hylids of Plio-Pleistocene divergence have provided evidence for allopatric speciation driven by reinforcement mechanisms, characterized by highly asymmetric F1-viability in experimental crosses . Neither F1 and F2-hybrids nor backcrosses could be identified at a contact zone between two Hyla species with an estimated Pliocene divergence (3.6 Mya) , indicating a lack of current gene exchange. By contrast, phylogeographic studies [e.g. [16, 17]] have shown that clades of more recent divergence [1.33 Mya; ] form "wide hybrid zone(s) with a considerable genetic exchange between the two gene pools" , suggesting that reproductive barriers are still low or inexistent.
The observed post-mating barriers (hybrid breakdown) in Sicilian green toads certainly induce selection to act against hybridization. Therefore, we expect some pre-mating mechanisms to have evolved since the first contact of these allopatrically evolved lineages, which might also explain the absence of F1-hybrids. These two species have already been shown to differ in breeding phenology . Pre-mating isolation between closely related species (through mating calls and female choice) is widespread in anurans, sometimes even in absence of post-mating isolation . The nature of barriers might also affect the structure of hybrid zones, since, due to female choice, mtDNA is expected to introgress more readily than nuclear DNA under pre-zygotic isolation mechanisms (and even more so when the invading species is relatively rare compared to residents) . Further insights on the evolutionary interactions between B. siculus and B. balearicus populations might be gained by collecting bio-acoustic and mate-choice data, together with additional crossing experiments.
Whatever the exact causes, our data clearly show a virtual absence of gene flow at the present contact zone (corroborated with a hybrid breakdown in backcrosses and F2), meaning that the speciation process can be considered as close to have reached complete reproductive isolation after some 2.7 My of divergence. These field data contrast sharply with the results from experimental hybridization in anurans, which show that some lineages may still produce viable F1 offspring after ca. 20 My of divergence (see Background). Information on F1-hybridizability or viability gained under laboratory conditions may thus grossly overestimate the time required for genetic isolation and speciation to occur in anurans.
Samples from 323 specimens of green toads were collected during intense fieldwork between 2004 and 2007 at 24 localities across the Italian Peninsula and Sicily (Figure 1 and Table 1). Italian Peninsula was added to the sampling to allow us to understand the context of the B. balearicus invasion in Sicily and to compare Sicilian B. balearicus genotypes with those of mainland origin. A few samples came from scientific collections (MVZ, NME, ZFMK) or were collected throughout the years (e.g. road-kills). Tissue samples consisted of finger tips and muscles (road kill) from sub adult and adult toads, and tail tips from tadpoles. Most adults were released at the sampling sites, some vouchers were deposited in institutional collections [MVZ, NME, ZFMK, details Appendix 1 in 23], and tissues were stored in 98% ethanol and/or at -20°C.
Toad crosses were performed with mature adult males and females in a naturally reproductive state during the breeding period in spring. Females were stimulated to spawn by injection of 0.1 ml of a 0.9% NaCl-solution containing 500-1000 IU of human choriogonadotropin (Sigma).
A first crossing experiment was made in the laboratory in 2007 between a female B. balearicus (Si 41, pop. 11) and a male B. siculus (Si 11, pop. 22), from which 200 F1-hybrid tadpoles were raised. Ten randomly chosen offspring were genotyped and added to our sampling ("population 25"). Seven others crossing experiments were made in the laboratory in 2009 using two B. balearicus individuals (male: Si 336; female: Si337; Sicily, Marina S. Biagio, pop. 9), two B. siculus individuals (male: Si 334; female: Si 335; Sicily, Pergusa Lake, province of Enna, 37°31'00.29" N, 14°18'04.34"E, W of pop. 18 and 19) and two F1-hybrid individuals coming from the first crossing experiment (male: Cross 13; female: Cross 11). As no proper control for intraspecific crosses could be performed (due to lack of available females), we provide results from intraspecific matings of other green toad lineages (Table 3: Control 1: B. turanensis, Control2: B. pewzowi), which were raised under identical conditions.
Crossing pairs from 2009 were observed during amplexus every hour. When 10-15 cm of clutch strings had been deposited, couples were removed from the tank and separated, animals were rinsed to avoid sperm contamination and arranged in new cross combinations in another tank (see Table 3 for cross combinations).
Clutch was left untouched until hatching or until visible signs of dying eggs/embryos were found that had to be removed to avoid a chain-reaction of embryo suffocation. After one month, 100 tadpoles were randomly chosen among the surviving (if so) offspring, raised in shallow oxygenated aquaria and fed with Elodea plants and fish food (Tetramin) under identical conditions. After the second month of development, percentage of surviving tadpoles was determined by individual counts. All of these tadpoles were further raised until metamorphosis or death. Toadlets were fed with Drosophila and juvenile crickets, and mealworms with weekly additions of calcium and vitamins. Illumination of terrariums included sun-light spectrum fluorescent lamps including a natural-like UV fraction. Clutch, embryonic, tadpole and juvenile development were documented photographically (Figures 5 to 7).
DNA extraction and microsatellite data generation
DNA was extracted using the Qiagen DNeasy™kit. Microsatellite primer (Bcal μ10) described by Rowe et al.  and six pairs of microsatellites (BaturaC203, C205, C218, C223, D105, D5) developed by Colliard et al.  were selected for analysis, due to their level of polymorphism and applicability in both species. PCR-amplification, electrophoresis and allele scoring were performed as described in . Multiplexing of PCR products was performed. For D105 and D5 we used a mixture in a 1:3 ratio; Bcal μ10 was 5% diluted and mixed with C223 and C203 products in a 1:2:2 ratio. The third multiplexing involved products of C218 and C205 in a ratio of 1:2.
Mitotyping, sequencing and analyses of mitochondrial DNA and nuclear DNA
We used a mitotyping approach described in Colliard et al.  to detect the mtDNA haplotype group for all individuals. We also sequenced 577 bp of the mitochondrial control region (D-loop) [as described in ] in 162 individuals from across the study area (Figure 2).
An intron of alpha-Tropomyosine, situated between exons 5 and 6 was amplified, cloned and sequenced as described by Stöck et al. . We used homologous sequences from 21 individuals from throughout the sampling area plus one individual from a museum collection coming from a locality near pop. 12. Four of these 22 green toads came from the contact zone with assignment values varying from 68% to 92% to either parental species.
All mitochondrial and nuclear sequences were submitted to GenBank (Acc.-Numbers HM852594-HM852744, in part from ). Maximum likelihood (ML) phylogenies of mitochondrial and nuclear sequence alignment were generated using PhyML version 2.4.5  and using HKY models according to MrModeltest . In each case, we choose a BioNJ as a starting tree, and optimized topology, branch length and rate parameters. Other parameters were used as defaults of the program. We generated bootstrap values based on 1000 resampled data sets.
Genotype data analyses
To exclude genotyping errors due to null alleles, stuttering and allelic drop-out, Micro-Checker v.2.2.3  was used. Genotypic linkage disequilibrium between each pair of loci per population was tested using ARLEQUIN v. 3.0 . Significance was adjusted for multiple tests using Bonferroni corrections. Estimation of heterozygote deficits and its significance was assessed using FSTAT v. 2.9.4 . For each microsatellite marker and species, we estimated genetic diversity by calculating number of alleles (Na), and within-sample expected heterozygosity (Hs) (Additional file 1).
Isolation by distance was investigated with Mantel tests (FSTAT) by regressing pairwise FST/(1-FST) against the natural logarithm-transformed Euclidean geographic distances , for B. balearicus and B. siculus independently, selecting localities with ten or more samples. Two sets of analyses were performed, including or not the three populations at the contact zone which showed mitochondrial introgression. A principal component analysis was conducted with PCAGEN  to visualize pairwise differentiation among populations (FST), with 1000 randomizations of genotypes to test for significance of axes.
In order to better determine potential signs of introgression between the two species, we used the Bayesian clustering algorithm of STRUCTURE v. 2.2 . We used the admixture model and allowed for correlated allele frequencies between populations, as recommended by the authors for cases of subtle population structure. We tested a range of cluster numbers (K) from 1 to the number of localities per analysis, plus an additional three to enable us to potentially infer subtle structure. Each run, replicated 10 times, consisted of 105 iterations, after a "burn-in" of 104. To infer which K best fits our data, we applied the ad hoc ΔK statistic developed by Evanno et al. . We performed three analyses: (1) all Sicilian localities plus the lab-cross individuals; (2) all B. siculus localities, and (3) all B. balearicus localities.
Identification of hybrids
Four alternative approaches were used to identify hybrids. First, we used the cluster assignment value (STRUCTURE, parameters as above, K = 2), considering as hybrids all individuals with assignment values below 90%. Second, we performed two analyses using NEWHYBRIDS v.1.1 Beta3  to assign Sicilian toads to genotypic classes (parental, F1, F2, backcrosses). The method computes Bayesian posterior probability that an individual belongs to each of these different hybrid classes while simultaneously estimating allelic frequencies for parental species. Runs were repeated several times with varying lengths of the "burn-in" and number of sweeps, as recommended in the program manual. The first analysis was based on all Sicilian individuals. The second focused on populations where hybrids were shown to occur or were likely to do so (pop. 12 to 16 and 18), with addition of few pre-assigned non-hybrid individuals from "pure" populations (pop. 9 to 11 and 17), using the z option as recommended by the authors . Third, we identified individuals as hybrids if they were assigned by STRUCTURE to one clade but contained mtDNA from the other (= cyto-nuclear disequilibrium). Finally, we compared the microsatellite allele composition of Sicilian B. siculus and B. balearicus populations far apart from the contact zone (i.e., presumably pure) with those closer to the contact zone (potentially admixed). Alleles found in presumably pure populations that were exclusively present in one species were considered as "diagnostic". Finding diagnostic alleles of one clade in individuals assigned to the other clade by STRUCTURE was considered a sign of past hybridization (Additional file 3).
Research has been carried out according to approved guidelines under the following permits: Regione Siciliana, Assessorato Agricoltura e Foreste Prot. No 89884, Palermo, Italy; Bundesamt für Veterinärwesen BVET Nr. 791/09, Bern, Switzerland; and Authorisation No. 1798, Service de la consommation et des affaires vétérinaires, Canton de Vaud, Epalinges, Switzerland.
million years: Mya: million years ago
Museum of Vertebrate Zoology, University of California, Berkeley, USA
Naturkundemuseum Erfurt, Germany
Zoologisches Forschungsinstitut und Museum Alexander Koenig, Bonn, Germany.
Coyne JA, Orr HA: Speciation. 2004, Sunderland, Sinauer, MA
Johns GC, Avise JC: A comparative summary of genetic distances in the vertebrates from the mitochondrial cytochrome b gene. Mol Biol Evol. 1998, 15 (11): 1481-1490.
Sasa MM, Chippindale PT, Johnson NA: Patterns of postzygotic isolation in frogs. Evolution. 1998, 52 (6): 1811-1820. 10.2307/2411351.
Wilson AC, Maxson LR, Starich VM: Two types of molecular evolution. Evidence from studies of interspecific hybridization. PNAS. 1974, 71: 2843-2847. 10.1073/pnas.71.7.2843.
Blair WF: Evolution in the genus Bufo. 1972, Austin and London
Malone JH, Fontenot BE: Patterns of reproductive isolation in toads. PLoS ONE. 2008, 12 (3): e3900-10.1371/journal.pone.0003900.
Noor MAF, Feder JL: Speciation genetics: evolving approaches. Nature. 2006, 7: 851-861.
Hoskin CJ, Higgie M, McDonald KR, Moritz C: Reinforcement drives rapid allopatric speciation. Nature. 2005, 437: 1353-1356. 10.1038/nature04004.
Mcdonnell LJ, Gartside DF, Littlejohn MJ: Analysis of a narrow hybrid zone between 2 species of Pseudophryne (Anura-Leptodactylidae) in Southeastern Australia. Evolution. 1978, 32 (3): 602-612. 10.2307/2407726.
Parris MJ: High larval performance of leopard frog hybrids: Effects of environment-dependent selection. Ecology. 2001, 82 (11): 3001-3009.
Gartside DF, Littlejohn MJ, Watson GF: Structure and dynamics of a narrow hybrid zone between Geocrinia laevis and Geocrinia victoriana (Anura, Leptodactylidae) in Southeastern Australia. Heredity. 1979, 43 (Oct): 165-10.1038/hdy.1979.72.
Littlejohn MJ, Watson GF: Hybrid zones and homogamy in Australian frogs. Ann Rev Ecol Syst. 1985, 16: 85-112. 10.1146/annurev.es.16.110185.000505.
Lamb T, Avise JC: Directional introgression of mitochondrial-DNA in a hybrid population of tree frogs - the influence of mating-behavior. PNAS. 1986, 83 (8): 2526-2530. 10.1073/pnas.83.8.2526.
Pfenning KS: Facultative mate choice drives adaptive hybridization. Science. 2007, 318: 965-967. 10.1126/science.1146035.
Yamazaki Y, Kouketsu S, Fukuda T, Araki Y, Nambu H: Natural hybridization and directional introgression of two species of Japanese toads Bufo japonicus formosus and Bufo torrenticola (Anura: Bufonidae) resulting from changes in their spawning habitat. J Herpetol. 2008, 42 (3): 427-436. 10.1670/07-186.1.
Santucci F, Nascetti G, Bullini L: Hybrid zones between two genetically differentiated forms of the pond frog Rana lessonae in southern Italy. J Evol Biol. 1996, 9 (4): 429-450. 10.1046/j.1420-9101.1996.9040429.x.
Canestrelli D, Nascetti G: Phylogeography of the pool frog Rana (Pelophylax) lessonae in the Italian peninsula and Sicily: multiple refugia, glacial expansions and nuclear-mitochondrial discordance. J Biogeogr. 2008, 35 (10): 1923-1936. 10.1111/j.1365-2699.2008.01946.x.
Stöck M, Moritz C, Hickerson M, Frynta D, Dujsebayeva T, Eremchenko V, Robert Macey J, Papenfuss TJ, Wake DB: Evolution of mitochondrial relationships and biogeography of Palearctic green toads (Bufo viridis subgroup) with insights in their genomic plasticity. Mol Phylogenet Evol. 2006, 41: 663-689. 10.1016/j.ympev.2006.05.026.
Knoop T, Merilä J: The postglacial recolonization of Northern Europe by Rana arvalis as revealed by microsatellite and mitochondrial DNA analyses. Heredity. 2009, 102: 174-181. 10.1038/hdy.2008.91.
Verardi A, Canestrelli D, Nascetti G: Nuclear and mitochondrial patterns of introgression between the parapatric European treefrogs Hyla arborea and H. intermedia. Ann Zool Fenn. 2009, 46 (4): 247-258.
Szymura JM, Barton NH: Genetic analysis of a hybrid zone between the fire-bellied toads, Bombina bombina and B. variegata, near Cracow in Southern Poland. Evolution. 1986, 40 (6): 1141-1159. 10.2307/2408943.
Yanachukov A, Hofman S, Szymura JM, Mezhzherin SV, Morozov-Leonov SY, Barton NH, Nürnberger B: Hybridization of Bombina bombina and B. variegata (Anura, Discoglossidae) at a sharp ecotone in western Ukraine: Comparisons across transects and over time. Evolution. 2006, 60 (3): 583-600.
Stöck M, Sicilia A, Belfiore NM, Buckley D, Lo Brutto S, Lo Valvo M, Arculeo M: Post-Messinian evolutionary relationships across the Sicilian channel: Mitochondrial and nuclear markers link a new green toad from Sicily to African relatives. BMC Evol Biol. 2008, 8: 56-10.1186/1471-2148-8-56.
Colliard C, Sicilia A, Moritz C, Perrin N, Stöck M: New polymorphic microsatellite markers and development of mitotyping primers for West Mediterranean green toad species (Bufo viridis subgroup). Mol Ecol Resour. 2009, 9: 1138-1140. 10.1111/j.1755-0998.2009.02600.x.
Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data. Genetics. 2000, 155: 945-959.
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-2620. 10.1111/j.1365-294X.2005.02553.x.
Anderson EC, Thompson EA: A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002, 160: 1217-1229.
Goudet J: 1999, [http://www2.unil.ch/popgen/softwares/pcagen.htm]
Canestrelli D, Cimmaruta R, Nascetti G: Phylogeography and historical demography of the Italian treefrog, Hyla intermedia, reveals multiple refugia, population expansions and secondary contacts within peninsular Italy. Mol Ecol. 2007, 16: 4808-4821. 10.1111/j.1365-294X.2007.03534.x.
Bonfiglio L, Manganoa G, Marraa AC, Masinib F, Paviac M, Petrusob D: Pleistocene Calabrian and Sicilian bioprovinces. Geobios. 2002, 35: 29-39. 10.1016/S0016-6995(02)00046-3.
Avise JC, Neigel JE, Arnold J: Demographic influences on mitochondrial-DNA lineage survivorship in animal populations. J Mol Evol. 1984, 20 (2): 99-105. 10.1007/BF02257369.
Sullivan BK: Mating system variation in Woodhouse's toad (Bufo woodhousii). Ethology. 1989, 83 (1): 60-68. 10.1111/j.1439-0310.1989.tb00519.x.
Sullivan BK, Ryan MJ, Verrell PA: Female choice and mating system structure. Amphibian biology. Edited by: Heatwole H, Sullivan BK. 2005, Baulkham Hills, NSW, Australia: Sons SB. Baulkham Hills, 469-517.
Canestrelli D, Verardi A, Nascetti G: Genetic differentiation and history of populations of the Italian treefrog Hyla intermedia: lack of concordance between mitochondrial and nuclear markers. Genetica. 2007, 130 (3): 241-255-10.1007/s10709-006-9102-9.
Lampert KP, Rand AS, Mueller UG, Ryan MJ: Fine-scale genetic pattern and evidence for sex-biased dispersal in the tungara frog, Physalaemus pustulosus. Mol Ecol. 2003, 12 (12): 3325-3334. 10.1046/j.1365-294X.2003.02016.x.
Monsen KJ, Blouin MS: Genetic structure in a montane ranid frog: restricted gene flow and nuclear-mitochondrial discordance. Mol Ecol. 2003, 12 (12): 3275-3286. 10.1046/j.1365-294X.2003.02001.x.
Currat M, Ruedi M, Petit RJ, Excoffier L: The hidden side of invasions: Massive introgression by local genes. Evolution. 2008, 62 (8): 1908-1920.
Petit RJ, Excoffier L: Gene flow and species delimitation. Trends Ecol Evol. 2009, 24 (5): 386-393. 10.1016/j.tree.2009.02.011.
Barton NH, Gale KS: Genetic-analysis of hybrid zones. Hybrid Zones and the Evolutionary Process. Edited by: Harrison R. 1993, New York: Oxford University Press, 13-45.
Flindt R, H H: Nachweis einer natürlichen Bastardisierung von Bufo calamita und Bufo viridis. Zool Anz. 1967, 178: 419-429.
Born G: Biologische Untersuchungen II. Weitere Beiträge zur Bastardisierung zwischen den einheimischen Anuren. Archive für Mikroskopie und Anatomie. 1888, 27: 192-271. 10.1007/BF02955420.
Zavadil V, Roth P: Natural hybridization between Bufo viridis and Bufo bufo in the Doupovské hory hills (northeast Bohemia, Czech Republic) with general comments on hybridization of european green toads and common toads. (Hrsg): Herpetologia Bonnensis. Edited by: Böhme, W, Bischoff W, Ziegler T. 1997, 401-404.
Drago A: Atlante climatologico della Sicilia. 2002, Assessorato Agricoltura e Foreste, Regione Sicilia (CD-ROM), 2
Turelli M, Orr HA: The dominance theory of Haldanes rule. Genetics. 1995, 140 (1): 389-402.
Haldane JBS: Sex ratio and unisexual sterility in hybrid animals. J Genet. 1922, 12 (2): 101-109. 10.1007/BF02983075.
Kawamura T, Nishioka M, Ueda H: Inter- and intraspecific hybrids among Japanese, European and American toads. Sci Rep Lab Amphib Biol Hiroshima Univ. 1980, 4: 1-125.
Ueda H: Offspring of sex-reversed males in Bufo viridis. Sci Rep Lab Amphib Biol, Hiroshima Univ. 1990, 10: 155-164.
Odierna G, Aprea G, Capriglione T, Castellano S, Balletto E: Cytological evidence for population-specific sex chromosome heteromorphism in Palaearctic green toads (Amphibia, Anura). J Biosciences. 2007, 32 (4): 763-768. 10.1007/s12038-007-0076-2.
Vines TH, Köhler SC, Thiel M, Ghira I, Sands TR, MacCallum CJ, Barton NH: The maintenance of reproductive isolation in a mosaic hybrid zone between the fire-bellied toads Bombina bombina and B. variegata. Evolution. 2003, 57 (8): 1876-1888.
Hofman S, Spolsky C, Uzzell T, Cogalniceanu D, Babik W, Szymura JM: Phylogeography of the fire-bellied toads Bombina: independent Pleistocene histories inferred from mitochondrial genomes. Mol Ecol. 2007, 16 (11): 2301-2316. 10.1111/j.1365-294X.2007.03309.x.
Szymura JM: Analysis of hybrid hones with Bombina. Hybrid Zones and the Evolutionary Process. Edited by: Harrison, R. 1993, 261-289.
Kruuk LEB, Gilchrist JS, Barton NH: Hybrid dysfunction in fire-bellied toads (Bombina). Evolution. 1999, 53 (5): 1611-1616. 10.2307/2640907.
MacCallum CJ, Nürnberger B, Barton NH, Szymura JM: Habitat preference in the Bombina hybrid zone in Croatia. Evolution. 1998, 52 (1): 227-239. 10.2307/2410938.
Hofman S, Szymura JM: Limited mitochondrial DNA introgression in a Bombina hybrid zone. Biol J Linn Soc. 2007, 91 (2): 295-306. 10.1111/j.1095-8312.2007.00795.x.
Sumida M, Kotaki M, Islam MM, Djong TH, Igawa T, Kondo Y, Matsui M, De Silva A, Khonsue W, Nishioka M: Evolutionary relationships and reproductive isolating mechanisms in the Rice Frog (Fejervarya limnocharis) species complex from Sri Lanka, Thailand, Taiwan and Japan, inferred from mtDNA gene sequences, allozymes, and crossing experiments. Zool Sci. 2007, 24: 547-562. 10.2108/zsj.24.547.
Blair WF: Isolating mechanisms and interspecific interactions in anuran amphibians. Q Rev Biol. 1964, 39: 334-344. 10.1086/404324.
Chan KMA, Levin SA: Leaky prezygotic isolation and porous genome: rapid introgression of maternally inherited DNA. Evolution. 2005, 59: 720-729.
Rowe G, Beebee T, Burke T: A further four polymorphic microsatellite loci in the natterjack toad Bufo calamita. Conserv Genet. 2000, 1 (4): 371-373. 10.1023/A:1011567502714.
Guindon S, Gascuel O: A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003, 52: 696-704. 10.1080/10635150390235520.
Nylander JAA: MrModeltest v2. Edited by: Pdbt. 2004, Evolutionary Biology Centre, Uppsala University
Van Oosterhout C, Hutchinson WF, Derek P, Wills M, Shipley P: MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004, 4 (3): 535-538. 10.1111/j.1471-8286.2004.00684.x.
Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol Bioinform. 2005, 1: 47-50.
Goudet J: FSTAT (Version1.2): A computer program to calculate F-statistics. J Hered. 1995, 86 (6): 485-486.
Rousset F: Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics. 1997, 145 (4): 1219-1228.
We sincerely thank M. Lo Valvo, F. Lillo, and R. Termine for contributing samples, and K. Ghali for help in the lab. We further thank T. Broquet, J. Goudet and A. Horn for advice and references, and J. van Rooyen, M. Vences, and two anonymous reviewers for their helpful comments on the manuscript.
Funding: Swiss NSF grant 31003A_129894 to NP.
The fieldwork was conducted by AS, CC, GFT and MSt, and the molecular work (genotyping, mitotyping and sequencing) by AS, CC and MSt. The crossing experiments were performed by CC and MSt, and the statistical analyses and data interpretation by CC, MSt, MA and NP. CC drafted the paper, which was edited by MSt and NP, then completed by AS, GFT and MA. All authors read and approved the final manuscript.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.