First Evidence of Cryptic Species Diversity and Population Structuring of Selaroides leptolepis in the Tropical Western Pacific

The yellowstripe scad, Selaroides leptolepis (Carangidae), is an important fish commodity in the Tropical Western Pacific (TWP). It has a latitudinal Pacific range from south of Japan down to northern Australia, with the highest concentration in Southeast Asia. However, its TWP fishing grounds have long been a hotspot of unsustainable exploitations, thus threatening the remaining wild populations. Despite the species’ commercial significance, there is limited understanding of its genetic structure and diversity. Herein, the genetic structure of S. leptolepis was examined using mitochondrial COI and CytB sequences. Both markers denoted significant genetic structuring based on high overall FST values. Hierarchical analysis of molecular variance (AMOVA), maximum likelihood (ML) phylogenetic trees, and median-joining (MJ) haplotype networks strongly supported the occurrence of two allopatrically distributed lineages. These comprised of a widespread Asian lineage and an isolated Australian lineage. Within-lineage distances were low (K2P < 1%) whereas across-lineage distances were remarkably high (K2P > 6%), already comparable to that of interspecific carangid divergences. Haplotype sequence memberships, high genetic variations, and the geographic correlation suggested that the Australian lineage was a putative cryptic species. Historical demographic inferences also revealed that the species experienced rapid expansion commencing on the late Pleistocene, most likely during the end of the Last Glacial Maximum (∼20,000 years ago). The present study encouraged the application of lineage-specific management efforts, as the lineages are experiencing different evolutionary pressures. Overall, accurate knowledge of the species’ genetic distribution is fundamental in protecting its diversity and assuring stock sustainability.

The yellowstripe scad, Selaroides leptolepis (Carangidae), is an important fish commodity in the Tropical Western Pacific (TWP). It has a latitudinal Pacific range from south of Japan down to northern Australia, with the highest concentration in Southeast Asia. However, its TWP fishing grounds have long been a hotspot of unsustainable exploitations, thus threatening the remaining wild populations. Despite the species' commercial significance, there is limited understanding of its genetic structure and diversity. Herein, the genetic structure of S. leptolepis was examined using mitochondrial COI and CytB sequences. Both markers denoted significant genetic structuring based on high overall F ST values. Hierarchical analysis of molecular variance (AMOVA), maximum likelihood (ML) phylogenetic trees, and median-joining (MJ) haplotype networks strongly supported the occurrence of two allopatrically distributed lineages. These comprised of a widespread Asian lineage and an isolated Australian lineage. Within-lineage distances were low (K2P < 1%) whereas acrosslineage distances were remarkably high (K2P > 6%), already comparable to that of interspecific carangid divergences. Haplotype sequence memberships, high genetic variations, and the geographic correlation suggested that the Australian lineage was a putative cryptic species. Historical demographic inferences also revealed that the species experienced rapid expansion commencing on the late Pleistocene, most likely during the end of the Last Glacial Maximum (∼20,000 years ago). The present study encouraged the application of lineage-specific management efforts, as the lineages are experiencing different evolutionary pressures. Overall, accurate knowledge of the species' genetic distribution is fundamental in protecting its diversity and assuring stock sustainability.

INTRODUCTION
The yellowstripe scad, Selaroides leptolepis, is an economically important fish commodity in the Tropical Western Pacific (TWP). This species is highly exploited in Southeast Asian countries, particularly Malaysia, Indonesia, the Philippines, and in the United Arab Emirates (Kempter, 2015; Figure 1). The S. leptolepis is commonly consumed fried, steamed, sun-dried, and even prepared as components of surimi and burger patties (Yu and Siah, 1998;Arfat and Benjakul, 2012). Aside from being a good source of dietary protein, it also contains other functional biomaterials such as protein hydrolyzates and histamine (Klompong et al., 2009;Huang et al., 2010). Since the market price of S. leptolepis is relatively more affordable than other fish groups (Kempter, 2015), it is highly patronized by local consumers. The species' commercial significance, however, also makes it a vulnerable target of extensive exploitations.
The S. leptolepis is the only representative of the genus Selaroides (Family Carangidae). Its closest relatives include Selar, Alepes, Hemicaranx, Chloroscombrus, and Caranx of the Tribe Carangini (Reed et al., 2002). Taxonomically, S. leptolepis is often misidentified as its Carangid relative-Selar crumenophthalmus (Bloch et al., 1793) due to morphological similarities. S. leptolepis differs by the absence of papilla on the lower pectoral girdle, absence of teeth on the upper jaw, and prominence of the longitudinal yellow stripe on its body (Nakabo, 2002). As a demersal species, S. leptolepis is commonly found in inshore waters shallower than 50 m (Allen and Erdmann, 2012). It is widely distributed throughout the tropical and subtropical Indian and Pacific waters (Abdussamad et al., 2013). In the TWP, it has a latitudinal distribution range from the south of Japan (Hata et al., 2017), down to the Indo-Malay-Philippine archipelagos (IMPA) (Mat Jaafar et al., 2012), and northern Australia (Dell et al., 2009).
Considered as one of the globally most important fishing regions, the TWP accounts for almost half of the world's marine fisheries production (Food and Agriculture Organization [FAO], 2020) and houses the top capture fish producers of the world, namely China, Indonesia, Vietnam, and Japan (Food and Agriculture Organization [FAO], 2020). Its IMPA region, in particular, is also considered a megadiverse area, with several studies documenting the presence of genetically distinct fish populations in the region (Rohfritsch and Borsa, 2005;Salini et al., 2006;Hubert et al., 2012;Thomas et al., 2014). The high biodiversity index of TWP had been attributed to its geologic history, oceanographic dynamics, capacity for biomass support, and vicariance histories (Carpenter and Springer, 2005;Gaither and Rocha, 2013). Unfortunately, some TWP fishing grounds have been experiencing high percentages of unsustainable fisheries. Stock depletions have already been reported in several fishing grounds within the TWP (Teh et al., 2007;Guanco et al., 2009;Tangke et al., 2018;Fauziyah et al., 2020). Records also revealed a steady decline in the global yield of S. leptolepis after reaching peak production in 2014 (Food and Agriculture Organization [FAO], 2021). The conservation of S. leptolepis and other demersal fishes is important since these groups constitute a significant portion of the TWP capture fisheries (Food and Agriculture Organization [FAO], 2020).
In fisheries management, the reduction of genetic diversity on the remaining natural populations remains a major problem. This loss of genetic diversity translates into reduced population fitness and its inability to adapt against evolutionary pressures (Kenchington, 2003). To address these problems, the incorporation of molecular techniques in conservation studies had been practiced. These techniques can reveal valuable genetic information, including a species' pattern of distribution and the demographic history of its populations. Clear information on the genetic structures allows us to define management zones and assess risks to overexploitation. Only a few records have documented the genetic distribution of S leptolepis, and all were done on small regional scales (Kempter, 2015;Mat Jaafar et al., 2020). There is no available information on its extant genetic structure and diversity in the TWP scale, hence this study.
The objective of this study was to investigate the intraspecies genetic diversity, structure, and demographic history of S. leptolepis in TWP. Sequence data from the mitochondrial DNA Cytochrome Oxidase I (COI) and Cytochrome B (CytB) regions were utilized to infer genetic structures. The mtDNA genome is conserved across animal lineages, contains few duplications, mutates rapidly, and is relatively easy to isolatemaking it a good marker for intraspecies genetic studies (Desalle et al., 2017;Zhang et al., 2020). The COI and CytB have moderate evolutionary rates than other mtDNA genes (Kochzius et al., 2010;Kartavtsev, 2011) and their fragments contain adequate informative phylogenetic information (Liu et al., 2021). The coupled use of COI and CytB had not only been applied to barcode animal taxa and infer deep intraspecies divergences (Baker et al., 1995;Smith et al., 2011;Çiftçi et al., 2013;Joshi et al., 2019), but also to detect provisional cryptic populations awaiting thorough screenings (Asgharian et al., 2011;Hubert et al., 2012;Sienes et al., 2014;Thongtam Na Ayudhaya et al., 2017). Findings from this study serve as baseline information for the fishery management and sustainable utilization of S. leptolepis given its economic relevance and understudied status.

Sampling and DNA Extraction
Selaroides leptolepis individuals were obtained from landing sites and wet markets from Central Philippines (n = 132) and Southern Taiwan (n = 43) between 2018 and 2021. Identification was based on morphological characters (i.e., longitudinal yellow stripe from upper part of eye to caudal peduncle, toothless upper jaw, and lower pectoral girdle absent of papilla; Nakabo, 2002;Motomura et al., 2017;Koeda and Ho, 2019). Approximately 1 g of muscle from the right caudal peduncle of each fish were stored in 2 ml tubes with 95% ethanol prior to genomic extraction. Donated muscle tissues from museum collections were also included and subjected to genomic extractions. These samples originated from Taiwan (n = 6), Philippines (n = 1), Malaysia (n = 1), and Australia (n = 16). More information regarding these samples can be found in Supplementary Table 1. Extraction was carried out using ReliaPrep TM gDNA Tissue Miniprep System (Promega) or Genomic DNA Extraction Kit 2.0 (Yeastern Biotech) following manufacturers' protocols.

Polymerase Chain Reaction Amplification and Sequencing
The COI gene portion was amplified using published primers, namely, FishF1 (5 -TCA ACC AAC CAC AAA GAC ATT GGC AC-3 ) and FishR1 (5 -TAG ACT TCT GGG TGG CCA AAG AAT CA-3 ) (Ward et al., 2005); whereas the CytB regions were amplified using CytbF (5 -GGC TGA TTC GGA ATA TGC AYG CNA AYG G-3 ) and CytbR (5 -GGG AAT GGA TCG TAG AAT TGC RTA NGC RAA-3 ) (Kochzius et al., 2010). Amplifications were performed in 20 µl reaction volumes containing 10 µl Ampliqon Taq DNA Polymerase Master Mix, 1 µl of genomic DNA, 0.5 µl of each primer, and 8 µl of ultrapure water. Polymerase chain reaction (PCR) for both markers was carried out with an initial denaturation at 94 • C for 2 min, followed by 35 cycles of 94 • C denaturation for 30 s, 52 • C annealing for 40 s, 72 • C extension for 1 min, and final extension of 72 • C for 10 min. The quality of the PCR products was evaluated in 1.5% agarose gel. PCR amplicons were either sent to Macrogen Inc. (Seoul, South Korea) or Genomics (New Taipei, Taiwan) for bidirectional sequencing using the mentioned respective PCR primers. Obtained sequences were quality checked and assembled in Sequencher 5.4.6 (Ann Arbor, MI, United States).

Data Analyses
Additional sequences publicly stored in GenBank were also incorporated in the analyses (Supplementary Table 2
For phylogenetic trees and haplotype networks generations that were not sensitive to sample sizes, individuals from sparsely represented geographic regions (n > 10) were incorporated to explore the global relationship among all the available sequences. These regions included SJPN, SCHS, JIDN, SIND, and ABGF for COI, and SJPN, SCHN, CVNM, PMYS, and ABGF for CytB (Table 1 and Supplementary Table 2). The inclusion of these sequences was relevant, as these regions represented individuals from the species' holotype collection site (i.e., JIDN, Cuvier and Valenciennes, 1833) and marginal distribution range (i.e., Southern Japan, Indian Ocean; Prabhu, 1956;Abdussamad et al., 2013;Hata et al., 2017). Maximum likelihood (ML) phylogenetic trees were created in MEGA X under substitution models K2P (COI) and K2P + G (CytB), with S. crumenophthalmus as the outgroup. PopART 1.7 (Leigh and Bryant, 2015) was used to generate the median-joining (MJ) haplotype networks.

Demographic History
Neutrality tests and effective population size change estimations were done for the n >10 sample groups. These were performed to infer historical demography and evolution neutrality. Deviation from the neutrality model was calculated using Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) in Arlequin 3.5.2.2. These indices indicate whether populations underwent expansions. Further inferences on historical demography were carried out with a mismatch distribution analysis. Demographic parameters such as tau (τ), θ 0 , θ 1 , sum of squared deviation (SSD), and Harpending's raggedness index (Hri) were also calculated in Arlequin 3.5.2.2. Graphical figures showing the pairwise comparison between the frequency of individuals (y-axis) with the corresponding number of pairwise differences (x-axis) were generated using DnaSP v6.12. Changes in effective population size (N e ) across time were inferred using Bayesian skyline plot analysis (Drummond et al., 2005) implemented in BEAST 2.6.5 (Bouckaert et al., 2019). XML files were initially prepared with BEAUti 2.6.5 (Bouckaert et al., 2019). The HKY + G nucleotide substitution model was selected for both markers to take into account possible sitespecific variations (Hill and Baele, 2019), and applied a strict clock mutation rate of 1 × 10 −8 per site per year as suggested for reef fishes (Stewart Grant et al., 2012;Delrieu-Trottin et al., 2017). Independent Markov chain Monte Carlo (MCMC) analyses were ran for 100 million generations with a burn-in of 10 million and sampled every 1,000 iterations. If necessary, runs were repeated until combined ESS >200 values were attained, and consensus of these parameter values were visualized in Tracer 1.7.2 (Rambaut et al., 2018).

Genetic Diversity
A total of 306 and 205 sequences were generated for COI and CytB, respectively. Fifty-one haplotypes were identified from the 517-bp COI fragment (GenBank accession numbers: MZ520638-MZ520664), whereas 49 haplotypes were identified from the 527-bp CytB (GenBank accession numbers: MZ55565 8-MZ555703). Thirty-eight singletons were found for COI and 40 for CytB. Overall, 69 and 65 polymorphic sites were identified for COI and CytB, of which 45 and 40 were parsimony informative, and 22 and 25 were singleton variables, respectively. Global genetic diversities were high (COI global Hd = 0.7825; CytB global Hd = 0.7364) while nucleotide diversities were low (COI global π = 0.0119; CytB global π = 0.0106). Corresponding nucleotide composition for COI and CytB were 23.21 and 22.69% adenine, 29.15 and 28.63% thymine, 29.68 and 33.22% cytosine, and 17.96 and 15.46% guanine. The genetic indices for each group were presented in Table 2.

Genetic Structure and Phylogeographic Relationships
Overall, pairwise F ST values for both markers showed significant genetic variation (p < 0.05) across populations. For COI, AUST, and CPHL significantly differed from the rest of the groups.  Table 3. Hierarchical AMOVA also showed the highest among-clusters F CT variability when AUST was separated from the others (COI = 0.8847; CytB = 0.9491), though the variation was statistically insignificant at p > 0.05. Additional AMOVA results from other hypothetical combinations are shown in Table 4. Combining all the other available sequences, the presence of two diverging lineages was detected from the network topologies (Figure 2) and phylogenetic trees (bootstrap > 75%) (Figure 3). Herein, these lineages were referred to as the Asian and Australian lineages. The Asian lineage covers all individuals excluding AUST, while the Australian lineage is comprised only of AUST individuals. The MJ haplotype networks for each marker showed a deep divergence on the Asian and Australian lineage in a reciprocally monophyletic network. Such type of network is characterized by the presence of more than one lineage, usually separated by numerous mutational steps (Jenkins et al., 2018). Pairwise K2P distances across the two lineages were COI = 6.77% and CytB = 6.64%. Within-lineage K2P differences were COI = 0.60% and CytB = 0.22% for the Asian lineage, and COI = 0.14% and CytB = 0.23% for the Australian lineage. Corresponding K2P differences across each collection region were presented in Table 3.
Further analyses without AUST were carried out to detect variations within the Asian lineage. F ST values revealed significant differences across Asian populations (COI = 0.6511; CytB = 0.4406). Only the COI ML tree showed a further divergence of the Asian lineage into two main sublineages with 76% bootstrap support. The first sublineage consisted of STWN and most of CPHL individuals, as well as a few PMYS, EMYS, SIND, SJPN, and ABGF individuals. On the other hand, the second sublineage included STWN, SECH, CVNM, PMYS, EMYS, JIDN, and CPHL individuals. These two Asian sublineages had a pairwise difference of K2P = 0.99%. The COI network displayed the divergence by the separation of its two dominant haplotypes-H01 and H15. H01 includes individuals from STWN, SECH, CPHL, PMYS, EMYS, and the additional representatives from JIDN, SIND, SCHS, and CVNM. Meanwhile, H15 is composed exclusively of CPHL individuals. For CytB, the dominant haplotypes were H01 and H02. H01 is majorly represented by individuals from CPHL and with few samples from STWN and SJPN, while H02 is comprised mostly of STWN and few CPHL, PMYS, and CVNM individuals (  p < 0.01) ( Table 2). Total mismatch distribution for COI and CytB reflected different demographic signatures. COI showed a multimodal pattern whereas CytB was bimodal (Figure 4). SSD value for COI revealed an insignificant difference (p > 0.05) from a predicted growth expansion model. Raggedness values (Hri) for both markers also showed statistical insignificance (p > 0.05), implying the samples had a relatively good fit to a population expansion model. Corresponding mismatch indices for the lineages were shown in Figure 4. Using COI mismatch parameters θ 0 and θ 1 (Harpending, 1994;Marini et al., 2021), the estimated effective female population size for S. leptolepis after expansion (θ 1 ) was approximately 5,000 times higher than prior (θ 0 ).
Demographic scenarios supporting the recent population expansion of S. leptolepis were presented in the Bayesian skyline plots (Figure 5). Both COI and CytB revealed patterns of a long history of constant population size, followed by a slight decline (bottleneck), and a subsequent demographic expansion. The fastest increase of its effective population N e happened between 21 and 10 thousand years ago (KYA). It showed a relatively stable effective population starting 2,000 years ago until the present. At the lineage level, the Asian lineage also displayed the rapid expansion signature, whereas the Australian lineage depicted slow population growths (Supplementary Figure 1).

Genetic Diversity
High overall haplotype and low nucleotide diversities were recorded for both markers, with numerous unique haplotypes or singletons present in the MJ haplotype networks. These singletons directly radiate from a largely shared haplotype, indicating few mutational step differences. On one hand, high haplotype diversity was detected when numerous unique sequences were present within the overall population. On the other hand, low nucleotide diversity was reflected when these nucleotide compositions were closely similar. This type of genetic pattern is usually attributed to a recently experienced expansion (Grant and Bowen, 1998), which is likewise supported by our historical demography findings. In such cases, individuals evolve into different haplotypes with minimal differences; and these haplotypes may either evolve directly or indirectly from an ancestral haplotype (Chanthran et al., 2020). This kind of genetic pattern has also been recorded in other TWP carangids (Rohfritsch and Borsa, 2005;Jamaludin et al., 2020;Mat Jaafar et al., 2020;Torres and Santos, 2020). The most recurrent and widespread haplotype is considered the oldest and most successful in traversing across sampling locations (Posada and Crandall, 2001;Mat Jaafar et al., 2012). This indicates that the second Asian sublineage might be the ancestral S. leptolepis lineage and, the IMPA waters could be its potential geographic origin.
The separation of the two lineages can probably be related to habitat discontinuity, oceanographic barriers, the species' dispersal capacity, and past geologic events. For instance, additional collections from Eastern Java, Indonesia corresponded with the widely distributed Asian haplotype-H01. This sampling locality is known to exhibit continuous coral and rocky reef bottoms that extend to the Nusa Tenggara region (Fahmi et al., 2021). The Eastern Indian Ocean and the Timor Sea separate Eastern Java, Indonesia from North Australia, which is the southernmost range of the species in TWP. These two regions have a proximate distance of ∼1,400 km and reach a maximum depth of ∼3,300 m through the Timor Trough, implying that water dynamics between the two neighboring localities might have acted as barriers which prohibited their genetic exchange. The small islands in the area are also known to serve as a gateway for the strong intrusive Pacific waters flowing westwards to the Indian Ocean (Gordon, 2005). Being an inshore, demersal fish, these suggest that S. leptolepis has a weak migration capability to overcome deep water and strong current conditions despite proximity. In contrast, genetic homogeneity was displayed within the Australian group (K2P ≤ 0.2%), which consisted of individuals collected from Queensland and Western Australia (Supplementary Table 1). Found on opposite sides of the continent, this reflects an extant gene flow and signifies a high dispersal capability for S. leptolepis in shallow coastal regions. All in all, oceanographic hindrances and physical limitations explained the separation of Asian and Australian lineages. Other examples of fishes showing homogenous North Australian groupings with sharp discontinuities along the Timor Sea include Lutjanus erythropterus (Salini et al., 2006), Decapterus russelli (Rohfritsch and Borsa, 2005), and Pristipomoides multidens (Ovenden et al., 2002). Additional samples from the Lesser Sunda and nearby island localities will be ideal to establish a more comprehensive narrative of the S. leptolepis distribution.

Shallow Structuring Across Asia
A broadscale geographic homogeneity was observed in the Asian samples. Within this lineage, K2P divergences were <1% on both markers. However, despite the evident genetic differences in our findings, the magnitude of these differences was low and inadequately supported (i.e., insignificant F ST s, low bootstrap support on CytB, and mixed haplotype and sublineage memberships). The discrepancy between the phylogenetic results of COI and CytB, where further substructuring was apparent in COI, might be an artifact of the level of sensitivity of the mtDNA marker. This undetected subclustering implies that such divergence is shallow. Moreover, the limited inference on phylogenetic trees is complemented with the use of networks. Since genetic diversity is usually low at the population level, this leads to indecisive tree resolutions and further overlooks other important evolutionary information. To effectively visualize reticulated relationships such as hybridization and recombination, networks are used. The implicit, sequence-based MJ haplotype network may suggest the possible occurrence of extant unsampled sequences or extinct sequences through median vectors (Bandelt et al., 1999;Kong et al., 2016 ; Figure 2). The mixed memberships of the dominant Asian haplotypes imply sympatric distribution and the absence of a clear genetic structure in the Asian lineage.
The absence of genetic structuring may partly be explained by the spawning patterns and pelagic larval duration of S. leptolepis, as well as the hydrogeographic features of the SCHS. S. leptolepis exhibits two annual spawning events succeeded by recruitment pulses which usually peak from March to June when the temperature is warmer (Prabhu, 1956;Guanco et al., 2009). In addition, observations from its close relative, S. crumenophthalmus, reveal an 18-day post-hatching time before reaching a nursing flexion stage (Welch et al., 2013), which supports the notion that the pelagic larval time of S. leptolepis might be close to 2-3 weeks. Our S. leptolepis Asian collection grounds also lie on the periphery of the SCHS, which is bounded to its adjacent areas by shallow waters and straits. The water movement of the SCHS's top layer is mainly influenced by seasonal monsoons, with its water circulation directly affected by two oppositely headed winds that sweep over the area at different times of the year (Huang et al., 1994). Consequently, the spawning times S. leptolepis temporally coincide with the timings of the opposing monsoons. With all these being said, the buoyancy of its ichthyoplankton, assisted with oceanographic circulators and ample pelagic larval durations, could be key factors in the successful genetic dispersal in the region.
Other notable clusterings of some haplotypes are shown in the ML trees and MJ networks. For instance, majority of the Malaysian haplogroup of Mat Jaafar et al. (2020) is now recognized as part of the widespread second Asian sublineage. This implies that most of the Malaysian S. leptolepis stocks are mainly comprised of direct ancient lineage descendants and thus support its IMPA origination. Additionally, SJPN sequences, which are collected from the northern-most TWP range of the species, appear to be more closely related to CPHL. This suggests a possible influence of the upward-moving Kuroshio current in dispersing the S. leptolepis to this northmost distribution. Also, the COI haplotype network showed that SIND (H49, H50) and ABGF (H51) samples formed a unique, close cluster. This suggests that S. leptolepis from the Indian Peninsula and ABGF region probably is a genetically distinct cluster. These samples are linked closely to the Peninsular Malaysian samples of Mat Jaafar et al. (2020), denoting that S. leptolepis from localities near the Malacca Strait might be genetically closer to its Indian Ocean conspecifics than those from the majority of IMPA. Additional samples are necessary to quantify their degree of genetic variability. This type of delineation was also observed in Seriola nigrofasciata (Mat Jaafar et al., 2012), Lutjanus lutjanus (Bakar et al., 2018), Uranoscopus cognatus (Mohd Yusoff et al., 2021), and Penaeus semisulcatus (Halim et al., 2021).

Demographic History
The negative Tajima's D and Fu's Fs indicated the presence of excessively rare haplotypes and imply a recent population expansion, thus rejecting the neutral evolution hypothesis. Fu's Fs has been considered as a more superior test, giving more reliability in inferring population growth (Ramos-Onsins and Rozas, 2002). A multi/bi-modal mismatch distribution suggested that the S. leptolepis populations are in equilibrium. However, relying solely on this graphical inference does not automatically warrant its respective history. Multiple and bimodality in mismatch analysis also happen in the presence of genetically distinct lineage in the samples. Its first peak in the graph represented the close intra-clade pairwise differences, and the next successive peaks illustrated a more ancient inter-clade pairwise difference. Therefore, lineages were suggested to be prior segregated to avoid possible violations in coalescent theory assumptions (Jenkins et al., 2018). In contrast, the unimodal patterns in each lineage depicted recent demographic expansion on the populations. Statistically insignificant values for SSD and Hri support such population expansion.
The Bayesian skyline plots suggested recent population expansion for S. leptolepis that initiated around the late Pleistocene era. The species recorded the most rapid increase in its effective population size during the Late Pleistocene (20 KYA) up to the early Holocene (10 KYA). During the Pleistocene, the TWP region, historically known as Sundaland, experienced glaciation and deglaciation processes which caused sea-level and temperature fluctuations that configured ocean dynamics. A contraction of the S. leptolepis effective populations was detected to coincide during these fluctuating periods . The formation of ice on continents and poles lowered sea levels, thus reducing the available space and ocean food supply for marine populations, and probably leading to depopulation. The abrupt rise in effective population size synchronously commenced during the Last Glacial Maximum's deglaciation timetable (∼20 KYA). This deglaciation caused a rapid rise in sea level which opened expansion opportunities. Once habitat conditions became acceptable, ancient S. leptolepis populations could have moved to these newly filled coastal regions. With the help of its high dispersal abilities in shallow environments, rapid population expansion, and new habitat colonization could have been easily achieved by the species. Overall, these collective circumstances might have influenced the present-day distribution patterns of S. leptolepis in the TWP coastal margins. Other relative carangids that exhibited similar expansion timelines in the region include D. maruadsi (Niu et al., 2019;Jamaludin et al., 2020), D. macrosoma, and D. macarellus (Arnaud et al., 1999).

Implications for Management and Conservation
The occurrence of two geographically isolated S. leptolepis lineages in TWP suggested at least two genetically distinct stocks were present in its waters. Their high level of K2P divergences was already at comparable levels for interspecies differentiation in Carangidae. The type locality for the original description of S. leptolepis was Java, Indonesia (Cuvier and Valenciennes, 1833). Specimens from Eastern Java grouped with the widespread Asian haplogroup suggested that this lineage might be the originally described S. leptolepis, and the Australian lineage is another putative species. Species complexes were considered recently diverged; therefore, their morphological differentiation is believed to have develop later due to new environment adaptations (Fahmi et al., 2021). A comprehensive and adequate collection on the S. leptolepis full range would also help us understand its global structure and detect the presence of any possible intermediate populations.
Selaroides leptolepis is of economic importance particularly in regions of Southeast Asia. However, this is also often coupled with a high exploitation rate. These Asian genetic stocks experience high fishery pressures driven by human consumption demands, which lead to localized fishery depletions. On the other hand, its Australian counterpart experiences relatively lesser pressures, as this species is not a top targeted fish commodity in the Australian region (Gunn, 1990). This species is usually documented as part of bycatch by fishery reports in this region (Blaber, 1993;Dell et al., 2009). Reports have indicated that >85% of Australia's fish stocks are well-managed and are at sustainable levels (Mobsby and Curtotti, 2018;Piddocke et al., 2020). This suggests that the two newly uncovered S. leptolepis lineages are experiencing different levels of fishery-induced pressures. Moreover, persisting localized unsustainable exploitations on the Asian stocks can cause fragmented isolation of small S. leptolepis populations. A lowered genetic diversity in these small S. leptolepis populations also means reduced fitness and higher risks against genetic degradation and drift. Small, isolated populations are also highly vulnerable to inbreeding, which reduces offspring number and viability. The S. leptolepis populations' ability in adapting to their constantly changing environment will also be restricted if genetic diversity is persistently decreased. If the highly targeted Asian lineage ultimately depletes, chances of replenishing its gene pool through migration will be unlikely due to the discussed isolating mechanisms.
Delimiting the S. leptolepis cryptic species complex will not only provide an advantage for its taxonomic recognition but can also aid in the formulation of better conservation measures. The challenge to fully delimit this S. leptolepis complex might rely on the combinatory use of genetic and non-genetic approaches. Exclusively associated characters in the morphology, reproductive traits, or habitat preferences will be helpful features for the species' field-based diagnostics. Otherwise, the absence of strong characters will make field-based differentiation impossible and even further complicate its management. Separating them into different units will foster more suitable stock-specific management approaches, especially since these groups are experiencing different exploitation pressures. Subdividing them will define their stock geographic boundaries, which can lead to more precise estimates on its fishery indices (e.g., recruitment, growth, and mortality) (Ovenden et al., 2009).
This study defined the genetic structure and the presence of a cryptic species of S. leptolepis in the TWP, wherein this species is of economic importance. Their isolated distribution, demographic history, and absence of inbetween populations warrant the clear separation; and this signifies the need for other diagnostic characters, whether morphological, habitat or behavioral, to disintegrate the species complex. A full distributional range survey coupled with robust genetic approaches [e.g., single nucleotide polymorphisms (SNP)] will reveal the global structuring and evolutionary history of this taxon. Regarding management and conservation, we recommend a lineage-specific approach since stocks face different environmental and fishery pressures.
A transnational management scheme can be designed for the widely distributed lineage. Most importantly, the integration of insights from genetic studies and other scientific information can foster the best management plan for this species in the future.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: MZ520638-MZ520664 and MZ555658-MZ555703).

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because no live samples were used in the study. All specimens were acquired from local fish markets.

AUTHOR CONTRIBUTIONS
LCH performed the experiments and data analyses. H-CL helped to supervise the study. LCH and H-CL drafted the manuscript. All authors conceived and designed the study, reviewed the drafts, and approved the final version of the manuscript.

FUNDING
This work was supported by the Department of Science and Technology-Science Education Institute, Philippines and Ministry of Science and Technology, Taiwan (Grant No. 109-2611-M-110-003).