Effects of Late Pleistocene Climatic Fluctuations on the Phylogeographic and Demographic History of Japanese Scad (Decapterus maruadsi)

The Late Pleistocene-Holocene climate fluctuations have greatly influenced the phylogeographic structure and historical dynamics of many marine organisms in the western Pacific marginal seas. Here, we investigated the impact of Pleistocene glacial-interglacial cycles on the phylogeographic structure and demographic dynamics of Decapterus maruadsi, an economically important fish along the coast of the East China Sea (ECS) and northern South China Sea (NSCS). We obtained 430 concatenated sequences (Cyt b + control region, 1548–1554 bp) of D. maruadsi, including 246 newly sampled from the ECS and 184 previously determined from the NSCS. Genetic structure and phylogenetic analysis demonstrated a lack of significant population structure among 16 populations. Moreover, there was no significant differentiation among populations from Chinese coastal waters and northern Vietnam. Neutrality tests, unimodal mismatch distributions, Bayesian skyline plots, and the star-like haplotype networks all indicated a recent demographic expansion for D. maruadsi population during the Late Pleistocene-Holocene, explaining the low genetic diversity in D. maruadsi along the southeast coast of China. Notably, phylogenetic analyses and net genetic distances based on Cyt b jointly confirmed that 57 Cyt b haplotypes identified as D. maruadsi from the previously defined Sundaland-Rosario-Ranong clade actually represented D. russelli. These results not only reveal the complex effects of Pleistocene-Holocene climate fluctuations on the phylogeographic structure and demographic history of D. maruadsi but also provide useful genetic information for the management of genetic resources.


INTRODUCTION
Climate oscillations during the Pleistocene greatly altered the environment of marginal seas of the western Pacific (Wang and Sun, 1994;Liu et al., 2011), including the East China Sea (ECS) and South China Sea (SCS), further influencing the spatial geographic distribution and historical dynamics of many coastal marine organisms (Hewitt, 2000;Hewitt, 2004). During the glacial period, the shallow continental shelf of the Bohai Sea, Yellow Sea, and ECS was exposed and the coastline migrated about 1200 km from the coast of Bohai Bay to Okinawa Trough (Wang, 1999;Shen et al., 2011), while the SCS formed a semienclosed sac-shaped gulf and exposed approximately 0.7 million km 2 of continental shelf (Wang and Sun, 1994;Clark et al., 2009). The coastline retreat resulted in the formation of a large land bridge from the Taiwan Strait to the Taiwan-Ryukyu (Kimura, 2000), restricting dispersal routes of marine organisms and isolating populations (Rex, 1997;Alvarado Bremer et al., 2005;Ruban, 2007). The rising sea level during the post-glaciation period resulted in the reconnection of two glacial sea basins (Okinawa Trough and the SCS), providing more habitat and an opportunity for secondary contact between isolated populations. The intermittent isolation and connection of ocean basins during repeated glacial-interglacial cycles also caused pronounced fluctuations in the ocean current pattern, upwelling intensity, and oceanic salinity, which further influenced paleo-productivity (He et al., 2013;Li et al., 2018). Moreover, Yangtze River, Pearl River, and other river systems have the potential to affect paleo-productivity in nearshore coastal waters of the ECS and SCS (Jian et al., 1999a;Jian et al., 1999b). All of these factors could significantly affect the distribution and abundance of marine organisms across the southeast coast of China and may have the profound genetic imprints on their populations (Liu et al., 2007b;He et al., 2014).
Given the substantial geographic variation and complex historical processes in the ECS and SCS, phylogeographic structure and historical demography might differ among taxa in response to coastal environmental changes during glacial periods caused by climate oscillations (Ni et al., 2014). Based on previous research, three evolutionary patterns have been identified for marine organisms along the coast of China: (1) genetic homogeneity with no lineage formation, as observed in the mud crab Scylla paramamosain (He et al., 2010), small yellow croaker Larimichthys polyactis (Wu et al., 2012), and black sea bream Acanthopagrus schlegelii (Zhao et al., 2021); (2) past lineage diversification with a current sympatric distribution and genetic homogeneity, as observed in the large yellow croaker Larimichthys crocea , whitespotted conger Conger myriaster (Zou et al., 2020), and spotted scat Scatophagus argus (Peng et al., 2021); and (3) past lineage diversification with an allopatric distribution and genetic heterogeneity, as observed in the redlip mullet Chelon haematocheilus (Liu et al., 2007b), lemonfish Plectorhinchus flavomaculatus (Han et al., 2008), mitten crab Eriocheir sensu stricto (Xu et al., 2009), and mottled spinefoot Siganus fuscescens (Ravago-Gotanco and Juinio-Meñez, 2010). Taken together, species-specific phylogeographic patterns and historical dynamics are expected in the ECS and SCS.
The Japanese scad Decapterus maruadsi (Temminck and Schlegel, 1844) is widely distributed across the marginal sea of the western Pacific and is an important economic fish in China, especially in the ECS and the northern SCS (NSCS) (Zheng et al., 2014;Liu et al., 2016). The annual catch of D. maruadsi in the ECS and NSCS reached approximately 5.6 × 10 5 tons from 1996 through 2020, with a maximum catch of about 6.3 × 10 5 tons in 2006(Chinese Fishery Statistical Yearbook, 1997. A recent fishery resources survey has shown that D. maruadsi can nearly always be classified as overfished in the ESC and NSCS, including in the coastal waters of southern Zhejiang (annual exploitation rate = 0.61) (Cui et al., 2020) and the Beibu Gulf (annual exploitation rate = 0.63-0.78) (Geng et al., 2018). In addition, global climate change is fundamentally altering marine environments at an alarming rate (Doney et al., 2009). Continuous overfishing and environmental changes have substantially altered the phenotypic characteristics and population structure of D. maruadsi, including miniaturization, precocious puberty, and structure simplification (Zheng et al., 2014;Geng et al., 2018). Thus, the population structure of D. maruadsi should be closely monitored in the long-term to facilitate the sustainable utilization and scientific management of fishery resources.
According to the migration patterns and ecological characteristics, the ECS D. maruadsi can be assigned to three migratory populations (i.e., West Bank of Kyushu, western ECS, and Minnan to eastern Guangdong) (Zheng et al., 2003), while the NSCS D. maruadsi can be partitioned to seven ecological populations (i.e., Minnan-Taiwan shoal, the area off Jiazi in eastern Guangdong, the Pearl River estuary, the sea off of the Pearl River estuary, western Guangdong, Hainan Qinglan, and Beibu Gulf) (Deng and Zhao, 1991). Additionally, the ECS and NSCS D. maruadsi exhibit different migratory capacities and morphological characters (Jiang et al., 2012;Wang et al., 2021). Accordingly, it is reasonable to speculate that genetic differentiation exists between the ECS and NSCS D. maruadsi if gene flow is overcome by selective pressures and genetic drift. It is, however, noteworthy that the ECS and NSCS are a relatively homogeneous and contiguous environment characterized by the near absence of physical barriers, and ocean currents further increase the connectivity. Long-distance gene flow is possible between the ECS and NSCS D. maruadsi with potential migratory ability, and this is expected to constrain population differentiation. Therefore, the level of genetic differentiation between the ECS and NSCS D. maruadsi is not clear.
Previous population genetic studies of D. maruadsi have revealed high genetic homogeneity between two populations of Fujian coastal waters based on AFLP  and mitochondrial DNA control region (CR) sequences (Niu et al., 2012) and among 11 populations in the NSCS based on cytochrome b (Cyt b) (Niu et al., 2018) and CR (Niu et al., 2019). Notably, Jamaludin et al. (2020) proposed that D. maruadsi from the central Indo-West Pacific could be divided into two distinct genetic clades (Northern Vietnam and Sundaland-Rosario-Ranong) based on Cyt b ( Figure 1A). The Northern Vietnam clade contained two northern Vietnam populations (Cat Ba and Nghe An), the Sundaland-Rosario-Ranong clade consisted of nine Sundaland populations, and the central Vietnam population (Khanh Hoa) was an admixed group. However, the genetic distance (0.029-0.051) between the Northern Vietnam clade and Sundaland-Rosario-Ranong clade was significantly higher than that among populations within clades (0.001-0.003) (Jamaludin et al., 2020), emphasizing the need for further molecular studies to validate these two genetic clades. Furthermore, no study to date has directly evaluated differences between populations from Chinese coastal waters and the central part of the Indo-West Pacific region.
In this study, we use Cyt b and CR sequences to investigate the phylogeographic structure and demographic history of D. maruadsi in the ECS and NSCS. The goals of the study were (a) to determine whether repeated glacial-interglacial climate cycles during the Pleistocene influenced the phylogeographic structure of D. maruadsi, (b) to characterize demographic history, and (c) to evaluate the molecular phylogenetic relationship between the Northern Vietnam clade and Sundaland-Rosario-Ranong clade defined by Jamaludin et al. (2020). Our results provide insight into the evolutionary history of D. maruadsi along the coastal areas of the ECS and NSCS, which is critical for species conservation.

Ethics Statement
The experimental animal protocols in the present study were reviewed and approved by the Animal Experimental Ethics Committee of Guangdong Ocean University, China.

ECS Sample Collection and Sequencing
A total of 246 D. maruadsi individuals were collected from eight geographic populations (Taizhou, Ningde, Changle, Pingtan, Putian, Xiamen, Dongshan, and Penghu) in the ECS from July 2013 to September 2018 ( Figure 1 and Table 1). After morphological identification, muscle tissues were obtained from individual fish and preserved in 95% ethanol at -20°C until DNA extraction.

NSCS Samples and Sequence Integration
Cyt b (722 bp) and CR (826-832 bp) sequences of 184 D. maruadsi individuals sampled from eight localities across the NSCS determined in our previous studies were also evaluated (Niu et al., 2018;Niu et al., 2019). Thus, a total of 430 Cyt b and CR sequences of 16 geographic populations from the southeast coast of China were included in analyses. Detailed information on these samples is provided in Figure 1 and Table 1. After sequencing, the sequenced peak maps of all samples firstly were adjusted manually by Chromos software (http://technelysium. com.au/wp/chromas/) to remove the low-quality sequences. Then, partial sequences of Cyt b and CR were manually edited and checked using DNAStar 7.1 (DNASTAR, Inc.). Multiple alignments were generated using ClustalX 1.81 (Thompson et al., 1997). The sequences of Cyt b and CR were manually trimmed to the same length using MEGA 7.0 (Kumar et al., 2016). Sequence alignment found that the trimming Cyt b (722 bp) and CR (826-832 bp) sequences was located at 20-741 base and 14-843 base of their full-length (1141 bp for Cyt b and 843 bp for CR, GenBank accession number: MT818506), respectively. The trimming sequences and their concatenated sequence set (5'-Cyt b + CR-3') were used for subsequent analyses.

Genetic Diversity and Population Structure Analysis
Population genetic diversity and differentiation among 16 geographic populations (430 individuals) of D. maruadsi collected from the ECS (246 individuals) and the NSCS (184 individuals) were assessed. All analyses described below were performed using Arlequin 3.5 (Excoffier and Lischer, 2010) and considered all sites including insertions/deletions (indels) data. Haplotype diversity (h), nucleotide diversity (p), and other molecular diversity indices were evaluated, while number of haplotypes (Hn) and number of private haplotypes (Hp) were calculated. Pairwise F ST values between populations were evaluated with 10,000 permutations and all P values were adjusted by the Benjamini-Hochberg false discovery rate (FDR) (Benjamini and Hochberg, 1995). Analysis of molecular variance (AMOVA) was also performed with 10,000 permutations to characterize the population structure and genetic variation at two hierarchical levels: (1) all populations were classified into one group according to pairwise F ST values and (2) 16 populations were assigned to the ECS group and the NSCS group according to sample source to detect differentiation and geographical barriers (e.g., Taiwan Strait).

Phylogenetic Analysis
The evolutionary models for each marker separately and the concatenated sequence set including and excluding outgroups were evaluated by the lowest Bayesian information criterion (BIC) score using jModelTest v2.1.10 (Darriba et al., 2012). Bayesian inference tree was constructed under Hasegawa-Kishino-Yano model (HKY + I for Cyt b, HKY + I + G for CR, and HKY + I + G for Cyt b + CR) using MrBayes v3.2.6 (Ronquist et al., 2012) with three independent runs. Each run included four chains and each chain was run for 100 million generations with a sampling frequency of 100 generations and a diagnostic frequency of 1000 generations until the average standard deviation of split frequencies fell below 0.01 (Hall, 2016). The first 10% of trees were discarded as burn-in while subsequent trees were used to produce a phylogram. The tree was visualized using FigTree (v1.4.4; http://tree.bio.ed.ac.uk/ software/figtree/). Furthermore, the haplotype network was generated with the TCS network using PopART 1.7 (Leigh and Bryant, 2015). The net genetic distance was calculated under the Tamura-Nei model (TN93 + G, G = 0.75) model using MEGA 7.0. It is worth mentioning that all sites with insertions/deletions were incorporated in the phylogenetic analysis.

Demographic Analysis
Historical demography of D. maruadsi for the pooled population were inferred using three different approaches. First, we conducted neutrality tests by calculating Tajima's D (Tajima, 1989a;Tajima, 1989b) and Fu's Fs (Fu, 1997) with 10,000 permutations using Arlequin 3.5. Second, nucleotide mismatch distributions with 10,000 bootstrap replicates were also implemented in Arlequin 3.5 to detect historical population dynamics (Rogers and Harpending, 1992). The sum of squared deviations (SSD) and Harpending's raggedness index (Hri) were calculated to examine the departure between observed and expected distributions. The parameter t gained from the mismatch distribution was used to estimate the time of the population expansion with the equation t = t/2u, where u represents a neutral mutation rate of the entire sequence in each generation calculated as u = 2mk, where m is the nucleotide mutation rate and k is the length of the nucleotide sequence (Schneider and Excoffier, 1999). Finally, Bayesian skyline plot (BSP) was implemented under the HKY model (HKY + I for Cyt b, I = 0.76; HKY + I + G for CR, I = 0.84, G = 0.87; and HKY + I + G for Cyt b + CR, I = 0.85, G = 0.86) using BEAST v1.7.5 to identify the changes in effective population size (Ne) over time (Drummond et al., 2012). These analyses were run using the parameters as followed: 1.2 billion generations for Cyt b, 1.3 billion generations for CR, and 1.5 billion generations for Cyt b + CR, respectively. The effective sample sizes (ESS) of all runs were at least 200, and the first 10% were discarded as burn-in. These independent runs were visualized using Tracer v1.7.1 (Rambaut et al., 2018). As the divergence rate of Cyt b reported as 2% per million years (myr) in multiple bony fishes (Bermingham et al., 1997;Bowen et al., 2001;Durand et al., 2002) and that of CR used as 5-10% per myr for D. maruadsi in the previous studies (Niu et al., 2019), we took 1%, 3.75%, and 2.5% per myr as the putative nucleotide mutation rates for Cyt b, CR, and concatenated sequences (Cyt b + CR), respectively.

Reassessment of Molecular Phylogenetic Relationships Between the Northern Vietnam Clade and Sundaland-Rosario-Ranong Clade
To confirm the Northern Vietnam clade and Sundaland-Rosario-Ranong clade defined by Jamaludin et al. (2020) and investigate genetic differentiation between populations from the southeast coast of China and the central part of the Indo-West Pacific region, 80 Cyt b haplotypes (GenBank accession numbers: KX212002-KX212079, KY440109-KY440110) previously identified as D. maruadsi from two northern Vietnam populations (Cat Ba and Nghe An), one central Vietnam population (Khanh Hoa), and nine Sundaland populations were included in a phylogenetic analysis with our Cyt b sequence data. The Cyt b sequences of congeneric species Decapterus russelli from the NSCS (GenBank accession number: OM479430; the sample was collected in our previous fishery survey and was identified when fresh, and the sequence was also uploaded by our research team, ensuring the accuracy of the identification result) and the coastal waters of India (GenBank accession number: MN711693, mitochondrion complete genome), Omani Sohar (GenBank accession number: KX512731) (Damerau et al., 2018), and Indonesian Ambon (GenBank accession number: LC646650) (Kimura et al., 2022), Decapterus macarellus (GenBank accession number: KM986880) , and Decapterus macrosoma (GenBank accession number: KF841444)  from the SCS were downloaded from NCBI as outgroups. Thus, a collective data set of Cyt b sequences containing 124 sequences of taxa previously identified as D. maruadsi (44 defined in our study plus 80 from NCBI) and six sequences of three congeneric species (outgroups) was formed. These sequences were truncated to target the same Cyt b fragments (606 bp) and were then used to construct a phylogenetic tree under the HKY + G model and TCS network and to calculate net genetic distances under the TN93 + G model (G = 0.29). Other parameters were the same as the methods described above.

Genetic Diversity and Population Structure
A total of 430 Cyt b and CR sequences were obtained from 16 geographic populations along the ECS and the NSCS. Partial Cyt b sequences (722 bp) had 39 polymorphic sites with no insertions/deletions, and partial CR sequences (826-832 bp) contained 73 polymorphic sites and 10 insertions/deletions (Tables S1, S2). These variant sites defined 44 Cyt b haplotypes (CH1-CH44) and 134 CR haplotypes (CR1-CR134), respectively. Two mtDNA fragments were combined to obtain a concatenated sequence of 1548-1554 bp (Cyt b + CR) for each individual for subsequent analyses. As listed in Table 1, 430 concatenated sequences (Cyt b + CR) had 112 polymorphic sites and 179 haplotypes, designated as H1-H179 (GenBank accession numbers: OM728655-OM728833). Details for all D. maruadsi haplotypes obtained in this study were shown in Table S3. The overall h and p values for the 16 populations were 0.954 and 0.00302, respectively, indicating low levels of genetic diversity in D. maruadsi along the southeast coast of China. Furthermore, high haplotype diversities (0.892-0.987) in concurrence with low nucleotide diversities (0.00199-0.00371) were also detected in each population. The average h and p values for the ECS group were 0.940 and 0.00277, respectively, and those of the NSCS group were 0.970 and 0.00333, respectively, indicating that the level of genetic diversity was lower in the ECS group than in the NSCS group. The proportion of haplotypes and proportion of private haplotypes were also lower in the ECS group (Hn/N, 43.09%; Hp/N, 31.30%) than in the NSCS group (Hn/N, 54.89%; Hp/N, 35.87%). All the results of genetic diversity based on each marker separately (Tables S1, S2) were consistent with concatenated sequences (Cyt b + CR).
Pairwise F ST values showed non-significant differentiation (FDR-adjusted P values > 0.05) among 16 geographic populations (Tables 2, S4, S5). The AMOVA with all populations set as a single group revealed that the majority of the molecular variance (> 99.50%) was partitioned within populations, while a minor fraction of variance (< 0.50%) existed among populations ( Table 3). When all populations were divided into the ECS group and the NSCS group, more than 99.00% of genetic variation was explained by individuals within populations and only a small fraction of variation (< 1.00%) was attributable to the difference between groups ( Table 3). These results indicated that no obvious genetic structure was found among 16 populations or between the ECS group and the NSCS group in D. maruadsi.

Phylogenetic Relationships
The phylogenetic trees constructed using 44 Cyt b haplotypes ( Figure 2A) and 134 CR haplotypes ( Figure 2B) revealed that the haplotypes of each population were scattered throughout the trees and no pronounced clusters corresponded to the sampling sites across the southeast coast of China. Interestingly, three main haplogroups (labeled A, B, and C) with low node support values (0.56 and 0.57) were detected in the Bayesian inference tree based on concatenated sequences (Cyt b + CR) ( Figure 2C). Haplogroups A, B, and C included 47, 85, and 47 haplotypes comprising 104, 230, and 96 individuals, respectively. Three haplogroups were broadly sympatric across the southeast coast of China and there were no appreciable differences in their frequencies of along the ECS and NSCS. The net genetic distances among three haplogroups were 0.0007-0.0027.

Demographic History
Tajima's D and Fu's F S values were highly significantly negative for the pooled population ( Table 4), suggesting that there were population expansion events for D. maruadsi along the southeast coast of China. The nucleotide mismatch distributions of the pooled population were unimodal (Figure 4), further indicating a demographic expansion of D. maruadsi population. Both SSD and Hri values were low and not significant (Table 4), demonstrating a lack of significant deviation between the observed and expected distributions and supporting sudden expansion models. Based on the expansion parameter t and the putative nucleotide mutation rates (1%, 3.75%, and 2.5% per myr for Cyt b, CR, and the concatenated sequences, respectively), the expansion times of D. maruadsi population were dated back to the Late Pleistocene to Holocene period ( Table 4).
Bayesian skyline plots revealed a more elaborate demographic history of D. maruadsi along the southeast coast of China ( Figure 5). The pooled population showed no pronounced demographic changes from 183-262 kya (coalescence times) to 18-20 kya and experienced a rapid population expansion from 18-20 kya to 3 kya (3-20 kya for Cyt b and Cyt b + CR, 3-18 kya for CR). From 3 kya to present, the effective population size remained relatively stable for the pooled population. After expansion, the effective population sizes of D. maruadsi based on Cyt b, CR, and Cyt b + CR increased about 550-, 350-, and 500-fold, respectively.

Molecular Phylogenetic Relationships of the Northern Vietnam Clade and Sundaland-Rosario-Ranong Clade
Phylogenetic relationships between the Northern Vietnam clade and Sundaland-Rosario-Ranong clade were further analyzed based on a collective Cyt b data set containing 124 sequences (606 bp) of named D. maruadsi and six sequences of three outgroup taxa (D. russelli, D. macarellus, and D. macrosoma). As illustrated in Figure 6A, the phylogenetic tree revealed two distinct clades (labeled I and II) with high node confidence values (1.00). The TCS network ( Figure 6B) also showed two divergent clades (labeled I and II) separated by a number of mutational steps (25 steps), corresponding to the clades in the phylogenetic tree. Within clade I, 44 D. maruadsi haplotypes defined in this study (labeled red, CH1-CH44) clustered together with 23 haplotypes from two northern Vietnam populations (Cat Ba and Nghe An) and one central Vietnam population (Khanh Hoa) (blue mark, GenBank accession numbers: KX212002, KX212004, KX212051-KX212058, KX212060-KX212065, K X 2 1 2 0 6 8 -K X 2 1 2 0 6 9 , a n d KX 2 1 2 0 7 2 -K X 2 1 2 0 7 6 ) . 1350. However, the net genetic distances within clade I and II were 0.0042 and 0.0087, respectively, which were much lower than those between clades or species mentioned above. Based on these results, we conclude that Cyt b haplotypes within clade II belong to D. russelli, rather than D. maruadsi. In addition, we did not detect significant genetic differentiation among populations from Chinese coastal waters and three Vietnam populations within clade I.

Demographic Expansion
Neutrality tests, mismatch distribution, BSP analysis, and haplotype network all indicated a recent demographic expansion event for D. maruadsi along the southeast coast of China since the last glacial period, which may be related to a series of unique oceanographic and ecological changes from the Late Pleistocene to Holocene. During the glaciation, the East Asian monsoon prevailed over the southeast coast of China, leading to strong sea water mixing, improving shallow and midwater nutrient transport, and accelerating nutrient recycling in the marine ecosystem (Steinke et al., 2011;Huang and Tian, 2012). These changes eventually resulted in a substantial increase in primary productivity. Based on Marine Isotope Stages (MIS), previous studies have shown that the mass accumulation rates (MARs) of total phytoplankton in the NSCS were 62.5 × 10 -12 mg·cm -2 ·kyr -1 and 109.5 × 10 -12 mg·cm -2 ·kyr -1 in MIS 3 (25-59 kya) and MIS 2 (11-25 kya), respectively (He et al., 2013), strongly indicating that paleo-productivity was significantly higher during MIS 2 than during past periods. Also, the total content of all marine phytoplankton biomarkers, used as a total productivity indicator, revealed a higher productivity from the ESC during MIS 2-3 (11-59 kya) than MIS 4 (59-73 kya) (Xing et al., 2008;Vats et al., 2020;Vats et al., 2021). Therefore, it is postulated that the expansion of D. maruadsi from the Last Glacial Maximum (LGM) to the early Holocene is related to increased primary productivity. After the LGM, sea levels have risen 120-140 m Waelbroeck et al., 2002), providing substantial niche space for many marine organisms (including D. maruadsi) and accommodating extensive terrigenous nutrients rapidly from the shallow continental shelves of Chinese coastal waters. The sediment from the estuarian regions of the Pearl River and Yangtze River also showed that high levels of terrigenous nutrients were transported into the ECS and the NSCS ecosystems through the sources of fluvial runoff during deglaciation (Jian et al., 1999b;Liu et al., 2007a;Qu and Huang, 2019;Vats et al., 2020). These factors further explain the rapid growth of D. maruadsi from the deglaciation to the early Holocene. Over thousands of years of rapid growth, the effective population size of D. maruadsi increased approximately 350-to 550-fold (each marker separately and the concatenated sequences for the pooled population). Along with the increase in the effective population size, intraspecific competition should theoretically increase after the early Holocene. Moreover, the MARs of total phytoplankton in the NSCS were 22.2 × 10 -12 mg·cm -2 ·kyr -1 in MIS 1 (0-11 kya) (He et al., 2013), revealing a decrease in paleo-productivity during the Holocene. Also, compared with the level of paleoproductivity from the LGM to early Holocene, biomarkers showed that for the ECS, total productivity and individual diatom, dinoflagellate, and coccolithophorid productivity were all lower during the Holocene (Xing et al., 2008;Vats et al., 2020;Vats et al., 2021). The increasing intraspecific competition, along with the decrease of paleo-productivity, may have limited recent D. maruadsi population growth and maintained a relatively stable effective population size within the last few thousand years.

Phylogeographic Structure
The previous studies have demonstrated that the alternation of historical glaciation and deglaciation had profound impacts on phylogeographical structures of many marine fishes in the western Pacific (Liu et al., 2007b;Ni et al., 2014;Gao et al., 2020). During the Late Pleistocene glacial periods, the sea level of the western Pacific dropped approximately 120 m Waelbroeck et al., 2002), and the shallow coastal continental shelves of ESC and SCS were mostly exposed, leading to a large loss of habitat for D. maruadsi and the formation of a land bridge. Surviving D. maruadsi may be isolated, with separate populations in semi-closed SCS, Okinawa Through, as  well as in the additional potential glacial refugia. Long-term glacial geographic isolation may have contributed to the formation of lineages. However, in this study, only Bayesian inference tree based on concatenated sequences (Cyt b + CR) revealed three main haplogroups with low node support values (0.56 and 0.57), which was similar to the phylogeographic pattern of Lateolabrax japonicus (Liu et al., 2006). Pairwise F ST values and AMOVA indicated genetic homogeneity with no lineage formation among 16 populations along the southeast coast of China. This point was also supported by the phylogenetic analysis based on Cyt b, CR, and the concatenated sequences (no obvious clustering of haplotypes corresponding to localities was observed). In addition, the phylogenetic tree and haplotype network based on Cyt b sequences showed a lack of differentiation between 16 populations along the southeast coast of China and three populations within the Northern Vietnam clade (belonging to the NSCS) defined by Jamaludin et al. (2020). In general, D. maruadsi along the coast of the ECS and the NSCS exhibited no significant genetic structure. Similar genetic structure has been found in the cutlassfish Trichiurus japonicus (He et al., 2014) and amphibious mudskipper Periophthalmus modestus (He et al., 2015) between the ECS and the NSCS. Avise et al. (1987) concluded that mtDNA phylogenetic continuities in the absence of current spatial separation could emerge in species having relatively extensive and recent historical interconnections through gene flow. This appears to hold true for D. maruadsi. After the last glaciation, rising sea levels opened the Taiwan Strait, which provided the opportunity for gene flow among D. maruadsi populations along the southeast coast of China. Moreover, D. maruadsi, especially the ECS population, possesses great dispersal capabilities, either as pelagic larvae, juveniles, and/or adults (Zheng et al., 2003). Adult D. maruadsi exhibit highly migratory connectivity across the species range (e.g., overlapping spawning time and spawning ground, longdistance spawning, overwintering, and preying migration), and the planktonic eggs and larvae could travel great distances on many annual or semi-annual ocean currents (Deng and Zhao, 1991;Niu et al., 2019). For instance, the southward currents (i.e., China Coastal Current and northeast monsoon drift) may carry the ECS sourced individuals toward the NSCS, whereas the northward currents (i.e., southwest monsoon drift and Taiwan Warm Current) may carry the NSCS sourced individuals toward the ECS (Li et al., 2000;Qiao, 2012). Previous studies have suggested that if there are no firm and longstanding physical  barriers between sampled locations, a high level of gene flow is possible for marine species with high dispersal ability, such as mud crab S. paramamosain (He et al., 2010), small yellow croaker L. polyactis (Wu et al., 2012), Chinese beard eel Cirrhimuraena chinensis Kaup , and black sea bream A. schlegelii (Zhao et al., 2021). In this study, we detected a high level of gene flow among 16 populations, confirming that the dispersal ability of D. maruadsi is an important determinant of the no significant genetic structure among all surveyed populations. One more point is that our BSP analysis showed that D. maruadsi recolonized the exposed continental shelf of the ECS and NSCS from glacial refugia with rising sea levels, which illustrated the species did not reach an equilibrium due to the lack of sufficient evolution time (Slatkin, 1993). Given all that, genetic homogeneity of D. maruadsi along the coast of the ECS and NSCS could mainly be explained by a combination of historical biogeography and the dispersal ability of the species.
Although the genetic differentiation between the ECS and NSCS group was not significant, a relatively higher level of genetic diversity was detected in the NSCS group than in the ECS group and the groups exhibited different migratory capacities and population compositions. The ECS group exhibits annual long-distant migration, including spawning migration, feeding migration, and overwintering migration ( Figure 1B) (Zheng et al., 2003). For example, the western ECS population can migrate from the center of Taiwan Strait (24°N) to Jeju Island (30°N) coastal regions (Zheng et al., 2003). However, the NSCS group only shows north-south or deepshallow short-distant migration across the coastal waters ( Figure 1B) (Deng and Zhao, 1991). Furthermore, recently studies have indicated that the dominant body lengths in the ECS group and NSCS group vary, e.g., in the spring, the dominant body lengths for larvae and adults of D. maruadsi in the ECS are 50-110 mm and 200-290 mm (Jiang et al., 2012), while those of the species in the NSCS are 101-110 mm and 181-200 mm (Wang et al., 2021), respectively. Overall, the population genetic structure (genetic homogeneity) of D. maruadsi along the coast of the ECS and NSCS was inconsistent with patterns of genetic diversity and biological characteristics (differences), emphasizing the importance of identifying genetic resources and management units of this species with caution and accounting for this information in the future studies.

Fisheries Management Implications
It is of vital significance to correctly understand the evolutionary history, genetic diversity, population structure, and ecological distribution of a species, in particular of a population with wide distribution, for the management and diversity conservation of important economic fishes (Utter, 1991;Grant et al., 2017). D. maruadsi is widely distributed in the coastal waters of China, and its effective population size may not be small. However, a low levels of genetic diversity was detected in D. maruadsi along the southeast coast of China, especially the ECS population. Moreover, recent fishery resources survey showed that D. maruadsi was once one of the most abundant fishes, but it has suffered a considerable catch reduction in some coastal waters of China, such as Yellow Sea, Bohai, and Beibu Gulf (Chinese Fishery Statistical Yearbook, 1997Geng et al., 2018), possibly due to overfishing, environmental changes, or other factors. Genetic characterization and current status of D. maruadsi resource indicate that its germplasm resources has presented downdraft phenomena. Therefore, in order to ensure the sustainable utilization of fishery resources, D. maruadsi should be closely monitored in the long-term, and conservation and management of its germplasm resources should also be strengthened to protect the genetic diversity of this fish as much as possible. In this study, D. maruadsi along the coast of the ECS and NSCS exhibited no significant genetic structure, indicating that the ECS and NSCS populations could be considered as a management unit for conservation. But we think that the ECS and NSCS D. maruadsi cannot be entirely classified into a management unit, due to the potential mismatch among genetic structure, genetic diversity, and biological characteristics, and weak genetic evolution signal. Given the maternally inherited characteristics of mtDNA, the population structure of D. maruadsi need to further investigate based on higher resolution nuclear genetic markers, such as microsatellites and single nucleotide polymorphisms, which would be of significant help to determine more accurate and refined management units and then design an effective management policy.

Reidentification of D. russelli and D. maruadsi
Many previous studies showed that Cyt b sequences have sufficient resolution and can effectively distinguish many different group of fishes including Carangidae, which is very similar with that of COI and has been used to identify many related species (Reed et al., 2001;Kartavtsev et al., 2007;Bektas and Belduz, 2008;Xia et al., 2008;Khosravi et al., 2020). Based on phylogenetic analysis of Cyt b sequences, we detected two highly divergent clades (I and II), which were monophyletic with high confidence (1.00). The genetic distance (0.0594) between the two clades was about 14-and 7-fold higher than that within clade I (0.0042) and clade II (0.0087), respectively. In fact, the genetic distances based on Cyt b between the Northern Vietnam clade (corresponding to our clade I) and Sundaland-Rosario-Ranong clade (corresponding to our clade II) (0.0290-0.0510) computed by Jamaludin et al. (2020) were also significantly greater than that among populations within clades (0.001-0.003), consistent with the ten-fold rule between species and genera (Ward et al., 2005), which confirmed that each clade was a valid species. The genetic distances between two clades were fairly comparable to that between D. maruadsi and D. russelli (0.0245 and 0.0310) (Hou et al., 2020;Zhang et al., 2020) and the low values reported between species within the family Carangidae (0.0213-0.1680) (Reed et al., 2001;Bektas and Belduz, 2008) or other marine fishes (0.0370-0.1300) (Johns and Avise, 1998), but were higher than the accepted threshold to delimit species (0.02) according to COI barcode sequences (Hebert et al., 2003). Accordingly, we can infer that the two clades defined by Jamaludin et al. (2020) are actually a cluster of two closely related species within the genus Decapterus, rather than significant geographic structure within D. maruadsi. Furthermore, we confirmed that 23 Cyt b haplotypes of the Northern Vietnam clade recorded in Jamaludin et al. (2020) were indeed D. maruadsi because they clustered together with 430 individuals of D. maruadsi collected in our study from the southeast coast of China to form clade I. However, 57 Cyt b haplotypes of the Sundaland-Rosario-Ranong clade defined by Jamaludin et al. (2020) clustered together with four D. russelli haplotypes (outgroups) and were assigned to clade II, indicating that they are D. russelli rather than D. maruadsi.
Both D. russelli and D. maruadsi are economically important pelagic fish in the family Carangidae, and are often mistakenly identified due to the overlapping distributions and similar morphological characteristics. D. russelli and D. maruadsi are widely distributed along the coastal waters of the Indo-West Pacific and western Pacific, respectively, with extensive overlap in the SCS (Smith-Vaniz, 1984;Liu et al., 2016; FishBase, https:// www.fishbase.se/search.php; GBIF, https://www.gbif.org/; TaiBNET, https://taibnet.sinica.edu.tw/). Moreover, most of the countable characteristics, such as dorsal fin, pectoral fin, pelvic fin, anal fin, gill rakes, scute, and vertebrae, largely overlap between D. russelli and D. maruadsi (Table S6) (Cheng and Zheng, 1987;Chen and Zhang, 2015;Wu and Zhong, 2021), making it impossible to distinguish the two species. There are several major morphological diagnostic characteristics to distinguish D. russelli and D. maruadsi, including color of fins, shape of the posterior end of the maxilla, scales before dorsal fin, etc (Table S7). For example, D. maruadsi has a white spot above the anterior part of the second dorsal fin but D. russelli does not; the caudal fins, dorsal fin, and pectoral fin of the latter are yelloworange or scarlet, while those of the former are blue-gray or silver (Table S7) (Cheng and Zheng, 1987;Chen and Zhang, 2015;Wu and Zhong, 2021). However, these different diagnostic characteristics may become difficult to distinguish due to long distance transport and prolonged oxidation, which may result in species misidentification. Collectively, the accurate identification of D. russelli and D. maruadsi based on only morphology is difficult, and future studies should therefore pay more attention to this issue.

CONCLUSIONS
The present study reveals that Late Pleistocene-Holocene climate fluctuations had profound impacts on the population dynamics and phylogeographic structure of D. maruadsi in the ECS and NSCS. With respect to historical dynamics, D. maruadsi experienced recent population expansions due to increasing paleo-productivity and niche space during the LGM to early Holocene. Concerning genetic structure, there was no significant differentiation among populations from the ECS and the NSCS, which probably reflected widespread and recent historical interconnections during the post-glaciation. It is noteworthy that the population genetic structure (genetic homogeneity) of D. maruadsi was inconsistent with genetic diversity and biological characteristics, suggesting that it is necessary to identify genetic resources and management units carefully. However, further studies with larger sample range and more molecular markers are required to resolve the apparent mismatch among population genetic structure, genetic diversity, and biological characteristics, which would be of significant help to determine more accurate and refined population structure and management units. In addition, 57 Cyt b haplotypes of the Sundaland-Rosario-Ranong clade should belong to the species D. russelli rather than D. maruadsi, suggesting that morphological characters are insufficient for accurate identification of the two fishes. Overall, this study will provide useful reference and genetic information for supervision of germplasm resources, formulation of effective management policy, and taxonomic identification and evolutionary relationship research of species in the genus Decapterus.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, OM728655-OM728833.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Experimental Ethics Committee of Guangdong Ocean University, China.