Maintenance of phenotypic plasticity is linked to oxidative stress in spadefoot toad larvae

www.oikosjournal.org OIKOS Oikos 1 © 2022 The Authors. Oikos published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Bob Wong Editor-in-Chief: Dries Bonte Accepted 6 February 2022 doi: 10.1111/oik.09078 1–11 2022: e09078 Phenotypic plasticity allows organisms to improve the match between their phenotype and heterogeneous environments. Theoretical models have argued that costs of maintaining the sensory and response machinery necessary for adaptive phenotypic plasticity are important determinants to the evolution of plasticity. Despite recurrent arguments invoking putative metabolic costs associated with maintenance of cellular machinery, no studies have yet attempted to quantify it from a molecular standpoint. Here we experimentally examine physiological differences across genotypes (sibships) of spadefoot toad larvae with different degrees of plasticity in response to predator cues. We observed marked differences across sibships in developmental, growth and morphological responses to predators, and tested whether increased plasticity was associated with oxidative stress or immune suppression. We observed that more plastic sibships experienced higher antioxidant enzymatic activity when reared in the absence of predator cues, i.e. not expressing their plastic responses. The degree of plasticity was also associated with higher lipid peroxidation and slightly greater granulocyte-tolymphocyte ratio. Higher antioxidant activity in highly plastic sibships suggests that maintenance of phenotypic plasticity may be linked to energy demanding metabolic processes. Our findings suggest that having the potential to produce plastic responses may incur oxidative and immunological costs. In the long term, such maintenance costs may erode individual fitness and even constrain the evolution of plasticity. To our knowledge, this is the first empirical evidence indicating the existence of a physiological cost to the maintenance of phenotypic plasticity. Keywords: antipredator responses, developmental plasticity, eco-physiology, immune response, oxidative stress, phenotypic plasticity Introduction Phenotypic plasticity can favour adaptation and diversification in heterogeneous environments (Charmantier et al. 2008). The evolution of plasticity, however, may be limited by intrinsic physiological costs associated with the ability to perceive environmental cues and arm an appropriate phenotypic response, i.e. by its ‘maintenance costs’ Maintenance of phenotypic plasticity is linked to oxidative stress in spadefoot toad larvae Pablo Burraco, Miguel A. Rendón, Carmen Díaz-Paniagua and Ivan Gomez-Mestre P. Burraco (https://orcid.org/0000-0002-9007-2643), M. A. Rendón (https://orcid.org/0000-0002-3735-7391), C. Díaz-Paniagua (https://orcid.org/00000002-6717-0173) and I. Gomez-Mestre (https://orcid.org/0000-0003-0094-8195) ✉ ([email protected]), Doñana Biological Station (CSIC), Seville, Spain. PB also at: Inst. of Biodiversity, Animal Health and Comparative Medicine, Univ. of Glasgow, Glasgow, UK. Research 2 (DeWitt et al. 1998, Auld et al. 2010, Ledon-Rettig et al. 2010). Such costs refer solely to the maintenance itself of the mechanisms needed to detect and respond to environmental inputs, and are different from ‘production costs’, which are those paid by organisms during the actual production or expression of alternative phenotypes (DeWitt et al. 1998, Teplitsky et al. 2005). The genetic and sensory machinery regulating phenotypic plasticity and their maintenance is still incomplete and therefore further empirical data is needed (Beldade et al. 2011, Kelly et al. 2011, Reuben and Touchon 2021). Theoretical models show that maintenance costs of plasticity would hinder the evolution of plasticity so that plasticity would only be favoured when its benefits outweigh its costs, especially in highly heterogeneous environments (Scheiner and Berrigan 1998, Gomez-Mestre and Jovani 2013, Chevin and Lande 2015). Various studies in the last two decades have tried to quantify maintenance costs of plasticity, although only a few of them have succeeded at empirically detecting such costs (DeWitt 1998, Agrawal et al. 2002, Relyea 2002a). A meta-analysis conducted on 27 studies across animals and plants reported 536 separate traits for which maintenance plasticity costs were examined (Van Buskirk and Steiner 2009). This meta-analysis concluded that fitness costs of plasticity were relatively mild (28.6% of the total negative fitness selection coefficients). Costs only seemed to be large when plasticity was induced by stressful environmental conditions causing a reduction in fitness (Van Buskirk and Steiner 2009). Hence, costs of plasticity may appear to be of little importance or at least hard to detect (Van Buskirk and Steiner 2009, Murren et al. 2015). Nevertheless, species with high investment in structures key to display plastic responses such as large brains or complex immune responses imply considerable maintenance costs (Snell-Rood 2012). Maintenance costs of plasticity have usually been sought evaluating their potential impact on direct or indirect measures of fitness, such as growth, developmental rate, survival or fecundity (Scheiner and Berrigan 1998, Dorn et al. 2000, Van Kleunen and Fischer 2007, Auld et al. 2010). Variation in such traits is likely under considerable selection and therefore it is perhaps not surprising that associated costs have been hard to find. Such broad fitness-related phenotypic consequences of plasticity maintenance may become quickly buffered by selection, whereas more subtle costs may persist. Auld et al. (2010) suggested that ‘for maintenance costs, plastic individuals must invest resources in maintaining the molecular physiological ‘machinery’ needed to detect, monitor and respond to environmental conditions’. This expectation of costly energy allocation or physiological toll associated with the maintenance of sensory machinery is commonly found in studies of plasticity, and yet, to the best of our knowledge, no study has explicitly attempted to quantify maintenance costs from a physiological perspective. If the maintenance of highly plastic genotypes requires an increased metabolic effort, they are likely to experience increased oxidative stress and require mechanisms to reduce the damage caused by excess reactive oxygen species (Pamplona and Costantini 2011, Halliwell and Gutteridge 2015). Indeed, intraspecific variation in metabolic rate has been suggested to explain phenotypic responses to environmental change, likely due to among-individual differences in resource allocation under unstable environments (Norin and Metcalfe 2019). In animals, plastic responses are regulated by a neuroendocrine cascade whose activation implies production costs in terms of altered metabolism, life-history traits or fitness (Gervasi and Foufopoulos 2008, Wingfield and Romero 2010, Love et al. 2014). If plasticity maintenance also required increased metabolic activity and the up regulation of catabolic processes, we would expect highly plastic genotypes to have to build up their antioxidant capacity to buffer the redox costs of elevated metabolism (Monaghan et al. 2009, Costantini et al. 2011, Pamplona and Costantini 2011). Moreover, maintenance of energetically expensive processes may also divert resources from growth, energy storage or the immune system (McEwen et al. 1998, Burraco et al. 2021), potentially generating a tradeoff between plasticity and body condition or immunocompetence. Hence, the evolution of plastic phenotypes may be constrained by the existence of physiological costs, but also by the need to buffer the costs of producing a plastic response (i.e. production costs). Here we used spadefoot toad larvae to test for an association between the degree of plasticity in response to predation risk and oxidative stress, antioxidant capacity, immune status, growth and development. We hypothesised that if maintaining adaptive antipredator plasticity is costly, highly plastic sibships would pay a physiological toll in terms of oxidative damage or compromised immune status, even if such plastic potential is not expressed. Also, we expected plastic genotypes to experience higher metabolic rates and therefore to upregulate buffering mechanisms such as the activity of antioxidant enzymes. To test these hypotheses, we collected eggs from 20 clutches of western spadefoot toads Pelobates cultripes and assessed their degree of plasticity by raising larvae to metamorphosis in the presence or absence of predator cues. We simultaneously and separately raised larvae from those same sibships under standardized conditions in the absence of predator cues, and examined among-sibship variation in oxidative stress, immune status, fat reserves or body mass. Methods We collected portions of 20 full sibships of the western spadefoot toad Pelobates cultripes from three localities in southwestern Spain: six from Doñana Natural Park (Huelva), nine from Doñana Biological Reserve (Huelva) and five from Sierra Norte Natural Park (Seville). We included sibships from different localities to increase genetic variation among-sibships and hence increase our ability to detect costs (Scheiner and Berrigan 1998), not to test hypotheses regarding geographic variation. The experimental design consisted in a two-step process. First, we quantified the degree of phenotypic plasticity across the 20 sibships in response to predator presence by 3 exposing larvae to the presence or absence of predator chemical cues (dytiscid beetle larvae). Second, we quantified several redox and immunological parameters for individuals of those same sibships under benign conditions, i.e. in the absence of predators. Finally, we combined the data from both experiments to test for an association between the degree of plasticity of each sibship and their constitutive physiological demands (Fig. 1). Eggs from each sibship were kept until hatching in 10-l plastic buckets filled with 5 l of carbon-filtered dechlorinated tap water in a walk-in climatic chamber at 18°C. Once hatched, 60 tadpoles from each sibship were haphazardly separated into 20 groups of 3 tadpoles and assigned these groups to two different treatments – control or added predator cues – to quantify their reaction norms. Similarly, 10 tadpoles from each sibship were individualized to quantify physiological parameters possibly associated with maintenance costs. Tadpoles were put in 3-l buckets (168 mm diameter × 184 mm high) filled with 2.7 l dechlorinated tap water in climatic chambers set at 24°C with a 12:12 light:dark photoperiod and fed rabbit chow ad libitum. We also collected Dytiscus circumflexus larvae as natural predators of P. cultripes tadpoles within all locations used in the experiment. Determination of reaction norms To assess the degree of plasticity in developmental rate, growth rate and morphology, we experimentally exposed tadpoles from each sibship to chemical cues from predators. Our experimental design consisted of a single factor with two treatments: absence of predator cues (‘control’) and added predator cues (‘predator’). We set 10 replicates per sibship × treatment combination (each replicate holding three tadpoles) for a total of 400 experimental units, randomly arrayed throughout three climatic chambers set at 24°C. The experiment started when tadpoles reached 25 Gosner stage (free feeding and active swimming; Gosner 1960). Water was renewed once a week. All tadpoles were fed rabbit chow ad libitum. In the predator treatment, we added water borne predator cues mixed with alarm cues from conspecific tadpoles twice a week. To obtain predator kairomones we maintained 12 Dytiscus larvae individually in 1-l buckets filled with 0.8 l dechlorinated water. We fed each Dytiscus one P. cultripes tadpole every other day. Twice a week, we filtered and pooled the water from each individual predator container to obtain a homogeneous mix and avoid biases due to predator identity. All experimental containers in the ‘predator’ treatment received 40 ml of this pool of water containing predator kairomones and alarm cues, whereas the ‘control’ containers received 40 ml of clean water. Two months after we initiated the experiment, we randomly chose one tadpole per container (for a total of ten tadpoles – all at Gosner stage 30 – per sibship per environment) to photograph it laterally for morphometric analysis of plasticity in body shape. Similarly, we estimated developmental plasticity considering only the first tadpole reaching metamorphic climax (42 Gosner stage, i.e. forelimbs emerged) from each container (ten containers per sibship per environment). Upon reaching the target developmental stage, we moved each individual to a different 1-l container filled with 50 ml of water and soaked tissue paper to provide cover until metamorphosis was completed (46 Gosner stage). Metamorphs were then weighed to the nearest 0.1 mg using a high precision balance. We used body mass at metamorphic climax for estimating growth plasticity as the linear daily gain in weight (in mg) over development. We estimated plasticity as the difference between the trait value measured in response to predator cues (i.e. time to metamorphosis, growth or morphology), and the trait value under control conditions (predator cues absent). In response Trait Env 1Env 2 Fitness Plasticity Trait Plasticity Stress physiology 1) Fitness costs of maintaining phenotypic plasticity 2) Stress-physiology parameters as cost indicators Stress ph ysiology Env 1Env 2 Env 1 G 1G 2G 3G 4 G 1 G 2 G 3 G 4 G 1 G 2 G 3 G 4 Figure 1. Schematic representation of 1) the traditional approach to the evaluation of maintenance costs of plasticity 2) the physiological approach used in this study, in which we tested for maintenance costs of antipredator plasticity in spadefoot toad Pelobates cultripes larvae. To that end, we ran two simultaneous experiments in which we determined the degree of plasticity for various traits in 20 different sibships in response to predator cues, and at the same time the baseline of several stress-related physiological parameters. Finally, we combined the data collected in each experiment to test for an association between the degree of plasticity and markers of physiological stress. 4 to predators, amphibian larvae reduce their activity rate, often resulting in longer larval periods that may result in larger size at metamorphosis (Relyea 2002b). In addition, predator cues also induce marked morphological changes such as deeper tail fins, shift of the tail insertion towards a more anterior position, or enhanced distal coloration (Touchon and Robertson 2018). Morphometric analyses We applied geometric morphometrics to describe shape variation in tadpoles across treatments. We photographed the left side of each tadpole and scaled all photos using a grid. We delimited the shape of tadpoles by digitizing 9 landmarks and 13 sliding semilandmarks (Supporting information) with tpsDig2 software (Rohlf 2010a, b). Landmarks were chosen according to their ability to capture the overall body shape of tadpoles while satisfying statistical restrictions associated with geometric morphometrics (Rufino et al. 2006). In order to control for postural changes in shape unrelated to treatment we corrected landmark position with a quadratic function using the unbend option in tpsUtil (Rohlf 2010a). We performed generalized Procrustes analysis [40] using the package geomorph ver. 3.0 in R ver. 3.5.2. Procrustes ANOVA was performed on shape variables to test if the interaction effect between treatment and sibship was significant, including log centroid size – Cs – as a covariate. Because there was an allometric component of shape variation (test for the association between Cs and shape: p = 0.001), we calculated residuals from a linear regression of shape on log(CS). We then used principal components analysis (i.e. relative warps, abbreviated as RWs) to determine allometry-free body shape variation among specimens. We included in our analyses the first four relative warps, which explained 31.0, 19.0, 11.0 and 9.4% of the total morphometric variance, respectively. These warps explained variation in common morphological features previously described for amphibian larvae (Rufino et al. 2006, Reuben and Touchon 2021). Determination of physiological parameters We estimated physiological parameters relative to metabolism, immune system and body mass on an independent set of tadpoles from the same sibships for which we had determined their reaction norms. We established ten replicates per sibship, rearing tadpoles individually in 3-l containers under identical conditions of temperature, feeding regime and photoperiod for two months. After this period, and with tadpoles at Gosner stage 30, we euthanized all individuals to determine their body mass, fat body content, activity of four antioxidant enzymes, lipid peroxidation, oxidative stress and leukocyte counts to assess immune state. For each tadpole, we blotted dry the excess water and weighed it to the nearest 0.1 mg. Tadpoles were euthanized with MS-222 (2 g l−1; 304506-5G) immediately prior to blood extraction for blood cell counts, and to dissection of fat bodies. Finally, each tadpole was eviscerated and then snap frozen and stored at −80°C until assayed for oxidative stress. Fat reserves The major triglyceride storage in amphibian larvae is located in the abdominal area forming fat bodies, which are essential for metamorphic success (Scott et al. 2007). We dissected and weighed fat bodies to the nearest 0.1 mg. Leukocyte count The immune state is susceptible to stress, which shifts the proportion of white blood cells (Davis et al. 2008). We estimated relative frequency of leukocytes and their abundance in relation to erythrocytes through flow cytometry (Uchiyama et al. 2005, Guava Easy Cyte Plus), a technique validated for Pelobates cultripes (Burraco et al. 2017a). See the Supporting information for further details. We differentiated four cell types according to gating strategy (Burraco et al. 2017a): erythrocytes, and three types of leukocytes: lymphocytes, granulocytes and monocytes. Because of the low number of monocytes in amphibian blood (Davis et al. 2008, Burraco et al. 2017a) we discarded them for further analysis. We used the granulocyte-to-lymphocyte ratio as well as the absolute count of granulocytes and lymphocytes (cells per μl of blood) to assess the immune state of individuals as higher values of this ratio are indicative of physiological stress (Rufino et al. 2006, Davis et al. 2008). Oxidative stress Eviscerated individuals were homogenized in a buffered solution (1:4; homogenate:solution) to inhibit proteolysis (100 mM Tris–HCl with 0.1 mM EDTA, 0.1% triton X-100, pH 7.8 and 0.1 mM PMSF; Burraco et al. 2013, 2017b) using a Miccra homogenizer (Miccra D-1). Homogenates for oxidative stress assessment were centrifuged at 20 817 g for 30 min at 4°C and supernatants were aliquoted into 0.6 ml tubes and stored at −80°C. We quantified the activity of four antioxidant enzymes: catalase (CAT), superoxide dismutase (SOD), glutathione reductase (GR) and glutathione peroxidase (GPx). We quantified oxidative damage in lipids by measuring the levels of malondialdehyde (MDA), and cellular redox status through estimates of the reduced-to-oxidized glutathione ratio (GSH:GSSG). See the Supporting information for details. Statistical analyses All statistical tests were conducted in R ver. 3.5.2 (<www.rproject.org>). Parametric assumptions were tested using Kolmogorov–Smirnov tests for normality and Breusch– Pagan tests for homoscedasticity. Before testing for maintenance costs of plasticity, we checked for among-sibship variation in plasticity (i.e. divergence of their reaction norms), and overall physiological differences in their oxidative and immune status, with the help 5 of linear models. We used linear mixed models (lme4 package) to check for differences in the degree of plasticity among sibships in response to predator cues testing for a ‘treatment by sibship’ interaction. Sibship was always included as a random factor in the models. MDA concentrations were logtransformed to fit parametric assumptions. We tested for maintenance costs of plasticity using general linear model selection. We conducted model selection attending to Akaike information criterion (AIC; Akaike 1973) with a correction for finite sample sizes (AICc) (Burnham 2002, Grueber et al. 2011). We only included plasticity values in the process of model selection for which we detected differences among sibships (i.e. significant ‘sibship by treatment’ interaction). Likewise, we only included physiological parameters that significantly differed among sibships. We excluded GR values in the model selection to avoid redundancy since its values were collinear with SOD and GPx (GR–SOD correlation: R2 = 0.11, p-value < 0.001, GR–GPx correlation: R2 = 0.17, p-value < 0.001). All independent variables included in the global model were standardized with the function scale to allow better comparison of the estimates. We generated a global model for each type of plastic response (developmental, growth and morphological) against those physiological variables that differed among sibships: plasticity SODCAT GPxMDA granulocyte-to-lymphocyteratio ~+++ ++ bbody mass We then used the function dredge implemented in the MuMIn package ver. 1.15.6, to generate a derived set of submodels. The analysis resulted in 64 models that were restricted to the top 2-AICc models using the function get.models implemented in the same package (Burnham 2002). We calculated the model average of the top 2AICc models with the function model.avg, recording the estimates and relative importance of each variable. For each selected model, we determined adjusted-R2, and indicated the ∆-values (AIC-differences) against the global and null models. Finally, we used linear model selection (following the same procedure indicated above) to test for differences among sibships in physiological parameters within control conditions, i.e. to check for possible physiological differences among sibships, linked to ‘sibship quality’ rather than to developmental plasticity. Results Determination of reaction norms Survival was high in both the control (77.0%) and the predator cue treatment (76.5%). Sibships differed in their degree of developmental plasticity in response to predator cues, as indicated by a significant ‘treatment × sibship’ interaction for time to metamorphosis (df = 19; χ2 = 114.75; p < 0.001; Table 1, Fig. 2a). Similarly, we found differences among sibships in how their growth was altered by predator cues (df = 19; χ2 = 35.49; p < 0.001, respectively; Table 1, Fig. 2b). Sibships also differed in their morphological plasticity, as the ‘sibship × treatment’ interaction was significant for all first four relative warps (Table 1, Fig. 2c–f). We found no evidence for among-population differences in plasticity (estimated as ‘population × treatment’ interaction) for either developmental plasticity (df = 2, χ2 = 0.63, p = 0.729), growth plasticity (df = 2, χ2 = 2.80, p = 0.246) or morphological plasticity (RW1: df = 2, χ2 = 0.38, p = 0.830; RW2: df = 2, χ2 = 5.80, p = 0.748; RW3: df = 2, χ2 = 0.58, p = 0.750; RW4: df = 2, χ2 = 3.28, p = 0.194). Determination of constitutive physiological differences among sibships We raised the larvae for two months (98% survival), during which time they developed synchronously and at the same developmental stage we euthanized them to determine a series of physiological parameters. To determine which physiological variables could be associated with differences in plasticity, we tested for among-sibship variation in those parameters. We found that body mass, the activity of four antioxidant enzymes (SOD, CAT, GPx and GR) and lipid peroxidation (MDA), as well as the granulocyte-to-lymphocyte ratio significantly differed among sibships (all p < 0.001, Supporting information). In contrast, we did not find significant differences among sibships in fat reserves, GSH, GSG-to-GSSG ratio or in the absolute count of lymphocytes or granulocytes (Supporting information). Physiological maintenance costs of plasticity The best models often indicated associations between specific physiological parameters and the degree of plasticity Table 1. Effect of predator cues on development, growth rate and morphology (first four relative warps, abbreviated RW) in Pelobates cultripes newly metamorphosed individuals (46 Gosner stage). A significant treatment-by-sibship interaction indicates differences among sibships in the degree of plasticity for a particular trait. Trait Effect df χ2p-value n Development predator × sibship 19 114.75 < 0.001 328 predator 1 5.76 0.017 328 clutch 19 95.98 < 0.001 328 Growth rate predator × sibship 19 35.49 < 0.001 307 predator 1 5.83 0.016 307 clutch 19 35.02 < 0.001 307 Morphology (RW1) predator × sibship 19 8.16 0.043 362 predator 1 3.13 0.077 362 clutch 19 8.11 0.004 362 Morphology (RW2) predator × sibship 19 24.94 < 0.001 362 predator 1 10.41 0.001 362 clutch 19 23.45 < 0.001 362 Morphology (RW3) predator × sibship 19 8.10 0.044 362 predator 1 10.80 0.001 362 clutch 19 15.83 < 0.001 362 Morphology (RW4) predator × sibship 19 20.26 0.001 362 predator 1 0.01 0.928 362 clutch 19 19.03 < 0.001 362 6 across sibships (in developmental, growth or morphology). Such associations varied depending on the trait in question (Table 2). We found that maintaining developmental and growth plasticity in response to predator cues incurred physiological consequences. Sibships showing plasticity in developmental time showed increased GPx activity (Table 2, Fig. 3A). Plasticity in growth was weakly but positively associated with changes in the granulocyte-to-lymphocyte ratio (Table 2, Fig. 3B). This relationship was strongly influenced by the G:L ratio data from one sibship (marked with an asterisk on Fig. 3B), but which presented values well within the range of G:L ratio observed for this species. Morphological plasticity represented in relative warps RW1, RW3 and RW4 also had associated variation in several physiological parameters. Plastic shifts in shape along RW1 were associated with increased GPX, SOD and CAT activities (Table 2, Fig. 3C and D). Plasticity in RW3 was associated with GPX levels and body mass (Table 2, Fig. 3E). Finally, plasticity in RW4 was associated with higher levels of GPx and MDA, and lower body mass (Table 2, Fig. 3F). We tested the relationship between among-sibship differences in development, growth or morphology and variation in physiological variables for larvae reared in the absence of predator cues. For those plastic traits that had associated maintenance costs, model selection only showed significant positive relationships between growth rate and body mass, and between RW3 and GPx, whereas the relationship with other physiological variables was negligible (Supporting information). -0.02 0.00 0.02 no pred pred -0.02 -0.01 0.00 0.01 0.02 no pred pred 200 300 400 no pred pred 0.0030 0.0035 0.0040 0.0045 0.0050 no pred pred -0.02 0.00 0.02 no pred pred (a) Days to metamorphosis (b) Growth (g/day of development) (c) Morphology (RW1) (d) Morphology (RW2) (e) Morphology (RW3) Morphology (RW4) -0.02 -0.01 0.00 0.01 0.02 no pred pred (f) Figure 2. Reaction norms of 20 Pelobates cultripes sibships responding to predator presence in terms of developmental (a), growth (b) and morphological (relative warps 1–4; c, d, e and f, respectively) plasticity. Transformation grid in the upper and lower part of c, d, e and f indicate the direction of variation of each landmark for extreme positive and negative values, respectively. 7 Discussion We find that the ability of spadefoot toad larvae to alter their development, growth and morphology in response to predation risk is associated with constitutively elevated activity of antioxidant enzymes, increased lipid peroxidation and granulocyte-to-lymphocyte ratio, as well as reduced growth rate. These results suggest that the sole maintenance of antipredator phenotypic plasticity in spadefoot toad larvae is associated with enhanced metabolism (denoted by higher activity of antioxidant enzymes), as well as to immune and oxidative costs. The activity of antioxidant enzymes was positively associated with predator induced plasticity in morphology and developmental rate. In particular, we found potential for larval plasticity to be strongly and positively associated with activity of glutathione peroxidase (GPx). One of the main catalytic functions of GPx is to reduce to water the hydrogen peroxide produced during catabolism. Hydrogen peroxide plays several important roles at the cellular level, especially in terms of ageing regulation, since increased intracellular H2O2 can induce cell death (Giorgio et al. 2007). High antioxidant capacity is often associated with increased metabolic effort (Pamplona and Costantini 2011, Halliwell and Gutteridge 2015, Castiglione et al. 2020). In amphibians, H2O2 seems to play a key role in mediating developmental transitions that involve extensive tissue remodelling (Johnson et al. 2013, Prokić et al. 2019). Additionally, GPx activity can be upregulated in individuals that have undertaken activities demanding considerable physiological efforts such as migration (Jenni-Eiermann et al. 2014), developmental acceleration (Gomez-Mestre et al. 2013) or high reproductive investment (Casagrande and Hau 2018). Therefore, increased antioxidant enzyme activity with increased plasticity may potentially indicates that more plastic sibships incur greater basal metabolic costs (i.e. even when they are not actively modifying their phenotype) and possibly that they may have the ability to buffer the overproduction of free radicals due to such enhanced metabolism. We also found lipid peroxidation (MDA levels), and, to a lesser extent, granulocyte-to-lymphocyte ratio, to be positively associated with predator-induced morphological Table 2. Estimates, standard errors (SE) and relative importance of the variables after model averaging of the top 2AICc models. Saturated models included all variables that differed significantly among sibships: superoxide dismutase (SOD), catalase (CAT), glutathione peroxidase (GPx), malondialdehyde (MDA), granulocyte-to-lymphocyte (G:L) ratio and body mass. We also tabulate the proportion of the variance explained by the best model (R2), which includes all the variables that were restricted after model averaging, as well as the delta values (differences in AIC) of this model to the saturated model and to the null model (i.e. only including the intercept). Development Estimate Unconditional SE Relative importance (Intercept) 25.29 14.42 GPx 32.73 14.79 1 best model adjusted-R2∆ to saturated model ∆ to null model ~ GPx 0.17 5.90 2.81 Growth Estimate Unconditional SE Relative importance (Intercept) 2.10×10-04 8.09×10-05 G:L ratio 1.12×10-04 1.02×10-04 0.67 best model adjusted-R2∆ to saturated model ∆ to null model ~ G:L ratio 0.15 7.68 2.24 Morphology – RW1 Estimate Unconditional SE Relative importance (Intercept) 3.75×10-03 1.89×10-03 GPx 5.98×10-03 0.001945 1 SOD 0.006876 0.002233 1 CAT 0.006643 0.002234 1 best model adjusted-R2∆ to saturated model ∆ to null model ~ GPx + SOD + CAT 0.48 4.69 10.64 Morphology – RW2 Estimate Unconditional SE Relative importance (Intercept) −5.35×10-03 1.53×10-03 best model Adjusted-R2∆ to saturated model ∆ to null model null model – – – Morphology – RW3 Estimate Unconditional SE Relative importance (Intercept) 4.15×10-03 8.85×10-04 GPx 0.0029398 0.0009138 1 body mass 0.0010003 0.0009138 0.29 best model adjusted-R2∆ to saturated model ∆ to null model ~ GPx + body mass 0.32 5.63 5.90 Morphology – RW4 Estimate Unconditional SE Relative importance (Intercept) 1.84×10-04 1.07×10-03 GPx 0.0038485 0.0013774 1 MDA 0.0065325 0.0020321 1 body mass −0.002325 0.0019435 0.28 best model adjusted-R2∆ to saturated model ∆ to null model ~ GPx + MDA + body mass 0.45 4.70 9.48 8 40 60 80 -0.01 0.00 0.01 0.02 Morphological plasticity (RW1) GPx activty (mU/mg total protein) 20 40 60 80 -0.01 0.00 0.01 0.02 Morphological plasticity (RW1) SOD activty (U/mg total protein) 40 60 80 100 0.000 0.005 0.010 0.015 Morphological plasticity (RW3) GPx activty (mU/mg total protein) 25 50 75 -100 0100 200 Developmental plasticity GPx activty (mU/mg total protein) 100 0.25 0.30 0.35 0.40 0 5e-04 1e-03 Growth plasticity Granulocyte-to-lymphocyte ratio (A) (B) (C)(D) (E) -0.010 -0.005 00.005 0.010 9 11 13 Morphological plasticity (RW4) Lipid peroxidation (nmol MDA/mL ) (F) * Figure 3. Values in the x-axis show the variations in developmental (A), growth (B) and morphological plasticity (C–F) in larvae from 20 different sibships of Pelobates cultripes exposed to predator cues. Plasticity values were calculated as the difference between the trait value measured in response to predator cues and the value in control conditions (predator cues absent). Values in the y-axis show significant variations in constitutive levels of either oxidative stress or immune related parameters measured in two-months old Pelobates cultripes larvae from the same 20 sibships used for determining plastic responses (Table 2). The asterisk in (B) indicates a datapoint with high leverage. Glossary: glutathione peroxidase (GPx), superoxide dismutase (SOD), malondialdehyde (MDA). 9 and growth plasticity, respectively. MDA is a proxy for the extent of lipid peroxidation resulting from oxidative stress (Hulbert et al. 2007). Excess lipid peroxidation usually modifies membrane composition and affects biological macromolecules like DNA, compromising cell stability and inflicting irreversible damages leading to accelerated ageing (Hulbert et al. 2007). Highly plastic sibships in our experiment experienced increased lipid peroxidation even if growing in benign conditions with ad libitum food and lack of predator cues, i.e. ‘control’ conditions. More plastic sibships may hence be paying an energetic toll in terms of excess production of free oxygen radicals due to enhanced metabolism (De Block and Stoks 2008, Murphy 2009) simply for retaining the potential to plastically alter their phenotype. Increased oxidative stress may translate into deteriorated body condition, reduced fertility, shorter lifespan and ultimately, lower fitness (Monaghan et al. 2009, Buttemer et al. 2010, Aitken et al. 2016). Similarly, we found that sibships with high plasticity in growth experienced a greater granulocyte-to-lymphocyte ratio, which is usually increased under stress and is mediated by glucocorticoid secretion in vertebrates (Davis et al. 2008), including amphibians (Burraco et al. 2017a). Increased granulocyte-to-lymphocyte ratio (or neutrophil/heterophil-tolymphocyte ratio) is also linked to higher disease susceptibility (Bhat et al. 2013) and to poor body condition (Lobato et al. 2005, Gomez et al. 2008). Therefore, growth plasticity seems to involve slight but quantifiable immunological costs in amphibian larvae responding to predation risk. These differences in physiology among sibships specifically reflect differences in their potential for predator induced developmental plasticity. When we tested for among-sibship variation in physiological parameters within the control environment, we found that differences in growth rate were only slightly associated with MDA levels, whereas developmental rate or changes in morphology showed minor or no association with physiological parameters. This indicates that trait plasticity itself, and not just constitutive trait differences among sibships, is associated with physiological maintenance costs. The existence of maintenance costs of plasticity has been proposed as one of the main causes for within-population variation in the degree of plasticity (Callahan et al. 2008, Gomez-Mestre and Jovani 2013, Lande 2014) and for the loss of plasticity in lineages evolving in isolated or stable environments compared to their plastic ancestors (Kulkarni et al. 2011, 2017, Luquet et al. 2011). However, only few studies have empirically detected such costs (Van Buskirk and Steiner 2009, Auld et al. 2010). The associations between oxidative stress and immunological status and the degree of phenotypic plasticity here detected in spadefoot toad larvae are congruent with what is known from the physiological consequences of producing plastic shifts in amphibian larvae (i.e. production costs of plasticity). For instance, accelerated development in response to pond drying is achieved at the expense of increased oxidative stress (Gomez-Mestre et al. 2013, Burraco et al. 2017c) or reduced post-metamorphic immune function (Gervasi and Foufopoulus 2008). These known physiological production costs of plasticity are consistent with the fact that maintaining such plastic potentiality may cause small but perceptible physiological costs requiring buffering mechanisms. Maintenance costs of plasticity can thus be relevant to understanding differential survival within and among populations, if individuals with genotypes of high potential plasticity experience greater oxidative stress as observed at the sibship level. This can be particularly important in a context of rapid environmental change (Gunderson and Stillman 2015), which is expected to negatively impact fitness-related parameters in ectotherms as a consequence of induced developmental and growth plastic responses (Burraco et al. 2020). Therefore, a dismissal of maintenance costs as an important factor in the evolution of adaptive plasticity may be premature, in agreement with many theoretical studies (Chevin et al. 2010, Gomez-Mestre and Jovani 2013, Chevin and Lande 2015). A more direct evaluation of the physiological variation linked to the maintenance of adaptive plasticity across taxa is needed to re-evaluate the importance of such costs and buffering mechanisms in shaping the evolution of plasticity. Acknowledgements – We thank F. Miranda and O. García for assistance in the Ecophysiology laboratory. We also thank M. Rojas, L. Asencio and C. Pérez for assistance in animal husbandry. Funding – This study was funded by Ministerio de Economía y Competitividad (grants CGL2012-40044 and CGL201459206-P). PB was supported by fellowship F.P.U.-AP2010-5373 from the Ministerio de Educación and by a Marie Curie Fellowship METAGE-797879. Laboratory facilities were provided by ICTSRBD. Author contributions Pablo Burraco: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Writing – original draft (lead); Writing – review and editing (lead). Miguel Ángel Rendón: Data curation (supporting); Formal analysis (equal); Writing – review and editing (supporting). Carmen Díaz-Paniagua: Conceptualization (supporting); Investigation (supporting); Resources (supporting); Supervision (supporting); Writing – review and editing (supporting). Ivan Gomez-Mestre: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (supporting); Writing – review and editing (supporting). Data availability statement Data are available from the Dryad Digital Repository: <https://doi.org/10.5061/dryad.p2ngf1vqw> (Burraco et al. 2022). Supporting information The supporting information associated with this article is available from the online version.