Thermal–metabolic phenotypes of the lizard Podarcis muralis differ across elevation, but converge in high-elevation hypoxia

To identify how organisms are reacting to rapidly changing environmental conditions, we must integrate patterns of interaction among physiology, behavior and abiotic factors. In addition to evidence for in situ adaptation (Bradshaw and Holzapfel, 2008; Catullo et al., 2019) and plasticity (Charmantier et al., 2008; Gunderson and Stillman, 2015; Kelly, 2019), habitat tracking is a common response (Pauchard et al., 2016). As environments warm, species can move toward higher elevations or poleward to track the emergence of their preferred thermal conditions (Dirnböck et al., 2011). Movement upslope, however, also presents multiple novel environmental conditions, including varying temperature regimes, changing plant communities and, of interest for this study, reduced oxygen availability (Storz et al., 2010; Jacobsen, 2020). By examining how behavior and physiology interact in response to changes in both temperature and oxygen availability, we can better anticipate phenotypic shifts and broader biogeographic patterns as global change marches on (Huey et al., 2012; Clusella-Trullas and Chown, 2014).

A multi-dimensional phenotypic approach is imperative because traits do not evolve or acclimate to environmental variation in isolation (Forsman, 2015; Muñoz and Losos, 2018; Bodensteiner et al., 2021a). The thermal dependence of physiological and performance traits is most often characterized with thermal performance curves (TPCs), which describe the relationship between core body temperature and a performance dimension (e.g. sprint speed). From these curves, we can extract key parameters, such as the thermal limits of performance, the optimal temperature for performance and total thermal breadth of performance (Huey and Stevenson, 1979; Taylor et al., 2021). The components of the TPC do not scale linearly across levels of organization (from molecules to organisms), nor are the effects of variation among physiological traits purely additive (Angilletta, 2009; Schulte et al., 2011). Instead, different parameters of the TPC and other thermally dependent traits can shift following a variety of patterns (reviewed in Bodensteiner et al., 2021a). Some traits, such as thermal preference, optimal sprinting temperature and heat tolerance, are correlated within individuals in some species, but not others (Van Damme et al., 1989; Huey et al., 2012). Other traits may evolve independently, such as upper and lower critical limits (Kellermann et al., 2012; Muñoz et al., 2014, but see Bodensteiner et al., 2021a). Importantly, thermal traits interact with other physiological pathways (e.g. metabolism and aerobic performance) to maintain whole-organism function (Gangloff and Telemeco, 2018; Ern, 2019). The thermodynamic effect on molecular processes increases metabolic rate, which can lead to a mismatch between oxygen supply and demand at temperature extremes (Pörtner, 2002; Pörtner et al., 2017). Correspondingly, interactions among these mechanisms can shape thermal performance and tolerance in a context-dependent fashion (Gangloff and Telemeco, 2018). In response to hypoxia, for example, hematological parameters affecting blood-oxygen carrying capacity shift in lizards, impacting both aerobic metabolism and thermal performance (Lu et al., 2015; González-Morales et al., 2017; Gangloff et al., 2019).

We present here an experiment designed to characterize the complex relationships among thermal–metabolic traits within individuals, describe how these traits differ among individuals from populations with different environmental characteristics, and quantify whether and how traits shift in response to novel high-elevation hypoxia exposure. We focused on common wall lizards (Podarcis muralis) from replicate populations in different elevational zones. This is a primarily lowland species, with historical populations extending to about 2000 m above sea level (ASL). Recent observations demonstrate upslope migration, up to 2300 m ASL in the Pyrénées (Pottier, 2017; F.A., personal observation), with further expansion expected as the environment continues to warm. Upslope colonization, however, will also expose lizards to novel oxygen conditions not experienced in their ancestral range. Given their large population sizes with a high degree of genetic variation (Michaelides et al., 2015; Beninde et al., 2018) and the range expansions and contractions they have experienced as a result of shifting climates over their evolutionary history (Gassert et al., 2013; Yang et al., 2021), we expect P. muralis to possess the raw genetic material needed for adaptive responses, thus offering a potential model for characterizing how phenotypic pathways respond to changing environments (Kelly, 2019). Our study combines field observations and laboratory experiments to address three primary questions. (1) Do low- and high-elevation lizards differ in their thermal–metabolic phenotypes when measured at their elevation of origin? (2) How do lizards shift thermal–metabolic phenotypes in response to novel hypoxia at extreme high elevation? (3) Do patterns of acclimation in response to captive and novel oxygen conditions differ between low- and high-elevation lizards?

To address these questions, we quantified nine physiological traits: thermal preference (T_pref), selected temperature range (T_sel), critical thermal minimum (CT_min), critical thermal maximum (CT_max), hematocrit (packed red blood cell density), hemoglobin concentration, running endurance, resting metabolic rate (RMR) and maximal post-exhaustion metabolic rate (MMR). We predicted that some of these traits will differ between lizards inhabiting different elevations when measured at their elevation of origin and some traits will shift in response to novel oxygen conditions (see detailed trait-specific predictions in Table 1). We additionally described phenotypes in multivariate space and tested for differences between groups, predicting that lizards from high elevation, which are pre-exposed to moderate levels of hypoxia, will shift less under exposure to more extreme oxygen reduction. We also estimated thermoregulatory patterns for lizards at both elevations, and predicted that lizards at high elevation will be more effective thermoregulators to compensate for their colder habitats (Caldwell et al., 2017; Ortega and Martín-Vallejo, 2019). Our study quantifies spatial variation in trait expression and plastic potential of traits in the ecological context of range expansion. Given that traits interact within individuals and that different trait combinations can lead to the same functional outcome, such a multivariate approach is essential in leveraging trait-based approaches to predict patterns at broader scales (Beissinger and Riddell, 2021).

Field and lab protocols were conducted under permit from the Direction régionale de l'environnement, de l'aménagement et du logement Midi-Pyrénées (Arrêté Préfectoral No: 2017-s-02 du 30 mars 2017), under current ethical committee approval (APAFIS#16359), and in accordance with Directive 2010/63/EU. The lizard carcass used for copper model validation was from OH, USA, under Ohio Division of Wildlife Wild Animal Permit (23-014) and Ohio Wesleyan University IACUC (protocol Spring_2020-21-04_Gangloff_B).

This experiment used adult male common wall lizards [Lacertidae: Podarcis muralis (Laurenti 1768)] from within their ancestral elevational range in the Departments of Ariège and Hautes-Pyrénées in southern France. We used a small lasso attached to an extendible pole to collect animals from low-elevation populations (425–473 m ASL; N=25) on 3–4 September 2018 and from high-elevation populations (1557–1812 m ASL; N=24) on 10–11 September 2018 (see Table S1 for sample size details). On the day of capture, we weighed lizards using a digital scale (AMIR, ±0.01 g; range 3.18–8.50 g, mean 5.75 g) and measured snout–vent length (SVL) using digital calipers (CD-6, Mitutoyo; SVL ±0.01 mm; range 50.90–68.42 mm, mean 60.21 mm). Lizards from low- and high-elevation populations did not differ in mass (t_95.4=−0.83, P=0.410) but lizards from low-elevation populations were ∼4% larger in SVL (low-elevation mean 61.42 mm; high-elevation mean 58.96 mm; t_86.6=−3.23, P=0.002). The day of capture, we transported animals to a laboratory at approximately the same elevation, either low elevation (436 m ASL) or high elevation (1735 m ASL; Table S1). Lizards were maintained under common conditions in each lab in groups of 3–6 individuals in plastic enclosures, with water provided ad libitum, and a basking platform that also served as a shelter. Ambient temperature varied between 15 and 20°C and light was provided 14 h day⁻¹ with fluorescent bulbs. During daytime hours, incandescent heat bulbs supplied a temperature gradient of ca. 25–40°C for 6 h day⁻¹ at 1 h intervals.

After 2–3 days of acclimation, we began experiments at the elevation of origin. Lizards were not fed until after the first set of measurements to ensure a post-absorptive state (Van Damme et al., 1991; Pafilis et al., 2007). We conducted measures on two consecutive days, fed and fasted lizards again, and then completed measures over several days to ensure that animals would be nourished, but post-absorptive, during measurements. All trials were conducted during daylight hours (09:30 h–20:00 h). Within each elevation cohort, we divided lizards into three groups and conducted a single assay of thermal physiology or performance on a given day: (1) thermal preference trial, (2) CT_max, (3) CT_min, (4) RMR, running endurance and MMR (as in Gangloff et al., 2019). Because of time constraints, only a subgroup of lizards could be measured for metabolic rates and endurance (Tables S1 and S2).

After initial measures, lizards were divided randomly into two treatment groups within each elevation cohort: one treatment remained at the elevation of origin to serve as a control and one treatment was transported to the Observatoire Pic du Midi de Bigorre at extreme high elevation, currently above the range limit of P. muralis (2877 m ASL; Table S1; Pottier, 2017). Lizards in the two treatment groups did not differ in size within either elevation cohort (low elevation: t_47.5=−1.33, P=0.19; high elevation: t_44.9=0.87, P=0.39). Lizards were maintained under standard husbandry conditions at low, high or extreme high elevation for 15–27 days and then all traits were re-measured. Previous work demonstrates that, within 3 weeks of transplantation to extreme high elevation, adult male lizards exhibit changes in blood biochemistry, sprint performance, running endurance, MMR and body condition (Gangloff et al., 2019). Our experimental design thus allows a test for differences between lizards from low- and high-elevation populations measured at their elevation of origin, changes in thermal biology over time in captivity, and plastic responses in lizards originating from different elevations to extreme high elevation. Following experiments, all animals were returned to their site of capture.

To measure core body temperature during the thermal tolerance experiments, a temperature probe (Type T, 36 gauge, Omega Engineering, ±0.1°C) connected to a digital temperature logger (HH806AU, Omega) was placed approximately 1 cm into the cloaca of each lizard and secured with medical tape. Lizards were placed in a plastic container in which they could move freely and acclimate to room temperature (∼20°C). The plastic container was then placed into a polystyrene foam cooler layered with ice (for CT_min trials) or heated with a 100 W light bulb (for CT_max trials). The lizard was then systematically cooled or warmed by ∼1°C min⁻¹. Once the lizard reached 15°C (for CT_min trials) or reached the panting threshold (Hertz et al., 1979) or 30°C (for CT_max trials), it was flipped onto its back and stimulated to right itself by prodding at the base of its tail and thighs with blunt tweezers. The flipping procedure was repeated every 0.5°C change in temperature until the lizard was unable to right itself after 15 s, with that temperature being defined as CT_min and CT_max.

We estimated the preferred body temperature (T_pref) and the preferred temperature range or the set-point range (T_sel) by placing lizards within a laboratory thermal gradient during their active hours (Taylor et al., 2021). T_sel refers to the range of the central 50% of body temperatures (interquartile range) measured from lizards that have been placed in a thermal gradient and allowed to choose where to sit, and T_pref is the mean of this range (Huey, 1982; Hertz et al., 1993; Taylor et al., 2021). The thermal preference arenas were constructed from smooth corrugated plastic (91 cm×60 cm×39 cm) with four identical lanes in each arena. We created a stable thermal gradient (20–40°C) with two ceramic heating bulbs. Lizards acclimated in the arena for 30 min before we recorded core body temperature every 5 min for 3 h with a thermocouple (40 gauge) that was inserted into the cloaca following the above protocol (Muñoz and Losos, 2018).

We measured RMR, running endurance and MMR consecutively in a subset of animals from one population at each elevation (Tables S1 and S2). Lizards were first placed in a 250 ml opaque plastic metabolic chamber and acclimated to 32°C (approximating the mean of field body temperatures recorded in this study and within the preferred temperature range of this species in this region; Trochet et al., 2018) in an incubator (ExoTerra Model PT-2445, Hagen Inc.) for 30 min. We then utilized pull-mode respirometry to continuously measure gas exchange for 30 min using a FoxBox-C Field O₂ and CO₂ Analysis System (Sable Systems, Inc.). Air was circulated at a rate of 500 ml min⁻¹ through the metabolic chamber, dried of water vapor with Drierite (W. A. Hammond Drierite Co., Ltd), and then measured for both O₂ and CO₂ content, corrected for barometric pressure. We calculated oxygen consumption (V̇_O2) and carbon dioxide production (V̇_CO2) rates following the equations of Lighton (2008) and extracted the lowest values averaged over a 10 min window for our value of RMR (as in Kouyoumdjian et al., 2019) using ExpeData software (v.1.7.30, Sable Systems, Inc.). After measuring RMR, we ran lizards to exhaustion in a sand-lined arena (82×31 cm) that was kept warm (∼35–38°C) with submerged heating cables to maintain lizard body temperature near the preferred value. One of us (E.J.G.) continuously chased the lizard with a soft paintbrush until further stimulation elicited no locomotor response for 5 s (Gleeson and Dalessio, 1989; Gangloff et al., 2019). Time (s) from the start of the trial to exhaustion, recorded with a stopwatch, was the measure of running endurance. Immediately after (within 10 s), we replaced lizards in a metabolic chamber and measured gas exchange for an additional 10 min. We extracted the highest rates of gas exchange over a 15 s interval as our measure of MMR (V̇_CO2,peak), excluding readings for the first 3 min to account for the time required to flush a chamber of this size at this flow rate. Measurement of gas exchange would ideally commence instantaneously; our protocol nonetheless captures the peak or near-peak in post-exercise aerobic metabolic rate (Hailey et al., 1987; Gangloff et al., 2019), which represents an important aspect of the energetic costs of locomotor activity (Gleeson and Hancock, 2002). Analysis of V̇_CO2 and V̇_O2 did not qualitatively differ; because our readings for CO₂ were consistently more reliable than for O₂ (due to sensor drift), we present data on V̇_CO2 only.

At the time of capture and then again during the second measurement period, we collected blood samples (25–30 μl) from the retro-orbital sinus using a heparinized glass capillary tube (MacLean et al., 1973), which was then stored on ice until processing. We quantified hematocrit and hemoglobin concentration following Gangloff et al. (2019). We measured hemoglobin in samples run in duplicate on three plates, with each plate including a pooled sample to provide estimates of intra- and inter-plate variation (CV of 5.4% and 9.1%, respectively).

The operative temperature (T_e) is the stable temperature of an organism when they are not behaviorally thermoregulating (Bakken, 1992). We deployed 20 electroform copper lizard models containing an iButton temperature logger (Model DS1921G-F5, Maxim Integrated) at one low and one high elevation site following Muñoz et al. (2014). Copper models were placed haphazardly in the environment (excluding roads and waterways) the night before sampling and perches were targeted to include microhabitats observed to be occupied by lizards from this and previous sampling efforts. Models recorded temperature every 10 min (07:00 h–19:00 h) on sampling days. Copper model accuracy was assessed by measuring the equilibrium temperature of models and a lizard carcass across a relevant range of temperatures (15–35°C in 5°C increments). Thermocouples were attached inside the model and inside the lizard's cloaca. We then allowed both to equilibrate at each temperature (∼1.5 h) in a programmed incubator. We then used the equation from a regression of lizard temperature on the model temperature (R²=0.998; y=0.779x+5.008) to correct recorded field temperatures from the copper models. On 4 September 2018 (Aubert, low elevation) and 11 September 2018 (La Mongie, high elevation) we sampled field body temperature (T_b) during daylight hours. Lizards were captured via methods detailed above and a T_b measure was taken by rapidly inserting a thermocouple connected to a digital thermometer (Omega HH801) ∼1 cm into the lizard's cloaca. We identified sex via external morphology and marked lizards before release with non-toxic paint to prevent resampling. We used a t-test to compare T_b between low- and high-elevation lizards.

We utilized linear mixed models to assess the relative influence of origin elevation, response to captive conditions and response to translocation on nine dependent variables: T_pref, T_sel, CT_min, CT_max, hematocrit, hemoglobin concentration, running endurance, RMR (V̇_CO2) and MMR (V̇_CO2,peak). Models included a single categorical factor composed of six levels, with levels describing the combination of elevation of origin (low versus high elevation), test location (elevation of origin versus extreme high elevation) and time point (first or second measure). We then tested our specific a priori hypotheses using linear contrasts of estimated marginal means using the emmeans package (Lenth et al., 2018), correcting for multiple comparisons with the Šidák method (see Table S2). Because our experiment was designed to first test all lizards at their elevation of origin and then the response to translocation, we did not have a full-factorial design. This modeling approach, therefore, decreases the number of parameters estimated and allows tests of our hypotheses. We included a random intercept for each individual and included a covariate of body size (mass for metabolic rate measures, SVL for all other dependent variables). If this covariate was not an important predictor (P>0.10), we removed this from final models. To meet the assumption of normally distributed residuals, we log₁₀-transformed T_sel, endurance, RMR and MMR before analysis. Models were implemented with the lme4 package (Bates et al., 2015) in R (http://www.R-project.org/). We visually assessed distributions of model residuals and determined the relative importance of fixed effects using type III sums of squares with corrected denominator degrees of freedom (Kenward and Roger, 1997). All data figures were created with ggplot2 (Wickham, 2011).

We analyzed the multivariate phenotype in a unified analysis that allowed us to simultaneously quantify the within-individual correlations of traits, differences in traits among lizards from different elevations, and shifts in phenotype in response to the hypoxia of extreme high elevation. We included in this analysis the same nine traits listed above. Individuals from a single low- and high-elevation population were measured at both their elevation of origin and extreme high elevation. As in the univariate analyses, T_sel and endurance were log₁₀-transformed before analysis. For measures of RMR and MMR, we calculated the rate of oxygen consumption per gram and then log₁₀-transformed these values before analysis. We utilized a non-parametric multivariate analysis of variance (NP-MANOVA) with residual randomization in permutation procedure (RRPP; Collyer and Adams, 2018), implemented in R (http://www.R-project.org/) and following Telemeco and Gangloff (2020). Because this approach requires a complete dataset (all measures for all individuals at a given measurement time), we included only those individuals for which we measured metabolic rates (see above; 46 measures of 9 variables on 23 individuals; Table S1). We were missing 6.3% (26 of 414) of individual trait measures. Following the guidelines of Telemeco and Gangloff (2020), we imputed these values using multiple regression models including all other measured traits. We then created a model that included a categorical factor combining elevation of origin, test location and time point to compare the multidimensional phenotype among different groups, with significance determined from 999 iterations of the residual randomization procedure. We extracted least-squares means in multidimensional space to compare phenotypes among groups, conducted pairwise tests for differences in all group combinations, and extracted least-squares means and 95% confidence intervals for each of the traits included in the multivariate response matrix.

T_pref, T_sel and CT_min did not differ between treatment groups and we found no evidence of shifts in these traits in response to novel oxygen conditions (Table S2, Fig. S1). When combining measurements made at the two time points, lizards from high elevation were 2.7°C more heat tolerant than lizards from low elevation (Table S2, Fig. S1C; Fig. 1A). Hematocrit decreased across time in lizards kept at the elevation of origin and was higher at the second time point in lizards transplanted to extreme high elevation than in lizards kept at their elevation of origin (Table S2, Fig. S1D; Fig. 1B). Hemoglobin concentration demonstrated similar, though non-significant, trends across time and in response to transplant to extreme high elevation (Table S2, Fig. S1E). Running endurance decreased during time in captivity for lizards kept at their elevation of origin, a result mostly driven by lizards originating from high elevation (Fig. S1G). Larger lizards had high RMR, but RMR did not differ among measurement groups (Table S2, Fig. S1H). Mass also positively influenced MMR and these measures differed among measurement groups. When measured at their elevation of origin and when combining all measurements made at both time points, low-elevation lizards exhibited higher MMR compared with high-elevation lizards (Table S2, Fig. S1I; Fig. 1C).

The NP-MANOVA with RRPP demonstrates that the multivariate thermal–metabolic phenotype clearly differs among measurement groups (Fig. 2; F_5,40=2.4, P=0.001). Lizards from both low- and high-elevation populations shifted phenotypes in response to captivity, even when held at their elevation of origin (low elevation: P=0.013; high elevation: P=0.003), driven primarily by decreases in both hematocrit levels and running endurance (Table 2). Transplanting lizards to extreme high elevation altered the thermal–metabolic phenotype in low-elevation lizards (P=0.046), but not high-elevation lizards (P=0.146), primarily driven by increases in hematocrit levels and hemoglobin concentration (Table 2). Multivariate phenotypes of lizards transplanted to extreme high elevation did not differ between lizards originating from low- or high-elevations (P=0.273; Table S3). Least-squares means of individual traits corroborate our univariate analyses (Fig. 3; Table S4).

We measured T_b in N=90 low-elevation and N=111 high-elevation lizards. Low-elevation lizards exhibited higher T_b compared with high-elevation lizards (32.5°C versus 29.9°C; t_193.4=5.56, P<0.001; Table S5; Fig. 4). However, thermoregulatory efficiency did not differ between the groups (Table S5). As values of d_b and d_e increase, both thermoregulatory accuracy and thermal quality of the environment decrease, respectively. For example, the mean d_e values between our high- and low-elevation sites only differed by 0.15, indicating relatively little difference in thermal environments available on the sampled days.

The interdependence of physiological traits, and their complex interaction with different environmental variables, challenges our ability to predict phenotypic responses in changing or novel environments and, by extension, the effects at a biogeographical scale (Beissinger and Riddell, 2021). Our study serves to quantify important spatial variation in multivariate trait expression and the plastic potential of these traits in a novel environment that may be soon colonized by a montane vertebrate. Whereas lizards from low- and high-elevation sites occupy distinct regions of trait space when measured at the elevation of origin, they converged in trait space under acute exposure to a novel oxygen environment (Fig. 2). Lizards from both low- and high-elevation populations exhibited plastic shifts in the multivariate thermal–metabolic phenotype in response to captive conditions, with low-elevation lizards changing more dramatically than high-elevation lizards in response to novel conditions of reduced oxygen availability (Figs 2 and 3). Aerobic capacity, an important determinant of behaviors and fitness (Gangloff and Telemeco, 2018; Seibel and Deutsch, 2020), is nonetheless limited in novel low-oxygen conditions (Table S2; Fig. 1C). This suggests that shifts in blood composition pertinent to oxygen transport – increases in hemoglobin concentration and hematocrit – only partially compensate for reduced oxygen availability. Describing this variation in phenotypic expression, especially the structure of trait correlation in the multidimensional thermal–metabolic phenotype, is a necessary first step in relating phenotypic expression to broader biogeographical patterns. Future work should be directed toward uncovering the genetic and physiological pathways that create these patterns of phenotypic variation and relating these underlying mechanisms to response potential.

To understand the potential for phenotypic change in response to novel conditions, we must first characterize existing physiological variation between populations. In our experiment, we found that critical thermal limits and thermal preferences were insensitive to variation in oxygen availability across an elevational gradient (Table S2). Nonetheless, we detected differences in CT_max and MMR between lizards measured at their elevations of origin. Our results of lower MMR but no effect of reduced oxygen availability in high-elevation lizards may indicate that the ‘starting point’ of physiological traits (i.e. due to local adaptation, developmental plasticity or other processes linked to population of origin that may lead to differentiation) may impact the trajectory of plastic response in these physiological phenotypes (West-Eberhard, 2003; Aubret and Shine, 2010; Phillips et al., 2016). In this case, individuals from the high-elevation populations exhibit MMR intermediate to those of lizards measured at low and extreme high elevation, suggesting that the variation observed may be the result of immediate oxygen environment. Thermal tolerance limits were not affected by the reduction of oxygen availability found at our extreme high-elevation treatment; correspondingly, these traits are unlikely to be oxygen limited as P. muralis tracks habitat upslope (Gangloff and Telemeco, 2018). Additionally, we found that high-elevation lizards are more heat tolerant than their low-elevation counterparts, a pattern documented in other lizard taxa (Llewelyn et al., 2016; Hodgson and Schwanz, 2019). The ‘hotter is better’ hypothesis supplies a compelling explanation for this counter-gradient pattern (Angilletta et al., 2010). Biochemical reactions accelerate as temperature increases (until an upper limit is reached) and, correspondingly, maximum performance, growth and metabolism should increase with higher body temperatures (Hamilton, 1973; Huey and Bennett, 1987). Counter-gradient variation (Conover and Schultz, 1995; Pettersen, 2020) in upper physiological limits, as we observed here, may be attributable to a greater basking need to compensate for limited thermal resources (Michniewicz and Aubret, 2010; Hodgson and Schwanz, 2019). This compensation may drive thermal physiological and behavioral evolution for increased thermal limits and preferences in cold climate reptiles relative to warm climate reptiles (Hertz and Huey, 1981; Huey et al., 2009; Llewelyn et al., 2017).

A wide variety of vertebrates, including humans, demonstrate shifts in hematocrit and hemoglobin concentration in response to short-term (days, weeks) exposure to hypoxic environments (Storz et al., 2010; Schlittler et al., 2020). This response increases oxygen diffusion and perfusion capacity by increasing the concentration of substrate in the blood to which oxygen can bind. This response, however, carries a physiological cost due to increased blood viscosity (Birchard, 1997; Dunlap, 2006). Previous work with P. muralis demonstrated that adult lowland lizards from both sexes increased hematocrit and hemoglobin concentration in response to both moderate and more extreme hypoxia (Gangloff et al., 2019; Kouyoumdjian et al., 2019). Notably, however, this shift was not sustained in the longer-term (6 weeks) and did not sufficiently compensate for performance decrements (Gangloff et al., 2019). Similarly, lizards in the present study exhibited plastic shifts in these blood parameters related to oxygen-carrying capacity, but in this case did not suffer decrement in running endurance or maximum metabolic rate at extreme high elevation. Further, running endurance and these blood parameters are correlated on the first axis of variation in the thermal–metabolic phenotype, suggesting a functional link. Plasticity in blood parameters, therefore, may be able to compensate for reduced oxygen availability, but only on relatively short time scales (several weeks) and, even then, only partially. This finding aligns with previous work showing that such short-term shifts in response to oxygen environments can be maladaptive in the long term (reviewed in Storz et al., 2010; Storz, 2021). In further support of this, we found no differences in hematrocrit or hemoglobin concentration between low- and high-elevation lizards, or in the plastic response of these measures. While this plastic response might be beneficial for individual animals under brief exposure, it is unlikely to facilitate range expansion and colonization of higher-elevation habitat (Hoffmann and Sgro, 2011; Munday et al., 2017; Catullo et al., 2019; Kelly, 2019). Nonetheless, hematocrit and running endurance acclimated in response to time in captivity, even at the elevation of origin (Table S2, Fig. S1F,G), similar to effects found previously with this species (Gangloff et al., 2019). Our findings shine a light on the importance of repeated measures and accounting for potential trait change of organisms while in captivity, even under consistent common garden conditions.

A multivariate approach to characterizing phenotypes can reveal important patterns of differentiation between populations and the response to specific environmental conditions. For example, Swaegers et al. (2020) described the transcriptomic response of damselflies (Ischnura elegans) from different latitudes to cool and warm temperatures. They found that the first and second axes of variation in whole-body RNA abundance described local adaptation and response to immediate temperature, respectively. Our multivariate approach similarly revealed important patterns not evident from the univariate analyses alone. Most clearly, we observed a separation of lizards from low- and high-elevation populations along the second axis of variation (Fig. 2). Importantly, these differences are consistent in measures made at the initial time point and those made after several weeks in captivity. This axis of variation most strongly contrasts lizards with different MMRs, consistent with our univariate results (Table S2; Fig. 1C). This difference in trait values may be due to genetically canalized differences, developmental differences in response to environmental cues, or the differences in the immediate conditions experienced during measurement. A reciprocal transplant experiment (measuring low-elevation lizards at high elevation and vice versa) is necessary to determine the relative influences of immediate environmental constraints and permanent mechanisms such as developmental plasticity or genetic differentiation. Selected temperature range also loads moderately high on this second axis of variation, potentially suggesting co-adaptation of thermoregulatory behavior and physiology, resulting in a narrowing of the TPC (Angilletta, 2009).

Whereas the second axis of variation in the multivariate phenotype clearly separates lizards originating from different elevations and measured at their elevation of origin, the first axis describes plastic responses to captivity. Most strongly, levels of hematocrit were reduced across time in captivity in animals kept at their native elevation. Running endurance also loaded moderately on the first axis, demonstrating that time in captivity reduces running endurance. In this case, changes in diet composition and water availability in captivity relative to the wild may have caused shifts in blood parameters, activity levels and performance, even as lizards were kept at the same elevation. Physiological phenotypes of lizards from both low- and high-elevation populations converged at extreme high elevation, suggesting a common response to this novel condition. Shifts in multivariate space in response to the novel oxygen environment occur along both the primary and secondary axes of variation. This belies simple interpretation of these axes as describing differences between populations (due to local adaptation or developmental differences) and within-individual shifts in traits (plasticity). That the multivariate phenotypes emerging in response to this novel condition did not differ in lizards originating from different elevations suggests that populations currently at high elevations – even those approaching the observed limit for the species in this area (Pottier, 2017) – are not necessarily more capable of dealing with these potential abiotic restrictions. This result in adult lizards parallels that observed in earlier life stages, specifically that developing embryos from low- and high-elevation P. muralis dams did not differ in their response to extreme high elevation (Kouyoumdjian et al., 2019). Future work directed toward determining whether the observed shifts in the thermal–metabolic phenotype are the result of immediate environmental constraints or active plasticity is crucial to determine fitness consequences and assess colonization potential (Forsman, 2015). Our results underscore the value of determining the patterns of variation in the multivariate phenotype among populations and potential plastic shifts in response to novel environmental conditions as an essential first step in predicting phenotypic responses and biogeographic shifts resulting from ongoing climate change.

The authors are grateful to J. Souchet, A. Josserand, E. Darnet and The Dø for field and logistical support. We thank the Observatoire Midi-Pyrénées staff at Pic du Midi and Y. Ardrit for their hospitality.