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