Demographic history and population genetic analysis of Decapterus maruadsi from the northern South China Sea based on mitochondrial control region sequence

Sample collection

Two hundred forty-one adult D. maruadsi were collected between March 2010 and May 2014 from 11 populations in the northern SCS (Fig. 1 and Table 1) namely, the Shantou inshore water (ST, N = 37), the Zhuhai inshore water (ZH, N = 20), the Yangjiang inshore water (YJ, N = 10), the Maoming inshore water (MM, N = 30), the Zhanjiang inshore water (ZJ, N = 10), the Wenchang inshore water (WC, N = 28), the Sanya inshore water (SY1, N = 25), the Sanya offshore water (SY2, N = 10), the Baimajing inshore water (BMJ1, N = 26), the Baimajing offshore water (BMJ2, N = 18), and the Beihai inshore water (BH, N = 27). After morphological identification, muscle tissues were taken from individual fish and stored in 95% ethanol at −20 °C until DNA extraction.

DNA extraction, PCR, and sequencing

Total genomic DNA was extracted from muscle tissue using the phenol-chloroform protocol (Huang et al., 2002). Polymerase chain reactions (PCRs) for mtDNA control regions were performed in 50 µl containing 39.5 µL dH₂O, 5 µL 10 × EasyTaq buffer, 1 µL dNTP mixture (10 mM), 1 µL forward primer (10 µM, 5′-AACTCCCAAAGCTAGGATTCT-3′), 1 µL reverse primer (10µM, 5′-TTTGTGCTTGCGGGGCTTT-3′) (Niu et al., 2012), 2.5U EasyTaq DNA polymerase (TransGen Biotech Co., Ltd., Beijing, China), and 2 µL template DNA (50–100 ng). PCR reactions were performed in a Veriti 96 Well Fast Thermal Cycler (Applied Biosystems, Foster City, CA, USA) under the following parameters: initial denaturation at 94 °C for 5 min followed by 35 cycles of 45 s denaturation at 94 °C, 45 s annealing at 55 °C, and 1 min elongation at 72 °C. The amplification cycle finished with a final elongation at 72 °C for 10 min. All sets of PCR included a negative control reaction tube in which all reagents were included except the template DNA. The PCR products were separated by electrophoresis on 1.2% agarose gel and the Trans 2K DNA Marker (TransGen Biotech Co., Ltd., Beijing, China) was used for identifying the target gene size. Successfully amplified PCR products were sent to GenScript Biotech (Nanjing, China) for sequencing.

Data analysis

Sequences were edited and aligned with the DNASTAR software package (Allex, 1999). All sequences were trimmed using the MEGA 6.06 software (Tamura et al., 2013). Molecular diversity indices, including haplotype diversity (h), nucleotide diversity (π), the mean number of pairwise differences (K), indels, and polymorphic sites, were calculated using the software Arlequin 3.5 (Excoffier & Lischer, 2010). The models of the nucleotide sequence evolution were evaluated by jModelTest v2.1.10 (Darriba et al., 2012), and the Hasegawa-Kishino-Yano (HKY) nucleotide substitution model with a gamma distribution (0.811) and a proportion of invariable sites (0.792) was chosen for the evolutionary analyses.

A Bayesian phylogeny was inferred for all haplotypes under the HKY+G+I model by MrBayes 3.2.6 software (Ronquist et al., 2012). The analysis was run for 10,000,000 generations until the average standard deviation of split frequencies fell below 0.01. The first 25% of burn-in trees were discarded and the subsequent trees were combined to produce a phylogram. Ultimately, a Bayesian inference tree was created and viewed using FigTree 1.4.3 software. The congener, Decapterus macarellus (GenBank accession number: KM986880) and Decapterus macrosoma (GenBank accession number: KF841444) mtDNA control regions were used as the outgroup. The genealogical relationship among the haplotypes was further evaluated using the media-joining network constructed with Network 5.0.0.3 software (Bandelt, Forster & Rohl, 1999).

Hierarchical analysis of molecular variance (AMOVA) with 10,000 permutations was performed to test the population genetic structure of D. maruadsi in the northern SCS. To quantify the genetic divergence, we computed pairwise F_ST (fixation index) values between populations with 10,000 permutations. An exact test of population differentiation based on haplotype frequencies was also examined with 100,000 steps in the Markov Chain. AMOVA, F_ST , and the exact test were all calculated using the software Arlequin 3.5.

Historical expansion patterns of D. maruadsi for all populations pooled together were calculated using three different approaches. First, Tajima’s D and Fu’s Fs neutrality tests with 10,000 permutations were calculated in Arlequin 3.5 to test the selective neutrality of the D. maruadsi control region. Secondly, nucleotide mismatch distribution was performed in Arlequin 3.5 and DnaSP 5.10.01 (Librado & Rozas, 2009) to detect the population expansion history. The fit between the unimodal observed and expected distributions was evaluated by the test statistics of goodness-of-fit, including the sum of squared deviation (SSD) and Harpending’s raggedness index (Hri). The statistical significance was estimated with 1,000 bootstrap replicates. Finally, the Bayesian skyline plot was implemented in BEAST 1.7.5 under the HKY+G+I nucleotide substitution model and strict molecular clock model to estimate the change in effective population size across time (Drummond et al., 2012). The Markov Chain Monte Carlo (MCMC) analysis was run for 1 billion generations, yielding high effective sample sizes (ESS) of at least 200, of which the first 10% were discarded as a burn-in. The analysis results were visualized using Tracer 1.7 (Rambaut et al., 2018).

Molecular rate and its calibration

The molecular rate of control region seems to vary among different bony fishes, such as 3% per myr for rainbow fishes (Zhu et al., 1994), 6.5–8.8% per myr for cichlid fishes (Sturmbauer et al., 2001), and 5–10% per myr for Arctic charr (Brunner et al., 2001). However, there has been no a molecular rate specialized for control region of the fish from Carangidae family. D. maruadsi has a small body size, high metabolic rate, and short generation time, suggesting that a rapid molecular clock may be more appropriate for the D. maruadsi control region. Thus, a per lineage mutation of 3.75% per myr calculated from the reported control region divergence rates (5–10% per myr) was used for the Bayesian skyline plot analysis of D. maruadsi in the present study (Brunner et al., 2001; Liu, 2006; Li et al., 2018).

Several researchers have suggested that it is necessary to calibrate the molecular mutation rate in the analysis of historical population demography (Ho et al., 2005; Ho et al., 2008; Crandall et al., 2012). If the molecular rate is not calibrated on timescales of less than approximately 1–2 Myr before the present, the population expansion time will be overestimated. According to the research results of Ho et al. (2008) and Crandall et al. (2012), 3–10 × faster mutation rates were used to calibrate the population expansion time in the present study. Even so, the evolutionary time estimated based on molecular rate provides an approximate time frame for demographic history.

Genetic diversity and haplotype patterns

The partial length sequences of the control region were 828–832 bp. The calculated genetic variations were summarized in Table 1. Sequence alignment revealed 65 transitions, 23 transversions, 11 indels, and 91 polymorphic sites. The h of the total population was 0.953, with a range of 0.905–0.980. The π and K ranged from 0.00269 (ZJ) to 0.00849 (SY1) and from 2.233 (ZJ) to 7.053 (SY1), respectively. This showed that the SY1 population had the highest level of genetic diversity, while the ZJ population had the lowest. The overall π and K were 0.00443 and 3.691, respectively, indicating low levels of genetic diversity in the D. maruadsi population in the northern SCS.

Ninety-one polymorphic sites were observed in the 241 sequences, and defined 94 haplotypes (GenBank accession number: MK294641–MK294734), of which 20 (21.28%) were shared among populations, with a relative frequency ranging from 0.83% to 15.77% (Table 2). For example, H21 (38 specimens, 15.77%) was found in all populations and showed the highest relative frequency among all haplotypes. The main haplotype, H7 (20 specimens, 8.30%), was shared by all populations except the ZJ population. Haplotypes H29 (nine specimens, 3.73%), H57 (18 specimens, 7.47%), H63 (14 specimens, 5.83%), and H74 (16 specimens, 6.64%) occurred in six, eight, seven, and seven populations, respectively. These shared haplotypes also showed different relative frequencies in different populations. The remaining 74 haplotypes (78.72%) were all unique to a single population (private haplotypes), and most of them (71, 95.95%) were singletons.

Phylogenetic relationship

The Bayesian inference tree (Fig. 2) for the 94 haplotypes of D. maruadsi showed that haplotypes in every population were scattered throughout the tree and no clustering was observed that corresponded to the sampling localities. A haplotype network (Fig. 3) revealed the presence of six central dominant haplotypes (H7, H21, H29, H57, H63, and H74) at high frequencies across 6–11 populations. All other low-frequency haplotypes radiated in a star-like fashion out from the dominant haplotypes by a few mutational steps. The sharing of interior haplotypes and the endemic distribution of terminal haplotypes indicated the genetic homogeneity among the populations and recent population expansions. This was consistent with the Bayesian inference tree of the 94 haplotypes.

Genetic structure

Pairwise F_ST values between populations ranged from −0.053 (between YJ and BMJ2) to 0.058 (between ZJ and SY2) and were statistically non-significant (Table 3). Hierarchical AMOVA results (Table 4) revealed most of the genetic variations (>99%) were attributable to the variations among individuals within populations and only a small amount of the genetic variations (<1%) were explained by variations among geographical groups and among populations within groups. The F_CT, F_ST, and F_SC values of hierarchical AMOVA analysis were all low and not statistically significant. The results demonstrated that no significant genetic differentiation existed among populations. The finding could be further supported by the results of the exact tests among the population which revealed no significance (exact P value = 0.0986 ± 0.0517).

Demographic history

Both Tajima’s D statistics (Tajima’s D =  − 2.201, p = 0.0005) and Fu’s Fs test (Fu Fs = −25.180, p = 0.0002) were negative and highly significant in pooled samples. The mismatch distributions of D. maruadsi for all populations pooled together showed a unimodal pattern under the sudden expansion model (Fig. 4A) and constant population size model (Fig. 4B). Under the sudden expansion model, the observed distributions were included in the 90%, 95%, and 99% confidence intervals obtained by simulated distribution. The SSD values (SSD = 0.001, p = 0.840) and Hri (Hri = 0.007, p = 0.936) were small and statistically non-significant, which showed no significant deviation between the observed distributions and expected distributions. Under the constant population size model, a poor goodness-of-fit was found between the observed and expected distributions. Both the neutrality test and mismatch distribution analyses suggested a recent population demographic expansion for D. maruadsi populations in the northern SCS. However, a low migrants per generation (M < 50) was detected under the spatial expansion model. This may lead to a larger variance of mismatch distribution and incorrectly estimating the parameters for spatial expansion (Ray, Currat & Excoffier, 2003; Excoffier, 2004). Therefore, the spatial expansion model was not chosen to analyze the range expansion history of D. maruadsi.

Bayesian skyline plot (Fig. 5) revealed a pattern of glacial decline and postglacial growth in the effective population size based on the mutation rate of 3.75% per myr. The population decline happened about 18,000–37,000 years before present (YBP) that coincided with the LGM. Following the LGM and continuing to the present (0–18,000 YBP), D. maruadsi has undergone a rapid episode of population growth, which resulted in an about 100-fold increase in the maximum effective population size. The results suggested that D. maruadsi in the northern SCS experienced demographic expansion after the LGM. Use of the 3–10 ×faster mutation rates shortened the timing of population expansion and the calibrated time was 1,800–6,000 YBP, which further confirmed that D. maruadsi in the northern SCS experienced a recent population expansion during the Holocene.

Artificial statistical bias adjacent geographical populations

To detect artificial statistical bias due to a low sample size from YJ, ZJ, and SY2 populations, some adjacent geographical populations were combined into three geographic groups for genetic differentiation analysis, including the western Guangdong group (WGD: ZH, YJ, MM, ZJ), Sanya group (SY: SY1, SY2), and Baimajing group (BMJ: BMJ1, BMJ2). Genetic homogeneity was demonstrated by AMOVA (Table 5), pairwise F_ST (F_ST =  − 0.0189–0.00674, P = 0.235 − 0.992), and exact tests (exact P value = 0.178 ± 0.0758) among three geographic groups and the other three populations (ST, WC, BH). This indicated that low sample sizes for YJ, ZJ, and SY2 have not influenced the results of this study.

Demographic history

This study revealed D. maruadsi in the northern SCS experienced a rapid demographic expansion after a period of low effective population size during the Holocene. The results suggested better conditions (e.g., the enlargement of habitat and the significant increase of paleo-productivity) after the LGM driving the population expansion of D. maruadsi.

During the LGM, the sea-level dropped 100–120 m in the South China Sea (Wang & Sun, 1994). The dropping of the sea-level exposed shallow coastal continental shelves and turned the SCS into a semi-closed land-locked sea (about one-fifth smaller than it is now) (Wang & Sun, 1994; Wang, 1999), greatly reducing the habitat and population abundance of many marine fishes in the continental shelf (He et al., 2015; Zhu, Cheng & Rogers, 2016). D. maruadsi widely distributes in the northern shallow continental shelf of the SCS within 180 m and its spawning grounds are located mainly in shallow nearshore waters less than 120 m (Deng & Zhao, 1991). With the falling sea-level during the LGM, most of its spawning grounds were exposed and gradually moved from shallow-water to the deeper water. Thus, D. maruadsi should have undergone intensive spawning ground loss and change, which may have caused a decrease in the effective population size, abundance, and distribution. The surviving D. maruadsi had to contract into the potential glacial refugia (semi-closed SCS). However, the subsequent rising sea-level of the postglacial-Holocene provided new niche space for D. maruadsi, leading to sustained and rapid growth in the effective population size of this species. Similar patterns of population dynamics have been revealed for some marine fishes in response to the postglacial rising sea-level in the northern SCS (Liu et al., 2006; Xiao et al., 2014; Gao et al., 2018).

The paleo-productivity of the northern SCS was mainly influenced by the fluctuant East Asian monsoon and Kuroshio Current (Wang, 1999; Yu & Qiu, 2016). After the LGM, the East Asian summer monsoon gradually strengthened, which caused high precipitation and increased land-source influx into the northern SCS (Jian, Wang & Kienast, 1999). Accordingly, high levels of terrigenous nutrients from the Asian continent were transported into the northern SCS ecosystems through the Pearl River, Han River, Nanliu River, Red River, and other sources of fluvial runoff (Jian et al., 1999; Gao, 2013; Qu & Huang, 2019). Moreover, the enhanced summer monsoon could drive upwelling containing high nutrient input into the northern SCS (Jian, Wang & Kienast, 1999). Obviously, the enhanced summer monsoon during the deglaciation period increased the nutrient input of the northern SCS, which led to higher paleo-productivity than that of the LGM. During the LGM, the Kuroshio Current was forced eastward due to the formation of a land bridge between Taiwan and the Ryukyu Islands, which prohibited its flow into the SCS (Ujiié et al., 2003). By contrast, the postglacial Kuroshio Current branches into the northern SCS through the Luzon Strait along with the rise of the sea-level (Nan, Xue & Yu, 2015). The Kuroshio intrusion induces upwelling and meso- and micro-scale eddies, strengthening the vertical mixing of the water layer and nutrient supplies to the upper ocean (Yu & Qiu, 2016). Additionally, the Kuroshio intrusion could transport large amounts of nutrients into the northern SCS (Xu et al., 2018), promoting the increase of paleo-productivity in the northern SCS. The combination of the above factors strongly prove that the northern SCS had significantly higher paleo-productivity during deglaciation (> 30 gC/(m²⋅10³a)) than the LGM (< 25gC/(m²⋅10³a)) (Kuhnt, Hess & Jian, 1999; Zhao & Wang, 1999). The significant increase of paleo-productivity in the northern SCS after the LGM has a close correlation to the recent rapid population growth of D. maruadsi.

All in all, LGM-Holocene climate oscillations gave rise to drastic changes in the East Asian monsoon circulation, sea-level, ocean current patterns, and upwelling density of the western Pacific marginal seas (Wang, 1999; Jian, Wang & Kienast, 1999; Ujiié et al., 2003). These changes resulted in enlargement of nearshore habitat and significant increase of primary productivity, which were believed to be the main cause of the rapid demographic expansion of D. maruadsi. Similar phenomena have been reported in other fish species, such as Pampus echinogaster (Li et al., 2018) and Trichiurus japonicus (He et al., 2014). However, it is not yet possible to discriminate the relative impact of changes in habitat or paleo-productivity based on the existing mtDNA data.

Population structure

In this study, low levels of genetic diversity and high levels of genetic homogeneity were detected in the D. maruadsi population from the northern SCS. Contemporary genetic diversity and structure is influenced by demographic history (Hewitt, 1996; Hewitt, 2000). The rapid population growth after a period of low effective population size could enhance the retention of new mutations, but it might not accumulate enough nucleotide diversity in the short term (Avise, Neigel & Arnold, 1984; Rogers & Harpending, 1992). Therefore, the population decline and the subsequent rapid population expansion may lead to the loss of diversity and explain the low levels of genetic diversity in D. maruadsi in the northern SCS. The recent rapid population expansion also indicates that D. maruadsi in the northern SCS is far from equilibrium due to the lack of sufficient evolutionary time (Slatkin, 1993), which should be an important historical reason for genetic homogeneity among D. maruadsi populations. A similar genetic structure and explanation was found in E. japonicus (Liu et al., 2006) and Pennahia argentata (Han et al., 2008) along the coastal waters of China.

Several other factors aside from the historical process may have also led to the lack of significant genetic differentiation among D. maruadsi populations in the northern SCS. First, the planktonic eggs and larvae of D. maruadsi (Deng & Zhao, 1991) could disperse on ocean currents to facilitate genetic exchange over great distances. There are many annual or semi-annual ocean currents in the northern SCS. The northward currents, such as the southwest monsoon drift, South China Sea warm current, the coastal current of eastern Guangdong, and the coastal current of eastern Hainan (Qiao, 2012), are beneficial for the northward dispersal of the southern groups (e.g., SY1, ZJ, and MM populations). The southward currents, including the northeast monsoon drift, Guangdong coastal current, and the Kuroshio branch in South China Sea (Li, Zhao & Wasiliev, 2000), promote populations from the eastern Guangdong to migrate south. Furthermore, D. maruadsi pelagic larvae on both sides of Leizhou Peninsula and Hainan Island could cross the Qiongzhou Strait for genetic exchange along the residual current in the Qiongzhou Strait (Gao, 2013). The planktonic eggs and larvae of D. maruadsi may travel great distances on ocean currents in the northern SCS, which may be sufficient to produce genetic homogeneity. Secondly, adult D. maruadsi from the northern SCS spawn throughout the year, and the spawning ground spreads from the Taiwan shoal to Beibu Gulf (Deng & Zhao, 1991). The overlapping spawning time and spawning ground could maintain genetic connectivity amid geographic separation.

The lack of significant genetic structure implies that D. maruadsi in the northern SCS should be considered a single panmictic population. However, some significant biological differences are found among the different populations of D. maruadsi in the northern SCS. The average fork length and body weight of D. maruadsi in Beibu Gulf and its mouth are 145 mm (dominant fork length = 110–160 mm) and 42 g (dominant body weight = 20–60 g), respectively; while those in the eastern Hainan Island are 158.7 mm (dominant fork length = 130–220 mm) and 102.5 g (dominant body weight = 70–140 g), respectively (Jia, Li & Qiu, 2005). Additionally, different peak spawning periods are found among the various geographical populations of D. maruadsi in the northern SCS (Deng & Zhao, 1991). For example, the peak spawning season for the Beibu Gulf population is from January to March, from May to June for the population from southeast Hainan Island, and from February to April for the populations from the Pearl River estuary to east Guangdong (Jia, Li & Qiu, 2005). The potential mismatch between the population’s genetic structure (genetic homogeneity) and biological characteristics (significant difference) reminds us to identify and manage D. maruadsi populations carefully. However, owing to maternal inheritance, the mtDNA control region sequence provides no information about male migration and has some limitations in respect of fishery management (Waples, Punt & Cope, 2008). Thus, future integrative studies based on nuclear markers with biparental inheritance (such as microsatellites and SNPs), and biological and ecological characteristics are needed for more accurate identification of the population structure and MUs of D. maruadsi in the northern SCS.