Population genetic structure of Marbled Rockfish, Sebastiscusmarmoratus (Cuvier, 1829), in the northwestern Pacific Ocean

Abstract Sebastiscusmarmoratus is an ovoviviparous fish widely distributed in the northwestern Pacific. To examine the gene flow and test larval dispersal strategy of S.marmoratus in Chinese and Japanese coastal waters, 421 specimens were collected from 22 localities across its natural distribution. A 458 base-pair fragment of the mitochondrial DNA (mtDNA) control region was sequenced to examine genetic diversity and population structure. One-hundred-six variable sites defined 166 haplotypes. The populations of S.marmoratus showed high haplotype diversity with a range from 0.8587 to 0.9996, indicating a high level of intrapopulation genetic diversity. Low non-significant genetic differentiation was estimated among populations except those of Hyogo, Behai, and Niiigata, which showed significant genetic differences from the other populations. The demographic history examined by neutrality tests, mismatch distribution analysis, and Bayesian skyline analysis suggested a sudden population expansion dating to the late Pleistocene. Recent population expansion in the last glacial period, wide dispersal of larvae by coastal currents, and the homogeneity of the environment may have important influences on the population genetic pattern. Knowledge of genetic diversity and genetic structure will be crucial to establish appropriate fishery management of S.marmoratus.


Introduction
The Marbled Rockfish, Sebastiscus marmoratus (Cuvier, 1829), valued for its high nutritional value and palatability (Zhu et al. 2011), is widely distributed in the coastal areas of the Northern Pacific Ocean, especially in China, Korea, and Japan (Higuchi and Kato 2002;Jin 2006;Nakabo 2013). In recent years, the number of S. marmoratus dramatically decreased due to overfishing, pollution, and habitat destruction , potentially influencing its genetic diversity and population structure.
The mitochondrial DNA (mtDNA) control region has been shown to be particularly effective in detecting population genetic structure and diversity, owing to its high polymorphism, maternal inheritance, high mutation rate, and nonrecombinant DNA (Bowen and Grant 1997;Whitehead et al. 2003;Dowling et al. 2008). Zhang et al. (2016) used mtDNA variation to determine that currents and larva dispersion with drifting seaweed influenced the phylogeographic pattern and genetic homogeneity of Sebastiscus schlegelii (Hilgendorf, 1880).
While S. marmoratus has been widely studied with respect to taxonomy (Hansen and Karlsbakk 2018), genetics (Deng et al. 2015;Cai et al. 2017;Xu et al. 2017), culture (Yin and Qian 2017;Watanabe et al. 2018), and breeding habitat (Chen et al. 2016;Guo et al. 2016), the genetic structure of its populations along the Chinese and Japanese coasts is not known. In view of the many genetic studies based on mtDNA Liu et al. 2018), we selected mtDNA markers to analyse the population genetics of S. marmoratus.
The goals of this study were to estimate genetic diversity, to characterize genetic structure, and to reconstruct the evolutionary relationships of S. marmoratus in its distribution range. Failure to characterize population units can lead to overfishing and severe decline (Waples 1998). Elucidation of S. marmoratus population genetic structure is crucial for its conservation management. The wide distribution of S. marmoratus throughout the NW Pacific, along with its short-distance migration life history, makes it an ideal candidate for investigating how the complex geological history of the Northwestern (NW) Pacific shapes intra-species diversity of the fish fauna.

Sample collection
From June 2009 to August 2015, we collected 421 wild S. marmoratus from 22 locations in coastal China and Japan, 10-24 specimens per site (Fig. 1, Table 1). Muscle tissue samples were preserved in 95% ethanol for subsequent DNA extraction.

DNA extraction, amplification, and sequencing
Genomic DNA was extracted from muscle tissue by proteinase K digestion followed by a standard phenol-chloroform technique. Fragments of the mtDNA control region were amplified with primers referenced from Han et al. (2008): DL-S (5'-CCC ACC ACT AAC TCC CAA AGC-3'), DL-R (5'-CTG GAA AGA ACG CCC GGC ATG-3').
Polymerase chain reactions (PCR) were carried out in 25 μL of reaction mixture containing 10-100 ng template DNA, 0.1 μL (5 U/μL) Taq DNA polymerase (Takara Co., Dalian, China), 1.5 μL (10 pmol/μL) of each forward and reverse primer, and 2 μL (200 μmol/L) deoxy-ribonucleoside triphosphate (dNTP). The PCR amplification was conducted in a Biometra thermal cycler under the following conditions: 2 min initial denaturation at 95 °C; 40 cycles of 60 s at 94 °C for denaturation, 45 s at 52 °C for annealing, and 60 s at 72 °C for extension; and a final extension at 72 °C for 8 min. The PCR product was purified with a Gel Extraction Mini Kit (Watson BioTechnologies Inc., Shanghai, China). The purified product was used as the template DNA for cycle sequencing reactions performed using BigDye Terminator Cycle Sequencing Kit (v. 2.0, PE Biosystems, Foster City, CA, USA), and bi-directional sequencing was conducted on an Applied Biosystems Instrument Prism 3730 automatic sequencer (Sunny Biotechnology Co. Ltd, Shanghai, China) with both forward and reverse primers. The primers used for sequencing were the same as those used for PCR amplification.

Data analysis
All sequences were edited and aligned manually by DNAStar software (DNAStar Inc., Madison, WI, USA) using default settings and were manually corrected. The genetic diversity indices of S. marmoratus, including haplotype diversity (h), nucleotide diversity (π), mean number of pairwise differences (k), and number of polymorphic sites were calculated by ARLEQUIN v. 3.5 (Excoffier and Lischer 2010).
Nucleotide sequence evolution models were evaluated using likelihood-ratio tests implemented by Modeltest v. 3.06 (Posada and Crandall 1998). The neighbor-joining (NJ) tree of the haplotypes was rooted with the out-group Sebastes schlegelii (Zhang et al. 2016) using MEGA v.5.0 and evaluated with 1000 bootstrap replicates (Tamura et al. 2011) to reconstruct phylogenies of haplotypes. Among-site heterogeneity was corrected with the shape parameter of gamma distribution (γ = 0.697). The GenBank accession number of S. schlegelii is JX241455. Pairwise genetic divergence among populations were tested by the fixation index Fst (Excoffier et al. 1992), and the significance of the Fst was evaluated by 10,000 permutations for each pairwise comparison in ARLEQUIN v. 3.5 (Excoffier and Lischer 2010). The P values were adjusted by Bonferroni correction (Rice 1989). The Fst and P-value heatmaps with dendrograms were created with the R project 3.5.1 (www.r-project.org).
Characterization of population subdivisions and population structure were conducted using a hierarchical analysis of molecular variance (AMOVA) of different gene pools (Excoffier et al. 1992). In addition to separate total population into the same gene pool analysis, AMOVA analyses were carried out on populations from the Chinese and Japanese coast; North China coast, South China and South Japan coast; North Yellow Sea, South Yellow Sea, East China Sea, South China Sea, and the Japanese coast. The haplotypes were assessed with 1000 permutations in AMOVA.
The Tajima D and Fu's F S tests were examined for neutrality (Tajima 1989;Fu 1997). Historical demographic expansions were also investigated by examination of the frequency distribution of pairwise differences between sequences (mismatch distribution) based on three parameters: θ 0 , θ 1 (θ before and after the population growth), and τ (time since expansion expressed in unit of mutational time) (Rogers and Harpending 1992). Historical pure demographic and range expansions were further investigated by the mismatch distributions using ARLEQUIN v. 3.5 (Excoffier and Lischer 2010). Unimodal distribution patterns reflect recent demographic or range expansion with a high level of migration between neighboring demes, while multimodal patterns indicate relatively stationary populations (Rogers and Harpending 1992;Ray et al. 2003). The sum of square deviations (SSD) and Harpending's raggedness index (HRI) were used to test goodness-of-fit of the observed unimodal mismatch distribution to that expected under the sudden expansion model. The time since population expansion was estimated using the equation τ = 2μt, where μ is the mutation rate for the entire DNA sequence under study, and t is the time since expansion. We used the sequence divergence rate of 5%-10%/MY (Brunner et al. 2001) for the control region sequences.
Bayesian skyline analyses, implemented in BEAST v. 1.7.4 (Drummond and Rambaut 2007), were performed to estimate changes in effective population size through time, which can indicate past demographic changes by comparison with current patterns of genetic diversity within a population (Drummond et al. 2005). To check for convergence, we executed multiple independent runs for 300,000,000 iterations under an HKY+I+G nucleotide substitution model and a strict molecular clock, with individual parameters estimated from the data with a piecewise constant skyline model of 10 groups. Genealogies and model parameters were sampled every 10,000 generations with the first 10% discarded as burn-in. Trace plots were inspected to assess mixing, convergence, and stationary distribution of the MCMC process in Tracer v. 1.5 (Rambaut and Drummond 2009). The effective population sizes were checked and confirmed as >200 for each parameter in order to avoid autocorrelation of parameter sampling.

Genetic diversity
A 458 bp segment of the mtDNA control region was amplified, and 106 polymorphic sites were detected, including 89 transitions and 17 transversions. A total of 166 haplotypes were identified based on the sequence variation in 421 individuals from 22 locations. Among these, 84 haplotypes were shared. The most common haplotypes, Hap4 and Hap5, were both shared by 40 individuals. Haplotype sequences were deposited in GenBank under accession numbers KY703229-KY703394. The estimated nucleotide diversity (π) and haplotype diversity (h) for the locations are shown in Table 2. The mean value of π was 0.0220±0.0110 with highest in Kochi (0.0250±0.0132) and the lowest in Hyogo (0.0098±0.0056). The mean value of h was 0.9560±0.0035 with the highest in Awaji (0.9996±0.0243) and the lowest in Hyogo (0.8587±0.0337).

Population structure
An unrooted phylogenetic tree was reconstructed by neighbor-joining analysis using 166 haplotypes with the best nucleotide substitution mode (HKY+I+G) rooted with the outgroup S. schlegelii. There were no significant genealogical branches or clusters corresponding to sampling localities (Fig. 2). The relationships among 166 haplotypes were represented on the minimum spanning tree (MST). The minimum spanning network was generally star-like with several common and ancestral haplotypes shared by most populations (Fig. 3). The MST was connected and indicated recent population expansion. The Hyogo population showed an obvious haplotype branch and others exhibited no unique haplotype corresponding to geographic populations.    (Table 4). In general, most of the pairwise Fst values among populations showed non-significant differences after sequential Bonferroni correction. However, significant genetic differences were obtained among Hyogo, Behai, and Niigata populations and between Hyogo, Behai, and Niigata and the other populations. The largest difference was seen between Niigata and Hyogo (Fst = 0.677, P < 0.05). Some pairwise Fst estimates were negative, indicating that within-population variation was greater than that between populations. The global AMOVA showed about 13.59% of the genetic variation to be among populations. Other grouping methods also indicated most genetic variation to be within populations, rather than among groups and populations.

Population historical demography
Tajima's D (D = -1.027; P > 0.05) and Fs test (Fs = -23.917; P < 0.01) results were negative, indicating departure from selective neutrality (Table 5). Non-significant and low values of SSD and HRI were found for each population and for the overall population, suggesting a sudden expansion model. The sudden expansion model of mismatch distribution was unimodal and a valid goodness-of-fit was observed between observed and expected distributions (Fig. 4), indicating strong demographic expansion. The τ value, which reflects the location of the mismatch distribution crest, provides a rough estimate of the time of initiation of rapid population expansion. According to t = τ/2μ, based on τ values (τ = 12.305) and divergence rate of 5-10% per site per Myr (Brunner et al. 2001), the pure population expansion occurred 268,000-448,000 years ago. The ratio of estimated effective female population size after expansion to that before expansion (θ 1 /θ 0 ) was 56.84 (Table 4). Results indicated that S. marmoratus underwent     colonization and recent population expansion events along the Chinese and Japanese coasts during the Pleistocene. The Bayesian skyline plot (Fig. 6) indicated an historic occurrence of a continual gradual increase in the effective size of S. marmoratus populations, dating to about 430,000 years BP at the end of the Pleistocene. These results are consistent with a process of historical expansion of S. marmoratus populations, as indicated by negative D and Fs values and mismatch distributions.

Discussion
Inbreeding depression and other genetic problems impacted by human behavior can be monitored by assessing genetic diversity under natural conditions (Ryman 1991;Smith et al. 1991). The adaption of marine organisms to their surroundings and their evolutionary potential can be affected by genetic diversity (Templeton 2010). We found that, despite high haplotype diversity (h = 0.9560±0.0035) of S. marmoratus in the northwestern Pacific Ocean, its nucleotide diversity was low (π = 0.0220±0.0110). The high mutation rate of the D-loop region may be a factor in this phenomenon (Wan et al. 2004). Haplotype diversity with low nucleotide diversity may indicate population reduction or the existence of a genetic bottleneck and may result in extinction under environmental pressure. It can be observed in a population experiencing rapid expansion from a low effective population size, assuming adequate time for the increase in haplotypes through mutation but inadequate time for accumulation of large sequence differences (Lowe et al. 2004). The retention of new mutations in the population can be enhanced by rapid population growth (Grant and Bowen 1998). The phenomenon of high haplotype diversity and low nucleotide diversity has been reported in organisms such as Glyptocephalus stelleri (Xiao et al. 2010), Trachurus japonicus (Song et al. 2013) and Circus spilonotus (Nagai et al. 2018) that have undergone a rapid severe population reduction.
Compared with anadromous and freshwater fishes, marine species are generally expected to show a low degree of genetic differences among geographic regions owing to their high dispersal potential through planktonic drifting of eggs, larvae, or adults and the absence of physical barriers (Palumbi 1994;Hewitt 2000;Liu et al. 2018). Our AMOVA results and the neighbor-joining analysis did not show significant genetic structure among geographic populations. Ecological characteristics and marine currents may play important roles in shaping the contemporary phylogeographic pattern of marine fishes. For example, the rockfish S. schlegelii is typical of fish that congregate in drifting seaweed during early development (Ikehara 1977;Safran and Omori 1990;Zhang et al. 2016). Sebastiscus marmoratus, a species of rockfish with life history similar to S. schlegelii, is believed to exhibit the same behavior, dispersing with drifting seaweed during November and December and in the following year from February to April (Mitchell and Hunter 1970;Wu et al. 1999). The Kuroshio Current is one of the strongest currents in the world and can accelerate gene flow from the southern East China Sea to the coastal waters of Japan (Liu et al. 2007). Inflows from the Yellow Sea enter the Bohai Sea along the west coast of Korea via the Yellow Sea warm current and the China Coastal Current (Jin et al. 2010). Waters also exchange between the warm Yellow Sea and Kuroshio currents. These strong currents might transport S. marmoratus larvae via drifting seaweed and promote exchange throughout its range.
Recent research reveals that the currently most common unintentional pathway for the transport of marine organisms is the ballast water of commercial vessels (Ruiz and Hines 1997;Wonham et al. 2000). As human activity becomes more frequent and extensive, trade between countries is strengthened, and commercial vessels traverse large area. Ballast water is usually taken from the harbor in one port and subsequently discharged in another port (Carlton and Geller 1993). Diverse organisms including protist, diatoms, invertebrate larvae, and copepods are collected and survive the voyage to the next port (Carlton 1985;Smith et al. 1999). Corrosion or other damage to protective grates or ballasting of water by gravitation, may provide access to the ballast for larger organisms such as post-larval fish (Springer and Gomom 1975;Wonham et al. 2000). Hansen and Karlsbakk (2018) reported S. marmoratus caught by bait in Norwegian waters and thus shown to be actively foraging, a strong indication that S. marmoratus may thrive in unfamiliar conditions.
We may conclude that transportation via ballast water may be a source of genetic homogeneity of S. marmoratus. It has been transported among ports and wharves along the NW Pacific Ocean and into rocky coastal areas with the release of ballast water, where it can easily survive (Hansen and Karlsbakk 2018). With adaptation to the same or similar environment, a number of invasive populations of S. marmoratus have been reported (Wonham et al. 2000). It is also possible that environmental factors such as salinity and temperature have brought about adaptive evolution of S. marmoratus. However, this is undetectable by the molecular markers we used in this study. It has been suggested that a genotyping-by-sequencing technique could reveal the occurrence of local adaptation .
Using genetics to understand biogeography is important to determine patterns influencing distribution of geographically distant populations. Genetic diversity, genetic distribution patterns, and effective population size were also influenced by paleogeological changes and fluctuations as well as life history and marine environment factors (Hewitt 1996). Following sea level falls in the glacial period (Wang 1999;Lambeck et al. 2002), S. marmoratus may have experienced a population contraction with the loss of genes of those dying out and the majority of survivors migrating to more suitable environments. A single branch from the star-like network representing the Hyogo population suggests a likely founder effect (Ramachandran et al. 2005).
Significant genetic differences were revealed between Hyogo and other populations based on the star-like network tree and Fst value analysis, which suggests that the deep and semi-open area of inland waters might have an impact on the geographic isolation. Genetic differences among Hyogo, Behai, and Niigata populations and between Hyogo, Behai, and Niigata and the other populations were primarily significant, and possibly relate to convergence evolution (Wilkens and Strecker 2003) and the formation of a refuge (Consuegra et al. 2002). In addition, mutation-drift disequilibrium may exist among these populations, which are in an unstable state of genetic mutation-drift (Lacy 1987). Further molecular marker studies are required to evaluate this proposition.

Conclusions
Climate fluctuations caused by glacial-interglacial alternation, early life-history, and ecological characteristics, combined with transport via ballast water may play important roles in the extensive gene flow among populations and the current genetic distribution pattern of S. marmoratus. Information provided by the current study will facilitate its comprehensive management. Future studies should be based on informative nuclear markers to provide additional information on genetic structure and differentiation of populations of S. marmoratus.