Evolutionary history of inshore oceanic island land snails diversified in shell colour

Oceanic islands provide an excellent opportunity to study the mode and tempo of phenotypic evolution of terrestrial organisms. Many studies have focused on oceanic islands far from the mainland. Oceanic islands near the mainland may provide distinct perspectives on phenotypic evolution, but a comprehensive understanding is still lacking. To address this gap, this study aimed to reveal when a land snail species inhabiting a volcanic archipelago within 40 km of the mainland diverged and how their shell colours evolved.


| INTRODUC TI ON
Islands provide a natural laboratory for ecological and evolutionary studies (Warren et al., 2015;Whittaker et al., 2017). In particular, oceanic island systems are best suited for studies on how traits have evolved for adaptation (Losos & Ricklefs, 2009). Since oceanic islands have never been connected to the mainland, terrestrial organisms are limited to those migrated by overseas colonisation.
Their abundance is markedly different between the oceanic islands, depending on the island's geographical distance from the mainland and the dispersal ability of the organisms. Consequently, the species composition could differ on each island (Aguilée et al., 2021).
Occasionally, surprising phenotypes can be found because island species are freed from the mainland interspecific interactions. For example, the morphology and colour of the shells of Mandarina snails on the Ogasawara Islands in Japan have been diversified because fewer predators are present (Chiba, 2002). The phenotype has also been diversified because of character displacement from closely related species (Chiba, 1999;Chiba & Davison, 2007;Davison & Chiba, 2006). Due to these reasons, many studies on phenotypic evolution have been carried out on oceanic islands far from the mainland, such as Galápagos (Grant & Grant, 2011;Kraemer et al., 2019) and Hawaii (Blonder et al., 2016;Zhang et al., 2021).
However, a unique phenotypic evolution has been observed on oceanic islands near the mainland (Melo et al., 2011). The colonisation rate should be higher in the islands close to the mainland (MacArthur & Wilson, 1967). Many predators and competitors can migrate to the coastal islands (Whittaker et al., 2008). The direction of trait evolution on these islands should differ from the oceanic islands far from the mainland (Lapiedra et al., 2018;Pringle et al., 2019;Reynolds et al., 2019). Therefore, inshore oceanic islands can provide new insights into the adaptive evolution of phenotypes.
The Izu Peninsula and Izu Islands in Japan provide an excellent system for studying the phenotypic evolution of terrestrial organisms on oceanic islands close to the mainland ( Figure 1). The Izu Peninsula was formed by submarine volcanic activity that occurred approximately 2.0 million years ago (Ma). The proto-Izu Peninsula was formed in parts by 1.35 Ma by uplifting of the ocean floor above the sea surface (Kano & Ito, 2016;Tani et al., 2011). The proto-Izu Peninsula was an island during its formation, but it is now a part of the Japanese mainland because of the Izu-Honshu collision that occurred due to the northward movement of the Philippine Sea Plate at approximately 1.0 Ma (Huchon & Kitazato, 1984). After this collision, the current Izu Peninsula was formed by many volcanic activities until 0.002 Ma (Hayakawa & Koyama, 1992;Koyama et al., 1995).
By contrast, the Izu Islands are located on the south-eastern coast of mainland Japan. These islands were formed because of volcanic activity that occurred less than 1.0 Ma (Tani et al., 2011). Oshima, the Izu Island closest to the mainland, is only 25 km from the mainland, whereas Aogashima, a remote Izu Island, is approximately 250 km from the mainland. There are approximately 10 islands between these two islands, which vary in size and elevation. Consequently, the species composition of reptiles, for instance, differs among the Izu Islands (Hasegawa, 2003;Hasegawa & Moriguchi, 1989).
Because of these compositional differences, predators and prey have evolved uniquely on these islands (Hasegawa, 1994(Hasegawa, , 2003. For example, phenotypic traits such as the body colour and hind leg length in Okada's five-lined skink Plestiodon latiscutatus can be influenced by different selective pressures from island-specific predators (Brandley et al., 2014;Landry Yuan et al., 2021). The feeding behaviour of the large Japanese field mouse, Apodemus speciosus, has shifted because of island-specific vegetation after migration to islands (Takechi & Hayashi, 2012). These geological and biological backgrounds enabled us to examine the inshore oceanic islandspecific evolutionary history of the phenotype on the Izu Islands.
Euhadra peliomphala simodae is one of the most diversified land snails in Japan in terms of shell colour (Hayashi & Chiba, 2004;Ito et al., 2021;Ito & Konuma, 2020). This snail inhabits the southern part of the Izu Peninsula and five Izu Islands (Figure 1; Oshima, Toshima, Niijima, Shikine and Kozu). Mainland snails tend to have bright shells, whereas peripheral island snails show various shell colours ranging from bright to dark (Hayashi & Chiba, 2004;Ito et al., 2021;Ito & Konuma, 2020). In addition, colour diversity differs among islands (Hayashi & Chiba, 2004). This may be explained by mark-recapture studies, which indicated that environmental changes after migration to the islands resulted in a difference in shell colour diversity Ito & Konuma, 2020). This study also demonstrated that predation pressure from the primary predator, A. speciosus, did not cause shell colour diversity, although predation pressure differed between the mainland and the islands . In other words, the island-specific environment, not predation pressure, could explain the shell colour diversity (Hayashi & Chiba, 2004;Ito et al., 2021). However, it remains unclear which island environmental factors determined snail shell colour diversity.
This study aimed to clarify the environmental factors of the islands that explain the diversity of shell colouration in E. p. simodae using a phylogenetic generalised linear mixed model (PGLMM).
We focused on four variables, three of which were environmental factors, and one was an index of historical demographic changes.
The first variable was the island area. Larger islands lead to a more heterogeneous environment, and it is more likely that small populations would be isolated within the islands (Gavrilets & Losos, 2009;MacArthur & Wilson, 1967). Indeed, correlations between oceanic island areas and the species diversity of snail communities have been reported (Kraemer et al., 2022;Parent & Crespi, 2006). The second variable was island elevation. Similar to island areas, elevation is also expected to increase environmental heterogeneity and contribute to species richness (Chiba et al., 2009;Parent & Crespi, 2006) and phenotypic diversity of the snail community (Kraemer et al., 2022).
The third variable is the distance from the surrounding landmasses, which affects the rate of immigration and extinction of the population (MacArthur & Wilson, 1967). The fourth variable is the rate of change in the effective population size because the founder effect also causes a bottleneck and can generate genetic drift and phenotypic differences (Carson & Templeton, 1984;Mayr, 1954). A previous study revealed the molecular phylogeny; however, the complex relationship between each island was unclear because the resolution of a molecular marker was low (Hayashi & Chiba, 2004). Therefore, we first estimated the fine-scale phylogenetic relationships and population structure of E. p. simodae using double-digest restriction site-associated DNA sequencing (ddRAD-seq). We then estimated the divergence time for each island population. Finally, using phylogeny with divergence times, we revealed the contribution of the above four factors to the diversification patterns of shell colour in land snails on oceanic islands.

| Sampling
A total of 117 snails were collected from Izu Peninsula (n = 42), Oshima (n = 5), Toshima (n = 16), Niijima (n = 24), Shikine (n = 14) and Kozu (n = 16) from 2016 to 2019 ( Figure 1). The sampling sites covered an area inhabited by E. p. simodae. Previous phylogenetic studies revealed that E. p. simodae is closely related to E. p. peliomphala (Hayashi & Chiba, 2000Shimizu & Ueshima, 2000). Therefore, E. p. peliomphala was used as an outgroup in this analysis. Euhadra eoa eoa was also used as an outgroup because it sympatrically inhabits the Izu Peninsula with E. p. simodae. The ddRAD-seq library was constructed using these samples (Table S1). The snail shell colour was quantified using a digital image after dividing each snail into the shells and soft parts (Appendix 1).

| ddRAD library construction and de novo assembly
Genomic DNA (gDNA) was extracted using the Qiagen DNeasy Blood & Tissue Kit (Qiagen), following a custom protocol. RNase A (Qiagen) was added as part of the extraction protocol. We then constructed two ddRAD libraries (Peterson et al., 2012) and sequenced them for 100 bp single-end reads using an Illumina Hiseq 4000 (Illumina). The raw sequenced Illumina reads were demultiplexed into individual reads. Then, we conducted quality filtering and de novo assembly using an 'ipyrad' on Python3 (Eaton, 2014;Eaton & Overcast, 2020). Because the number of consensus sequences was relatively low in three individuals from the Izu Peninsula, they were removed during single nucleotide polymorphisms (SNP) calling. The detailed methods for library construction and de novo assembly are in Appendix 2. Summaries of de novo assembly and SNP calling are presented in Tables S1 and S2 respectively.

| Phylogenetic and population structure analysis
Two approaches were used to estimate the phylogenetic relationships among E. p. simodae in each population. First, a maximum likelihood (ML) tree was constructed with 100 replicated bootstraps using RAxML HPC2 (Stamatakis, 2006). Analysis was conducted at the San Diego Supercomputer Center using the CIPRES Science Gateway (Miller et al., 2010). Second, a lineage tree was constructed based on multi-species coalescence with singular value decomposition (SVD) quartets using PAUP* software (Chifman & Kubatko, 2014;Swofford & Sullivan, 2009). All estimated trees were drawn using FigTree v1.4.4.
Population structure analysis was performed using STRUCTURE software (Pritchard et al., 2000) via the ipyrad-analysis toolkit after removing the outgroups. Additional STRUCTURE analysis was conducted using the Izu Island populations to obtain further insight into their population structures. A Markov chain Monte Carlo (MCMC) simulation was used to estimate each parameter of the STRUCTURE.
In total, 300,000 MCMC simulations were performed on five parallel chains, and initial 30,000 simulations were discarded as burn-ins.
Each analysis was performed with K = 1, 2, … , 15. In the second round of STRUCTURE, the analysis was conducted until K = 10 .
Principal component analysis (PCA) was conducted using the 'adegenet' library in R (Jombart, 2008;Jombart & Ahmed, 2011;R Core Team, 2022). The scaled allele frequency was calculated using the same SNPs dataset used in the first round of STRUCTURE. The mean allele frequencies replaced the gap sites in each SNP site because PCA could not be calculated when gap sites were present.
However, this operation may result in artefacts when there are many missing sites. Therefore, we used an SNPs dataset with a relatively low percentage of missing sites (14.2%) to avoid underestimation as far as possible.
The co-ancestry relationship in each snail was calculated using the fineRADstructure (Malinsky et al., 2018). Default settings were used to calculate the co-ancestry matrix, assign each snail to estimated haplotypes and construct a co-ancestry tree. The estimated co-ancestry tree was drawn following a modified R script using R (R Core Team, 2022).

| Population dynamics estimation
Approximation Bayesian computation (ABC) was used to estimate the effective population size, divergence time and mutation rate (Marjoram & Tavaré, 2006). The ABCtoolbox was used for the ABC analysis (Wegmann et al., 2010). The fastsimcoal2 software (Excoffier et al., 2013;Excoffier & Foll, 2011) was used to simulate the SNP alignments. Arlsumstat (Excoffier & Lischer, 2010) was used to estimate four statistical genetic summaries (Table S3); Pairwise F ST (Excoffier et al., 1992), Tajima's D (Tajima, 1989), Fu's F S (Fu, 1997) and pairwise differences (Tajima, 1983(Tajima, , 1993. To construct the basic topology of the coalescent simulation, a species phylogenetic tree was estimated using SVD quartets. The population of each island was treated as a species unit. The northern and southern parts of Niijima were treated separately because the genetic compositions of northern and southern Niijima were well differentiated (see Results). However, the genetic composition of northern Niijima was slightly nested in the southern population (see Section 3). The presence of gene flow between them was evaluated by comparing 10 scenarios using ABC before estimating the divergence time. Generation time, effective population size of each branch and minor allele frequency were assumed to be the parameters in the ABC analysis. Ten thousand simulations were conducted for each scenario. We then selected plausible scenarios from the top 10% of the simulations based on the rejection algorithm and the Bayes factor (Kass & Raftery, 1995).
Subsequently, divergence time was estimated using ABC. The divergence time, effective population size of each branch and the mutation rate were simulated from the prior distribution. Each simulation calculated two additional parameters: the changing rate in the effective population size after branching and the difference in the current effective population size. The rate of change was calculated as the ratio of the effective population size after branching to that before branching (Nadachowska-Brzyska et al., 2013). Differences in effective population sizes were calculated by subtracting the current effective population sizes from each other. Five base pairs were generated at 50,000 loci based on randomly selected values from the prior distribution. A total of 500,000 coalescent simulations were performed, and statistical summaries were calculated for each simulation. Finally, the top 1000 parameters were obtained based on the rejection algorithm. The divergence time can be converted into the number of generations used for the coalescent simulations by assuming that the generation time is 2.5 years (Appendix 3). Except for one parameter, the divergence time between Shikine and southern Niijima, a non-informative prior distribution was used (Tables S4 and   S5). This was assumed to have diverged 1134 years ago (Appendix 4).
Several pilot tests were run before the final estimation to determine prior distributions. The 'abc' package was used (Csilléry et al., 2012) in R for the model selection and parameter estimation.

| Analysis of shell colour diversity
We added measurements of snail shell colour from the Izu Peninsula was used on Python3 for GMMs (Pedregosa et al., 2011), and R was used for the binomial test (R Core Team, 2022).
A phylogenetic generalised linear mixed model (PGLMM) was used to examine whether the four independent variables contributed to shell colour diversity (de Villemereuil & Nakagawa, 2014).
Each population's variance in shell colour was used as an index of colour diversity. This was used as the response variable, and Gaussian regression was applied. For this analysis, we calculated and used the variance of every population on each island because the phylogeny and population structure were not divided within an island (see Section 3). Then, the variance within each island was modelled as a random effect, which followed a Gaussian distribution with mean 0 and variance 2 . As independent variables, this study used island area, elevation, the average distance from the surrounding landmasses and the changing rate of effective population sizes (Appendix 5; Table S6). Although previous studies have demonstrated that island age affects the species richness and functional diversity of land snails (Kraemer et al., 2022;Parent & Crespi, 2006), we did not use island age because the Izu Islands are relatively young (Tani et al., 2011) and detailed geological data are unknown. The species phylogeny with the estimated divergence times was converted to a variance-covariance matrix using the 'ape' library on R (Paradis & Schliep, 2019;R Core Team, 2022) and added to our model. In this analysis, the Izu Peninsula was removed because it is part of the mainland, and it is difficult to determine the precise area of E. p.
simodae. MCMC simulations with non-informative prior distributions were used to estimate the parameters. We standardised the four independent variables using the z-score to avoid numerical instability and improve the mixing of each chain (Gilks & Roberts, 1996).
Then, three Markov chains using 5000 MCMC simulations that were thinned at a rate of five following the initial burn-in of 2000 iterations were run to obtain the posterior distribution of each parameter. We confirmed the Gelman and Rubin statistics to be lower than 1.1 and convergence (Gelman & Shirley, 2011). The 'brms' package was used in R for the MCMC simulations (Bürkner, 2017(Bürkner, , 2018(Bürkner, , 2021 R Core Team, 2022).

| Phylogeny and population structure
Phylogenies estimated from the ML and SVD quartets were almost consistent ( Figure 2 and Figure S1). The populations on each island, except the Izu Peninsula and Niijima, formed a monophyletic group.
The snail population of Niijima is divided into northern and southern regions. Dark and bright shell colours did not show reciprocally monophyletic relationships within any island clade. Population structures showed similar phylogenetic trends. The genetic groups comprised six categories (Figure 3 and Figure S2). The Izu Peninsula is comprised of four genetic groups. However, the Izu Islands consisted of two genetic groups, one of which includes Oshima, Toshima and Kozu, and the other includes Niijima and Shikine. Toshima and Oshima slightly include Niijima and Shikine's composition. Oshima also slightly includes groups of the Izu Peninsula. When the snail population on the Izu Peninsula was excluded from the STRUCTURE estimation, the genetic groups were comprised of five categories ( Figure 3 and Figure S3)

| Population dynamics and divergences
The species tree estimated from the SVD quartets was concordant with the ML tree and STRUCTURE results for clades with strong bootstrap supports. Toshima, Oshima and Kozu were clustered in the same clade. In contrast, Niijma and Shikine belonged to a different clade ( Figure S6). The divergent histories of the Niijma and Shikine populations were selected for the two scenarios with a similar Bayesian posterior probabilities of 23% (Figure 4 and Figure S7; Table S7). One scenario was the no-gene flow scenario (Model 1), and the other was the unidirectional gene flow scenario from northern Niijima to southern Niijima after the divergence between southern Niijima and Shikine (Model 4). However, the Bayes factor of these two scenarios compared to other scenarios assuming recent gene flow (Models 2 and 3) was lower than 2.0, indicating that the evidence of differences between these two scenarios was weak (Table S8).
The divergence times for all populations are shown in Figure 5.  (Table S9). When comparing the current effective population sizes, the 95% BCI of the difference included zero; therefore, no significant differences were observed (Table S10).

| Shell colour diversity against three environmental factors and a historical event
The frequency distribution of dark and bright colours differed among populations ( Figure. 5 and Tables S12). Most populations showed skewed shell colour frequencies in dark or bright colours ( Figure 5 and Table S13). More than 70% of the snails on the Izu Peninsula and Kozu had bright shell colours, whereas Oshima, southern Niijima and Shikine were biased towards dark snails at 69%, 75% and 96% respectively. In contrast, Toshima and northern Niijima had similar colour frequencies. The PGLMM showed that the posterior probability of the coefficient of elevation was 93.4% positive for shell colour diversity, although the 95% BCI included zero ( Figure 6 and Table S14).
The posterior probability of coefficients on island area, the average distance from surrounding landmasses and changing rate of effective population size were relatively less positive (65.8%, 39.7% and 33.0% respectively), and their 95% BCI included zero.

| Evolutionary history of Euhadra peliomphala simodae
We found that the divergence into each island occurred between 0.21 and 0.74 Ma (95% BCI: 0.02-1.43; nodes 2-4 in Figure 5). The current Izu Islands were formed by Quaternary volcanic activity since  (Hayashi & Chiba, 2000;Shimizu & Ueshima, 2000), but they may also tend to diversify more between mainland and island populations than other animals. These biological characteristics may explain why the divergence time of land snails was older than that of other animals on the Izu Islands.
Population structure analysis revealed that land snails are genetically grouped on each island. In addition, the probability of gene flow between ancient northern Niijima and Shikine is extremely low. These results suggest that the gene flow is relatively limited between inter-island populations. However, land snails possess a passive dispersal ability across the ocean, which can explain their colonisation of oceanic islands. For example, some snails exhibit salt tolerance by sealing themselves with an epiphragm at the shell aperture to protect their soft body. They can survive for over 20 days, even if they are swept out to the sea, for example, by a storm (Ożgo et al., 2015). Some land snails on the Izu Islands may have been dispersed to approximately 1000 km by ocean currents (Hirano et al., 2019). Furthermore, other factors, such as wind and bird migration, are related to passive dispersal in oceanic systems (Cowie & Holland, 2006). For example, micro-land snails on an oceanic islands can be transported by surviving in the digestive tracts of birds (Wada et al., 2012). Although adult snails of E. p.
simodae (approximately 30 mm) are larger than micro snails, dispersion via bird predation may be conceivable because juveniles and eggs are of the same size as micro snails (less than 10 mm), and the Izu Islands are close to the mainland. Therefore, it may be easy to limit the gene flow between island populations due to immobility, but they are physiologically tolerant to long passive transportation.
Thus, they can form endemism on the remote islands in the early stages.
A similar result, in which the divergence time of land snails was older, was reported in a study on an inshore oceanic island of the Canary Islands (Greve et al., 2010). The Canary Islands are 100-500 km away from the mainland, which is a relatively small distance.  In contrast, the divergence time of land snails is in line with that of other animals on the Galápagos Islands (Parent et al., 2008;Phillips et al., 2020) and Hawaii Islands (Cowie & Holland, 2008;Holland & Hadfield, 2004). These animals diversified during the island formation. The distances to the mainland are 950-1200 km for the Galápagos Islands and 3000 km for the Hawaiian Islands. Therefore, the distance from the source landmass may explain the difference in divergence time between immobile animals and other animals, although the present study alone does not clarify this.
Our results also showed that multiple genetic groups were distributed within a single island, but only in Niijima. The genetic structure of southern Niijima was almost the same as that of Shikine, an island approximately 2.5 km from Niijima. In contrast, the genetic structure of northern Niijima was half-mixed with the northern Niijima and southern Niijima (Shikine) components. This can be attributed to the secondary contact. The estimated historical demography supported scenarios with a high probability, in which, after the formation of south Niijima at 886 AD (Sugihara et al., 2001;Tsukui et al., 2006), the population in Shikine migrated into southern Niijima, and the gene flow between northern and southern Niijima still persists. Nevertheless, our results do not completely dismiss a scenario that does not involve gene flow. In this case, with no gene flow, some loci might have been derived from incomplete lineage sorting (ILS) because the genetic structure of northern Niijima was nested with that of southern Niijima (Shikine).
The divergence between the two major insular clades occurred at 1.35 Ma (95% BCI: 0.58-2.51; node 5 in Figure 5). The confidence interval includes the times of Izu Island formation, the proto-Izu Peninsula rising above sea level and the formation of the proto-Izu Peninsula (Kano & Ito, 2016;Tani et al., 2011). Thus, the cause of divergence between the two major insular clades remains unclear. In contrast, the divergence between the peninsular and insular assemblages occurred at 2.11 Ma (95% BCI: 0.92-2.96; node 6 in Figure 5).

This timing was considered a period when both the proto Izu
Peninsula and the Izu Islands had not yet formed (Kano & Ito, 2016;Tani et al., 2011). Therefore, divergence may have occurred in Cynops pyrrhogaster, in the current southern and central parts of the Izu Peninsula diversified at 1.9-5.1 Ma (Tominaga et al., 2015). The Izu-Honshu collision of the Tanzawa Block via the northward movement of the Philippine Sea Plate occurred at 2.5-6.8 Ma (Ishikawa et al., 2016). Hence, these geological events might have led to the above divergences. However, it is poorly understood whether landmasses such as the proto or the current Izu Peninsula and Izu Islands were lifted above the sea level. Thus, it must be determined when and where these landmasses exist from the paleontological and geological approaches. This clarifies the cause of diversification in the Izu Peninsula and Izu Islands in the early period.

| Shell colour diversity of Euhadra peliomphala simodae
The posterior probability of the coefficient of elevation was 93.4% positive in the PGLMM results. This suggests that higher elevation islands tend to have higher colour diversity. This result supports previous studies demonstrating that island elevation contributes to the species diversity and functional morphological diversity of land snails (Chiba et al., 2009;Kraemer et al., 2022;Parent & Crespi, 2006). An ecological study using E. p. simodae suggested that environmental factors cause natural selection that works on shell colour . Thus, an environmental factor that covaries with elevation is a candidate factor that affects shell colour. Temperature, precipitation, frogs and ultraviolet radiation are environmental factors that covary with elevation. For example, in respect to temperature, some studies using land snails have suggested that shell colour plays a role in thermoregulation (Heath, 1975;Knigge et al., 2017). Bright colours may help maintain optimum body temperature when exposed to high temperatures. In contrast, dark colours could have the advantage of quickly warming up at low temperatures. Furthermore, it is important to consider biological factors that may be correlated with elevation, such as vegetation and intensity of interspecific interactions. For example, a previous study found that vegetation varied with elevation on some Izu Islands (Kojima & Fujiwara, 2008). In addition, Parent and Crespi (2006) showed that elevation has an in- colour diversity due to bottleneck (Carson & Templeton, 1984;Mayr, 1954)  tion with a small number of individuals (Carson & Templeton, 1984;Mayr, 1954). A previous study using a partial mitochondrial DNA sequence also showed a similar result: wherein the genetic diversity in the Izu Peninsula was higher than that in the Izu Islands (Hayashi & Chiba, 2004). This is contrary to the population dynamics estimation. Therefore, we could not rule out the possibility of genetic drift. Regardless of whether a bottleneck occurred, our PGLMM results suggest that the island area, the average distance from the surrounding landmasses and changing rate of the effective population size did not explain shell colour diversity.
The frequency of shell colour in snails is highly diverse on the Izu Peninsula and in northern Niijima, Oshima and Toshima; however, the shell colour profile differs depending on the island. Many snails on the Izu Peninsula exhibit a bright shell colour through stabilising selection during adulthood . However, stabilising selection does not work among juvenile snails , implying that only a small number of dark-shelled snails persisted on the Izu Peninsula.
Thus, the population on the Izu Peninsula would have persisted with large colour diversity. In contrast, the frequency of dark-shelled snails was higher on the three islands of northern Niijima, Oshima and Toshima than on the Izu Peninsula. When a population expands into a novel environment, a phenotype with a low frequency that has previously persisted may sometimes increase due to selection (Colosimo et al., 2005). This implies that dark-shelled snails increased on the three islands of northern Niijima, Oshima and Toshima, while maintaining a large colour diversity. In addition, both northern Niijima and Toshima showed a symmetric dimorphism in snail shell colour.
In other words, the frequency of dark snails increased by as much as that of bright snails. In northern Niijima, disruptive selection works on snail shell colour (Ito & Konuma, 2020); therefore, darker snails may have the advantage of surviving as brighter snails in Toshima as well as northern Niijma. These two island populations were genetically positioned in different clades. Hence, the parallel evolution of symmetric colour dimorphism may have occurred because of natural selection in E. p. simodae on the two islands.
The shell colour profile differed on three islands (populations), southern Niijima, Shikine and Kozu, where shell colour did not diversify. Southern Niijima and Shikine had dark-shelled snails, whereas Kozu had bright-shelled snails. One possible explanation for the skewed colour distribution is a shift in natural selection after the migration (Hunt & Rabosky, 2014). For example, the adaptive peak in the shell morphology of Mandarina land snails found on the Ogasawara Islands, Japan, could rapidly shift to another peak owning to environmental changes (Chiba, 1996). However, the mechanisms underlying the natural selection working shell colour on these three islands are currently unknown. It has been observed that natural selection acts on shell colour in the Izu Peninsula and Niijima of the Izu Islands Ito & Konuma, 2020), so it is plausible that the same may occur for these islands. A possible alternative explanation is that a phenotypic shift is not dependent on natural selection because a mathematical model and fossil record analysis revealed that genetic drift could cause a dramatic shift in phenotypic frequency (Hunt & Rabosky, 2014;Lande, 1976). This study could not determine whether a bottleneck occurred during the migration into the islands, as previously discussed. Thus, our study was limited to suggesting the cause of shell colour shift up to the current study. Unravelling the function of shell colour will lead us to elucidate the functional mechanism of the phenotypic shift in shell colours on the island.

| CON CLUS ION
The genetic composition of Euhadra peliomphala simodae differed between the islands. The only exception was Niijima, where snails from the northern and southern parts of the island exhibited distinct genetic compositions. We also observed that E. p. simodae began to diversify approximately 2.0 Ma. Until about 0.2 Ma, the formation of each island population was complete. These divergence times are greater than those of other animals on the Izu Islands. It was suggested that both proximity to the mainland and the immobility of land snails were the causes of early divergence. Moreover, island elevation positively contributed to shell colour diversity. In contrast, island area, distance from surrounding landmasses and the changing rate in effective population size did not correlate with colour diversity. In addition, populations with similar shell colours were not always phylogenetically close to each other, suggesting parallel evolution. Therefore, we conclude that environmental factors that covary with island elevation would specify the colour profiles of different snails on each island via natural selection.

ACK N O WLE D G E M ENTS
We thank Daishi Yamazaki for providing a sample of the outgroup.
We also thank Osamu Kagawa, Takahiro

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Demultiplexed fastq files on ddRAD-seq are available from the NCBI (DRA015787; Junji Konuma is an associate professor at Toho University, Chiba, Japan. He studies the mechanisms of ecological specialisation using the field research, morphometrics and genomics of the Damaster beetle.
Authors' contributions: All authors conceived the ideas; SI designed the methodology, collected and analysed the data and led the manuscript's writing. All authors contributed critically to the drafts and gave final approval for publication.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information can be found online in the Supporting Information section at the end of this article.