Divergent northern and southern populations and demographic history of the pearl oyster in the western Pacific revealed with genomic SNPs

Abstract In the open ocean without terrain boundaries, marine invertebrates with pelagic larvae can migrate long distances using ocean currents, suggesting reduced genetic diversification. Contrary to this assumption, however, genetic differentiation is often observed in marine invertebrates. In the present study, we sought to explain how population structure is established in the western Pacific Ocean, where the strong Kuroshio Current maintains high levels of gene flow from south to north, presumably promoting genetic homogeneity. We determined the population structure of the pearl oyster, Pinctada fucata, in the Indo‐Pacific Ocean using genome‐wide genotyping data from multiple sampling localities. Cluster analysis showed that the western Pacific population is distinct from that of the Indian Ocean, and that it is divided into northern (Japanese mainland) and southern (Nansei Islands, China, and Cambodia) populations. Genetic differentiation of P. fucata can be explained by geographic barriers in the Indian Ocean and a local lagoon, and by environmental gradients of sea surface temperature (SST) and oxygen concentration in the western Pacific. A genome scan showed evidence of adaptive evolution in genomic loci, possibly associated with changes in environmental factors, including SST and oxygen concentration. Furthermore, Bayesian simulation demonstrated that the past population expansion and division are congruent with ocean warming after the last glacial period. It is highly likely that the environmental gradient forms a genetic barrier that diversifies P. fucata populations in the western Pacific. This hypothesis helps to explain genetic differentiation and possible speciation of marine invertebrates.


| INTRODUC TI ON
A fundamental question regarding ecological and evolutionary processes of marine species is how differentiation and diversification occur in the open ocean. In marine environments, ocean currents advect pelagic organisms from one locality to another, allowing wide distributions. Marine animals with long pelagic larval periods can potentially disperse hundreds or thousands of kilometers (Shanks, 2009;Shanks, Grantham, & Carr, 2003). This long-range dispersal enhances connectivity of local populations, maintains gene flow, and large effective population size, thereby reducing genetic differentiation among populations (Palumbi, 1994;Ward, Woodwark, & Skibinski, 1994).
Dispersal of marine organisms is, however, confined by geographic factors such as ocean currents and terrain boundaries, and by environmental conditions such as sea surface temperature (SST), oxygen concentration, and salinity. Steady ocean currents restrict the direction of gene flow between localities, resulting in spatial patterns of population structure (Thomas et al., 2015;Xuereb et al., 2018). Marine population can also be dissociated by land exposure.
In particular, marine regression that accompanied the last glacial maximum (LGM) around 23,000-18,000 years ago, affected extant population structures (Hewitt, 2000). Furthermore, there is growing evidence that environmental gradients also serve as barriers that prevent genetic flow between local populations (Benestan et al., 2016;Bernatchez et al., 2019;Sandoval-Castillo, Robinson, Hart, Strain, & Beheregaray, 2018;Sylvester et al., 2018). Larvae delivered by ocean currents can settle only new localities that provide suitable environmental conditions. Therefore, two populations adapted to different environmental conditions are expected to be genetically isolated. Environmental heterogeneity also drives differentiation in genomic loci related to local adaptation (Guo, Li, & Merilä, 2016).
In order to understand the impact of environmental factors on differentiation of marine populations, the western Pacific is a fascinating region in which geographic conditions are particularly well studied (Gallagher et al. (2015) for review). The steady Kuroshio Current carries warm water from tropical to subtropical and temperate areas around the Nansei (Ryukyu) and Japanese Archipelagos ( Figure 1). Along the Nansei Archipelago, islands are located at less than 100-km intervals, except for a 200-km gap between Kume and Miyako Islands. This continuous chain of islands, connected by the Kuroshio Current (mean maximum surface velocity reaches 1.2 m/s (Yang et al., 2015)), allows long-range distribution of littoral and shallow-water marine organisms from the Nansei Archipelago to Kyushu, Shikoku, Honshu, and adjacent islands (hereafter we refer to these islands as the "Japanese mainland" for simplicity). In addition, there has been no drastic geological change or ocean current variation since the last glacial period (LGP) in the western Pacific (Gallagher et al., 2015;Voris, 2000). Physical geography, including a continuous chain of islands and ocean currents may have contributed to genetic homogeneity of marine populations in this region.
In fact, some intertidal gastropod species exhibit uniform populations in the western Pacific (Kojima et al., 2006;Yamazaki et al., 2017). It is noteworthy that an intertidal gastropod species Monodonta confusa, having only a 3-day pelagic larval duration (PLD), shows homogeneity from the Nansei Islands to the Japanese mainland (Sasaki, 1985;Yamazaki et al., 2017). This suggests that the Kuroshio Current can transport larvae sufficiently to maintain longrange population connectivity. In contrast, genetic discontinuity between the Nansei Islands and the Japanese mainland is evident in other animals with various life cycles and PLDs, such as the ostracod Vargula hilgendorfii (free-living adult, but poor swimming ability) (Ogoh & Ohmiya, 2005), the soft shore barnacle Fistulobalanus albicostatus (PLD = 5-7 days) (Chang et al., 2017), intertidal gastropods Mondosa perplexa and M. labio (PLD = 3 days for M. perplexa) (Sasaki, 1985;Yamazaki et al., 2017), and a mudskipper Boleophthalmus pectinirostris (PLD = 1 month) (Chen, Hong, Chen, Wang, & Zhang, 2015). These observations demonstrate that genetic structure varies among species and reflects more than just pelagic duration. One possible explanation for genetic structure under extensive oceanographic connectivity is that organisms have capabilities specific to environmental conditions; therefore, population structure is restricted by environmental gradients. Despite a number of population genetics studies that have been carried out previously (Chang et al., 2017;Chen et al., 2015;Kojima et al., 2006;Ogoh & Ohmiya, 2005;Sasaki, 1985;Yamazaki et al., 2017), there are no empirical studies that explain population structure of marine organisms in the western Pacific in light of environmental factors, based on genome-wide variations.
The pearl oyster, Pinctada fucata, is a good model to understand genetic differentiation in the western Pacific because of its wide distribution. P. fucata is a marine bivalve commonly cultured for pearl production, mainly in Asian countries. Since this species is one of the most valuable aquaculture species, assessment of the genetic structure of natural populations should be highly beneficial for conservation of this resource. Pinctada species are typically found in the Indo-Pacific and western Atlantic tropical and subtropical areas (Cunha, Blanc, Bonhomme, & Arnaud-Haond, 2011), while P. fucata is distributed from low-latitude tropics to north temperate regions around Japan, representing the northern distributional limit of the genus. P. fucata spends 2-3 weeks in its pelagic larval stage during early development (Fujimura, Wada, & Iwaki, 1995). Then it settles on rocky substrates in coastal areas facing the open ocean. Its broadcast spawning strategy, as well as its broad distribution, explains its high-distance dispersal capability (Yu & Chu, 2006).
To date, the population structure of P. fucata in the western Pacific remains controversial. An early biogeographic study demonstrated that the P. fucata population in the Nansei area differs from that in the Japanese mainland in frequency of shell color morphs, shell shape, and allozyme variation (Wada, 1984). A DNA marker study using the inter-simple sequence repeat (ISSR) showed that one PCR fragment was frequently found in individuals from the Nansei Islands, China, Cambodia, and Myanmar, while it was absent from the Japanese mainland (Masaoka & Kobayashi, 2006). Contrarily, no genetic marker was found to distinguish individuals from Japan, China, Vietnam, and Thailand using allozyme loci (Atsumi, Komaru, & Okamoto, 2004). Molecular phylogenetic analyses using complete sequences of 18S and 28S ribosomal RNA and ITS regions did not show significant differentiation among Japanese, Chinese and Cambodian populations (Masaoka & Kobayashi, 2002. In addition, artificial fertilization experiments demonstrated that individuals from western Pacific countries produce normal veliger larvae, indicating that they are not isolated reproductively (Atsumi et al., 2004). In summary, population structure may be possible, but genetic analyses with more high-resolution molecular markers and more comprehensive  sampling are necessary to resolve the population structure of P. fucata in the Indo-Pacific Ocean.
In this study, we first conducted restriction site-associated DNA (RAD) sequencing and retrieved genome-wide SNPs from P. fucata specimens throughout the Indo-Pacific region to resolve the population structure. Next, in order to identify possible barriers that limit gene flow in the western Pacific, we tested correlations between genetic distance and various environmental factors, including SST, oxygen, salinity, phosphate, nitrate, and carbon dioxide concentrations. Using outlier analysis, genomic loci and genes putatively under non-neutral evolution were investigated. Finally, we integrated the above results, Bayesian simulations, and admixture analysis to reconstruct the demographic history of P. fucata that formed extant populations after the LGP. The results presented here will provide a striking example of genetic differentiation of organisms with high dispersal ability in the open ocean.

| Sample collection
We collected pearl oysters, Pinctada fucata, from five regions in East

| Data processing and genotyping
After sequencing, raw reads with restriction enzyme recognition sites and barcode sequences at both ends were retained with bcl2fastq software (Illumina). Reads were then demultiplexed using reaper and tally, parts of Kraken software tools (Davis, van Dongen, Abreu-Goodger, Bartonicek, & Enright, 2013), based on custom barcodes to give focused genomic sequence data derived from each individual. Then, adapters and sequences with low quality (average quality below 30 with a 3-base sliding window) were trimmed with Trimmomatic (version 0.36) (Bolger, Lohse, & Usadel, 2014).

| Population structure analyses
We analyzed population structure using model-free approaches.
Population structure of P. fucata was further investigated with a model-based approach using a Bayesian clustering method implemented in fastSTRUCTURE (Raj, Stephens, & Pritchard, 2014). We tested K = 2 to K = 10 clusters using the simple prior model. Then we chose the number of clusters (K = 4) that minimized marginal likelihood using the "chooseK.py" script provided in the fastSTRUCTURE package. We also conducted discriminant analysis of principal components (DAPC) in R using the adegenet library (Jombart, 2008). In order to determine to which cluster each specimen from the western Pacific belonged, samples from MMR were excluded and the number of clusters K = 3 was applied to the DAPC. The number of principal components was chosen to optimize the α-score estimated by "optim.a.score" module in adegenet.
The Bayesian program BayeScan was used to search for loci with extreme F ST values that should be explained by selection (Foll & Gaggiotti, 2008 OutFLANK estimates a likelihood based on a trimmed distribution of F ST to infer the neutral distribution of F ST (Whitlock & Lotterhos, 2015). Highest and lowest outlier trimming was set to 0.05 and the minimum heterozygosity required for inclusion was set to 0.1, as the default parameters. For each analysis, loci that satisfied a threshold of the false discovery rate (FDR) or q-value of 0.1 were retained.
Finally, loci detected by all three methods were considered outliers.
Genomic regions around outlier loci were surveyed to find coding sequences that might be under selection. We searched gene models located within 1 kb of outlier loci in the P. fucata reference genome version 2.0 (Takeuchi et al., 2016). We used a 1-kb window for physical linkage between SNPs and gene models based on the observation that oysters show a low level of linkage disequilibrium (LD) (Gutierrez et al., 2017). Putative gene models that are subject to selection were searched with BLASTP against the UniProtKB database using an Evalue threshold <1e −5 . UniProt IDs were retrieved from the BLASTP results and used for Gene Ontology annotation (Ashburner et al., 2000) using PANTHER (Mi, Muruganujan, Ebert, Huang, & Thomas, 2019).

| Assessment of isolation by environmental factors
We assumed that environmental gradients shape P. fucata popula- concentration, an annual average over 19 years was used. These ecological values were scaled for further analyses. Geographic distances between localities were calculated using the Imap package in R (Wallace, 2012).
Since collinearity between variables can affect RDA (Dormann et al., 2013), correlation among variables, including environmental factors and geographic distance, was tested. First, we assessed correlation between geographic distance and annual mean differences in environmental factors among localities ( Figure S2). We selected localities located along the Kuroshio Current (JPN_1, 2, 3,7,8,9,10,13,NAN_18,19,20,and 21) because direct distances between JPN localities in Pacific coast and in Japan Sea coast do not represent the distances that marine organisms actually migrate. Simple Mantel tests for distance and ecological factors were assessed using the R ecodist package (Goslee & Urban, 2007) and Pearson correlation with 1,000 permutations. Because distance and environmental factors other than CO 2 are significantly correlated (see Figure S2 for detail), distance was omitted from RDA. Effect of isolation-by-distance was tested independently (see Note S1 for detail). Second, correlations between each pair of environmental factors were checked by estimating Pearson's correlation coefficient in R. Results showed that sea surface temperature (SST) and oxygen concentration are strongly correlated throughout the year ( Figure S3). It is reasonable to expect a negative correlation between SST and oxygen concentration because oxygen solubility in water increases at lower temperatures. Thus, we prepared two datasets: (a) without oxygen and (b) without SST in order to analyze SST and oxygen independently. In addition, because there were correlations between seasons for SST, oxygen, and salinity, we separated the datasets for each season.
Consequently, we analyzed eight environmental datasets. For genotyping data, we prepared two datasets: (a) all SNPs and (b) non-neutral SNPs identified using outlier analyses in order to assess which environmental factors are related to the SNPs under selection. We conducted analysis of variance (ANOVA) with 1,000 permutations to assess the global significance of the model and to evaluate the influence of each environmental factor. To control FDR in the multiple testing, q-values were calculated using the Benjamini-Hochberg (BH) method (Benjamini & Hochberg, 1995).
Using RDA, we also investigated SNP loci putatively related to local adaptation, based on the method of Forester, Lasky, Wagner, and Urban (2018). We conducted RDAs with the all-SNP dataset for eight environmental datasets as mentioned above. Two RDA axes were taken into account for SNP search because at most two environmental factors were significantly correlated with genetic variations in each dataset (see the Results section for detail). Loadings were calculated for each SNP locus, and loci that showed loading ± 3SD from the mean loading of RDA axes were regarded as associated with local adaptation.

| Modeling scenario for establishment of population structure
In the western Pacific, northern and southern populations were identified (see Results). Next, we focused on how the two populations separated after the last glacial period (LGP). One premise is that P. fucata originally inhabited the southern region, and then northern population appeared afterward, because (a) the genus Pinctada is mainly distributed in tropical and subtropical areas at present (Cunha et al., 2011), showing that Pinctada is a warm-water taxon, (b) fossils of P. fucata are absent during the LGP in Japan (Masuda & Noda, 1976) and their shells are found in middens created about 6,000 years ago in the Japanese mainland (Yamamoto, 1992), indicating that the northern population appeared after the LGP, and (c) the northward direction of the Kuroshio Current has been stable since the LGP (Gallagher et al., 2015), making it impossible for pearl oyster larvae to have passively migrated from north to south. In order to understand how the two populations become established, we hypothesized two scenarios: expansion (scenario 1) and migration (scenario 2). In scenario 1, the ancestral population in the southern area expanded to reach the northern area. Then the southern and northern populations become divided, assuming that neither reduced its population size significantly during the separation. In scenario 2, the size of the southern population is consistent with that of the ancestral population, and a small number of individuals incidentally migrated and settled in the northern area. Unlike scenario 1, scenario 2 assumes that the present northern population was derived from a limited number of founders, experiencing a founder effect.
We investigated these hypotheses using approximate Bayesian computation (ABC) analysis (Beaumont, 2010) implemented in DIYABC (Cornuet et al., 2014). ABC is a flexible method for drawing statistical inference about complex evolutionary scenarios based on genetic data. ABC compares summary statistics of observed and simulated data under evolutionary scenarios, and then infers parameters without explicit likelihood calculations (Beaumont, 2010). We generated 1 million simulated datasets representing demographic scenarios. The simulated datasets were compared to empirical data to assess the likelihood of each scenario. Summary statistics were calculated for the proportion of zero values, the mean and variance of nonzero values, and the mean of the complete distribution for each population. In addition, summary statistics for F ST and Nei's distances between populations were used.
The posterior probability of each scenario was calculated using a direct approach (Cornuet et al., 2008) using the 500 closest simulated datasets.

| Inference of population admixture
We inferred introgression of one population into the other using TreeMix version 1.13 (Pickrell & Pritchard, 2012). The maximumlikelihood tree was constructed using an SNP window size of 100 and the MMR population as an outgroup. The number of tested migration events varied from 1 to 5. The model fit was quantified by the covariance for each migration event. The plotting function in R was used to visualize the results.

| Population structure
Principal component analysis (PCA) of all samples showed a distinct population structure of P. fucata in the region (Figures 2 and S1).
The first two principal components explained 11.7% and 7.6% of the total variance, respectively. The population at MMR was clearly dis-
Assuming that genes physically close to these loci are under non-neutral selection, we identified 42 gene models located within 1 kb of non-neutral loci in the reference genome. Among those, 35 gene models were annotated based on results of a BLASTP search against the UniProt database (Table S3). GO annotation showed that 34 genes were classified into six categories of molecular function ( Figure 5b). Catalytic activity (10 genes, 47.6%) and binding (4 genes, 19.0%) occupied two-thirds of the class. The predominant category in biological process is cellular process (6 genes, 25.0%), followed by metabolic process (5 genes, 20.8%), localization (3 genes, 12.6%) and biological regulation (3 genes, 12.5%). In the cellular component class, most genes were categorized as organelle (7 genes, 35.0%).

| Isolation by environmental factors
Although northern and southern populations are distinct from each other in the western Pacific, there is no visible terrain boundary separating the populations. We next investigated environmental variables that may disrupt gene flow between the southern and northern populations.
To test the correlation between genetic variation and environmental variables, a series of RDA was performed ( Figure 6, Table S4). When analyzing all SNP datasets (Figure 6a-h), genetic variation was significantly associated with SST in spring, summer, and winter (Figure 6a Based on the RDA results, SNP loci possibly associated with adaptation to local environments were investigated (Table S5). In total, 1,049 unique loci were identified and 11 of them were mutually detected by three outlier analyses, BayeScan, PCAdapt, and OutFLANK.
We searched genomic regions near the 11 loci (≤1 kb) and found 4 gene models that may be related to adaptation (Table S3).

| Population history
Next, we addressed the demographic history of P. fucata that has shaped current population structure after the LGP in the western

F I G U R E 2 Distinct population structure of Pinctada fucata in the Indo-Pacific area. Principal component analysis (PCA) of all samples. Each dot indicates one individual
Pacific. Using DIYABC, we tested two hypotheses: expansion (scenario 1) and migration (scenario 2) ( Figure 7a). Prior scaled parameters, including effective population sizes (N) and timing of population events (t) were set as shown in Table 2, with uniform distributions.
Our results showed that the posterior probability (PP) and the confi-

| D ISCUSS I ON
Population structure of marine invertebrates is determined by biological attributes such as pelagic larval duration and physiology, as well as by seascape features, including geographic setting and environmental variables.
Terrain boundary is one of the apparent barriers that prevents gene flow among marine populations, resulting in vicariance and allopatric differentiation. For example, the Central American Isthmus, which emerged in the late Pliocene (2.8 million years ago), separated the Pacific and the Atlantic Oceans (Coates, McNeill, Aubry, Berggren, & Collins, 2005). Marine organisms diverged into genetically distinct sister clades on both sides of the Isthmus (Lessios, 2008). In the Indo-Pacific area, a large land mass, corresponding to the Sunda Shelf, was exposed by low sea level during glacial periods of the Pleistocene (Voris, 2000). The appearance of this land mass severed the genetic connection between Indian and Pacific marine populations with long-lived pelagic larvae, such as starfish and reef fish species (Benzie, 1999;Gaither et al., 2011). The pearl oyster population in MMR may also have been differentiated from the western Pacific population by this terrain boundary (Figures 2,   3, and 4).
The pearl oyster population at Kamikoshiki Island (JPN_17) exemplifies small-scale, but distinct allopatric differentiation by a geographic barrier. Pearl oysters inhabit a coastal lagoon called Lake Namako-ike on the shore of the island. The lake is completely isolated from the open sea by a coastal gravel bar with no tidal inlets. The geological barrier, established in mid-Holocene about 5,000-6,000 years ago (Aramaki, Yamaguchi, & Tanaka, 1976;Woodruff, Donnelly, & Okusu, 2009), genetically isolates the Kamikoshiki population from all others, resulting in unique genetic constitution (Figures 2, 3, and 4a).
Except for the population in Kamikoshiki Island, genetic divergence among localities in the Japanese mainland is low (Figure 3a).
One possible explanation is that the pearl oysters have been artificially delivered from one place to another for aquaculture purpose, causing genetic homogenization in the JPN region. However, our sampling localities include islands (JPN_1, 11, and 14) more than 100 km from pearl farming areas (Figure 1b). To the best of our knowledge, there are no records of artificial transfers of pearl oysters in these areas; therefore, local pearl oyster populations of isolated islands cannot have been contaminated by human activities. Thus, the genetic uniformity of population in JPN localities is attributed to natural dispersal of the pearl oyster throughout the region by the Kuroshio and Tsushima Currents (Figure 1a), irrespective of artificial intervention. Genetic uniformity negates the possibility of isolation-by-distance to explain the P. fucata population in the JPN region ( Figure S4). Similarly, the southern population is distributed across a wide area ranging from the South China Sea (KHM and CHN) to the East China Sea (NAN). The wide distribution of the population is maintained by the long-lived pelagic larvae, together with steady ocean currents. Population uniformity from the Gulf of Thailand to the South China Sea was demonstrated in various marine animals such as the seahorses, Hippocampus kuda and H. trimaculatus (Lourie, Green, & Vincent, 2005), the acorn barnacle, Chthamalus malayensis (Tsang et al., 2012), and the limpet, Cellana toreuma (Wang, Ganmanee, Shau-Hwai, Mujahid, & Dong, 2016). A genetic connection between the South China Sea and Nansei Islands populations is also evident from studies on the toothed top shell, Monodonta labio (Yamazaki et al., 2017). These studies support the idea that the South China Sea Current connects KHM and CHN, and that the Kuroshio and Taiwan Currents connect the South China Sea and the NAN area ( Figure 1a), maintaining the southern population of P. fucata.

| Environmental barriers and adaptation to local conditions
As P. fucata shows potential to inhabit much of the coastal habitat of the western Pacific, it is speculated that southern and northern populations could also be demographically connected by the Kuroshio Current. In fact, the presence of wild P. fucata population in the JPN region confirms that they can migrate from the NAN area to JPN via the Kuroshio Current, at least in the past. In the present ocean, however, the southern and northern populations are clearly differentiated (Figures 2, 3, and 4). In the absence of a geographic barrier, we hypothesize that existing population patterns of P. fucata can be explained by their capacity to adapt to local environmental gradients.
As pearl oysters adopt a sessile adult lifestyle, dispersion occurs in the pelagic larval stage from fertilization to settlement, taking 2-3 weeks. Maturation of gonads and spawning of P. fucata occurs throughout the year, with a small peak of activity in the winter in the Nansei area (Wada, Komaru, Ichimura, & Kurosaki, 1995). Therefore, the environmental barrier that limits dispersal of larvae from south to north should persist regardless of the season. Among the environmental factors studied here, SST showed the strongest correlation    with the genetic distance between the southern and northern populations ( Figure 6), a correlation that is evident in all four seasons.
The results indicate that the SST gradient is one of the major factors affecting the distribution of P. fucata populations in the western Pacific. In fact, seawater temperature directly affects metabolic and physiological activity of pearl oysters (Numaguchi, 1994;Yukihira, Lucas, & Klumpp, 2000). It is noteworthy that the optimum seawater temperature for P. fucata in China is significantly higher than that for individuals in Japan (Numaguchi & Tanaka, 1986;Wang, Zhu, Wang, Luo, & Liu, 2012). Furthermore, their growth is considerably affected by water temperature, especially in early life stages such as fertilization, development, and settlement (Numaguchi & Tanaka, 1986; O Thermal tolerance is associated with physiological oxygen supply and demand. Oxygen deficiency is caused by increasing oxygen demand at high temperatures, and by decreased ventilation and circulation at low temperatures (Pörtner, 2001(Pörtner, , 2002. Therefore, ectothermic marine organisms may modify functional properties and density of mitochondria to maximize ATP production efficiency at ambient temperature. In our outlier study, several genes with metabolic functions are identified. For instance, a gene encoding a subunit of succinate dehydrogenase (pfu_aug2.0_10.1_00001. t1), located in the mitochondrial matrix, participates in the electron transport chain, as well as the TCA cycle (Hederstedt & Rutberg, 1981;Yankovskaya et al., 2003). Inverted formin-2 or INF2 homolog (pfu_aug2.0_1964.1_15330.t1) is involved in mitochondrial fission (Korobova, Ramabhadran, & Higgs, 2013). Notably, this gene is detected not only by all 3 outlier approaches, but also by a genotype-environment association approach using RDA (Table S3). It is possible that gene expression level or activity of these proteins is optimized in response to environmental conditions.
Along with lower seawater temperatures, changes in mitochondrial properties and higher concentrations of oxygen in the environment may cause increasing oxidative stress by F I G U R E 7 Population history scenarios of Pinctada fucata in the western Pacific tested using DIYABC. (a) In scenario 1, the ancestral population (N A ) was increased at t 2 , and separated into two populations (N S and N N ) without a reduction in the sum of the population size at t 1 . Total population size is stable during the subdivision, and no bottleneck event was occurred. In scenario 2, the southern population (N S ) size is stable, while the northern population (N N ) originated from a small number of individuals from the ancestral population (N NA ) derived from N S at t 2 , and expanded at t 1 . A bottleneck in the northern population is assumed in scenario 2. (b) Posterior probability and confidence interval of the scenarios estimated using DIYABC  (Abele & Puntarulo, 2004;Guderley, 2004;Viarengo, Canesi, Garcia Martinez, Peters, & Livingstone †, 1995). ROS strongly oxidize all cell components including DNA, proteins, and lipids. The enzymatic defense system against cytotoxic effects of ROS includes superoxide dismutase, catalase, and glutathione peroxidase (Manduzio, Rocher, Durand, Galap, & Leboulenger, 2005). In the present study, however, we did not find any of these genes among adaptation-related candidates. On the other hand, there is a gene involved in DNA repair, single-strand selective monofunctional uracil-DNA glycosylase (pfu_aug2.0_1138.1_28083.t1).
In summary, the present outlier analysis detected a number of genes that may be related to local adaptation. It should be noted that the gene candidate list includes gene models of unknown relevance and that some outliers could be false positives. Meanwhile, information about outliers in the genome suggests potential candidates, the biological function of which should be analyzed in future studies.
Experimental approaches such as RNA-seq to detect differentially expressed genes in different conditions, gene knockdown, or genome editing will be essential to identify genes related to adaptation.

| Reconstruction of population history
Extant populations are derived from ancestral populations that experienced environmental changes in the Pleistocene (Hewitt, 2000).
Therefore, in order to understand how contemporary population structure was established, it is necessary to consider historical seascape features.
Based on the present population genetics study discussed above, we reconstructed the demographic history of P. fucata ( Figure 9). During the last glacial period, the SST around the Japanese mainland was 5-6°C lower than at present. The Kuroshio Current was thought to have turned westward at a lower latitude than today (Ujiié, Tanaka, & Ono, 1991;Ujiié, Ujiié, Taira, Nakamura, & Oguri, 2003), but see Lee et al. (2013) and Gallagher et al. (2015) for a review. During this period, the population may not have been able to inhabit the JPN area (Figure 9a). After the last glacial period, temperatures gradually increased and peaked ~6,000 years ago, 2-3°C higher than at present in Japan (Schöne et al., 2004). SST warming enabled the P. fucata population to extend northward and to reach the Japanese mainland during this period ( Figure 9b). As the SST gradually decreased to its present level, individuals able to tolerate lower water temperatures remained in the JPN area and formed the northern population (Figure 9c). Environmental gradients, especially SST, caused a genetic discontinuity between the modern northern and southern populations. Meanwhile, artificial fertilization experiments indicate that speciation has not occurred between the populations (Atsumi et al., 2004).
Population patterns along the Kuroshio Current in the western Pacific vary among organisms. In particular, diversity of coral species diminishes drastically in the Japanese mainland compared to F I G U R E 9 Possible population history of Pinctada fucata in the western Pacific. (a) In the last glacial period, the distribution of P. fucata was restricted to the South, due to low SST. The light gray area indicates exposed land, assuming that the sea level was 120 m below present. Note that this map is based on modern terrain height and that tectonic change is not considered. (b) The ancestral population expanded to the north, accompanied by increased SST. Tropical species such as hermatypic corals were also distributed to 35°N (dotted line). (c) As the SST gradually decreased to present levels, demographic flow from the South to the North became restricted by this environmental barrier. The northern population was retained, presumably due to adaptation to the local environment the Nansei Archipelago (Veron & Minchin, 1992). In other words, the Nansei area is the northern limit of distribution for many coral species, irrespective of their pelagic larval duration. Marine shallow-water species show common population structure, being divided into southern (Nansei Islands) and northern (Japanese mainland) populations (Chang et al., 2017;Ogoh & Ohmiya, 2005;Yamazaki et al., 2017). In this study, we demonstrated the population structure of a marine bivalve species. Despite their long-dispersal potential, P. fucata is genetically differentiated due to an environmental barrier in the western Pacific. Evidence from this study suggests that the present population distribution was shaped by past environmental changes and by the adaptive capability of the animal. This perspective can be generalized to explain varied population patterns of marine organisms that may have different thermal tolerances and evolutionary consequences. Fossil records indicate that 6,000 years ago, tropical fauna, including hermatypic corals and mollusk species were distributed at higher latitudes (~35°N) (Hoshino, 1967;Matsushima, 1984;Yabe & Sugiyama, 1932) (Figure 9b). Some tropical species retreated to the south as the SST decreased to present levels, while others colonized the JPN area to form northern populations, presumably due to tolerance of lower temperatures.
Our genome-wide analysis detected genetic introgression of the northern into the southern population ( Figure 8). Gene flow, opposite the Kuroshio Current, might be due to artificial introductions. It is known that pearl oysters were brought from the Japanese mainland to Amami Island (NAN_18) and Okinawa Island (NAN_19) for aquaculture purposes in the past. Although pearl farming using P. fucata is not continuing at present, hybridization has occurred in these areas and genetic attributes have been retained in the NAN population.
In terms of conservation of P. fucata, an important aquaculture species for pearl farming, our results provide essential information about population structure. In order to conserve genetic diversity of the species, artificial admixture for pearl farming should be planned carefully based on natural population structure. In particular, animals from the southern population should not be introduced for pearl farming in the Japanese mainland, or vice versa, to avoid losing their genetic uniqueness. As the northern and southern populations may be adapted different environmental conditions, it is beneficial for sustainable pearl farming to conserve genetic diversity that might be able to overcome climatic changes in the future.

| CON CLUS IONS
Using genome-wide genotyping data, we determined the population structure of Pinctada fucata in the Indo-Pacific area. Population isolation by geographic barriers is evident in the Indian Ocean (MMR) and a local lagoon in Kamikoshiki Island (JPN_17), whereas environmental gradients can shape large-scale population differentiation in the western Pacific (northern and southern populations). Genomic loci that may be associated with adaptation to the local environment were identified. The present study highlights demographic history affected by environmental changes that explain the extant marine population in the western Pacific. In addition, the population structure demonstrated by the present study should provide essential information for genetic conservation of the pearl oyster.

ACK N OWLED G EM ENTS
This work was supported by grants from the Project to Advance Institutional Bio-oriented Technology Research, NARO (special project on advanced research and development for next-generation technology), and by internal funds from the Okinawa Institute of Science and Technology (OIST). We thank Dr. Steven D. Aird for editing the article.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
SNP dataset is available on DRYAD (https ://doi.org/10.5061/dryad. qnk98 sfbv) (Takeuchi, 2019).