Tadpoles of hybridising fire-bellied toads (B. bombina and B. variegata) differ in their susceptibility to predation

Abstract

The role of adaptive divergence in the formation of new species has been the subject of much recent debate. The most direct evidence comes from traits that can be shown to have diverged under natural selection and that now contribute to reproductive isolation. Here, we investigate differential adaptation of two fire-bellied toads (Anura, Bombinatoridae) to two types of aquatic habitat. Bombina bombina and B. variegata are two anciently diverged taxa that now reproduce in predator-rich ponds and ephemeral aquatic sites, respectively. Nevertheless, they hybridise extensively wherever their distribution ranges adjoin. We show in laboratory experiments that, as expected, B. variegata tadpoles are at relatively greater risk of predation from dragonfly larvae, even when they display a predator-induced phenotype. These tadpoles spent relatively more time swimming and so prompted more attacks from the visually hunting predators. We argue in the discussion that genomic regions linked to high activity in B. variegata should be barred from introgression into the B. bombina gene pool and thus contribute to gene flow barriers that keep the two taxa from merging into one.

About the Authors:

Radovan Smolinský

Roles Data curation, Investigation, Methodology, Project administration, Writing – review & editing

Current address: Department of Biology, Faculty of Education, Masaryk University, Brno, Czech Republic

Affiliation: Research Facility Studenec, Institute of Vertebrate Biology, Czech Academy of Sciences, Brno, Czech Republic

Vojtech Baláž

Roles Funding acquisition, Investigation, Methodology, Resources, Writing – review & editing

Affiliation: Faculty of Veterinary Hygiene and Ecology, University of Veterinary and Pharmaceutical Sciences, Brno, Czech Republic

Beate Nürnberger

Roles Conceptualization, Formal analysis, Funding acquisition, Methodology, Software, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Research Facility Studenec, Institute of Vertebrate Biology, Czech Academy of Sciences, Brno, Czech Republic

ORCID logo https://orcid.org/0000-0002-6133-830X

Introduction

Adaptation to local environments is commonplace in nature [1] and drives the evolution of novel ecotypes. Whether adaptive phenotypic divergence plays an important role in the origin of new species has been the subject of much discussion [2–7]. In this case, adaptation to different ecological conditions should entail a marked fitness reduction for ecotypes outside the habitat to which they are adapted and it should contribute to reproductive isolation [7]. Of particular interest are case studies in which the selective agents and the traits under selection have been identified [7], because then the mechanisms that link adaptation to the emergence of reproductive barriers can be studied and their importance assessed. For example, some populations of Rhagoletis pomonella fruit flies have recently switched from hawthorn to apple hosts and have adapted their phenology and fruit odour preference accordingly [8, 9]. In stickleback fish (Gasterosteus aculeatus), the repeated evolution of distinct trophic ecotypes within lakes is mirrored in the intricate, adaptive divergence of cranial morphology [10]. The intertidal snail Littorina saxatilis produces shells with different tickness and aperture in wave-exposed versus crab-rich sections of the shoreline, which provide better protection against dislodgement [11] and predation [12], respectively. These adaptations contribute to reproductive isolation, because each ecotype is less fit when transplanted into the respective other habitat [13, 14], because F1 hybrids are at a disadvantage [15, 16] and/or because the adaptations cause assortative mating [17–19]. A genomic dissection of such divergent traits can not only pinpoint the genomic regions under selection but also investigate genetic linkage between ecologically selected and mating traits [19] or, for example, explore the role of inversions in novel adaptations [20, 21]. Ecological divergence has long been assumed to generate barriers to gene flow in the hybrid zone of European fire-bellied toads [22], for which genomic resources have recently been developed [23]. Here we seek experimental evidence for adaptive divergence in tadpole traits as a first step in establishing a link between ecology and reproductive isolation.

In contrast to the very young taxa referenced above, the red-bellied toad (Bombina bombina) and the yellow-bellied toad (B. variegata) are more anciently derived [MRCA ~ 3.2 mya, 23]. Nevertheless, they interbreed wherever they meet and produce a wide range of recombinants in typically narrow (2–7 km wide) hybrid zones [24]. Classic cline theory [25] suggests that the mean fitness of intermediate hybrid populations is reduced by 42% compared to that of pure populations on either side [26]. There is experimental evidence for endogenous selection against hybrid embryos and tadpoles [27]. Profound ecological divergence is also strongly implicated. B. bombina reproduces in semi-permanent lowland ponds, such as flooded meadows and oxbows that fall dry in winter [e.g. 28]. B. variegata lays eggs in more ephemeral habitat (‘puddles’) and sometimes also in small steams [29], typically at higher elevations. Consequently, hybrid zones tend to be located at ecotones. There is a long list of taxon differences that presumably adapt both taxa to their current habitat [reviewed in 22]. For example, adult B. variegata have a relatively stronger skeleton, longer legs [30] and thicker skin [31], traits which presumably aid their frequent movements over land [32, 33], to track the ever-shifting mosaic of suitable breeding sites. Ponds are not only more predictable year-on-year but also allow for larger breeding aggregations. With their well-developed vocal sacs [34], male B. bombina produce much louder mating calls [35] and form choruses that can be heard from afar. A B. variegata male would likely be at a disadvantage if it joined such a chorus.

Here we focus on tadpoles, because they experience the strongest habitat contrast. Faced with the ever present risk of desiccation, B. variegata tadpoles hatch from relatively larger eggs, grow faster and metamorphose earlier [30, 36]. They also spend relatively more time swimming and feeding [37, 38]. The less active behaviour of B. bombina tadpoles should protect them against visually hunting predators in ponds [37]. These taxon differences reflect a general pattern across anurans. More permanent aquatic habitats have a greater abundance and diversity of predators [39–41]. Anuran species that reproduce in them have typically quiescent tadpoles [40, 42–44] that take longer to reach metamorphosis than those that breed in ephemeral sites [42, 43, 45]. More active tadpoles species are typically more susceptible to predation [42, 46–50]. This robust pattern has been interpreted as a trade-off [42, 45, 51, 52]: tadpoles can either be highly active and maximise food intake for fast growth and development in temporary habitats or avoid predators through quiescent behaviour in more permanent sites, but not both. Tadpole activity and growth rate covary as expected in evolutionary contrasts [53] but evidence for a causal relationship is mixed at the intra-specific level [54–59]. However, the trade-off between rapid growth and predator avoidance is almost universal [60] irrespective of whether activity is the mechanism that links the two. We therefore expect that B. variegata tadpoles are at greater risk of predation than B. bombina.

Nevertheless, B. variegata tadpoles do sporadically encounter predators such as newts, dragonfly and dytiscid larvae, water scorpions and salamander larvae [29]. Thus, tadpoles would maximise their chance of reaching metamorphosis by adjusting their phenotype to the local balance of risk. Indeed, a wide range of anurans display phenotypic plasticity in response to predators. Almost invariably, tadpoles reduce their activity level when they perceive chemical cues that indicate predator presence [45, 61, 62]. This response can occur within minutes [63]. Longer term exposure to such cues can cause an increase in refuge use [64, 65], alter larval morphology [62], prompt the development of colour spots on the tail [66, 67], and slow growth and development [68, 69]. These effects are seen across venues from the laboratory to the field [70, 71].

The commonest and most closely analysed morphological response to predators is the development of a relatively deeper tail fin [43, 44, 55, 62], usually in combination with a shorter tail and a shallower body. A deeper tail may facilitate rapid escape (faster burst swimming speed) or lure attacks away from the more vulnerable body. While there is evidence that natural variation in tail shape affects tadpole swimming performance [50, 72], studies to date have given stronger support to the lure hypothesis [73, 74]. Irrespective of the particular mechanism, there is evidence that predator-induced phenotypes are eaten less than predator-naïve ones [75–78] and that survival is correlated with measured traits that show a plastic response [78–80]. A recent study of Rana temporaria [81] makes a particularly strong case for adaptive phenotypic plasticity in both behaviour and morphology, including an assessment of the variability in predation risk within and between ponds and over time that could bring about such a flexible strategy.

Both B. bombina and B. variegata respond to the presence of predators with reduced activity [37, 38] and by developing a relatively deeper tail fin [82, see also 62]. Exposure to chemical predation cues slows their growth and development under laboratory conditions [82]. While there is at present no evidence that the induced phenotypes are better protected from predators, they should be used for a fair comparison of relative predation risk, because tadoles would develop with predators in nature. A previous study using predator-naïve tadpoles, while showing the expected result, could therefore not fully settle the issue [37]. In a second laboratory experiment with predator-induced tadpoles, a strong population effect overwhelmed any taxon differences [82]. In both cases, significant effects were observed only in the earlier of two sets of predation trials, which suggests that relative susceptibility may change over time.

Here, we carried out laboratory experiments that improve on the design of the previous two studies by using predator-induced phenotypes throughout, by increasing the level of replication within each taxon and by estimating relative predation risk over a four-week period, from hatching to the time when the first B. variegata tadpoles start to metamorphose. Dragonfly larvae were used as predators, because they are ubiquitous in ponds and also found in puddles. We also assessed the contributions of activity, body size and tail shape to the observed differences in predation risk. Our experiments confirmed our expectation that pure B. variegata suffer a higher predation rate than pure B. bombina due to a higher activity level. In the discussion, we argue that elevated activity contributes to the gene flow barrier across the hybrid zone which keeps the two gene pools from merging into one.

Materials and methods

Animal collections, pairing scheme and toad care

Adult toads were collected from Moravia (Czech Republic) where the Carpathian lineage of B. variegata meets the Southern lineage of B. bombina [83]. The experiment was carried out in 2017. Starting in early May, aquatic habitats were visited regularly (14 known B. bombina sites and 19 known B. variegata sites). Due to an unusually cold spring, toads did not appear in appreciable numbers until early June and were then promptly collected to ensure that females had not yet laid eggs. Collection sites (Table 1) were located at least 6 km from known areas of hybridisation and deemed 'pure' based on the taxon-specific spot pattern on the toads' bellies [22]. For each taxon, minimum and maximum distances between sites were 1.5 and 43.6 km (B. bombina) and 0.7 and 14.9 km (B. variegata). B. bombina sites were ponds with a minimum width of 5 m and a maximum depth greater than 1 m. They all contained predators (incl. Aeshna cyanea, Natrix natrix, fish) and, with one exception (site E), they had abundant aquatic vegetation. B. variegata sites were either water-filled wheel ruts (max. depth < 30 cm) or in one case (site 2) the shallow end of an abandoned swimming pool. They were devoid of aquatic vegetation and three of nine sites contained no predators at all (see S1 Table for details).

[Figure omitted. See PDF.]

Table 1. Collection sites and pairing scheme.

https://doi.org/10.1371/journal.pone.0231804.t001

Random numbers were drawn to pair males and females within sites whenever possible. However due to unequal sex ratios per site, four B. variegata pairs involved toads from separate, but nearby sites (Table 1, maximum distance between sites per pair: 0.96 km). Each pair was housed in a 52 x 35 cm plastic box with 7L of well water, a polystyrene island and three pieces of plastic aquarium plants as support for egg batches. As a precaution against the potential spread of Batrachochytrium infection, all waste water from the entire experiment was heated to 100° C and treated with UV (Eheim UV reflex UV800) before disposal to kill any fungal spores.

On 20 June, a belly swab was taken from each toad to test for Batrachochytrium dendrobatidis (Bden) and B. salamandrivorens (Bsal) infection. The former pathogen is widely distributed in the Czech Republic, but so far no population declines have been observed [84]. The qPCR followed the protocol of Blooi et al. [85] using a Roche Light Cycler480. No evidence of either infection was found. We note that samples in the same qPCR run but from an unrelated study scored positive for Bden infection. The adults were subsequently returned to their exact collection sites.

Predators

We used Aeshna cyanea larvae as predators. This dragonfly is widely distributed across Europe and was present in most of the habitats we surveyed (see above). Two hundred larvae were collected from one site near Jihlava (Vysočina County) and were kept in groups of 15 in 50L plastic boxes. They were fed every other day with Tubifex worms before the predation trials started and with tadpoles thereafter.

Oviposition

To induce reproduction, 250 I.U. of human chorionic gonadotropin, dissolved in 0.2 ml of amphibian Ringer solution, were injected into the lymph sac of an outer thigh of both males and females on 6 June. Most pairs produced egg clutches overnight. Over the following two days, 30–40 eggs per family were transferred into the rearing set-up (see below). The remaining offspring were raised in the oviposition boxes, with weekly water changes. They were used as food for dragonflies in order to condition the water in the rearing set-up. The adults were moved to new boxes of the same size. Tadpoles hatched on 12–13 June. In the following, we refer to 13 June as day 0 of the experiment.

Tadpole rearing

We constructed two gravity flow circuits (Fig 1). For each of them, five plastic tubs (52 x 37 x 16 cm) were interconnected with silicone tubing such that water could be pumped from a 20 L bucket at floor level into the top tub and then descend stepwise via the other four tubs before returning to the bucket. Inflows and outflows in each tub were positioned diagonally opposite to maximise mixing. Each tub held 27L of water at a depth of 14 cm. In order to maintain uniform water quality, the pumps (Eheim StreamON+ 2000) were operated 24 hours per day at a low flow rate and per circuit 7 L of water were replaced daily. Throughout the experiment, water temperature was 22–24°C and the air temperature was 24–26°C. All water changes for adults and tadpoles were carried out with well water heated to 22°C and sterilised with UV (Eheim UV reflex UV800).

[Figure omitted. See PDF.]

Fig 1. Schematic of one gravity flow circuit.

At the start of the experiment, each of the 32 compartments contained 10 hatchlings from a single family. A dragonfly was regularly fed non-experimental tadpoles in the bottom tub (P) so that chemical predation cues were present throughout tadpole development. Blue lines indicate silicone tubing and arrows show the direction of flow.

https://doi.org/10.1371/journal.pone.0231804.g001

The upper four tubs in each circuit were subdivided with fibreglass insect mesh into eight rearing compartments (size: 18 x 11.5 x 16.5 cm, 64 compartments in total). A plastic aquarium plant was added to each compartment for structural complexity. At the start of the experiment, each compartment held 10 hatchlings from a single family. Families (nine per taxon) were distributed over either three or four compartments in total and either one or two compartments in each circuit. Exactly half of the compartments in each circuit were assigned to each taxon. Placement of families within circuits was decided by random draws. Tadpoles were fed boiled nettle (Urtica dioica) leaves ad libitum and uneaten food was removed every day.

In order to simulate predator presence, we let chemical cues of tadpole predation circulate from day 0 onwards. For this, we used the fifth tub in each circuit to feed 1–2 non-experimental tadpoles (depending on size) to one dragonfly larva every other day. Tadpoles of both taxa were used for this in equal proportions.

Predation trials

Predation trials started on day 2 when the tadpoles were active. They were carried out in two plastic boxes subdivided with fibreglass insect mesh into four sections (26 x 26 x 18 cm), each with two plastic aquarium plants (akin to Myriophyllum), placed into diagonally opposite corners. They contained a 1:1 mixture of well water and water from the rearing circuits. For each trial, four tadpoles (2 B. bombina, 2 B. variegata) were placed into one section and allowed to acclimatize for 30 min. A dragonfly that had not been fed on the previous day was then added and the trial was stopped when the first tadpole had been eaten (maximum duration: 45 min). The identity of the missing tadpole was determined from digital photographs of the animals taken before the trial (see below). The survivors were returned to their rearing compartments. The trials were recorded with a Microsoft LifeCam VX-2000 video camera suspended over each plastic box.

Predation trials were carried out on 18 days as follows: two consecutive experiment days were followed by one ‘off-day’ (1 exception); this pattern was repeated nine times until day 29 post hatching, when the first B. variegata tadpoles entered metamorphosis. There were 20 trials per day between 10:00 and 15:00 hours. A custom Perl script was used to haphazardly draw the composition of each trial per day such that (a) four different families were represented, (b) one tadpole per taxon came from circuit 1 and the other from circuit 2 and (c) each rearing compartment contributed at least one tadpole on a given day. However towards the end of the experiment, a few rearing compartments were empty and the contribution of others had to be increased accordingly. No individual tadpole or dragonfly was used more than once per experiment day.

Morphometrics

Digital images were taken by placing a tadpole into a glass cuvette and photographing it from the side with a Canon EOS 400D camera (lens: Canon EF 28-90mm). tpsDIG2 (F. James Rohlf) software was used to take the following measurements from these images: total length, tail length, body height and maximum tail height.

Behaviour

From the video footage of the predation trials, we recorded the behaviour of a focal tadpole during one 30 s interval per trial for half of the trials on a given day. A custom Perl script provided the pseudo-random start point of the observation interval and allowed behaviour (swimming, feeding, resting) to be recorded with separate key strokes. From one trial to the next, we alternated between a B. bombina and B. variegata tadpole. Taxon identification was based on size in the early trials. Morphological differences were used as soon as they became apparent. Footage from the first two days was not used, because the taxa could not be reliably distinguished. We discarded all observation intervals that overlapped with a predator attack. We could not determine whether the focal tadpole survived the trial, because that would have increased the data collection time roughly 10 to 20-fold. Moreover, it was difficult to track an individual over several minutes, especially when it temporarily disappeared together with others under a plastic plant. We also recorded for 68 predator attacks (four on each of 17 experiment days) whether they were initiated by and directed towards an actively swimming tadpole.

Permits

Collection permits were granted by the relevant regional offices (South Moravia and Zlín) of the Czech Ministry for the Environment. The experimental protocol was approved by the ethics committee of the Institute of Vertebrate Biology, CAS, permit no. 15/2016.

Statistical analyses

Estimates of predation risk.

We estimated predation risk per rearing compartment (= the unit of observation, containing offspring of one family) as the number of tadpoles eaten divided by the number of valid trials in which that compartment was represented. This figure subsumes the ontogeny of tadpoles in terms of changes in size, body shape, behaviour and, possibly, experience in predators avoidance. Recall that surviving tadpoles were returned to their rearing compartment and were likely picked again in later predation trials. Just as in nature and in the classic mesocosm studies [46, 55, 86], in which tadpole cohorts interacted with predators for a substantial part of larval development, older tadpoles were likely to have survived previous predator encounters. Specific to our design is that we staged these encounters and collected additional data about each one of them.

Especially with a small number of trials, our estimate of predation risk may be biased if tadpoles from a certain compartment happened to be assigned to trials with, say, particularly well defended congeners. There is a range of statistical methods to estimate an individual's skill from a limited number of contests with opponents of varying strengths (as e.g. in chess), starting with the model of Bradley and Terry [87]. However, none of these match our trial design of a single winner (or rather: looser) and a 'tie' for the other three.

We therefore checked for biases due to a limited number of trials as follows. Each trial k involves tadpoles from four (of 64) compartments. The risk context rcik of the tadpole from compartment i (1 ≤ i ≤ 4) is the mean of the predation risks rj of the other three compartments that are represented:

For a given compartment m (1 ≤ m ≤ 64), the average over its associated rcik estimates in valid trials is the mean risk context, rcm. If our estimates of predation risk, rm, are unbiased, then the estimates of risk context from the observed dataset should be very similar to those from a much larger set of trials. To test this, we haphazardly drew compartments for another 2550 trials with the same script as before and used these ‘mock trials’ to re-compute risk contexts using the rm estimates from the actual experiments.

We used linear mixed models to investigate the effect of taxon on predation risk and included ‘circuit’ and ‘family’ terms to properly represent the non-independence of compartments. Circuit was treated as a fixed effect, because it had only two levels. We did not include a population term in the model, because populations were unevenly represented and some families had parents from different sites. Any differences between populations should appear as family effects. For these, we fitted not only random intercepts but also allowed the effect of ‘circuit’ to vary among families (random slopes). Model were fitted with the R module lme4 [88, 89]. We used the ‘REML = FALSE’ option in order to compare models via a likelihood ratio test.

Morphological correlates of survival.

Logistic models were used to investigate the effects of taxon, tadpole size and tail shape on survival per trial. The input was generated by sorting all valid trials in chronological order and picking from them alternately (a) a surviving tadpole at random or (b) the non-surviving tadpole. For the entire dataset, the four morphological measurements were log-transformed to improve normality and entered into a principal component analysis (PCA). We treat the first principal component as a measure of tadpole size. From a regression of log(tail height) on PC1 we extracted the residuals as a measure of relative tail height (≈ tail shape). To ensure shared allometry between the taxa, we inspected pairwise plots of the four measured variables [90]. For the focal tadpole per trial, PC1 and residual log(tail height) were expressed as a deviation from the trial mean. All analyses were carried out in R v.3.4.4 [91].

Results

Taxon differences in predation risk

We assessed relative predation risk over a four-week interval. Our experiments started two days after hatching when tadpoles were uniformly at Gosner [92] stage 23. On the last experiment day, B. bombina tadpoles were mostly at stages 35–36 whereas the commenest B. variegata stages were 39–40 with some individuals entering metamorphosis (stages 42+). In total, 360 trials were conducted. Of these, no tadpole was eaten in 97 trials, more than one was eaten in five trials and in three trials the two taxa were not evenly represented. Discarding these, 255 valid predation trials remain. The proportion of valid trials in which a B. variegata tadpole was eaten was 0.65. The taxon difference was particularly strong for 5–12 day old tadpoles, when around 77% of the tadpoles eaten were B. variegata (Fig 2A). There was also a slightly higher risk for tadpoles from circuit 1 as opposed to circuit 2 (proportion = 0.56). This effect was largely restricted to B. bombina. The number of trials in which no tadpole was eaten was smaller in the first half of the experiment (tadpole age: 2–14 days, n = 31) than in the second half (tadpole age: 15–29 days, n = 66, Fig 2B). The mean durations of valid trials were 9.9 min (fist half) and 13.8 min (second half).

[Figure omitted. See PDF.]

Fig 2. Relative predation risk and number of valid trials per day.

A. The proportion of tadpoles consumed that were B. variegata. Tadpole age (days after hatching) is plotted on the x-axis. Each data point gives the proportion for one experiment day. B. The number of valid trials (out of 20) per experiment day.

https://doi.org/10.1371/journal.pone.0231804.g002

Predation risk was computed per rearing compartment as the number of tadpoles eaten during valid predation trials over the number valid trails in which this compartment was represented. The observed risk contexts (as defined in the Materials and Methods) closely matched those re-computed for a much larger number of hypothetical trials ( compartment = 159.9): for 56 of 64 compartments the observation deviated by less than 10% from the expectation (max. deviation 14%) and in all cases it was well within one standard deviation of the prediction (S1 Fig). We therefore use predation risk per rearing compartment as the response variable in the following analyses (Fig 3).

[Figure omitted. See PDF.]

Fig 3. Predation risk per rearing compartment arrayed by family.

B. variegata: families 1–4 and 6–10, B. bombina: families 11–13 and 15–20.

https://doi.org/10.1371/journal.pone.0231804.g003

The full linear mixed model (Table 2, model 1) included taxon, circuit and family (random effect, nested within taxon). We allowed ‘circuit’ to affect each family differently (random slopes). As expected, there was an increased risk for B. variegata (β = 0.2069 ± 0.0357) and a slightly decreased in risk in circuit 2 (β = -0.0814 ± 0.0376). Family effects (intercepts and slopes) accounted for 13.9% of the variance. The reduced model 2 without the ‘taxon’ term gave a much worse fit to the data, resulting in a highly significant likelihood ratio test (χ2 = 19.12, 1 df, p-value = 1.228 x 10−05). There is thus clear evidence that B. variegata tadpoles were at relatively greater risk of predation by Aeshna cyanea larvae. Nevertheless, a large amount of unexplained variation remains (Fig 3).

[Figure omitted. See PDF.]

Table 2. Linear mixed models of predation risk per rearing compartment.

https://doi.org/10.1371/journal.pone.0231804.t002

Tadpole size and shape

Throughout the experiment, B. bombina was on average smaller with a shallower tail (Fig 4A and 4B). But after correction for body size, B. bombina's tail fin was relatively higher than that of B. variegata (Fig 4C). Note that this is a difference between two types of predator-induced phenotypes. It does not tell us about the strength of the plastic response in either taxon, because predator-naïve tadpoles were not part of this study. Our shearing approach to correct for body size is validated by the shared allometry in both taxa (S2 Fig). If size and/or tail shape are functionally important, we would expect to see within-taxon correlations with survival as well. The logistic regression analysis (Table 3, see Materials and Methods for details) confirmed the significantly lower survival of B. variegata and also showed a strong interaction of tadpole size (PC1) and taxon. There were no significant effects associated with relative tail height.

[Figure omitted. See PDF.]

Fig 4. Morphological development by taxon.

Means and standard deviations by taxon for each of three descriptors of morphology. A. PC1 = 1. principal component ≈ body size, B. log(tail height), C. residual(log(tail height)) = residuals of a regression of log(tail height) on PC1 ≈ size-corrected tail height.

https://doi.org/10.1371/journal.pone.0231804.g004

[Figure omitted. See PDF.]

Table 3. Coefficients of the logistic model of survival.

https://doi.org/10.1371/journal.pone.0231804.t003

Fig 5A illustrates that larger B. bombina tadpoles survived better than smaller ones (t = 3.5993, df = 85.589, p-value = 0.0005). But this effect is not consistent across the whole dataset. B. variegata tadpoles are overall larger and more at risk of predation, whereas smaller size is associated with greater risk in B. bombina. There is therefore no evidence that a general preference of Aeshna larvae for larger prey causes higher mortality in B. variegata tadpoles. The analogous plot for relative tail height hints at an advantage for higher tail fins (Fig 5B), but the present dataset provides no statistical support for this hypothesis.

[Figure omitted. See PDF.]

Fig 5. Means of body size and tail shape by taxon and survival.

Plotted are the fitted means from the logistic regression analysis for PC1 (A) and residual log(tail height) (B). Error bars indicate one standard error of the mean.

https://doi.org/10.1371/journal.pone.0231804.g005

Tadpole activity

During the 30 s observation intervals, both taxa spent most of their time either resting or, seemingly, feeding. No food was provided in these trials. So, the tadpoles' persistent scraping of surfaces with their mouthparts, especially during the first half of the experiment, probably provided little, if any nutrition. Swimming typically occurred in short bursts of a few seconds separated by inactive intervals. Ninety-two percent of recorded predator attacks were directed at actively swimming tadpoles.

The distributions of the time spent swimming were strongly skewed to the right in both taxa (Fig 6). As expected, B. variegata spent more time swimming than B. bombina (means: 4.27 and 1.78 s, respectively; unpaired Wilcoxon rank sum test, W = 2234.5, p-value = 0.0096). Focal tadpoles were chosen based on their clear visibility during the 30 s observation interval. Their survival status is unknown (see Materials and Methods). We can therefore not investigate within-taxon correlations between activity and survival.

[Figure omitted. See PDF.]

Fig 6. Distribution of the time spent swimming by taxon.

Each observation interval lasted for 30 s. Number of intervals: 75 and 74 for B. variegata and B. bombina, respectively.

https://doi.org/10.1371/journal.pone.0231804.g006

Discussion

In our experimental setting, tadpoles of the puddle breeder B. variegata are at greater risk of predation from dragonfly larvae than those of B. bombina that develop in semi-permanent ponds. This result confirms our expectation based on similar species comparisons [42, 43]. It also agrees with the findings of Kruuk and Gilchrist [37] who conducted predation trials with larger groups of predator-naīve B. bombina and B. variegata tadpoles that could only adjust their behaviour. Here, we allowed tadpoles to develop a more complete predator-induced phenotype. Nonetheless, the survival difference between the taxa persisted and it was particularly strong during days 5–12 (Fig 2A). In previous mesocosm experiments with pure and hybrid Bombina and free-ranging dragonflies, the overall predation rate declined sharply from day 15 onwards even though tadpoles continued to be attacked [93, see similar results for Rana, 94]. Thus, the days with the highest total predation risk overall may coincide with a particularly strong disadvantage for B. variegata.

The increase in no-predation trials in the second half of the experiment is consistent with the mesocosm observations. The overall proportion of such trials (0.27) was high in relation to comparable studies. In our previous Bombina study [82] using the same design (4 tadpoles, 1 dragonfly), all trials resulted in predation. But the trial compartments were only half the size and so presumably led to more predator-prey encounters. McCollum and Van Buskirk [75] exposed Hyla versicolor tadpoles (two per trial, 1.5 times larger boxes) to predators for up to 500 min. Both tadpoles survived in 7% of cases. Our trials may have led to a roughly similar outcome, had they been extended from 45 to 500 min. The fact that all tadpoles managed to avoid predation in a quarter of cases suggests that we did not create a completely ‘hopeless’ and thus unrealistic situation but instead one that even the more susceptible taxon could negotiate successfully. Many predator attacks were unsuccessful, because tadpoles could escape via sprints or agile movements through the vegetation.

Many tadpoles were used in more than one trial. Later trials may therefore have involved individuals that benefited from experience gained and/or were inherently better at avoiding predators. We believe that this resembles the situation in nature more closely than the alternative of shielding tadpoles from predator encounters until the day of their one experimental test. As pointed out by Melvin and Houlahan [95], typical tadpole mortality in nature (approx. 95% according to their review) may be non-random with regard to traits under study. The generally much more benign conditions in the laboratory may thus produce experimental cohorts that are not representative of the population of origin. Moreover, we specifically set out to integrate predation risk over a substantial part of larval ontogeny with all its dramatic changes. We consider any gain in experience to be part of this process. Nonetheless, this close-up view of predator-prey interactions is a piece in a mosaic that should be complemented with experiments in semi-natural settings and using additional predator species.

If B. bombina and B. variegata were reproductively isolated, these results would address the question why B. variegata reproduces only rarely in semi-permanent, predator-rich ponds [96]. Whether or not B. bombina tadpoles are present in these habitats may not matter, because there is little evidence for competition among tadpoles in ponds [41, 43, 46]. Predators tend to keep tadpole numbers sufficiently low so that competition is negligible. B. bombina then simply provides a benchmark for a successful 'pond strategy'. B. variegata's greater susceptibility to predation shows that its tadpoles fall short of this benchmark even when they display the predator induced phenotype. Other traits such as relatively fewer eggs per clutch [36] may also limit B. variegata's reproductive success in ponds.

Hybridisation shifts the focus away from a comparison between distinct types and towards selection on traits in recombinants. These genotypes predominate in Bombina hybrid zones, while F1 individuals are very rare, if present at all [24, 26]. A hybrid genome consists of blocks of varying sizes that are inherited from one or the other species. Within the hybrid zone and with each round of recombination, these blocks get smaller, while the influx of pure genomes from the periphery introduces new large blocks [97]. The chance of, say, a B. variegata-derived block to cross into the B. bombina gene pool depends on how it affects the fitness of its carriers along the way. A genetic variant that reduces fitness in ponds is likely to be eliminated from this habitat together with the entire B. variegata-block in which it resides. It is in this sense that such variants act as barriers to gene flow [98].

This reasoning extends to adaptively diverged traits, no matter whether the taxon difference is based on one locus or on several unlinked, additive loci: any maladapted B. variegata variant would increase the probability that the block in which it resides will be eliminated from ponds [99]. Smooth clines of six quantitative traits mapped across the hybrid zone near Pešćenica, Croatia, indeed suggested an largely additive genetic basis [30]. For highly integrated compound traits, the gene flow barrier may be further strengthened if intermediate hybrids have reduced fitness due to mismatched phenotypes. This has been demonstrated in F2 crosses of sticklebacks (limnetic x benthic ecotypes), where numerous quantitative trait loci underlying the trophic adaptations are known: F2 hybrids with certain intermediate genotypes expressed mismatched phenotypes and grew more slowly [100]. Overall, the contribution of ecological divergence to the maintenance of distinct taxa should increase as the list of traits under selection gets longer [101, 102].

Prey activity alerts visually hunting predators and correlates strongly with tadpole predation rates in inter-specific comparisons (see Introduction) and also within species [75, 103]. Our findings confirm this pattern. Because the taxon difference in activity is seen in animals raised under uniform laboratory conditions, it should have a genetic basis. We do not know how many loci are involved but we expect that each B. variegata-derived genetic variant that crosses the hybrid zone and increases activity in a B. bombina genetic background is selected against. It thus reduces gene flow at genetically linked loci.

For the sake of argument we consider two alternative scenarios. There may be a hybrid genotype with an activity level that increases fitness over that of either pure taxon. This type of hybrid might then contribute to an increase rather than a reduction in gene flow. However, we do not find it likely that activity has an intermediate optimum that has not yet been reached in either taxon. It could also be that a high-activity variant is closely linked to a universally favourable mutation that arose in B. variegata. The mean fitness of blocks including both variants might then be higher than that of their B. bombina counterparts, generating a conduit of gene flow in this genomic region. But the effect of the activity variant would still be to slow the introgression of such a block into the B. bombina gene pool and the ‘pull’ of the favourable variant would cease as soon as recombination separates the two [99]. On balance, these considerations do not alter our prediction that higher activity most likely contributes to the reduction in gene flow across the hybrid zone.

There was no evidence that relatively higher tail fins caused better survival in either taxon. The large variances that rendered differences between group means (taxon, survival) non-significant (Figs 4C and 5B) may reflect heritable variation among families in the amplitude and type of morphological plasticity [104]. This trait warrants further study using more detailed morphometrics [72]. Finally, our analyses disproved the hypothesis that an overall preference of Aeshna cyanea larvae towards larger prey contributed to B. variegata's greater mortality.

Is higher predation risk in B. variegata a direct consequence of adaptation to ephemeral habitat? Population studies have shown that successful metamorphosis can be restricted in a given year to less than 50% of sites in which eggs have been laid (depending on annual rainfall) and that desiccation is the main cause of reproductive failure [105–107]. This is reflected in the growth rate of B. variegata tadpoles which was the fastest in a set of 15 European and North American anurans [68]. In contrast to the situation in ponds, intra-specific competition should be commonplace for B. variegata as tadpoles grow and puddles shrink [108]. The typically high activity of tadpoles growing in this type of habitat supports this [42, 109]. In the laboratory, large B. variegata tadpoles displace smaller conspecifics from food patches with their powerful tail movements [110]. Not only is the abundance and diversity of predators that B. variegata encounters lower than in ponds [37, see also 39–42], sometimes these early-successional sites are predator-free. This suggests a hierarchy of risk: rapid development takes priority whether or not there are predators. They add to the already considerable unpredictability of reproductive success. Bet-hedging may be a way to sidestep an unresolvable trade-off [111, 112]. Adult B. variegata can live up to at least 15 years in nature [29, but see: 113] and females are known to lay eggs repeatedly over their lifespan [29, 112].

There may be aspects of B. variegata’s anti-predator strategy that remain to be discovered. B. variegata frequently develops dark tail spots that may divert attacks away from the body, whereas B. bombina's tail fin is always uncoloured [114]. Moreover, the animals may not have been able to express some behavioural strategies in our experimental setting. Nonetheless, our results show that high activity in Bombina tadpoles is associated with higher predation risk that causes relatively higher mortality in the more active taxon, B. variegata, even when the tadpoles express a predator-induced phenotype. Population genetic theory predicts that B. variegata alleles underpinning high activity will be selected against in typical B. bombina breeding habitat and should therefore act as barriers to gene flow across the hybrid zone.

In the context of ecological speciation (cf. Introduction), the remaining question is whether these taxon differences reflect selection towards different phenotypic optima. Given that quiescent behaviour benefits B. bombina, what are the benefits of high activity for B. variegata? If indeed rapid growth and development in a competitive environment require active tadpoles, then the classic trade-off applies: activity mediates two mutually exclusive strategies [42, 45, 46, 51, 52]. In Bombina, this could be tested not just by comparing the pure taxa, but also by investigating phenotypic and genetic covariances in natural hybrids and by mapping the component traits in the Bombina genome [23].

Supporting information

[Figure omitted. See PDF.]

S1 Fig. Risk context per rearing compartment.

https://doi.org/10.1371/journal.pone.0231804.s001

(PDF)

S2 Fig. Pairwise plots of morphological traits.

https://doi.org/10.1371/journal.pone.0231804.s002

(PDF)

S1 Table. Ecological features of the collection sites.

https://doi.org/10.1371/journal.pone.0231804.s003

(PDF)

Acknowledgments

We thank Oliver Gast and Zuzana Hiadlovská for help in the field and with the experiments, respectively. Jarrod Hadfield, Lumír Gvoždík and Crispin Jordan offered advice on the experimental design and statistics. We thank Günther Gollmann, Jacek M. Szymura and Josh Van Buskirk for comments on an earlier version of the manuscript.

Citation: Smolinský R, Baláž V, Nürnberger B (2020) Tadpoles of hybridising fire-bellied toads (B. bombina and B. variegata) differ in their susceptibility to predation. PLoS ONE 15(12): e0231804. https://doi.org/10.1371/journal.pone.0231804

References

1. Hereford J. A quantitative survey of local adaptation and fitness trade-offs. Am Nat. 2009;173: 579–588. pmid:19272016

2. Schluter D. Ecology and the origin of species. Trends Ecol Evol. 2001;16: 372–380. pmid:11403870

3. Schluter D. Evidence for ecological speciation and its alternative. Science. 2009;323: 737–741. pmid:19197053

4. Hendry AP. Ecological speciation! Or the lack thereof? Can J Fish Aquat Sci. 2009;66: 1383–1398.

5. Nosil P. Ecological Speciation. Oxford: Oxford University Press; 2012.

6. Bierne N, Gagnaire P-A, David P. The geography of introgression in a patchy environment and the thorn in the thigh of ecological speciation. Curr Zool. 2013;59: 72–86.

7. Faria R, Renaut S, Galindo J, Pinho C, Melo-Ferreira J, Melo M, et al. Advances in ecological speciation: an integrative approach. Mol Ecol. 2014;23: 513–521. pmid:24354648

8. Filchak KE, Roethele JB, Feder JL. Natural selection and sympatric speciation in the apple maggot fly, Rhagoletis pomonella. Nature. 2000;407: 739–742. pmid:11048719

9. Linn C, Feder JL, Nojima S, Dambroski HR, Berlocher SH, Roelofs W. Fruit odor discrimination and sympatric host race formation in Rhagoletis. Proc Natl Acad Sci USA. 2003;100: 11490–11493. pmid:14504399

10. McGee MD, Schluter D, Wainwright PC. Functional basis of ecological divergence in sympatric stickleback. BMC Evol Biol. 2013;13: 277. pmid:24380474

11. Le Pennec G, Butlin RK, Jonsson PR, Larsson AI, Lindborg J, Bergström E, et al. Adaptation to dislodgement risk on wave-swept rocky shores in the snail Littorina saxatilis. PLOS ONE. 2017;12: e0186901. pmid:29059225

12. Johannesson B. Shell morphology of Littorina saxatilis Olivi: The relative importance of physical factors and predation. J Exp Mar Biol Ecol. 1986;102: 183–195.

13. Schluter D. Adaptive Radiation in Sticklebacks: Trade-Offs in Feeding Performance and Growth. Ecology. 1995;76: 82–90.

14. Janson K. Selection and migration in two distinct phenotypes of Littorina saxatilis in Sweden. Oecologia. 1983;59: 58–61. pmid:25024147

15. Linn CE, Dambroski HR, Feder JL, Berlocher SH, Nojima S, Roelofs WL. Postzygotic isolating factor in sympatric speciation in Rhagoletis flies: Reduced response of hybrids to parental host-fruit odors. Proc Natl Acad Sci U S A. 2004;101: 17753–17758. pmid:15591346

16. Hatfield T, Schluter D. Ecological speciation in sticklebacks: environment-dependent hybrid fitness. Evolution. 1999;53: 866–873. pmid:28565618

17. Feder JL, Opp SB, Wlazlo B, Reynolds K, Go W, Spisak S. Host fidelity is an effective premating barrier between sympatric races of the apple maggot fly. Proc Natl Acad Sci USA. 1994;91: 7990–7994. pmid:11607491

18. Rundle HD, Nagel L, Boughman JW, Schluter D. Natural selection and parallel speciation in sympatic sticklebacks. Science. 2000;287: 306–308. pmid:10634785

19. Bay RA, Arnegard ME, Conte GL, Best J, Bedford NL, McCann SR, et al. Genetic coupling of female mate choice with polygenic ecological divergence facilitates Stickleback speciation. Curr Biol. 2017;27: 3344–3349.e4. pmid:29056455

20. Feder JL, Roethele JB, Filchak KE, Niedbalski J, Romero-Severson J. Evidence for inversion polymorphism related to sympatric host race formation in the apple maggot fly, Rhagoletis pomonella. Genetics. 2003;163: 939–953. pmid:12663534

21. Meyers PJ, Doellman MM, Ragland GJ, Hood GR, Egan SP, Powell THQ, et al. Can the genomics of ecological speciation be predicted across the divergence continuum from host races to species? A case study in Rhagoletis. Philos Trans R Soc B Biol Sci. 2020;375: 20190534. pmid:32654640

22. Szymura JM. Analysis of hybrid zones with Bombina. In: Harrison RG, editor. Hybrid zones and the evolutionary process. New York: Oxford University Press; 1993. pp. 261–289.

23. Nürnberger B, Baird SJE, Čížková D, Bryjová A, Mudd AB, Blaxter ML, et al. A dense linkage map for a large repetitive genome: discovery of the sex-determining region in hybridising fire-bellied toads (Bombina bombina and B. variegata). BioRxiv [Preprint]. 2020 bioRxiv 328633 [posted 9 Oct 2020; cited 21 Oct 2020]: [51 p.]. Available from: https://www.biorxiv.org/content/10.1101/2020.10.06.328633v1

24. Yanchukov A, Hofman S, Szymura JM, Mezhzherin S, Morozov-Leonov SY, Barton NH, et al. 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: 583–600. pmid:16637503

25. Barton NH, Gale KS. Genetic analysis of hybrid zones. In: Harrison RG, editor. Hybrid zones and the evolutionary process. Oxford: Oxford University Press; 1993. pp. 13–45.

26. Szymura JM, Barton NH. The genetic structure of the hybrid zone between the fire-bellied toads Bombina bombina and B. variegata: comparison between transects and between loci. Evolution. 1991;45: 237–261. pmid:28567861

27. Kruuk LEB, Gilchrist JS, Barton NH. Hybrid dysfunction in Fire-Bellied Toads (Bombina). Evolution. 1999;53: 1611–1616. pmid:28565554

28. Philippi D, Gollmann G. On the status of the fire-bellied toad, Bombina bombina, in Lobau (Vienna, Donau-Auen National Park). Acta ZooBot Austria. 2014;150/151: 25–32.

29. Gollmann B, Gollmann G. Die Gelbbauchunke: von der Suhle zur Radspur. Bielefeld: Laurenti Verlag; 2002.

30. Nürnberger B, Barton NH, MacCallum C, Gilchrist J, Appleby M. Natural selection on quantitative traits in the Bombina hybrid zone. Evolution. 1995;49: 1224–1238. pmid:28568525

31. Czopkowa G, Czopek J. The vascularization of the respiratory surfaces in Bombina bombina. Bull Acad Pol Sci Cl II. 1955;3: 313–318.

32. Hartel T. Movement activity in a Bombina variegata population from a deciduous forested landscape. North-West J Zool. 2008;4: 79–90.

33. Kovar R, Brabec M, Vita R, Bocek R. Spring migration distances of some Central European amphibian species. Amphib-Reptil. 2009;30: 367–378.

34. Boulenger GA. On two European species of Bombinator. Proc Zool Soc Lond. 1886;1886.

35. Lörcher K. Vergleichende bioakustische Untersuchungen an der Rot- und Gelbbauchunke, Bombina bombina (L.) und Bombina v. variegata (L.). Oecologia. 1969;3: 84–124. pmid:28309207

36. Rafińska A. Reproductive biology of the fire-bellied toads, Bombina bombina and B. variegata (Anura: Discoglossidae): egg size, clutch size and larval period length differences. Biol J Linn Soc. 1991;43: 197–210.

37. Kruuk LEB, Gilchrist JS. Mechanisms maintaining species differentiation: predator mediated selection in a Bombina hybrid zone. Proc R Soc Lond B. 1997;264: 105–110.

38. Reichwaldt E. Der Einfluss von Prädatoren auf die Kaulquappen hybridisierender Unken (Bombina, Discoglossidae, Amphibia): Entwicklung, Morphologie und Verhalten. Diploma Thesis, University of Munich. 1999.

39. Babbitt KJ, Baber MJ, Tarr TL. Patterns of larval amphibian distribution along a wetland hydroperiod gradient. Can J Zool. 2003;81: 1539–1552.

40. Richter-Boix A, Llorente GA, Montori A. A comparative study of predator-induced phenotype in tadpoles across a pond permanency gradient. Hydrobiologia. 2007;583: 43–56.

41. Werner EE, Skelly DK, Relyea RA, Yurewicz KL. Amphibian species richness across environmental gradients. Oikos. 2007;116: 1697–1712.

42. Woodward BD. Predator-prey interactions and breeding pond use of temporary pond species in a desert anuran community. Ecology. 1983;64: 1549–1555.

43. Smith DC, Van Buskirk J. Phenotyic design, plasticity, and ecological performance in two tadpole species. Am Nat. 1995;145: 211–233.

44. Relyea RA. Morphological and behavioural plasticity of larval anurans in response to different predators. Ecology. 2001;82: 523–540.

45. Lawler SP. Behavioural responses to predators and predation risk in four species of larval anurans. Anim Behav. 1989;38: 1039–1047.

46. Morin PJ. Predation, competition, and the composition of larval anuran guilds. Ecol Monogr. 1983;53: 119–138.

47. Azevedo-Ramos C, Van Sluys M, Hero J-M, Magnusson WE. Influence of tadpole movement on predation by odonate naiads. J Herpetol. 1992;26: 335–338.

48. Chovanec A. The influence of tadpole swimming behaviour on predation by dragonfly nymphs. Amphib-Reptil. 1992;13: 341–349.

49. Relyea RA. The Relationship between predation risk and antipredator responses in larval Anurans. Ecology. 2001;82: 541–554.

50. Dayton GH, Saenz D, Baum KA, Langerhans RB, DeWitt TJ, Lindström J. Body shape, burst speed and escape behavior of larval Anurans. Oikos. 2005;111: 582–591.

51. Werner EE, Anholt BR. Ecological consequences of the trade-off between growth and mortality-rates mediated by foraging activity. Am Nat. 1993;142: 242–272. pmid:19425978

52. Cressler CE, King AA, Werner EE. Interactions between behavioral and life‐history trade‐offs in the evolution of integrated predator‐defense plasticity. Am Nat. 2010;176: 276–288. pmid:20629540

53. Richardson JML. The relative roles of adaptation and phylogeny in determination of larval traits in diversifying Anuran lineages. Am Nat. 2001;157: 282–299. pmid:18707290

54. Skelly DK, Werner EE. Behavioral and life-history responses of larval American toads to an odonate predator. Ecology. 1990;71: 2313–2322.

55. McCollum SA, Van Buskirk J. Costs and benefits of a predator-induced polymorphism in the gray treefrog Hyla chrysoscelis. Evolution. 1996;50: 583–593. pmid:28568914

56. Laurila A, Kujasalo J, Ranta E. Predator induced changes in life history in two anuran tadpoles: effects of predators on diet. Oikos. 1998;83: 307–317.

57. Relyea RA, Werner EE. Quantifying the relation between predator-induced behavior and growth performance in larval anurans. Ecology. 1999;80: 2117–2124.

58. Steiner UK. Linking antipredator behaviour, ingestion, gut evacuation and costs of predator-induced responses in tadpoles. Anim Behav. 2007;74: 1473–1479.

59. Orizaola G, Richter‐Boix A, Laurila A. Transgenerational effects and impact of compensatory responses to changes in breeding phenology on antipredator defenses. Ecology. 2016;97: 2470–2478. pmid:27859081

60. Lima SC, Dill LM. Behavioral decisions made under the risk of predation: a review and prospectus. Can J Zool. 1991;268: 619–640.

61. Anholt BR, Werner E, Skelly DK. Effect of food and predators on the activity of four larval ranid frogs. Ecology. 2000;81: 3509–3521.

62. Van Buskirk J. A comparative test of the adaptive plasticity hypothesis: Relationships between habitat and phenotype in Anuran larvae. Am Nat. 2002;160: 87–102. pmid:18707501

63. Dijk B, Laurila A, Orizaola G, Johansson F. Is one defence enough? Disentangling the relative importance of morphological and behavioural predator-induced defences. Behav Ecol Sociobiol. 2016;70: 237–246.

64. Teplitsky C, Plénet S, Joly P. Tadpoles’ responses to risk of fish introduction. Oecologia. 2003;134: 270–277. pmid:12647168

65. Kurali A, Pásztor K, Hettyey A, Tóth Z. Resource-dependent temporal changes in antipredator behavior of common toad (Bufo bufo) tadpoles. Behav Ecol Sociobiol. 2018;72: 91.

66. McCollum SA, Leimberger JD. Predator-induced morphological changes in an amphibian: predation by dragonflies affects tadpole shape and color. Oecologia. 1996;109: 615–621.

67. Touchon JC, Warkentin KM. Fish and dragonfly nymph predators induce opposite shifts in color and morphology of tadpoles. Oikos. 2008;117: 634–640.

68. Van Buskirk J. The costs of an inducible defense in anuran larvae. Ecology. 2000;81: 2813–2821.

69. Relyea RA. Getting out alive: how predators affect the decision to metamorphose. Oecologia. 2007;152: 389–400. pmid:17356812

70. Van Buskirk J. Natural variation in morphology of larval Amphibians: Phenotypic plasticity in nature? Ecol Monogr. 2009;79: 681–705.

71. Michel MJ. Spatial dependence of phenotype-environment associations for tadpoles in natural ponds. Evol Ecol. 2011;25: 915–932.

72. Arendt J. Morphological correlates of sprint swimming speed in five species of spadefoot toad tadpoles: Comparison of morphometric methods. J Morphol. 2010;271: 1044–1052. pmid:20730918

73. Van Buskirk J, Anderwald P, Lüpold S, Reinhardt L, Schuler H. The lure effect, tadpole tail shape, and the target of dragonfly strikes. J Herpetol. 2003;37: 420–424.

74. Johnson JB, Burt DB, DeWitt TJ. Form, function, and fitness: Pathways to survival. Evolution. 2008;62: 1243–1251. pmid:18266988

75. Van Buskirk J, McCollum SA. Functional mechanisms of an inducible defence in tadpoles: morphology and behaviour influence mortality risk from predation. J Evol Biol. 2001;13: 336–347.

76. Teplitsky C, Plenet S, Léna J-P, Mermet N, Malet E, Joly P. Escape behaviour and ultimate causes of specific induced defences in an anuran tadpole. J Evol Biol. 2005;18: 180–190. pmid:15669975

77. Kraft PG, Wilson RS, Franklin CE, Blows MW. Substantial changes in the genetic basis of tadpole morphology of Rana lessonae in the presence of predators. J Evol Biol. 2006;19: 1813–1818. pmid:17040378

78. Innes-Gold AA, Zuczek NY, Touchon JC. Right phenotype, wrong place: predator-induced plasticity is costly in a mismatched environment. Proc R Soc B Biol Sci. 2019;286: 20192347. pmid:31795869

79. Van Buskirk J, McCollum SA, Werner EE. Natural selection for environmentally induced phenotypes in tadpoles. Evolution. 1997;51: 1981–1990. pmid:28565118

80. Van Buskirk J, Relyea RA. Selection for phenotypic plasticity in Rana sylvatica tadpoles. Biol J Linn Soc. 1998;65: 301–328.

81. Van Buskirk J. Spatially heterogeneous selection in nature favors phenotypic plasticity in anuran larvae. Evolution. 2017;71: 1670–1685. pmid:28346658

82. Vorndran IC, Reichwaldt E, Nürnberger B. Does differential susceptibility to predation in tadpoles stabilize the Bombina hybrid zone? Ecology. 2002;83: 1648–1659.

83. Fijarczyk A, Nadachowska K, Hofman S, Litvinchuk SN, Babik W, Stuglik M, et al. Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Mol Ecol. 2011;20: 3381–3398. pmid:21749513

84. Baláž V, Vojar J, Civiš P, Šandera M, Rozínek R. Chytridiomycosis risk among Central European amphibians based on surveillance data. Dis Aquat Organ. 2014;112: 1–8. pmid:25392037

85. Blooi M, Pasmans F, Longcore JE, Spitzen-van der Sluijs A, Vercammen F, Martel A. Duplex real-time PCR for rapid simultaneous detection of Batrachochytrium dendrobatidis and Batrachochytrium salamandrivorans in Amphibian samples. J Clin Microbiol. 2013;51: 4173–4177. pmid:24108616

86. Werner EE, McPeek MA. Direct and indirect effects of predators on two Anuran species along an environmental gradient. Ecology. 1994;75: 1368–1382.

87. Bradley R, Terry M. Rank analysis of incomplete block designs. I. The method of paired comparisons. Biometrika. 1952;39: 324–345.

88. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67: 1–48.

89. Winter B. Linear models and linear mixed effect models in R with linguistic applications. 2013. Available: https://arxiv.org/abs/1308.5499

90. McCoy MW, Bolker BM, Osenberg CW, Miner BG, Vonesh JR. Size correction: comparing morphological traits among populations and environments. Oecologia. 2006;148: 547–554. pmid:16604370

91. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2018. Available: https://www.R-project.org/

92. Gosner KL. A simplified table for staging anuran embryos and larvae with notes on identification. Herpetologica. 1960;16: 183–190.

93. Fischer N. Performance of pure, natural hybrid and F1 Bombina tadpoles in semi-natural enclosures: Morphology, behaviour, growth, survival and development time. Diploma Thesis, University of Munich. 2006.

94. Eklöv P, Werner EE. Multiple predator effects on size-dependent behavior and mortality of two species of Anuran larvae. Oikos. 2000;88: 250–258.

95. Melvin SD, Houlahan JE. Tadpole mortality varies across experimental venues: do laboratory populations predict responses in nature? Oecologia. 2012;169: 861–868. pmid:22298111

96. Van Buskirk J. Habitat partitioning in European and North American pond-breeding frogs and toads. Divers Distrib. 2003;9: 399–410.

97. Barton NH. Multilocus clines. Evolution. 1983;37: 454–471. pmid:28563316

98. Barton NH, Bengtsson B-O. The barrier to genetic exchange between hybridising populations. Heredity. 1986;56: 357–376. pmid:3804765

99. Sachdeva H, Barton NH. Introgression of a block of genome under infinitesimal selection. Genetics. 2018;209: 1279–1303. pmid:29895560

100. Arnegard ME, McGee MD, Matthews B, Marchinko KB, Conte GL, Kabir S, et al. Genetics of ecological divergence during speciation. Nature. 2014;511: 307–311. pmid:24909991

101. Rice WR, Hostert EE. Laboratory experiments on speciation: what have we learned in 40 years? Evolution. 1993;47: 1637–1653. pmid:28568007

102. Nosil P, Harmon LJ, Seehausen O. Ecological explanations for (incomplete) speciation. Trends Ecol Evol. 2009;24: 145–156. pmid:19185951

103. Anholt BR, Werner EE. Interaction between food availability and predation mortality mediated by adaptive behavior. Ecology. 1995;76: 2230–2234.

104. Touchon JC, Robertson JM. You cannot have it all: Heritability and constraints of predator-induced developmental plasticity in a Neotropical treefrog. Evolution. 2018;72: 2758–2772. pmid:30346031

105. Barandun J, Reyer H-U. Reproductive Ecology of Bombina variegata: Development of Eggs and Larvae. J Herpetol. 1997;31: 107–110.

106. Köhler SC. Mechanisms of partial reproductive isolation in a Bombina hybrid zone in Romania. Ph.D. Thesis, University of Munich. 2003.

107. Hartel T, Nemes S, Mara G. Breeding phenology and spatio-temporal dynamics of pond use by the Yellow-Bellied Toad (Bombina variegata) population: The importance of pond availability and duration. Acta Zool Litu. 2007;17: 56–63.

108. Newman RA. Effects of density and predation on Scaphiopus couchi tadpoles in desert ponds. Oecologia. 1987;71: 301–307. pmid:28312260

109. Relyea RA. Competitor-induced plasticity in tadpoles: Consequences, cues, and connections to predator-induced plasticity. Ecol Monogr. 2002;72: 523–540.

110. Lemcke C. Verwandtenerkennung bei Kaulquappen der Gelbbauchunke, Bombina variegata. Diploma Thesis, University of Munich. 2002.

111. Böll S. Ephemere Laichgewässer: Anpassungsstrategien und physiologische Zwänge der Gelbbauchunke (Bombina variegata) in einem Lebensraum mit unvorhersehbarem Austrocknungrisiko. Ph.D. Thesis, Julius-Maximilians-Universität, Würzburg. 2002.

112. Barandun J, Reyer H-U, Anholt B. Reproductive ecology of Bombina variegata: aspects of life history. Amphib-Reptil. 1997;18: 347–355.

113. Hantzschmann AM, Gollmann B, Gollmann G, Sinsch U. The fast–slow continuum of longevity among yellow-bellied toad populations (Bombina variegata): intrinsic and extrinsic drivers of variation. PeerJ. 2019;7: e8233. pmid:31871841

114. Thiesmeier B. Fotoatlas der Amphibienlarven Deutschlands, Österreichs und der Schweiz. Bielefeld: Laurenti Verlag; 2019.