Statistical analyses
We tested for differences in lizard foraging body temperatures between
snake and snake-free islands as well as between the two periods of
1981-1991 and 2018-2019 by building a linear mixed-effects model using
the “lme4” package (De Boeck et al. 2011) in R (R Core Team 2019),
with air temperature in the shade as a covariate, snake presence and
period as fixed effects, and island, month and year as random
intercepts. We used the “MuMIn” package (Bartoń 2020) to calculate
R2 values for the model.
To compare the thermal dependency of maximal locomotory speeds between
Hachijo-Kojima and Kozu P. latiscutatus on the racetrack, and
between E. quadrivirgata and Kozu P. latiscutatusin open field, we modelled thermal performance curves (Angilletta 2009).
For this, we fit raw data for each group to four different non-linear
models (Tables S1 and S2 in Supporting Information; Angilletta 2006). We
then chose the best fitting model overall in each of the two comparisons
based on yielded Akaike’s information criterion (AIC), assuming a
difference larger than 2 to be significant (Burnham and Anderson 2002),
relative AIC weights (wi ), residual sums of
squares (RSS), as well as the biological relevance of yielded curve
shapes, where an expected shape would see a gradual increase in
performance with body temperature up until an optimum after which
performance drops rapidly to zero (Angilletta 2006; Angilletta 2009).
From the thermal performance curves, we extracted values pertinent to
the characterization of the thermal dependency of locomotion (Logan et
al. 2013). These included the thermal optimum for sprinting
(Topt ), maximal performance
(Pmax ), and the range of body temperatures at
which 95% of Pmax is attainable. The latter
estimate is termed the B95 performance breadth,
and is interpreted as the thermal optimum range (Hertz et al. 1983).
In addition to thermal physiology, we aimed to determine if lizards from
Kozu, where the snake predator is present, also show adaptations in limb
morphology relative to those from Hachijo-Kojima. This is a relevant
morphological dimension to compare considering that lizards that have
evolved longer hindlimbs run faster at both population and species level
(Bonine and Garland 1999; Aerts et al. 2000). We thus tested for
differences in hind leg and hind foot lengths, both as ratios of SVL,
between lizards from the two islands with two-sampled t-tests. We then
applied a Bayesian model to determine relationships between
morphological traits and maximum sprint speeds for both lizard
populations while accounting for multi-collinearity in morphological
traits. Here we used individual maximum sprint speed measurements rather
than values extracted from modelled thermal performance curves. To
remove the effect of temperature on sprint speed, we estimated relative
sprint speed by regressing log-transformed body temperature against
log-transformed sprint speed (Lleonart et al. 2000), where we used the
standard major axis regression as the regression. We subsequently used
the obtained residuals as estimates of relative sprint speed. We
conducted the analysis using the “smart” package in R (Warton et al.
2012; R Core Team 2019). We constructed a hierarchical model consisting
of four equations to estimate the contribution of each measurement to
relative sprint speed (Zipkin et al. 2020; Figure S2). We defined the
core of this model with the following equation:
\begin{equation}
\begin{matrix}\text{rS}P_{i}\ \sim\ Normal\left(\beta_{1,k}\bullet\text{HL}_{i}+\beta_{2,k}\bullet\text{HF}_{i}+\beta_{3,k}\bullet\text{Tail}_{i}+\beta_{4,k}\bullet\text{SVL}_{i}+ep_{1,k},\ \text{\ σ}_{1}^{2}\right)\#(1\\
\end{matrix})\nonumber \\
\end{equation}\begin{equation}
i=1,2,\ldots N,\ k=1,2\nonumber \\
\end{equation}where \(\text{rS}P_{i}\) represents the value of the \(i\)th relative
sprint speed, \(HL_{i}\) represents the value of \(i\)th hind leg
length, \(HF_{i}\) represents the value of \(i\)th hind foot length,\(\text{Tai}l_{i}\) represents the value of \(i\)th tail length,\(\text{SV}L_{i}\) represent the value of \(i\)th SVL,\(\beta_{1,k},\ \beta_{2,k},\ldots,\ \beta_{4,k}\) are the coefficients
of each variable, \(ep_{1,k}\) is the fixed effect on the k th
population, \(\sigma_{1}^{2}\) is the variance, \(N\) is the number of
individuals, \(k\) is the index for the island (Hachijo-Kojima or Kozu).
This model assumes that the contribution of each measurement to relative
sprint speed varies among populations. Since the autocorrelation between
each measurement should be considered only by this model, we defined the
following equations to consider these correlations:
\begin{equation}
\begin{matrix}\text{HL}_{i}\ \sim\ Normal\left(\beta_{5,k}\bullet SVL_{i}+ep_{2,k},\ \sigma_{2}^{2}\right)\ \#\left(2\right)\\
\end{matrix}\nonumber \\
\end{equation}\begin{equation}
\begin{matrix}\text{HF}_{i}\sim Normal\left(\beta_{6,k}\bullet SVL_{i}+ep_{3,k},\ \sigma_{3}^{2}\right)\ \#\left(3\right)\\
\end{matrix}\nonumber \\
\end{equation}\begin{equation}
\begin{matrix}\text{Tail}_{i}\sim Normal\left(\beta_{7,k}\bullet SVL_{i}+ep_{4,k},\ \sigma_{4}^{2}\right)\ \#\left(4\right)\\
\end{matrix}\nonumber \\
\end{equation}where \(\beta_{5,k},\ \beta_{6,k},\) and \(\beta_{7,k}\) are the
coefficients of SVL, \(ep_{2,k},\ ep_{3,k}\), and \(ep_{4,k}\) are fixed
effects, and\(\ \sigma_{2}^{2},\ \sigma_{3}^{2}\),\(\ \sigma_{4}^{2}\)are the variances showing individual heterogeneity.
To uncover the functional links between morphology, performance and
fitness (Arnold 1983) for lizards, we estimated selection gradients and
shapes of natural selection on hind leg length. We determined if lizards
with longer hind legs showed higher rates of survival from predation
using tail break distance as a proxy for predator escape success. We
assumed individuals with longer tail break distances relative to SVL to
be more capable of escaping predatory events, and therefore more likely
to survive. We thus examined the survival component of fitness related
to hind leg length and sprint speed by comparing tail break distance as
a function of hind leg length for two snake-free islands (Miyake and
Hachijo-Kojima) and two islands with snakes (Kozu and Mikura). To
correct for the effect of SVL on tail break distance and hind leg
length, we calculated their relative values
(\(\text{rTBD}\ \text{and}\text{\ rHLL}\), respectively) as ratios of
SVL. We then applied a second-order polynomial regression to estimate
selection gradients (Equation 5), and a P-spline function (Equation 6)
to visualize the shape of natural selection (Lande and Arnold 1983;
Blows and Brooks 2003; Gimenez et al. 2006; Ito and Konuma 2020):
\begin{equation}
\begin{matrix}f(\text{rTBD}_{i})=\ \mu+\ \beta\cdot\text{rHLL}_{i}+\frac{1}{2}\cdot\gamma\cdot\text{rHLL}_{i}^{2}\ \#\left(5\right)\\
\end{matrix}\nonumber \\
\end{equation}\begin{equation}
\begin{matrix}f(\text{rTBD}_{i})=\ \mu+\beta_{1}\cdot\text{rHLL}_{i}\ +\cdots+\beta_{p}\cdot\text{rHLL}_{i}^{P}\ +\sum_{k=1}^{K}{b_{k}\left(\text{rHLL}_{i}\ -\kappa_{k}\right)_{+}^{P}}\#\left(6\right)\\
\end{matrix}\nonumber \\
\end{equation}where \(\mu\) is the intercept, \(\beta\) is the linear selection
gradient representing directional selection, \(\gamma\) is the
non-linear selection gradient representing disruptive selection when\(\gamma>0\) or stabilizing selection when \(\gamma<0\), \(P\) is
the degree of freedom in the P-spline (which we set to 3), \(\eta\) is a
vector of parameters\(\left(\beta_{1},\cdots\beta_{P},\ b_{1},\ \cdots b_{K}\right)\),\(\left(\text{rHLL}_{i}-\kappa_{k}\right)_{+}^{P}\) is either\({(\text{rHLL}_{i}-\kappa_{k})}^{P}\) when
(\(\text{rHLL}_{i}-\kappa_{k})\ \geq\) 0 or 0 otherwise, and\(\kappa_{k}\) is k ’s fixed knots, with\(\kappa_{1}<\ \kappa_{2}<\cdots<\kappa_{k}\). Following Ruppert
(2002) and Ito and Konuma (2020), we set the number of knots to 35
where\(\ K=min(\frac{I}{4},\ 35)\), and used \(k/(K+1)\) quantiles
for all values of relative hind leg length, with \(k\) between 1 and 35.
We assumed that the coefficient \(b\) for\(\left(\text{rHLL}_{i}-\kappa_{k}\right)_{+}^{P}\) has a normal
distribution with mean 0 and variance \(\sigma_{b}^{2}\) (Gimenez et al.
2009). Since relative tail break distance was a positive continuous
variable, we applied the exponential function as the link function and
assumed it followed a gamma distribution.
To estimate each parameter, we used Markov Chain Monte Carlo (MCMC)
simulations in a Bayesian framework. We used non-informative
distributions as prior distributions for all parameters, the uniform
distribution for the interval [0, 30] for the prior distribution of
the hyper-parameters\(\left(\sigma_{1}^{2},\ \sigma_{2}^{2},\ldots,\ \sigma_{4}^{2}\right)\),
and used a normal distribution with a mean of 0 and variance of
1002 as the prior distributions for all
gradients\(\ {(\beta}_{1,k},\beta_{2,k},\ldots,\beta_{7,k})\), fixed
effects \((ep_{1,k},\ ep_{2,k},\ldots,ep_{4,k})\), intercept \(\mu\),
linear selection gradient β, non-linear selection gradient\(\gamma\) and the parameter \(\varphi\) in the gamma distribution. We
used an inverse-gamma distribution with the two parameters both at 0.001
for the prior distribution of\(\ \sigma_{b}^{2}\). We obtained the
posterior distributions for each parameter by generating three Markov
chains using 2,000 MCMC simulations sampled at a rate of five times
following an initial burn-in of 1,000 iterations. We standardized all
measurements using z-scores in order to avoid numerical instabilities
and to improve mixing each chain (Gilks & Roberts 1996). The posterior
distributions for all the parameters were summarized by 95% Bayesian
confidence intervals. We used Gelman and Rubin statistics to evaluate
the MCMC convergence (Gelman & Shirley 2011) and confirmed all
parameters underwent acceptable convergence. We performed MCMC
simulations using PyStan (https://github.com/stan-dev/pystan; accessed
15 December 2019) in Python, and deposited the source code for the
Bayesian model in $$$$$$$$ (e.g. GITHUB, ZENODO).
RESULTS
Foraging body temperatures , snake predator presence, and
recent climatic warming
Lizards on islands with the snake predator forage at warmer body
temperatures than those on islands without (Fig. 1), with respective
mean foraging body temperatures of 35.4 ± 0.1°C (n = 223; 135 on
Kozu, 42 on Mikura, 14 on Niijima, 14 on Oshima, and 18 on Toshima) and
32.2 ± 0.2°C (n = 219; 85 on Hachijo-Kojima and 134 on Miyake).
After controlling for data collection period and measured air
temperature in the shade with the linear mixed-effects model, we found
lizard foraging body temperatures to be 2.9°C hotter (p< 0. 01) on average on snake predator islands than on
snake-free islands (Table S3).
We observed a warming trend across islands, with annual average daily
means having increased from 16.9 ± 0.2°C in 1981-1991 to 18.3 ± 0.2°C in
2018-2019 (Fig. 2). We also found mean foraging body temperatures to
have increased from 35.2 ± 0.2°C in 1981-1991 (n = 183; 111 on
Kozu, 42 on Mikura, 12 on Oshima and 18 on Toshima) to 36.4 ± 0.2°C in
2018-2019 (n = 40; 24 on Kozu, 14 on Niijima and 2 on Oshima) on
islands where the snake predator is present, and from 31.7 ± 0.2°C in
1981-1991 (n = 134, all on Miyake) to 32.6 ± 0.3°C in 2018-2019
(n = 73, all on Hachijo-Kojima) on islands where it is absent
(Fig. 2). We found overall lizard foraging body temperatures in
2018-2019 to have, on average, increased by 1.3°C (p < 0.001)
relative to 1981-1991 (Table S3).