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.