Mojave desert hydric physiology

Technical supplement

This is the analysis notebook for the working paper Phylogeny vs. physiology: how a guild of desert lizards maintain water balance. The project examines cutaneous evaporative water loss (CEWL) among a sample of lizards of 9 distinct species in the Mojave desert before and after acclimation to a controlled environment.

Comments appear like this.

To-do’s and placeholders for missing content appear like this.

Executive summary

This supplement presents descriptive and inferential analyses of cutaneous evaporative water loss (CEWL) among eight lizard species in the Mojave desert. CEWL was measured at capture — reflecting physiology embedded in field conditions — and after laboratory acclimation, reflecting intrinsic physiology with hydric stress removed. The difference in log(CEWL) between states quantifies physiological flexibility. Three parallel Bayesian phylogenetic mixed models were fit, one per phenotypic state, each modeling log(CEWL) as a function of blood plasma osmolality, ambient temperature, and vapor pressure deficit (VPD), with species-level random effects and phylogenetic covariance.

Key findings:

• Temperature and VPD credibly predict log(CEWL) at capture (temp: +19% per °C; VPD: −46% per kPa unit) but neither is credible at acclimation, suggesting these environmental factors influence CEWL under field conditions but not intrinsically.
• Osmolality is not credibly associated with log(CEWL) at any phenotypic state.
• Species-level variation in log(CEWL) is substantially larger at acclimation (\(\hat{\sigma}_s \approx 0.73\)) than at capture (\(\hat{\sigma}_s \approx 0.37\)), consistent with field conditions compressing apparent species differences.
• Phylogenetic signal (\(\lambda\)) is highly uncertain across all three models; posteriors span nearly the full [0, 1] interval and the data do not strongly support or refute phylogenetic structuring of any CEWL trait.
• Among species, only Chuckwallas showed a clearly elevated mean log(\(\Delta\)CEWL); flexibility estimates are otherwise imprecise and broadly consistent with zero change.

The supplement is organized as follows. The Descriptive Analysis section presents sampling information and exploratory summaries and graphics. The Statistical Modeling discusses model specification and presents summaries of fit and findings for each of the three phenotypic states: CEWL at capture, CEWL at acclimation, and physiological flexibility (log-scale change in CEWL from capture to acclimation). The Synthesis section compares estimates across all three models and examines phylogenetic signal.

Descriptive analysis

Primary measurements taken both at capture and after acclimation (2 days after capture) are:

• cewl: mean cutaneous evaporative water loss (in \(g/m^2h\)) measured in technical replicates
• osml: mean blood plasma osmolality (\(mmol/kg\)), or the concentration of solutes in the blood (proxy for hydraton)
• mass: mass (\(g\)) at time of measurement
• clo.temp: cloacal temperature (\(^\circ C\))
• temp: ambient temperature (\(^\circ C\)) at time of measurement
• vpd: ambient vapor pressure deficit (\(kPa\)) at time of measurement

Example rows are shown below.

Table: example data rows (continued below)

id	species	date	day	mass	osml	cewl	clo.temp
zt08	zt	2022-05-14	capture	15.1	356.5	6.48	25.5
zt08	zt	2022-05-16	acclimation	15.1	324	8.97	25.3
zt07	zt	2022-05-14	capture	14.7	394.3	6.66	26.84
zt07	zt	2022-05-16	acclimation	14.5	344.7	8.21	25.09
zt06	zt	2022-05-14	capture	11.2	374.3	8.85	26.94
zt06	zt	2022-05-16	acclimation	11	341.3	9.41	25.91

temp	vpd
25.12	2.29
25.7	2.25
25.34	2.29
25.68	2.27
25.04	2.26
26.87	2.52

For most species, about ten individuals were captured for the study, but sample sizes range from 2 to 24 individuals.

Table: sample sizes by species

sp.abbrv	sp.common	sp.sci	n
bg	western banded gecko	coleonyx variegatus	11
di	desert iguana	dipsosaurus dorsalis	24
cw	chuckwalla	sauromalus ater	14
lb	long-tailed brush lizard	urosaurus graciosus	10
sb	side-blotched lizard	uta stansburiana	17
tw	tiger whiptail	aspidoscelis tigris	10
zt	zebra-tailed lizard	callisaurus draconoides	7
mf	mojave fringe-toed lizard	uma scoparia	10
ln	long-nosed leopard lizard	gambelia wislizenii	2

Capture sites

Individuals were captured from 9 locales in the Mojave desert.

Figure: Map of study sites (blue) and capture locations (inset)

Capture locations were classified by microhabitat type. The capture site for individual TW10 was not classified.

Table: sample sizes by locale and microhabitat type

locale	floor	road	rock	sand	veg
cronese lake rd	-	-	-	11	1	-
desert studies center	-	1	-	2	-	-
horsethief springs	2	-	5	-	1	-
junk pile	1	-	-	-	-	-
kel-baker rd	-	8	-	-	-	-
kelso dunes	6	-	1	25	12	-
lava flow	2	-	21	-	-	1
salt creek hills	2	-	1	1	-	-
zzyzx rd	1	2	-	-	-	-

The microhabitat distribution of species sampled for the study is shown below. Most species were sampled from 1 or 2 primary microhabitats.

Figure: Microhabitat allocation of individuals sampled by species

Microhabitats are not considered in subsequent analysis because, although there are apparent marginal differences in CEWL and osmolality at capture when grouped by microhabitat, species is a confounding factor and after accounting for species the apparent marginal differences vanish, as shown in the figure below.

Figure: Distributions of CEWL and osmolality by microhabitat (top); Distributions of residual CEWL and osmolality by microhabitat after accounting for species (bottom).

Indeed, microhabitat does not explain a significant amount of variation in CEWL after adjusting for species.

Table: ANOVA of CEWL by species and microhabitat.

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
species	8	2169	271.2	19.96	1.284e-16
mh	4	6.441	1.61	0.1186	0.9756
species:mh	6	10.77	1.795	0.1321	0.9919
Residuals	85	1154	13.58	NA	NA

Moreover, there’s no obvious strong relationship between microhabitat temperatures and CEWL or osmolality at capture.

Figure: Microhabitat ambient temperature at capture site and CEWL and osmolality values after capture

Missing values

Several osmolality values are missing both at capture and after acclimation.

Table: counts of missing values

name	acclimation	capture
cewl	0	0
clo.temp	0	1
mass	0	0
osml	33	30
temp	0	0
vpd	0	0

The missing cloacal temperature value is from individual LB02. We’ll ignore this value since it was agreed on 8/20/25 that we’d use ambient temperatures rather than cloacal temperatures in statistical models.

The missing osmolality values are distributed across species as follows:

Table: counts of missing osmolality values by species

day	variable	bg	cw	di	lb	ln	mf	sb	tw	zt
acclimation	osml	11	5	1	8	0	0	8	0	0
capture	osml	11	1	0	7	0	0	11	0	0

Here is a list of the individuals with missing values:

• acclimation only: CW01, CW11, CW12, CW13, CW14, DI18, LB04
• capture only: CW02, SB01, SB02, SB08
• both: LB05, LB06, LB07, LB08, LB09, LB10, LB11, SB10, SB11, SB13, SB14, SB15, SB16, SB17, SB18

Missing osmolality measurements may become available for Chuckwallas at a later date per discussion on 8/20/25.

Sample sizes after removing missing values are rendered below.

Table: sample sizes by species after removing missing osmolality values

species	n	osml missing	osml recorded
bg	11	0	11
cw	14	8	6
di	24	23	1
lb	10	2	8
ln	2	2	0
mf	10	10	0
sb	17	6	11
tw	10	10	0
zt	7	7	0
Total	105	68	37

Cloacal and ambient temperatures are only moderately correlated (r ≈ 0.6 at capture, weaker after acclimation), and their ordering reverses between states: on average, cloacal temperature exceeds ambient at capture but falls below it after acclimation. Ambient temperature is used as the predictor in all models — it better captures the external hydric-thermal environment and avoids conflating the outcome physiology (thermoregulation affects cloacal temperature) with a predictor.

Distributions of primary measurements

The marginal distributions of individual measurements by species at capture and flexibility (log-scale change in CEWL) after acclimation are shown below.

Figure: Distributions of CEWL, osmolality, and mass by species and ambient measurement conditions at capture (top); distributions of physiological flexibility (log-scale change in CEWL) after acclimation (bottom).

At a glance, only Chuckwallas showed clearly elevated flexibility (increased CEWL after acclimation). Owing to the difference in scale of flexibility across species, CEWL is shown on the log scale.

Statistical modeling

A note on interpreting cross-sectional models. The capture and acclimation models describe cross-sectional patterns at each phenotypic state — they answer what predicts log(CEWL) at that moment, not what changed between states. Differences in estimates between these two models should be interpreted with care: they may reflect differences in measurement conditions and individual baseline variation rather than responses to acclimation per se. For inferences specifically about physiological flexibility (change in CEWL after acclimation), the flexibility model is the appropriate tool.

The framework for analysis is:

\[ \log(\text{CEWL}) \longleftarrow \text{osml} + \text{temp} + \text{vpd} + \text{species} \]

Individuals with missing osmolality measurements (discussed above) are removed from analysis, and covariates are centered for ease of interpretation when fitting models.

Phylogenetic mixed model specification

Denote the vector of CEWL measurements (or change in CEWL measurements) by \(Y\) and the covariate measurements (or change in covariate measurements) by:

\[ X = (1 \; \text{osml} \; \text{temp} \; \text{vpd}) \]

To enter species into the model, encode the species membership of each individual by the \(n \times S\) matrix:

\[ Z_S = \text{diag}(\mathbf{1}_{n_1}, \mathbf{1}_{n_2}, \dots, \mathbf{1}_{n_S}) \]

Here \(\mathbf{1}_{n_s}\) denotes a vector of ones of length \(n_s\), where \(n_s\) denotes the sample size for species \(s = 1, \dots, S\) (note \(\sum_s n_s = n\)).

Preliminary analysis using a model with uncorrelated random effects was used to test for species/osmolality interactions; no such interactions were significant at the 5% level, though there is suggestive evidence of an interaction for the acclimation and flexibility models. This interaction term is retained. Additional interactions were not identifiable and could not be tested.

Data: cdata
Models:
fit1_cap: log.cewl ~ osml + temp + vpd + (1 | species)
fit2_cap: log.cewl ~ osml + temp + vpd + (1 + osml | species)
         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
fit1_cap    6 37.121 51.026 -12.560    25.121                     
fit2_cap    8 37.442 55.982 -10.721    21.442 3.6788  2     0.1589

Data: adata
Models:
fit1_acc: log.cewl ~ osml + temp + vpd + (1 | species)
fit2_acc: log.cewl ~ osml + temp + vpd + (1 + osml | species)
         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
fit1_acc    6 91.699 105.36 -39.850    79.699                       
fit2_acc    8 89.778 107.99 -36.889    73.778 5.9214  2    0.05178 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Data: ddata
Models:
fit1_chg: log.cewl ~ osml + temp + vpd + (1 | species)
fit2_chg: log.cewl ~ osml + temp + vpd + (1 + osml | species)
         npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
fit1_chg    6 76.655 89.972 -32.327    64.655                       
fit2_chg    8 75.710 93.466 -29.855    59.710 4.9453  2    0.08436 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

To include an osmolality-species interaction term (see remark below), define \(Z_I = Z_S \circ \left(\mathbf{1}_S^\top \otimes x_1\right)\) and specify the linear mixed model:

\[ Y = X\beta + Z_S \gamma_S + Z_I \gamma_I + \epsilon \]

This is a random slope and random intercept model by species where \(\gamma_S\) captures species-specific differences in mean log(CEWL) and \(\gamma_I\) captures species-specific differences in the log(CEWL)–osmolality relationship. We assume \(\epsilon \sim N(0, \sigma^2 I_n)\) independent of the random effects, and allow correlated intercepts and slopes with within-species covariance:

\[ \Delta = \text{var}\left(\begin{array}{c}\gamma_{Sj} \\ \gamma_{Ij}\end{array}\right) = \left(\begin{array}{cc} \sigma_S^2 &\rho\sigma_S\sigma_I \\ \rho\sigma_S\sigma_I &\sigma_I^2\end{array}\right) \]

Phylogenetic relationships may induce trait dependence between species in proportion to the length of shared evolutionary history. Coleonyx variegatus (Western Banded Gecko) is excluded from all statistical models: its small body size makes blood draws impractical, so osmolality measurements are unavailable for this species. Under a Brownian motion model of evolutionary drift, the phylogeny induces the following expected correlations among an arbitrary continuous trait for the remaining eight species:

\[ V = \quad \begin{array}{l|cccccccc} & TW & DI & CW & LN & ZT & MF & SB & LB \\ \hline TW & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ DI & 0 & 1 & 0.688 & 0.504 & 0.504 & 0.504 & 0.504 & 0.504 \\ CW & 0 & 0.688 & 1 & 0.504 & 0.504 & 0.504 & 0.504 & 0.504 \\ LN & 0 & 0.504 & 0.504 & 1 & 0.562 & 0.562 & 0.562 & 0.562 \\ ZT & 0 & 0.504 & 0.504 & 0.562 & 1 & 0.847 & 0.712 & 0.712 \\ MF & 0 & 0.504 & 0.504 & 0.562 & 0.847 & 1 & 0.712 & 0.712 \\ SB & 0 & 0.504 & 0.504 & 0.562 & 0.712 & 0.712 & 1 & 0.755 \\ LB & 0 & 0.504 & 0.504 & 0.562 & 0.712 & 0.712 & 0.755 & 1 \end{array} \]

We incorporate a phylogenetic correlation structure for the random intercepts only, so that \(\gamma_S \sim N(0, \sigma^2_S V_\lambda)\) and \(\gamma_I \sim N(0, \sigma^2_I I_S)\), where \(V_\lambda = (1 - \lambda)I_S + \lambda V\) and \(\lambda \in (0, 1)\) captures the strength of phylogenetic signal in the trait after adjusting for covariates. Within-species slope–intercept correlation \(\rho = \text{Corr}(\gamma_{Sj}, \gamma_{Ij})\) is retained. Together these assumptions give the marginal covariance structure:

\[ \Sigma = \text{var}Y = \sigma^2_S Z_S V_\lambda Z_S^T + \rho\sigma_S\sigma_I\left(Z_S Z_I^T + Z_I Z_S^T\right) + \sigma^2_I Z_I Z_I^T + \sigma^2 I_n \]

The model is estimated in a Bayesian framework with normal priors for regression coefficients and half-normal priors for variance parameters. We used the one-to-one reparametrization \(\tau_S^2 = \sigma_S^2(1-\lambda)\) and \(\tau_P^2 = \sigma_S^2\lambda\), which expands the intercept covariance term and gives:

\[ \Sigma = \text{var}Y = \tau^2_S Z_S Z_S^T + \tau^2_P Z_S V Z_S^T + \rho\sigma_S\sigma_I\left(Z_S Z_I^T + Z_I Z_S^T\right) + \sigma^2_I Z_I Z_I^T + \sigma^2 I_n \]

Half-normal priors (specified as \(\text{Normal}(0, \sigma)\) with a lower bound of zero) were used for all variance parameters. Exponential priors were avoided because their sharp spike at zero creates a funnel-shaped posterior geometry when variance components are small, which produces divergent transitions and poor sampler efficiency. All three models share log(CEWL) as the outcome, so prior scales were calibrated to the log scale throughout. For the capture and acclimation models, the log(CEWL) outcome has mean \(\approx 2.4\) and SD \(\approx 0.5\)–\(0.8\), giving fixed-effect priors \(\text{Normal}(0, 2)\) for the intercept and \(\text{Normal}(0, 1)\) for slopes, and variance priors \(\tau_S \sim \text{HN}(0, 1)\), \(\tau_P \sim \text{HN}(0, 0.5)\), and \(\sigma \sim \text{HN}(0, 1)\). The flexibility model outcome (log-scale CEWL change, i.e., log(\(\Delta\)CEWL)) has mean \(\approx 0\) and SD \(\approx 0.6\), giving \(\text{Normal}(0, 1)\) for the intercept, \(\text{Normal}(0, 2)\) for slopes (to accommodate larger VPD effects), and \(\text{HN}(0, 0.5)\) for both \(\tau_S\) and \(\tau_P\). Note that \(\rho\) applies only to the iid intercept component \(\tau_S\); the phylogenetic component \(\tau_P\) is assumed independent of the slope.

Analysis of capture data

Sample size: n = 75 individuals across 8 species (after removing missing osmolality values).

Model diagnostic checks

Diagnostic MCMC trace and density plots show good chain mixing and consistent posteriors, indicating successful convergence.

Figure: Trace (left) and posterior density (right) plots for each model parameter.

Figure: Pairs plots for variance parameters.

Coefficient estimates

Guide to reading these results. All estimates are on the log(CEWL) scale; to convert to a multiplicative effect on raw CEWL, exponentiate: a slope of \(\beta\) means a 1-unit increase in the predictor is associated with an \(e^\beta\)-fold change in CEWL. For example, the capture model temperature slope of \(\approx 0.175\) means each 1°C increase is associated with \(e^{0.175} \approx 1.19\), or roughly a 19% increase in CEWL. The 95% credible interval (CI) is the range of values that the posterior assigns 95% probability to — if the CI excludes zero, the association is considered credible. Key variance parameters:

• \(\tau_S\) (species intercept iid SD): typical spread in baseline log(CEWL) across species, after removing phylogenetic structure
• \(\tau_P\) (phylogenetic SD): how much of species-level variation is attributable to shared evolutionary history; feeds into \(\lambda\) below
• \(\lambda\) (phylogenetic signal): proportion of total species variation explained by phylogeny (0 = no phylogenetic structure; 1 = trait perfectly tracks the phylogenetic tree)
• \(\sigma_I\) (slope SD): how much the osmolality–log(CEWL) slope varies across species
• \(\rho\) (intercept–slope correlation): whether species with higher baseline log(CEWL) also tend to have stronger (or weaker) osmolality responses
• \(\sigma_\epsilon\) (residual SD): unexplained individual-level variation after accounting for species and phylogeny

The table below shows parameter estimates.

Table: parameter estimates and 95% credible intervals

term	estimate	lwr	upr
Intercept	2.357	1.698	2.932
osml	0.002411	-0.002213	0.006907
temp	0.175	0.09579	0.2542
vpd	-0.6089	-1.045	-0.1613
lambda	0.4239	0.001234	0.995
sigma_S	0.5428	0.2779	0.9822
sigma_I	0.003911	0.0003943	0.01086
rho	-0.3305	-0.9655	0.7088
sigma_eps	0.2496	0.2093	0.3003

The figure below shows the posterior densities for each parameter from the MCMC draws, and two posterior predictive checks (PPCs): a scatter of individual predicted versus observed means, and a density overlay comparing replicated datasets from the posterior to the observed distribution. Together these check both mean-level accuracy and distributional shape.

Figure: Posterior densities (top); posterior predictive checks — predicted vs. observed means (bottom left) and posterior predictive density overlay (bottom right).

Owing to the random slopes and intercepts, each species has a distinct estimated baseline log(CEWL) when covariates are fixed at their means, and a distinct slope in osmolality. The table below shows the species-specific estimates.

Table: species-specific relationships after adjusting for phylogeny

	Intercept	osml
cw	2.83	0.0008697
di	2.092	0.005056
lb	2.486	0.002066
ln	2.24	0.002639
mf	2.186	0.00447
sb	2.429	0.002701
tw	2.417	0.0005989
zt	2.269	0.0006366

Displayed visually:

Figure: (Left) Species-specific intercepts and osmolality slopes (95% CI); dashed line shows population-average. (Right) Implied species-specific osmolality–log(CEWL) relationships at average temperature and VPD.

Summary of findings

• At mean conditions, mean log(CEWL) is estimated to be 2.36 (95% CI: 1.70, 2.93), corresponding to approximately 5.5–18.8 \(g/m^2h\) on the raw scale
• There is no credible association between log(CEWL) and osmolality (95% CI: −0.002 to 0.007)
• Temperature is credibly and positively associated with log(CEWL): each 1°C increase is associated with \(e^{0.175} \approx 1.19\)-fold higher CEWL (95% CI on slope: 0.096 to 0.254)
• VPD is credibly and negatively associated with log(CEWL) (95% CI on slope: −1.045 to −0.161), consistent with higher ambient evaporative demand reducing CEWL
• Species-specific variation in log(CEWL) is estimated at \(\hat{\sigma}_s\) = 0.37 (95% CI: 0.04, 0.80)
• Phylogenetic signal is highly uncertain (\(\lambda\): median 0.36, 95% CI ≈ 0 to 1.0); the posterior spans nearly the full unit interval and the data do not meaningfully constrain \(\lambda\)
• The slope/intercept correlation is imprecisely estimated (95% CI: −0.965 to 0.709), with a weakly negative central tendency suggesting species with higher baseline log(CEWL) may exhibit a shallower osmolality response

Analysis of acclimation data

Sample size: n = 72 individuals across 8 species (after removing missing osmolality values).

Model diagnostic checks

Diagnostic MCMC trace and density plots show good chain mixing and consistent posteriors, indicating successful convergence.

Figure: Trace (left) and posterior density (right) plots for each model parameter.

Pairs plots of the posterior draws for the variance parameters.

Figure: Pairs plots for variance parameters.

Coefficient estimates

Parameters are interpreted as in the capture section. The intercept here represents log(CEWL) in the absence of hydric stress and field ecological variation. A larger \(\sigma_S\) relative to the capture model would indicate that acclimation reveals greater between-species heterogeneity than is visible under field conditions.

The table below shows parameter estimates.

Table: parameter estimates and 95% credible intervals

term	estimate	lwr	upr
Intercept	2.374	1.424	3.222
osml	0.0006288	-0.005686	0.006696
temp	-0.001903	-0.1517	0.1445
vpd	0.2223	-0.5105	0.9406
lambda	0.3124	0.0005509	0.9795
sigma_S	0.921	0.5422	1.515
sigma_I	0.005362	0.001124	0.01398
rho	-0.3609	-0.962	0.5926
sigma_eps	0.3433	0.2865	0.4151

The figure below shows the posterior densities for each parameter and two posterior predictive checks.

Figure: Posterior densities (top); posterior predictive checks — predicted vs. observed means (bottom left) and posterior predictive density overlay (bottom right).

Owing to the random slopes and intercepts, each species has a distinct estimated baseline log(CEWL) at acclimation and a distinct slope in osmolality. The table below shows the species-specific estimates.

Table: species-specific relationships after adjusting for phylogeny

	Intercept	osml
cw	3.611	-0.003224
di	2.007	0.004936
lb	2.405	0.0007533
ln	2.372	0.0004921
mf	2.065	-0.001305
sb	2.579	0.000583
tw	2.393	-0.0001061
zt	1.907	0.002243

Displayed visually:

Figure: (Left) Species-specific intercepts and osmolality slopes (95% CI); dashed line shows population-average. (Right) Implied species-specific osmolality–log(CEWL) relationships at average temperature and VPD.

Summary of findings

• At mean conditions, mean log(CEWL) at acclimation is estimated to be 2.37 (95% CI: 1.42, 3.22)
• There is no credible association between log(CEWL) and osmolality (95% CI: −0.006 to 0.007)
• Neither temperature nor VPD is credibly associated with log(CEWL) at acclimation (temp 95% CI: −0.152 to 0.144; VPD 95% CI: −0.511 to 0.941) — in contrast to the capture model
• Species-specific variation in log(CEWL) is larger at acclimation (\(\hat{\sigma}_s\) = 0.73, 95% CI: 0.14, 1.41) than at capture (\(\hat{\sigma}_s\) = 0.37), suggesting field conditions partially compress apparent species differences
• Phylogenetic signal is highly uncertain (\(\lambda\): median 0.21, 95% CI ≈ 0 to 0.98)
• The slope/intercept correlation is imprecisely estimated (95% CI: −0.962 to 0.593), with a weakly negative central tendency

Analysis of flexibility

Sample size: n = 68 individuals across 8 species with complete osmolality measurements at both time points (used to compute log-scale flexibility).

Model diagnostic checks

Diagnostic MCMC trace and density plots show good chain mixing and consistent posteriors, indicating successful convergence.

Figure: Trace (left) and posterior density (right) plots for each model parameter.

Pairs plots of the posterior draws for the variance parameters.

Figure: Pairs plots for variance parameters.

Coefficient estimates

Here the outcome is log(ΔCEWL) — physiological flexibility — the log-scale change in CEWL from capture to acclimation. All three covariates are likewise computed as acclimation-minus-capture differences (Δosmolality, Δtemperature, ΔVPD), so slopes represent the association between a unit change in the covariate difference and log(ΔCEWL). A positive slope for a predictor means that individuals whose conditions changed more in that direction showed greater flexibility. The intercept represents expected mean flexibility at mean conditions. \(\lambda\) here quantifies how phylogenetically structured flexibility is, not CEWL level; a high \(\lambda\) would indicate closely related species tend to show similar magnitudes of flexibility.

The table below shows parameter estimates.

Table: parameter estimates and 95% credible intervals

term	estimate	lwr	upr
Intercept	0.1316	-0.507	0.7681
osml	-0.0003449	-0.00567	0.00623
temp	0.2608	0.1691	0.3528
vpd	-1.269	-1.871	-0.6681
lambda	0.4446	0.001027	0.9975
sigma_S	0.6076	0.337	1.022
sigma_I	0.004889	0.0004437	0.01269
rho	-0.3914	-0.9705	0.7319
sigma_eps	0.3393	0.2805	0.4119

The figure below shows the posterior densities for each parameter and two posterior predictive checks.

Figure: Posterior densities (top); posterior predictive checks — predicted vs. observed means (bottom left) and posterior predictive density overlay (bottom right).

Owing to the random slopes and intercepts, each species has a distinct estimated mean flexibility (log(ΔCEWL)) and a distinct slope in osmolality. The table below shows the species-specific estimates.

Table: species-specific relationships after adjusting for phylogeny

	Intercept	osml
cw	0.7798	-0.005721
di	0.02576	-0.002349
lb	0.06994	0.000662
ln	0.1485	1.946e-05
mf	-0.02002	0.0006076
sb	0.1839	0.0006314
tw	0.03236	0.001856
zt	-0.1576	0.001959

Displayed visually:

Figure: (Left) Species-specific intercepts and osmolality slopes (95% CI); dashed line shows population-average. (Right) Implied species-specific osmolality-log(ΔCEWL) relationships at average temperature and VPD.

The species-level intercepts from the flexibility model directly estimate mean log(\(\Delta\)CEWL) for each species — i.e., which species tend to show greater or lesser flexibility after acclimation, after accounting for measurement conditions. The plot below ranks species by estimated flexibility.

Figure: Species-specific mean log(ΔCEWL) with 95% credible intervals, ranked by estimate. Dashed red line shows the population average; dashed black line at zero indicates no change.

Summary of findings

• At mean conditions, mean log(\(\Delta\)CEWL) is estimated at 0.13 (95% CI: −0.51, 0.77); Chuckwallas are the clearest outlier with a species-level estimate of ~0.78
• There is no credible association between log(\(\Delta\)CEWL) and osmolality (95% CI: −0.006 to 0.006)
• Temperature is credibly and positively associated with log(\(\Delta\)CEWL): each 1°C increase is associated with \(e^{0.261} \approx 1.30\)-fold greater flexibility (95% CI on slope: 0.169 to 0.353)
• VPD is credibly and negatively associated with log(\(\Delta\)CEWL) (95% CI: −1.871 to −0.668), suggesting higher evaporative demand at measurement suppresses the acclimation response
• Species-specific variation in flexibility is uncertain (\(\hat{\sigma}_s\) = 0.40, 95% CI: 0.04, 0.79)
• Phylogenetic signal in flexibility is highly uncertain (\(\lambda\): median 0.40, 95% CI ≈ 0 to 1.0); with only eight species the data cannot meaningfully constrain \(\lambda\)
• The slope/intercept correlation is imprecisely estimated (95% CI: −0.971 to 0.732)

Synthesis

Cross-model parameter comparison

The tables below compare estimates across phenotypic states. The first places capture and acclimation side by side (both on the log(CEWL) scale); the second shows the flexibility model estimates (log(\(\Delta\)CEWL) scale). Cells show posterior mean and 95% credible interval.

Table: Posterior mean and 95% CI — capture and acclimation models (outcome: log(CEWL)).

Parameter	Capture	Acclimation
Intercept	2.357 (1.698, 2.932)	2.374 (1.424, 3.222)
Osmolality	0.002 (-0.002, 0.007)	0.001 (-0.006, 0.007)
Temperature	0.175 (0.096, 0.254)	-0.002 (-0.152, 0.144)
VPD	-0.609 (-1.045, -0.161)	0.222 (-0.511, 0.941)
Species SD	0.367 (0.042, 0.803)	0.728 (0.138, 1.408)
Phylogenetic SD	0.326 (0.016, 0.835)	0.444 (0.02, 1.115)
Slope SD	0.004 (0, 0.011)	0.005 (0.001, 0.014)
Intercept-slope corr.	-0.331 (-0.965, 0.709)	-0.361 (-0.962, 0.593)
Residual SD	0.25 (0.209, 0.3)	0.343 (0.287, 0.415)
Phylogenetic signal	0.424 (0.001, 0.995)	0.312 (0.001, 0.98)

Table: Posterior mean and 95% CI — flexibility model (outcome: log(ΔCEWL)).

Parameter	Flexibility (log(ΔCEWL))
Intercept	0.132 (-0.507, 0.768)
Osmolality	0 (-0.006, 0.006)
Temperature	0.261 (0.169, 0.353)
VPD	-1.269 (-1.871, -0.668)
Species SD	0.397 (0.036, 0.793)
Phylogenetic SD	0.376 (0.016, 0.92)
Slope SD	0.005 (0, 0.013)
Intercept-slope corr.	-0.391 (-0.971, 0.732)
Residual SD	0.339 (0.281, 0.412)
Phylogenetic signal	0.445 (0.001, 0.997)

Phylogenetic signal across phenotypic states

The figure below shows the full posterior distribution of \(\lambda\) from each model, making visible both the central tendency and the considerable uncertainty across all three phenotypic states.

Figure: Posterior distributions of phylogenetic signal (λ) for each phenotypic state. All three posteriors are broad, indicating the data are largely uninformative about the degree of phylogenetic structuring in any CEWL trait.

Summary

Temperature and VPD credibly predict log(CEWL) at capture — each 1°C increase is associated with a ~19% increase in CEWL, and higher VPD is associated with ~46% lower CEWL — but neither predictor is credible at acclimation or for flexibility. This pattern is consistent with field conditions driving CEWL variation at capture through behavioral thermoregulation and microclimate differences, while intrinsic CEWL (measured under controlled acclimation) varies primarily among species rather than tracking environmental conditions. Osmolality shows no credible association with CEWL at any state, suggesting hydration status (at the range observed here) may not be a strong proximal driver of cutaneous water loss.

Species-level variation in log(CEWL) is nearly twice as large at acclimation as at capture, pointing to field conditions obscuring inherent between-species differences — consistent with the interpretation that capture CEWL reflects both physiology and ecology, while acclimated CEWL is a cleaner physiological signal. Phylogenetic signal is highly uncertain across all three phenotypic states, with posteriors that are nearly flat across the full 0–1 interval; with only eight species, the data lack the resolution to reliably partition species-level variation into phylogenetic and non-phylogenetic components. Among species, Chuckwallas stand out with the highest estimated flexibility (mean log(\(\Delta\)CEWL) ≈ 0.78), though even this estimate carries a credible interval that includes zero.