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

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.

Total genomic DNA was extracted from muscle tissues using a standard phenol-chloroform extraction protocol (Huang et al., 2002). PCRs were implemented to amplify partial fragments of mitochondrial Cyt b and CR sequences using previously described primers, L14724 (5′-GAC TTG AAA AAC CAC CGT-3′) and H15915 (5′-CTC CGA TCT CCG GAT TTAC AAG AC-3′) for Cyt b (Xiao et al., 2001) and LF (5′-AAC TCC CAA AGC TAG GAT TCT-3′) and LR (5′-TTT GTG CTT GCG GGG CTT T-3′) for CR (Niu et al., 2012). PCR amplification was carried out in a reaction volume of 50 μL, including 39.5 μL of ddH₂O, 5 μL of 10× EasyTaq buffer (containing Mg²⁺), 200 nM each primer, 200 μM dNTPs, 2.5 U of Taq DNA polymerase (TransGen Biotech Co., Ltd., Beijing, China), and 40–100 ng of template DNA. PCRs of Cyt b and CR were performed according to our previously described procedures (Niu et al., 2018; Niu et al., 2019). Each PCR amplification contained a negative control, which included all reagents except the template DNA. PCR products were separated by electrophoresis on a 1% agarose gel along with the Trans 2K DNA Marker (TransGen Biotech Co., Ltd.) for identifying the target gene size and sent to Sangon Biotech Co. Ltd. (Shanghai, China) for sequencing.

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.

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 (π), 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).

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.

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 τ gained from the mismatch distribution was used to estimate the time of the population expansion with the equation t = τ/2u, where u represents a neutral mutation rate of the entire sequence in each generation calculated as u = 2μk, where μ 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.

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) (Zou et al., 2016), and Decapterus macrosoma (GenBank accession number: KF841444) (Li et al., 2016) 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.

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 π 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 π 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.

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.

The TCS networks constructed using Cyt b, CR, and concatenated sequences (Cyt b + CR) all showed star-like topologies structure ( Figure 3 ), which contained three (CH1, 15.12%; CH3, 66.98%; and CH6, 4.19%), five (CR1, 7.91%; CR2, 10.93%; CR3, 8.37%; CR5, 19.07%; and CR28, 3.72%), and five (H1, 6.51%; H3, 7.21%; H5, 16.05%; H8, 9.07%; and H55, 2.23%) dominant haplotypes at high frequencies, respectively. These high-frequency haplotypes were distributed throughout most of the sampling area and were located in the centre of the haplotype networks. All other low-frequency haplotypes were separated from these dominant haplotypes by only one or low number of mutational steps. Overall, phylogenetic analysis further indicated that there was no obvious genetic structure but a highly diverse within the species D. maruadsi across the southeast coast of China.

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 τ 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.

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, KX212068–KX212069, and KX212072–KX212076). Interestingly, within clade II, 57 Cyt b haplotypes from nine Sundaland populations and Khanh Hoa population (green mark, GenBank accession numbers: KX212003, KX212005–KX212050, KX212059, KX212066–KX212067, KX212070–KX212071, KX212077–KX212079, and KY440109–KY440110) previously identified as D. maruadsi clustered together with four D. russelli haplotypes (outgroups) from the NSCS and the coastal waters of India, Omani Sohar, and Indonesian Ambon (brown mark, GenBank accession numbers: OM479430, MN711693, KX512731, and LC646650, respectively). The net genetic distances based on Cyt b were as follows: clade I/II, 0.0594; clade I/D. macarellus, 0.1420; clade I/D. macrosoma, 0.1945; clade II/D. macarellus, 0.1256; clade II/D. macrosoma, 0.1573; and D. macarellus/D. macrosoma, 0.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.

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 mid-water 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 (Lambeck et al., 2002; 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 paleo-productivity 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.

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 (Lambeck et al., 2002; 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, long-distance 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 (Li et al., 2014), 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 deep-shallow 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.

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, 1997–2021; Geng 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.

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 yellow-orange 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.