Modelling of natural selection of shell colour in mainland
As in a previous study (Ito & Konuma, 2020), we used the survival rate
as an index of fitness. First, second-order polynomial regression was
used to estimate the selection gradients with Equation 1 (Lande &
Arnold, 1983; Blows & Brooks, 2003; Gimenez et al. , 2009b). To
visualise the survival rate according to shell colour, penalized spline
regression (P-spline) was used with Equation 2 (Gimenez et al. ,
2006a; b).
\begin{equation}
\begin{matrix}\text{logit}\left(S_{i,t}\right)=\ \mu+\ \beta\cdot col_{i}+\frac{1}{2}\cdot\gamma\cdot col_{i}^{2}+\ \epsilon_{i}+\ B_{t}\ \#\left(1\right)\\
\end{matrix}\ \ \nonumber \\
\end{equation}\begin{equation}
\begin{matrix}\text{logit}\left(S_{i,t}\right)=\ f\left(\text{co}l_{i}\right)+\epsilon_{i}\ +\ B_{t}\ \#\left(2\right)\\
\end{matrix}\ \ \nonumber \\
\end{equation}\begin{equation}
f\left(\text{co}l_{i}\right)=\mu+\beta_{1}\cdot col_{i}\ +\cdots+\beta_{p}\cdot col_{i}^{P}\ +\sum_{k=1}^{K}{b_{k}\left(\text{co}l_{i}\ -\kappa_{k}\right)_{+}^{P}}\text{\ \ }\nonumber \\
\end{equation}\begin{equation}
t=1,2,\ldots,12;\ i=1,2,\ldots,\ 1714\ \left(\text{in\ the\ Izu\ Peninsula}\right)\nonumber \\
\end{equation}In these equations, \(logit(x)\) means \(log[x/(1-x)]\).\(\mu\) is the overall average for the survival rate on the logit scale.\(\beta\) is the linear selection gradient, which represents the
strength and direction of directional selection. \(\gamma\) is the
nonlinear selection gradient, which represents the strength of the
disruptive selection when \(\gamma>0\) or that of the stabilizing
selection when \(\gamma<0\). \(\epsilon_{i}\) represents the
individual heterogeneity, which follows a normal distribution with mean
0 and variance \(\sigma_{\varepsilon}^{2}\) (Royle, 2008). \(B_{t}\) is
the monthly effects on the logit scale for month \(t\). \(P\) is
the degree of freedom in the P-spline, which was set to \(P=3\) in
this study. \(\left(\text{co}l_{i}-\kappa_{k}\right)_{+}^{P}\) is
either \({(col_{i}-\kappa_{k})}^{P}\) when
(\(\text{co}l_{i}-\kappa_{k})\ \geq 0\) or 0 otherwise.\(\kappa_{k}\) represents k ’s fixed knots, with\(\kappa_{1}<\ \kappa_{2}<\cdots<\kappa_{k}\). The numbers of
knots \(K\) is decided according to \(K=min(\frac{I}{4},\ 35)\), so we
set to \(K\ =\ 35\) (Ruppert, 2002). We calculated \(k/(K+1)\)quantiles by using all values of shell colour, with \(k\) varying
between 1 and 35 (Ruppert, 2002). We assumed that the coefficient\(b_{k}\) of \(\left(\text{co}l_{i}-\kappa_{k}\right)_{+}^{P}\)follows a normal distribution with mean 0 and variance\(\sigma_{b}^{2}\) (Gimenez et al. , 2009a). The survival rate of
this species shows monthly fluctuations (Ito & Konuma, 2020), so we
considered the monthly effects \(t\) (\(t\ =\ 1,\ 2,\ \ldots,\ 12\)).
We constructed a multistate state-space model to estimate the parameters
in Equation 1 and 2. Natural selection was estimated while excluding the
predation effects by using two stages for the “dead” state: “dead
from predation” and “dead from other causes” (Schaub & Pradel, 2004;
Marescot et al. , 2015). This model was constructed by extending a
previous multistate model (Ito & Konuma, 2020). In this model, we
defined the state matrix with 12 states (Equation. S1): (1) survival
within the study site as a juvenile, (2) dead from a factor other than
predation, such as temperature, within the study site as a juvenile, (3)
dead from predation within the study site as a juvenile, (4) survival
within the study site as an adult, (5) dead from a factor other than
predation within the study site as an adult, (6) dead from predation
within the study site as an adult, (7) survival outside the study site
as a juvenile, (8) dead from a factor other than predation outside the
study site as a juvenile, (9) dead from predation outside the study site
as a juvenile, (10) survival outside the study site as an adult, (11)
dead from a factor other than predation outside the study site as an
adult, and (12) dead from predation outside the study site as an adult.
In the matrix, we defined the survival rate of juveniles and adults, the
predation rate of juveniles and adults, the site fidelity rate, and the
transition rate from juveniles to adults. The survival rate was applied
to Equation 1 or 2 separately according to our purpose. The predation
rate was assumed to be constant regardless of shell colour and month
from comparing the causes of mortality among marked snails (see the
results). The transition rate was assumed to have monthly fluctuation,
and the site fidelity rate was considered constant according to a
previous study (Ito & Konuma, 2020). The transition rate was modelled
with a random effects model.
In the observation matrix, we defined seven observations states
(Equation. S2): (1) first capture or recapture of a living juvenile
snail, (2) recovery of a juvenile snail that died from causes other than
predation, (3) recovery of a juvenile snail that died from predation,
(4) first capture or recapture of a living adult snail, (5) recovery of
an adult snail that died from causes other than predation, (6) recovery
of an adult snail that died from predation, and (7) an undiscovered
snail or unrecovered marked snail. In this matrix, the recapture rate
was for live individuals, and the recovery rate was for marked dead
individuals. These parameters were defined as differing between survey
months, although they were the same for juveniles and adults. In
addition, these parameters could vary according to shell colour (Ito &
Konuma, 2020). We modelled the effect of shell colour using a GLMM
framework, in which the link function was a logit function. In addition,
the same individual heterogeneity \(\varepsilon_{rp,i}\) was defined in
these parameters. This heterogeneity parameter is normally distributed
with mean 0 and variance \(\sigma_{\text{rp}}^{2}\). The probabilities
model obtained from the state and observation matrices were defined with
categorical distributions (Equations S3-S4; Kery & Schaub, 2012).