Introduction

Competition for food, at both intra- and interspecific levels, is an intrinsic characteristic of animals’ life and a major mechanism of natural selection (Schoener 1974, 1982; Levine 1976). As a consequence, animals have evolved the ability to perceive the presence of competitors and to plastically modify their morphology and behavior accordingly (Horat and Semlitsch 1994; Relyea 2002; Pfennig et al. 2007). For example, under strong interspecific competition, tadpoles of the wood frog, Lithobates sylvaticus, develop a relatively large body, a wide mouth, and a short and shallow tail (Relyea 2002). These morphological changes are thought to improve tadpoles’ competitive ability, because large tadpoles, with their longer guts, are more efficient at converting ingested food into metabolic energy (Lindgren and Laurila 2005) and can thus sustain higher levels of feeding activity. Indeed, an increase in feeding activity is the most common behavioral response to competitors (Stephens and Krebs 1986).

Feeding activity depends on the balance between two opposing necessities: finding food and staying safe from predators. Under strong competition, the tradeoff is expected to shift to higher activity levels because the per capita food decreases, whereas the perceived safety increases (Ramirez et al. 2021). For this reason, individuals are expected to respond plastically to the presence of competitors by increasing their feeding activity. Within populations, however, there might be differences among individuals in their response to competitors. For example, all individuals might show a similar plastic response, but some of them might be more active than others independent of the number and types of competitors. These behavioral differences that are repeatable across time and/or contexts are called “animal personality” and are usually described in terms of cross-context repeatability (Réale et al. 2007; Dingemanse et al. 2010). Individuals, however, might also show context-dependent variation. For example, some individuals might respond more to conspecific than to heterospecific competitors, whereas others do the opposite or do not discriminate between the two categories. This among-individual variation in the plastic response (i.e., variation in mean change of behavior across contexts) is called “behavioral plasticity” or “individual-by-environment variation,” \(I\times E\) (Dingemanse and Wolf 2013; Stamps 2016; Houslay et al. 2018).

Personality and \(I\times E\) describe behavioral variation at the within-population level, but their effects reverberate at the higher levels of biological organization, with important ecological and evolutionary consequences (Wolf and Weissing 2012). To understand these effects, behavioral variation might be analyzed at multiple hierarchical levels, simultaneously. For example, Schirmer et al. (2020) studied inter-individual differences in boldness in two syntopic rodent species and found that personality, by promoting niche partitioning, plays an important ecological role in favoring species coexistence. Von Merten et al. (2020) adopted a multiple hierarchical approach to test the pace-of-life syndrome hypothesis. They compared the patterns of covariation between life history traits and behaviors at the individual and species levels to infer whether the same evolutionary mechanisms are involved.

In this paper, we adopt a multi-hierarchical approach, to study how differences in the intensity of interspecific competition affected behavioral variation at the individual, population, and species levels. We chose as study models two species of brown frogs, Rana latastei and R. dalmatina, whose ecology and recent history make them highly suitable for our goals. Within the European brown frogs, they form a monophyletic group (Veith et al. 2003; Yuan et al. 2016), but have different distribution ranges: R. dalmatina is widespread all over the Italian peninsula and Central and Eastern Europe, whereas R. latastei is endemic mostly to Northern Italy, where it often lives in syntopy with R. dalmatina (Ficetola and De Bernardi 2005; Sillero et al. 2014). Despite these differences in range, the populations of both species show low genetic variation, which supports the hypothesis that their ranges are the consequence of post-glacial expansions from single refugia (Ficetola et al. 2007; Vences et al. 2013). The difference in range extension causes the two species to interact asymmetrically: most populations of R. latastei are sympatric with R. dalmatina (Ficetola and De Bernardi 2005), whereas only few peripheral populations of R. dalmatina are sympatric with R. latastei. According to the sink-source hypothesis (Kirkpatrick and Barton 1997; Galipaud and Kokko 2020), local adaptations in peripheral populations are inhibited by gene flow from more central populations. For this reason, the effects of interspecific competition are expected to be stronger in R. latastei than in R. dalmatina.

We collected tadpoles of both species from either syntopic or allotopic populations, and raised them either in syntopy or allotopy. We then repeatedly tested tadpoles in open field trials (OFTs) in the presence of a caged conspecific, a caged heterospecific, or an empty cage. By adopting a “character state” approach (Houslay et al. 2018), we analyzed behavioral variation at individual, population, and species levels and addressed the following questions: (i) do species differ in their response to conspecific and heterospecific stimuli? As explained above, our hypothesis was a stronger response in R. latastei than in R. dalmatina. (ii) Does the syntopic-allotopic condition affect tadpoles’ behavior? Our hypothesis was a stronger response in syntopic R. latastei tadpoles. (iii) Do the patterns of among-individual variation in personality and plasticity differ among species and populations? In Fig. 1, we show three ways in which these patterns might differ, but we made no explicit predictions about them. As a corollary to these questions, we also investigated the association between behavior and both body size and growth. In tadpoles, body size is known to affect the plastic response to predators (Castellano and Friard 2021) and we asked (iv) whether it affects also the response to interspecific competitors. Finally, since the causal link between size and behavior runs in both directions (e.g., animals that invest more in feeding activity are expected to grow faster), we asked (v) whether there is an association between growth and behavior and how it varies between species.

Fig. 1
figure 1

Intra- and inter-individual variation in a hypothetical flexible trait. Each individual (blue, red, green) is tested three times (small circles) in the three environmental conditions (A, B, and C). Large circles are the individual environmental means. On the top-left of each panel, the variance–covariance structure of the random factor (individual) is shown. In (a), individuals differ in “personality” but show similar patterns of individual plasticity, variances and covariances are homogeneous, and between-environment correlations did not differ from one; in (b), individuals differ in personality, but differences vary between environments (i.e., there is \(I\times E\)), variances and covariance are no longer homogeneous, and between-environment correlations are lower than one; (c) as in (b) with a stronger \(I\times E\) effect

Full size image

Materials and methods

On 16th April 2020, we collected with hand nets a total of 120 tadpoles of R. latastei and R. dalmatina from three breeding sites, located in Special Areas of Conservation (SAC, Directive 92/43/ECC) in Northwest Italy. Sixty tadpoles (30 for each species) were collected in “Lanca di San Michele” (IT1110024 SAC), where the species breed in syntopy, whereas the other sixty were caught in ponds where only one species breeds: “Confluenza Po-Varaita” (IT1160013 SAC, R. latastei) and “Po morto di Carignano” (IT1110025 SAC, R. dalmatina). All tadpoles were at an early developmental stage (stages 24–25, Gosner 1960), with size ranging from 12 to 15 mm, in R. latastei, and from 11 to 13, in R. dalmatina. Since tadpoles belonged to relatively large populations (number of spawning females > 30) and since they were captured at different points in the breeding sites, they were unlikely to be siblings.

From the two populations of both species, we randomly chose 24 tadpoles and placed them into separate plastic containers (33.5 × 19 × 12 cm) in about 5.5 l of water. The 96 containers (each with one tadpole) were then arranged, in groups of eight, into 12 fiberglass troughs (217 × 40 × 15 cm) (LAMAR Udine s.r.l). Six troughs contained tadpoles of both species (syntopic treatment), with two tadpoles from each of the four populations. The other six troughs contained either one of the species, with four tadpoles from each population. The troughs were placed, in groups of two (one syntopic and one allotopic treatment), outdoor in a lawn under a shelter of 50% knitted shade cloth material, to avoid full-sun exposition. To maintain the physical isolation of tadpoles, but to allow homogeneous water flow through the containers, we cut two windows (16 × 6 cm) into the large sides of the containers and sealed them with 1-mm plastic mesh. Tadpoles were fed fish vegetable flakes ad libitum until the end of the experiment (14th May 2020), when they were returned to their native ponds, to complete larval development.

From April 30th to May 11th, we repeatedly video-recorded the activity of the 96 tadpoles inside experimental arenas. We used 12 tanks (60 × 40 × 15 cm) half filled with well water. At 80 cm above each tank, we placed a Raspberry Pi v2.1 8 MP camera, connected to a Raspberry 3 single board B3 + computer. The 12 Raspberries were connected via internet with a laptop computer, which controlled the recording activity of the 12 cameras by means of a custom-designed software written in Python 3. Technical details of the hardware and the software source codes can be found in the Github repository at https://github.com/olivierfriard/raspberry_video-recording_coordinator.

In a recording trial, a tadpole was placed into the arena, which contained one of the following acute stimuli: (i) an empty cage (ECT, empty cage treatment), (ii) a cage with a conspecific tadpole (CCT, caged conspecific treatment), or (iii) a cage with a heterospecific tadpole (CHT, caged heterospecific treatment). Tadpoles used as acute stimuli were from the syntopic population and, during the experiment, they were kept in allotopy in separate aquaria. The tadpole was first let to acclimatize inside a plastic cage, for about 5 min; then, the cage was lifted and the tadpole was free to move. Recordings had a 1280 × 720 resolution and a 10-Hz frame rate and lasted for 40 min, after which the tadpole was returned to its rearing container. Tadpoles were recorded once a day for a total of nine times (three times in each of the three acute treatments). The recordings were organized in nine daily sessions, during which all 96 tadpoles were recorded with the same acute stimulus. The daily session timeline was divided in three rounds. In the first, we tested tadpoles with ECT (30th April), and then with CCT (1st May) and with CHT (2nd May). After 2-day rest, we started the second and then the third round, adopting the same acute stimulus sequence (Bell et al. 2009). To improve cue effectiveness, the night before a daily session, we changed the water and positioned the acute stimuli in the tanks. Water was not changed during a daily session.

We analyzed 863 (out of 864) recorded videos with the semi-automatic tracking software DORIS v. 0.0.19 (https://github.com/olivierfriard/DORIS.program), an open-source program in Python, which uses the OpenCV library for image processing and a user-friendly graphical interface (GUI) to set the input parameters of the analysis. The program saves, for each video, a table with frame-by-frame Cartesian coordinates of the tracked objects. From these tables, we computed three variables: the shortest straight-line distance between the focal tadpole and the center of the cage; the inter-frame speed, which is the Euclidean distance between the tadpole positions in frames f and (f + 1) multiplied by the video frame rate; and the activity state, which is a binary variable with value “1” (“moving state”) if the inter-frame speed is greater than or equal to 2 cm/s or “0” (“resting state”), if it is less then 2 cm/s. We used the latter variable to identify movement bouts and to measure their number and duration. In this case, we applied a high-pass filter and considered only bouts longer than 0.5 s. From these variables, we derived the six parameters to describe tadpole behavior during the entire recording trial: (i) the average distance from the cage (CD); (ii) the mean speed (mSPEED) and (iii) its standard deviation (sdSPEED); (iv) the activity index (IND), defined as the proportion of frames with tadpoles in a “moving state”; (v) the mean duration of movements (mMOV); and (vi) the number of movements (nMOV).

At the beginning and at the end of the recording experiment, we photographed the tadpoles, by placing them in a Petri dish, lined with graph paper. Pictures were taken with a Raspberry Pi v2.1 8 MB camera on a Raspberry Pi model 3B +. From the pictures, we measured tadpole total length (from the tip of the head and the tip of the tail), using the open-source software ImageJ (Schneider et al. 2012). To minimize observer bias, blinded methods were used when the morphological data were collected. By plausibly assuming constant growth rate during the 2-week experiment, we computed the expected body size of tadpoles at the nine recording days, to be used as a predictor in the successive analyses (see below).

Statistical analyses

Between-species and among-individual behavioral (co)variation

We used Pearson’s bivariate correlation analyses to describe the pattern of correlations between the six behavioral parameters, separately in R. dalmatina and R. latastei tadpoles. With the exception of cage distance (CD), all variables were highly inter-correlated (see “Results”). For this reason, we carried out a principal component analysis on the correlation matrix of the five highly inter-correlated parameters from the entire dataset (i.e., including tadpoles of both species) and used the first component as a descriptor of tadpole activity (i.e., the amount of movements, AM) and the second component as a descriptor of the moving pattern (i.e., the type of movements) (see “Results”). The normalized cage distance (nCD) was considered separately and used as a descriptor of tadpoles’ spatial behavior. AM, TM, and nCD were then entered as dependent variables in successive analyses.

We carried out two series of general linear mixed effect models. In the first, we considered, as fixed factors, the predictors responsible for both inter-individual variation (species identity; allotopic or syntopic origin; the allotopic-syntopic ontogenetic treatment) and intra-individual variation (acute treatment, trial order, and body size). Trial order was included to control for possible carryover effects of repeating testing (Dochtermann 2010), whereas body size was included to control for its possible spurious effects on behavior (in fact, size differed between species) and for the effects that the increasing in size during the experiment could have eventually had on individual behavioral variation. In these analyses, we first run the full models with all the two- and three-way interactions; then, we derived a reduced model that included all the main effects and only the statistically significant interactions. The random part of these models included two categorical variables: the troughs (nested within the ontogenetic treatments), which accounted for uncontrolled differences between the experimental units (i.e., blocks); and the individual tadpoles (nested within the species, the population, and the ontogenetic treatment), which accounted for individual components of the residual variance. For both random variables, we assumed a homogeneous variance and covariance structures (random intercepts mixed models). Visual inspection of their residuals from all models was conformed to the assumption of residual normality.

The second series of analyses focused on the pattern of intra-individual variation. We used multivariate mixed models, separately on the R. dalmatina and R. latastei samples, to estimate repeatability scores and variance and covariance components of the behavioral traits in the three acute treatments (Dingemanse and Dochtermann 2013). Unlike the first series of models, which considered acute treatment as a three-level, categorical predictor, the second series of models adopted a multivariate approach and used the acute treatment to split the behavioral variables into three dependent variables (one for each acute treatment), using the trial order as a replicate-pairing criterion (Houslay et al. 2018). Moreover, the models of this second series did not include the predictors that might have affected tadpole personality (i.e., the allotopic/syntopic origin and treatment), and they included only the trial order, to statistically control for its effect on within-individual variation. These analyses were carried out on a balanced dataset, obtained by replacing the only missing value (see above) with its overall mean (i.e., with zero). These multivariate mixed models adopted a “character state” approach (Houslay et al. 2018; Mitchell and Houslay 2021), which was preferred over the more familiar “reaction norm” approach (Dingemanse et al. 2010) because the three acute treatments could not be a priori aligned along an ordinal axis. From these models, we estimated among-individual variances and all cross-treatment covariances (i.e., between CCT-ECT, CCT-CHT, and ECT-CHT). Variances were used to test for among-individual differences in behavior (personality) (Nakagawa and Schielzeth 2010), whereas covariances were used to test for among-individual variation in behavioral plasticity (\(I\times E\)). We used the Deviance Information Criterion (DIC) to compare the fits of the null models, which assumed a \(3\times 3\) (co)variance matrix with homogeneous values (that is, for any two treatments, \(i\) and \(j\), \(VAR\left[i\right]=VAR\left[j\right]=COVAR[i,j]\)) with the fits of an unconstrained matrix (that is, \(VAR\left[i\right]\ne VAR\left[j\right]\ne COVAR[i,j]\)). Results were considered to support the \(I\times E\) hypothesis if the unconstrained models fitted better than their equivalent null models (Houslay et al. 2018). All models assumed homogeneous residual variances across contexts.

We used the MCMCglmm package in R (Hadfield 2010) to fit all models. The package adopts Bayesian methods for inference based on Monte Carlo Markovian chain (MCMC) simulations. In all models, we derived non-informative priors from inverse-Wishart distributions (Hadfield 2010; Meuthen et al. 2018). In models with homogeneous variance, uninformative priors were obtained by setting the expected variance V = 1 and the degree of belief v = 0.002, whereas in models with heterogeneous variance by setting a covariance matrix with V = diag(N)*0.002 and v = N + 1 (where N is the number of parameters to be estimated and diag(N) is the identity matrix of order N). In the MCMC simulations, we set the following parameters: number of iterations = 330000; duration of the burn-in period = 30000; thinning sample interval = 50. Estimates of both the expected values of fixed effects and the variances of random factors were inferred from the modes of their posterior distributions. For each parameter, the program uses the 2.5 and 97.5 percentiles of its posterior distribution to calculate the credible interval, that is, the range that has a 95% probability to contain the true value (Faraway 2016). To calculate the null hypothesis probability that a parameter is zero (pMCMC), the program first computes the frequency of simulations with estimates either greater or lower than zero and considers pMCMC as two times the lowest of these frequencies. We considered an effect to be statistically significant if its pMCMC < 0.05. To test for differences in repeatability between species and acute treatments, we calculated differences of their posterior values, computed the credible intervals of these differences, and rejected the null hypothesis of no differences if the credible interval did not include zero.

The trace and distribution of all parameters were checked visually for autocorrelation and sampling stationary (Faraway 2016). Models were run repeatedly to check for result consistency and the Gelman-Rubin criterion was used to check for convergence of the MCMC processes.

Body size, growth, and behavior

We used two series of nested linear mixed effect models to analyze the factors accounting for variation in tadpoles’ body size. The first series of models considered initial body size as a dependent variable and compared three nested models that included, progressively, the trough (as a random factor), the species, the allotopic/syntopic origin, and the allotopic/syntopic treatment (as fixed factors). These models were run on the entire sample of 96 tadpoles. The second series considered growth (the increase in size during the experiment) as a dependent variable and used the same predictors of the first series with the inclusion of a behavioral variable (the mean amount of movements). These analyses used only 90 tadpoles, because the final size of six tadpoles was not measured. In both series, we compared nested models by considering the change in deviance with a \({\chi }^{2}\) distribution. The analyses were carried out using the package lme4 in R (Bates et al. 2015).

Results

Tadpole behavior

Table 1 shows the descriptive statistics of the behavioral variables in R. dalmatina and R. latastei tadpoles under the three acute treatments. With the exception of CD, the other five variables were always highly inter-correlated (Table 2) and in all bivariate tests, the correlation coefficients were always higher in R. dalmatina than in R. latastei tadpoles (Table 2).

Table 1 Descriptive statistics of the six descriptors of tadpoles’ behavior and of the first (Amount of movement) and second (Type of movements) principal components (with the exclusion of CD) in R. dalmatina (D) and R. latastei (L) samples
Full size table
Table 2 Between-movement descriptor correlations. Correlation coefficients are measured separately for R. dalmatina (above the diagonal) and R. latastei tadpoles (below the diagonal). *P < 0.5; **P < 0.01
Full size table

Table 3 shows the canonical loadings of the first and second principal components of the five highly inter-correlated variables (CD was excluded from the analysis). The first component explains about 76% of the variance and, since it is positively associated with all the original variables, it is a size factor that describes the amount of movements. The second component explains 14% of variance. Since it is influenced positively by the mean movement duration (mMOV) and negatively by the number of movements, it is a shape factor that describes the type of movement. Tadpoles that made many short movements scored lower on this component than those that made few but longer movements. In the successive analyses, we used these two components and CD as the dependent variables of GLMM models to analyze behavioral differences at the species, population, and individual levels.

Table 3 Canonical loading of the first and second components extracted from the correlation matrix of the five descriptors of tadpoles’ activity
Full size table

The amount of movements

The first model used as dependent variable the amount of movements (first principal component). Tadpoles are known to increase their activity in the presence of competitors. We, therefore, tested whether the increase in activity differed between species and between the allotopic-syntopic origins and the allotopic-syntopic treatments (Table 4). Independent of the acute treatments, R. dalmatina tadpoles moved more than R. latastei tadpoles (Fig. 2). However, while the amount of movements of R. dalmatina tadpoles did not vary markedly with the acute treatments, those of R. latastei tadpoles did: in the ECT, R. latastei tadpoles moved much less than R. dalmatina tadpoles, their movements increased slightly than those of R. dalmatina tadpoles in the CCT (posterior mode = 0.11, pMCMC = 0.389) and increased even more markedly in the CHT (posterior mode = 0.48, pMCMC = 0.001).

Table 4 Summary of the fixed effects on the amount of movements estimated from the posterior distributions of Bayesian GLMMs, using, as the reference level, R. dalmatina tadpoles from allotopic population, under allotopic treatment and tested in the empty-cage control experiments. The model included, as random factors, tadpole identity and troughs (results not shown). 95% credible intervals (CI) are shown in parentheses. The estimated level of significance (pMCMC) of each effect is indicated by the asterisks (*P < 0.05; **P < 0.01; *** P < 0.001). D: R. dalmatina, L: R. latastei
Full size table
Fig. 2
figure 2

Individual variation in the amount of movements of Rana dalmatina (red) and R. latastei (blue) tadpoles in three acute treatments: an empty cage (ECT), a caged conspecific (CCT), and a caged heterospecific (CHT). Small, transparent dots are the by-treatment individual means. Large markers are the means of the allotopic (circles) and syntopic (squares) populations

Full size image

Both the syntopic-allotopic ontogenetic treatment and the syntopic-allotopic origin of tadpoles did not show any net effect on the amount of movements. However, the syntopic-allotopic origin showed a marginally significant interaction with the species. In fact, the difference in the amount of movements between syntopic and allotopic R. latastei tadpoles was slightly larger than that observed in R. dalmatina tadpoles (posterior mode = 0.46, pMCMC = 0.026).

Independent of the environmental conditions, the amount of movements increased with tadpole body size and with trial order (Table 4). While the effect of body size did not vary with respect to the species and the environmental conditions, the effect of the trial order varied in relation to the acute treatments: with respect to ECT, the amount of movements decreased in CCT and slightly increased in CHT.

The type of movements

This variable describes how tadpoles organize their movements. Tadpoles that move less but for longer score high on it. The variable was affected by the species, by the acute treatment and by the syntopic-allotopic origin, but not by the syntopic-allotopic ontogenetic treatment (Table 5). Overall, it was higher in R. dalmatina than in R. latastei tadpoles (Fig. 3) and the interspecific differences increased in the syntopic population. In fact, syntopic-native R. latastei tadpoles made more and shorter movements than allotopic-native conspecific. In ECT, tadpoles of both species scored higher than in CCT and in CHT. In CCT, both species showed a similar decrease, whereas, in CHT, only R. latastei plastically changed the pattern of movements and tadpoles from the syntopic population changed more than those from the allotopic population.

Table 5 Summary of the fixed effects on the type of movements estimated from the posterior distributions of Bayesian GLMMs, using, as the reference level, R. dalmatina tadpoles from allotopic population, under allotopic treatment and tested in the empty-cage control experiments. The model included, as random factors, tadpole identity and troughs (results not shown). 95% credible intervals (CI) are shown in parentheses. The estimated level of significance (pMCMC) of each effect is indicated by the asterisks (*P < 0.05; **P < 0.01; *** P<0.001). D: R. dalmatina, L: R. latastei
Full size table
Fig. 3
figure 3

Individual variation in the type of movements of Rana dalmatina (red) and R. latastei (blue) tadpoles in the empty (ECT), caged conspecific (CCT), and caged heterospecific (CHT) treatments. Small, transparent dots are the by-treatment individual means. Large markers are the means of the allotopic (circles) and syntopic (squares) populations

Full size image

While body size did not show any effect on the type of movements, the trial order showed a negative association, which, in both species, was stronger in ECT than in CCT and CHT.

The spatial behavior

The tadpoles’ average cage distances were affected neither by the syntopic-allotopic origins nor by the ontogenetic treatment, but only by the interaction between species and acute treatments. Specifically, R. latastei tadpoles stayed closer to the cage when it contained a conspecific than when it was empty or with a R. dalmatina tadpole, whereas the R. dalmatina tadpoles showed the opposite pattern, staying closer to the cage when it was empty or with a R. latastei tadpole than when it contained a conspecific (Fig. 4, Table 6). Differently of what observed in the amount and in the type of movements, body size and trial order did not show any significant effect on cage distances.

Fig. 4
figure 4

Individual variation in the observed cage distances of Rana dalmatina (red) and R. latastei (blue) tadpoles in the three acute treatments (empty cage, ECT; caged conspecific, CCT; and caged heterospecific, CHT). Small, transparent dots are the by-treatment individual means. Large markers are the means of the allotopic (circles) and syntopic (squares) populations

Full size image
Table 6 Summary of the fixed effects on cage distances estimated from the posterior distributions of Bayesian GLMMs, using, as the reference level, R. dalmatina tadpoles from allotopic population, under allotopic treatment and tested in the empty-cage control experiments. The model included, as random factors, tadpole identity and troughs (results not shown). 95% credible intervals (CI) are shown in parentheses. The estimated level of significance (pMCMC) of each effect is indicated by the asterisks (*P < 0.05; **P < 0.01). D: R. dalmatina, L: R. latastei
Full size table

Within-individual variation and between-individual correlation

We used multivariate mixed models to analyze among-individual behavioral variation, separately for the two species (see “Methods”). With only one exception (cage distances in R. latastei), models with a constrained variance structure performed better than those with an unconstrained variance structure (Table 7), providing evidence that repeatabilities did not differ significantly within species, among treatments. The repeatability of the amount of movements was slightly higher in R. dalmatina than in R. latastei, but the posterior mode of their differences did not differ significantly from zero (posterior mode of the difference between R. dalmatina and R. latastei = 0.045, 95% CI: − 0.113 to 0.189). The repeatability of the shape of movements was slightly higher in R. latastei than in R. dalmatina, but, even in this case, the difference was not significant (posterior mode of the differences between R. dalmatina and R. latastei = − 0.005, 95% CI: − 0.149 to 0.121). With respect to the amount and the type of movements, cage distances showed significantly lower repeatability in both species. In R. dalmatina, repeatability did not differ from zero, whereas in R. latastei, it was low but significantly higher than zero. In this species, the unconstrained model fitted better than the constrained alternative, although the credible intervals of repeatability in the three treatments largely overlapped.

Table 7 Repeatabilities of behaviors in the three acute treatments, with their 95% credible intervals. DIC is the Deviation Information Criterion of the constrained model, whereas ΔDIC is the difference between the constrained and the unconstrained DIC values. For negative ΔDICs, repeatability was computed on the constrained model, which assumes homogeneous variance of acute treatments (empty-cage, ECT; caged conspecific, CCT; and caged heterospecific, CHT) and between-treatment correlations equal to + 1. For positive ΔDIC, repeatability was computed for each acute treatment, from the unconstrained model
Full size table

Body size, growth, and behavior

Previous analyses have shown a significant positive effect of body size on the amount of movements. However, the causal arrow that links body size to behavior runs in both directions, because if differences in the amount of movements reflect differences in feeding activity, then the higher the activity of a tadpole the higher its expected growth rate. To test for this prediction, we analyzed differences in body size between species, populations, and treatments and we investigated the pattern of association between individual behavior and growth rate.

At the beginning of the experiment, R. latastei tadpoles were significantly larger than R. dalmatina tadpoles (\({\chi }_{1}^{2}\) = 18.093, P < 0.001), and the differences in size were larger between allotopic than between syntopic populations (\({\chi }_{1}^{2}\) = 6.401, P = 0.011) (Table 8). During the experiment, R. dalmatina tadpoles grew faster than R. latastei tadpoles (\({\chi }_{1}^{2}\) = 49.537, P < 0.001) and, independent of the species, tadpoles from the syntopic population grew faster than those from the allotopic populations (\({\chi }_{1}^{2}\) = 9.719, P = 0.002). In contrast, independent of both species and origin, tadpoles of the syntopic treatment grew as fast as those of the allotopic treatment (\({\chi }_{1}^{2}\) = 0.005, P = 0.944). Independent of species, origin, and treatment, tadpoles that moved more during the recording experiments were also those that grew faster (\({\chi }_{1}^{2}\) = 5.413, P = 0.020) and the effect did not differ between the two species (\({\chi }_{1}^{2}\) = 0.609, P = 0.406) (see also Fig. 5).

Table 8 Body size of tadpoles at the beginning and at the end of the experiment. Growth is the difference between final and initial sizes (A, allotopic; S, syntopic)
Full size table
Fig. 5
figure 5

Relationship between tadpoles’ growth (increase in size from the beginning to the end of the experiment) and amount of movement (cross-context mean of the first principal component of tadpoles’ activity). Tadpoles of R. dalmantina moved more and grew faster than those of R. latastei. Independently of the species, however, the increase in size was positively associated with the amount of movements

Full size image

Discussion

R. latastei and R. dalmatina have a similar ecology and evolutionary history, but differ greatly in their range extension: R. latastei is endemic to Northern Italy, whereas R. dalmatina has undergone a successful post-glacial expansion and today is widespread in the low-plain territories of much of Central and Western Europe (Ficetola et al. 2007; Vences et al. 2013). As a consequence, R. latastei is sympatric to R. dalmatina in most of its range, whereas R. dalmatina is sympatric to R. latastei only in a small portion of it (Fig. S1). Since, in widespread species, gene flow from central populations is expected to dampen local adaptations in peripheral populations (Kirkpatrick and Barton 1997; Nicolaus et al. 2015), we made two predictions. First, we predicted R. dalmatina to be less responsive to heterospecifics than R. latastei. Our results were consistent with this prediction. In fact, at the population level, we found that R. dalmatina tadpoles were less flexible than R. latastei tadpoles and did not modify their behavior in the presence of heterospecifics. Second, we predicted R. latastei tadpoles from syntopic populations and treatments to be more sensitive to heterospecifics than those from their allotopic counterparts. Our results were only partially consistent with this prediction, in that we found a significant effect of the syntopic/allotopic population, but no effect of the syntopic/allotopic treatment. Besides these predictions, we found that, in both species, the pattern of variation at the species level mirrors that at the individual level (as exemplified in Fig. 1a) and there was no clear evidence for among-individual variation in plasticity (that is, no \(I\times E\)). In both species, among-individual variation in the amount of movements was positively associated with variation in body size and growth rate, suggesting that it might be a good descriptor of tadpoles’ feeding activity.

The ability of tadpoles to increase their activity in relation to the perceived increase of intraspecific and interspecific competition is well known (Anholt and Werner 1995; Relyea 2002, 2004). Tadpoles of many species do not simply respond to a decrease of resources, but to competitor-specific cues, which they use to exhibit competitor-specific, adaptive behaviors. For example, wood frogs can discriminate between intraspecific and several interspecific competitors and adjust their behavior accordingly (Relyea 2002). In tadpoles of the spadefoot toads, Spea bombifrons and S. multiplicata, intraspecific competition, in allopatry, favors resource polymorphism, with individuals of each species adopting either a carnivorous or a detritivorous diet, whereas, in sympatry, interspecific competition favors trophic monomorphism and divergence between species, with tadpoles of S. multiplicata adopting a detritivorous diet and those of S. bombifrons adopting a carnivorous diet (Pfennig et al. 2007). Our results not only confirm the importance of competitor cues to regulate tadpoles' behavior, but they also provide evidence that the behavioral responses to these cues varied between species. While the analyses of spatial behavior suggest that both species could discriminate between conspecific and heterospecific stimuli, the analyses of the amount and type of movements show that only R. latastei tadpoles flexibly modify their behaviors in response to conspecific and heterospecific cues. The amount of movements of R. latastei tadpoles was low in control treatments (ECT), increased in the presence of conspecifics (CCT), and increased even more in the presence of heterospecifics (CHT). The type of movements showed the opposite trend. Interestingly, the amount of movements of R. latastei tadpoles was also affected by their allotopic-syntopic origin, but not by the allotopic-syntopic treatment: tadpoles from the syntopic population tended to be more active than those from the allotopic population. Overall, these results provide strong evidence that R. latastei (but not R. dalmatina) tadpoles behave flexibly to intra- and interspecific cues, as predicted by the source-sink hypothesis. They also suggest that flexibility might be influenced by the environment experienced in the early stages of development (Ferrari et al. 2005; Urszan et al. 2018). However, since the experiment had no replicates for the allotopic-syntopic conditions of natural populations, the role of the environment in the development of the flexible response to interspecific cues needs further investigation.

Results also showed that there is more to behavioral variation than that predicted by the source-sink hypothesis. In particular, at the species level, R. dalmatina tadpoles were more active than R. latastei tadpoles and their movements, although similar in number, were longer in duration. Moreover, between-species differences depended on the type of stimulus: with respect to controls, in the presence of caged heterospecific, between-species differences in the amount of movements decreased, whereas differences in the type of movements increased. In both cases, differences were mainly explained by R. latastei tadpoles, which, as we mentioned above, behaved more flexibly than R. dalmatina tadpoles.

Moving is costly, because it consumes energy and increases predation risk, but it positively impacts feeding rate (Werner and Anholt 1993). Tadpoles’ activity thus is found to increase in presence of competitors (Relyea 2002; Cabrera-Guzman et al. 2013) and to decrease in presence of predators (Relyea 2001; Van Buskirk 2001; Urszan et al. 2018; Castellano and Friard 2021). Our finding that R. dalmatina tadpoles were more active than R. latastei tadpoles under all three acute treatments might be the consequence of a different tradeoff between the costs and the benefits of moving (Nathan et al. 2008). Tadpoles, however, differ not only in the quantity but also in the type of movements: R. latastei tadpoles tended to stay longer at the bottom of the tank and make short movements, while R. dalmatina tadpoles were more likely to swim through the water column and make longer movements. We speculate that these differences in the amount and type of movements might reflect some fine-scale differences in the ecology and life history of the two species (Biro and Stamps 2008; Bolnick et al. 2011). On hatching, R. dalmatina tadpoles are much smaller than R. latastei tadpoles (Castellano S. unpublished results) and this might make advantageous for them to invest more in feeding. Consistently with this hypothesis, we observed that at the beginning of our experiment, R. dalmatina tadpoles were smaller than R. latastei tadpoles, but they grew faster and they were of similar size when the experiment terminated. Intriguingly, the positive association between the type of movements and growth at the between-species level was also observed at the within-species level, supporting the hypothesis that the same proximate mechanism is at work at both levels of biological variation. A candidate mechanism for this association is metabolic rate (Campos-Candela et al. 2019). Biro and Stamps (2010) suggested that individuals with higher metabolic rates might be expected to grow faster and to be more active and bolder during feeding than individuals with lower metabolic rates. In tadpoles, metabolic rate is responsible for variation in growth and development rates (Beck and Congdon 2000) and it is plausibly associated with feeding efficiency and activity. For example, tadpoles are known to increase their feeding efficiency by developing larger mouths (Relyea 2002) and longer guts (Lindgren and Laurila 2005). In this way, they can plastically increase the capacity of their metabolic machinery and, thus, the ability to sustain costly behaviors. If this explanation is correct, then behavioral differences observed between R. latastei and R. dalmatina tadpoles might be viewed as different solutions along the fast-slow continuum, with R. dalmantina showing a faster life history than R. latastei.

These interspecific differences in larval life history might also explain why R. latastei was more plastic than R. dalmatina. In fact, the environmental conditions that favor “slow” pace-of-life strategies are also those that promote the evolution of plasticity (Wright et al. 2019). At the species level, “fast” species with a short life span were found to evolve less plastic responses and to invest less in learning than “slow” species with a longer life span (Kokko and Sutherland 2001). At the within-population level, fast-type individuals, with a bold and aggressive personality, were found to behave in a less plastic way than slow-type individuals (Nicolaus et al. 2015; Spiegel et al. 2015). For this reason, the observed differences in behavioral plasticity between R. latastei and R. dalmatina tadpoles might reflect their different larval life history: R. latastei tadpoles have evolved “slower” personality traits, which make them less active than R. dalmatina tadpoles, but more sensitive to the environmental cues of interspecific competition. Both the “pace-of-life” and the “source-sink” hypotheses predict R. dalmatina to be less plastic than R. latastei. Although the two hypotheses are not mutually exclusive, their relative importance can be tested empirically. In fact, while the “pace-of-life” hypothesis predicts a generalized low plasticity in R. dalmatina, the “source-sink” hypothesis predicts R. dalmatina tadpoles to be insensitive to cues relevant only in the sink, but not in the source, populations. An empirical test to discriminate between the two hypotheses could thus be to compare the plastic response of R. dalmatina and R. latastei tadpoles to predators. While the “pace-of-life” hypothesis predicts R. dalmatina tadpoles to be bolder and less sensitive to them than R. latastei tadpoles, the “source-sink” hypothesis predicts similar responses, if the predator is widespread over the range of both species.

So far, we have considered interspecific differences at the species level. But how these differences are related to the pattern of variation at the individual level? As shown in Fig. 1, when there is \(I\times E\) interaction, changes at the population level can no longer predict those at the individual level. For example, the populations of the spadefoot toads, Spea bombifrons and S. multiplicata, released from interspecific competition, increase their niche width, not because individuals become more generalist, but because they specialize and diversify their feeding strategies (Pfennig et al. 2006, 2007). Three-spine sticklebacks show a similar effect when released from competition with juvenile trout Oncorhynchus clarki, but not with prickly sculpin Cottus asper (Bolnick et al. 2010). In this case, individuals increase but do not diversify the range of resources they exploit, with no effects on the population niche width (Bolnick et al. 2010). Our results provide no evidence for an \(I\times E\) interaction and, consequently, they show consistent behavioral variation at both the individual and population levels. In the case of cage distances, all individuals in the population responded in a similar way to environmental cues (i.e., no personality effects), whereas in the case of the amount (both species) and the type of movements (only R. latastei), individuals differed, but differences were consistent across acute treatments. Moreover, we found no evidence that selection that promoted plasticity at the species level reduced among-individual differences at the population level. In the amount of movements, we observed a slightly higher repeatability in R. dalmatina, whereas in the type of movements, we observed the opposite effect, but in both cases differences were not statistically significant.

In conclusion, as predicted by the source-sink hypothesis, R. dalmatina tadpoles were less sensitive to the presence of R. latastei tadpoles than vice versa. In both species, the pattern observed at the species level largely mirrored that observed at the individual level, either because personality has no effects or because its effects are consistent across different contexts. Besides these predicted differences in context-dependent behavioral variation, however, the two species also showed unpredicted, context-independent differences, suggesting that there was more to interspecific variation than the simple effect of selection by heterospecific competitors. We hypothesize that interspecific differences are the consequences of different larval life history strategies, along the fast-slow continuum. With respect to R. latastei, R. dalmatina might have evolved faster larval phenotypes and this, by prioritizing growth over survival, might have further constrained the evolution of plastic responses to heterospecific competitors.