Context dependent variation in corticosterone and phenotypic divergence of Rana arvalis populations along an acidification gradient

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Context dependent variation in corticosterone and phenotypic divergence of Rana arvalis populations along an acidification gradient © The Author(s) 2022 Published version Mausbach, Jelena; Laurila, Anssi; Räsänen, Katja Mausbach, J., Laurila, A., & Räsänen, K. (2022). Context dependent variation in corticosterone and phenotypic divergence of Rana arvalis populations along an acidification gradient. BMC Ecology and Evolution, 22, Article 11. https://doi.org/10.1186/s12862-022-01967-1 2022 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 https://doi.org/10.1186/s12862-022-01967-1 RESEARCH ARTICLE Context dependent variation incorticosterone andphenotypic divergence ofRana arvalis populations alonganacidification gradient Mausbach Jelena1,2*, Laurila Anssi3 and Räsänen Katja1,2,4* Abstract Background: Physiological processes, as immediate responses to the environment, are important mechanisms of phenotypic plasticity and can influence evolution at ecological time scales. In stressful environments, physiological stress responses of individuals are initiated and integrated via the release of hormones, such as corticosterone (CORT). In vertebrates, CORT influences energy metabolism and resource allocation to multiple fitness traits (e.g. growth and morphology) and can be an important mediator of rapid adaptation to environmental stress, such as acidification. The moor frog, Rana arvalis, shows adaptive divergence in larval life-histories and predator defense traits along an acidification gradient in Sweden. Here we take a first step to understanding the role of CORT in this adaptive divergence. We conducted a fully factorial laboratory experiment and reared tadpoles from three populations (one acidic, one neutral and one intermediate pH origin) in two pH treatments (Acid versus Neutral pH) from hatching to metamorphosis. We tested how the populations differ in tadpole CORT profiles and how CORT is associated with tadpole life-history and morphological traits. Results: We found clear differences among the populations in CORT profiles across different developmental stages, but only weak effects of pH treatment on CORT. Tadpoles from the acid origin population had, on average, lower CORT levels than tadpoles from the neutral origin population, and the intermediate pH origin population had intermediate CORT levels. Overall, tadpoles with higher CORT levels developed faster and had shorter and shallower tails, as well as shallower tail muscles. Conclusions: Our common garden results indicate among population divergence in CORT levels, likely reflecting acidification mediated divergent selection on tadpole physiology, concomitant to selection on larval life-histories and morphology. However, CORT levels were highly environmental context dependent. Jointly these results indicate a potential role for CORT as a mediator of multi-trait divergence along environmental stress gradients in natural populations. At the same time, the population level differences and high context dependency in CORT levels suggest that snapshot assessment of CORT in nature may not be reliable bioindicators of stress. Keywords: Acidification, Adaptive divergence, Amphibians, Corticosterone, Environmental stress, Evolutionary physiology, Phenotypic plasticity © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Open Access BMC Ecology and Evolution *Correspondence: [email protected]; [email protected] 1 Department of Aquatic Ecology, Eawag, Ueberlandstrasse 133, 8600 Duebendorf, Switzerland Full list of author information is available at the end of the article Page 2 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 Background Environmental change, be it natural or anthropogenic, is often associated with the exposure of individuals and populations to abiotic and biotic environmental stressors, which can lead to strong natural selection [1]. At short evolutionary time scales, environmental stress can lead to “rapid evolution” [2], raising the questions how natural selection acts on multiple interacting traits and what are the mechanisms of rapid adaptation? [1] In order to understand eco-evolutionary responses of populations to stressful and fast changing environments, we need to understand how environmental and genetic effects jointly act on the organismal phenotype [1, 3–5]. A major source of environmental responsiveness of organisms is physiological plasticity, which determines the immediate responses and the ability of individuals to acclimate to environmental stress [6–10]. Physiological responses are, hence, expected to be under natural selection [8, 11–13], whereby genotypes with optimal combinations of stress responses and energy metabolism for a given ecological context should be favoured [14, 15]. Although the role of physiology in adaptation has received attention in ecophysiology [16] and evolutionary physiology [17, 18], it has yet to be fully integrated across fields (i.e. as eco-evolutionary physiology of contemporary populations). In vertebrates, a candidate pathway in integrated stress responses arises via the glucocorticoid hormones corticosterone (CORT) and/or cortisol [19]. In addition to stress responses, these glucocorticoids are involved in general metabolic processes and a range of gene expression networks (metabolism, growth, tissue repair, reproduction, immune function; reviewed in [19, 20]), and can therefore be under strong natural selection. Consequently, glucocorticoids are major mediators of the phenotype and of special interest in the context of ecoevolutionary physiology [3, 19, 20]. CORT is secreted after the activation of the hypothalamic–pituitary–adrenal (HPA) axis which, among other functions, is one main physiological pathway responsible for stress responses in vertebrates [19, 20]. In the shortterm, elevated glucocorticoid levels can allow energy mobilization in stressful situations (e.g. via fat catabolism and decrease in digestion; [19]). However, chronically elevated CORT levels can be costly and cause reproductive malfunction, cellular damage and immunosuppression [19, 21–25]. Therefore, under long-term exposure to stress, populations face a trade-off: natural selection should prevent detrimental effects of elevated CORT levels, yet maintain the ability to respond adaptively to temporally varying stressors (e.g. predation attempts or extreme temperatures). In general, natural selection on CORT levels may act on both “supportive” (i.e. maintaining ability to respond by elevating CORT levels) and “protective” (i.e. reducing negative effects of elevated CORT) processes [14]. Thereby selection may favour either higher (supportive process) or lower (protective process) CORT levels in energetically challenging or stressful environments, and act on baseline and/or stress induced CORT levels in a context dependent manner (i.e. depending on costs versus benefits; [14]). In amphibians, environmental stress activates the hypothalamic-pituitary-interrenal (HPI) axis, leading to the secretion of CORT [26–29]. In tadpoles, CORT levels are influenced by many different stressors, such as acidity, predators and parasites [30–32], and CORT can influence many fitness traits, from growth rates and immune function to traits related to resource acquisition and predator defense [31]. CORT levels can also vary strongly across the developmental stages, generally peaking at metamorphosis [33]. At early to mid-larval stages elevated CORT levels may decrease growth and development rates [34–36], as well as reduce body length and increase tail depth [30, 31]. From late to mid-larval stages elevated CORT levels may instead accelerate development [26, 37]. CORT related performance trade-offs [31, 35, 36] and geographic variation in glucocorticoid levels has been demonstrated [14, 30, 34], indicating the potential for divergent natural selection through CORT. However, how CORT profiles and CORT—trait associations vary across divergent environments in natural populations is poorly known. Environmental acidity, both natural and anthropogenic, is stressful for a range of organisms [38, 39], including amphibians [reviewed in 40], and can be a potent agent of natural selection. Moor frog, Rana arvalis, populations along a pH gradient in Sweden show phenotypic divergence in multiple tadpole traits [41, 42]. Specifically, laboratory studies show that under common garden conditions, tadpoles from acid origin populations develop slower, grow faster and are larger at metamorphosis, and have deeper tails, than tadpoles from neutral origin populations [41, 43]. This divergence is mediated through a combination of maternal and direct genetic effects [43, 44] and is a response to both acidity and predator induced divergent selection [41]. However, the physiological underpinnings of this multi-trait divergence are unknown. Here we aim to increase the understanding of the role of glucocorticoids in adaptive divergence of natural populations. Specifically, we study the effects of acid stress on CORT levels, and associations between CORT and functionally relevant traits (life-history and morphology), in R. arvalis tadpoles from three divergent populations. In a fully factorial laboratory experiment, we reared tadpoles from an acidic (Tottatjärn, TT), a neutral (Rud, RD) Page 3 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 and an intermediate (Bergsjön, BS) pH origin population in two contrasting pH treatments: Acid pH (physiologically stressful) and Neutral pH (physiologically benign). We compared tadpoles from these population-treatment combinations at three developmental stages (mid-larval stages G32 and G38, and metamorphosis, G42, [45]). We made the following predictions. First, if acidic pH is stressful and CORT is an indicator of stress, tadpoles in the Acid pH treatment should show elevated CORT levels relative to the Neutral pH treatment. Second, if there has been divergent selection on either baseline (e.g. due to differential metabolic demands) or stress induced CORT [14], tadpoles from the acidic and neutral pH origin population should differ in their CORT profiles. This phenotypic divergence among the populations could be in the form of Genotype × Environment interactions (i.e. seen as differential CORT responses of populations to the pH treatments) and/or differences in mean CORT levels (independent of pH treatment). Finally, if CORT is a key mediator of multi-trait adaptive divergence, CORT levels should correlate with larval life-history traits and tail morphology (latter representing a typical inducible predator defence trait, [41]). In particular, we expected elevated CORT levels to correlate with larval developmental time, body size and tail morphology. Results Tadpoles from the three phenotypically divergent populations (TT: Acid origin, RD: Neutral origin, and BS: Intermediate origin) were reared in either Acid (target pH 4.3) or Neutral (target pH 7.5) pH in the lab. We tested i) how tadpoles from the three populations differ in CORT profiles and ii) what are the CORT—trait relationships in multivariate space for larval development time (days from G25 to a given tadpole stage), tadpole size (mass, g) and tadpole morphology (Fig. 1). With regard to the core hypotheses for CORT, Population main effects would be indicative of genetic divergence in response to selection on baseline CORT levels (involved in organismal metabolism in absence of stress) or chronic stress induced CORT levels [14], pH main effects and higher CORT in the Acid pH treatment would be indicative of stress induced CORT elevation, and Population x pH interaction effects would indicate among population differences in chronic CORT stress responses. Notably, given chronic exposure (weeks to months) of tadpoles, no difference between the benign (Neutral) and stressful (Acid) pH treatment may indicate that CORT levels have returned to baseline levels (e.g. to reduce detrimental effects of chronically elevated CORT) (see Discussion). Sampling and rearing were conducted in two Blocks (A: morning sampling/warmer temperature; B: afternoon sampling/cooler temperature) and measurements were taken at three larval stages: G32, G38 and G42 (see “Methods” for details). The intended number of replicates for each population—treatment combination was eight individuals, but the following population treatment combinations at G42 had N = 7 (TT4B and TT7B) and N = 9 (BS7B) (as detailed in Methods). Therefore, a total of 287 individuals were included in the statistical analyses. The data was analysed using univariate and multivariate AN(C)OVAs and interpreted based on differences in LS means (univariate analyses) and Hypothesis-Error (HE) plots (MANOVAs) where relevant (see Methods for details). Only final models are presented here. Corticosterone G32 and G38 tadpoles had, on average, lower CORT levels than the G42 metamorphs (Fig.2a–c, Additional file1: Table1.1A). However, univariate analyses of log(CORT) across developmental stages revealed a significant pH treatment × Block × Stage interaction (Additional file1: Table1.1A)—indicating that CORT variation was context dependent. To examine this three-way interaction Fig. 1 Morphological traits of Rana arvalis tadpoles measured at G32 and G38: body length (BL), body depth (BD), tail length (TL), maximum tail depth (TD) and tail muscle depth (TMD). BL was taken from mouth to the base of the hind leg, BD was taken where it was longest orthogonal to BL, TL was taken from base of the hind leg to tail tip, TD was measured where it is deepest and TMD was taken orthogonal to the “spine” right at its base Page 4 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 further, we next conducted models separately within each of the developmental stages (G32, G38 and G42). At G32, there were significant Population, Block and pH treatment × Block effects (Table1): Acid origin (TT) tadpoles had, on average, lower CORT levels than Neutral (RD) and Intermediate (BS) origin tadpoles (Tukey test, Table1, Fig.2a). Moreover, CORT levels were higher in the Acid than the Neutral pH treatment in the B block, whereas there was no difference between the pH treatments in the A block (Fig.2a, Tukey test Table1). At G38, only the Block effect was significant, with tadpoles in the A block having higher CORT levels than tadpoles in the B block (Fig.2b, Table1). At G42 no statistically significant effects on CORT were found (Fig.2c, Table1). Fig. 2 Mean ± SE of CORT levels (a–c), and LS means ± SE of log(developmental time, days from G25) (d–f) and log (mass in g) (g–i), across three larval developmental stages (G32, left panel, G38, middle panel, and G42, right panel) in Rana arvalis. Tadpoles from the Acid (TT), Neutral (RD) and Intermediate (BS) origin population were reared in Acid (4) or Neutral (7) pH treatment across two rearing Blocks (A: morning sampling/warmer; B: afternoon sampling/colder). Sample size was N = 8 except for the following cases: a) TT4B, TT7A and TT7B: N = 7; b) RD4B and BS4B: N = 6, BS4A and TT4B: N = 7; c) BS4B and TT4B: N = 6, TT7B: N = 7; f) TT4B and TT7B: N = 7, BS7B: N = 9; g) BS4B: N = 7; i) TT4B and TT7B: N = 7, BS7B: N = 9 (see Methods for details) Page 5 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 The multivariate phenotype CORT andlife history traits Mid‑larval stage G32 MANOVAs—Population, pH treatment and Block had significant main effects, but no significant interactive effects, on the joint variation between CORT, developmental time and tadpole mass at G32 (Additional file1: Table1.3). Block explained most of the variation in this multivariate space (eta2: 54%), followed by pH treatment (21%) and Population (18%) (Additional file1: Table1.3; for canonical HE analyses see Additional file2: Table2.1 and Fig.2.1A). On average, TT tadpoles developed slower and were larger than RD tadpoles, with BS tadpoles being intermediate—though life-history trait variation was pH treatment and Stage dependent (Fig.2d and g; for Univariate ANOVAs see Additional file1: Table1.4). HE-plots from the MANOVAs—At G32, there was a strong negative association between CORT and developmental time across Blocks (Fig. 3a, block ellipsoid): individuals with higher CORT levels (A block) developed faster than those with lower CORT levels (B block). There was also a negative association between CORT and developmental time across populations (Fig. 3a, pop ellipsoid): TT individuals had lower CORT levels and developed slower, whereas RD and BS individuals had higher CORT levels and developed faster. The HE plots further indicated that the pH treatment effect in the MANOVA (Additional file1: Table1.3) was primarily due to tadpoles developing slower in the Acid (4) than the Neutral (7) treatment (Fig.3c and f, pH ellipsoid), but there was no relationship between CORT and developmental time across pH treatments (Fig.3a). There was a subtle negative association between CORT and body mass across the Blocks (Fig.3b, block ellipsoid): tadpoles from the A block tended to have higher CORT levels but be smaller than those from the B block. There was a stronger negative association between CORT and tadpole mass across populations: individuals with lower CORT levels were larger (TT population) than individuals with higher CORT levels (RD and BS tadpoles) (Fig.3b, pop ellipsoid). There was no clear association between CORT and tadpole size in relation to the pH treatments. Finally, there was a subtle positive relationship between developmental time and mass of tadpoles across Blocks (Fig.3c): tadpoles in the A block developed faster (fewer days from G25 to G32) but were smaller, whilst tadpoles in the B block developed slower and were larger. In contrast, there was a strong positive association between developmental time and mass of tadpoles across Populations (Fig.3c, pop ellipsoid): TT tadpoles developed slower and were larger whereas RD and BS tadpoles developed faster and were smaller. Interestingly, the pH treatment (Fig.3c, pH ellipsoid) reversed this development time-mass Table 1 Univariate models of log(CORT) within G32, G38 and G42 stages Results of univariate linear models on log(CORT) for G32, G38 and G42, respectively. Results are presented for Rana arvalis tadpoles from the Acid (TT), Neutral (RD) and the Intermediate (BS) origin population, reared in Acid (4) or Neutral (7) pH treatment and two Blocks (A: morning/warmer, B: afternoon/colder). These are final models following removal of non-significant threeor two-way interactions. Statistically significant effects (p < 0.05) are shown in bold. Different letters in posthoc tests denote significantly different LS means (Tukey tests), indicating that 4A and 4B, and 7A and 7B, differ from each other (i.e. letters do not overlap), whereas there are no differences between the pH treatments within each block Factors Developmental stage G32 G38 G42 df F p df F p df F p Population 2 5.71 0.005 2 2.60 0.081 2 1.55 0.218 pH treatment 1 0.02 0.880 1 0.06 0.800 1 1.45 0.233 Block 1 15.15 < 0.001 1 152.51 < 0.001 1 2.10 0.151 Population × pH 2 0.34 0.712 2 1.36 0.262 2 0.52 0.596 pH × Block 1 4.84 0.031 – – – – – – Residual SE and (df) 0.46 (85) 0.22 (83) 0.32 (84) Posthoc (Tukey) LS mean ± SE Population: TT 1.21 ± 0.09 b BS 1.63 ± 0.08 a RD 1.72 ± 0.08 a Block: A 1.89 ± 0.07 a B 1.15 ± 0.07 b pH x Block: 4A 1.83 ± 0.09 a 7A 1.95 ± 0.10 a 4B 1.30 ± 0.10 b 7B 1.00 ± 0.10 b LS mean ± SE Block: A 2.09 ± 0.03 a B 1.50 ± 0.03 b Page 6 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 relationship with individuals that developed slower (i.e. Acid, 4, treatment) being smaller than those that developed faster (i.e. Neutral, 7, treatment)—reflecting stressful conditions in the acid treatment. Mid‑larval stage G38 MANOVAs—At G38, Block, Population and pH treatment had strong and significant main effects and a significant Population x pH interaction on joint variation of CORT, developmental time and mass (Additional file1: Fig. 3 HE plots from the CORT and life history MANOVAs at larval stage G32 (a–c upper panel), G38 (d–f middle panel) and G42 (g–i lower panel) in Rana arvalis. All response variables were log transformed. Hypothesis ellipses that are outside of the Error ellipse (in red) indicate significant effects. The ellipses depicted are pop = Population, ph = pH treatment, pop:ph = Population-pH treatment interaction, block = rearing block. The solid dots indicate fixed effect means for Population (TT: Acid origin, RD: Neutral origin, BS: Intermediate origin), pH treatment (4: Acid, 7: Neutral) and Block (A morning sampling/warmer, B afternoon sampling/colder) Page 7 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 Table1.3). Block explained most of the variation in this multivariate space (eta2: 66–67%), followed by pH treatment (46%), Population (17–18%) and Population x pH treatment (7%) (Additional file1: Table1.3). HE plots from the MANOVAs—There was a strong negative association between CORT and developmental time across Blocks at G38 (Fig.3d): individuals with high CORT levels (A block) developed faster than those with lower CORT levels (B block). Likewise, there was a negative association between CORT and body mass across the Blocks (Fig.3e, block ellipsoid): tadpoles from the A block had higher CORT levels and were smaller, whereas tadpoles from the B block fell to the opposite end of the axis. At G38, CORT was correlated with developmental time across Blocks (Fig.3d, block ellipsoid): tadpoles that had higher CORT levels (A block) developed faster than those that had lower levels (B block). There was a strong negative relationship between developmental time and mass of tadpoles in relation to pH (Fig.3f, pH ellipsoid): tadpoles developed slower but were smaller in the Acid pH treatment and developed faster and were larger in the Neutral pH treatment. Jointly with univariate analyses (Additional file1: Table1.4), the HE plots (Fig.3a–f, and below) indicated that the pH and population effects were primarily driven by effects on development time and mass. There was a significant association between development time and mass across populations (Fig.3f, pop ellipsoid), with TT tadpoles developing slower and being larger than RD tadpoles. Metamorphosis G42 MANOVA—At G42, only Population and Block had a significant main effect on the joint variation of CORT, developmental time and mass (Additional file1: Table1.3). Variance partitioning showed that Population explained 18–20% and Block 17% of this variation (Additional file1: Table1.3). TT metamorphs were substantially larger in both pH treatments than RD and BS metamorphs and individuals reached metamorphosis somewhat slower in the B block (Fig.2f and i; for univariate ANOVAs see Additional file1: Table1.4). HE plots from the MANOVAs As for G32 and G38, there was a negative association between CORT and developmental time at G42 across the populations and Blocks. The Population ellipsoid (Fig.3g and h) indicated that individuals with lower CORT levels (mostly TT tadpoles) were larger and developed slower than those with higher CORT levels (mostly BS and RD tadpoles). Individuals with higher CORT levels (i.e. A block) developed faster to metamorphosis than those with lower CORT levels (B block) (Fig.3g, block ellipsoid). There was a strong positive relationship between developmental time and mass across Populations (Fig.3i, pop ellipsoid): TT metamorphs developed slower and were larger whereas RD and BS metamorphs developed faster and were smaller. There was a weak positive association between developmental time and mass at G42 across Blocks (Fig.3i, block ellipsoid), with tadpoles that developed faster (A block) being smaller and tadpoles that developed slower being larger (B block). The patterns of developmental time and body mass trait means seen in the HE plots (Fig.3) for G32, G38 and G42 were confirmed in univariate models (Fig.2, Additional file1: Table1.1, Table1.4). CORT: morphology relationships Discriminant analyses of principal components (DAPC)— A multivariate DAPC including CORT and log transformed morphological traits (body length, body depth, tail length, tail depth and tail muscle depth) across G32 and G38 showed a clear phenotypic separation of the two developmental stages (Additional file3: Fig.3.1, 3.2A, B & Table3.1). As G32 is more reflective of functionally relevant mid-larval stage morphology and showed more variance across individuals in our data set (see Additional file3), we next conducted a separate DAPC on G32 tadpoles only. CORT: morphology relationships atG32 In the DAPC for G32 tadpoles, LD1 explained 64.8% of the variance, with tail muscle depth (64.6%), tail depth (18.4%) and CORT (12.2%) loading strongest on this axis. LD2 explained 18.2% of the variance, with tail length (37.9%) and body depth (29.3%) loading strongest (Fig.4, Additional file3: Fig.3.1, Fig.3.2C & Table3.1). Visual inspection indicated that LD1 reflects mostly population level variation, with TT tadpoles having lower CORT and relatively deeper tail muscles and tails, RD tadpoles having higher CORT and relatively shallower tail muscles and tails, and BS tadpoles being intermediate (Fig.4). LD2, on the other hand, reflected mostly variation related to the pH treatments, with tadpoles in the Acid treatment having a shallower body and longer tail and tadpoles in the Neutral treatment having a relatively deeper body and longer tail (Fig.4). MANOVAs—To investigate the multivariate relationships between CORT and tadpole morphology at G32, we next conducted a MAN(C)OVA with CORT and body length, body depth, tail length, tail muscle depth and tail depth as response variables, and mass as a covariate. All traits were log transformed. This analysis found significant Population, pH treatment, Block and mass main effects, but no significant Population x pH interaction effect (Additional file1: Table1.3). Partitioning of Page 8 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 variance indicated that strongest effects on the multivariate phenotype were by mass (94%), Block (45%), Population (28–31%) and pH treatment (21%) (Additional file1: Table1.3) (This ranking held for Pillai’s, Wilk’s Lambda as well as Hotelling Lawley’s test statistics. Details of the canonical analyses and HEplots can be found in Additional file 2. For LSmeans of individual morphological traits, and univariate ANOVAs see Additional file1: Fig.1.1 and Additional file1: Table1.4, respectively). Based on the HE plots, there was a subtle negative relationship between CORT and tail length across blocks: tadpoles with higher CORT (A block) had shorter tails than tadpoles with lower CORT (B block) (Fig. 5 and Additional file2: Fig.2.3). There was a subtle negative relationship between CORT and both tail length and tail depth: tadpoles with higher CORT (RD) had shorter and shallower tails. The relationship between CORT and tail muscle depth was strongly associated to population: tadpoles with higher CORT levels (RD) had shallower tail muscles compared to those with lower CORT levels (TT). No other morphological traits were associated with variation in CORT. In summary, we found among population divergence in multivariate space for CORT, life-history traits and morphology of R. arvalis tadpoles (see Fig.6). Specifically, tadpoles from an acid origin population (TT) had lower CORT levels, developed slower and were larger, and had relatively deeper tails and tail muscles during the mid-larval stage G32 relative to the neutral and intermediate origin population (RD and BS). Many effects were strongly affected by block effects, indicating strong phenotypic plasticity. Most intriguingly, we found negative associations between CORT levels and developmental time (at G32, G38 and G42), as well as CORT and tadpole tail length (at G32), across the rearing blocks: tadpoles from the A block had higher CORT levels and developed faster across all three developmental stages, they were also larger but had a shorter tail than those in the B block. Discussion We found clear multivariate phenotypic divergence among tadpoles of three R. arvalis populations that were reared in acid versus neutral pH in the lab. In accordance with our previous studies [42], tadpoles from the acidic (TT) and the neutral pH origin population (RD) (two ends of an acidification gradient [42],) were more divergent, and tadpoles from the intermediate pH origin (BS) population were more variable or intermediate in their phenotype. Most intriguingly, we found among population divergence in CORT levels: TT tadpoles had, on average, lower CORT levels (especially at mid-larval stage G32) than RD tadpoles, with BS tadpoles being intermediate. Variation in CORT was, however, highly context dependent. As expected, CORT levels at metamorphosis (G42) were much higher than during mid-larval stages (G32 and G38). However, the effects of pH treatment on CORT were weak: only G32 tadpoles in one of the two rearing Blocks (B block) showed higher CORT levels in the Acid than the Neutral treatment. In contrast, block effects were generally strong both on CORT and other tadpole traits, and can reflect variation due to circadian rhythm and/or temperature (discussed below). Finally, our multivariate analyses of CORT—trait associations showed that higher CORT levels were related to faster development (at all stages) and to relatively shorter and shallower tails and shallower tail muscles (at G32). Corticosterone The most striking and novel finding in our study is divergence in CORT profiles between the three R. arvalis populations. Given common garden rearing (and individuals originating from multiple families within each population), these results suggest genetic divergence in CORT levels, although a contribution of maternal effects is also possible [e.g. 44, 46, see below]. The lower CORT levels of TT tadpoles indicate that acidity mediated selection may have favoured the general downregulation of CORT to reduce CORT induced costs under chronic stress [19, 21–25], as a so called “protective” mechanism [14]. Alternatively, it is possible that selection has favoured lower baseline CORT levels through selection on other traits, Fig. 4 DAPC on stage G32 tadpoles for CORT and morphological traits (body length, body depth, tail length, tail depth, tail muscle depth) in Rana arvalis. Tadpoles from an Acid origin (TT), Intermediate origin (BS) and a Neutral origin (RD) population were reared in an Acid (here: 4) and a Neutral (here: 7) pH treatment and two Blocks (A and B). All response variables were log transformed for this analysis. LD1 represented mostly variation in tail depth (TD), tail muscle depth (TMD) and CORT, LD2 represents mainly tail length (TL) and body depth (BD) (See main text and Additional file 3 for details) Page 15 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 reach normality (where relevant) and to assure all traits were on similar scales in multivariate analyses. Data was checked for normality before and after log transformation. In all linear models, normality was visually assessed using QQ plots and by checking the distribution of residuals. All data sets were analysed by J. Mausbach and had information on population treatment combination. In a few cases, individuals appeared as outliers in some traits (one extreme value in CORT and TD). These individuals were retained in the analyses because they did not influence statistical significance and were more likely to represent biologically extreme values than measurement error or other unwanted variation. Univariate linear models Data was analysed using linear models (ANOVA Type III in nlme package using options(contrast = c("contr. sum","contr.poly")) and drop1(model,. ~ .,test = "F") [96] or the “car” package [97]). The full models included fixed factors of Population (3 levels), pH treatment (2 levels), Developmental stage (3 levels), Block (2 levels) and relevant twoto four-way interactions. We always started with a full model containing all fixed effect interactions. To then reduce model complexity, we sequentially removed non-significant interactions (P ≥ 0.05, starting with four-way interactions) using backwards selection of the linear models, with always maintaining all main effects and experimentally meaningful interactions in the model. For example, the Population x pH treatment interaction was always retained as this tested a key hypothesis and was an essential part of the fully factorial study design. Next, CORT was analysed within each of the developmental stages (G32, G38 and G42) separately. In these models, main effects of Population, pH, Block and their interactions were included. Where relevant and not statistically confounded, tadpole mass was included as a covariate. Non-significant covariate interactions and covariate main effects were sequentially removed and only final models are reported here. As the three study populations differ substantially in larval mass (Table1 and Additional file1: Table1.3), and population and mass would be statistically confounded, mass was not added as a covariate in final univariate statistical models on CORT. Instead, to test whether tadpole size affected CORT levels models with log(mass) as covariate were run within each study population (Additional file 1: Table 1.2). This ANCOVA was done only within the G32 stage. We report Means ± SE of the data and tests for relevant pairwise differences in LSmeans using post hoc Tukey tests (R package: lsmean, FSA, ggplot2). Multivariate phenotype CORT andLife history In order to assess covariance between CORT, developmental time (days from G25 until G32, G38 or G42) and tadpole mass, MANOVAs were run within G32, G38 and G42 stages. These analyses included log(CORT), log(developmental time) and log(mass) as response variables. Fixed factors of Population (3 levels), pH treatment (2 levels) and Block (2 levels), and their twoto threeway interactions, were included as predictors. We always started with a full model containing all fixed effect interactions and reduced model complexity by sequentially removing non-significant interactions (P ≥ 0.05, starting with four-way interactions) using backwards selection. Only final models are reported here. For final models, we calculated the partial variance eta2 [98], using the heplot package in R [99, 100]. We report Wilk’s, Pillai’s and Hotteling Lawley’s eta2 with ranking of all partial variances. However, we only present test statistics for Wilk’s tests (MANOVA type III, contrasts = list(topic = contr.sum, sys = contr.sum)). For visual presentation of MANOVAs we used the package heplot (MANOVA type III) in R [99, 100] that plots ellipsoids of hypotheses (H) and the error (E) of a given model. Significance is indicated by H ellipsoids reaching out of the E ellipsoid [99, 100]. The HE plots were run both on canonical models (candisc package [101]), to visualize the overall MANOVA results in one plot, and as simple HE plots derived from the MANOVA. As we were specifically interested in the mid-larval stage G32, the canonical analyses were conducted only for G32 tadpoles. All MANOVAs were followed by univariate linear models (ANOVA type III, using options(contrast = c("contr. sum","contr.poly")) and drop1(model,. ~ .,test = "F")). We report means ± SE of the univariate data and test for relevant pairwise differences of LSmeans using post hoc Tukey tests (R package: lsmean, FSA, ggplot2 [102–104]). Morphology andCORT relationship Visual representation morphology and CORT—To assess covariance visually across the full phenotype at mid-larval stages G32 and G38, a Discriminant Analysis of Principle Components (DAPC) was conducted (R adegenet package 2.0, [105, 106]). In this multivariate statistical approach, the variance of the data is partitioned into a between-group and within-group component, in order to maximize the discrimination between the groups [105]. The data was first analysed using a principal component analysis (PCA) and afterwards clusters identified using discriminant analysis (DA) [105]. In our data set, combinations of developmental stage (G32 and G38), Population (TT, BS, RD), pH treatment (Acid and Neutral) and Block (A and B) were used to Page 16 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 assign individuals into groups. These analyses were conducted sequentially for two DAPCs: I. The 1st DAPC included log CORT, tail length, tail depth, tail muscle depth, body depth and body length at G32 and G38. II. The 2nd DAPC included log CORT, tail length, tail depth, tail muscle depth, body depth and body length at G32. For all DAPCs, contributions of the Loadings (e.g. LD1, LD2) were calculated by dividing the respective LD through the sum of “eigenvalues”. Only LD1 and LD2 were considered further as the data were sufficiently described by these two (> 80%) and all variable contributions that were above 10% for those LDs are reported. An ordination and loading plot was used to illustrate the grouping and most relevant contribution of variables. Multivariate tadpole morphology and CORT at G32—In order to assess the covariation between CORT and multivariate morphology (body depth, body length, tail length, tail depth, tail muscle depth), a MANCOVA was run at G32. In these models all traits were log transformed. Fixed effects of Population (3 levels), pH treatment (2 levels) and Block (2 levels) as main effects, log(mass) as a covariate and relevant two to four-way interactions were included as explanatory variables. Model selection, partitioning of variance, and univariate model testing and visualization, was conducted as for the CORT—life history MANOVA (see above). Visual representation of the MANCOVA conducted with HE plots using the heplot package and manual [99, 107]. Abbreviations CORT: Corticosterone; EIA: Enzyme Immuno Assay; TT: Tottatjärn (acid origin population, breeding pond pH 4.0); BS: Bergsjön (intermediate origin population, breeding pond pH 6.1); RD: Rud (neutral origin population, breeding pond pH 7.0); 4 Acid: pH treatment (target pH 4.3); 7 Neutral: pH treatment (target pH 7.5); BL: Body length; BD: Body depth; TL: Tail length; TD: Maximum tail depth; TMD: Tail muscle depth; RSW: Reconstituted soft water; temp.: Temperature; dev. Time: Developmental time (days from G25); treatm.: Treatment. Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s1286202201967-1. Additional file1: Additional Results for univariate and multivariate AN(C) OVAs Additional file2: Canonical models and additional HE plots of MAN(C) OVAs Additional file3: Additional DAPC results Additional file4: Validating the hormonal sampling and assay methods Acknowledgements We thank Maike Demski, Robin Haglund, Frida Sjösten, Nadine Tardent, Tanja Trüb, Mehdi Khadraoui, Catia Chavez and Gunilla Engström for help during field and/or laboratory work and N. Tardent and Nora Weissert for morphology measurement on ImageJ. We thank the land owners of our study sites for access to their land for field work. We thank Tobias Deschner and Roisin Murtagh for conducting the LCMS validations for Corticosterone for our study at the Endocrinology Laboratory at Max Planck Institute for Evolutionary Anthropology in Leipzig and Christoffer Bergvall at Uppsala University and Johan Meyr at the Swedish Agricultural University for lending out laboratory equipment. We thank Pablo Buracco and Ivan Gomez-Mestre for feedback on hormonal extraction methods and EIA analysis and Wolfgang Goymann for advice on the endocrinological part of this study. Additionally, we thank Wolfgang Goymann and Jukka Jokela for commenting on the concepts and analysis, and anonymous reviewers as well as Maren Vitousek for constructive comments on earlier versions of this manuscript. Authors’ contributions JM and KR jointly designed the study. JM performed the experiments, analysed the data and wrote the first draft of the manuscript. KR acquired the SNF grant for this study, assisted during data collection, guided data analyses and commented and refined earlier versions of this manuscript. AL gave input to the design of the experiment, was the host at the facilities to carry out the experiment, gave feedback during analysis and commented and refined earlier versions of this manuscript. All authors read and approved the final manuscript. Funding This study was financed by grants from Swiss National Science foundation (SNF) (to KR, Number: 31003A_166201). The SNF has no role in the design, analysis, or reporting of the study but warrants scientific funding based on a rigorous peer review process. Availability of data and materials The datasets generated and/or analysed during the current study are available upon publication at Dryad. https:// doi. org/ 10. 5061/ dryad. tx95x 69zp. Declarations Ethics approval and consent to participate The experiments were conducted under permissions from the Ethical committee for animal experiments in Uppsala County (“Uppsala djurförsöksetiska nämnd”, Uppsala Animal Experimental Ethics Board), https:// djur. jordb ruksv erket. se/ amnes omrad en/ djur/ olika slags djur/ forso ksdjur/ etisk tgodk annan deavd jurfo rsok.4. 78507 16f11 cd786 b52d8 00021 46. html (permit number 5.8.18–01518/2017). The County board of Västra Götaland permitted the egg collections for the experiments (permit number 522–6251-2017). Consent for publication Not applicable. Competing interests The authors declare that they have no competing interests. Author details 1 Department of Aquatic Ecology, Eawag, Ueberlandstrasse 133, 8600 Duebendorf, Switzerland. 2 Institute of Integrative Biology, ETH Zurich, Universitätstrasse 16, 8092 Zurich, Switzerland. 3 Animal Ecology/Department of Ecology and Genetics, Evolutionary Biology Centre, Uppsala University, Norbyvägen 18D, 75236 Uppsala, Sweden. 4 Department of Biological and Environmental Science, University of Jyväskylä, Survontie 9C, 40014 Jyväskylä, Finland. Received: 29 August 2020 Accepted: 26 January 2022 References 1. Hoffmann AA, Parsons PA. Extreme environmental change and evolution. Cambridge: Cambridge University Press; 1997. Page 17 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 2. Bijlsma R, Loeschcke V. Environmental stress, adaptation and evolution: an overview. J Evol Biol. 2005;18:744–9. 3. Badyaev AV. Stress-induced variation in evolution: from behavioural plasticity to genetic assimilation. Proc Biol Sci. 2005;272:877–86. 4. Houle D, Govindaraju DR, Omholt S. Phenomics: the next challenge. Nat Rev Genet. 2010;11:855–66. 5. Pigliucci M. Phenotypic plasticity: beyond nature and nurture. Baltimore: Johns Hopkins University Press; 2001. 6. Bonier F, Martin PR. How can we estimate natural selection on endocrine traits? Lessons from evolutionary biology. Proc Biol Sci. 2016;283:20161887. 7. Denver RJ. Environmental stress as a developmental cue: Corticotropinreleasing hormone is a proximate mediator of adaptive phenotypic plasticity in amphibian metamorphosis. Horm Behav. 1997;31:169–79. 8. Schoenle LA, Zimmer C, Vitousek MN. Understanding context dependence in glucocorticoid–fitness Relationships: The role of the nature of the challenge, the intensity and frequency of stressors, and life history. Integr Comp Biol. 2016;58:777–89. 9. Taff CC, Vitousek MN. Endocrine flexibility: optimizing phenotypes in a dynamic world? Trends Ecol Evol. 2016;31:476–88. 10. Zimmer C, Taff CC, Ardia DR, Ryan TA, Winkler DW, Vitousek MN. On again, off again: acute stress response and negative feedback together predict resilience to experimental challenges. Funct Ecol. 2019;33:619–28. 11. Bonier F, Martin PR, Moore IT, Wingfield JC. Do baseline glucocorticoids predict fitness? Trends Ecol Evol. 2009;24:634–42. 12. Mormede P, Terenina E. Molecular genetics of the adrenocortical axis and breeding for robustness. Domest Anim Endocrinol. 2012;43:116–31. 13. Stedman JM, Hallinger KK, Winkler DW, Vitousek MN. Heritable variation in circulating glucocorticoids and endocrine flexibility in a free-living songbird. J Evol Biol. 2017;30:1724–35. 14. Vitousek MN, Johnson MA, Downs CJ, Miller ET, Martin LB, Francis CD, Donald JW, Fuxjager MJ, Goymann W, Hau M, Husak JF, Kircher BK, Knapp R, Schoenle LA, Williams TD. Macroevolutionary patterning in glucocorticoids suggests different selective pressures shape baseline and stress-induced levels. Am Nat. 2019;193:866–80. 15. Romero LM, Beattie UK. Common myths of glucocorticoid function in ecology and conservation. J Exp Zool A Ecol Integr Physiol. 2021 16. Blaustein AR, Gervasi SS, Johnson PT, Hoverman JT, Belden LK, Bradley PW, Xie GY. Ecophysiology meets conservation: understanding the role of disease in amphibian population declines. Philos Trans R Soc Lond B Biol Sci. 2012;367:1688–707. 17. Feder ME, Bennett AF, Huey RB. Evolutionary physiology. Annu Rev Ecol Syst. 2000;31:315–41. 18. Storz JF, Bridgham JT, Kelly SA, Garland T. Genetic approaches in comparative and evolutionary physiology. Am J Physiol Regul Integr Comp Physiol. 2015;309:R197–214. 19. Sapolsky RM, Romero LM, Munck AU. How do glucocorticoids influence stress responses? Integrating permissive, suppressive, stimulatory, and preparative actions. Endocr Rev. 2000;21:55–89. 20. Guillette LJ, Crain DA, Rooney AA, Pickford DB. Organization versus activation—the role of endocrine disrupting contaminants (EDCS) during embryonic development in the wild. Environ Health Perspect. 1995;103:157–64. 21. Hodges K, Brown J, Heistermann M. Endocrine Monitoring of Reproduction and Stress. In: Thompsons KV, Kirk Baer C, editors. Kleiman DG. Wild mammals in captivity. Principles and Techniques for Zoo Management. Chicago: The University of Chicago Press; 2010. p. 447–68. 22. Kindermann C, Narayan EJ, Hero JM. Urinary corticosterone metabolites and chytridiomycosis disease prevalence in a free-living population of male Stony Creek frogs (Litoria wilcoxii). Comp Biochem Physiol A Mol Integr Physiol. 2012;162:171–6. 23. McEwen BS, Wingfield JC. The concept of allostasis in biology and biomedicine. Horm Behav. 2003;43:2–15. 24. Monfort S. Non-invasive endocrine measures of reproduction and stress in wild populations. In: Holt W, Pickard A, Rodger J, Wildt D, editors. Reproductive science and integrated conservation. Cambridge: Cambridge University Press; 2002. p. 147–65. 25. Sands J, Creel S. Social dominance, aggression and faecal glucocorticoid levels in a wild population of wolves Canis lupus. Anim Behav. 2004;67:387–96. 26. Denver RJ. Stress hormones mediate environment-genotype interactions during amphibian development. Gen Comp Endocrinol. 2009;164:20–31. 27. Hill R, Wyse G, Anderson M. Animal physiology. Sunderland: Sinauer Associate; 2008. 28. McEwen BS, Wingfield JC. Allostasis and allostatic load. In: Fink G, editor. Encyclopedia of stress. New York: Academic Press; 2007. 29. Rollins-Smith LA. Neuroendocrine-immune system interactions in amphibians - implications for understanding global amphibian declines. Immunol Res. 2001;23:273–80. 30. Chambers DL, Wojdak JM, Du P, Belden LK. Pond acidification may explain differences in corticosterone among salamander populations. Physiol Biochem Zool. 2013;86:224–32. 31. Middlemis Maher JM, Werner EE, Denver RJ. Stress hormones mediate predator-induced phenotypic plasticity in amphibian tadpoles. Proc R Soc Lond B Biol Sci. 2013;280:20123075. 32. Marino JA, Holland MP, Middlemis Maher JM. Predators and trematode parasites jointly affect larval anuran functional traits and corticosterone levels. Oikos. 2014;123:451–60. 33. Glennemeier KA, Denver RJ. Developmental changes in interrenal responsiveness in anuran amphibians. Integr Comp Biol. 2002;42:565–73. 34. Dahl E, Orizaola G, Winberg S, Laurila A. Geographic variation in corticosterone response to chronic predator stress in tadpoles. J Evol Biol. 2012;25:1066–76. 35. Glennemeier KA, Denver RJ. Role for corticoids in mediating the response of Rana pipiens tadpoles to intraspecific competition. J Exp Zool. 2002;292:32–40. 36. Glennemeier KA, Denver RJ. Small changes in whole-body corticosterone content affect larval Rana pipiens fitness components. Gen Comp Endocrinol. 2002;127:16–25. 37. Denver RJ, Boorse GC, Glennemeier KA. Endocrinology of complex life cycles: Amphibians. In: Pfaff D, Arnold A, Etgen A, Fahrbach S, Moss R, Rubin R, editors. Hormones, brain and behaviour, vol. 2. San Diego: Academic Press; 2002. p. 469–513. 38. Schindler DW, Mills KH, Malley DF, Findlay DL, Shearer JA, Davies IJ, Turner MA, Linsey GA, Cruikshank DR. Long term ecosystem stress— the effects of years of experimental acidification on a small lake. Science. 1985;228:1395–401. 39. Petrin Z, Englund G, Malmqvist B. Contrasting effects of anthropogenic and natural acidity in streams: a meta-analysis. Proc Biol Sci. 2008;275:1143–8. 40. Räsänen K, Green D. Acidification and its effects on amphibian populations. In: Heatwole H, Wilkinson JW, editors. Amphibian Decline. Diseases, Parasites, Maladies and Pollution. Amphibian Biology. Chipping Norton, Australia: Surrey Beatty and Sons; 2009. p. 3244–67. 41. Egea-Serrano A, Hangartner S, Laurila A, Räsänen K. Multifarious selection through environmental change: acidity and predatormediated adaptive divergence in the moor frog (Rana arvalis). Proc Biol Sci. 2014;281:20133266. 42. Hangartner S, Laurila A, Räsänen K. Adaptive divergence of the moor frog (Rana arvalis) along an acidification gradient. BMC Evol Biol. 2011;11:366. 43. Hangartner S, Laurila A, Räsänen K. Adaptive divergence in moor frog (Rana arvalis) populations along an acidification gradient: inferences from QST–FST correlations. Evolution. 2012;66:867–81. 44. Hangartner S, Laurila A, Räsänen K. The quantitative genetic basis of adaptive divergence in the moor frog (Rana arvalis) and its implications for gene flow. J Evol Biol. 2012;25:1587–99. 45. Gosner KL. A simplified table for staging anuran embryos and larvae with notes on identification. Herpetologica. 1960;16:183–90. 46. Räsänen K, Laurila A, Merilä J. Maternal investment in egg size: environmentand population-specific effects on offspring performance. Oecologia. 2005;142:546–53. 47. Haase CG, Long AK, Gillooly JF. Energetics of stress: linking plasma cortisol levels to metabolic rate in mammals. Biol Lett. 2016;12:20150867. 48. Bókony V, Lendvai AZ, Liker A, Angelier F, Wingfield JC, Chastel O. Stress response and the value of reproduction: are birds prudent parents? Am Nat. 2009;173:589–98. Page 18 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 49. Hau M, Ricklefs RE, Wikelski M, Lee KA, Brawn JD. Corticosterone, testosterone and life-history strategies of birds. Proc Biol Sci. 2010;277:3203–12. 50. Kulkarni SS, Denver RJ, Gomez-Mestre I, Buchholz DR. Genetic accommodation via modified endocrine signalling explains phenotypic divergence among spadefoot toad species. Nat Commun. 2017;8:993. 51. Reeder DM, Kramer KM. Stress in free-ranging mammals: integrating physiology, ecology, and natural history. J Mammal. 2005;86:225–35. 52. Vitousek MN, Johnson MA, Donald JW, Francis CD, Fuxjager MJ, Goymann W, Hau M, Husak JK, Kircher BK, Knapp R, Martin LB, Miller ET, Schoenle LA, Uehling JJ, Williams TD. HormoneBase, a populationlevel database of steroid hormone levels across vertebrates. Sci Data. 2018a;5:180097. 53. Schjolden J, Backström T, Pulman KG, Pottinger TG, Winberg S. Divergence in behavioural responses to stress in two strains of rainbow trout (Oncorhynchus mykiss) with contrasting stress responsiveness. Horm Behav. 2005;48:537–44. 54. Trenzado CE, Carrick TR, Pottinger TG. Divergence of endocrine and metabolic responses to stress in two rainbow trout lines selected for differing cortisol responsiveness to stress. Gen Comp Endocrinol. 2003;133:332–40. 55. Atwell JW, Cardoso GC, Whittaker DJ, Campbell-Nelson S, Robertson KW, Ketterson ED. Boldness behavior and stress physiology in a novel urban environment suggest rapid correlated evolutionary adaptation. Behav Ecol. 2012;23:960–9. 56. Béziers P, San-Jose LM, Almasi B, Jenni L, Roulin A. Baseline and stressinduced corticosterone levels are heritable and genetically correlated in a barn owl population. Heredity (Edinb). 2019;123:337–48. 57. Narayan, EJ. Non-invasive reproductive and stress endocrinology in amphibian conservation physiology. Conserv Physiol. 2013;1. 58. Wikelski M, Cooke SJ. Conservation Physiology. Trends Ecol Evol. 2006;21:38–46. 59. Narayan EJ, Forsburg ZR, Davis DR, Gabor CR. Non-invasive Methods for Measuring and Monitoring Stress Physiology in Imperiled Amphibians. Front Ecol Evol. 2019;7:431. 60. MacColl AD. The ecological causes of evolution. Trends Ecol Evol. 2011;26:514–22. 61. Briscoe Runquist RD, Gorton AJ, Yoder JB, Deacon NJ, Grossman JJ, Kothari S, Lyons MP, Sheth SN, Tiffin P, Moeller DA. Context dependence of local adaptation to abiotic and biotic environments: a quantitative and qualitative synthesis. Am Nat. 2020;195:412–31. 62. Bachmann JC, Van Buskirk J. Adaptation to elevation but limited local adaptation in an amphibian. Evolution. 2021;75:14109. 63. Sachs LM, Buchholz DR. Insufficiency of Thyroid Hormone in Frog Metamorphosis and the Role of Glucocorticoids. Front Endocrinol (Lausanne). 2019;10:287. 64. Pfister HP, King MG. Adaptation of the glucocorticosterone response to novelty. Physiol Behav. 1976;17:43–6. 65. D’Agostino J, Vaeth GF, Henning SJ. Diurnal rhythm of total and free concentrations of serum corticosterone in the rat. Acta Endocrinol (Copenh). 1982;100:85–90. 66. Coe CL, Levine S. Diurnal and annual variation of adrenocortical activity in the squirrel monkey. Am J Primatol. 1995;35:283–92. 67. Dupont W, Bourgeois P, Reinberg A, Vaillant R. Circannual and circadian rhythms in the concentration of corticosterone in the plasma of the edible frog (Rana esculenta L.). J Endocrinol. 1979;80:117–25. 68. Pancak MK, Taylor DH. Seasonal and daily plasma-corticosterone rhythms in American toads Bufo americanus. Gen Comp Endocrinol. 1983;50:490–7. 69. McEwen BS, Brinton RE, Sapolsky RM. Glucocorticoid Receptors and Behavior: Implications for the Stress Response. In: Chrousos GP, Loriaux DL, Gold PW, editors. Mechanisms of Physical and Emotional Stress. Adv Exp Med Biol. Boston, MA:Springer;1988. p.35–45. 70. Gangloff EJ, Holden KG, Telemeco RS, Baumgard LH, Bronikowski AM. Hormonal and metabolic responses to upper temperature extremes in divergent life-history ecotypes of a garter snake. J Exp Biol. 2016;219:2944–54. 71. McEwen BS, Wingfield JC. What is in a name? Integrating homeostasis, allostasis and stress. Horm Behav. 2010;57:105–11. 72. Romero LM, Dickens MJ, Cyr NE. The Reactive Scope Model - a new model integrating homeostasis, allostasis, and stress. Horm Behav. 2009;55:375–89. 73. Florencio M, Burraco P, Rendón MA, Díaz-Paniagua C, Gomez-Mestre I. Opposite and synergistic physiological responses to water acidity and predator cues in spadefoot toad tadpoles. Comp Biochem Physiol A Mol Integr Physiol. 2020;242:10654. 74. Rich EL, Romero LM. Exposure to chronic stress downregulates corticosterone responses to acute stressors. Am J Physiol Regul Integr Comp Physiol. 2005;288:R1628–36. 75. Mizoguchi K, Yuzurihara M, Ishige A, Sasaki H, Chui DH, Tabira T. Chronic stress differentially regulates glucocorticoid negative feedback response in rats. Psychoneuroendocrinology. 2001;26:443–59. 76. Lema SC. Hormones, developmental plasticity, and adaptive evolution: endocrine flexibility as a catalyst for ‘plasticity-first’ phenotypic divergence. Mol Cell Endocrinol. 2020;502:110678. 77. Vitousek MN, Taff CC, Hallingen KK, Zimmer C, Winkler DW. Hormone and Fitness: Evidence for Trade-Offs in Glucocorticoid Regulation Across Contexts. Front Ecol Evol. 2018. https:// doi. org/ 10. 3389/ fevo. 2018. 00042. 78. Crespi EJ, Williams TD, Jessop TS, Delehanty B. Life history and the ecology of stress: how do glucocorticoid hormones influence life-history variation in animals? Funct Ecol. 2013;27:93–106. 79. Harkey GA, Semlitsch RD. Effects of Temperature on Growth, Development, and Color Polymorphism in the Ornate Chorus Frog Pseudacris Ornata. Copeia. 1988;4:1001–7. 80. Crossin GT, Love OP, Cooke SJ, Williams TD. Glucocorticoid manipulations in free-living animals: considerations of dose delivery, life-history context and reproductive state. Funct Ecol. 2016;30:116–25. 81. Denver RJ. Proximate mechanisms of phenotypic plasticity in amphibian metamorphosis. Am Zool. 1997;37:172–84. 82. Denver RJ, Middlemis-Maher J. Lessons from evolution: developmental plasticity in vertebrates with complex life cycles. J Dev Orig Health Dis. 2010;1:282–91. 83. Altwegg R, Reyer HU. Patterns of natural selection on size at metamorphosis in water frogs. Evolution. 2003;57:872–82. 84. Husak JF. Measuring Selection on Physiology in the Wild and Manipulating Phenotypes (in Terrestrial Nonhuman Vertebrates). Compr Physiol. 2015;6:63–85. 85. Glandt D, Jahle R. Der Moorfrosch/The Moor Frog. Germany, Bielefeld: Laurenti Verlag; 2008. 86. Räsänen K, Soderman F, Laurila A, Merilä J. Geographic variation in maternal investment: Acidity affects egg size and fecundity in Rana arvalis. Ecology. 2008;89:2553–62. 87. APHA: Standard methods for the examination of water and wastewater. 1985, American Public Health Association, Washington, DC, 16 88. Alvarez D, Nicieza AG. Effects of temperature and food quality on anuran larval growth and metamorphosis. Funct Ecol. 2002;16:640–8. 89. CarmonaOsalde C, OlveraNovoa MA, RodriguezSerna M, FloresNava A. Estimation of the protein requirement for bullfrog (Rana catesbeiana) tadpoles, and its effect on metamorphosis ratio. Aquac. 1996;141:223–31. 90. Ramlochansingh C, Branoner F, Chagnaud BP, Straka H. Efficacy of tricaine methanesulfonate (MS-222) as an anesthetic agent for blocking sensory-motor responses in Xenopus laevis tadpoles. PLoS One. 2014;9:e101606. 91. Cakir Y, Strauch SM. Tricaine (MS-222) is a safe anesthetic compound compared to benzocaine and pentobarbital to induce anesthesia in leopard frogs (Rana pipiens). Pharmacol Rep. 2005;57:467–74. 92. Smith MS, Booth NJ, Peterson BC, Stephens WS, Goudie CA, Simco BA. Analysis of Short-Term Cortisol Stress Response in Channel Catfish by Anesthetization with Metomidate Hydrochloride and Tricaine Methanesulfonate. J Aquat Anim Health. 2015;27:152–5. 93. Hernández SE, Sernia C, Bradley AJ. The effect of three anaesthetic protocols on the stress response in cane toads (Rhinella marina). Vet Anaesth Analg. 2012;39:584–90. 94. Archard G and Goldsmith AR. Euthanasia methods, corticosterone and haematocrit levels in Xenopus laevis: Evidence for differences in stress?. Anim Welf. 2010;19. Page 19 of 19 Jelenaetal. BMC Ecology and Evolution (2022) 22:11 • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year • At BMC, research is always in progress. Learn more biomedcentral.com/submissions Ready to submit your research Ready to submit your research ? Choose BMC and benefit from: ? Choose BMC and benefit from: 95. Burraco P, Arribas R, Kulkarni SS, Buchholz DR, Gomez-Mestre I. Comparing techniques for measuring corticosterone in tadpoles. Curr Zool. 2015;61:835–45. 96. Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1–148; 2020. 97. Fox J, Weisberg S. An R Companion to Applied Regression. 3rd ed. Thousand Oaks, CA: Sage; 2019. 98. Langerhans RB, DeWitt TJ. Shared and unique features of evolutionary diversification. Am Nat. 2004;164:335–49. 99. Fox J, Friendly M, Monette G. heplots: Visualizing Tests in Multivariate Linear Models. R package version 1.3–5. 2018. 100. Friendly M. HE plots for Multivariate General Linear Models. J Comput Graph Stat. 2007;16:421–44. 101. Friendly M, Fox J. candisc: Visualizing Generalized Canonical Discriminant and Canonical Correlation Analysis. R package version 0.6–5. 2013. 102. Length RV. Least-Squares Means: The R Package lsmeans. J Stat Softw. 2016;69:1–33. 103. Ogle DH, Wheeler P, Dinno A . FSA: Fisheries Stock Analysis. R package version 0.8.30. 2020. 104. Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016. 105. Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5. 106. Jombart T, Ahmed I. adegenet 1.3–1: new tools for the analysis of genome-wide SNP data. Bioinformatics. 2011;27:3070–1. 107. Friendly M, Sigal M. Graphical methods for multivariate linear models in psychological research: An R tutorial. Quant Method Psychol. 2017;13:20–45. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.