High-throughput SNPs dataset reveal restricted population connectivity of marine gastropod within the narrow distribution range of peripheral oceanic islands

Molecular studies based on the high resolution genetic markers help us to grasp the factor shaping the genetic structure of marine organisms. Ecological factors linking to life history traits have often explained the process of genetic structuring in open and connectable oceanic environments. Besides, population genetic divergence can be affected by fragmented habitat, oceanic current, and past geographical events. In the present study, we demonstrated the genetic differentiation of marine gastropod Monodonta sp. within a narrow range of peripheral oceanic islands, the Ogasawara Islands. Genetic analyses were performed not only with a mitochondrial DNA marker but also with a high-throughput SNPs dataset obtained by ddRAD-seq. The results of the mtDNA analyses did not show genetic divergence among populations, while the SNPs dataset detected population genetic differentiation. Population demographic analyses and gene flow estimation suggested that the genetic structure was formed by sea level fluctuation associated with the past climatic change and regulated by temporal oceanographic conditions. These findings provide important insights into population genetic patterns in open and connectable environments.

After larval settlement and metamorphosis, M. confusa matures to a shell height of approximately 10 mm, which takes 1 to 2 years. Around the Japanese archipelago, a previous molecular study on M. confusa detected a significant genetic population structure, but the level of differentiation was smaller than that of closely related species 46 .
For Monodonta sp. distributed in the Ogasawara Islands, previous phylogenetic analyses based on standard molecular markers showed no differentiation within the Ogasawara Islands 35 . However, given the relatively short planktonic larval duration of Monodonta, it is possible that Monodonta sp. is genetically differentiated across the archipelagos. To clarify the genetic structure of marine species within the Ogasawara Islands, it is necessary to carry out a high resolutive molecular approach.
Here, we aimed to determine the level of genetic differentiation and patterns of genetic structuring within a narrow range of peripheral oceanic islands using marine organisms Monodonta sp. as our model. In addition, we tested the hypothesis that genetic differentiation occurred due to historical climatic events and regulated temporal oceanographic conditions. Genetic analyses were performed using a mitochondrial DNA (mtDNA) marker, and a high-throughput SNPs dataset obtained by ddRAD-seq. Based on the dataset obtained by these two types of genetic markers, we discuss the genetic structure, temporal population connectivity, and timing of genetic divergence in open and connectable marine environments of peripheral oceanic islands.

Results
Genetic structure estimated by mtDNA datasets. The obtained length of the COI sequence alignment was 513 base pairs. Monodonta sp. had nine haplotypes from 59 individuals of four sampling islands. The haplotype diversity and nucleotide diversity of each population are shown in Table 1. The haplotype diversity ranged from 0.53 (locality 2: Yomejima Island) to 0.70 (locality 1: Kitanoshima Island), and the nucleotide diversity ranged from 0.00136 (locality 2) to 0.00195 (locality 1). The AMOVA analysis did not show the existence of population genetic differentiation in Monodonta sp. (Φ ST = 0.02, P = 0.23; Table 2). Similarly, no genetic differentiation was detected in any combinations of pairwise F ST values among the four populations (Table 3a). Table 1. Sampling localities and population genetic diversity indices (COI gene) of Monodonta sp. nN: number of individuals used data of nDNA obtained by ddRAD-seq analyses; mtN: number of individuals used mitochondrial population genetic analyses (COI); nH: number of haplotypes (COI); HD: Haplotype diversity (COI); ND: Nucleotide diversity (COI); IC: inbreeding coefficient estimated by using nDNA.  Genetic analyses based on the SNPs datasets. After demultiplexing the raw datasets, we obtained 675,370-10,388,536 reads for each individual. A total of 1731 filtered loci (9.72% of missing sites) and a total alignment length of 158,328 bp (9.20% of missing sites) were obtained (details in Table S1). The results of the STRU CTU RE analysis are shown in Fig. 3, and the maximum value of ΔK was at K = 2. These two genetic clusters were mainly divided by the Mukojima Islands [Kitanoshima Island (locality 1) and Yomejima Island (locality 2)] and the southernmost Hahajima Island (locality 4). On Chichijima Island (locality 3), which is located centrally in the Ogasawara Islands, the two genetic clusters were mainly mixed. A scatterplot of the PCA also showed that individuals fell into two genetic clusters (Fig. S1). This trend was similar to the results of the aforementioned STRU CTU RE analysis. In PCA plots [ Fig. S1a (PC 1 and PC 2), b (PC 1 and PC 3)], the northern Mukojima Islands (locality 1 and 2) and several individuals of Chichijima Island and Hahajima Island came together in a single cluster. Other individuals of Chichijima Island and Hahajima Island formed a distinct cluster. The above two clusters did not appear in PC2 and PC3 (Fig. S1c). The BA3-SNPs analysis indicated the gradient of gene flow among populations (Fig. 4), which did not necessarily correspond to the geographic distance between islands. Migrations from Yomejima Island and Hahajima Island to Chichijima Island tended to occur relatively easily (> 0.10), while those from Chichijima Island to other islands were less likely to occur (< 0.04). Migrations between Yomejima Island and Hahajima Island were infrequent. The inbreeding coefficient exceeded 0.10 at each locality except at Kitanoshima Island, and was about 0.5 at Chichijima Island ( Table 1).
The AMOVA estimated by using these SNPs dataset showed the presence of population genetic differentiation (Φ ST = 0.219, P < 0.0001; Table 2). Pairwise F ST values between the populations were higher than the results from mtDNA (Table 3). Except for pairwise F ST values between Chichijima Island and Hahajima Island, genetic differentiation was detected between the population.
The estimating the timing of genetic differentiation by fastsimcoal2 approach showed that the northern (Kitanoshima Island and Yomejima Island) and the southern (Chichijima Island and Hahajima Island) islands

Discussion
The present study tested the level of population genetic differentiation of an intertidal gastropod species distributed within a narrow range of peripheral oceanic islands using two types of genetic markers. A genome-wide single nucleotide polymorphisms (SNPs) dataset obtained by double-digest RAD sequencing (ddRAD-seq) was able to detect genetic differentiation among populations of Monodonta sp. endemic in the Ogasawara Islands.
In the case of our study, the mitochondrial DNA (mtDNA) did not provide sufficient information on genetic polymorphisms compared to the above high resolution genetic marker. For molecular ecological studies, mtDNA has been commonly used to demonstrate the population genetic structure and evolutionary history of various taxa. However, mitochondrial-specific traits such as frequent introgression may bias the investigation of population genetic structuring and dynamics 47 . Incomplete lineage sorting also makes the interpretation of the genetic results difficult. The discordance often occurs between mtDNA and nuclear DNA (nDNA) due to their different inheritance patterns 48,49 and mtDNA seems to lose the ancestral polymorphisms faster than nDNA 50,51 . In contrast, genetic traces left in mtDNA are important for understanding the complex evolutionary history between organelle and nuclear genomes 52 . By using both organellar DNA and genome-wide datasets, we can understand the population genetic dynamics in detail. From the above context, the present mtDNA results provided useful information on the differences in genetic diversity between Monodonta sp. on the Ogasawara Islands and its closely related continental species, M. confusa. The mtDNA of M. confusa displays a higher level of genetic   54 . These results indicate that endemic marine species underwent a strong bottleneck when colonising the Ogasawara Islands, which is considered a fundamental pattern in peripheral oceanic islands 28,40 .
Populations of Monodonta sp. accumulated genetic differentiation within a relatively narrow range from Kitanoshima Island to Hahajima Island (< 120 km). The genetic components of Monodonta sp. are divided between the Mukojima Islands (Kitanoshima Island and Yomejima Island) and the southern islands (Chichijima Island and Hahajima Island). The genus Monodonta is known to undergo a planktonic larval phase during its life cycle. In the case of M. confusa, which is the closest continental relative to Monodonta sp. endemic to the Ogasawara Islands, the planktonic duration is 3 days 45 . Besides, trochid snail species related to Monodonta have relatively short PLD (Gibbula umbilicalis and Osilinus lineatus ≤ 7 days) 55 . Dispersal ability during the planktonic larval phase is thought to be one of the most important factors for genetic differentiation among populations in the marine environment. However, a theoretical study indicated that population genetic connectivity is not necessarily determined by planktonic larval dispersal duration 56 . Isolated island populations have experienced significant founder effects and bottlenecks when colonising and could not receive migrants from the source continental populations. The genetic drift is thus important for the genetic differentiation of island endemic species 57,58 .
The geographical distance at which genetic differentiation occurs depends on the oceanographic conditions and specific characteristics of species 59 . Although the Kuroshio Current and its countercurrent are effective factors in the dispersal of marine organisms distributed in the North Pacific ocean 60 , the detailed system of oceanic currents around the Ogasawara Islands has not been sufficiently studied 61,62 . In seed dispersal plants Pandanus boninensis (Pandanaceae) distributed in the Ogasawara Islands, molecular studies with gene flow estimations demonstrated that migration mainly occurs from the south islands (Hahajima Island) to the north islands (Mukojima Island) via Chichijima Island and also from the Mukojima to the Chichijima Islands 62 . The result of gene flow estimation of Monodonta sp. also showed that the opportunities for migration to Chichijima Island from Yomejima and Hahajima Islands tended to occur relatively frequently, while reverse migration occurred less frequently. Although migration via oceanic currents vary among taxa around the Ogasawara Islands, there may be currents heading toward Chichijima Island from the south and north islands. In addition, frequent migration to Chichijima Island may also be due to the habitat preference of Monodonta sp. While Monodonta sp. is common in the rocky intertidal zone of the Ogasawara Islands, their inhabitation is less common in environments that are exposed to strong waves 63 . Chichijima Island may hold many suitable habitats for Monodonta sp. since it is the largest of the Ogasawara Islands and has a long coastline and large bay. Furthermore, the Chichijima Islands are composed of three large islands ( Fig. 1; Otoutojima Island, Anijima Island, and Chichijima Island), and each island can weaken the effects of waves on the intertidal zone. In fact, the inbreeding coefficient of Chichijima Island is at a high level. The presence of such suitable habitats may assist in the settlement of larvae that migrate from other islands to the Chichijima Islands. Monodonta sp. is not distributed to islands surrounded by coastlines with strong wave exposure, such as Nishinoshima Island (approximately 130 km from Chichijima Island) and Kita-iwoto Island (approximately 150 km from Hahajima Island) 64,65 . This suggests that planktonic larvae of Monodonta sp. are unable to settle in an environment attacked by strong waves or it may be difficult to disperse over distances of 100 km.
The life cycle of Monodonta sp. on the Ogasawara Islands is not well understood, but the knowledge of M. confusa can be applied to estimate the results of population demographic estimation [42][43][44] . According to these studies, M. confusa matures in 1 to 2 years. If we apply the above generation time to interpreting the population demographic estimation of Monodonta sp. on the Ogasawara Islands (Fig. 5), the bottlenecked ancestral population size have recovered 454,514-909,028 years ago (t1). Therefore, the time of divergence of Monodonta sp. from its continental sister species is considered to be before t1, and there is no significant discrepancy between t1 and   68 . Following LGM, because the expansion of the intertidal zone was completed and the sea level began to rise, genetic differentiation of Monodonta sp. may have been initiated. Although the present pattern of genetic structure, which includes two different genetic components in the two most external islands and a mixture of the two in the central islands, seem to fit the stepping stone model, it is not confirmed by gene flow estimation. Therefore, this pattern was established by genetic differentiation between the north and south islands due to the above discussed past climate event and is thought to be regulated by the oceanic current moving to Chichijima Island. In addition, population genetic differentiation of endemic species may be promoted by the disadvantages of dispersals in remote environments like oceanic islands isolated from the continental landmass [69][70][71] . Larval dispersal is considered to have several advantages, including colonising new areas and decreasing the likelihood of inbreeding and extinction due to genetic exchanges between populations 69,71-73 . Dispersal also fundamentally has various unavoidable risks. Even if Monodonta sp. maintains similar dispersal ability of its continental relatives, the dispersed larvae may not be able to settle in suitable habitats on isolated oceanic islands. If few suitable habitats are available, dispersed larvae may fail to settle successfully, thereby increasing mortality 74 . For the intertidal snail species genus Littorina, direct developing species are present on isolated and remote islands, while there are no closely related species with a pelagic larval stage (The paradox of Rockall) 75 . The possibility that isolated situations may promote genetic differentiation due to the disadvantages of larval dispersal is an interesting issue to be solved by future ecological studies on larval behaviours.
The marine species inhabiting the Ogasawara Islands are mainly characterised by taxa that once successfully colonised these islands after long-distance travel. Endemic phenomena in oceanic islands are affected by the distance of the island and time elapsed since colonisation 27 . Around the Ogasawara Islands, several marine species were confirmed to be endemic [33][34][35][36][37][38] . However, the level of genetic connectivity and differentiation within the Ogasawara Islands is poorly understood in other marine species. In the case of coral species with far-reaching dispersal ability, a low level of genetic differentiation among the three populations on the Ogasawara Islands was detected by microsatellite markers 76 . In contrast, Monodonta sp., with a relatively short larval phase, exhibited a genetic structure within this narrow range of archipelagos. Similarly, other endemic marine species with short dispersal ability may also show genetic differentiation within a similar range. To understand the connectivity among populations of the Ogasawara Islands, comparative genetic studies using multiple taxa are needed. The present study indicated that it is necessary for marine conservation on the Ogasawara Islands to focus on both the north and south island populations which differ genetically. Besides, the genetic structure of marine species may be more differentiated than previously thought within the Ogasawara islands. High resolution genetic analyses of various marine taxa will provide invaluable insights into marine conservation strategies in highly endemic water regions such as oceanic islands 77,78 . By clarifying the genetic structure and population demographic history of various species, we can set up a marine reserve system at an appropriate geographic scale 79,80 .
In conclusion, we detected the population genetic divergence of intertidal gastropod with a planktonic larval phase, even within a narrow range of the Ogasawara Islands. This was caused by historical climatic events and maintained by temporal oceanographic conditions. To accurately grasp the degree of genetic differentiation and population history, high-throughput genetic data is useful. In the genomic era in which we live, unexpected genetic differentiation is bound to be detected. These findings provide important insights for understanding the population genetic patterns of open and connectable environmental situations observed in remote oceanic islands.

Methods
Sample preparation. In total, 59 samples of Monodonta sp. were collected from the four localities of the Ogasawara Islands ( Fig. 1; Table 1). Our sampling range covered the distribution area of Monodonta sp. on the Ogasawara Islands. To prepare the samples, a part of the foot muscle was dissected from each individual and stored in 99.5% ethanol for subsequent molecular analyses.
Mitochondrial DNA sequencing and population genetic analyses. Total DNA was obtained from tissue samples using the NucleoSpin Tissue (TaKaRa, Shiga Pref., Japan) according to the manufacturer's instructions. Fragments of the COI gene were amplified using the primers CoxAF (5′-CWA ATC AYA AAG ATA TTG GAAC-3′) and CoxAR (5′-ATA TAW ACT TCW GGG TGA CC-3′) 81 . Polymerase chain reactions (PCR) were performed under the following conditions: 94 °C for 3 min followed by 5 cycles at 94 °C for 30 s, 45 °C for 30 s, 72 °C for 1 min, followed by 35 cycles at 94 °C for 30 s, 52 °C for 30 s, and 72 °C for 1 min, with a final extension at 72 °C for 5 min. The products were then purified using ExoSAP-IT (Amersham Biosciences, Buckinghamshire, UK). Cycle sequencing was performed using the PCR primers with the BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, CA, USA) and the products were directly sequenced from both directions using an ABI 3130xl automated sequencer (Applied Biosystems). The validity of the sequences was verified with the software package 4Peaks 82 , and the forward and reverse sequences were assembled using  84 . Sequence data are available on Genbank (LC316340-LC316346, LC671974-LC672025). Two genetic diversity indices (i.e., haplotype diversity and nucleotide diversity) were calculated for each population using Arlequin v. 3.5 85 . The population genetic structures were estimated using analysis of molecular variance (AMOVA) 86 . Pairwise F ST between populations was calculated by Arlequin v. 3.5 with 1000 permutations. To visualize the geographical distribution pattern of haplotypes, haplotype networks were reconstructed using a median-joining network 87 implemented in PopART 88 . ddRAD-seq, SNPs detection, and genetic analyses. To prepare a library of double digest restriction site-associated DNA sequencing (ddRAD-seq), 29 genetic samples of Monodonta sp. were used (Table 1). RNase was added to the total DNA and it was digested by two restriction enzymes (EcoRI and MspI). P1 and P2 adapters were ligated to DNA fragments according to the ddRAD-seq protocol 89 . The ligated samples were multiplexed and purified with the NucleoSpin gDNA Clean-up kit. Pippin Prep (Sage Science, MA, USA) was used to collect approximately 450 base pair (bp) DNA fragments. The DNA fragments were amplified in eight single PCR reactions, and the products were cleaned using the NucleoSpin gDNA Clean-up kit. The constructed DNA library was sent to Oregon State University's Center for Genome Research and Biocomputing and sequenced using Illumina HiSeq 4000 single-end sequencing, yielding maximum read lengths of 100 bp. Raw sequence data were deposited in the DDBJ Sequence Read Archive (DRA013313). The raw Illumina sequence reads were demultiplexed and processed with the ipyrad pipeline 90 . If the obtained reads were low in quality or contained adapter sequences, it was trimmed and set to a minimum length of 35 bp using ipyrad's parameter setting. Reads were clustered at 85% sequence similarity. The minimum number of samples was set to 26/29 (0.896%). Other parameters followed the default setting of the ipyrad pipeline and output datasets for subsequent genetic analysis.
To grasp genetic population structure, we estimated individual genotypes using STRU CTU RE v. 2.3.4 according to the ipyrad analysis toolkit (https:// ipyrad. readt hedocs. io/ en/ latest/ index. html). The number of preassigned genetic clusters (K) was assumed to range from 1 to 10, and seven independent runs were performed for each K value. Each run included 50,000 burn-in iterations and 100,000 iterations. We also performed a principal component analysis (PCA) with GenoDive v. 3.03 91 and PCAs were plotted by PAST 4.04 92 . Temporal gene flow rates and inbreeding coefficient among four populations were estimated by BA3-SNPs 93 . The delta values (− m, − a, and − f) selected were 0.3250, 0.5500, and 0.0750, respectively. The number of runs were 30,000,000 MCMC iterations including 6,000,000 burn-ins (20%) and 1000 sampling intervals. Tracer v. 1.6.0 was used to verify the Bayesian convergence and parameter values throughout generations 94 . The AMOVA and pairwise F ST between populations were calculated with Arlequin v. 3.5.
The above results of genetic analyses showed that the populations of the northern (Kitanoshima Island and Yomejima Island) and the southern (Chichijima Island and Hahajima Island) islands are genetically differentiated (see "Results" section), and we then estimated its divergence timing by fastsimcoal2 with simplified demographic models (Fig. 5) 95 . To grasp the demographic history of Monodotna sp. in the Ogasawara Island, we calculated the two time series: the recovery timing of ancestral population size from the bottleneck that occurred when arriving at the Ogasawara island (t1) and the timing of genetic differentiation (t2) between the northern (Kitanoshima Island and Yomejima Island) and the southern (Chichijima Island and Hahajima Island) islands. The pilot run with a broad range of parameters was conducted and the parameter settings were determined according to the results of pilot run (Table S2). We calculated the site frequency spectrum using easySFS (https:// github. com/ isaac overc ast/ easyS FS). We carried out 100 independent runs of 100,000 coalescent simulations in each model. Then, we used a simulation with the highest likelihood values for parameter estimation.