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).