Physiological mechanisms of adaptive developmental plasticity in Rana temporaria island populations

Adaptive plasticity is essential for many species to cope with environmental heterogeneity. In particular, developmental plasticity allows organisms with complex life cycles to adaptively adjust the timing of ontogenetic switch points. Size at and time to metamorphosis are reliable fitness indicators in organisms with complex cycles. The physiological machinery of developmental plasticity commonly involves the activation of alternative neuroendocrine pathways, causing metabolic alterations. Nevertheless, we have still incomplete knowledge about how these mechanisms evolve under environments that select for differences in adaptive plasticity. In this study, we investigate the physiological mechanisms underlying divergent degrees of developmental plasticity across Rana temporaria island populations inhabiting different types of pools in northern Sweden. In a laboratory experiment we estimated developmental plasticity of amphibian larvae from six populations coming from three different island habitats: islands with only permanent pools, islands with only ephemeral pools, and islands with a mixture of both types of pools. We exposed larvae of each population to either constant water level or simulated pool drying, and estimated their physiological responses in terms of corticosterone levels, oxidative stress, and telomere length. We found that populations from islands with only temporary pools had a higher degree of developmental plasticity than those from the other two types of habitats. All populations increased their corticosterone levels to a similar extent when subjected to simulated pool drying, and therefore variation in secretion of this hormone does not explain the observed differences among populations. However, tadpoles from islands with temporary pools showed lower constitutive activities of catalase and glutathione reductase, and also showed overall shorter telomeres. The observed differences are indicative of physiological costs of increased developmental plasticity, suggesting that the potential for plasticity is constrained by its costs. Thus, high levels of responsiveness in the developmental rate of tadpoles have evolved in islands with pools at high but variable risk of desiccation. Moreover, the physiological alterations observed may have important consequences for both short-term odds of survival and long term effects on lifespan.


Background
Developmental plasticity occurs when an environmental input induces a lasting alteration in an organism's phenotype. Developmental plasticity is generally a nonreversible process and can modulate acclimation capacity (reversible plasticity) in later developmental-stages [1]. Evolutionary benefits of adaptive developmental plasticity have been broadly assessed by both theoretical studies [2][3][4] and empirical approaches [5,6], reviewed in [7]. Theoretical studies argue that developmental plasticity can increase population viability and facilitate the maintenance of genetic variation in two main ways: reducing the effect of genetic drift by moderating bottlenecks and shielding genetic variants from selection [8,9]. This shielding effect may slow down the response to selection [10], but at the same time increases the odds of population persistence under environmental heterogeneity and preserves genetic variation hence granting a greater adaptive potential [9,10]. Because plasticity confers a higher capacity for surviving new or rapidly changing environments, it can facilitate divergence among populations [11][12][13] and ultimately foster speciation [3,8], as both empirical and theoretical studies have shown.
Environmental heterogeneity affects the degree of adaptive developmental plasticity, as predicted by theory [9,14] and demonstrated empirically [15,16]. Thus, heterogeneous environments are expected to select for highly plastic genotypes, whereas homogeneous environments would tend to reduce plasticity, especially if there are maintenance costs associated to such plasticity [17][18][19]. Populations diverging in the extent of their adaptive plastic responses as a consequence of selection under different environmental regimes would require qualitative or quantitative divergence in the mechanisms underlying phenotypic expression. The expression of alternative phenotypes relies on changes in gene expression [19,20] and/or changes in the regulation of physiological pathways [21]. Thus, understanding the physiological mechanisms regulating developmental plasticity is key to understand the origin of evolved differences among populations in contrasting environments [22][23][24].
Developmental plasticity is particularly critical for amphibians because they are typically species with low vagility and high philopatry to highly variable habitats [25]. The majority of amphibians exhibit ancestral aquatic reproduction [26] and breed in temporary water bodies where the larvae develop until metamorphosis. A major larval-stage risk is pool drying and amphibian larvae can often detect fluctuations in water level and accelerate development in response to decreased water levels to achieve an early metamorphosis [27,28]. Theory predicts that populations exposed to more fluctuating hydroperiods would exhibit greater developmental plasticity since it would allow tadpoles to adjust the timing of metamorphosis accelerating development when ponds dry up or extending the larval period if ponds persist for longer. However, the evolution of developmental plasticity might be constrained by costs of maintaining and/or producing such responses and also limited by factors such as reduced genetic variation or the impossibility of developing appropriate phenotypic responses at all points throughout development [29][30][31]. The populations of Rana temporaria on the Swedish islands are a paradigmatic example of the relationship between environmental heterogeneity and degree of adaptive plasticity [17]. These populations were established between 70 and 800 years ago [32] by frogs migrating from the mainland [33]. The islands have rocky pools where the frogs breed. These pools vary greatly and consistently in depth and size among islands, and consequently vary in average pond duration, ranging from permanent to ephemeral [17]. Over the last decade, several studies have analyzed the adaptive divergence in developmental rate among these R. temporaria populations. These studies have found signs of developmental canalization so that populations occupying islands with only ephemeral pools show overall faster developmental rates than populations from islands with permanent pools [34] while also showing reduced within-population genetic variation for developmental rate [34]. Moreover, island populations show marked differences in plasticity of their developmental rate [35] according to pond duration during the breeding season. In addition, [36] found evidence of costs of developmental plasticity, although these were only noticeable for the most plastic populations [36]. Faster developing populations were also found to express higher levels of thyroid hormone receptors (alpha and beta) and genes associated with higher metabolic activity [20]. This is congruent with mechanisms of developmental acceleration found in other species, namely increased thyroid hormone and corticosterone levels [37].
Developmental acceleration in amphibians is mediated by neuroendocrine pathways, especially by the hypothalamic pituitary adrenal axis (HPA) [38]. The HPA-axis is activated by external environmental inputs such as decreased water height [39], ultimately resulting in increased corticosterone production, which together with thyroid hormone synergistically activate or repress nuclear target genes causing increased metabolic rate, and accelerated morphogenesis [38]. Prolonged corticosterone secretion enhances metabolism and development [40,41] but it can be at the cost of overproduction of toxic substances called reactive oxygen species (ROS) [42] that can inflict considerable cellular damage. Similarly, [43] found that bird species with faster development and smaller body size produced higher antioxidant enzymes. In our study, we expect that R. temporaria populations exposed to divergent regimes of pool drying will show adaptive differences in developmental rate. These among-population differences in developmental plasticity would be explained by altered gene expression associated with adaptive changes in metabolic activity and rate of morphogenesis, which ultimately would involve divergent corticosterone and antioxidant levels.
High ROS production also results in DNA damage, of which telomere shortening is of great importance because of its association with life-history trade-offs and lifespan [44,45]. Telomeres are non-coding tandem repeat sequences of the terminal regions of the chromosomes with high G-C strand asymmetry [46]. Telomeres are determinants of cell senescence and are also involved in chromosome stability by avoiding chromosome fusion [47]. Telomere replication occurs via reverse transcriptase telomerase that adds telomeric repeats to 3′ overhang. However, telomere ends shorten over time after each cell division until critical telomere length is reached and initiates apoptosis, or programmed cell death [48]. Enhanced telomere abrasion has been described as a consequence of compensatory growth and developmental acceleration [49], mainly as a direct consequence of oxidative damage. In the wild, telomere shortening studies are principally focused on ageing [50,51] and body condition quantification [52,53]. However, only a few evolutionary studies include telomeres as a mechanism under selection even though telomere shortening correlates with numerous variations in life-history traits [54], and most of these studies have been conducted on humans [55][56][57].
Here we examine whether adaptive differences in developmental rate among R. temporaria populations are associated with changes in corticosterone levels, oxidative stress, and telomere length. We expected highly plastic populations to increase their corticosterone levels to a greater extent in response to pond drying than less plastic populations. We also hypothesized that populations with high developmental plasticity may pay a cost of maintaining such plastic ability in terms of higher constitutive levels of corticosterone and oxidative stress, and likely shorter telomeres, compared with less plastic populations.

Study area and field sampling
We studied the physiological consequences of adaptive divergence in developmental rate of Rana temporaria tadpoles from six islands located in the Gulf of Bothnia (Umea, Sweden) in a 10 km section of the coastline. The size range of the islands is between 9 and 38 ha. Frogs breed in water filled pools created by rocky depression on these islands. The pools on the islands differ in their water permanence such that some have ephemeral pools, other permanent pools, and others have a mixture of permanent and ephemeral pools. There is no relationship between predator abundance and life history traits among pools on these islands [32], and consequently pond duration seems to be the main environmental factor determining larval development in this system.

Experimental setup
One week after egg sampling tadpoles reached the freefeeding stage and started to swim actively (Gosner stage 25). At that time, 12 tadpoles per clutch (six tadpoles per treatment) for a total of 216 experimental units were individually transferred to 1 L containers filled with 750 mL of reconstituted soft water, where each tadpole was haphazardly assigned to an experimental treatment, i.e. constant water level or simulated pool drying. We renewed water every 4 days and tadpoles were fed ad libitum with lightly boiled spinach on each day of water change [60]. Containers were placed on shelves in a random order with respect to treatment and island in a walk-in climate chamber. Temperature and light were set to 22°C and to 18 h: 6 h of light: dark. The experiment included a water level factor with two levels: constant water and simulating pool drying conditions. In the simulating pool drying treatment we decreased the initial water volume of 750 mL (10 cm) by 33% at each water change starting on day 5 until day 25, keeping the water volume constant at 66 mL (1 cm depth) afterwards. Water temperature did not differ significantly between treatments despite the differences in water level. The experiment ended when tadpoles reached Gosner stage 42 (front legs visible). We monitored tadpoles twice a day (09.00 and 21.00) to check for metamorphs. At this stage tadpoles were photographed to determine their length using ImageJ (version 1.47 t), and were weighed on a high precision balance to the nearest 0.1 mg. Clutch size from each island population varied between 2 and 4 (see Additional file 1).

Tissue collection
Tadpoles at 42 Gosner were euthanized by immersion in a buffered solution of MS-222. A portion of muscle from the tadpoles' tail was removed with a surgical blade and preserved at −20°C for telomere analyses. The rest of the tail was snap frozen in liquid nitrogen and preserved at −80°C until corticosterone assays were conducted. The rest of the body was also snap frozen for oxidative stress assays.

Corticosterone assay
Corticosterone content was determined in tails (50-60 mg) collected from tadpoles at 42 Gosner Stage. Tissue was homogenized with an Ultraturrax TP18/10 (Hanke & Kunkel; IKA-Werk) during 30 s and the hormone was extracted following an organic phase extraction with 1 mL ethyl acetate during 30 min at 4°C and continuous shaking. Samples were then centrifuged at 5000 rpm during 15 min and a known volume of the supernatant was taken and evaporated in a speedVac. Dried elutes were re-suspended in a final volume of 120 μL of the assay buffer provided with the enzymoinmunoassay kit supplemented with ethanol (<5% of final volume) to aid the re-suspension of steroids. We assayed each sample in duplicate (50 μL per sample) for corticosterone determination through specific enzymoinmunoassays (Arbor Assays, K014-H1/H5). The corticosterone antibody has low cross-reactivities to cortisol (0.38%), 11-desoxycorticosterone (12.3%), or progesterone (0.24%). The efficiency of the extraction was checked by spiking several aliquoted samples with 100 pg of exogenous corticosterone prior to extraction and comparing them to the non-spiked aliquotes. Recovery of exogenous corticosterone was never lower than 96%. The lowest point in the corticosterone standard curve was 2.29 pg/mL. To test for assay precision and variability, we determined the coefficient of variation (CV%) for intraand inter-assay variation. Intra-assay variation was 8.74%. Inter-assay variation was 13.23% in the highest point of the standard curve and and 6.86% in the lowest point (n = 8 assays in both cases). The intra-sample correlation coefficient was 0.97. Corticosterone concentration was calculated from the %B/B0 curve by using the 4PLC fitting routine and following the online tool from [61].

Oxidative stress assay
We quantified the activity of three antioxidant enzymes (catalase, glutathione peroxidase, and glutathione reductase). We also quantified malondialdehyde (MDA) concentration, a product formed during lipid peroxidation, total glutathione (GSH t ) and the ratio of oxidized to reduced glutathione (GSH/GSSG ratio). After evisceration, tadpoles were individually homogenized with a Miccra homogenizer (Miccra D-1) at 35,000 rpm in a buffered solution to inhibit proteolysis (1:4; w:v; [62]). The homogenates were centrifuged at 14,000 rpm for 30 min at 4°C and the resulting supernatants were aliquoted into six 0.6 mL tubes and cryopreserved at −80°C. We determined total protein content by standard Bradford's procedure [63]. The coefficient of variation in total protein determinations between duplicated samples was on average 4.27%.
Catalase (CAT) catalytic activity was indirectly quantified following [64]. According to this protocol, potassium permanganate (KMnO 4 ) was used as an oxidizing agent that reacts with the catalase substrate hydrogen peroxide (H 2 O 2 ) giving a red color that can be read at 480 nm 5 min after KMnO 4 was added. We used commercial catalase (SIGMA -60,634) for standard curves preparation. The coefficient of variation between duplicated samples was on average 3.26%. The intra-sample correlation coefficient was 0.95. Glutathione peroxidase (GPX) activity was quantified following the protocol developed by [65]. An excess of glutathione reductase (GR) reduces continually with NADPH the oxidized glutathione (GSSG) producing a constant level of reduced glutathione (GSH). We estimated NADPH oxidation spectrophotometrically at a wavelength of 340 nm. The coefficient of variation between duplicated samples was on average 2.12%. The intra-sample correlation coefficient was 0.88. GR activity was determined assessing the decrease in absorbance at 340 nm due to NADPH oxidation [66], and assays showed an average coefficient of variation between duplicated samples of 2.74%. The intra-sample correlation coefficient was 0.96. Lipid peroxidation was determined according to [67] by measuring MDA concentration. MDA is one product of the lipid peroxidation that reacts with thiobarbituric acids, generating a red product that absorbs at 535 nm. We quantified MDA concentration in nmol / mL by measuring optical density of each sample and then subtracting the blank values and comparing with the calibration curve. The coefficient of variation between duplicated samples was on average 3.93%. The intra-sample correlation coefficient was 0.95. Total glutathione (GSH t ) was determined following [68]. Homogenates were diluted (1:10, w:v) and homogenized in a stock buffer solution (0.01 M PBS and 0.02 M EDTA). Three working solutions were prepared from the same stock buffer as follows: (a) 0.03 mM of NADPH, (b) 6 mM 5,5′-Dithiobis(2-nitrobenzoic acid) (DNTB), and (c) 50 units of GR/mL. Solution a and b were mixed using a ratio of 7:1 (A:B) and then 160 μL of this mixture was added to 40 μL of supernatant. After 15 s, 20 μL of the solution C were added and absorbance was read at 405 nm after 30 and 60 s. Total concentration of glutathione was determined by comparing changes in absorbance between the two readings, according to a standard curve generated by serial dilutions of glutathione from 1 nM to 0.031 nM. The repeatability was 90.3% (n = 10 samples).

Relative telomere length quantification
Genomic DNA for telomere measurements was isolated using a high-salt DNA extraction protocol on a portion of tail muscle. All determinations were made from muscle tissue to avoid confounding differences between tissues in telomere length [69,70]. Relative telomere length assays were performed using quantitative PCR (qPCR) following the protocol developed by [70] and calculated as the ratio of telomere repeats to a single-copy gene (T/S ratio; [71]). As a single copy gene we used glyceraldehyde-3-phosphate dehydrogenase (GAPDH, GenBank accession no. AF255390). Forward and reverse sequences of primers used to amplify GAPDH were 5-AACCAGCCAAGTACGATGACAT-3′ (GAPDH-F) and 5′-CCATCAGCAGCAGCCTTCA-3′ (GAPDH-R), respectively. Forward and reverse primers of the target gene were 5′CGGTTTGTTTGGGTTTGGGTTTGGGTTTG GGTTTGGGTT-3′ (Tel1b) and 5′-GGCTTGCCTTAC CCTTACCCTTACCCTTACCCTTACCCT-3′ (Tel2b), respectively. We performed qPCR in two separate plates for GAPDH and telomere genes by adding 20 ng of genomic DNA from each sample. The set of primers used (Tel1b/ Tel2b and GAPDH-F/GAPDH-R) was at an initial concentration of 900 nM containing 10 μL of Brilliant SYBR Green QPCR Master Mix (Roche) in a final volume of 20 μL. PCR protocol consisted of 10 min at 95°C followed by 30 cycles of 1 min at 56°C and 1 min at 95°C for telomere fragment amplification, and 10 min at 95°C followed by 40 cycles of 1 min at 60°C and 1 min at 95°C for GADPH fragment. We conducted qPCRs on a LightCycler 480 (Roche) and we tested the efficiency of each plate by performing a control standard curve by serially diluting a pool of samples from the different treatments and islands in triplicate (160, 40, 10, 2.5 and 0.66 ng of DNA per well). We calculated the cycle threshold (C t ) value of each sample for each plate. Threshold is the basal level of fluorescence. C t is defined as the number of cycles needed to detect a signal above the threshold. All samples were run in duplicate and relative telomere length was calculated following the formula [72]: where E telomere is the qPCR efficiency of telomere fragment; E GAPDH is the qPCR efficiency of the GAPDH fragment; ΔC t telomere is the C t deviation of controlsample of the telomere (target gene) fragment; ΔC t GAPDH is the C t deviation of controlsample of reference of GADPH (gene of reference) fragment. The coefficient of variation for duplicated samples was 0.65% on average for GAPDH assays and 1.36% for telomere assays. The intra-sample correlation coefficient was 0.90 and 0.86 for the telomere and GAPDH assays, respectively.

Statistical analyses
Statistical analyses were conducted in R, version 3.3.1 (R Development Core Team). We checked whether residuals followed normal distributions conducting Kolgomorov-Smirnov tests ("lillie.test" in "nortest" package, version 1.0-2). We also tested for homoscedasticity using Barlett's tests with the function "bartlett.test" ("car" package; version 2.0-22). We fitted linear and generalized mixed effect models using "lmer" function for normal data and "glmer" for nonnormal data (package "lme4"; version 1.1-7). We used water level (constant or simulated pool drying) and island habitat as fixed factors, and included clutch nested within island as random factors in the models. We used likelihood ratio tests to determine the significance of each predictor. Days to metamorphosis, corticosterone, GPX, GR, GSH, and telomere data were log-transformed to meet parametric assumptions. We estimated growth rate as log (body mass at 42 Gosner stage) log (days to reach 42 Gosner stage). Intraclass correlation coefficients were calculated using "ICCest" function ("ICC" package; version 2.3.0).

Results
Survival was very high throughout the experiment (93.05%) and we found no significant differences between treatments, clutches, or island habitats. Simulated pool-drying induced accelerated development in R. temporaria larvae (χ 2 = 58.05, P < 0.001; Fig. 1a), but the degree of response differed among island habitats as indicated by a significant island habitat by water level interaction (χ 2 = 6.73, P = 0.035; Fig. 1a). Tadpoles from islands with only ephemeral pools showed the highest developmental plasticity since they reduced the time to metamorphosis by 8.13% on average, compared to tadpoles from islands with both pool-drying regimes and with only permanent pools, which accelerated their development by 4.40% and by 4.57% on average, respectively. Pool drying decreased larval growth rate by an average of 12.03% (χ 2 = 42.39, P < 0.001; Fig. 1b). However, we did not find differences in growth rate among island habitats (χ 2 = 0.49, P = 0.784) or in the interaction between island habitat and water level (χ 2 = 1.99, P = 0.390; Fig. 1b). Also, pool drying decreased snout-to-vent length of the emerging metamorphs (χ 2 = 51.01, P < 0.001), but we did not find differences in length among island types (χ 2 = 0.31, P = 0.856) or in the interaction between island type and water level (χ 2 = 2.15, P = 0.341).
Simulated pool desiccation induced endocrine responses in R. temporaria tadpoles, indicated by increased corticosterone levels by an average 25.48% (χ 2 = 8.71, P = 0.003; Fig. S1 in Additional file 2). However, we did not find significant differences between island habitats or their interaction with water level (all P > 0.110; Fig. S1 in Additional file 2).
We also found alterations in antioxidant responses of the different populations in relation to water level. CAT activity decreased by 9.14% on average in tadpoles responding to pool drying (χ 2 = 3.84, P = 0.049; Fig. 2a). Island habitat also affected CAT activity (χ 2 = 6.56, P = 0.038; Fig. 2a) so that it was 11.21% and 16.86% lower in tadpoles from islands with ephemeral pools than in tadpoles from islands with both types of pools or with permanent pools, respectively. CAT levels were not affected by the interaction between pool drying and island habitat (χ 2 = 4.37, P = 0.112; Fig. 2a). Pool drying seemed to induce a slight increase in GPX levels, but the effect was not significant (χ 2 = 2.84 P = 0.092; Fig. S2 in Additional file 2). Neither island habitat nor its interaction with pool drying significantly affected GPX levels (χ 2 = 3.25, P = 0.197, and χ 2 = 0.01, P = 0.996, respectively; Fig. S2 in Additional file 2). Pool drying did not alter GR levels (χ 2 = 2.371, P = 0.124; Fig. 2b). However, we found marginally non-significant differences between island habitats (χ 2 = 5.78, P = 0.056; Fig. 2b). Tadpoles from islands with only ephemeral pools and from islands with both types of pools showed lower constitutive GR levels than tadpoles from permanent pools (6.26% and 7.84% lower GR levels on average, respectively). We did not find differences in the interaction between pool drying and island habitat in GR levels (χ 2 = 1.18, P = 0.555; Fig. 2b).
We found that alterations in antioxidant enzymatic activity were not necessarily associated with cellular oxidative damage, since MDA values were not affected by pool drying, and showed no differences among island habitats, or their interaction (all P > 0.109; Fig. S3 in Additional file 2). Decreased water level resulted on average in 17.31% lower GSH t values (χ 2 = 9.49, P = 0.002). Levels of GSH t was unaltered by either island habitat or the island habitat by water level interaction (χ 2 = 1.90 P = 0.386 and χ 2 = 4.69, P = 0.096, respectively; Fig. S4 in Additional file 2). The ratio of reduced to oxidized glutathione (GSH/GSSG ratio) was not affected by pool drying (χ 2 = 3.16, P = 0.076; Fig. S5 in Additional file 2), and showed no differences among island habitats (χ 2 = 2.54, P = 0.281; Fig. S5 in Additional file 2), or their interaction (χ 2 = 0.36, P = 0.833; Fig. S5 in Additional file 2). Relative telomere length was not affected by water level treatment (χ 2 = 0.01, P = 0.923; Fig. 3), but it varied significantly among island habitats (χ 2 = 6.078, P = 0.048; Fig. 3). Tadpoles from ephemeral pools showed on average 22.31% shorter relative telomere length than tadpoles from permanent pools, and 19.89% shorter relative telomere length than tadpoles from semi-permanent pools. The interaction between relative telomere length and water level was not significant (χ 2 = 1.46, P = 0.482; Fig. 3).
Relative telomere length and GSH t, did not correlate with duration of larval period or growth rate (all P > 0.110). We found, however, a slight negative correlation between GSH/GSSG ratio and larval period (R 2 = 0.089; P = 0.019).

Discussion
Reaction norms can evolve under selection in divergent environmental regimes if plasticity confers adaptive value to individuals developing alternative phenotypes, but also if costs of plasticity select for less plastic genotypes in less variable environments [2,3,17]. Here we found evidence for divergence among Rana temporaria populations in their ability to accelerate development in response to decreased water level in accordance with the predominant pool drying regimes. The studied populations have been estimated to experience a small degree of neutral genetic differentiation [33], although Fst estimates could have been biased downwards due to high heterozygosity of the markers used [73,74] and hence population differentiation may be greater than previously thought. Populations from islands with only ephemeral pools showed the greatest capacity for developmental acceleration in response to pool drying. This indicates that selection on the regulation of developmental rate under different flooding regimes has resulted in increased plasticity in these populations (genetic accommodation), not in canalized fast development or genetic assimilation [2,75]. We suspect that such a steep reaction norm is derived within the species, but phylogenetic reconstructions of the ancestral state of plasticity for these populations would be required to resolve the polarity of the changes in our Rana temporaria populations. Selection for rapid larval development can result in canalized fast development and loss of plasticity [20,76]. In our study, populations inhabiting islands with ephemeral pools have evolved greater developmental plasticity. Our results differ somewhat from a previous study on this system [17], where islands with both types of pools showed higher plasticity. This discrepancy can be a consequence of broader than expected interannual variability in pool duration [17], among-population variation in costs of plasticity maintenance [30], and/or stochastic effects of sampling different genotypes over different years.
The mechanisms underlying developmental acceleration in amphibians are well known, and rely on activation of the hypothalamic-pituitary-axis [27,38]. This neuroendocrine activation in amphibians results in higher thyroid hormone and corticosterone levels, which enhance cell replication rate [54] and morphogenesis [38]. We observed that larvae from all populations increased corticosterone levels to a similar extent (around 25%) when facing decreased water levels. Such upregulation of corticosterone explains the ability of these populations to accelerate development, but does not reflect the observed among-population differences in their degree of plasticity. In contrast, the activity of antioxidant enzymes varied significantly among populations, which would suggest that populations under divergent environments evolved metabolic differences resulting in different levels of ROS production. In this case, highly plastic Rana populations showed lower constitutive levels of both CAT and GR enzymes than less plastic populations. CAT transforms hydrogen peroxide to water and oxygen whereas GR reduces glutathione disulfide to the sulfhydryl form of glutathione, both processes being essential in protecting cells from oxidative damage. Low enzymatic levels are associated with low production of ROS, but also with exhaustion of the enzymes as a result of oxidative stress via toxic substances or ROS production [77][78][79]. In addition, [20] found an increase in the gene transcript of catalase in tadpoles facing pool drying. Levels of catalase transcript were higher in populations that developed faster, supporting the idea of exhausted levels of these enzymes in populations permanently exposed to pool drying conditions. Thus, lower activities of CAT and GR might indicate an enzymatic inactivation caused by an excess of ROS produced during developmental acceleration [45,80]. However, the absence of among-population differences in MDA or GSH levels do not fully support these conclusions since it indicates a lack of intense lipid peroxidation or antioxidant responses. Therefore, an alternative explanation might be that selection favored individuals that maximized ROS production when they experienced higher metabolic rates, hence showing lower antioxidant activity [81,82].
Telomere length varied among populations adapted to different pool-drying regimes by evolving different extents of developmental plasticity. This is a novel and intriguing finding and it may help to understand the implications of developmental plasticity to lifespan and fitness. Telomere shortening occurs after each cell replication so that telomeres usually shorten with age [83,84]. However, telomere shortening is also linked to increased metabolism [54], particularly to oxidative stress originated during intense cell respiration [54,85]. In wild populations, telomere shortening has also been described as a predictor of mortality [86], reproductive costs [87] or the impact of infections [88]. Therefore, telomere shortening can result a reliable indicator of the biological age of individuals [86]. In our experiment, antioxidant responses observed in Rana populations from islands with only ephemeral pools suggest a high ROS production derived from intense metabolism, which is associated with the attrition of telomeres. Thus, based on our oxidative stress and telomere length determinations, we would expect accelerated development to bear the consequence of reduced lifespan and hence possibly reduced fitness. Among-population differences in telomere length have been associated to trans-generational effects of male age at reproduction due to a progressive elongation of telomeres in sperm with age [89][90][91]. Analogous maternal age effects have not been found [90]. Parental exposure to stressful conditions is also relevant to inheritance of telomere length, although these processes remain poorly understood [92]. In our study, shorter telomeres in tadpoles from islands with only ephemeral pools might be related to early age of first reproduction of males, which could be a long-term consequence of accelerated development against pool drying. However, further empirical studies will elucidate underlying mechanisms in telomere inheritance. Evaluations of telomerase activity and long-term studies testing the effects of parental telomere shortening on life-history traits of offspring will help to understand telomeric dynamic across generations.
Physiological differences among populations found in this study suggest that increased developmental plasticity may be associated with metabolic alterations that compromise the health and lifespan of individuals, as populations differed in antioxidant enzymes activity and length of terminal chromosomal regions. Reduced individual lifespan could have cascading demographic effects on population viability, although this remains to be explored.

Conclusions
In sum, populations evolving in contrasting environments showed divergent levels of developmental plasticity and associated oxidative stress and telomere length variation, despite the slight neutral genetic differentiation previously described. These results emphasize the importance of including physiological measurements in the study of phenotypic plasticity, in order to be able to understand the underlying mechanisms of particular evolutionary and ecological processes.

Funding
This study was funded by Ministerio de Economía y Competitividad (grant CGL2012-40044). P.B. is supported by fellowship F.P.U.-AP2010-5373 and by the travel scholarship Est13/00128 from Ministerio de Educación. We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Availability of data materials
All data analysed during this study are included in the Additional file 1.