2 MATERIALS AND METHODS

2.1 Study area

Beibu Gulf (17°–22°N, 105°–110°E) is in the northwestern part of the South China Sea and has a long coastline that belongs to China and Vietnam. It is a semi-closed gulf of approximately 128,000 km2 that ranges from Leizhou Peninsula, Qiongzhou Strait, and Hainan Island to Vietnam, and extends to the Guangxi coast in the north (Ma et al., 2010). The bottom of Beibu Gulf is flat and slopes from the northwest to the southeast, and the water depth in the gulf is typically less than 100 m (average depth, 42 m). The surrounding climate of this gulf is subtropical and monsoonal. Moreover, this gulf contains numerous estuaries from which rivers discharge nutrients. Beibu Gulf is a traditional fishing ground and main source of fishery products for coastal areas because of its high productivity and rich biological diversity, which benefit from its geographic location and climatic conditions (Z. Chen, Xu, Qiu, Lin, & Jia, 2009).

2.2 Sample collection

All fish specimens were collected from fishery surveys carried out by the South China Sea Fisheries Research Institute; these surveys were conducted by the commercial fishing vessel “Beiyu60011” in the northern South China Sea using bottom trawler nets in September 2018. Eleven sampling sites from Beibu Gulf were investigated in the study area, which covered over 45,000 km2 (Fig. 1). In total, 11–16 specimens of each sampling site were used for the DNA analyses after morphological identification. The dorsal fin was removed from each specimen and preserved in absolute ethanol at −20°C.

2.3 DNA extraction, amplification, and sequencing

Genomic DNA was extracted from dorsal fin tissue using the TIANamp Marine Animals DNA Kit (TIANGEN, China) following the manufacturer’s protocol. The concentration for use as a PCR template was adjusted to an A260 of approximately 0.05–0.2. Fragments of the mtDNA COI gene were amplified from total genomic DNA using the polymerase chain reaction. The primers FishF1 and FishR1 (Ward, Zemlak, Innes, Last, & Hebert, 2005) were used for COI PCR amplification. Each 50 µl PCR consisted of 31.25 µl ddH2O, 5 µl PCR buffer, 5 µl CoralLoad concentrate, 4 µl of 25 µM MgCl2, 1 µl of 10 µM dNTPs, 0.5 µl of 25 µM solution of each primer, 2.5 µl DNA template, and 0.25 µl TopTaq DNA polymerase (QIAGEN, Germany). The PCR conditions for amplification were: an initial step of 2 min at 95°C; 35 cycles of 0.5 min at 95°C (denaturation), 0.5 min at 54°C (annealing), and 1 min at 72°C (extension); followed by 7 min at 72°C (final extension) on a 2720 Thermal Cycler (Applied Biosystems, USA). PCR products were visualized on 1.5% agarose gels and the most intense products were selected for sequencing. Products were bidirectionally sequenced on an ABI 3730XL automated sequencer with both forward and reverse primers following manufacturer’s instructions.

2.4 Genetic analysis

The sequences were assembled and edited in Bioedit (Hall, 1999) and aligned using the CLUSTALW multiple algorithm under default options. Ambiguous sequences were trimmed after alignment. The authenticity of our COI sequences was verified by a BLAST search in GenBank (BLASTn, megablast algorithm) and compared with the highest match (99%–100%). Molecular diversity from COI sequences was measured using DnaSP 5.0 (Librado & Rozas, 2009) with the following variables: number of haplotypes (H), polymorphic sites (S), haplotype diversity (h), and nucleotide diversity (π).
To visually represent the relationships among the mtDNA haplotypes of the sampled E. cardinalis individuals, we performed a network analysis using HAPLOVIEWER, which turns trees built from traditional phylogenetic methods into haplotype genealogies (Salzburger, Ewing, & Von Haeseler, 2011). The phylogenetic tree used for this graphical representation was obtained by employing a maximum likelihood approach in PhyML 3.0 (Guindon et al., 2010) using a GTR model with four gamma-distributed rate categories as the substitution model.
Population genetic differentiation was inferred by estimating pairwise ΦST values among 11 sampling sites using the Tamura–Nei model of nucleotide substitution in Arlequin 3.5 (Excoffier & Lischer, 2010). This model was assessed as most suitable for our data using jModeltest 2.1 (Darriba, Taboada, Doallo, & Posada, 2012). The significance of ΦST values was tested by 10,000 permutations. To estimate migration rates between sampling sites, we performed maximum likelihood analysis using the coalescence method in Migrate 3.2.1 (Beerli & Felsenstein, 1999). The estimated parameters were 𝒫ij = θiMij, where 𝒫ij is the number of effective migrants from i to j, θ is mutation-scaled population size, and Mij is mij/μ (where mij is the immigration rate from population i into j, and μ is the mutation rate per generation).

2.5 Inferring the historical demography

Signatures of population demographic changes (bottlenecks or expansions) in E. cardinalis were first examined by Tajima’ D and Fu’s Fs (Fu, 1997; Tajima, 1989) statistics with 10,000 permutations in Arlequin 3.5 to determine whether E. cardinalis in the Beibu Gulf data conformed to or departed from neutral theory model expectations because of factors such as a population bottleneck or expansion. For neutral markers, significant negative Tajima’s D and Fu’s Fs values can be expected under population expansion.
Next, we calculated the mismatch distribution among haplotypes under the sudden expansion model. This measure quantifies the smoothness of the observed mismatch distribution. We employed parametric bootstrapping (1,000 replicates) as implemented in Arlequin 3.5 to test the goodness of fit of the observed mismatch distribution to that expected under the sudden expansion model using the sum of squared deviations and raggedness index (R index) (Excoffier, 2004; Ray, Currat, & Excoffier, 2003; Schneider & Excoffier, 1999). Populations that underwent substantial expansion are expected to exhibit unimodal mismatch distributions with a low R index.
In addition, we used two methods to estimate the time of E. cardinalis population expansion in Beibu Gulf. First, the expansion time was directly estimated from the mismatch distribution with the statistic τ (tau) and translated into absolute time in years (t) using the equation τ = 2ut, where u is the mutation rate of the sequence and is calculated as u = 2µk; k is the number of nucleotides in the sequence and µ is the mutation rate of the mtDNA gene per generation (Rogers & Harpending, 1992). Following Cantatore et al. (1994), mutation rates between 1% and 3% per million years were selected for our mitochondrial analysis.
Subsequently, we inferred historical demography from effective population size estimates over time using the Bayesian skyline plot (BSP) method in BEAST 1.8.0 (Alexei J. Drummond & Rambaut, 2007; A. J. Drummond, Rambaut, Shapiro, & Pybus, 2005). An HKY model with among-site rate heterogeneity across all branches that assumed a strict molecular clock was used for this calculation. Markov chains were run for 2.5×107 generations and sampled every 1,000 generations, with the first 2,500 samples discarded as burn-in. Three replicates were run and combined to separately analyze the COI datasets. Other parameters were set as default values. TRACER 1.5 was used to visualize the posterior probabilities of the Markov chain statistics and to calculate a statistical summary of the genetic parameters.