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

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.

Neutrality tests and demographic parameters like Tajima’s D (Tajima, 1989), Fu’ Fs (Fu, 1997), Fu’s F* and Li’s D* (Fu & Li, 1993), R2 (Ramos-Onsins & Rozas, 2002), mismatch distribution, and tau (τ), were calculated in DnaSP v6 (Rozas et al., 2017) using the clusters reported by BAPS. The significance of these tests was assessed using the coalescence simulations implemented in DnaSP, with 1000 replicates and the standard neutral model (SNM) that assumes a large constant population size (Hudson, 1990). Genetic diversity indices like the number of haplotypes (H), haplotype and nuclear diversity (h, π) (Nei, 1987), and the number of polymorphic sites (S), were also calculated in DnaSP.

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⁻⁸ 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⁷ 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 marginal seas. Also, most of the nodes included haplotypes with a single individual. BAPS suggested the presence of three clusters (log(ml) = −8545.8301), herein Lineages A, B, and C; in which the frequency of lineage B and C were similar in all marginal seas except in the South China Sea, and lineage A was more abundant in the South China Sea compared with the other locations (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.