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.

The primary research question is:

To what extent do osmolality, ambient temperature, vapor pressure deficit, body mass, species differences, and shared phylogeny explain variation in cutaneous evaporative water loss at capture and after acclimation?

Comments appear like this.

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

Descriptive summaries

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.

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 after 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

Which temperature to use?

It is worth noting that cloacal and ambient temperatures are related. Some exploration was carried out in consideration of whether both should factor in the analysis.

Figure: Cloacal and ambient temperatures at capture and after acclimation (top); changes in value after acclimation (bottom)

Contrary to expectation: (a) the two are only moderately correlated and (b) the correlation is weaker after acclimation, as shown in the top panel. This also makes visible that the ordering of temperature measurements reverses after acclimation – on average, cloacal temperatures are higher than ambient temperatures at capture but lower than ambient temperatures after acclimation. The bottom panels show that while ambient temperatures did not appear to change systematically after acclimation, cloacal temperatures decreased on average.

Per discussion on 8/20/25, ambient temperature is of greater interest as a possible predictor of cutaneous evaporative water loss. We may wish to revisit this at a later date.

Distributions of primary measurements

The marginal distributions of individual measurements by species at capture and changes after acclimation are shown below.

Figure: Distributions of CEWL, osmolality, and mass by species and ambient measurement conditions at capture (top); distributions of changes after acclimation (bottom).

At a glance, only Chuckwallas showed a dramatic change (increase, in this case) in CEWL due to acclimation, yet they experienced little change in mass or osmolality (though the latter changes exhibit a relatively large variance). Owing to the difference in scale of change among Chuckwallas compared with other species, change in CEWL is shown above (and later modeled) on the log scale.

Other observations regarding marginal distributions of primary measruements:

• ambient conditions at times of measurement are typically 24-26 degrees C and 2.0-2.4 kPa
• a few species (LN, TW, LT, ZT) exhibited changes in osmolality, but they were of relatively small magnitude (20-40 mmol/kg)
• ambient temperatures fluctuated by 1-2 degrees between measurements taken after capture and measurements taken after acclimation
• vapor pressure deficits fluctuated by 0.2-0.3 kPa between measurements taken after capture and measurements taken after acclimation
• changes in mass are negligible

Statistical modeling

The framework for analysis is:

\[ \text{CEWL} \longleftarrow \text{osml} + \text{mass} + \text{temp} + \text{vpd} \]

Marginal relationships between the primary covariates and CEWL are shown below. Mass is log-transformed owing to large species differences in typical body size.

Figure: Marginal relationships between change in CEWL and changes in other primary measurements; logratio changes in CEWL shown at top; and arithmetic difference in CEWL shown at bottom

Considering that changes in mass are negligible, do we want to consider instead compare change in CEWL directly with mass at capture?

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

Bayesian phylogenetic mixed model

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 \; \log(\text{mass}) \; \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 = \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\)).

Assuming the relationship between CEWL and each covariate does not depend on species yields the linear mixed model:

\[ Y = X\beta + Z \gamma + \epsilon \]

Here the random vector \(\gamma\) captures species differences in mean CEWL. Species-specific intercepts are \(\beta_0 + \gamma_s\) for species \(s = 1, \dots, S\). The random vector \(\epsilon \sim N(0, \sigma^2 I_n)\) represents individual variation, and it is assumed that \(\epsilon \perp \gamma\).

Preliminary analysis using a model with uncorrelated random effects was used to test for species/covariate interactions; no such interactions were significant. These results are not shown here, but codes to reproduce the inference are included in the source code for this notebook.

If \(\gamma \sim N(0, \sigma^2_\gamma I_S)\) species random effects are uncorrelated, but phylogenetic relationships may induce trait dependence between species in proportion to the length of shared evolutionary history. For the species in this study, an ultrametric (i.e., time-calibrated) tree was constructed from the phylogeny in Pyron et al. (2013) and is shown below.

Under a Brownian motion model of evolutionary drift, this structure induces the following expected correlations among an arbitrary continuous trait:

\[ 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.696 & 0.528 & 0.528 & 0.528 & 0.528 & 0.528 \\ cw & 0 & 0.696 & 1 & 0.528 & 0.528 & 0.528 & 0.528 & 0.528 \\ ln & 0 & 0.528 & 0.528 & 1 & 0.59 & 0.59 & 0.59 & 0.59 \\ zt & 0 & 0.528 & 0.528 & 0.59 & 1 & 0.864 & 0.74 & 0.74 \\ mf & 0 & 0.528 & 0.528 & 0.59 & 0.864 & 1 & 0.74 & 0.74 \\ sb & 0 & 0.528 & 0.528 & 0.59 & 0.74 & 0.74 & 1 & 0.781 \\ lb & 0 & 0.528 & 0.528 & 0.59 & 0.74 & 0.74 & 0.781 & 1 \end{array} \]

To account for phylogeny, the covariance structure of the random effects is specified with the phylogenetic correlation matrix \(V\) shown above so that \(\gamma \sim N(0, \sigma^2_\gamma V_\lambda)\) where \(V_\lambda = (1 - \lambda)I_S + \lambda V\). This gives the marginal covariance structure:

\[ \Sigma = \text{var}Y = \sigma^2_\gamma Z V_\lambda Z^T + \sigma^2 I_n \]

The model is estimated in a Bayesian framework with normal priors for regression coefficients and exponential priors for variance parameters. To estimate \(\sigma_\gamma\) and \(\lambda\), the variance structure was reparametrized:

\[ \Sigma = \text{var}Y = \sigma^2_S ZZ^T + \sigma^2_P Z V Z^T + \sigma^2 I_n \]

Tighter priors were needed to ensure that \(\sigma^2_S, \sigma^2_P\) were well-identified by the data; a weak prior was used for \(\sigma^2\).

Model fit to capture data

Table: parameter estimates and 95% credible intervals

term	estimate	lwr	upr
Intercept	12.68	10.6	14.74
log.mass	2.072	0.6908	3.472
osml	0.008218	-0.0145	0.03077
temp	2.495	1.551	3.4
vpd	-10.33	-15.55	-5.222
lambda	0.016	4.634e-06	0.1154
sigma_g	2.461	1.749	3.363
sigma	2.98	2.507	3.574

Figure: Posterior densities (left); observed vs. predicted CEWL values (right).

Findings at capture:

• There are significant positive associations between CEWL and each of (log) mass and ambient temperature
• There is a significant negative association between CEWL and vapor pressure deficit
• There is a significant species effect
• There is roughly as much estimated interspcies as intraspecies variation in CEWL
• Phylogenetic signal is negligible (\(\lambda \in (0, 0.12)\))

Owing to the random intercepts, each species has a distinct estimated baseline CEWL when covariates are fixed at their means. The table below shows the species-specific estimates.

Table: species-specific relationships after adjusting for phylogeny

	Intercept	log.mass	osml	temp	vpd
cw	18.16	2.072	0.008218	2.495	-10.33
di	7.833	2.072	0.008218	2.495	-10.33
lb	15.82	2.072	0.008218	2.495	-10.33
ln	10.22	2.072	0.008218	2.495	-10.33
mf	9.716	2.072	0.008218	2.495	-10.33
sb	16.11	2.072	0.008218	2.495	-10.33
tw	12.67	2.072	0.008218	2.495	-10.33
zt	11.81	2.072	0.008218	2.495	-10.33

The plot below shows population-level marginal trends in each variable according to the model.

Figure: Posterior means with 95% credible intervals in each variable while holding other terms fixed at population averages.

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

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

Pairs plots of the posterior draws for the variance parameters indicate that the phylogenetic signal is reasonably well-identified by the data, though there may still be a slight ridge.

Figure: Pairs plots for variance parameters.

Model fit to change after acclimation

Under development

Table: parameter estimates and 95% credible intervals

term	estimate	lwr	upr
Intercept	0.1221	-0.4908	0.7177
log.mass	0.6008	-2.139	3.33
osml	-0.002425	-0.005564	0.000749
temp	0.2478	0.158	0.3389
vpd	-1.38	-1.988	-0.7722
lambda	0.3791	0.0002392	0.9986
sigma_g	0.5301	0.3	0.9229
sigma	0.3549	0.2958	0.4304

Figure: Posterior densities (left); observed vs. predicted CEWL values (right).

Findings: …

Owing to the random intercepts, each species has a distinct estimated baseline CEWL when covariates are fixed at their means. The table below shows the species-specific estimates.

Table: species-specific relationships after adjusting for phylogeny

	Intercept	log.mass	osml	temp	vpd
cw	0.7887	0.6008	-0.002425	0.2478	-1.38
di	-0.0289	0.6008	-0.002425	0.2478	-1.38
lb	0.1215	0.6008	-0.002425	0.2478	-1.38
ln	0.123	0.6008	-0.002425	0.2478	-1.38
mf	-0.03077	0.6008	-0.002425	0.2478	-1.38
sb	0.2144	0.6008	-0.002425	0.2478	-1.38
tw	0.0067	0.6008	-0.002425	0.2478	-1.38
zt	-0.2233	0.6008	-0.002425	0.2478	-1.38

The plot below shows population-level marginal trends according to the model.

Figure: Posterior means with 95% credible intervals in each variable while holding other terms fixed at population averages.

Diagnostic MCMC trace and density plots show…

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

Pairs plots of the posterior draws for the variance parameters suggest…

Figure: Pairs plots for variance parameters.