Sampling around major landing sites and local markets was conducted in both East and Peninsular Malaysia where 134 samples of various sizes were collected randomly from five wild populations of T. jarbua. Populations were provisionally divided into five groups according to region: 1) Kuala Selangor (KS, N = 31) of west Peninsular which is surrounded by the Straits of Malacca; 2) Kuantan, Pahang (KN, N = 30) of east Peninsular which is adjacent to the South China Sea; 3) Mukah, Sarawak (MH, N = 21) of East Malaysia which is surrounded by the South China Sea; 4) Sandakan (SN, N = 28) and 5) Tawau (TW, N = 24) of East Malaysia which are surrounded by the Sulu Sea and the Celebes Sea, respectively (Fig. 1). Samples were collected within the period of April 2015 to August 2018. Approximately 20 mg of muscle tissue from each fish sample was removed and immediately preserved in 95% ethanol and stored at -20°C until genetic analysis was performed.

Genomic DNA was extracted using 10% Chelex Resin following the protocol of Hyde et al. (2005). Approximately 680 bp of the COI-5’ gene was amplified using the FishF1 or FishF2 forward primers and FishR1 or FishR2 reverse primer pairs (Ward et al. 2005). Polymerase Chain Reaction (PCR) amplification of approximately 1000 bp from the 5’-end of the Cyt b gene was performed using the primer pairs Glu31 and Thr33 (Liu et al. 2015) and internal primers: Glu231 (5’-CTT ACA GGC CTC TTT CTG GCC AT- 3’) and Thr233 (5’- TTT GAG CTA CTA ATG CAG TAT- 3’) were designed for this study.

PCR was performed using a Mastercycler epgradient S thermalcycler (Eppendorf, Hamburg, Germany) and 25μl reaction mixtures consisting of 12.5μl exTEN 2X PCR master mix (1^st BASE, Selangor, Malaysia), 9.5 μl of sterile distilled water, 1μl each of forward and reverse primers, and 1μl of DNA template. PCR cycling conditions were as follow: initial denaturation for 1 min at 96°C, 36 cycles of denaturation at 95°C for 30 s, annealing for 30 s at 44°C (CO1) or 48°C (Cyt b), elongation for 1 min at 72°C, and final elongation for 10 min at 72°C. The amplicons were checked for correct length via electrophoresis on a 1% agarose gel (90V for 25 min). PCR products were sent to Apical Scientific Sdn. Bhd. (Selangor, Malaysia) for purification and DNA sequencing.

Multiple sequence alignment was first performed separately for each gene region using the CLUSTAL X (Thompson et al. 1994) program implemented in BIOEDIT ver. 7.0.5 (Hall et al. 2011). The sequences were subsequently trimmed and aligned manually prior to concatenation of COI and Cyt b sequences. Analyses performed in this study were based on the final truncated length of 1446-bp concatenated sequences. All haplotype sequences were deposited in Genbank under the accession numbers MN529663–MN52993.

Unique haplotypes were quantified and the genetic diversity, nucleotide diversity, and pairwise distance were calculated using DNASP v. 4.0 (Rozas et al. 2003). The level of gene flow among populations (Nm) based on Hudson et al. (1992) was also calculated in DNASP v. 4.0. Analysis of molecular variance (AMOVA) was performed using ARLEQUIN v.3.5 (Excoffier and Lischer 2010) for the four, hypothetical, region-based groupings (Selangor, Pahang, Sabah and Sarawak) to investigate the partition of genetic variation among regions (F_CT), among populations within regions (F_SC), and within populations (F_ST). The significance of the F-statistics for population comparisons was assessed using 1000 permutations. The Tamura Nei plus gamma rate model (TN93+G) was selected by MEGA v. 7.0 (Kumar et al. 2016) as the best-fitting substitution model based on the Bayesian information criterion. A Maximum Likelihood (ML) tree was reconstructed in MEGA 7.0 to show the level of divergence and relationships among haplotypes of T. jarbua. The confidence level at each node was assessed by 1000 bootstrap replications. This tree was compared against the median-joining network generated using the program’s default settings of NETWORK 4.5.0.2 (Bandelt et al. 1999).

In addition, a neutrality test of the pairwise differences among all populations was performed to infer historical demographic and deviation of sequence variation from evolutionary neutrality. Deviations from neutrality were evaluated using Fu’s Fs (Fu 1997) and Tajima’s D (Tajima 1989) via DNASP. Statistical tests and confidence intervals for D and F’s were based on a coalescent simulation algorithm. A large negative value of Fu’s Fs or the Tajima’s D rejecting the null hypothesis of neutrality indicates population expansion(s). The demographic changes were also examined using the mismatch distribution analysis (Rogers and Harpending 1992) in ARLEQUIN with 1000 permutations. The Harpending’s raggedness index (Harpending et al. 1993) and the sum of squared deviations (SSD) between observed and expected mismatch for each of the populations under the model of constant population size were analyzed according to Schneider and Excoffier (1999). This method quantifies the smoothness of the observed mismatch distribution and a non-significant result indicates an expanding population. The spatial expansion hypothesis (both raggedness index and SSD) was tested using a parametric bootstrap approach with 1000 replicates.

The 1446 bp concatenated COI (631 bp) and Cyt b (815 bp) sequences were analyzed for 134 individuals obtained in five different locations (Fig. 1) from East Malaysia and Peninsular Malaysia. The nucleotide composition was 23.0% adenine, 28.9% thymine, 32.0% cytosine and 16.1% guanine. The higher A+T content (51.9%) compared to G+C content (48.1%), is common among fishes (Appendix 1). There were 202 polymorphic sites of which 31 (15.4%) were singleton variable sites, and 171 (84.6%) were parsimony informative (Appendix 2). A total of 149 mutations with 66 transitions, 8 transversions and 75 substitutions were found within the dataset (data not shown).

A total of 83 putative haplotypes were derived from the 134 individuals sequenced with 79 of them being unique haplotypes (95.18%) and four were shared haplotypes (4.82%). The dominant haplotype of Malaysian populations is Hap5 (KS, KN, TW, SN, MH, TAI) while other shared haplotypes are Hap3 (KS, KN, TW, MH, IND), Hap8 (KS and KN) and Hap45 (TW and SN). The population from KS recorded the highest total number of haplotypes (22) of which 19 were unique haplotypes, while Tawau recorded the lowest number of haplotypes (15) with 12 unique haplotypes. The nucleotide diversity (π) of T. jarbua populations in this study ranged from 0.0288 ±0.0158 (mean ±SD) to 0.3434 ±0.1722 while haplotype diversity (h) ranged from 0.9167 ±0.0482 to 0.9952 ±0.0165 (Table 1). The MH population recorded the highest π and h. A high h and low π indicate that the populations studied were moderate in genetic diversity.

A ML tree was reconstructed based on the 83 haplotypes of this study and four COI + Cyt b sequences from Hainan, Taiwan, India and Philippines which were downloaded from National Center for Biotechnology Information (NCBI) (Appendix 3). The mtDNA concatenated dataset defined the haplotypes into three major clades with no significant clusters corresponding to sampling localities (Fig. 2). Apart from Clade I which consists of eight haplotypes solely from Sarawak, the other four clades include mixtures of haplotypes from various localities without any obvious geographical structuring among them. Clade II consists of haplotypes from MH and KS populations while Clade III consists of the haplotype from the Philippines and haplotypes from Malaysia except MH. Clade III is the most geographically inclusive with haplotypes from India, Hainan, Taiwan, Philippines and all five Malaysian populations.

Figure 2. 

Maximum likelihood haplotype tree reconstructed based on the concatenated mtDNA dataset. The bootstrap values higher than 50% are shown near the nodes.

The general topology of the median-joining network (Fig. 3) corresponded with the ML tree (Fig. 2) with three major clusters identified. In the network presented, shared haplotypes occupy the central area, while the unique haplotypes branched out from the center. This formation provides a star-like profile, which indicates population expansions. Hap3 and Hap5 are the dominant haplotypes in cluster III. The distribution frequency of all the 83 haplotypes in T. jarbua populations is presented in Appendix 4. Hap5 recorded the highest distribution frequency with 21 individuals. It is the only common haplotype shared between all five populations from Malaysia and the representative from Taiwan. Hap3 with 16 individuals was found in India and in all locations sampled in Malaysia except for Sandakan.

Figure 3. 

Haplotypes median-joining network corresponding to the ML tree with three observed clusters. The star-like profile observed in cluster III indicates the presence of sudden expansion.

Pairwise F_ST comparisons between populations in Malaysia were significant at the 95% confidence level except for the comparison between TW and SN (Table 2). Populations of MH and KN showed the greatest pairwise differentiation (F_ST = 0.5353; p < 0.05) while SN and TW showed the least differentiation (F_ST = 0.0452; p > 0.05). The pairwise nucleotide divergence among populations (Table 3) showed the same trend as the F_ST values and was not correlated with geographical distance. The overall gene flow (Nm) estimated among populations was low at 0.82. The sequence divergence was calculated using the Kimura two parameter (K2P) distance model for both genes (Table 3). The greatest genetic differences (COI: 0.019 and Cyt b: 0.029) were observed between MH-KN, MH-SN and MH-TW. The T. jarbua populations displayed a low level of conspecific divergence within 2% (COI).

The genetic structure of the T. jarbua populations analysed by AMOVA showed little (39.52%) genetic differentiation among regions but high (62.13%) variation within populations (Table 4). This indicates that the populations were not genetically differentiated among regions and the genetic variation was mainly from within the population level. There is essentially no genetic structuring (-0.14% variation) among populations within region.

The overall Tajima’s D value was negative with an insignificant p-value, indicating deviation from evolutionary neutrality. Similarly, the Fu’s Fs test which is based on the distribution of haplotypes, revealed negative but significant p-values for all five populations studied, indicating an excess of rare haplotypes or rare mutations in the population compared to what is expected under a neutral model of evolution. Following the results of Fu’s Fs test, the hypothesis of neutral evolution was rejected.

In the present study, all populations demonstrated bimodal and ragged shaped patterns which points to the population having remained largely constant in size and that the lineage was widespread (Rogers and Harpending 1992). The scatterplot of the bimodal illustrations is shown in Figure 4. The results of mismatch distribution are contradictory to the results of the neutrality analysis. Hence, to further test the validity of the neutrality test results, we calculated the raggedness index and SSD under the demographic expansion model as shown in Table 5. P-values of SSD between the observed and expected mismatch distributions were all statistically insignificant (p > 0.10), indicating the presence of non-equilibrium and a population expansion event in T. jarbua. Besides that, the studied populations showed non-significant raggedness index (p > 0.10) indicating the data has relatively good fit to a model of population expansion (Harpending 1994). A ragged distribution suggests that the lineage was widespread (Excoffier et al. 1992; Rogers and Harpending 1992; Rogers 1995).

Figure 4. 

Pairwise number of difference (mismatch distribution) analysis was conducted using the constant population size model to observe the population size changes. The observed frequencies were represented by red dotted line. The frequency expected under the hypothesis of population expansion model was depicted by continuous green line. a Kuala Selangor b Kuantan c Mukah d Sandakan e Tawau f all populations.

Species identification was confirmed by morphological observation and DNA sequence data in which the intraspecific COI divergence was within the 2% threshold value (Hebert et al. 2003; Ward 2009; Srivathsan and Meier 2012). In general, a relatively high haplotype diversity (0.92–0.99) and low nucleotide diversity (0.03–0.34) were observed for the populations in this study. The combination of high haplotype diversity and low nucleotide diversity is common in pelagic marine fishes (Liu et al. 2015). This is likely due to rapid demographic expansion from a small effective population size, assuming there is sufficient time for the number of haplotypes to increase through mutation but insufficient for accumulation of large sequence differences (Grant and Bowen 1998; Avise 2000; Lowe et al. 2004). The inference of population expansion is further supported by the star-like patterns in Figure 3. Our results agree with an earlier study of T. jarbua populations from the Taiwanese waters which exhibited similarly high population haplotype diversity ranging from 0.86 to 1.00 as inferred from COI and Cyt b sequence data (Liu et al. 2015). In a study on T. jarbua populations from the Gulf of Aden, Lavergne et al. (2014) reported similarly high genetic diversities ranging from 0.216 to 0.698. The high number of haplotypes (83 haplotypes) in the present study is likely due to the high mutation rate of the mtDNA genes. Most of the haplotypes are unique to its region which may indicate the presence of different founding populations in the studied localities (Teixeira et al. 2011). The haplotype diversity of a relatively rapid evolving genome within a population often approaches 1.0 as many individuals will tend to have unique haplotypes (Freeland et al. 2011).

A population’s genetic structure is affected by genetic drift, local adaptation, and gene flow. In a marine environment, the development of population structure is greatly influenced by factors that affect dispersal, such as ocean currents, historical variance, and geographic distance coupled with differences in dispersal ability and habitat discontinuity (Saarman et al. 2010). Population structure inferred from mtDNA markers displays less genetic divergence in the pelagic and moderately pelagic species due to their potential to undertake long-distance migrations in oceanic waters (Jaafar 2014).

The haplotype tree (Fig. 2) revealed three major lineages but geographic structuring among the five populations is not distinct. In general, haplotypes specific to certain geographic regions did not form monophyletic groups, but appeared to be randomly distributed across the haplotype tree. Hap5 and Hap3 which recorded the highest distribution frequency are likely the ancestral haplotypes among the populations sampled. Recent haplotypes were evolved directly or indirectly from the ancestral haplotypes. The existence of two ancestral points indicates that T. jarbua in Malaysia probably exist from two different sources. According to the coalescent theory, common haplotypes at the center of a network are inferred to be ancestral, while tip haplotypes at the periphery are derived or descendant from ancestral haplotypes (Akib et al. 2015). The occurrence of star-like patterns radiating from these major haplotypes suggests that T. jarbua populations have undergone significant population size expansions in the relatively recent past (Forster et al. 2001; Akib et al. 2015).

F_ST values are often used to infer gene flow, in which a lower F_ST value indicates low genetic divergence and higher gene flow. F_ST values below 0.05, as observed between SN and TW populations, indicate negligible genetic divergence, probably due to active exchange of genetic material between populations through breeding. Furthermore, the pairwise divergence between these populations is not statistically significant. According to Wright (1965), F_ST of 0–0.05 is described as little differentiation, 0.05–0.15 as moderate differentiation, 0.15–0.25 as great differentiation and values greater than 0.25 as very great differentiation. All populations studied showed moderate to very great pairwise differentiation except for TW-SN. The overall gene flow recorded was rather low (Nm = 0.82) which suggests limited genetic connectivity among the five populations.

Populations from the same region, i.e., TW and SN of Sabah, were the least genetically variable (Tables 2 and 3), which is likely due to the close geographical distance between these populations. The theory of fish migration across adjacent drainage systems due to flooding, which follows the one-dimensional stepping stone model that allows migration to adjacent population (Song et al. 2013), may apply in the case of T. jarbua. The significantly higher genetic differentiation between populations of KS-SN (pairwise F_ST = 0.2681) compared to KN-SN (pairwise F_ST = 0.0702) may be attributed to distance and physical barrier. Some genetic exchange can be expected since the Straits of Malacca connects the Andaman Sea and South China Sea via the narrow Tebrau strait. It is likely that mixing between the two bodies of water is very limited which supports the F_ST value obtained. However, the higher genetic variation between populations of MH-SN (pairwise F_ST = 0.5038) as compared to MH-KS (pairwise F_ST = 0.3310) implies that geographical distance is not the only driving factor of genetic variation among populations of T. jarbua in the Malaysian waters, similar to the results observed in the wider Gulf of Aden (Lavergne et al. 2014) where populations of adjacent locations showed low genetic connectivity despite the absence of a geographic barrier. Populations bordering a common origin such as the South China Sea (KN, MH, TW, and SN) may have evolved independently of each other over time, but there might have been insufficient time for genetic divergence to accumulate in these populations.

Another interesting finding of this study is the occurrence of shared haplotype between the populations from Peninsular and East Malaysia, India, Hainan, Philippines and Taiwan. Common haplotypes between localities and mixed haplotypes of different lineages in some populations in the current study can be explained by the biogeographical history of Southeast Asia (historically known as the Sundaland). Southeast Asia is believed to have experienced simultaneous glaciation and consequent deglaciation along with its associated decrease and increase of seawater levels during the Pleistocene period, which greatly influenced continental and oceanic configuration (Voris 2000). The shared haplotypes between Malaysian populations and those from as far as India suggests that the range of population expansion after glacial retreat was not restricted to the South China Sea but also extended into the Indian Ocean (Liu et al. 2015). Lavergne et al. (2014) also reported high connectivity between populations in the Gulf of Aden and South China Sea due to the unique sharing of COI haplotypes between both regions. The haplotype sharing and their consequent gene flow may also be attributed to breeding migration, mutation, pelagic larvae, and sharing of common ancestors (Frankham 1996).

The MH population is the most genetically distinct with the highest between-group mean distances, haplotype and nucleotide diversity among the five populations. Geographical isolation of allopatric populations restricts gene flow between two populations, which in turn allows the evolution of a genome adapted to local condition (Hall 1993). Cluster I (MH) is estimated to form after the separation of Peninsular Malaysia from the Borneo Island due to the rise in the depth of the Sunda River between 40 to 100 m. This gradual separation was suspected to have caused accumulative genetic drift. According to Halliday (1993), genetic drift is likely to occur, particularly in small populations that are isolated from the main population and it may become the major source of genetic variation between some populations.

Among the four populations, MH is genetically closest to KS. Geological evidence suggests that the river systems of Sarawak were historically interconnected with most major river systems of Peninsular Malaysia via the Sunda River during Pleistocene glaciation (about 10000 years ago), thus allowing gene flow among these drainages (Kamarudin and Esa 2009). Gene flow from populations in the Straits of Malacca to those in Sarawak has been reported in several studies including Ryan and Esa (2006), Azhar and Hassan (2015), Samani et al. (2016) and Lau et al. (2018). Meanwhile, populations of SN-KN, which are separated by the South China Sea, showed high genetic connectivity (pairwise F_ST = 0.0702). This could probably be explained by the high migration ability of T. jarbua (>1000 km, Liu et al. 2015), human-mediated transfer through ballast waters (Liu et al. 2019) or past glaciation events. Furthermore, the high similarity in the sequence data (Table 3) perhaps indicates remnants of identical haplotypes from both populations, and that they were essentially similar at one time before the separation (Inger and Chin 2002). Eventually, sea level rise during the last Pleistocene period caused Borneo to be separated from mainland Asia (Peninsular Malaysia), which we suggest, resulted in shelf submergence and subsequent genetic differentiation between grunters from KN and SN. Pleistocene sea level fluctuations could also explain the incomplete divergence of grunters between East and Peninsular Malaysia. Similar evidence of a close genetic relationship between fishes of Borneo and mainland Asia in relation to their biogeographical history was discussed by several other authors (Pin et al. 2001; Nadiatul et al. 2011; Tan et al. 2012; Song et al. 2013).

Historical demographic expansions were determined by analysing the frequency distributions of pairwise differences between sequences (Rogers and Harpending 1992; Ray et al. 2003; Excoffier 2004). Neutrality tests with Tajima’s D and Fu’s Fs statistics estimate the deviation from neutrality, which is based on the expectation of a constant population size at mutation-drift equilibrium. Here, a negative Tajima’s D signifies an excess of low frequency polymorphisms relative to expectation, indicating population size expansion or positive selection (Tajima 1983). The negative and significant Fu’s Fs statistical value provides strong evidence for past population expansion, and rule out the possibility of genetic hitching or background selection, and evolutionary forces that produce a pattern similar to population expansion (Fu and Li 1993; Fu 1997; Okello et al. 2005). The T. jarbua populations displayed a genetic pattern typical of a population that has undergone a recent population expansion due to its two common haplotypes (Hap3 and Hap5) present across the range while the rest of the haplotypes are unique. The range expansion was a recent phenomenon and may not have achieved the migration-drift equilibrium, as shown by the lack of phylogeographical structure. Neutrality test statistics were in overall negatively significant and not consistent with a population at drift-mutation equilibrium.

The mismatch distribution is generally displayed as a multimodal pattern for populations showing demographic equilibrium. In contrast, a unimodal pattern depicts populations which have experienced recent demographic expansion (Rogers and Harpending 1992). In the results, all localities presented a multimodal pattern proving recent expansion. The hypothesis that the observed data fit the sudden expansion model was tested using the SSD and the raggedness index. Here, non-significant values for SSD signifies that the observed data do not deviate from that expected under the model of expansion. Non-significant raggedness index also indicates population expansion. Our observations of non-significant values in goodness-of-fit distribution for all populations suggest that population expansion occurred recently (Rogers and Harpending 1992).