Threadfin Porgy (Evynnis Cardinalis) Haplotype Pattern and Genetic Structure in Beibu Gulf, South China Sea

Threadfin porgy (Evynnis cardinalis) is one of the important commercial fishing targets of bottom trawl fishery in the northern South China Sea. It is mainly threatened by overexploitation and listed as endangered species in the IUCN Red List. To investigate the demographic history and genetic structure of E. cardinalis population, partial sequences of the mitochondrial cytochrome c oxidase subunit I (COI) gene were obtained from 162 individuals collected from Beibu Gulf, South China Sea. In total, 44 different haplotypes were identified, and the dominant haplotype was found in all sampling sites. Across the dataset, nucleotide diversity was low, whereas haplotype diversity was high. Low pairwise comparisons of ΦST and high gene flow among sampling sites revealed a genetically homogeneous population structure in Beibu Gulf, indicating a single panmictic stock of E. cardinalis in this area. The star-like haplotype network, unimodal mismatch distribution, and significantly negative Tajima’s D and Fu’s Fs values indicated recent population demographic expansion of E. cardinalis. The mismatch distribution and Bayesian skyline plot results indicated that E. cardinalis from Beibu Gulf might have experienced colonization and demographic expansion due to sea level fluctuations during the late Pleistocene.


INTRODUCTION
Marine fish are generally deemed to have high dispersal potential because of their high moving capability at both larval and adult stages, and the absence of obvious physical barriers to dispersal (Caley et al., 1996;Hellberg, 2009). Theoretically, a species with the higher dispersal capability, the lower the genetic structure between populations. Therefore, in marine fish, especially migratory fish, the signal of population differentiation is weak and difficult to observe because of high levels of gene flow (Gandra et al., 2020). Understanding the genetic diversity of commercially important species is critical to implement protection policies and management regulations (Araki and Schmid, 2010). Genetic diversity (both within and between populations) greatly influences the adaptive potential of species to environmental changes, and ultimately determines their long-term resilience to ecological disturbances (Pauls et al., 2013). Additionally, knowledge of genetic structure is useful for understanding the migration routes, areas, and seasons of fish spawning. Such information can help fisheries managers define the spatiotemporal scales over which they can implement effective stock management and conservation plans (Bradbury et al., 2008). Therefore, from the resource management perspective, understanding the levels of gene flow among populations and patterns of genetic diversity is a fundamental issue. The effectiveness of marine protected areas depends on both their ability to self-recruit (reproductive potential) and the spillover of adults and export of larvae to nearby fished areas (Harrison et al., 2012;Le Port et al., 2017). Over the last 30 years, population genetic studies have become an essential tool for stock management and conservation of coral reef, estuarine, and coastal populations because population genetic studies are useful for estimating genetic diversity and the ability to survive anthropogenic activities such as overfishing, habitat degradation, eutrophication, invasive species, and pollution (Ryman et al., 1995;Ruzzante et al., 1998;Gill and Kemp, 2002;Gaither et al., 2010;Machado et al., 2020). Mitochondrial DNA markers (mtDNA) are widely employed to detect population structure in marine species because they have large number of copies, high mutation rates, generally maternal inheritance, and almost nonexistent recombination. The primary advantages of mtDNA are the inheritance pattern and nonexistent recombination: speciation events over multiple generations can be traced through maternal cloning, and male dispersal does not homogenize the population (Prugnolle and de Meeus, 2002). These factors make mtDNA markers particularly suitable indicators of the genetic structure of marine populations with high gene flows, such as zooplankton and migratory fishes (Wang et al., 2013;Xu et al., 2019;Machado et al., 2020).
Threadfin porgy, Evynnis cardinalis (Lacepède 1802), occurs in the Indo-West Pacific from Japan, Korea, and China to Vietnam and Indonesia (Chen and Qiu, 2005). The species is mainly distributed at depths of 30-60 m, but also can occur to 100 m depth (Iwatsuki and Carpenter, 2014). E. cardinalis can be found in various types of bottoms, but is more common close to coral reefs or rough bottoms. Small individuals are very common in shallow, such as sheltered bays, while larger fish often live in deeper water (Eggleston, 1974). E. cardinalis is one of the main commercial fishing targets of bottom trawl fishery in the northern South China Sea and is thought to have three geographical stocks in the northern South China Sea: the Taiwan Strait, South China Sea, and Beibu Gulf stocks (Zhang et al., 2020). According to the surveys data, more than 90% catches of the total of sparids is E. cardinalis in Beibu Gulf. E. cardinalis is a migratory fish that is captured all year round with significant seasonal differences among different fishing areas. In Beibu Gulf, E. cardinalis undergoes seasonal migration and is found toward the northeastern shallow area of the gulf in late autumn and early winter. Then, spawning takes place in northwest Weizhou Island in early spring, and the recruits disperse widely in the northern area of the gulf and begin to disperse to the south in late spring (Iwatsuki and Carpenter, 2014). However, this species exhibits late maturity and longevity, which predispose it to impacts from heavy exploitation. In the northern South China Sea, the Beibu Gulf stock of E. cardinalis declined by 58% from 2001 to 2005, whereas the Taiwan Strait stock declined by 62% from 1993 to 2008. Therefore, a recent IUCN Red List assessment reported E. cardinalis as endangered, and this species is mainly threatened by overexploitation (Iwatsuki and Carpenter, 2014). Previous studies of E. cardinalis in Beibu Gulf have focused on feeding habits, growth and mortality, ecological distribution, phylogeny, and stock density (Chen and Qiu, 2003;Chen and Qiu, 2005;Zhang et al., 2007;Zhang et al., 2014;Cai et al., 2017;Zhang et al., 2020). However, no information is available to date on the population genetic structure of E. cardinalis in the northern South China Sea, and this prevents reliable stock assessments and protection policy formulation. The objective of this study was to provide a population genetic analysis using a portion of the mitochondrial cytochrome c oxidase subunit I (COI) gene to assess the population genetic diversity pattern and historical demography, and estimate E. cardinalis expansion time in Beibu Gulf.

Study Area
Beibu Gulf (17°-22°N, 105°-110°E) is located in the northwestern part of the South China Sea and has a long coastline that belongs to China and Vietnam. It is an approximately 128,000 km 2 of semi-closed gulf that ranges from Leizhou Peninsula, Qiongzhou Strait, and Hainan Island to east Vietnam coast, and extends to the Guangxi coast in the north (Ma et al., 2010). The bottom of Beibu Gulf is flat and deepens from the northwest to the southeast, with a depth of 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 important source of fishery products for coastal areas because of high productivity and rich biodiversity, which benefit from its unique geographic location and climatic conditions (Chen et al., 2009).

Sample Collection
All E. cardinalis 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. We set eleven sampling sites in the study area of Beibu Gulf, which covered over 45,000 km 2 ( Figure 1). In total, 11-16 specimens of each sampling locations were used for the DNA analyses after morphological identification. The dorsal fin or muscle were removed from each specimen and preserved in absolute ethanol at −20°C.

DNA Extraction, Amplification, and Sequencing
Genomic DNA were extracted from the dorsal fin or muscle tissue using the TIANamp Marine Animals DNA Kit (TIANGEN, China). The concentration used as PCR template was adjusted to an A260 of approximately 0.05-0.2. Fragments of the mtDNA cytochrome c oxidase subunit I gene were amplified from total genomic DNA by the polymerase chain reaction. The primer sets FishF1 and FishR1 (Ward et al., 2005) were used for PCR amplification. Each 50 µl PCR tube contained of 5 µl PCR buffer, 4 µl of 25 μM MgCl 2 , 2.5 µl DNA template, 5 µl CoralLoad concentrate, 0.5 µl of 25 µM solution of each primer, 1 µl of 10 µM dNTPs, 0.25 µl TopTaq DNA polymerase, and 31.25 µl ddH 2 O (QIAGEN, Germany). The procedure for PCR amplification were in turn as follows: an initial step of 95°C for 3 min; 35 cycles with 95°C for 30 s (denaturation), 54°C for 30 s (annealing), and 72°C for 1 min (extension); then followed by 5 min at 72°C (final extension) on a 2,720 Thermal Cycler (Applied Biosystems, United States). PCR products were visualized on 1.2% agarose gels and the most intense bands were selected for sequencing. All PCR products were bidirectionally sequenced on an ABI 3730XL automated sequencer with both forward and reverse primers.

Genetic Analysis
The authenticity of all COI sequences was first verified by BLAST search in GenBank (BLASTn, megablast algorithm) and compared with the highest match (99-100%).
Then all sequences were assembled in Bioedit (Hall, 1999) and aligned using the CLUSTALW multiple algorithm under default options. Ambiguous sequences were trimmed after alignment. Molecular diversity from COI sequences was measured using DnaSP 5.0 (Librado and Rozas, 2009) with the following variables: number of haplotypes (H), polymorphic sites (S), haplotype diversity (h), and nucleotide diversity (π).
In order to visually represent the relationships between the mtDNA haplotypes of the sampled E. cardinalis individuals, we performed a haplotype network analysis using HAPLOVIEWER, which turns traditional phylogenetic tree into haplotype genealogies network (Salzburger et al., 2011). The phylogenetic tree used for this visual 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 estimated by pairwise Φ ST values among 11 sampling sites using the Tamura-Nei model of nucleotide substitution in Arlequin 3.5 (Excoffier and Lischer, 2010). The model was assessed as most suitable for our data using jModeltest 2.1 (Darriba et al., 2012). The significance of the pairwise Φ 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 and Felsenstein, 1999). The estimated parameters were P ij θ i M ij , where P ij is the number of effective migrants from i to j, θ is mutation-scaled population size, and M ij is m ij /μ (where m ij is the immigration rate from population i into j, μ is the mutation rate per generation).

Inferring the Historical Demography
Signatures of population demographic patterns (bottlenecks or expansions) in E. cardinalis were first examined by Tajima (1989) and Fu (1997) 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 Fu's Fs and Tajima's D values can be expected under population expansion.
Then, we calculated the mismatch distribution among sampling sites under the sudden expansion model. This measure quantified the smoothness of the observed mismatch distribution. We used parametric bootstrapping (1,000 replicates) as implemented in Arlequin 3.5 to test the goodness of fit of the observed mismatch distribution to the expected under the sudden expansion model using raggedness index (R index) and the sum of squared deviations (Schneider and Excoffier, 1999;Ray et al., 2003;Excoffier, 2004). Populations that underwent substantial expansion are expected to display unimodal mismatch distributions with a low R index.
Furthermore, we used two methods to estimate the population expansion time of E. cardinalis in Beibu Gulf. First, the population expansion time was directly calculated from the mismatch distribution of τ (tau) statistic and converted to the absolute time in years (t) using the equation τ 2 ut, where u is the mutation rate of the sequence and is calculated as u 2 µk; k is the nucleotides number of the sequence and µ is the mutation rate of the mtDNA gene per generation (Rogers and 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 (Drummond et al., 2005;Drummond and Rambaut, 2007). This calculation used a strict molecular clock with an HKY model assumed with among-site rate heterogeneity across all branches. Markov chains were run for 2.5 × 10 7 generations and sampled every 1,000 generations, with the first 2,500 samples discarded as burn-in. Three replicates were performed and combined to analyze COI datasets, respectively. The rest of parameters were set as default values. TRACER 1.5 was employed to visualize the posterior probabilities of the Markov chain statistics and to estimate a statistical summary of the genetic parameters.
The haplotype network analyses of the mtDNA from E. cardinalis in Beibu Gulf showed a typical star-like pattern, with the most common haplotype, H1, located at the center of the haplotype network and surrounded by many low-frequency haplotypes that were divergent from H1 by only one or two mutations (Figure 2). Most of the haplotypes had low-frequency, and 27 haplotypes were uniquely present at only one sampling site. Only 11 out of the 44 haplotypes were observed at more than two sampling sites. However, several of the abundant haplotypes were present at multiple sampling sites. For example, the dominant widespread haplotype H1 distributed in all 11 sampling sites ( Figure 3); the largest distance between two sampling sites harboring this haplotype was more than 300 km. H6 and H8 were distributed at eight and seven sampling sites, respectively. Moreover, between six and ten haplotypes were found per sampling site, but there was no clear pattern of geographical variation (Figure 3). The highest number of haplotypes (H 10) was observed at sites 8 and 10 ( Table 1).
The average genetic differentiation, Φ ST , of E. cardinalis within Beibu Gulf was 0.031 (range, −0.0657-0.0671). All pairwise Φ ST values among the sampling sites were small and not significant    Table 2). The effective number of migrants per generation among sampling sites ranged from 5.98 to 45.76. The rates of migration from S7 to S1 were the lowest, whereas those from S1 to S3 were the highest ( Table 3).

Demographic History
The Fu's Fs and Tajima' D tests of E. cardinalis in Beibu Gulf were significantly negative , p < 0.01; Tajima's D −2.446, p < 0.01); this indicated that E. cardinalis may have experienced population expansion ( Table 4). The mismatch distribution appeared to be unimodal (Figure 4), which was consistent with the expected distribution under a sudden expansion model (R index 0.055, p 0.26).
On the basis of the mutation rate of 1-3% of mitochondrial genes per million years (Cantatore et al., 1994) and τ value of all data (1.531, 95% confidence interval: 1.281-1.941), the expansion time for E. cardinalis in Beibu Gulf was estimated to have occurred from approximately 62-21 ka.
The BSP analysis indicated that the E. cardinalis haplotypes in Beibu Gulf coalesced approximately 40 ka when a mutation rate of 2% of mitochondrial genes per million years was used for analysis. However, the BSP pattern revealed a coalescence time of 30-70 ka when using a mutation rate of 1-3% per million years ( Figure 5). Despite these different coalescence times, the BSP patterns showed similar tendencies and indicated that a steady population expansion took place between 5 and 30 ka.

DISCUSSION
Studies of population connectivity with mitochondrial markers provide critical information on gene flow and genetic relationships between neighboring populations (Turner et al., 2004;García et al., 2008). Many studies showed that mitochondrial markers are highly effective for revealing marine fish genetic diversity and population connectivity (Lavergne et al., 2014;Gao et al., 2019;Machado et al., 2020). In this study on E. cardinalis, mtDNA sequence analysis of specimens from Beibu Gulf revealed no significant genetic differentiation among sampling sites, with low Φ ST values indicating genetic homogeneity. In contrast to freshwater fish, marine fish are generally expected to show relatively low genetic differentiation across their distribution. It is mainly attributed to genetic exchange being maintained by adult mobility throughout Beibu Gulf during reproduction, and through the passive dispersal of eggs and larvae due to the lack of noticeable physical barriers in pelagic oceans (Grant and Bowen, 1998;Hellberg, 2009;Machado et al., 2020). The dominant widespread haplotype H1 was found in all 11 sampling sites, which also indicated high dispersal potential of E. cardinalis at planktonic egg, larval, or adult stages in Beibu Gulf.
Previous studies suggested that E. cardinalis breed once a year in Beibu Gulf (Chen and Qiu, 2003;Hou et al., 2008). E. cardinalis gonads begin to develop in November and begin to spawn from December to February. The population concentrated in the northern Beibu Gulf during spawning. In early spring, parent fish mainly occur in the northeast of the gulf after spawning, and juveniles concentrate in shallow nearshore of this area in late spring. Thereafter, juveniles gradually migrate southwest and widely disperse in deep waters of Beibu Gulf in summer or early autumn (Zhang et al., 2020).
In addition, the dispersal pattern of E. cardinalis was also impacted by circulation in Beibu Gulf. In spring, the density gradient and monsoon wind drive the ocean current from northeast to southwest in the gulf. The surface current velocity reaches 30 cm/s, and the current in the middle layer is approximately 5-10 cm/s (Gao et al., 2017). The direction of the spring currents roughly coincides with the migration route of E. cardinalis. Therefore, the seasonal migration and ocean current may be responsible for gene exchange between different locations, and E. cardinalis shows low levels of genetic differentiation in Beibu Gulf. If we refer to the biological description of a single stock as given by Ihssen et al. (1981), ''a stock is an intraspecific group of randomly mating individuals with temporal and spatial integrity,'' then the lack of distinct spatial boundaries and genetic substructure (low Φ ST values) revealed by genetic analyses indicated that E. cardinalis in Beibu Gulf belong to a single stock.
The presence of a single stock in Beibu Gulf indicates that geographical isolation might block gene exchange between the Beibu Gulf stock and the other 2 E. cardinalis stocks, the Taiwan Strait and South China Sea stocks. In Beibu Gulf, the circulation, Hainan Island, and Leizhou Peninsula could act as barriers that impede free dispersal of E. cardinalis into this gulf from other areas of the South China Sea (Gao et al., 2017). Our findings and similar investigations conducted elsewhere demonstrated that marine fish that inhabit coastal waters usually constitute a single panmictic stock. For example, Rodrigues et al. (2008) revealed that Cynoscion acoupa from northern Brazil belongs to a single stock, even though it covers at least 1,260 km of coastline. Gao et al. (2019) reported a high level of genetic homogeneity in the Pholis fangi population around Yellow Sea and Bohai Sea, and suggested it should be considered as a single panmictic stock. Hoolihan et al. (2006) also reported a homogeneous distribution of Spanish mackerel (Scomberomorus Frontiers in Environmental Science | www.frontiersin.org August 2021 | Volume 9 | Article 726533 commerson) throughout the Gulf of Oman, Arabian Gulf, and Arabian Sea on the basis of mtDNA analyses. In addition, mtDNA sequence regions are particularly suitable for inferring historical processes that might be responsible for the contemporary geographical distribution of marine species because they are more prone to genetic drift than nuclear genes and have a relatively small effective population size (Slatkin and Hudson, 1991;Avise, 1994). In our study, the haplotype network of E. cardinalis from Beibu Gulf exhibited a star-like and unstructured pattern with a predominance of scattered. The dominant haplotype (carried by 45% of the specimens) was at the center of the haplotype network and surrounded by many haplotypes that diverged from the dominant haplotype by only few mutations. Most surrounding haplotypes were unique to each sampling site and showed few differences among them (Figure 2). Similar star-like haplotype networks have been observed for other species in different coastal areas: Terapon jarbua, which consists of a panmictic stock from the Socotra Archipelago to the Hadhramout Coast along the wider Gulf of Aden (Lavergne et al., 2014); and Pogonias courbina, which did not display distinct structure along the coast of the southwestern Atlantic Ocean (Machado et al., 2020).
A star-like haplotype network pattern, low nucleotide diversity, and high haplotype diversity are often considered consequences of recent population expansion linked to the Pleistocene environmental changes (Craig et al., 2007;Pereira et al., 2009;Liu et al., 2011). The recent demographic expansion of E. cardinalis from Beibu Gulf is also supported by the unimodal mismatch distribution and significantly negative Fu's Fs and Tajima's D values. The population expansion of E. cardinalis in Beibu Gulf, which was directly estimated from the mismatch distribution, started 62-21 ka before present, corresponding to the late Pleistocene. BSP analysis indicated steady population expansion that started around 30 ka. Both two methods of estimated period of population expansion are consistent with the environment changes of the northern South China Sea in the Pleistocene.
E. cardinalis is mainly distributed from 30 to 60 m depth, and spawns in coastal habitats and shallow shorelines. Therefore, the E. cardinalis distribution is closely related to historical sea level fluctuations. When sea level was 120 m lower than the present level during the last glacial maximum of the Pleistocene, the northern South China Sea included Beibu Gulf, which was part of the South China continent, Hainan Island, and Taiwan Island were connected to mainland China, and the entire South China Sea was separated from the Indian Ocean to form a semi-closed basin (Wang and Sun, 1994;Voris, 2000). Similar to other terrestrial species, E. cardinalis may have moved and survived in a potential glacial refuge during this period, such as the semiclosed South China Sea (Hewitt, 1999). In the late Pleistocene, the sea level was still 30 m lower than the present level, but the glaciation began to disappear and the sea water gradually poured into Beibu Gulf through the Qiongzhou Strait (Lu et al., 2003). An initial population of E. cardinalis may have immigrated to Beibu Gulf from neighboring areas after it was filled with sea water and sufficiently deep. This initial panmictic stock quickly colonized the empty novel environment under the founder and priority effects, and might have experienced rapid population expansion when favorable conditions occurred (Shulman et al., 1983;Boileau et al., 1992).

CONCLUSION
A homogeneous population structure with low genetic diversity and star-like haplotype pattern was revealed for E. cardinalis in Beibu Gulf of the northern South China Sea. Pairwise comparisons of Φ ST and the effective number of migrants between all sampling sites in Beibu Gulf indicated a single panmictic stock of E. cardinalis. This is mainly attributed to the free genetic exchange of E. cardinalis during reproduction. However, E. cardinalis from the eastern South China Sea cannot disperse into Beibu Gulf because of the circulation patterns and ocean geographical features such as the semi-closed sea of Beibu Gulf. This molecular evidence revealed that E. cardinalis from Beibu Gulf might have experienced colonization and population expansion during the late Pleistocene due to sea level fluctuations.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: National Center for Biotechnology Information (NCBI) BioProject database under accession numbers: MW881382-MW881425.

ETHICS STATEMENT
The animal study was reviewed and approved by South China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences.