Genetic variation of the small yellow croaker (Larimichthys polyactis) inferred from mitochondrial DNA provides novel insight into the fluctuation of resources

1. Introduction

The small yellow croaker, Larimichthys polyactis, is an economic fish of the order Perciformes, mainly inhabiting estuaries and coastal waters (Wu et al., 2012; Xiao et al., 2009). It is known as famous seafood product in China, Japan, D. P. R. Korea and R. O. Korea (Zhang and Cheng, 2005). Since the 1950s, however, the resources and catch of L. polyactis have undergone huge fluctuations. In the 1950s, L. polyactis was in the fishery boom period, whose average annual production could reach 120 000 t. However, in the 1960s and 1970s, it has gradually decreased due to overfishing. With the continuous increase of fishing pressure, the resource of much marine fish has declined severely in the 1980s. Some protection policies have been developed to protect the fish stock since the 1990s (Lin et al., 2008), playing a key role in the recovery of L. polyactis resources, which have kept its market landing at a high level in recent years (Lin et al, 2011, Han et al., 2019). However, the individual miniaturization younger-age trend of L. polyactis is becoming more and more serious (Tang and Zhou, 1999). To adapt to the environment and fishing pressure, the following changes have been detected: the generation time of the species has shortened from 6 years old to 2 years old, the distribution area of the spawning stock has extended, the environmental characteristics of the spawning have ground changed obviously, the relative individual fecundity has increased significantly, and the egg diameter have decreased (Lin et al., 2009). To protect and fully utilize the resources of L. polyactis, many researchers have focused on the studies of its resources (Lin and Cheng, 2004), physiology (Zheng et al., 2020), biology (Xu and Cheng, 2009) and morphology (Lin et al., 2010).

The summer closed fishing in China was formulated based on the study of fish stock (Chen et al., 2009). This policy certainly played a role on protecting resources of L. polyactis. However, L. polyactis spawns before May, and the fishing targets were reported to almost be young individuals, which indicated it may suffer overfishing before the 4-month moratorium (Yan et al., 2019). Besides, because of high fishing stress after the open of fishing, the conservation effect accumulated in summer was almost exhausted (Lin and Cheng, 2004).

Genetic diversity is an important component of biodiversity, and could generally reflect effective population size of species (Giovannoni et al., 1990). High genetic diversity within species usually suggested a large population size and can provide its evolutionary potential for long-term survival. Some studies have analyzed the genetic diversity of L. polyactis populations collected a dozen years ago based on different markers such as random amplified polymorphic DNA (RAPD), amplified fragment length polymorphism (AFLP) and mitochondrial DNA (Meng et al., 2003; Han et al., 2009; Xiao et al., 2009; Kim et al., 2012; Li et al., 2013). However, since L. polyactis populations continued to be under high fishing pressure, it is necessary to analyze its current genetic diversity status.

Molecular markers, important method for studying the mechanism of resource restoration, have been used widely in research of L. polyactis. Mitochondrial DNA (mtDNA) is a type of molecular marker which can give insights into contemporary levels of gene flow (Moritz et al., 1987). The control region of the mtDNA has the fastest mutation rate, and it has been widely used in conservation genetics and population genetics because of the characteristics of clear mechanism, simple structures and low molecular weights (Wirgin et al., 2000; Tokuyama et al., 2020). In this study, we analyzed the genetic diversity and genetic differentiation of L. polyactis based on mitochondria DNA control region. This work, in addition to offering a wealth of biological novelties which can indicate genetic structure of L. polyactis populations, also can provide the foundation for resources protection and policy setting.

2. Materials and methods

2.1 Study area and sampling

A total of 168 L. polyactis specimens collected from 7 sites (the Yellow Sea: YT, RS, QD, LYG, YC; the East China Sea: ZS, WZ) were used in this study (Fig. 1, Table 1). The muscle tissue of these individuals was stored in alcohol for total genomic DNA extraction using phenol/chloroform method (Sambrook et al., 1989).

Figure  1.  Sampling sites of Larimichthys polyactis. The dashed red circle represents the North Yellow Sea and Bohai Sea overwintering group, and the dashed green circle represents the South Yellow Sea and East China Sea overwintering group (Xu and Cheng, 2009).

DownLoad: Full-Size Img PowerPoint

Table  1.  Sampling information and genetic diversity parameter of Larimichthys polyactis from different locations

Site	Collection date	Sample number (N)	Number of polymorphic site (n)	Haplotype number	Haplotype diversity (h)	Nucleotide diversity (π)
YT	Dec. 2019	24	43	21	0.989±0.015	0.018±0.010
RS	Apr. 2019	24	40	23	0.958±0.012	0.015±0.008
QD	Aug. 2019	24	39	24	1.000±0.012	0.016±0.009
LYG	Aug. 2019	24	43	21	0.989±0.015	0.015±0.008
YC	Oct. 2020	24	44	22	0.989±0.017	0.019±0.011
ZS	Aug. 2019	24	36	23	0.996±0.013	0.011±0.006
WZ	Nov. 2020	24	42	21	0.989±0.015	0.017±0.009
All samples		168	110	134	0.993±0.002	0.016±0.009

 | Show Table

DownLoad: CSV

2.2 Polymerase Chain Reaction (PCR) and sequencing

The primers were designed to amplify the fragment of the mtDNA control region in this study were F2-CGGACGTCGGGGGTTAAAT and R2- ATGGGGAGCAACCACAAGAA. The PCR was performed in A300 Fast Thermal Cycler (LongGene Scientific Instruments Co., Ltd., China). PCR amplifications were carried out in volume of 25 μL, the reaction system containing deionized water (17.35 µL), dNTPs (2 µL, 10 mmol/L), 10×PCR buffer (2.5 µL, 10 mmol/L Tris–HCl pH 8.3, 50 mmol/L KCl, 1.5 mmol/L MgCl₂), forward and reverse primers (1 µL, respectively), Taq polymerase (0.15 µL), DNA template (1 µL). The amplification conditions as follows: 5 min denaturation at 94°C, 38 alternating cycles of 45 s at pre-denaturation 5 min (95°C), denaturation 45 s (94°C), annealing 45 s (50°C), extension 45 s (72°C), 38 cycles and extension 10 min (72°C). The amplification products were detected by 1% agarose gel electrophoresis. The products were sequenced commercially by Tsingke Biotech Co., Ltd. (China). All samples were sequenced with both directions to ensure the accuracy of these fragments.

2.3 Data analysis

DNASTAR software was used to align and edit all sequences (DNASTAR Inc., USA). Haplotypes based on sequence were obtained using DnaSP v.5.0 (Librado and Rozas, 2009). Phylogenetic relationships based on haplotypes were analyzed by MEGA v.5.0 (Tamura et al., 2011). Genetic diversity index such as polymorphic sites, nucleotide diversity (π), haplotype diversity (h) (Nei, 1987) was calculated using ARLEQUIN v.3.0 (Excoffier et al., 2007). The D test of Tajima and the Fs test of Fu indicating population demographic expansion condition were also carried out by ARLEQUIN v.3.0 (Tajima, 1989; Fu, 1997). Mismatch distribution analysis and neutrality tests were used to examine the historical dynamic of L. polyactis. F-statistics (F_st) and molecular variance (AMOVA) with Kimura-2-parameters model of substitution were calculated in ARLEQUIN v.3.0 in order to evaluate population structure (Kimura, 1980; Weir and Cockerham, 1984; Excoffier et al., 1992).

The effective population size was estimated using Lamarc v.2.0 (Kuhner, 2006). The parameter θ (θ=2N_eμ, where N_e was the effective population size, μ was the mutation rate) which indicated the effective population size of each population was estimated by Lamarc v.2.0 (Kuhner, 2006; Felsenstein et al., 1999). This software was also used to calculate gene flow (M=m/μ, m was migration rate) and parameter g of growth rate.

The Bayesian skyline plot (BSP) was generated with BEAST v.2.3.0 (Bouckaert et al., 2014) and Tracer v.1.7.1 (Rambaut et al., 2018). The strict molecular clock and stepwise skyline were selected as a model, and the sequence divergence rate for the control region of L. polyactis was approximately set to 3%−12% every million years.

3. Results

3.1 Sequence variation and genetic differentiation of L. polyactis

A total of 447-bp control region sequences were obtained for 168 samples from 7 geographic populations. There were 110 polymorphic sites, which defined 75 substitutions consisting of 74 transitions and 22 transversions. The average base composition content was 27.5% for A, 14.4% for C, 22.3% for G, 35.8% for T. The number of haplotypes was 134, which was more than the number of polymorphic sites, indicating that there was strong homoplasy among populations of L. polyactis. The diversity indices of each population were shown in Table 1. Besides, we analyzed the genetic diversity of the control region sequences of L. polyactis collected from the Yellow Sea (Qingdao) and the East China Sea (Zhoushan) based on the previous studies. The result showed that genetic diversity of two sites remained relatively constant in the present study compared with previous data (Table S1)

3.2 Phylogenetic relationships and population genetic structure

Clustering analysis of control region haplotypes was conducted by Neighbor-Joining method. The Neighbor-Joining tree was constructed using Larimichthys crocea as outgroup. From the tree topologies, we concluded that no obvious genealogical branch and geographical structure were found in control region haplotypes of 7 geographic populations, and haplotypes were scattered throughout the tree (Fig. 2). This result was further supported by the haplotype network (Fig. 3). The starburst structure also suggested a very recent origin for most haplotypes.

Figure  2.  The phylogenetic analyzes for Larimichthys polyactis investigated using Neighbor-joining tree. Larimichthys crocea sequence was used as an outgroup.

DownLoad: Full-Size Img PowerPoint

Figure  3.  Haplotype network showing genetic relationship among control region haplotypes for Larimichthys polyactis. The size of circles is proportional to haplotype frequency.

DownLoad: Full-Size Img PowerPoint

The genetic structure of L. polyactis populations was estimated based on pairwise F-statistics and AMOVA analyses. The F_st among different geographical populations was −0.008 (YT vs. WZ) to 0.050 (QD vs. YC). No genetic differentiation which consisted with the distribution pattern caused by geographically isolated distribution pattern was detected, revealing high genetic homoplasy of L. polyactis in the Yellow Sea and East China Sea. Three values of the pairwise F_st estimates were negative (Fig. 4) (YT vs. LYG, YT vs. WZ and LYG vs. ZS), indicating that the variation within populations was greater than that between population, and most of them were not statistically significant after correction for multiple tests.

Figure  4.  Population genetic structure analyse (F-statistics, F_st) of Larimichthys polyactis. * is significant genetic differentiation.

DownLoad: Full-Size Img PowerPoint

The result of AMOVA analysis demonstrated that the genetic variation among groups was 0.04%, while the 98.96% variation was detected within populations when the populations were divided into two groups (group from YT, RS, QD, LYG, YC; and group from ZS, WZ). The results implemented under two patterns of gene pools showed most of total genetic variances were attributed to differences among individuals within populations, indicating no population structure existed throughout the sample locations of L. Polyactis (Table 2).

Table  2.  AMOVA of Larimichthys polyactis populations of different locations

Source of variation	Variance components	Percentage variation	$\varPhi $-Statistics	P
1. One gene pool (YT, RS, QD, LYG, YC, ZS, WZ)
Among populations	0.04	1.02	$\varPhi_{\rm{st}} $=0.010	<0.001
Within populations	3.73	98.98	$\varPhi_{\rm{st}} $=0.010	<0.001
2. Two gene pool (YT, RS, QD, LYG, YC) (ZS, WZ)
Among groups	0.01	0.04	$\varPhi_{\rm{st}} $ = 0.010	<0.001
Among populations within groups	0.03	1.00	$\varPhi_{\rm{st}} $ = 0.010	<0.001
Within populations	3.72	98.96	$\varPhi_{\rm{st}} $ = 0.010	<0.001

 | Show Table

DownLoad: CSV

3.3 Population historical dynamics

Neutrality test was used to estimate population historical dynamics of L. polyactis. Tajima’s D (D =−1.355 6, P=0.084) and Fu’s Fs (Fs=−13.572 3, P=0.000) statistics were all significantly negative, indicating that the population departure significantly from mutation-drift equilibrium (Table 3). Therefore, we inferred that the recent large quantity expansion took place in population of L. polyactis.

Table  3.  Results of Tajima’s D and Fu’s Fs statistics based on different locations

Population	Tajima’s D			Fu’s Fs
	D	P		Fs	P
YT	−1.229 1	0.099		−10.934 6	0.000 0
RS	−1.357 3	0.079		−9.854 1	0.000 0
QD	−1.692 2	0.031		−11.136 2	0.000 0
LYG	−1.349 4	0.075		−12.157 4	0.000 0
YC	−1.013 2	0.139		−12.162 2	0.000 0
ZS	−1.151 3	0.118		−18.778 7	0.000 0
WZ	−1.010 1	0.145		−11.239 8	0.000 0
All samples	−1.355 6	0.084		−13.572 3	0.000 0

 | Show Table

DownLoad: CSV

The mismatch distribution for L. polyactis showed distinct unimodality, which supported the conclusion of population sudden expansion (Fig. 4). The BSP was also consistent with the result of neutrality test, and further uncovered the timing of expansion of L. polyactis (Fig. 5). The results indicated that the effective population of L. polyactis along the Yellow Sea and East China Sea have increased from about 75 000 to 40 000 years ago and leveled off by 30 000 years ago (Fig. 4).

Figure  5.  The mismatch distribution under sudden expansion model (solid line).

DownLoad: Full-Size Img PowerPoint

Figure  6.  The Bayesian skyline plots for Larimichthys polyactis populations.

DownLoad: Full-Size Img PowerPoint

3.4 Effective population size and gene flow

The effective population size, g and gene flow was estimated by the program Lamarc. Some discrepancies of effective population size were detected among different geographic populations. The lowest effective population size in geographic populations was found in QD population (3.42×10⁵−1.37×10⁶), while the largest size was detected in LYG population (8.83×10⁵−3.53×10⁶) (Table 4).

Table  4.  Estimates of effective population size (N_e), grow parameter (g) and migration among Larimichthys polyactis populations

Population	$\theta $ ($\theta $=Ne μ)	Ne	g		Gene flow (Ne×m)
				YT	RS	QD	LYG	YC	ZS	WZ
YT	0.091	7.58×105−3.03×106	345	−	64	58	207	0	371	110
RS	0.045	3.75×105−1.50×106	733	435	−	145	351	265	0	102
QD	0.041	3.42×105−1.37×106	197	296	169	−	215	245	0	85
LYG	0.106	8.83×105−3.53×106	418	467	0	62	−	158	143	210
YC	0.055	4.58×105−1.83×106	402	346	196	203	152	−	213	152
ZS	0.077	6.42×105−2.57×106	468	0	204	576	0	202	−	265
WZ	0.062	5.17×105−2.07×106	383	265	163	358	362	356	321	−
Note: − represents no data; μ, matation rate.

 | Show Table

DownLoad: CSV

The g of all populations was positive (Table 4), indicating L. polyactis populations have been growing. The gene flow among populations was asymmetric and high value of gene flow was found among most geographical populations, which supports the result that there was strong homoplasy in populations of L. polyactis (Table 4).

4. Discussion

The mtDNA fragments have been proven to be very effective and reusable for population genetic study and widely used in population genetics (Simons et al., 2001; Gao et al., 2020). It could be relatively accurate to calculate the genetic parameters of different populations and then to detect the genetic structure, genetic diversity and population historical dynamics of the populations (Simons et al., 2001; Song et al., 2014; Shan et al., 2020). Larimichthys polyactis is an economically important species in the northwestern Pacific Ocean. There are two main spawning grounds for L. polyactis in offshore area of China, the North Yellow Sea and Bohai Sea spawning ground and the South Yellow Sea and East China Sea spawning grounds. In this study, our sampling sites are near the main spawning grounds of L. polyactis, where the resource of this species is relatively concentrated. These sites are the most important and representative distribution areas (Lin et al., 2009). Mitochondrial DNA control region was used to analyze the current genetic diversity status and the genetic differentiation among these sites. Larimichthys polyactis from different geographic populations showed high haplotype diversity and low nucleotide diversity, which was a common phenomenon in marine fish (Ishikawa et al., 2001; Zheng et al., 2012). The result indicated that this species may have experienced a recent sudden expansion, which was supported by the results of neutrality test and mismatch distribution (Slatkin and Hudson, 1991). Small populations would lose partial low frequency of haplotype due to genetic drift, and the rapid growth of the population size within a short time may lead to decrease of genetic drift (Rogers and Harpending, 1992; Dann et al., 2013). Some rare haplotypes were thus preserved, resulting in high genetic diversity in population of L. polyactis.

The higher genetic diversity and larger effective population size of L. polyactis were detected in this study, which indicated that its resource had a tendency to gradually restore in recent years. Severe population declines can be recognized by decreased genetic diversity (Glenn et al., 1999). Our genetic assessment of L. polyactis provides a baseline from which L. polyactis populations can be appraised in the future. The average age of L. polyactis was 5.120 years old in the 1960s, whereas it was only 1.001 years old since the 1990s. The individual miniaturization younger-age trend is becoming more and more serious (Tang and Zhou, 1999). The control region sequences of L. polyactis collected from the coastal waters of China from May 2003 to April 2005 (GenBank: KF412650-KF412764) (Xiao et al., 2015) were downloaded for the comparative analysis. The results showed that h, π and effective population size in the present study were slightly higher than those data from Xiao et al (2009). The predators with larger size may have declined due to overfishing which might cause the increase of L. polyactis resources (Lin et al., 2011). Some laws and regulations protecting the spawning grounds and breeding period of L. polyactis have been formulated since the 1990s, which also effectively preserved its resources (Cheng et al., 2004). Moreover, in recent years, formulation of strict ban on fishing resulted in the gradually declined in the annual catch of L. polyactis (Chen et al., 2020). These measures have played an important role in the recovery of L. polyactis resources. Even so, recently, miniaturization and younger-age trend of L. polyactis were still acute (Lin et al., 2008). To restore the resource of L. polyactis, it is essential to promulgate more reasonable laws and regulations on fishery protection. The small mesh size of the net prevents the small size individuals escaping from the opening meshes. Increasing the mesh size of the net can release small size or juvenile individuals than the control net while keeping the large size catch (Zhang et al., 2014). Additionally, fishery spawning ground and adjacent sea area of L. polyactis should be protected as a priority. It is necessary to create fishery spawning ground protection area based on biological characteristics, fishery survey data and geographical environment (Wang et al., 2021).

The F_st is an effective parameter to estimate differentiation among populations (Teacher et al., 2013). In this study, the conventional population statistic F_st also revealed no significant genetic structure throughout the examined range of L. polyactis and detected strong gene flow existed among the populations. Ikeda (1964) found that recruitment phenomenon existed from north to south in L. polyactis populations among the East China and the Yellow seas, which was consistent with this result. Similarly, the AMOVA conducted on two patterns of gene pools was used to illustrate the species population genetic structure and showed that most variations existed within populations in both one gene pool and two gene pool patterns. The previous studies using mitochondrial DNA markers, microsatellite markers and Restriction-site Associated DNA Sequencing all showed no significant differentiation among L. polyactis populations in the coastal waters of China, which supported the result in this study (Xiao et al., 2009; Kim et al., 2012; Zhang and Cheng, 2005). Low genetic differentiation in L. polyactis populations can be attributed to its strong dispersal ability. The pelagic eggs may float with the ocean current, resulting in gene flow among populations (Lin et al., 2011). Moreover, for spawning, feeding and wintering, L. polyactis populations will migrate in a large range every year (Kim et al., 2006). Theoretically, if N_e×m>1, gene flow among populations can block genetic differentiation caused by genetic drift (Palumbi, 1994). High gene flow was detected among L. polyactis populations in this study, which also indicated that there was no significant genetic differentiation in populations. This may also account for the high genetic diversity. The assessment of population structure can be vital for maintaining a productive fishery (Giovannoni et al., 1990). There was no evidence for genetic subdivision of L. polyactis in this study. However, it is premature to conclude that L. polyactis is composed of a single panmictic stock based on only mtDNA. Therefore, recommendations for future evaluation and management research of fishery resources are studies based on fishery survey combined with molecular markers of higher coverage and larger data sets, such as whole-genome resequencing.

The results of mismatch distribution, neutrality test all showed that L. polyactis populations had a recent expansion, which was support by the result of BSP. Besides, we estimated the population expansion time of L. polyactis using the software of BEAST. The population expansion time of L. polyactis in offshore area of China dated to be approximately in 75 000−40 000 years ago. The population expansion time estimated in the previous study (about 49 300 years ago) was basically consistent with this study (Xiao et al., 2009). This time was in the late Pleistocene and the global climate has changed dramatically, which exerted great influence on survival distribution pattern of global organisms (Herbert et al., 2001; Gao et al., 2020). The effective population size of L. polyactis maybe had a sudden expansion after they experienced severe bottleneck effect because of environmental degradation. Therefore, L. polyactis populations may also be not at equilibrium yet between genetic drift and population migration because of sudden expansion (Palumbi, 1994).

5. Conclusions

The present study analyzed the variation in genetic diversity and population structure of L. polyactis based on mitochondrial DNA which can provide novel insight into the fluctuation of resources. In summary, L. polyactis populations had high genetic diversity and effective population size, which indicated that its resource had a tendency to gradually restore in recent years. No significant genetic differentiation was detected among populations. Besides, like many marine fish, the effective population size of L. polyactis had a sudden expansion recently. This study can provide a significant basis for the formulation of laws and regulations on fishery protection.