Parallel evolution of passive and active defence in land snails

Predator-prey interactions are major processes promoting phenotypic evolution. However, it remains unclear how predation causes morphological and behavioural diversity in prey species and how it might lead to speciation. Here, we show that substantial divergence in the phenotypic traits of prey species has occurred among closely related land snails as a result of adaptation to predator attacks. This caused the divergence of defensive strategies into two alternatives: passive defence and active defence. Phenotypic traits of the subarctic Karaftohelix land snail have undergone radiation in northeast Asia, and distinctive morphotypes generally coexist in the same regions. In these land snails, we documented two alternative defence behaviours against predation by malacophagous beetles. Furthermore, the behaviours are potentially associated with differences in shell morphology. In addition, molecular phylogenetic analyses indicated that these alternative strategies against predation arose independently on the islands and on the continent suggesting that anti-predator adaptation is a major cause of phenotypic diversity in these snails. Finally, we suggest the potential speciation of Karaftohelix snails as a result of the divergence of defensive strategies into passive and active behaviours and the possibility of species radiation due to anti-predatory adaptations.

distribution pattern 37,38 . Two of these snail species-Karaftohelix (Ainohelix) editha and Karaftohelix (Ezohelix) gainesi-on Hokkaido Island illustrate the extreme inter-and intra-specific levels of phenotypic variation of this group 38 . K. editha and K. gainesi have very different shell morphologies; therefore, these species previously belonged to different genera. However, these two species are current nearly indistinguishable both genetically and anatomically 38 . Because frequent hybridization has occurred between these two species, it is likely that divergence of the phenotypic traits of K. editha and K. gainesi evolved relatively rapidly due to natural selection 38 .
However, the factors driving this natural selection are not clear because K. editha and K. gainesi often inhabit the same region without divergence in habitat use 38 . Therefore, differences in habitat may not be the main selective pressure 38 . On the other hand, it is plausible that the phenotypic differences between land snails in this region have evolved in part due to predation pressure because predators that specialize on land snails, including some types of carabid beetles, are distributed throughout northeastern Asia 39 . Some species of carabid beetles were well known as the specialist predators of land snails 40,41 . For example, Damaster blaptoides and Acoptolabrus gehinii are distributed throughout most of Hokkaido Island as are K. editha and K. gainesi. It is likely that these carabid beetles prey on K. editha and K. gainesi because there are a few other large snail species on Hokkaido, but no species other than K. editha and K. gainesi are distributed throughout all of Hokkaido 42 .
Although a number of land snail predators are known 40,41 , there are not many malacophagous species on Hokkaido and southern far-east Russia. In addition, all of the malacophagous species other than carabid beetles are rare in these regions. Mammals and birds can also be predators of these land snails, but there is no evidence of such predation except by a chipmunk, Tamias sibiricus, in these areas. However, T. siviricus is omnivorous, and they prey on relatively low numbers of snails 43 ; therefore, this species may not be an important snail predator.
Malacophagous flat worms and a burying beetle, Phosphuga atrata, were also distributed in these regions but also relatively rare. In contrast, carabid beetles, especially D. blaptoides and A. gehinii, are very abundant on Hokkaido, and Coptolabrus smaragdinus and Acoptolabrus schrencki on southern far-east Russia. Thus, it is likely that the phenotypic divergence and speciation of land snails in northeast Asia have been driven by predation pressure from carabid beetles.
To test this hypothesis, we examined the morphology of several land snail species in northeastern Asia, exposed them to natural predators, and observed their defensive responses. Second, we estimated the evolutionary relationships among the different species and morphotypes using molecular phylogenetic approaches. Finally, we discussed the distinctive defence behaviours that Karaftohelix snails display and how these behaviour and defence-related morphologies have evolved and resulted in species radiation.

Materials and Methods
Samples. Nine genetically related bradybaenid species were collected and analysed. Five species of  Table 1). All bradybaenid land snails used in this study are large species (larger than 8 mm). We found that several congeneric species coexisted in locations on Hokkaido Island and in far-east Russia. Acusta despecta (Pulmonata: Bradybaenidae) was used as an outgroup for phylogenetic analyses because it has been shown that the genus Acusta is the sister genus of Karaftohelix 45,46 . Behavioural observations. In total, 55 adult snails from eight species (four island species and four continental species; Supplementary Table 1) were used for behavioural observations. The foot of each snail was given an artificial external stimulus by pushing on it with fine-tip tweezers, and the response behaviour of the snail was observed. All snails used for behavioural observations were cultured individually in plastic cases (15.5 cm × 11.0 cm × 4.5 cm) with wet tissue paper at the bottom in room temperature (20~25 °C) for three to seven days before observations. Some trials were recorded with a video camera (HDR-XR500 V; SONY, Japan).
In addition, ten adult K. editha and K. gainesi snails from Bibai (locality no. 19-3 in Fig. 2 Table 1). In each experiment, one beetle and one snail were put into a plastic case (15.5 cm × 11.0 cm × 4.5 cm) with horticultural soil at the bottom for 15 minutes under low-intensity light, and the behaviour of the snail in avoiding predatory attacks by the beetles was observed. All beetles were fed sufficient amounts of fish meat sausage three days before trials. Some trials were recorded with a video camera (HDR-XR500 V; SONY, Japan). Table 1). Nine shell morphological traits were measured from pictures of the shell (Fig. 3). Traits included: aperture height (AH), aperture width (AW), shell diameter (D), total shell height (H), shell height above the aperture (SH), spire width (SW), number of coils (NC), aperture area (AA) and total area in the shell (AT). The shell shape and size were analysed separately. A principal component analysis (PCA) was used for the analysis of shell shape; this was conducted with JMP software (SAS Institute, North Carolina) using mean values of the five lengths relative to the shell diameter (AH/D, AW/D, H/D, SH/D, SW/D and NC/D) from each species and from each locality. In addition, the mean relative aperture area (AA/AT) of each species from each locality was compared among the different behavioural traits. The shell diameter (D) of each species was also compared among the different behavioural traits for the size analysis. All measurements are shown in Supplementary Table 2 proteinase K, incubated at 60 °C for approximately 1 hour, extracted once with phenol/chloroform and precipitated with two volumes of ethanol. The DNA pellet was then rinsed with 70% ethanol, vacuum-dried for approximately 1 hour and dissolved in 50 μ L of distilled water.

Morphological analyses. Shell morphological analyses were conducted for 165 individuals of nine species (Supplementary
To estimate the phylogenetic relationships among the collected snails, sequenced fragments were sampled from two mitochondrial DNA regions (cytochrome oxidase subunit 1 (CO1) gene (~530 bp) and 16S ribosomal DNA (16S; ~700 bp)) and from two nuclear DNA regions (ribosomal internal transcribed spacer regions 1 and 2 (ITS; ~1200 bp) and external transcribed spacer region (ETS; ~380 bp)). Polymerase chain reaction (PCR) conditions and the primers used are shown in Supplementary Table 3. The PCR products were purified using Exo-SAP-IT (Amersham Biosciences, Little Chalfont, Buckinghamshire, UK). The sequencing cycle was carried out with both forward and reverse primers using ~80-100 ng of PCR product in the reaction and the BigDye TM Terminator v3.0 Cycle Sequencing Ready Reaction Kit (Applied Biosystems, California). The DNA sequences were electrophoresed on a 310 Genetic Analyser or a 3130 Genetic Analyser (both Applied Biosystems, California).
Phylogenetic analyses. Sequences were aligned using MUSCLE v3.8 47 , and the results were cleaned of problematic alignment blocks using GBLOCKS v0.91 48 with the default parameters. Gene trees were constructed using Bayesian inference (BI), maximum likelihood (ML), maximum parsimony (MP) and neighbour joining (NJ) models based on the combined dataset with all sequences (16S, CO1, ITS and ETS).
Prior to the ML and BI analyses, the appropriate models of sequence evolution were selected with Kakusan software version 4-4.0.2011.05.28 49,50 , and all the sequences were combined. BI analysis used MrBayes v3.1.2 51 . Tree spaces were explored using two concurrent runs with four simultaneous Markov chain Monte Carlo (MCMC) simulations for 10 million generations with sampling every 100 generations for the combined data set of all sequences (16S, CO1, ITS and ETS). The number of generations before stationarity of likelihood values was obtained and was estimated with TRACER v1.5 software 52 such that the effective sample sizes of all parameters were more than 190 after the burn-in. The heating parameters were set to 0.15. After discarding the first 10001 trees as the burn-in, the 50% majority rule consensus tree and the posterior probabilities of nodes in the tree were obtained.
The ML analysis was performed with TREEFINDER v2008 53 under the maximum likelihood criterion. The MP and NJ trees were reconstructed using MEGA v6.0 54 . Prior to the MP and NJ analyses, the 16S, CO1, ITS, and ETS sequences were combined using MEGA. Nodal support for the ML, MP and NJ analyses were assessed using bootstrap analyses 55 with 1000 replications.

Results
Behavioural observations. The observed behaviours were classified into two main categories-passive defence and active defence. Almost all individuals from all six species (K. editha, K. blakeana, K. takahidei, K. maackii, K. middendorffi and K. ussuriensis) retracted their soft body into their shells very quickly, which is a passive defence behaviour (Table 1; Fig. 4a; Supplementary Movie 1, 3, 4, 5). In contrast, no individuals from the other two species (K. gainesi and K. selskii) retracted their soft body into the shell. Rather, they became even more active than before the external stimulus and vigorously swung their shells. This motion was usually repeated several times at the same frequency (approximately one swing every three seconds; Table 1; Fig. 4b; Supplementary Movie 2, 6). Two individuals of K. blakeana did not show a quick response and finally retracted their soft body into the shell, and one individual of K. gainesi and two individuals of K. selskii showed a different behaviour in which they created bubbles around their soft body (indicated as "other behaviours" in Table 1; Supplementary Movie 7). The behaviour of K. editha, K. blakeana and K. gainesi was also observed tentatively in the wild, and the same behaviour observed in the laboratory was seen in the wild.
These extremely different behaviours were clearly associated with snail species. All species were separated into two groups: "passive defence species" (K. editha, K. blakeana, K. takahidei, K. maackii, K. middendorffi and K. ussuriensis) and "active defence species" (K. gainesi and K. selskii). When the snails received external stimulus, the passive defence species retracted their soft body into the shell, and the active defence species showed aggressive behaviours.
These behaviours were also observed when the malacophagous carabid beetles attacked the snails (n = 10 for each of two snail species, K. editha and K. gainesi; Supplementary Table 1). The K. editha eventually escaped from the predator by retracting deep inside the shell (Supplementary Movie 8). In contrast, K. gainesi escaped from the predator by flipping away or even knocking the predator over with its shell (Supplementary Movie 9, 10).

Morphological analyses.
More than 91% of the variation in shell morphology among the individual snails was explained by two principal components (PC1 and PC2; Table 2; Fig. 5a). The AW/D, AH/D, SW/D and NC/D had large loading values (more than 0.8 or less than − 0.8). These four traits are related to the relative aperture size; thus, PC1 can be interpreted as a factor explaining the relative aperture size. Two other factors, H/D and SH/D, had high loading values on PC2 (more than 0.8). Therefore, PC2 can be interpreted as a factor explaining the relative shell height.
Morphologically, passive and active defence species were clearly separated from each other based on PC1 scores (Steel-Dwass test, P < 0.001). In contrast, PC2 scores were not significantly different between the behavioural groups (Steel-Dwass test, P > 0.05). These results clearly indicated that the relative aperture size was much larger in active defence species than in passive defence species. The PCA results were confirmed more directly by comparing the relative aperture area (AA/AT) between passive and active defence species (Fig. 5b; Steel-Dwass test, P < 0.001).
In addition, the shell diameter (D) was significantly larger in active defence species than in passive defence species (Fig. 5c; Steel-Dwass test, P < 0.001). This might indicate that the two different defensive behaviours were associated with different shell sizes, although the shell size of K. selskii was barely larger than the passive defence species.
Phylogenetic analyses. In the molecular phylogenetic analyses, 74 individuals of the nine species as well as the outgroup taxa were analysed, and 60 haplotypes were detected. All analyses (BI, ML, MP and NJ) resulted in nearly identical topologies. The inferred phylogenetic relationships among the haplotypes are shown in Fig. 6.
Two major clades were identified, and these clades corresponded to island and continental species. The monophyly of the island clade was strongly supported (Bayesian posterior probability (BPP) = 1.00, bootstrap support value (BV) for ML, MP and NJ analyses = 65, 88 and 94%, respectively). The monophyly of the continental clade was also well supported in all trees except for the ML tree (BPP = 0.90, BV for MP and NJ = 98 and 99%, respectively). In the ML analyses, the relationships among the continental species and populations were unclear.
The island clade was further subdivided into approximately four groups (subclades I-a, I-b, I-c and I-d; Fig. 6). Each subclade was represented by only a single species except for subclade I-a, which included K. editha and K. gainesi. Only one individual of K. editha from Yubari (Ke-22-1) was not included in subclade I-a.
In the continental clade, individuals collected from the different regions tended to have very distinctive genotypes even within the same species (e.g., K. maackii and K. middendorffi). However, obvious subclades were not  recognized, and the relationships among haplotypes were not clearly resolved due to incongruence of the topologies among the BI, ML, MP and NJ models.
The passive and active defence species were both separated into islands and continental clades indicating that the divergence of passive and active defence species occurred independently on the island and on the continent (Fig. 6).

Discussion
Two alternative anti-predator behaviours, passive and active defence, were documented among closely related Karaftohelix snails. Although further studies are needed, the results of feeding experiments suggest that these two alternative behaviours have the same function-avoiding predation by malacophagous carabid beetles. The passive defence snails use their shell as a "shield" to defend their soft body from the predator's attack, whereas the active defence snails use their shell as a "club" to hit the predators and knock them over. This study is one of only a few to report on land snails using their shell for active defence by swinging it against a predator. One example of a similar but more obscure behaviour in a Japanese bradybaenid land snail, Acusta despecta, has been described when these snails are attacked by the larvae of fireflies [56][57][58] .
The hypothesis that predation pressure led to morphological divergence in these snails seems to be reasonable because the alternative behaviours demonstrated by the different taxa are associated with differences in shell shape and size. The shell morphology analyses indicated that the relative aperture size of the shell was strongly associated with behavioural differences (Fig. 5a,b). In addition, the shell diameter might be larger in the active defence species than in the passive defence species when coexisting species are compared (Fig. 5c). For the passive defence snails, which retract as an anti-predatory behaviour, a shell with a narrower relative aperture prevents the predator from inserting its head in the shell 41,57,58 . The outer lip of the adult shell of all the passive defence species is also markedly thickened, suggesting that this characteristic is effective in protecting the shell from being broken by the predator when the soft body is deeply retracted into the shell 41,57,58 .
In contrast, active defence snails swing their shell as an anti-predatory behaviour, and a larger relative aperture might allow development of strong muscle to swing the shell around the soft body. In addition, a large relative aperture size relates to a relatively large body size-this can help shake off the predator and can even damage the predator. Aperture size is positively correlated with foot size 59 and thus with muscle mass. Thus, these two different defensive strategies are incompatible because there is a fundamental functional trade-off between those that use the shell as a "shield" versus a "club". The morphological analyses support this idea because there were no intermediate morphotypes between species with both passive and active defence strategies.
The phylogenetic analyses clearly suggested that the passive and active defence species and the morphotypes related to these defensive strategies arose independently on the islands and the continent, although the divergence pattern of the island clade was more complex than continental clades (Fig. 6). This may suggest that the divergence of passive and active defence strategies and island speciation has been ongoing. This pattern-parallel evolution of similar adaptive traits in several independent regions-strongly implicates natural selection against predation pressure as the cause of the evolution of these traits [4][5][6][7]12 . It is unlikely that the morphological differences among the several species are due to major differences in habitat because there were no obvious differences in the local microhabitats occupied by species when they coexisted. The hypothesis that predation pressure led to speciation is a uncommon explanation of morphological divergence because adaptation to different microhabitats is a major factor underlying phenotypic divergence among species 4,6,8,14,15,60 .
The divergence in body size between coexisting passive and active defence species may promote the evolution of reproductive isolation between these species-especially between K. editha and K. gainesi on Hokkaido Island. Although the shell diameter of K. selskii is only slightly larger than the sympatric K. middendorffi, the size of the soft body is larger in the former than the latter because of the larger relative aperture size of the former. In bradybaenid snails, differences in body size cause reproductive isolation because of the presence of size-assortative mating 61 . Although further experimental approaches are needed, this study implies that divergence in defence strategies against the same predator can cause speciation in Karaftohelix. The patterns of Karaftohelix diversification indicates that an ecological radiation occurred among these land snails because they clearly share a common ancestry (Fig. 6) 45,46 . Although the exact ages of speciation of each species are difficult to estimate, divergence of passive and active defence species appears to have started 1-3 Ma within the islands and continent based on the evolutionary rate of bradybaenid 16SrRNA 62 . This suggests that there has been enough time for the populations to diverge into highly distinctive phenotypes.
The active defence species K. gainesi coexists with one or two congeneric passive defence species on the islands, and another active defence species K. selskii coexists with two or three congeneric passive defence species on the continent. However, no clear relationship is found between morphology and number of coexisting congeneric species. A large overlap of habitat and/or resource use among these snail species suggests that interspecific competition among sympatric species is weak. Differences in shell size among the sympatric passive defence species are unlikely to be caused by interspecific competition because of no difference in habitat use among these species. Although further analyses are needed, we speculate that the larger passive defence species is more advantageous to protect the shell from being broken by the beetle and the smaller passive defence species is more advantageous to prevent the beetle from inserting its head in the shell.
Therefore, predation is shaping the evolutionary change among these land snails. Morphological changes in relative aperture size represent an ecological trade-off and only one strategy can be employed.
This type of radiation does not follow the existing ecological models of radiation. The evolutionary pattern of the bradybaenid land snails observed here seems to follow the model of prey species divergence by "apparent competition" 63 . Similar examples of phenotypic divergence of prey driven by a small number of predator species have been shown in some previous studies on freshwater and land snails 64,65 . Although further studies are needed to clarify the genetic patterns of speciation, morphological variations and behavioural traits in these snail, the present findings shed light on ecological factors other than resource competition that are important forces to drive phenotypic divergence and species radiation.