The role of the isolation of the marginal seas during the Pleistocene in the genetic structure of black sea bream Acanthopagrus schlegelii (Bleeker, 1854) in the coastal waters of Japan

The black sea bream Acanthopagrus schlegelii (Bleeker, 1854) is a commercially important species in Japanese waters. Assessing its population structure is essential to ensure its sustainability. In the Northwestern Pacific, historical glacial and interglacial periods during the Pleistocene have shaped the population structure of many coastal marine fishes. However, whether these events affected the population of black sea bream remains unknown. To test this hypothesis and to assess the population structure of black sea bream, we used 1,046 sequences of the mitochondrial control region from individuals collected throughout almost the entire Japanese coastal waters and combined them with 118 sequences from populations distributed in other marginal seas of the Northwestern Pacific Ocean. As in other coastal marine fish with similar distribution, we also found evidence that the glacial refugia on the marginal seas prompted the formation of three lineages in black sea bream. These lineages present signatures of population growth that coincided with the interglacial periods of the Pleistocene. While the origin of Lineages B and C remains unclear, the higher relative frequency of Lineage A in the southernmost location suggests its origin in the South China Sea. The non-significant pairwise ΦST and AMOVA of Japanese populations and the presence of these three lineages mixed in Japanese waters; strongly suggest that these lineages are homogenized in both the Sea of Japan and the Pacific Ocean. Our results indicate that the black sea bream should be managed as a single stock in Japan until the strength of connectivity in contemporary populations is further addressed using non-coding nuclear markers.


INTRODUCTION
Collection sites of A. schlegelii specimens used in this study. Fully details of populations are described in Table 1.

Sample collection
A total of 1,046 individuals of black sea bream were collected at 13 different locations throughout almost the entire coastal waters of Japan. We also collected 30 individuals from the Miaoli County coast in Taiwan (Fig. 1, Table 1). For each individual, a piece of the pectoral or caudal fin was cut and preserved in 99% ethanol for DNA isolation.

DNA isolation and sequencing
Genomic DNA was isolated using the TNES-urea buffer (Asahida et al., 1996) followed by standard phenol-chloroform protocol. Fragments longer than 686 base-pairs of the mitochondrial DNA (mtDNA) control region was amplified by polymerase chain reaction (PCR) using the forward primer HDDloopF-54 (5 -CCTATTGCTCAGA GAAAAGGGATT-3 ) and the reverse primer HDDloopR-43 (5 -CCTGAAGTAACCAGATG-3 ), designed by Shi et al. (2015). For each individual, the PCR reaction was carried out in a total volume of 10 µL containing 6.75 µL of ultra-pure water, 1.0 µL of 10X Taq Buffer, 0.8 µL of dNTP, 0.05 µL of TaKaRa ExTaq DNA polymerase (TaKaRa, Shiga, Japan), 0.2 µL of each primer and 1 µL of template DNA. PCR was performed in a Master cycler Gradient 96-Well system (Eppendorf, Hamburg, Germany) with initial denaturation at 95 • C for 5 min followed by 32 cycles of 94 • C for 30 s, 56 • C for 30 s and 72 • C for 1 min 30 s; and a final extension at 72 • C for 10 min. PCR products were treated with ExoSAP-IT (Affymetrix/USB Corporation, Cleveland, OH) and then sequenced on a Genetic Analyzer (ABI3130x1, Applied Biosystems) using the BigDye v3.1 Terminator Sequencing Kit (Applied Biosystems) and the forward primer.

Data analyses
We pooled our sequences with additional fifty-nine sequences of the same DNA region reported in Shi et al. (2015) that correspond to three different locations in the coast of China (the East and South China Sea) (GenBank accession numbers: KJ586516-KJ586574). Sequences were aligned in Clustal X (Thompson, Higgins & Gibson, 1994), and collapsed into haplotypes after trimmed. Pairwise genetic divergence was calculated with the fixation index ST implemented in Arlequin v3.5 (Excoffier & Lischer, 2010) using the Kimura 2P substitution model (Kimura, 1980) and with the significance tested by 10,000 permutations and the Benjamini and Yukutieli's FDR approach (Benjamini & Yekutieli, 2001) for non-independent test. Additionally, hierarchical population structure was evaluated using the analyses of molecular variance (AMOVA) implemented in Arlequin. For this purpose, populations were grouped considering: (1) populations in Japanese waters, (2) populations in the marginal seas of the Sea of Japan, East China Sea, and the South China Sea; and the Pacific Ocean side, and (3) mitochondrial clusters (see BAPS methods and results).
To display the relationship between mitochondrial haplotypes, we constructed a TCS parsimony network using the program PopART ver. 1.7 (Leigh & Bryant, 2015). We also construct a maximum likelihood phylogenetic tree using the IQ-TREE software (Nguyen et al., 2015) with the ultrafast bootstrap approximation and the best model selected by ModelFinder (Kalyaanamoorthy et al., 2017). As the outgroup, we used Acanthopagrus berda (Genbank accession number: AM992253.1). Besides, we assessed the number of clusters using a bayesian approach implemented in BAPS v5.2 (Corrander, Marttinen & Mäntyniemi, 2006) with the linked loci model (Corander & Tang, 2007) and ten independent runs for K = 1:7.
The time of demographic expansion (T) was calculated from the relationship T = τ /2u, where u = 2µk and with µbeing the mutation rate and k the alignment length (686 bp here). We consider an average generation time of four years (Gonzalez, Umino & Nagasawa, 2008), and a mutation rate of 3.6 × 10 −8 per site per year (3.6%/Myr) reported for the mtDNA control region of teleost (Donaldson & Wilson Jr, 1999).
Past dynamics of the female effective population size were observed through Bayesian Skyline Plots (BSPs) calculated in BEAST 2.4.8 (Bouckaert et al., 2014). We used the strict molecular clock and the HKY substitution model. BSPs were run with three independent runs of 5 × 10 7 generations sampled every 5,000 iterations, and with the first 25% discarded as burn-in. All runs yielded an effective sample size (ESS) of more than 200 for the parameter of interest after burn-in. The skyline plots were generated in Tracer v.1.7.1 (http://tree.bio.ed.ac.uk/software/tracer/).

Genetic structure
Our alignment matrix of 686 base pairs yielded 597 different haplotypes with 205 variable sites. The pairwise ST between the 13 Japanese populations were only significant between Nagasaki (NG) and Hiroshima (HH) populations ( ST = 0.023, P = 0.006). When Japanese populations are compared with Chinese and Taiwanese population, the Zhenjian population (DJ) from the Northern East China Sea is only significantly different from the Hiroshima (HH) population ( ST = 0.041, P = 0.005). From the remaining 45 pairwise comparisons, 31 were significant and their ST ranged from 0.036-0.121. While some of the pairwise ST comparisons were significant, all of the values were still very low (Table 2).
Hierarchical AMOVA supported the genetic structure among groups based on marginal seas ( CT=0.00761, P = 0.0039); however, 98.5% of the variation was observed within populations. When only Japanese groups are considered, AMOVA indicates the absence of population genetic structure among Japanese populations ( ST = 0.002, P > 0.05) with 99.73% of the variation explained within populations. Similarly, the FST values were significant but very low (Table 3).
The haplotype network did not display signs of geographical structuring, supporting the high connectivity and low FST statistics between geographic populations and among the    Fig. 2A). However, our phylogenetic tree did not support these lineages nor showed clustering of haplotypes (Fig. S1). In the TCS network, the three lineages provided by BAPS displayed many start-like shapes with a center haplotype which suggested a population expansion. These shapes were more clear in Lineages B and C, and some haplotypes from Lineage A showed a closer relationship with haplotypes from Lineage B. Pairwise FST of sequences that included these lineages was high and strongly significant ( FST of 0.211, 0.411, and 0.540 with P = 0 for all comparisons) and the AMOVA among them had accumulated a variation of 41.29% with an FST of 0.412 at P = 0 (Table 3). The TCS network of haplotypes and the distribution of these lineages are shown in Fig. 2B.

Demographic history
The genetic diversity indices of each lineage are shown in Table 4. The haplotype diversity for each lineage was high and similar (h ranged from 0.957-0.995), and the nucleotide diversity was overall low but higher in Lineage A (π = 0.0097) than in the Lineage B (π = 0.0048) and C (π = 0.0043). Supporting the shape for each lineage within the haplotype network, we observed strongly significant values for all the neutrality test and R2 statistics suggesting events of past demographic expansion ( Table 4). The mismatch distribution also fits the model of demographic growth rather than the model of constant population size (Fig. 3A). The expansion time calculated based on tau (τ ) ( Table 4) were similar to those calculated from the Bayesian skyline plot (Fig. 3B). Overall, the expansion time for Lineage A seems to have occurred between 300-340kya, and for Lineage B and C between 120-140 kya.

DISCUSSION
Population genetic structure of marine fishes are assumed to be very low or inexistent because of absence of physical barriers and wide-range dispersal of larvae by the ocean currents (Ward, Woodwark & Skibinski, 1994). In the Northwestern Pacific, however, mitochondrial sequences of marine fishes shows signatures of moderate and strong genetic differentiation depending on the life-history trait of species (Hirase & Ikeda, 2014;Liu et al., 2007;Shen et al., 2011;Shui et al., 2009;Song et al., 2017;Wang et al., 2008;Yu, Kai & Kim, 2016). These signatures are largely attributed to the isolations of marginal seas which produced topography variations, isolation, sea level fluctuation, and passive dispersal of ocean currents. These isolations have occurred during the Pleistocene epoch (Song et al., 2017;Voris, 2000;Wang, 1999;Yan et al., 2015).
For black sea bream, BAPS bayesian clustering detected the presence of three different lineages which explained 41.29% of the variation in the AMOVA analyses. The presence of these lineages within each of our geographical locations ( Fig. 2A) also produced nonsignificant pairwise FST and AMOVA (Tables 2 and 3) of black sea breams inhabiting Japanese waters. In the Northwestern Pacific, glacial periods during the Pleistocene also played a key role in the formation of three lineages in the coastal scaled sardine Sardinella zunasi (Wang et al., 2008), based the mtDNA control region; the flathead mullet Mugil cephalus (Shen et al., 2011), based on the mtDNA COI; and for the manila clam Ruditapes philippinarum (Mao et al., 2011), based on the mtDNA COI. These results support that Pleistocene glaciation might also have shaped the divergence of black sea bream populations into three lineages and that interglacial periods homogenized these lineages around its current distributions.
The South China Sea appeared composed mostly by individuals of Lineage A, suggesting the origin of this lineage in this northernmost location of our collections, very likely due to the enclosure of this marginal sea during glacial periods of the Pleistocene (Wang, 1999). The relative frequency of Lineages B and C, however, are similar in the Pacific   Obs. Exp.

Const.
Obs. Exp. Ocean side, the Sea of Japan, and the East China Sea, making impossible to assign the correct origin of them and suggest that biological characteristics of this fish and physical factors of these waters might have homogenized them during Pleistocene interglaciations (Kimura, 1996;Kitamura et al., 2001). Supporting the formation of two lineages driven by Pleistocene glaciations in similar regions of Lineage B and C here, Bae et al. (2020) also found two lineages of the grass puffer Takifugu niphobles in the Japan Sea, East China Sea, Yellow Sea and the Pacific Ocean region, using the mtDNA COI. The frequency of the COI haplotypes in their samples allowed Bae et al. (2020) to assign these lineages to their correct origin. Thus, inferring the origin of lineage B and C would perhaps require the sequence of other genes with slower mutation rate than the mtDNA control region, as neither bayesian approach ( Fig. 2A) nor maximum likelihood analyses (Fig. S1) provide enough information to assign them to a unique marginal sea. The expansion time calculated the female effective population size of black sea bream lineages coincides with the rise of sea levels during interglacial periods around 300,000 and 150,000 years ago (Table 4 and Fig. 3B), which might have led this species to homogenize more intensively in the Northwestern Pacific.

Lineage
Further studies of populations genetics in black sea bream will benefit from the use of more polymorphic loci from the nuclear genome. These loci will be useful to assess the contemporary population connectivity of this fish. However, studies such as from Bae et al. (2020) suggest that the mito and nuclear patterns (with microsatellite loci) of coastal marine fishes in the Northwestern Pacific are similar.

CONCLUSION
In the present study, we assess the population genetic structure of black sea bream and test whether the glacial-interglacial periods of the Northwestern Pacific Ocean played a role in the divergence and distribution of its populations. To test this hypothesis, we used 1,046 sequences of the mitochondrial control region from black sea breams distributed in almost the entire coast of Japanese waters and combined them with 118 sequences from other marginal regions. We observed the presence of three different lineages of black sea bream with signatures of population expansion during interglacial periods of the Pleistocene. The sequences of these lineages are carried by black sea bream populations inhabiting Japanese waters suggesting that, although separated in the past, they have already homogenized. Our results highlight the role of historical climate fluctuations into the genetic structure of coastal marine fishes of the Northwestern Pacific and provide useful information for the fisheries management of black sea bream. Based on our results, we recommend managing black sea bream as a single stock in Japan until more studies using more polymorphic nuclear DNA markers clarify the absence or presence of any contemporary barrier.