A total of 1,224 specimens of the genus Trichiurus were collected from 17 fishing ports in the China Seas in the fall (September and October) of 2018 with longline, gill, and trawl nets by fishermen ( Figure 1A ; Table 1 ). The localities were classified into three seas. Six sampling localities were established off the coast of the YS: (1) Dondog, DD; (2) Dalian, DL; (3) Weihai, WH; (4) Qingdao, QD; (5) Lianyungang, LY; and (6) Nantong, NT. Four sampling localities were established off the coast of the ECS: (7) Zhoushan, ZS; (8) Wenzhou, WZ; (9) Ningde, ND; and (10) Quanzhou, QZ. Seven sampling localities were established off the coast of the SCS: (11) Shantou, ST; (12) Shenzhen, SZ; (13) Yangjiang, YJ; (14) Zhanjiang; ZJ; (15) Beihai, BH; (16) Wenchang, WC; and (17) Sanya, SY ( Table 1 ; Figure 1A ). The COI gene was amplified by polymerase chain reaction (PCR) using the primers Fish-F2 (5’-ACCTCTGTGTGTGGGGCTAC-3’) and Fish-R2 (5’-GTGATGCATTGGCTTGAAA-3’) (Gu et al., 2021). Each 50 µl PCR mixture contained 5 ng of template DNA, 5 µl of 10× reaction buffer, 4 µl of dNTP mix (10 mM), 5 pmol of each primer and 2 U of Taq polymerase (TaKaRa Biomedical Technology, Beijing). PCR was carried out on an MJ Thermal Cycler with one cycle of denaturation at 94°C for 3 min and 40 cycles of denaturation at 94°C for 30 s, annealing at 51°C for 30 s, and extension at 72°C for 1 min 30 s, followed by a 72°C extension for 10 min and storage at 4°C. The purified PCR products were sequenced using an ABI 377 automated sequencer (Applied Biosystems, Foster City, CA, USA). The chromatograms were checked with Chromas software (Technelysium), and the sequences were manually edited using BioEdit 6.0.7 (Hall, 1999).

The nucleotide sequences were aligned with Clustal X 1.81 (Thompson et al., 1997). The selection of the best-fitting nucleotide substitution models was performed using the Bayesian information criterion (BIC) in jModelTest 2.0 (Darriba et al., 2012). The most appropriate nucleotide substitution model was GTR + I + G. Phylogenetic analysis was used to identify the species in MEGA-X (Kumar et al., 2018). The intrapopulation genetic diversity levels were estimated using haplotype diversity (h) (Nei and Tajima, 1983) and nucleotide diversity (θ_π) indices (Jukes and Cantor, 1969) in DnaSP v5 (Librado and Rozas, 2009). The existence of a phylogeographic structure was examined following the method of Pons and Petit (1996) by calculating the two genetic differentiation indices (G _ST and N _ST) in DnaSP. G _ST, or genetic differentiation among populations, depends only on the frequencies of the haplotypes. N _ST, the degree of substructure, is influenced by both haplotype frequencies and the distances between haplotypes.

The demographic histories of hairtail fish were reconstructed using four different approaches. First, we performed Tajima’s D and Fu’s F_S neutrality tests (Tajima, 1989; Fu, 1997) in DnaSP. Under a population expansion model, significant negative values for Tajima’s D and Fu’s F_S neutrality tests were expected. Statistical significance was tested using 10,000 permutations. Second, the mismatch distribution (Rogers and Harpending, 1992) was estimated under the assumption of sudden expansion as implemented in Arlequin. The sum of squared deviations (SSD) between observed and expected mismatched distributions, the raggedness index (Rg), and the Ramos-Onsins and Rozas (R2) statistic were used as test statistics with 1,000 bootstrap replicates. In the third approach, we reconstructed historical demography using the coalescent-based Bayesian skyline plot (BSP) approach implemented in BEAST (Drummond et al., 2013). A stick clock model with a Bayesian skyline tree was used. We ran 10⁶ generations. The time to the most recent common ancestor (T_MRCA) was also calculated using BEAST. A mutation rate of 0.55% per million years (myr) has been estimated for the mtDNA COI gene in electric fish (Picq et al., 2014). This rate falls within the range of values previously reported for other freshwater fishes (Webb et al., 2004). The burn-in and plots for each analysis were visualized using Tracer v1.6 (Rambaut et al., 2013) to determine whether the convergence and suitable effective sample sizes were achieved for all parameters. This software was also used to generate skyline plots.

Finally, we performed an approximate Bayesian computation (ABC) with DIYABC v.2.0 software (Cornuet et al., 2014) to detect population size changes and then infer demographic parameters. To examine past transitions in effective population size, our study tested six population demographic scenarios. The scenarios were as follows:

In scenario A (the constant model), the effective population size was constant at N1 from the present to the past.

In scenario B (the decline model), the effective population size changed from Na to N1 at t, and Na was larger than N1.

In scenario C (the expansion model), the effective population size changed from Na to N1 at t, and Na was smaller than N1.

In scenario D (the isolation and constant model), the populations were isolated, and the effective population size was constant at N1 from the present to the past.

In scenario E (the isolation and decline model), the populations were isolated, the effective population size changed from Na to N1 at t, and Na was larger than N1.

In scenario F (the isolation and expansion model), the populations were isolated, the effective population size changed from Na to N1 at t, and Na was smaller than N1.

For all analyses, the reference table was built with 1,000,000 simulated datasets per scenario using the following summary statistics: one-sample statistics for the number of haplotypes, Tajima’s D, the mean number of pairwise differences, the variance of pairwise differences, and the number of segregating sites; and two-sample statistics for the mean of within-sample pairwise differences, the mean of between-sample pairwise differences, the number of segregating sites, and F_ST between samples. The posterior probabilities were compared by logistic regression. Additionally, the DIYABC program was used to estimate the population size of the ancestral populations of each species or metapopulation.

Population genetic connectivity was assessed using the Bayesian MCMC method implemented in MIGRATE-N 4.4.3 (Beerli, 2016). MIGRATE-N calculates the mutation-scaled population size (θ = 2Neμ for haploid mtDNA) and immigration rate (M = m/μ) for each area. Within T. japonicus, all sampling populations were separated into three metapopulations based on three oceanographic areas: the YS, ECS, and SCS. Four models of dispersal were first evaluated:

M1 is a full migration model with three population sizes and six immigration rates.

M2 is a source–sink model with three population sizes and three directional north-to-south immigration rates.

M3 is a source–sink model with three population sizes and three directional south-to-north immigration rates.

M4 is an island model where all areas share a single mean estimate of θ and exchange genes with all other areas at the same mean rate.

The models were ranked by log Bayes factors (LBFs), which compare the marginal likelihoods of models calculated using the thermodynamic integration method. The most useful information is found in the model ranked first. The effective number of immigrants per generation was calculated for haploid data with female transmission following the equation Nm = θ × M. The MIGRATE-N program was also used to estimate the emigration and immigration of each population to better understand migration among populations. In these analyses, we ignored some populations that included too few samples, for example, population ZJ of T. nanhaiensis, population BH of T. brevis, and populations ST, ZJ, and WC of T. japonicus ( Table 1 ). Finally, to examine whether the migrations of T. japonicus corresponded to feeding, our study also checked whether the stomach was empty and estimated the ratio of empty stomachs in each population.

A total of 1,224 hairtail fish were collected from 17 fishing ports along the China Seas ( Figure 1A ). Species were first identified by morphology. Two species groups were recognized: the T. lepturus complex, which has an anal opening positioned vertically at the 38th to 41st dorsal fin rays, and the T. russelli complex, which has an anal opening positioned vertically at the 34th to 35th dorsal fin rays ( Table 2 ) (Burhanuddin et al., 2002). Within the T. lepturus complex, T. japonicus has a longer tail, and T. lepturus has a whitish dorsal fin when fresh; in contrast, T. nanhaiensis has a yellowish-green dorsal fin (Hsu et al., 2009); the frontal bone of T. nanhaiensis is very smooth, the frontal bone of T. japonicus is slightly inverted, and the frontal bone of T. lepturus is obviously inverted and bulges in the upper part of the orbit and is accompanied by indentation ( Figures S1, S2 ; Yi et al., 2022). A total of 1,551 bp of the COI gene was sequenced. The nucleotide sequences were A+T rich (53.3%). To identify species, these sequences were analyzed by phylogenetic analysis and BLAST of the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov). The neighbor-joining (NJ) tree showed that these 1,224 COI sequences in the China Seas fell into three phylogroups ( Figure 1B ). Compared with the sequences in GenBank, these three phylogroups were identified as three species belonging to two species complexes: T. nanhaiensis and T. japonicus in the T. lepturus complex and T. brevis in the T. russelli complex. Trichiurus nanhaiensis and T. brevis are found only in the SCS, and the China Seas harbor T. japonicus ( Figure 1A ; Table 1 ).

In total, 127 COI haplotypes from 203 specimens of T. brevis were sampled from six populations in the SCS ( Figure 1A ; Table 1 ). Among the 127 haplotypes, ten were shared by two or more populations. Population SY had the most shared haplotypes ( Table 1 ). The most widespread haplotype was distributed among all six populations. According to N _ST (0.099) and G _ST (0.031), T. brevis displayed no significant phylogeographic structure or population differentiation ( Table 3 ).

In total, 66 COI haplotypes from 119 specimens of T. nanhaiensis were sampled from three populations in the SCS ( Figure 1A ; Table 1 ). Although the sample size and the number of haplotypes of T. nanhaiensis were smaller than those of T. brevis, the nucleotide diversity (θ_π) was larger than that of T. brevis ( Table 1 ). According to N _ST (0.149) and G _ST (0.038), T. nanhaiensis also displayed no significant phylogeographic structure or population differentiation ( Table 3 ). Among the 66 haplotypes, only six were shared by two or more populations. ST and SZ populations shared the most haplotypes (five shared haplotypes).

In total, 521 COI haplotypes from 902 specimens of T. japonicus were sampled from 16 populations in the YS, ECS, and SCS ( Figure 1A ; Table 1 ). Among the 17 sampling populations, only SZ did not contain T. japonicus. The haplotype and nucleotide diversities of T. japonicus were larger than those of the other species ( Table 1 ). Among the 521 haplotypes, 96 were shared by two or more populations. Population QD in the YS had the most private haplotypes, while population ZS in the ECS had the most shared haplotypes ( Table 1 ). The northernmost population (DD) and the southernmost population (SY) had nine shared haplotypes. Two shared haplotypes were distributed widely in twelve and fourteen populations. The populations in the YS, ECS, and SCS had 85, 73, and 45 shared haplotypes, respectively. A total of 60, 39, and 32 haplotypes were shared by the YS and ECS, YS and SCS, and ECS and SCS, respectively. There were 26 haplotypes distributed in the China Seas. According to N _ST (0.000) and G _ST (0.002), T. japonicus displayed no phylogeographic structure or population differentiation ( Table 3 ).

The significant negative values of Tajima’s D and Fu’s F_S supported a population expansion model for these three species ( Table 3 ). The results of the mismatch distribution analysis showed that T. japonicus, T. nanhaiensis, and T. brevis all underwent population expansions, although the distributions of pairwise differences seemed bimodal ( Figure 2 ; Table 3 ). The BSPs showed that all species displayed population expansion at approximately 0.2–0.4 mya ( Figure 3A ). Trichiurus japonicus first slightly expanded at approximately 3.0 mya (data not shown); its population size was stable at 3.0–0.4 mya, and then it expanded again. At approximately 0.25 mya, the population size declined, and then it expanded again to 0.22 mya. Finally, at approximately 0.08 mya, the rate of population expansion declined ( Figure 3A ). Trichiurus nanhaiensis and T. brevis displayed a similar pattern. The populations of these two species both expanded at 0.20 and 0.24 mya, and the rate of population expansion declined at 0.15 and 0.18 mya ( Figure 3A ). Additionally, the population size of T. nanhaiensis has declined recently ( Figure 3A ).

Assigning the populations of T. japonicus to three metapopulations, YS, ECS, and SCS, the significant negative values of Tajima’s D and Fu’s F_S supported a population expansion model for these three metapopulations ( Table 3 ). The results of the mismatch distribution analysis indicated a population expansion, and the distributions of pairwise differences seemed bimodal ( Figure 2 ; Table 3 ). The patterns shown in the BSP among these three seas were slightly different, but still showed similarities ( Figure 3B ). The population underwent two expansion events. First, the population sizes in the YS, ECS, and SCS all expanded by 0.50 mya. Second, the population sizes in the YS, ECS, and SCS all expanded by 0.15–0.23 mya. Finally, the rate of population expansion in the YS and ECS declined by 0.15 and 0.07 mya, respectively, but the rate in the SCS did not decline ( Figure 3B ). Additionally, before the second population expansion event, the YS suffered from declining pattern.

The results of a BEAST analysis suggested that the times to coalescence for T. japonicus, T. nanhaiensis, and T. brevis were sometime in the Late Miocene (T_MRCA = 8.419), Middle Pleistocene (T_MRCA = 1.485), and Middle Pliocene (T_MRCA = 3.988) ( Table 3 ). The times to coalescence for the three metapopulations of T. japonicus were 4.863 (YS), 4.305 (ECS), and 8.145 (SCS).

The DIYABC results differed among these three species. In T. nanhaiensis, the highest posterior probability was found for the isolation and constant scenario (Scenario D). Its posterior probability (0.8908, 95% CI: 0.6205–1.0000) was higher than that of other scenarios ( Table 4 ). In T. brevis, the highest posterior probability was also found for Scenario D. Its posterior probability was the highest (0.9976, 95% CI: 0.9791–1.0000). In T. japonicus, the highest posterior probability was found for the isolation and expansion scenario (Scenario F). Its posterior probability was the highest (1.0000, 95% CI: 0.9999–1.0000) ( Table 4 ). In metapopulations of T. japonicus, the YS and ECS displayed the highest posterior probabilities for the isolation and constant scenarios (Scenario D, 1.0000 and 0.9931, respectively), and the highest posterior probability in the SCS was observed for the isolation and decline scenario (Scenario E; 0.9992, 0.9947–1.000; Table 4 ).

Migration estimates suggest that T. japonicus complies with a full migration model among the China Seas (M1; Figure 4A ). M1 had the lowest log marginal likelihood of the four gene flow models tested (LBF = 0.00). Migrations among populations were also estimated ( Table S1 ). The results of the regression analysis displayed no significant relationship between immigration and emigration in each population (y = 0.1186x + 6.2647, R² = 0.0064, where x = emigration and y = immigration; Figure 4B ). In each population, populations WZ, ZS, and SY had the most immigration, and populations WH, ZS, and ND had the most emigrations ( Figure 5 ). The values of emigration minus immigration in each population showed that most individuals did not stay in populations WH, ND, and YJ, and most individuals tended to stay in populations WZ, QZ, and SY ( Figure 5 ).

In total, our study found that T. japonicus in the YS migrated southward to the ECS, and a few specimens migrated southward to the SCS ( Figure 6 ; Table S1 ). Additionally, our study also found that some specimens of T. japonicus migrated from the SCS northward to the ECS ( Figure 6 ). Furthermore, the results of the MIGRATE analysis showed that most specimens of T. nanhaiensis and T. brevis migrated northward from the SCS to the Taiwan Strait (population ST), but few specimens migrated southward along the coast ( Figure 6 ).

To determine whether the migrations of T. japonicus corresponded to feeding, our study also estimated the ratio of empty stomachs in each population ( Table 1 ). Our study found that the empty stomach ratio in most populations is less than 0.50. The population DL in YS displayed the highest empty stomach ratio (0.82), and the population QZ in ECS displayed the lowest empty stomach ratio (0.17). The mean empty stomach ratios in the three seas were 0.52 (YS), 0.36 (ECS), and 0.53 (SCS). The regression analysis displayed no significant relationship between immigration and emigration and the ratio of empty stomachs (R² = 0.0067 and 0.0038, data not shown). These results might be because we only determined whether the stomach was empty in each specimen. Moreover, our study divided that data into two groups, with ratios of empty stomachs of more than 0.50 and less than 0.50. In the group with empty stomach ratios of less than 0.50, the regression analysis displayed no significant relationship between immigration and emigration and the ratio of empty stomachs (R² = 0.0934 and 0.0056, Figure 7 ). In the group with a ratio of empty stomachs greater than 0.50, the regression analysis displayed a slightly significant relationship between immigration and emigration and the ratio of empty stomachs (R² = 0.1182 and 0.3449, Figure 7 ). When our study estimated immigration, emigration, and the ratio of empty stomachs in three seas, the regression analysis displayed a strong significant relationship between immigration and emigration and the ratio of empty stomachs (R² = 0.5207 and 0.9710, data not shown), although the sample size was only three (YS, ECS, and SCS).

A species is a necessary unit in all biological, ecological, and genetic studies. Clear species identification for fishery management purposes has been the subject of a major FAO program since the 1960s. However, the taxonomic status of the hairtail has not been resolved, especially within the genus Trichiurus. The species within the genus Trichiurus could be assigned to two groups: large-headed and short-tailed fishes (Hsu et al., 2009). The large-head group is also referred to as the T. lepturus species complex. Hsu et al. (2009) found three species within this species complex: T. japonicus, T. nanhaiensis, and T. lepturus. T. japonicus was distributed along the continental shelves in the northeastern Pacific, especially in the seas around Japan and mainland China. T. nanhaiensis was primarily distributed in the SCS, and T. lepturus was distributed in the Pacific and Atlantic. According to Hsu et al. (2009), only T. japonicus and T. nanhaiensis were distributed along the coast of the Asian continent. The findings of this study ( Figure 1A ; Table 1 ) and He et al. (2014) supported these distribution patterns. Although T. japonicus is distributed throughout all the China Seas, it is much less abundant in the SCS. The findings of this study revealed that T. brevis is more abundant in the SCS than other Trichiurus species ( Figure 1A ; Table 1 ).

The sea temperature in the SCS is higher than that in the ECS and YS. Accordingly, our study suggests that these three Trichiurus species have different habitat preferences along the coast of China. T. japonicus prefers colder waters on the continental shelves, while T. nanhaiensis prefers warmer waters on continental slopes and T. brevis prefers warmer waters on the continental shelves. These results corresponded to the results of Hou et al. (2021). Hou et al. (2021) found that the spawning of T. japonicus occurs along the continental shelf, T. nanhaiensis occurs along the continental shelf to the slope, and T. brevis occurs in shallow water. Although our study found that T. nanhaiensis and T. brevis migrated from the SCS to the southern end of the Taiwan Strait (see below: Migration for Feeding and Spawning), the Taiwan Strait is a barrier to the distribution patterns of other species (e.g., Liu et al., 2007; Xu et al., 2009; Liu et al., 2011; Hsu et al., 2021). The spatial and temporal distributions of Trichiurus species are determined by hydrological conditions, e.g., salinity and temperature (Hou et al., 2021).

Among these three Trichiurus species along the Chinese coastline, the genetic diversity (θ_π) of T. japonicus was higher than that of T. nanhaiensis and T. brevis, and the latter two species had similar genetic diversity ( Table 1 ). However, the T_MRCA of T. japonicus was higher than that of T. nanhaiensis and T. brevis, and that of T. brevis was higher than that of T. nanhaiensis ( Table 3 ). Genetic diversity did not correspond to the TMRCA values. Yi et al. (2021a) found that the size of the specimen distribution of a given species and its genetic diversity are positively related based on the family Apogonidae in the SCS. Our study found that T. nanhaiensis and T. brevis are found only in the SCS, and T. japonicus is found in all the China Seas ( Figure 1A ; Table 1 ). Thus, our study suggests that the genetic diversity and the size of the specimen distribution area have positive relationships.

The results of N _ST and G _ST showed that these three Trichiurus species lacked population structure ( Table 3 ). The distribution of the shared haplotypes also suggested that gene flows among populations were frequent. However, our study found that the results of the mismatch distribution analysis showed that the distributions of pairwise differences in these three Trichiurus species seemed bimodal ( Figure 2 ). Previous studies reported that the bimodal mismatch distribution reveals population divergences (Yan et al., 2020; Yi et al., 2021b). Some studies found that there were two spawning cohorts of T. japonicus in the ECS and SCS, from February–July and November–December (Sun et al., 2020a; Sun et al., 2020b). The results of the DIYABC analyses also showed that the populations of Trichiurus species met in isolation ( Table 4 ).

Additionally, the results of the T_MRCA displayed a strange pattern ( Table 3 ). The T_MRCA of the total of T. japonicus and that of the metapopulation in the SCS were close (8.419 and 8.15), but those in the metapopulations in the YS and ECS were 4.863 and 4.305, respectively. However, T. japonicus was mainly distributed in the YS and ECS, and the population sizes in the YS and ECS were higher than those in the SCS ( Figure 3B ). Our study suggests that during glaciations, YS and ECS were completely exposed, and the SCS became a refugium. Thus, the T_MRCA in the SCS and that of the total samples were similar and higher than those in the YS and ECS ( Table 3 ). This result also implied that some populations of T. japonicus retreated to the SCS during glaciations and settled in interglacial periods. Thus, although our study found no significant genetic differentiation of these three Trichiurus species along the Chinese coastline because their swimming abilities were strong, we considered that Trichiurus species might have been differentiated in the China Seas and were shaped by landform changes during glaciations, as in other previous phylogeographic studies (e.g., Shen et al., 2011; Yang et al., 2017; Qu et al., 2018). Sea level and temperature fluctuations promoted differentiation in the Northwestern Pacific Ocean (Liu et al., 2006; Shen et al., 2011; Qu et al., 2018).

Jaureguizar et al. (2004) proposed that the major aims of fish migrations are seasonal and reproductive. Sun et al. (2020a) also proposed that T. japonicus undergoes seasonal migration. Our study sampled specimens in the late fall and early winter. In this season, marine species might migrate southward, following warm waters through the China Coastal Current ( Figure 1 ). Our study found that these three Trichiurus species generally did not migrate southward through the China Coastal Current ( Figures 1 , 6 ). The migration pattern of T. japonicus was supported by the full migration model ( Figure 4A ). Regression analysis revealed no significant positive relationship between immigration and emigration in each population (y = 0.1186x + 6.2647, R² = 0.0064, where x = emigration and y = immigration; Figure 4B ). These results might have resulted in their strong swimming abilities.

According to the migrations between populations, our study inferred the migration routes of T. japonicus and found that most specimens migrated from the YS southward to the ECS and SCS ( Figure 6 ). Most specimens remain in the ECS, including the Taiwan Strait. Moreover, our study also found that some specimens migrated northward from the SCS to the Taiwan Strait. Our study also revealed that T. nanhaiensis and T. brevis migrated from the SCS to the Taiwan Strait ( Figure 6 ). This region, south of the Yangtze River and north of the Taiwan Strait, has been recognized as the wintering ground of T. japonicus (Wang et al., 2011; Sun et al., 2020a; Sun et al., 2020b). In this region, the waters are formed by the Kuroshio Current and diluted Yangtze waters. The former current has high salinity and temperature, and the latter has low salinity and abundant nutrients (Chen et al., 2009). The branch of the Kuroshio Current flows from the east of the Philippines through the Bashi Channel into the Taiwan Strait. Moreover, Han et al. (2013) also found that nutrient transport from the ECS is a major source of support for winter primary production on the northeastern SCS shelf through the Taiwan Strait. Thus, there is more biomass in this area than in the surrounding areas. The northern wintering ground extended to the Taiwan Strait and northeastern SCS. Jaureguizar et al. (2004) suggest that both seasonal migrations and reproductive movements are affected by temporal structure. Chiang et al. (2002) also proposed that water temperature is a limiting factor for phytoplankton biomass density, and that carbon biomass increases in the warm season. Accordingly, our study suggests that these migrations of the Trichiurus species were shaped by feeding ( Figure 7 ). Our study suggested that the ECS, including the Taiwan Strait, is a wintering ground for T. japonicus, and T. nanhaiensis and T. brevis winter in the northeastern SCS ( Figure 6 ).

In addition to the wintering grounds in the ECS, including the Taiwan Strait, our study found that some specimens of T. japonicus remained in population SY ( Figure 5 ). The water temperature in the SCS is higher than that in the YS and ECS. Hou et al. (2021) found that in late summer and early autumn, eggs of T. japonicus and T. nanhaiensis mainly occurred in the surface water mass of the SCS. Thus, our study suggests that population SY is a wintering and spawning ground for T. japonicus ( Figure 6 ). These results also supported the results of the T_MRCA ( Table 3 ). During the LGM, T. japonicus migrated from the YS and ECS to the SCS as a refugium, and some populations settled during the interglacial period. Accordingly, our study suggests that the migrations of Trichiurus species are mainly associated with feeding and spawning preferences and that the behavior is shaped by oceanographic conditions (e.g., Han et al., 2013; Liu et al., 2020; Hou et al., 2021).

Many studies of marine species have revealed a population expansion before or after the Last Glacial Maximum (pre-LGM or post-LGM; Delrieu-Trottin et al., 2017; Jenkins et al., 2018; Qu et al., 2018). Jenkins et al. (2018) also proposed that some species display stable population sizes pre- and post-LGM. Clark et al. (2009) identified the interval with LGM as between 0.027 and 0.019 mya. Our study found that these three Trichiurus species all showed a population expansion after the LGM ( Figure 3A ). Previous studies also displayed similar patterns (He et al., 2014; Lin et al., 2021). However, our results and those of He et al. (2014) displayed two different patterns. First, in the results of He et al. (2014), T. japonicus in the ECS experienced a genetic bottleneck, but in our results, T. japonicus in the YS experienced a genetic bottleneck, and that in the ECS did not ( Figure 3B ). Second, in the results of He et al. (2014), the population of T. japonicus in the SCS declined, but in our results, the population expanded ( Figure 3B ). Our study considered that these differences might have resulted from the sampling season. He et al. (2014) might have sampled in winter, while our study samples were collected in late fall and early winter. In winter, T. japonicus migrated from the YS to the ECS and SCS. Thus, the patterns in He et al. (2014) shifted from north to south.

The neutrality tests, mismatch distribution, and coalescent BSPs supported population expansion ( Table 3 ; Figures 2 , 3 ), but the approximate Bayesian computational approaches supported population declines or remained constant ( Table 4 ). The incongruent results between approximate Bayesian computational approaches ( Table 4 ) and other demographic analyses ( Table 3 ; Figures 2 , 3 ) might have resulted in a population differentiation (Heller et al., 2013; Yan et al., 2020; Yi et al., 2021b). Heller et al. (2013) proposed that the results of BSP analyses show false signals of population decline under biologically plausible combinations of population structure. Yan et al. (2020) and Yi et al. (2021b) also found that when populations met differentiation, the results of the BSP analyses displayed population expansion, but other analyses, e.g., mismatch distribution and Tajima’s D test, did not support the population expansion. Furthermore, the BSPs revealed that the population sizes of three Trichiurus species and three metapopulations of T. japonicus all displayed population expansion after the LGM and experienced a genetic bottleneck during the LGM ( Figure 3 ). However, our study did not find that during the last interglacial period, the population size was not higher than that during glaciation in BSP ( Figure 3 ; He et al., 2014). In other words, the population size was stable in the glacial and the last interglacial periods. Grant et al. (2012) found that the BSPs of simulated sequences showed a flat demographic curve before the LGM and suggested that the results of the BSP analyses lost the pre-LGM population history. Accordingly, our results also support limits in BSP analyses (Grant et al., 2012; Heller et al., 2013). Thus, the results of BSP analyses may need to be evaluated carefully.

Our study suggests that the Trichiurus species underwent population differentiation due to landform changes during the Pliocene and Pleistocene glacial cycles, as in previous phylogeographic studies (Kizaki and Oshiro, 1980; Ota, 1994; Wang, 1999; Voris, 2000; Yang et al., 2017; Qu et al., 2018). The complex geological history might cause incorrect assessment in BSP analyses. However, our study did find evidence supporting the population history of Trichiurus species based on all the results.

Among the three Trichiurus species, the BSP results revealed that the population size of T. nanhaiensis was lower than that of the other two species, and only this species recently displayed population decline ( Figure 3A ). However, the genetic diversity of T. nanhaiensis was not lower than that of the other species ( Table 1 ). Previous studies found that population differentiation could result in overestimated genetic diversity (Yan et al., 2020; Yi et al., 2021b). Our study found that the level of population differentiation of T. nanhaiensis was higher than that of the other two species (N _ST; Table 3 ). Thus, our study showed that the BSP results and genetic diversity in T. nanhaiensis might be overestimated, and the population size of T. nanhaiensis seriously declined ( Figure 3A ). Our study also found that people ate T. nanhaiensis more frequently in the South Sea (our observations) because this species is meatier than T. japonicus and T. brevis.

The BSP results also revealed different patterns among the three metapopulations of T. japonicus ( Figure 3B ). The estimated population size in the ECS was the highest, and that in the SCS was the lowest. These results corresponded to the annual marine capture fishery production of China (Kang et al., 2018). Thus, our study suggests that T. japonicus can be used as an indicator organism to monitor fishery resources in the China Seas. The ECS is the largest marginal sea in the northwestern Pacific, and there are more nutrient sources here than in other areas, including those from the Kuroshio Current and Yangtze River runoff (Chiang et al., 2002). Additionally, in our sampling season, T. japonicus migrated from the YS to the ECS as a wintering ground. Thus, the population size in ECS was higher. The YS is a semi-enclosed body of water, but the coastline is heavily populated, urbanized, and industrialized. Ji et al. (2019) proposed that the largehead hairtail in YS, including the Bohai Sea, experienced overexploitation. Moreover, in our sampling season, the populations in the YS migrated southward. Thus, the YS metapopulation declined significantly ( Figure 3B ). Furthermore, the population size in the SCS increased according to the BSP results ( Figure 3B ). Liu (2013) proposed that among the marginal seas of China, the South China Sea has the highest biodiversity. Our study also inferred a wintering and spawning ground in the SCS, and some populations migrated from the YS and ECS southward to the SCS ( Figure 6 ). As the BSP analyses have limitations, the BSP results would display false signals of population decline under combinations of population structure. Thus, although the BSP results showed population growth in the SCS, the SCS population size was the lowest ( Figure 3B ). Kang et al. (2018) and Du et al. (2020) found that resources experienced a decline in catch per unit effort (CPUE), the annual catch in the YS decreased after 2008, the annual increments declined in the ECS and SCS, and a rapid decrease in the annual capture of the SCS began in 2018. Thus, our study suggests that policies are needed for the resource management of Trichiurus species.