Genetic population structure and demographic history of the endemic fish Paralichthys olivaceus of the Northwest Pacific Ocean

Abstract The Northwest Pacific marginal waters comprising the South China Sea, East China Sea, Yellow Sea, and the Sea of Japan have unique geomorphic features. The Japanese flounder Paralichthys olivaceus, which is endemic to the Northwest Pacific, has high nutritional, economic, and ecological value. To allow the examination of the demographic history and population structure of the most common P. olivaceus species range over the five marginal seas (East China Sea, Yellow Sea, Bohai Sea, Northwest Pacific Ocean, and the Sea of Japan), the mitochondrial DNA control region of 91 individuals from six populations in China was sequenced. These sequences were combined with 233 sequences from four populations distributed in the Northwestern Pacific Ocean for analysis. Higher levels of nucleotide diversity (0.032 ± 0.016) and haplotype diversity (0.996 ± 0.001) were observed. The peripheral Fuqing population in the East China Sea had the relatively lowest genetic diversity and highest differentiation. Furthermore, when the results of the isolation by distance test, spatial analysis of molecular variation and geographic barrier analysis are also considered, there is a clear need to prioritize resource conservation and enhancement measures in this area. The phylogenetic trees, structure assignment test, and haplotypes network revealed no significant differences in the genealogical structure among ten populations. Mismatch distribution analysis, Bayesian skyline plots, and neutrality tests suggested that P. olivaceus experienced population expansion during the Pleistocene. Ocean currents and climate change play important roles in shaping the geographical distribution and genetic population structure of P. olivaceus.


| INTRODUC TI ON
Marginal water ecosystems in the Northwest Pacific are exposed to intense anthropogenic stresses, such as pollution and overfishing (Halpern et al., 2008;Ni et al., 2014;Yamashita et al., 2021;Zhang et al., 2020), and are affected by the complex ocean current patterns in the region. During the Last Glacial Maximum, the sea level in the Northwest Pacific decreased by 120-140 m. Marginal seas in the region were separated by land, forming many independent small sea areas with fragmented habitats (Lambeck et al., 2002;Voris, 2000;Zhao et al., 2022). The processes of sea-level rise and fall have greatly impacted the species formation and genetic structure of organisms in the marginal waters of the Northwest Pacific (Liu et al., 2007;Ni et al., 2014). Many researchers have used molecular biology methods to evaluate the effects of regression, sea immersion, and sea habitat fragmentation on organisms in this region (Bae et al., 2020;Tang et al., 2010). The dispersal ability in the early life cycle of marine organisms can determine the genetic structure of marine biological populations to some extent (Caley et al., 1996). Through the action of ocean currents, groups of marine organisms with long pelagic larval durations may have genetic connectivity with others in different geographic locations. This connectivity may contribute to the lack of obvious genetic structures among many marine organisms with large geographical ranges (Blanco-Bercial & Bucklin, 2016;Grant & Bowen, 1998;Hewitt, 2000).
The Japanese flounder Paralichthys olivaceus (Temminck & Schlegel, 1846) belongs to the order Pleuronectiformes and the family Paralichthyidae. Paralichthys olivaceus lives in bottom waters with warm temperatures and has important economic and ecological value (Hamidoghli et al., 2020;Kim et al., 2010;Sekino, 2006;Shigenobu et al., 2013). The Japanese flounder is the only member of Paralichthys along the coast of Asia, where it mainly feeds on crustaceans and small fish. It is widely distributed in the coastal areas of China, the Korean Peninsula, and Japan (Fujii & Nishida, 1997;Sekino et al., 2002;Zhang et al., 2019). The number of eggs laid by a 480 mm sexually mature female P. olivaceus is approximately 2 × 10 5 , whereas 600 mm sexually mature females lay approximately 4 × 10 5 (Ochiai & Tanaka, 1986). The eggs and larvae can float on ocean waves for 25-50 days (Ochiai & Tanaka, 1986;Xu et al., 2012); the long floating period of the fish larvae and juveniles enables P. olivaceus to diffuse over a long distance with the ocean current. In theory, this diffusion strengthens connectivity between various geographical groups. Some of the characteristics of P. olivaceus, such as its large spawning numbers, large eggs, long floating period of larvae, and spawning migration (Kim et al., 2010;Sekino et al., 2002), are suitable for studying the impact of recontact of marine fish on the genetic structure of the population after habitat fragmentation in the Quaternary ice age.
Many studies of the mitochondrial DNA (mtDNA) of P. olivaceus have been conducted in recent years (Ando et al., 2016;Sekino et al., 2002;Yamashita et al., 2021) to clarify the genetic population structure of this fish. A study of the populations of P. olivaceus in the Sea of Japan (Fujii & Nishida, 1997) revealed high variability in the mtDNA control region. Analysis of the sequences of the mtDNA ND2 and ND5 genes indicated the presence of separate stocks among the northern and southern waters of the Pacific coast of Tohoku, Japan (Shigenobu et al., 2013). The highest rates of base substitutions and insertion/deletion events in mtDNA have been detected in the first half of the control region (adjacent to the proline transfer RNA gene) (Fujii & Nishida, 1997;Saccone et al., 1987). Therefore, in this study, we selected the first half of the mtDNA control region as the amplification region. The most common P. olivaceus species range was used as the sampling area, as dense sampling across species ranges is necessary to understand phylogeographic dynamics. We analyzed the population structure and demographic history to describe the geographic distribution patterns of P. olivaceus in all marginal seas of the Northwestern Pacific.

| Sample collection and DNA isolation
Ninety-one P. olivaceus specimens were collected from six geographical locations across Chinese coastal waters from 2004 to 2008 ( Figure 1; Table 1). All individuals were identified based on their morphological characteristics (Nakabō, 2002;Nakabo & Doiuchi, 2013). A piece of muscle tissue was collected from each individual and stored at −20°C. Genomic DNA was extracted using standard phenol-chloroform extraction protocols with proteinase K treatment (Sambrook et al., 1989).

| DNA amplification and sequencing
The mtDNA control region was amplified using the forward and reverse primers 5′-GTT AGA GCG CCA GTC TTG TA-3′ and 5′-CCT GAA GTA GGA ACC AAA TGC-3′, as described by Sekino et al. (2002). Amplification was performed in a 30 μl reaction volume containing 10 μM of each primer, 10 ng template DNA, 100 μM of each dNTP, 3 μl 10× PCR buffer containing MgCl 2 , and 1.25 U of DNA polymerase (Ex Taq™, Takara Co). The cycling conditions were as follows: preliminary denaturation at 94°C for 4 min, followed by strand denaturation at 94°C (45 s), annealing at 58°C (45 s), and primer extension at 72°C (1 min) for 33 cycles; and final elongation at 72°C (10 min). The PCR products were evaluated using 1% agarose gel electrophoresis. Products with a single bright band were selected and purified. The purified fragments were sequenced using an ABI-3730 Automatic Sequencer (Applied Biosystems).

| Data analysis
We pooled our sequences with an additional 233 sequences reported by Sekino et al. (2002) and Ando et al. (2016), corresponding to four locations (Table 1) on the coast of Japan (GenBank accession numbers: LC129561-LC129707). All sequences were screened, edited, and aligned in Lasergene 7.1 (DNAStar; Burland, 2000). A haplotype network was constructed in HapView, to examine genealogical relationships among haplotypes. All sequences were sorted into complete datasets based on their locations for population genetic analysis. Genetic diversity indices including the number of polymorphic sites (s), numbers of haplotypes (N), mean number of pairwise differences (k), nucleotide diversity (π), and haplotype diversity (h) were calculated for each population using Arlequin v3.5.2.2 (Excoffier & Lischer, 2010).
Genetic divergences between pairs of populations were tested by determining the pairwise fixation index, Φ ST , in Arlequin using the Kimura 2P substitution model; significance was tested by performing 10,000 permutations. When multiple comparisons were performed, p values were adjusted using the sequential Bonferroni procedure (Rice, 1989). Analysis of molecular variance was conducted to study hierarchical population structure and the geographic pattern of population subdivision (Excoffier et al., 1992). In this analysis, populations were grouped consid- Ocean, and Sea of Japan. Meanwhile, spatial analysis of molecular variation (SAMOVA v2.0) was used to detect the geographic genetic structure. The number of groups ranged from 2 to 11, with 1000 permutations, and the maximum Φ CT value (genetic diversity between groups) was calculated to determine the best way to group the population. A genetic distance matrix (Sunde F I G U R E 1 Map of China sea coastal waters showing sampling locations and flow directions of coastal currents and ocean currents (from: Ni et al., 2014). Geographic location of sampling locations of the study area are shown in Table 1 by distance (IBD) patterns in RStudio using the "geosphere" (Hijmans et al., 2017) and "vegan" package (Oksanen et al., 2013) using Spearman's method and 10,000 permutations. Geographic barriers were computed and visualized using Barrier v2.2 based on the Φst pairwise comparison matrix and geographic distance.
This method applies the Monmonier's maximum distance algorithm to identify barriers to gene flow among sites, namely the zones in which differences between pairs of sites are the largest.
Phylogenetic analysis was conducted for haplotypes using the Harvester Web 0.6.94 (Earl & VonHoldt, 2012) based on the ΔK method (Evanno et al., 2005) and the K with the highest likelihood (Pritchard et al., 2000). A structure distribution plot was drawn using Distruct 1.1 (Rosenberg, 2004).
A neutrality test (Excoffier & Lischer, 2010), mismatch distribution analysis (Rogers & Harpending, 1992), and Bayesian skyline plot (Ho & Shapiro, 2011) analysis were used to analyze the demographic history of P. olivaceus. The Tajima's D and Fu's Fs values were calculated to evaluate neutrality using Arlequin (Excoffier & Lischer, 2010), and the significance of the obtained values was tested by generating 1000 random samples under the null hypothesis of selective neutrality. The mismatch distribution analysis was carried out for the inference of historical population dynamics. Bayesian skyline plots were generated using BEAST v2.6.6 (Drummond et al., 2012) based on the GTR + I + G + F model using a strict molecular clock; the mutation rate of the control region was considered as 3.6%/Myr (Donaldson & Wilson Jr, 1999)

| Population genetic structure
Phylogenetic trees of the haplotypes based on ML and BI analyses revealed no clusters corresponding to sampling locations or significant genealogical branches ( Figures S1 and S2). The results from the structure assignment test supported three clusters (K = 3 as the highest value), but the peak of the delta k value was low. No distinct population genetic structure could be obviously distinguished by this assignment ( Figure S3). The haplotype network showed that each geographical group had a mixed distribution pattern with no central haplotype, and the evolutionary relationship showed no genealogical branches (Figure 2; Table S1).

| Population genetic differentiation
Among the 45 pairwise comparisons, 12 were significant after sequential Bonferroni correction. Pairwise Φ ST values between the 10 populations ranged from −0.033 to 0.444 ( Table 3) An examination of the IBD patterns with Φ ST revealed genetic differentiation that increases linearly as a function of geographic distance, as well as significant differences in population differentiation across sampling sites (r = 0.37, p = .021, Figure 3). According to the results of the barrier analysis, the priority barrier was observed between FQ and other populations, implying that genetic differentiation changes abruptly, and these abrupt changes are associated with barriers ( Figure 4). Analysis of molecular variance revealed significant genetic differentiation among the nine populations F I G U R E 2 Phylogenetic relationships of Paralichthys olivaceus haplotypes represented in median-joining haplotype network based on mitochondrial region D-loop sequences. Each circle represents a unique haplotype, and the size of the circle is proportional to its total frequency. Each branch connecting different circles represents a single nucleotide change, and blue dot on the branches represent an additional nucleotide change. Colors denote sample geographic origins as indicated by the legend.
(Φ ST = 0.045, p < .001), among populations within five marginal seas (Φ SC = 0.065, p < .001; Table 4). The best partitioning of the genetic structure was obtained when samples were divided into two groups, based on the largest value of Φ CT (0.169). The first and second groups were composed of fish from the FQ and other populations, respectively. Results of the SAMOVA supported regional genetic subdivision as revealed by the geographic barrier analysis (Table 4).

| Demographic history
A unimodal distribution was observed in the mismatch distribution analysis ( Figure 5). The neutrality test statistics (Tajima's D and Fu's

| Genetic diversity
The mtDNA control regions of P. olivaceus in the Northwest Pacific Ocean showed high haplotype diversity (0.996 ± 0.001) and high nucleotide diversity (0.032 ± 0.016). The higher diversity may be that a TA B L E 3 Pairwise Φ ST estimation (below diagonal) and its associated probability (above diagonal) among ten populations of Paralichthys olivaceus based on mtDNA control region data. large stable population with a long evolutionary history or secondary contacts among differentiated lineages (Grant & Bowen, 1998;Srinivas et al., 2021;Zhong et al., 2020). When combined with the results of Liu et al., 2019 and other analysis results of this study, it is shown that P. olivaceus has a large stable population and carries out panmixia within the distribution range. Li and Wang (1995) found

| Genetic structure and genetic differentiation
In the phylogenetic analysis, the 223 haplotypes were not clustered by location or marginal sea. The haplotype network and structure assignment test also showed no obvious structure or central haplotype. Thus, P. olivaceus does not have an obvious geographical structure. The dispersal of larvae with ocean currents is an important cause of the limited genetic differentiation of marine fishes that have a geographically large distribution range (Strathmann et al., 2002). It is generally thought that the levels of genetic differentiation among marine fish populations are low. Previous studies have reported that extensive gene exchange occurs over a wide geographical range in marine fishes (De Queiroz-Brito et al., 2022;Grant & Bowen, 1998;Han et al., 2008;Song et al., 2010). The Φ ST between the ND and FQ populations (East China Sea) was the highest value among all populations (Table 3)

| Demographic history
The population expansion of P. olivaceus was analyzed using three methods, haplotype nucleotide mismatch analysis, neutrality tests, events affecting a population. Therefore, when the population accumulates large amounts of variation in a short time, Fu's Fs tends to give a large negative value (Su et al., 2001). Combined with haplotype nucleotide mismatch analysis, the results showed a unimodal distribution, indicating that the population recently experienced population expansion events. However, this is different from a demographic scenario assessed by genetic diversity. We hypothesized that it may be because P. olivaceus has a large enough population, and that the population expansion occurs in a small range, which has a small impact on the distribution range of the whole species.
The calculated population expansion time is approximately 310,000-610,000 years ago in the middle and late Pleistocene.
Climatic fluctuations during the Pleistocene led to contraction and expansion of the distribution range of many species, which typically impact genetic diversity (Avise, 2000;Hewitt, 2000;Li, 1997). In the Last Glacial Maximum, the Yellow Sea and the Bohai Sea became land, and the East China Sea was approximately 50% of its current size (Wang, 1990). The survival range of marine fish decreased sharply; therefore, the P. olivaceus population may have been isolated in one or more glacial refugia, and the Northwest Pacific Ocean population may have experienced bottleneck effects.
The marine environment is a very open setting. In marine fish, gene exchange between populations is likely affected by various marine environment-related factors, including geographical distance, ocean circulation, seawater temperature, and salinity (Han et al., 2012). Historical migration trends may be related to major sea-level cycles that occur at intervals of ∼100 kyr over the past ∼800 kyr, with maximum amplitudes of 120-140 m (Imbrie et al., 1992;Lambeck et al., 2002). Many studies have shown a weak genetic differentiation between the geographical populations of surface marine fish that can migrate long distances or swim. This is because of the free diffusion of floating eggs, fish larvae, and juveniles or adults, and few geographical or other obstacles in the open ocean environment, enabling widespread and extensive gene exchange in these marine fish (Canfield et al., 2022;Hewitt, 2000;Palumbi, 1994;Pérez-Rodríguez et al., 2021). However, P. olivaceus is a benthic fish, and its life history indicates no long distance migration habit. Therefore, the reason for panmixia among populations may be related to their early life habits. Active diffusion of fish larvae and juveniles as well as marine environmental factors, such as ocean circulation and climate change in the late Pleistocene, played important roles in shaping the systematic geographical pattern and population genetic structure of P. olivaceus.

| Genetic resource conservation
Based on our results, the FQ population exhibits the lowest genetic diversity and highest differentiation and should be prioritized for protection to avoid further loss of genetic diversity. The mtDNA control region may exhibit hypervariability and homoplasy, which cannot accurately reflect the overall population structure (Takeshima et al., 2005;Verma et al., 2016). To better manage, protect, and utilize P. olivaceus, larger numbers of polymorphic loci from the nuclear genome combined with multiple mtDNA (COI, ND2, ND5) regions should be collected to investigate its population genetics. These loci can be used to assess contemporary population connectivity and genetic differentiation. By combining traditional fishery investigations, physiological and ecological information, and molecular biology methods, the resource status of P. olivaceus can be determined. This study not only has basic theoretical significance for understanding the historical evolution of marine fish populations in the Northwest Pacific but also reveals how fish respond to climate and environmental changes, which is of great practical value for the rational development and formulation of effective species protection strategies. (lead); investigation (lead); project administration (lead); writing -review and editing (supporting).

ACK N OWLED G M ENTS
We thank Dr. Masashi Sekino of Fisheries Resources Institute, Japan Fisheries Research and Education Agency, for kindly providing valuable comments on the earlier draft of the manuscript.
We would like to thank Editage (www.edita ge.com) for their support with language editing. The present study was supported

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available