Diversifying of Two Pampus Species across the Indo–Paciﬁc Barrier and the Origin of the Genus

: Among marine species distributed in the Indian Ocean and the Paciﬁc, the Indo–Paciﬁc Barrier (IPB) has been found to be an important barrier of divergence of species distributed on both sides. Among the ﬁve species of the genus Pampus , only Pampus chinensis and P. cinereus are distributed across the western Paciﬁc and the Indian Oceans and have not been studied comparatively using extensive sampling and gene markers. Furthermore, the origin and history of genus Pampus remain unrevealed. We used thousands of nuclear loci based on target gene enrichment to explore genetic structure of P. chinensis and P. cinereus across the western Paciﬁc and Indian Oceans. We performed divergence dating and ancestral area reconstruction analysis and inferred the dispersal routes of the Pampus species. The results suggest that the IPB played an important role in the differentiation between populations among the two oceans for both P. chinensis and P. cinereus , dividing species into the Paciﬁc lineage and the Indian Ocean lineage. Low sea level in the late Pleistocene may be the main cause. The result also showed that the South China Sea was the center of origin of the genus Pampus , and dispersal routes of each species may be associated with the ocean currents. Our study provided new examples for the IPB effect on marine species.


Introduction
The Ocean, accounting for about 71% of the earth's surface area, has a total area of 360 million km 2 , which can be divided into a series of biogeographical regions, with the Indo-Polynesian province as the broadest one. In the tropical Indo-Pacific region, six biogeographical provinces (Western Indian Ocean, Red Sea, Indo-Polynesian, Hawaiian, Marquesas, Easter Island) are defined by 10% endemism [1]. The Indo-Pacific Barrier (IPB), which is located in the center of the Indo-Polynesian province, represented by the Indo-Malaysian archipelago, is considered as the major block to gene flow for marine species between the Indian Ocean and the western Pacific [2]. There are three major patterns reported: (1) sister species divided by the IPB, such as Chaetodon trifasciatus distributed in the Indian Ocean and C. lunulatus in the Pacific [3] and Coris cuvieri distributed in the Indian Ocean and C. gaimard in the Pacific [4]; (2) two separated lineages of the same species, such as the reef fish Cephalopholis argus [5], the bigscale soldierfish Myripristis berndti [6], the snapper Lutjanus fulvus [7] and the scalloped hammerhead shark Sphyrna lewini [8]; and (3) no population structure or ongoing high gene flow between the populations distributed in the two oceans, such as the moray eels [9], unicornfishes [10], some groupers [11] and goby Smilosicyopus fehlmanni [12]. Different patterns of genetic structure among different species might be due to different life history characters of the organisms, such as dispersal capability [9], spawning aggregation [11] or population history [12]. Previous phylogeographic studies focused on the coastal fishes along the Northwest Pacific, Bay of Bengal and the Arabian Sea have been seen on some marine species, such as the seahorse [13], the Chinese four-eyed sleeper [14], etc.; however, the species of genus Pampus have not been studies.
The genus Pampus (Perciformes: Stromateidae) comprises five species [15], that is, P. chinensis, P. cinereus, P. minor, P. argenteus and P. punctatissimus, all of which are pelagic fishes and multiple batch spawners [16,17]. Only two species of this genus, P. chinensis and P. cinereus, are widespread across the Indian Ocean and the Northwest Pacific [18,19]. Samples of P. cinereus collected from the Bay of Bengal have been reported to have close affinity to the P. cinereus in the South China Sea but are genetically distinct [20,21]. Pampus chinensis also was reported with differentiation between the Indian Ocean and the western Pacific based on the mitochondrial cytochrome c oxidase I (COI) gene [22]. This differentiation between populations may be influenced by the IPB. Most of these Pampus studies used mtDNA genes, which reflect only the history of the maternal lineage instead of biparental inheritance. Moreover, the origin of genus Pampus and the dispersal history of each species remain untested.
In this study, we collected thousands of nuclear loci based on a target-gene enrichment method [23] to (1) test how the IPB has affected the two cross-ocean species P. chinensis and P. cinereus; (2) examine whether P. chinensis and P. cinereus have diverged into two lineages; and (3) reveal the origin and dispersal routes of species of the genus Pampus.

Laboratory Protocols
Right pectoral fins or muscle tissue were sampled and preserved in absolute ethanol. Genomic DNA was extracted using the Ezup Column Animal Genomic DNA Purification Kit (Sangon, Shanghai, China). The quantity of purified DNA was measured using a NanoDrop 3300 Fluorospectrometric (Thermo Fisher Scientific, Wilmington, DE, USA).
Agarose gel electrophoresis was used to visualize the purified DNA. Then, 1000~1500 ng of purified DNA of each sample was used to prepare an Illumina pair-end library following the methods of Meyer and Kircher [24] and Fisher et al. [25] with some modifications. Then, a subset of purified DNA of was used for target-gene enrichment [23] with multi-120 bp biotinylated RNA baits, which were designed for a total of 17,675 single-copy conserved nuclear coding sequences of Oreochromis niloticus [26]. The products after gene enrichment were sequenced on an Illumina HiSeq X10 platform at Novogene, Inc. (Beijing, China).

Single Nucleotide Polymorphism Calling
Single nucleotide polymorphism (SNP) calling was performed within P. chinensis and P. cinereus, respectively. First, consensus sequences for each gene were generated as reference sequences using consensus.pl [31]. Then, BWA v0.7.16 [32] and SAMTOOLS v1.10 [33] were used to make an index for the consensus sequences with default parameters, and the trimmed reads were mapped to the reference sequences using BWA-MEM [34]. GATK v4.1.8.1 [35] was used to perform the SNP calling and quality filtering.

Population Structure
The SNP sites derived from Hardy-Weinberg equilibrium (HWE, p-value threshold 0.05) were excluded using VCFTOOLS v0.1.17 [36]. Briefly, individuals were grouped corresponding to their sampling ocean, the Indian Ocean and Pacific, respectively. Loci which were derived from the HWE for each group were recorded and filtered out. A custom Perl script was used to filter out low-coverage SNPs and to select only one SNP per targeted locus to avoid the problem of linkage disequilibrium (vcftosnps.pl). This set of SNPs was used as input for STRUCTURE analysis [37] to infer the population structure clustering. Parameters included 200,000 MCMC replicates and 20,000 burnin, groups (K) from 1 to 5, with 10 iterations. The results were compressed and uploaded to STRUCTURE HARVEST [38]. CLUMPP v1.1.2 [39] was used to find the optimal alignments among the 10 replicates of cluster analysis, and DISTRUCT v1.1 [40] was used for plotting. Moreover, another set of independent SNPs with no more than 20% missing data were used to infer population structure through the ADEGENET R package [41] using its implementation of discriminant analysis of principal components (DAPC). First, the xvalDapc function was used to perform cross-validation for deciding the numbers of principal components (PCs) to be retained based on the lowest root mean squared error (RMSE). Second, we chose the number of clusters with lowest Bayesian Information Criterion (BIC) and performed the dapc function with PCs chosen in the first step.
The mean pairwise Analysis of Molecular Variance (AMOVA) F ST and the mean nucleotide diversity (Π) for each population were calculated between populations using the populations program implemented in STACKS v2.59 [42,43]. Then the populations were grouped according to their distribution relative to ocean or sea areas, and the AMOVA F ST was calculated between groups.

Phylogenetic Analysis
The aligned DNA sequences of P. chinensis and P. cinereus were used to reconstruct the phylogenetic tree to examine their divergence. All sequences for P. chinensis and P. cinereus, one P. punctatissimus and two outgroup taxa (Peprilus medius and Psenopsis anomala), were checked for molecular clocklikeness to exclude any sites which did not pass the test using a custom Perl script [44]. A custom Perl script (decidata.pl) was used to retain the loci that contained at least one individual for each sampling site. RAXML v8.2.10 [45] was used to reconstruct the concatenate-based phylogenetic tree, using the rapid bootstrap analysis, a best-scoring maximum likelihood tree search method (option "-f a") and a GTR-GAMMA model with 1000 bootstrap replicates. Coalescent-based phylogenetic tree reconstruction was performed using ASTRAL-III v5.7.5 [46,47]. First, IQ-TREE v2.1.2 [48] was used to reconstruct gene trees with the best-fit model found in ModelFinder [49] ("-m MFP") and 1000 ultrafast bootstrap replicates [50]. Then, these gene trees were used as input for ASTRAL-III v5.7.5. The resulting trees were visualized using FIGTREE v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/, 25 November 2019).

Divergence Time Dating and Demographic History
To reduce the calculation time, up to three individuals per collection site were used, and loci containing at least 80% of these individuals were selected for divergence time dating using BEAST v2.6.3 [51,52]. The GTR+G+I nucleotide substitution model was applied, and the strict clock model was used for estimating the clock rate. The calibration point was the common ancestor of the genus Pampus and Peprilus that was set to within the interval of 38.5 million years ago (Mya) to 57.8 Mya [53][54][55]. The speciation rate was set to around 0.05 per million year [54], a log-normal distribution was employed and other priors were set as default. The MCMC chain length was set to 20 million with 2 million preburnin and logged every 2000 generations. TRACER and TREEANNOTATOR in the BEAST package were used to analyze the results and summarize time-trees with 10% burnin, and then the maximum lineage credibility (MCC) tree was visualized using FIGTREE v1.4.4.
For estimating the demographic history of P. chinensis and P. cinereus, SNPs as well as sequence data were used. Individuals were grouped according to the STRUCTURE and ADEGENET analysis results, and the concatenated sequences were used to perform Bayesian Skyline plot [56], implemented in BEAST with 10 million MCMC chains. The Tajima's D test and Fu's Fs test were calculated with 10,000 replicates using ARLEQUIN v3.5.2.2 [57], while a significant negative value of Tajima's D test or Fu's Fs test indicates recent population expansion after bottleneck.
One SNP per locus was selected by a custom Python script (select_snp.py). This subset of SNPs in VCF format were converted to ARLEQUIN format using PGDSPIDER v2.1.1.5 [58]. ARLEQUIN v3.5.2.2 was used to generate the folded site frequency spectrum (SFS) and 100 bootstrapped folded SFS. The folded SFS was used as input for FAST-SIMCOAL v2.7 [59,60], which simulated different evolutionary scenarios and estimated parameters of predefined models. We assumed 6 models for P. chinensis and 8 models for P. cinereus, and the best models for P. chinensis and P. cinereus were chosen according to the likelihood distribution and Akaike information criterion (AIC). In detail, the models for P. chinensis are: (1) assuming that the eastern Indian Ocean group (containing samples from the Bay of Bengal, BB) is derived from the Pacific group (containing samples from the western Pacific, PAC) and with migration between them; (2) same as model 1 but without migration; (3) assuming that the PAC was derived from the BB and with migration between them; (4) same as model 3 but without migration; (5) assuming that two groups were diverged from a common ancestor that lived in the western Pacific and migration persisted between them; (6) same as model 5 but without migration. Models for P. cinereus are: (1) assuming that the PAC was the ancestral group, that the BB was derived from the PAC, then the western Indian Ocean group (containing samples from the Arabian Sea, AS) was derived from the BB, and without migration between groups; (2) same as model 1, but migration persisted between PAC and BB, AS and BB; (3) assuming that the AS was the ancestral group, BB was derived from the AS, and the PAC was derived from the BB and without migration between groups; (4) same as model 3, but migration persisted between PAC and BB, AS and BB; (5) assuming that the PAC and BB diverged from a common ancestor that lived in the western Pacific, the AS was derived from the AS and without migration between groups; (6) same as model 5, but migration persisted between PAC and BB, AS and BB; (7) assuming that AS and BB diverged from a common ancestor that lived in eastern Indian Ocean, the PAC was derived from the BB and without migration between groups; (8) same as model 7, but migration persisted between PAC and BB, AS and BB. In total, 100 independent runs were performed for point estimation for each model, including 100 k coalescent simulation and 30 ECM cycles ("-L 30"), and the mutation rate was ignored ("-0 option"). The confidence intervals (CI) for the point estimations from the best models were performed with 10 runs per bootstrapped SFS, including parameter set the same as the point estimation.

Infer Ancestral Area of Genus Pampus
Two thousand trees output from BEAST were randomly chosen and used for RASP v4.2 [61] to reconstruct ancestral areas for the genus Pampus. Five states were defined which corresponded to the sea area of sample collection, including the Arabian Sea, the Bay of Bengal, the South China Sea, the East China Sea, and the Bohai Sea. The RASP4 built-in R package BIOGEOBEARS [62] was used to select the historical biogeography model including jump dispersal; then the ancestral area reconstruction was performed with the best-fitting model and 2 max areas included.

Results
An average of 8429 genes were enriched for the respective 99 Pampus individuals, with a minimum capture of 1640 genes and maximum captured of 13,565 genes (Table S1). After filtering the SNP calling results using the GATK pipeline, 22,020 and 65,871 SNPs were obtained for P. chinensis and P. cinereus, respectively.

Population Structure
For P. chinensis, 6579 SNPs were used in STRUCTURE analysis after HWE testing and filtering to retain one SNP per locus. The highest delta K value was obtained when K equaled 2 ( Figure S1a), which suggested that clustering into two groups was the best-fitting model to the observed data ( Figure 2a). The result of DAPC analysis using 5282 SNPs was the same as STRUCTURE result in that the lowest BIC appeared when clustered into 2 groups ( Figure S2a), that is, INO, BAN and SM in one group and SS and TS in the other (Figure 2b). The Π values of each population ranged from 0.102 to 0.141 (Table 2), with BAN having the lowest Π values. The mean pairwise F ST values between populations ranged from 0.041 to 0.206 (Table 3). The mean pairwise F ST values between populations within each ocean were smaller than cross oceans. The highest mean pairwise F ST values was observed between BAN and SS with 0.206. After grouping the populations according to their ocean area, the pairwise F ST values between the Pacific group and the Indian Ocean group was 0.131, indicating high genetic differentiation (Table 3).  For P. cinereus, 9488 SNPs were used in STRUCTURE analysis. The highest delta K value was obtained when K equaled 2 ( Figure S1b), suggesting that two clusters is the best supported model (Figure 3a). The result of DAPC analysis using 7027 SNPs was different from the STRUCTURE result in that the lowest BIC appeared when the data were clustered into three groups ( Figure S2b), counting IRAN as one group, INO and BAN as one group and SS and ES as another group (Figure 3b). The Π values for each population ranged from 0.049 to 0.103 (Table 2), with ES having the lowest Π values. The mean pairwise F ST values between populations ranged from 0.020 to 0.186 (Table 3). Pairwise F ST values between populations within oceans were smaller than cross oceans. The highest mean pairwise F ST value was observed between IRAN and ES s 0.186. After grouping the populations according to their region, the pairwise F ST values between the Pacific group and the western Indian Ocean group was 0.180, while 0.094 between the Pacific group and the eastern Indian Ocean group and 0.053 between the west Indian Ocean group and the east Indian Ocean group (Table 3).

Phylogenetic Partitions
After filtering, 7379 genes were used to reconstruct the phylogenetic tree. In the concatenate tree reconstructed using RAXML, samples of individuals of P. chinensis collected from the western Pacific were clustered into one lineage while those from the Indian Ocean were clustered into another. The pattern of divergence within P. cinereus was similar as in P. chinensis, except that P. cinereus collected from the Indian Ocean were clustered into two sublineages, the eastern Indian Ocean and the western Indian Ocean (Figure 4). The coalescent-based species tree reconstructed using ASTRAL-III exhibited the same partitioning pattern for P. chinensis and P. cinereus ( Figure S3). . Concatenated phylogenetic tree of Pampus chinensis and P. cinereus. Each species is mainly diverged into two lineages that are respectively the Pacific and Indian Ocean lineages. Note that two outgroups Peprilus medius and Psenopsis anomala were omitted from this figure, but P. punctatissimus was retained as the outgroup of P. chinensis. Bootstrap values were beside the node.

Divergence Time and Demographic History
A total of 1437 loci were concatenated and used for estimation of divergence time. In the result, the estimated first occurrence of genus Pampus was in the early Miocene,  Figure 5). For the five species of genus Pampus, the estimated first occurrence of P. chinensis was~2.05 Mya, P. cinereus was~5.53 Mya, P. punctatissimus was~3.76 Mya, P. argenteus was~4.63 Mya, and P. minor was~3.3 Mya. Particularly, for P. chinensis, the estimated first occurrence of the Indian Ocean and Pacific lineages were~1.6 Mya and~1.43 Mya, respectively. For P. cinereus, the Indian Ocean lineage was first occurred in~4.55 Mya, and the Pacific was first occurred in~2.14 Mya ( Figure 5).  (Table S2). In the FASTSIMCOAL2 analysis, a total of 5557 independent loci were used to generate folded SFS for P. chinensis. The best-fitting demographic model for P. chinensis was model 5, with the highest likelihood ( Figure S4a) and lowest AIC (the second-best model ∆AIC = 84.62), indicating that the Pacific and the eastern Indian Ocean groups were divergent from a common ancestor around 898,108 generations ago (95% CI 888,139-917,095). Since the gene flow was estimated backward-in-time, the forward gene flow should be as the same value but reversed in direction. In this context, estimated gene flow from the Pacific to the eastern Indian Ocean (4.58 × 10 −7 , 95% CI 4.43 × 10 −7 -4.66 × 10 −7 ) was stronger than that from the eastern Indian Ocean to the Pacific (1.19 × 10 −7 , 95% CI 1.12 × 10 −7 -1.23 × 10 −7 ; Figure 7a). For P. cinereus, a total of 7222 independent loci were used to generate folded SFS. The likelihood distribution of model 6 and model 8 overlapped ( Figure S4b); however, model 6 had a lower AIC than model 8 (∆AIC = 10.32). Hence, model 6 was chosen, indicating that the Pacific group and the eastern Indian Ocean group were divergent from a common ancestor 1,543,404 generations ago (95% CI 1,417,135-1,659,188), and after a few generations (1,141,612 generations ago, 95% CI 1,141,526-1,321,759). The western Indian Ocean group diverged from the eastern Indian group, and gene flow persisted within and between oceans. The estimated forward-in-time gene flow from the Pacific to the eastern Indian Ocean (9.22 × 10 −6 , 95% CI 0-9.79 × 10 −6 ) was stronger than that from the eastern Indian Ocean to the Pacific (4.46 × 10 −7 , 95% CI 0-1.29 × 10 −5 ). Moreover, the gene flow from the eastern Indian Ocean to the western Indian Ocean (1.05 × 10 −6 , 95% CI 0-3.11 × 10 −5 ) was stronger than the reverse direction (2.07 × 10 −7 , 95% CI 9.19 × 10 −8 -2.75 × 10 −5 ; Figure 7b).

Ancestral Area
The best-supported model for ancestral area reconstruction was Dispersal-Extinction-Cladogenesis + jump dispersal (DEC + j), which means the ancestor of Pampus could undergo all kinds of anagenetic (dispersal, extinction) and cladogenetic (narrow-, fullsympatry or allopatry, etc.) events [63,64]. In the analysis, 15 dispersal and vicariance events were inferred within the genus Pampus, and the results suggested that the South China Sea was the ancestral area for genus Pampus (Figure 8).

Intraspecific Diversification of P. chinensis and P. cinereus
Interoceanic differentiation was detected within both P. chinensis and P. cinereus. In the results of STRUCTURE analysis, P. chinensis clustered into two groups, and the same result was found for P. cinereus. The DAPC analysis showed similar results, except that P. cinereus was further split into three clusters. The clustering of these two species was similar, according to oceanic areas, suggesting that the IPB acted as a barrier between the Indian Ocean and the western Pacific for a long enough time to derive genetic differentiation. The same effect was observed in other marine taxa that are distributed among the Pacific and the Indian Ocean, such as Labroides dimidiatus, Plectroglyphidodon lacrymatus, Sargocentron diadema, etc. [65]. Some differences also were seen in the F ST statistics results. When we grouped the populations according to their oceanic areas, P. chinensis showed a large differentiation between the eastern Indian Ocean and the western Pacific, but P. cinereus showed a medium level of differentiation between the eastern Indian Ocean and the western Pacific, and large differentiation between the western Indian Ocean and the western Pacific, which indicated that P. cinereus may be further diverged due to the barrier of the Indian subcontinent.
The partitioning patterns also were observed in the phylogenetic tree, which explicitly divided into the Indian Ocean lineage and the western Pacific lineage (Figure 4) for each species. Interestingly, the sample of P. chinensis from the Strait of Malacca was diverged from other individuals of the Indian Ocean lineage (Figure 4), and STRUCTURE analysis indicated that its genetic variants were roughly half from the Indian Ocean and half from the western Pacific (Figure 2a). Since this sampling location was near the IPB, hybridization may have happened between the individuals from the eastern Indian Ocean and western Pacific to form a hybrid population; however, more individuals should be sampled from this location to provide robust evidence. P. cinereus of the Indian Ocean lineage was separated into the western Indian Ocean lineage and the eastern Indian Ocean lineage, which may be a consequence of the isolation between the Arabian Sea and the Bay of Bengal. Divya et al. [20] reported that Pampus species collected from the Arabian Sea and the Bay of Bengal have affinity to P. cinereus collected from the South China Sea. In their COI tree, these three populations were clustered into three lineages, which was the same as our topology. Combining the results based on mitochondrial DNA and nuclear DNA, respectively, we suggest that three lineages exist within P. cinereus. Divya et al. [21] described the Pampus species collected from the Arabian Sea as a new species, P. candidus, but we think further scrutiny should be made to split more species from P. cinereus lineages before distinct morphological characters are identified. On the other hand, P. liuorum, a controversial Pampus species, was defined as a synonym of P. cinereus in our previous study [15], but a recent study justified P. liuorum as a valid species using two mitochondrial DNA barcodes [66]. In their Cytb trees, P. cinereus and P. liuroum sampled in the western Pacific were separated in two lineages, which was different from our result, but their COI tree corroborated our result that P. cinereus and P. liuroum are not reciprocal in monophyletic. More comparison of specimens and integrated data sets should be collected to assess the status of P. liuroum. Nevertheless, our results suggest that the IPB is a major factor diversifying P. chinensis and P. cinereus, both of which consist of an Indian Ocean lineage and a Pacific lineage.

The Mechanism of Forming Different Intraspecific Lineages
Since the IPB has shown its barrier effect on P. chinensis and P. cinereus, exploring the mechanism of how P. chinensis and P. cinereus formed into different lineages is meaningful. Due to glaciation in the Pleistocene, the sea level was below its present level by as much as 120 m [67], and there have been at least eight 100 thousand-year cycles [68]. In these periods, the Sunda Shelf, which forms the IPB, was mostly exposed, restricting exchange between the western Pacific and the tropical Indian Ocean [2]. Estimated time-tree showed that the occurrence of each lineage of P. chinensis was in the mid-late Pleistocene, and lineages of P. cinereus were in the early Pleistocene or even late Pliocene. We hypothesize that divergence between oceans of both P. chinensis and P. cinereus were due to sea level fluctuations caused by repeated glaciation in the Pleistocene that isolated the populations in these two ocean areas for a long time. Combining the results, we propose these scenarios for species divergence: (1) for P. chinensis,~2 million years ago from the South China Sea, the common ancestor of the Indian Ocean and Pacific lineages started to disperse and then diverged. Temperature dropped sharply after the new populations colonized the eastern Indian Ocean and western Pacific, dropping the sea level below the then-present level and the IPB isolated the two new populations until the end of the ice age. Due to a long time of geographic isolation, the Indian and Pacific lineages were formed and subsequently restored a low-level of gene flow between the two oceans until present. (2) The divergence of P. cinereus was similar to the scenario of P. chinensis, as though it occurred in an earlier age. The common ancestor of the western and eastern Indian Ocean lineages and Pacific lineage dispersed from the South China Sea, and they first colonized in the eastern Indian Ocean and then the western Indian Ocean. Similarly, the mostly exposed Sunda Shelf interrupted gene flow between the Indian Ocean and Pacific and formed the Indian lineage and Pacific lineage. The genetic breaks caused by the IPB are also reported in other marine organisms [69]. However, in the Indian Ocean lineage, there are sublineages in the western Indian Ocean and the eastern Indian Ocean. We suggest this should not be affected by the glaciations but maybe the founder effect. Meanwhile, only both the eastern Indian Ocean lineages of P. chinensis and P. cinereus obtained significant negative Tajima's D values; combining the Bayesian skyline plots result, we suggest all lineages of P. chinensis and P. cinereus are experiencing a demographic plateau period since diverging expansion.
FASTSIMCOAL2 results estimated that divergences between the Indian Ocean and Pacific lineages of P. chinensis and P. cinereus were 898,108 generations ago and 1,543,404 generations ago, respectively. Previous studies reported that P. argenteus spawns twice a year [16], and all individuals in a population reach sexual maturity in two years [70]. Since the reproductive biological characteristics are similar within the genus [17], the generation times of P. chinensis and P. cinereus were set to 2 years. Therefore, the divergence between the Indian Ocean lineage and Pacific lineage of P. chinensis and P. cinereus should be around 1,796,216 years ago and 3,086,808 years ago, respectively. The estimated time frame overlaps with the result of our reconstructed time-tree. In other words, this result also suggests that the low sea level in the Pleistocene may have caused the diversification within both P. chinensis and P. cinereus between the western Pacific and the Indian Ocean.

The Origin and Dispersal Routes of Genus Pampus
The time-tree estimated that the common ancestor of all Pampus species lived~22.22 Mya, a period of explosive diversification in most modern marine fishes [13]. Our result suggests that dispersals and diversifications were carried on subsequently. The first species appearing was P. cinereus in the South China Sea, and then it began westward dispersal to the eastern Indian Ocean and then the western Indian Ocean and northward dispersal to the East China Sea. The second species diverging was P. argenteus in the South China Sea, and then it began dispersal to the Bohai Sea and then colonized the East China Sea. P. punctatissimus and P. minor appeared almost at the same time, while P. punctatissimus dispersed northward to the East China Sea and P. minor northward to the Taiwan Strait. The most modern Pampus species is P. chinensis that first appeared in the South China Sea and then began to disperse westward to the eastern Indian Ocean and northward to the Taiwan Strait ( Figure 9). The dispersal routes may have been facilitated by the ocean currents in the western Pacific and north Indian Ocean. We suggest that the South China Sea southern cyclonic gyre [71] would facilitate the individuals through the shallow sea strait and reach the eastern Indian Ocean, and then the Monsoon Current in the north Indian Ocean [72] would facilitate further westward dispersal. The Guangdong Coastal Current, South China Sea Warm Current, South China Sea Branch of Kuroshio and Luzon cyclonic gyre [71] would facilitate individuals through the Taiwan Strait and Luzon Strait, and then the Kuroshio Current [73] would facilitate further northward dispersal. The Indo-Australian Archipelago, also known as the Coral Triangle (CT), is a marine biodiversity center which contains numerous marine species [5]. Three hypotheses were proposed to explain the high biodiversity in the CT: center of origin, center of overlap and center of accumulation. Although Bellwood and Meyer [74] suggested the center of accumulation model could best describe the high biodiversity in CT, Briggs [75] proposed that the species accumulation that happened in the Eocene of the West Tethys and CT had its ability to produce species after the early Miocene. In our results, the occurrence and most probable origin of genus Pampus were estimated to be the early Miocene (~22.22 Mya) in the South China Sea, which is near the CT. Thus, our results provide new evidence for the center of origin model for CT and emphasize the CT as an important region for biodiversity conservation.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d14030180/s1, Figure S1: Mean log-likelihood Pr (X|K) and delta K, Figure S2: BIC values for each cluster, Figure S3: Coalescent-based species tree, Figure S4: Likelihood distributions for different models, Table S1: Numbers of enriched genes for each individual, Table S2: Values of Tajima's D test and Fu's Fs test.
Author Contributions: C.L. conceived this study and collected the samples. G.Y. and A.S. collected the samples. G.F. conducted the experiments, analyzed the data and wrote the manuscript. All authors discussed the results and edited the final manuscript. All authors have read and agreed to the published version of the manuscript.

Informed Consent Statement: Not applicable.
Data Availability Statement: The raw FASTQ reads are available in the NCBI repository with accession number PRJNA809874. All scripts used in this study were deposited in a GitHub repository (link to this repository: https://github.com/naFgG/Scripts_used_in_Pampus, 23 February 2022).