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

Halasan, Lorenzo C.; Geraldino, Paul John L.; Lin, Hsiu-Chin

doi:10.3389/fmars.2021.756163

ORIGINAL RESEARCH article

Front. Mar. Sci., 18 October 2021

Sec. Marine Biology

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.756163

This article is part of the Research Topic The Philippine Seas: Biodiversity and Ecological Impacts of Natural and Anthropogenic stressors in Tropical Reef Systems View all 14 articles

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

$\r\nLorenzo C. Halasan$ Lorenzo C. Halasan¹

Paul John L. Geraldino²

Hsiu-Chin Lin^1,3*

¹Department of Marine Biotechnology and Resources, National Sun Yat-sen University, Kaohsiung, Taiwan
²Department of Biology, University of San Carlos, Cebu, Philippines
³Doctoral Degree Program in Marine Biotechnology, National Sun Yat-sen University, Kaohsiung, Taiwan

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

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.

FIGURE 1

Figure 1. (A) Major collection sites from this study indicated by www.frontiersin.org (n > 10 individuals) and (n < 10 individuals). The map is courtesy of The World Bank Org (2021). (B) Fresh and (C) dried S. leptolepis sold at a local fish market in Central Philippines (CPHL).

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 isolate—making 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.

Materials and Methods

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). COI sequences included the Malaysian dataset published by Mat Jaafar et al. (2020) (n = 70).

Genetic Diversity and Structuring

Sequences were aligned in BioEdit 7.2.5 (Hall, 1999) using Clustal W option with default settings and trimmed to 517- and 527-bp lengths, respectively, for COI and CytB. These were then exported to DnaSP 6.12 (Rozas et al., 2017) for clustering and haplotype generation. Sequences were grouped according to collection regions, namely: South Japan (SJPN), South Taiwan (STWN), Southeast China (SECH), Central Vietnam (CVNM), South China Sea (SCHS), Central Philippines (CPHL), Peninsular Malaysia (PMYS), East Malaysia (EMYS), Java Indonesia (JIDN), Southern Indian Peninsula (SIND), Arabian Gulf (ABGF), and Australia (AUST) (Figure 1 and Table 1). Selection of the best-fit nucleotide substitution models through the Akaike information criteria (AIC) were carried out separately for COI and CytB using MEGA X (Kumar et al., 2018). Since some analyses were sensitive to the sample sizes, only groups with n >10 sequences were included in the following calculations: estimations for genetic diversity (Hd), nucleotide diversity (π), and analysis of molecular variance (AMOVA)—which were carried out in Arlequin 3.5.2.2 (Excoffier and Lischer, 2010). K2P sequence divergence was calculated using MEGA X. For COI, the n >10 groups included STWN, SECH, CPHL, PMYS, EMYS, and AUST, while CytB only included STWN, CPHL, and AUST (Table 1).

TABLE 1

Table 1. Sampling region and abbreviations, sample size (N), and the number of haplotypes.

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 (τ), θ₀, θ₁, 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 site-specific 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).

Results

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.

TABLE 2

Table 2. Polymorphic sites (PS), haplotype diversity (Hd), nucleotide diversity (π), and Tajima’s D and Fu’s Fs indices for cytochrome oxidase I (COI) and cytochrome B (CytB).

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. AUST exhibited the highest pairwise differentiation (F_ST > 0.95) followed by CPHL (F_ST > 0.62). For CytB, all three sampling regions (STWN, CPHL, and AUST) significantly differed from each other. Pairwise F_ST values across each sampling group were presented in 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.

TABLE 3

Table 3. Population pairwise distances using F_STs (upper diagonals) and Kimura two-parameter (%K2P) (lower diagonals) for cytochrome oxidase I (COI) (upper table) and cytochrome B (CytB) (lower table).

TABLE 4

Table 4. Analysis of molecular variance (AMOVA) results for hierarchical genetic subdivision for the percentage of variation and F-statistics of the cytochrome oxidase I (COI) and cytochrome B (CytB) genes.

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.

FIGURE 2

Figure 2. Median-joining (MJ) networks of (A) cytochrome oxidase I (COI) and (B) cytochrome B (CytB) haplotypes of S. leptolepis Tropical Western Pacific (TWP) populations. The sizes of nodes (circles) are proportional to the number of individuals. Black dots represent median vectors (unsampled haplotypes).

FIGURE 3

Figure 3. (A) Cytochrome oxidase I (COI) and (B) cytochrome B (CytB) maximum likelihood (ML) trees of S. leptolepis rooted with S. crumenophthalmus (KJ502072, AY050731) from GenBank. Bootstrap consensus tree for the topology was inferred from 1,000 replicates and >75% support labels are retained at tree nodes.

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 (Figure 2 and Supplementary Table 3).

Demographic History

Overall Tajima’s D were negative for COI (D_global = −1.2809; p < 0.1) and CytB (D_global = −1.4956; p < 0.05). Fu’s Fs revealed highly significant negative values for both markers (COI_global Fs = −17.1323, p < 0.01) (CytB_global Fs = −22.5609, 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 θ₀ and θ₁ (Harpending, 1994; Marini et al., 2021), the estimated effective female population size for S. leptolepis after expansion (θ₁) was approximately 5,000 times higher than prior (θ₀).

FIGURE 4

Figure 4. Mismatch distribution analysis for S. leptolepis mtDNA cytochrome oxidase I (COI) (left) and cytochrome B (CytB) (right) using DnaSP v6.12. Graphs of the mismatch distributions of (A,B) overall populations; (C,D) Asian lineage; and (E,F) Australian lineage. The X-axis shows the observed distribution of pairwise nucleotide differences while Y-axis shows the frequencies. Solid lines represent the observed frequency of pairwise differences and broken lines show the expected values under the sudden expansion model. Sum of squared deviation (SSD) and Harpending’s raggedness index (Hri) indices with corresponding p-values are also presented.

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

FIGURE 5

Figure 5. Bayesian skyline plot based on (A) cytochrome oxidase I (COI) and (B) cytochrome B (CytB) for the effective population size changes throughout time. X-axes represent time in thousand years ago (KYA), while Y-axes show the effective population size (N_e). The blue line represents the median, while the thick blue band represents 95% highest posterior density (HPD) intervals.

Discussion

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.

Genetic Structuring

Differentiation of Asia and Australia Groups

Our findings coherently supported the distinction of the Australian group. The high, significant genetic variations and correlation with geography supported an allopatric separation. The lineages’ pairwise K2P divergences (COI = 6.77%; CytB = 6.64%) were higher than intraspecific variances on carangids with detected cryptic populations: Atule mate COI < 4.82%; S. crumenophthalmus COI < 4.66%; Seriolina nigrofasciata COI < 4.32% (Mat Jaafar et al., 2012); Decapterus maruadsi CytB < 5.0% (Jamaludin et al., 2020). These values were also more than 10× higher than the mean intraspecific divergences of barcoded fishes: COI = 0.39% (Ward et al., 2005); COI = 0.34% (Thu et al., 2019); COI = 0.32% (Zhang, 2011); COI = 0.59% and CytB = 0.57% (Kochzius et al., 2010); and CytB < 1.0% (Li et al., 2018). Herein, it is suggested that the two lineages represented different putative species in an S. leptolepis complex. This magnitude of divergence was at comparable values with the interspecific differentiation in the confamilial Pseudocaranx dentex complex (Bearham et al., 2020).

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_STs, 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 sub-clustering 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 (200–40 KYA). 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 in-between 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).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank our following contacts, especially for providing tissue sample collections: A. Graham (CSIRO National Research Collections, Australia), H. Motomura (Kagoshima University Museum, Japan), M. Nakae (National Museum of Nature and Science, Japan), A. Murase (University of Miyazaki, Japan), and S. P. Huang (Academia Sinica, Taiwan). We also thank A. C. J. Hurboda, T. A. B. Itong, P. H. T. Apag, O. P. C. Halasan, T. N. A. Mat Jaafar, Y. H. Guo, S. Yo, and J. W. Kuo for rendering technical assistance.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.756163/full#supplementary-material

References

Abdussamad, E. M., Rohit, P., Koya, S. K. P., Habeeb Mohamed, O. M. M. J., and Jeyabalan, K. (2013). Carangids (Family: Carangidae) in the seas around Indian subcontinent with description of macro-taxonomic characters for the field identification of genera and species. Indian J. Fish. 60, 21–36.

Google Scholar

Allen, G. R., and Erdmann, M. V. (2012). Reef Fishes of the East Indies: Volumes I-III, Tropical Reef Research. Perth: Universitiy of Hawai’i Press.

Google Scholar

Arfat, Y. A., and Benjakul, S. (2012). Gelling characteristics of surimi from yellow stripe trevally (Selaroides leptolepis). Int. Aquat. Res. 4:5. doi: 10.1186/2008-6970-4-5