Phenology and plasticity can prevent adaptive clines in thermal tolerance across temperate mountains: The importance of the elevation‐time axis

Abstract Critical thermal limits (CTmax and CTmin) decrease with elevation, with greater change in CTmin, and the risk to suffer heat and cold stress increasing at the gradient ends. A central prediction is that populations will adapt to the prevailing climatic conditions. Yet, reliable support for such expectation is scant because of the complexity of integrating phenotypic, molecular divergence and organism exposure. We examined intraspecific variation of CTmax and CTmin, neutral variation for 11 microsatellite loci, and micro‐ and macro‐temperatures in larvae from 11 populations of the Galician common frog (Rana parvipalmata) across an elevational gradient, to assess (1) the existence of local adaptation through a PST‐FST comparison, (2) the acclimation scope in both thermal limits, and (3) the vulnerability to suffer acute heat and cold thermal stress, measured at both macro‐ and microclimatic scales. Our study revealed significant microgeographic variation in CTmax and CTmin, and unexpected elevation gradients in pond temperatures. However, variation in CTmax and CTmin could not be attributed to selection because critical thermal limits were not correlated to elevation or temperatures. Differences in breeding phenology among populations resulted in exposure to higher and more variable temperatures at mid and high elevations. Accordingly, mid‐ and high‐elevation populations had higher CTmax and CTmin plasticities than lowland populations, but not more extreme CTmax and CTmin. Thus, our results support the prediction that plasticity and phenological shifts may hinder local adaptation, promoting thermal niche conservatism. This may simply be a consequence of a coupled variation of reproductive timing with elevation (the “elevation‐time axis” for temperature variation). Mid and high mountain populations of R. parvipalmata are more vulnerable to heat and cool impacts than lowland populations during the aquatic phase. All of this contradicts some of the existing predictions on adaptive thermal clines and vulnerability to climate change in elevational gradients.


| INTRODUC TI ON
Climate change is promoting fast increases in both mean temperatures and the frequency of extreme heat events and temporal anomalies, which may jeopardize biodiversity worldwide (IPCC, 2014;Pacifici et al., 2015;Parmesan, 2006). Species basically rely on four strategies to cope with this crisis: evolutionary changes in their tolerance limits, thermal acclimation (phenotypic plasticity), shifts in behavior and activity timing, and shifting geographical ranges in order to track historical climates (Habary et al., 2017;Hoffmann & Sgrò, 2011;Walther et al., 2002). Therefore, the study of population variation and phenotypic plasticity in physiological traits, and the correct characterization of thermal microenvironments can be decisive to predict the consequences of global warming (Camacho et al., 2015;Garland et al., 1991;Huey et al., 2012;Somero, 2010).
Spatial variation in thermal physiology (e.g., Critical Thermal Limits, CT max and CT min ) in relation to latitude and elevation has been thoroughly described at the interspecific level being often associated with environmental thermal heterogeneity (Bozinovic et al., 2011;Pintanel et al., 2019Pintanel et al., , 2022Shah et al., 2017;Stevens, 1989;Sunday et al., 2019). Compared with longer range climatic variation of latitudinal gradients, elevational gradients change climate over shorter distances. That results in predictable changes in air temperatures between −6.5°C/km and − 3.5°C/km, due to adiabatic cooling, and an increase in thermal variability, associated with higher solar radiation and the lowering of air density (Hodkinson, 2005). These steeper climatic gradients may select for thermal adaptations to local conditions and thermal plasticity, and act as physiological barriers to gene flow potentially producing genetic differentiation, particularly in non-seasonal tropical latitudes (Janzen, 1967;Polato et al., 2018).
In turn, under moderate gene flow, populations living at divergent climates could introduce maladapted genotypes into each other, potentially decreasing the frequency of local adapted genotypes.
This would determine a reduction in the steepness of the slope of adaptive clines, as it has been predicted theoretically (Endler, 1977;Slatkin, 1973), and demonstrated empirically in temperate amphibians (Bachmann et al., 2020). Recent intraspecific studies have revealed adaptive clinal variation with elevation in both CT max and CT min (Bishop et al., 2017;Klok & Chown, 2003;Sørensen et al., 2005) with more variation in CT min than CT max (Muñoz et al., 2014). In contrast, a number of studies did not reveal such elevational trend in physiological traits related to thermal tolerances (Buckley et al., 2013;Gvoždík & Castilla, 2001;Senior et al., 2019;Slatyer et al., 2019;Slatyer & Schoville, 2016;Tonione et al., 2020). Yet, reliable data to support adaptive clinal variation is still scant because of the need and difficulty of integrating phenotypic (P ST ) and molecular divergence (F ST ) (Brommer, 2011;Leinonen et al., 2008).
Several biotic and abiotic factors may prevent adaptive differentiation in thermal physiology when analyzing mountain clines. Most animals experience climate at fine-scale patches, and microenvironmental temperatures actually faced by the organism can deviate greatly from recorded mean air temperatures obtained at larger spatial scales (Helmuth, 2009;Potter et al., 2013;Suggitt et al., 2011).
Spatiotemporal changes in microclimate may result, for example, from differences in topography and vegetation cover (Porter et al., 2002). Since microclimatic heterogeneity cannot be captured by downscaling regional climatic variation (Caillon et al., 2014;Diamond & Chick, 2018;Navas et al., 2013;Pincebourde et al., 2016), it becomes important to measure it. This is particularly important for assessing climatic tolerance in animals living at buffered microhabitats, like ponds (Ex. tadpoles, Gutiérrez-Pesquera et al., 2016;Katzenberger et al., 2018;Pintanel et al., 2022). This thermal heterogeneity at the microclimatic scales allows organisms for behavioral thermoregulation, which may preclude the evolution of physiological adaptations in performance by reducing the exposure of organisms to extreme temperatures (i.e., the Bogert effect; Huey et al., 2003Huey et al., , 2012Kearney et al., 2009;Buckley et al., 2015;Farallo et al., 2018;Muñoz, 2022). Otherwise, organisms may be non-active year-round, adopting dormant physiological states such as diapause, hibernation or estivation, in order to escape stressful extreme temperatures (Ragland & Kingsolver, 2008), or because breeding habitats or resources are temporarily unavailable (e.g., ice covered aquatic habitats, and pond drying). Temporal adjustments of activity can be the result of two additive components, one seasonal (or phenological) and one at a finer temporal scale (24-h) that can be dependent on season and particular weather conditions. The limit for these adjustments will be imposed by the energy demand (Kearney & Porter, 2017. Therefore, in absence of energetic constraints, populations can persist at their spatial locations despite change in thermal conditions, without changes in physiological traits. Thus, phenological adjustments in activity can modify the strength of directional selection over thermal tolerance limits through altitudinal gradients (Álvarez et al., 2012;Phillimore et al., 2010;Socolar et al., 2017). This is important because the physiological adjustment of thermal traits may be subjected to Spanish Ministry of Science and Innovation K E Y W O R D S local adaption, microclimate, niche conservatism, phenotypic plasticity, thermal tolerance, warming tolerance

T A X O N O M Y C L A S S I F I C A T I O N
Ecophysiology severe constraints. For instance, extreme heat stress can occur even at high elevations (Sunday et al., 2014), so populations living in mountain areas should face both heat and cold extremes, which leads to an unlikely solution ("master-of-all" superorganism, Remold, 2012).
All these factors can contribute to buffer the realized thermal variation (i.e., the range of effective temperatures experienced by individuals) along elevational gradients and, ultimately, they can promote the pervasiveness and retention of climatic niches; hence, organisms would maintain their thermal niches unchanged while moving along an elevation gradient. In this context, behavioral tracking of the microclimatic niche over space and phenology can allow for retention of the microclimatic niche rather than adapting to the new local conditions with changes in physiology (Farallo et al., 2020;Huey et al., 2003;Kearney et al., 2009;Muñoz, 2022).
Populations living in highly variable thermal environments would express greater plasticity in their thermal tolerances (Angilletta, 2009;Gunderson & Stillman, 2015;Chevin & Hoffmann, 2017;Mallard et al., 2020; but see, for deeper discussion, Gilchrist, 1995, Enriquez-Urzelai et al., 2020. Besides behavioral thermoregulation and phenological adjustments, thermal acclimation may also prevent directional selection on physiological thermal traits, which would constrain thermal adaptation to the new climatic conditions. In fact, organisms can retain plasticity enough to allow for rapid shifts in thermal tolerances, thus precluding thermal stress, death and, likely, the effects of natural selection (Chevin et al., 2010;Huey et al., 2012;Levins, 1969).
Here we examined elevational clinal variation in the upper (CT max ) and lower (CT min ) critical thermal limits and their plasticity in 11 populations of the Galician common frog, Rana parvipalmata ( Figure 1), across an altitudinal gradient from 40 to 1800 m a.s.l. We focus on the aquatic tadpole stage because the breeding aquatic habitat of many amphibian species exhibit low thermal heterogeneity. This limits tadpole ability for behavioral thermoregulation compared with terrestrial adult stages (see Feder & Hofmann, 1999) and, thus, determining that thermal selection on the aquatic phase could be a major driver of tadpole variation in critical thermal limits and its plasticity. In addition, recent research has identified maximum pond temperature as an important range-limiting factor for R. temporaria / Rana parvipalmata (Enriquez-Urzelai, Kearney, et al., 2019). Second, we analyzed population vulnerability to thermal stress by estimating warming and cooling tolerances (sensu Sunday et al., 2014;Gutiérrez-Pesquera et al., 2016). Since average temperatures usually decline with elevation, and thermal physiology limits are driven by peak environmental temperatures (Buckley & Huey, 2016;Gutiérrez-Pesquera et al., 2016;Overgaard et al., 2014;Pintanel et al., 2019), we posit the following two predictions: (1) greater heat impacts are expected at the lowlands Sunday et al., 2014), whereas higher cold impacts are forecasted for mountain populations Sunday et al., 2014) Table S2). We selected a total of 11 populations along an elevational transect between 40 and 1800 m a.s.l.
All of them belong to the T2 (eastern) lineage of Rana parvipalmata (Dufresnes et al., 2020) and were assigned to a maximum of five genetic clusters (Choda, 2014).
For each site, we haphazardly collected 5-7 recently fertil- Experiments) for how to report animal research in scientific publications (https://arriv eguid elines.org/arriv e-guide lines).

| Environmental data
To better characterize potential selective pressures that might lead to thermal local adaptation, we obtained micro-and macroenvironmental thermal data for all the sampled populations. To define the thermal profile at each site, we took into account the population's phenology, encompassing both the breeding and larval periods ( Figure 2; Tables S1-S5). In order to characterize the macroclimatic thermal environments, we used the 'extract' function in the R package raster (Hijmans & van Etten, 2014; R Core Team, 2019) to obtain temperature data from WorldClim 2 layers (the temperature of shaded air at around 2 m off the ground, 30 s or 1 km 2 spatial resolutions; records from 1970 to 2000) (Fick & Hijmans, 2017). For each population, we assessed monthly maximum (TMAX) and minimum (TMIN) temperatures restricted to the time period with presence of larvae in the ponds (Tables S2   and S3). In addition, we calculated seasonal temperature range (SR) as the difference between TMAX and TMIN (Supporting Information, Table S4). Since most animals experience climate at fine-scale patches (Porter et al., 2002;Suggitt et al., 2011), macroclimatic data may be poorer predictors of thermal physiology Hence, we gathered microclimatic data by placing HOBO Pendant temperature dataloggers that recorded temperature every 10-30 min at each pond (Table S5). All the Hobo data-loggers were placed underwater on the bottom of the ponds at a depth of 15-25 cm. One datalogger (Aliva population, 1420 m) was lost and we employed the information from a nearby population with similar elevation and pond characteristics (Pandébano, 1220 m). For each location, we calculated maximum daily temperature (tmax), minimum daily temperature (tmin), seasonal thermal range (sr) (sr = tmaxtmin) and average daily temperature range (dr) for the period when tadpoles are present in the ponds ( Figure 1, Table 1, Table S5).

| Estimation of critical thermal limits and warming and cooling tolerances
Previous research with tadpoles showed that there is not significant variation in CTs across the larval period and only later on, at the verge of metamorphosis climax (Gosner stage 42, Gosner, 1960), both CTs exhibit a strong decline (e.g., Floyd, 1983).
Therefore, to determine the CT max and CT min of R. parvipalmata populations, we haphazardly selected 32 larvae within Gosner stages 26-39 from each population pool (a mix of 5-7 families per population). Then, each population sample was split into two groups of 16 tadpoles for the estimation of CT max and CT min , respectively. Tadpoles were kept individually in 400 ml plastic cups and acclimated inside environmental chambers (FitoClima, Aralab) to a constant temperature of 20°C with a photoperiod of 12 L:12D for 4 days before conducting the tolerance assays. This is the time required for adult amphibians and in tadpoles (J. Turriago and M. Tejedo, unpublished data) to stabilize CT max after a large change in environmental temperature, such as field and laboratory conditions (Hutchison, 1961). Tadpoles were fed ad libitum with Purina rabbit chow. Oxygen saturation in the vessels was daily monitored with a laboratory multi-parameter probe (WTW CellOx® 325) and recorded values were always over 60%. Thermal tolerance limits (CT max and CT min ) were determined using the Hutchison's dynamic method (Gutiérrez-Pesquera et al., 2016;Lutterschmidt & Hutchison, 1997b). Tadpoles were weighed immediately before the beginning of the test and placed in individual 100 ml containers with dechlorinated tap water inside a thermal bath of 15 L (HUBER K15-cc-NR) previously stabilized to 20°C (acclimation temperature and start temperature) for five minutes. Afterwards, each animal was exposed to a constant heating / cooling rate (ΔT = 0.25°C min −1 ; CT max and CT min , respectively) until the endpoint was attained. The endpoint for both thermal limits was defined as the temperature at which tadpoles become motionless and failed to respond to external stimuli (10 consecutive gentle prods with a wooden stick applied in 2 s intervals). Since tadpoles were small in size, we assumed that body temperature was equivalent to water locations appear in Tables S1-S3. All environmental data were obtained only during the breeding period of each population Hutchison, 1997a

| Phenotypic plasticity in CT max and CT min
We studied temperature acclimation in both thermal limits (CT max and CT min ) of five populations, corresponding to low: Nueva   Figure 6). CT max and CT min were determined following the above protocols with a start temperature of 20°C for all acclimation temperatures.

| Molecular markers. DNA extraction
In order to examine the potential for local adaptation in thermal limits (CT max and CT min ), we conducted P ST -F ST comparisons across the 11 populations by using 11 microsatellite loci to assess neutral variation (Table S7). We chose microsatellites because these markers are suitable to define population structure at the fine scale (Camacho-Sánchez et al., 2020;DeFaveri et al., 2013;Lemopoulos et al., 2019;Saint-Pé et al., 2019). To reduce the likelihood of including closely related individuals, we obtained the material for genetic analyses from breeding adults, either as buccal swabs (Pidancier et al., 2003)

| Estimates of neutral genetic and phenotypic divergences (F ST and P ST )
We used the MICRO-CHECKER software (Oosterhout et al., 2004) to check for genotyping errors and null alleles. No evidence of scoring alleles and large alleles dropout was found, but RtU7 was excluded due to the possibility of null alleles. In addition, four markers (Rtempμ1, RtU4, RtμH, and RtμB) were discarded due to failed amplifications. The remainder six markers (Rtempμ2, Rtempμ4, BFG072, BFG093, BFG183, and BFG241) were quantifiable for the 11 experimental populations and therefore used to estimate F ST values. Exact tests for departure from Hardy-Weinberg equilibrium (HWE) and tests for linkage disequilibrium were performed for each population across all loci, and at each locus individually, using GENEPOP v2.1 (Raymond & Rousset, 1995).
Significance was evaluated using the Markov chain methods (Guo & Thompson, 1992) with 5000 dememorizations steps and 1000 batches of 10,000 interactions per batch for HWE, and 5000 interactions for linkage disequilibrium tests.
To test whether local adaptation or neutral divergence drive the elevational variation in CT max and CT min , we compared the genetic differentiation (F ST ) and the genetic divergence of quantitative traits (Q ST ) among populations, which allow to discern whether trait differentiation is due to genetic drift or natural selection (Leinonen et al., 2008). F ST estimates were calculated according to Weir and Cockerham (1984), using FSTAT 2.9.3.2 (Goudet, 2002). As a surrogate for Q ST data, we used a measure of phenotypic divergence of a trait (P ST ) (Brommer, 2011), which can be calculated by the equation: where 2 B is the phenotypic variance between populations, 2 w , the phenotypic variance within-population, and h 2 is the character heritability.
The constant c represents the proportion of the total variance due to additive genetic effects across populations (Leinonen et al., 2006).
To obtain the P ST values for each pair of populations, we used a linear mixed-effects model (LMM), with population defined as a random factor, using the LMe4 r-package: CT max ~ 1 + (1|Population), and CT min ~ 1 + (1|Population) (Bates et al., 2015). We used the error variance as a proxy for 2 w (within-population variance) and the intercept variation for 2 B (variance between populations). Thus, P ST estimates depend on the ratio c h 2 . Since these parameters are extremely challenging to obtain in the wild and usually unknown (Pujol et al., 2008), we considered a set of values to calculate P ST (Brommer, 2011). We constructed several matrices for the P ST values obtained for different values of c and h 2 . For each possible combination, the overall mean values (overall P ST ) and their 95% confidence intervals were calculated using a nonparametric bootstrap and compared with the upper limit of the confidence interval for overall F ST (Supporting Information S7-S9). The value of the c/h 2 ratio at which the lower confidence interval for P ST and the upper F ST estimates overlap, can be interpreted as a measure of the robustness of the difference between F ST and P ST estimates (see Brommer, 2011). The use of P ST presents several caveats, since non-additive genetic variance (epistasis or dominance effects), maternal effects or environmental factors and genotype-environment interaction, can lead to a distorted picture of additive genetic variation when studying only phenotypic variation in natural conditions (Brommer, 2011;Leinonen et al., 2008;Pujol et al., 2008). However, because experimental individuals were raised and analyzed under the same environmental conditions, we can assume a lower risk of unwanted effects.

| Statistical analyses
To test for geographic variation in CT max and CT min among populations of R. parvipalmata, we fitted two separate ANCOVA models for CT max and CT min , as dependent variables, including populations as categorical factor, and body mass as a covariate. Tukey HSD posthoc tests were conducted to identify which populations differed in their thermal limits.
To test for covariation between CT max and CT min , thermal variation and elevation, we assessed linear and quadratic regression models using elevation and the thermal data as independent variables and CT max , CT min , WT, wt, CT, ct as dependent variables. This allowed to determine the relationship between species' physiological limits and environmental thermal predictors, including elevation.

Paired t-Test and non-parametric test of Kolmogorov-Smirnov (KS)
were used to assess differences between estimates of vulnerability of exposure to extreme temperatures, determined using macroclimate (WT, CT) or microclimate (wt, ct) thermal data. Finally, we ran Mantel tests with 999 permutations in the R-package 'vegan' (Oksanen et al., 2018) to test the correlations between P ST and F ST pairwise population matrices.
In order to evaluate whether acclimation effects differ among populations, we used a two-way analysis of variance of CT max and CT min , with temperature (6, 13, 20, and 27°C) and population as fixed factors. We estimated the Acclimation Response Ratio (ARR) that measures the change in thermal tolerance relative to change in acclimation temperature (Claussen, 1977;Ruthsatz et al., 2022). In the case of ARR for upper thermal tolerance, ARR = [(highest CT max -lowest CT max ) /Δ°C]. Then it provides a metric of thermal plasticity capturing acclimation capacity, thus allowing standardized comparisons between populations and critical thermal limits, (e.g., whether beneficial acclimation is greater for CT min than CT max ). During CT min estimates at the lower acclimation temperatures (6 and 13°C), water often reached crystallization (exothermic freezing reaction) before the immobility endpoint was attained. Most of the tadpoles briefly exposed to the freezing point recovered activity within a few seconds and survived the 24-h period, indicating resistance to these extreme cold temperatures, but the actual CT min was inestimable (i.e., in these cases CT min was lower [cooler] than the observed freezing point of water). Thus, in order to avoid bias associated with either biased sample or a biased CT min estimates, we restricted the statistical analyses of CT min acclimation scope to the 20-27°C range. All the statistical analyses were performed in R 3.6.1 (R Core Team, 2019).
A preliminary analysis revealed that tadpoles from populations between 900 and 1700 m were smaller than tadpoles from low elevation (<700 m) and the highest elevation (1800 m) (ANOVA;  (Tables 1 and Table S1) although they were w not affected by either elevation or any or macro-and microclimatic predictor ( Figure 4; Table 2). Regarding vulnerability to heat stress, we also found contrasting differences between warming tolerance estimates based on macro-(WT; WorldClim) and micro-climate data (wt) (t = 2.74, df = 10, p = .02; D = 0.82, p < .01). These differences were remarkable in two ways. First, WT was consistently higher than wt, except for two populations at the lowest elevation ( Figure 5a).
The overall value of neutral differentiation (F ST ) of R. parvipalmata was 0.066 (95% CI: 0.058-0.075) and revealed significant genetic differentiation among the 11 populations (Table S8) We found significant divergence among populations in the level of plasticity of CT max and CT min to variation in acclimation temperature (population × acclimation interaction; Table 3). Acclimation to warm temperatures resulted in higher CT max and CT min (Table 4, Figure 6). Mid and high-altitude populations showed the highest phenotypic plasticity for CT min (ARR; mean ± 1SD, 0.32 ± 0.05, n = 3, Table 4), being twice greater in magnitude when compared with the ARRs of lowland populations (mean ± 1SD, 0.16 ± 0.01, n = 2; F I G U R E 3 Elevational variation in: Absolute maximum (a) and minimum temperatures (b) using the WorldClim database (TMAX and TMIN), and microenvironmental pond datalogger information (tmax and tmin). Absolute seasonal temperature range (c) (st = tmax-tmin) and mean daily temperature range (dr) (daily tmax-tmin) (d), based in microenvironmental pond datalogger information. Thermal information corresponds only to the larval period for each studied location (see Figure 1b). TA B L E 2 Results for the linear and quadratic regressions between temperature data, critical thermal limits elevation, warming tolerance (WT, wt) and cooling tolerance (CT, ct).

F I G U R E 4
Boxplot showing the variation of thermal tolerance limits along the elevational gradient. The first and third quartiles ("hinges") and the 95% confidence interval of the median ("notches") are shown. The dashed line indicates the mean CT values of R. parvipalmata for the overall data, and dots placed past the line edges indicate outliers. Table 4). Greater ARR indexes for CT max were obtained for the midelevation populations (Table 4).

| DISCUSS ION
Thermal limits (CT max , CT min ) of tadpoles varied between Rana parvipalmata populations across the elevation gradient. However, CT max and CT min did not correlate with elevation nor with any of the macro-and microclimate predictors, suggesting niche conservatism. We did find some indications of directional selection as the divergence in thermal physiology limits (P ST ) tended to be slightly greater than the neutral The prevalence of intraspecific adaptive clinal variation with elevation in thermal traits is still a contentious topic in current F I G U R E 5 Estimates of Warming Tolerance: CT max minus TMAX or tmax, for either WorldClim or microclimate pond maximum temperature estimates (datalogger), respectively (a), and Cooling Tolerances: TMIN or tmin, for either WorldClim or microclimate pond minimum temperatures (datalogger), respectively, minus CT min (b), for 11 populations of R. parvipalmata, considering only the larval period. Triangles: estimates based on air temperature from WorldClim. Squares: estimates based on water temperatures registered in ponds with dataloggers.

Elevation (m)
Warming Tolerance ( (Berven et al., 1979;Luquet et al., 2015;Richter-Boix et al., 2015), including the closely related Rana temporaria (Laugen et al., 2003;Lind et al., 2011;Muir et al., 2014). However, our findings were consistent with the absence of elevational variation in thermal sensitivity of locomotion and thermotolerance reported for post-metamorphic and adults of Rana parvipalmata in the same study system (Enriquez-Urzelai et al., 2018. In that case, the terrestrial stages of amphibians can thermoregulate via micro-habitat selection or activity timing, sheltering underground in crevices and rodent burrows to avoid peak temperatures (Enriquez-Urzelai et al., 2020). This is parallel to the known Bogert effect in many ecotherms, (Bogert, 1949;Buckley et al., 2015;Farallo et al., 2018;Huey et al., 2003;Muñoz, 2022;Muñoz & Losos, 2018). However, the scope for behavioral thermoregulation in water is much more limited than in the terrestrial environment due to its high specific heat capacity and conductivity that reduces spatial thermal heterogeneity. Amphibian tadpoles, although able to thermoregulate (Hutchison & Dupré, 1992), can be exposed to unavoidably thermal stress, particularly in sunlit ponds without canopy cover (Duarte et al., 2012). Exposure to sunlight is prevalent in the breeding habitats of mid-and high-elevation populations of R. parvipalmata, which are exposed to relatively stressful high temperatures with warming tolerances <8°C ( Figure 5a) and a wide range of both seasonal and short-term thermal variation ( Figure 3d).
Shifts in the reproductive period and microclimatic conditions that are conditioned by elevation may dampen temperature changes along the gradient and prevent the expected linear decrease in both critical thermal limits. Therefore, thermal conservatism in tolerance limits together the absence of local adaptation may simply be a consequence of a coupled variation of reproductive timing with elevation (the "elevation-time axis" for temperature variation, in contrast with the elevation axis). In montane areas, altitudinal variation in the timing of reproduction seems to be constrained by hydroperiod rather than temperature itself, with higher altitude populations delaying breeding until snow melting (Álvarez et al., 2012;Corn, 2003).
In addition, the time of spawning may be genetically determined in mountain populations (Álvarez et al., 2012;Phillimore et al., 2010;Wilczek et al., 2010). In this sense, it appears that during warm winters, when early snow melting occurs, frogs still delay reproduction until a threshold time is reached (Álvarez et al., 2012).
Phenotypic plasticity of thermal limits matched the observed variability in temperature through the elevation gradient. Populations from mid and high-elevations showed higher levels of plasticity than low-elevation populations, especially for CT min . These populations, particularly those from mid-elevation, are also exposed to greater seasonal and daily thermal ranges, which supports the idea that phenotypic plasticity in critical thermal limits can be a response to the increased environment thermal variability. Similar adaptive plasticity in thermal tolerances are shown by populations of two toad species exposed to more variable climates (Alveal-Riquelme et al., 2016;McCann et al., 2018). Furthermore, since mid-elevation populations showed the lowest warming and cooling tolerances, and both thermal limits showed no clear clines, our data suggest the existence of a trade-off between phenotypic plasticity and tolerance to environmental temperatures (Stillman, 2003 Iberia, likely expanded during the cold glacial cycles (Dufresnes et al., 2020). Thus, it is possible that these extremely low CT min are remnants of adaption to environmental conditions from that period, resulting in no current micro-geographic adaptive differentiation along elevation gradients. This hypothesis of 'evolutionary anachronism' (Janzen & Martin, 1982; see also Qu & Wiens, 2020, Moreira et al., 2021 was supported by our finding that lowland populations showed extreme cold tolerances (e.g., Color, 380 m) although they presented the higher minimum temperature recorded.
The contrasting estimates of warming tolerances derived from macro-and microclimate data highlight the importance of monitoring the microhabitat when assessing vulnerability to global warming (Baudier et al., 2015;Gutiérrez-Pesquera et al., 2016;Katzenberger et al., 2018;Enriquez-Urzelai, Kearney, et al., 2019;Pintanel et al., 2019Pintanel et al., , 2022; see also Sunday et al., 2014 and higher tmax were found in mid-elevation populations and not at the high-peak populations. This a priori unexpected outcome is likely the result of phenological shifts in these local populations and differences in habitat structure (i.e., canopy cover, topographical shadow) between low-, mid-, and high-elevation wetlands. In the study area, the breeding habitat of lowland (below 500-700 m) common frogs consists of very small and shallow waters scattered on a rather humanized landscape (e.g., track pools, ditches, and less often small temporary ponds), and located in valley bottoms with dense canopy cover, which prevents direct beam solar radiation. In contrast, breeding habitats in mid-altitude areas (700-1200 m) are small, shallow, temporary ponds, most often located on high plains and hills without canopy cover, and therefore, exposed to high levels of direct solar radiation and a low thermal buffering. Finally, although high altitude wetlands (1300-2100 m) can be affected by topographic shading, most often they lack canopy cover, which besides a thinner atmosphere leads to low thermal buffering. This, along with the change in reproductive phenology (the elevationtime axis), may reduce the actual differences in temperatures experienced by larvae at different elevations.
Phenological shifts and microgeographic variation in habitat structure can determine the thermal regimes experienced by populations along mountain gradients (see Muñoz & Losos, 2018 has reduced the thermal differences between populations, thus hindering physiological evolution (see also Enríquez-Urzelai et al., 2018;Muñoz, 2022).

| CON CLUS IONS
Directional changes in reproductive phenology and phenotypic plasticity can block local thermal adaptation by lowering selective pressures for population differentiation, (consistent with niche conservatism hypothesis; Muñoz & Losos, 2018). Therefore, although "phenological buffering" may override thermal stress under current conditions, it could hinder long-term adaptation to climate change, potentially compromising long-term population sensitivity (Buckley et al., 2015;Enriquez-Urzelai et al., 2018;Kearney et al., 2009). This agrees with the idea that phenotypic responses do not occur in selective vacuums, and therefore, any adjustment in one response can cause an evolutionary ripple in others (Muñoz, 2022).