Mitochondrial marker implies fishery separate management units for spotted sardinella, Amblygaster sirm (Walbaum, 1792) populations in the South China Sea and the Andaman Sea

Ethical statement

Since only dead specimens were sampled, no permit was required, and no ethical consideration were required for this study. In addition, A. sirm is not listed under the International Union for Conservation of Nature (IUCN) list of endangered or protected species.

Sample collection

A total of 179 specimens of A. sirm were collected from ten landing sites from vessels that used purse seine fishing gears in Peninsular Malaysia, Sabah, and Sarawak from 2015 to 2019. These sampling locations cover the southern region of the SCS, SS, CS, and AS. The species were morphologically identified by of 10 to 20 series of black spots down the flank, according to Whitehead (1985) (Table 1; Fig. 1).

In addition, a small clipping of A. sirm pectoral fin was collected during the sampling activity. Each clipping was fixed in a separate vial containing 95% ethanol and stored at −20 °C until further use.

The sampling locations were divided into four regions following Kasim et al. (2020). The details of each landing port were obtained via interviews with fishermens/boatmen and confirmed by consulting with the Department of Fisheries staff stationed at the landing sites. Notably, no samples were obtained from the Strait of Malacca (a waterway that links the AS and the SCS), as A. sirm is absent from this location (Whitehead, 1985). Meanwhile, sample collections in Ranong, Thailand, were facilitated by the Andaman Coastal Research Station for Development, Kasetsart University of Thailand. Two samples from Taiwan were contributed by collaborators from the National University of Taiwan and the National Museum of Marine Biology, Taiwan (see Fig. 1). Due to the low number of samples (N = 2), the data was only used in the phylogenetic analysis.

Extraction and polymerase chain reaction (PCR) amplification

Genomic DNA extraction of 179 samples was performed using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. The quality of extracted DNA was assessed using UV spectrophotometer Q3000 (Quawell, San Jose, CA) and diluted to a final concentration of 50 ng/µL. The mtDNA Cyt b was targeted using the following primers; F: WMA15-F (5′ACC GTT GTA ATT CAA CTA TAG AAA C 3′) and R: TruCytb-R (5′ CCG ACT TCC GGA TTA CAA GAC CG 3′) (Jérôme et al., 2003). The PCR amplification was conducted in a Thermal Cycler (Techgene, Irving, TX, USA) at a final volume of 25 µL, composed of 12.5 µL of 2x EasyTaq^® SuperMix (TransGen Biotech, Beijing, China), 0.2 µL of each primer, 2.0 µL of genomic DNA (50 ng/mL), followed by nuclease-free water to achieve the final reaction volume. The thermal cycling conditions were as follows: initial denaturation at 95 °C for 2 min, followed by 35 cycles of denaturation at 94 °C for 10 s, annealing at 56.8 °C for 10 s, and extension at 72 °C for 15 s. The quality of PCR products was visualised on 1.5% agarose gels stained with 1 to 3 µL of GelRedTM Nucleic Acid Gel Stain (Biotium Inc., Fremont, CA, USA). The unpurified PCR amplicons were sent to Repfon Technologies Sdn Bhd for purification and sequencing in forward and reverse directions in an automated sequencer (ABI3730xl, Applied Biosystems, Bedford, MA, USA).

Data analysis

The raw sequences for both forward and reverse sequences were edited and assembled after removing the primer sites.The final assembled sequences were aligned with ClustalW implemented in MEGA v7.0 with default setting (Kumar et al., 2018). The complete aligned data set was inspected for nucleotide variable sites, parsimony informative sites, number of haplotypes, haplotype distributions, transitions and transversions, and nucleotide frequencies in DnaSP v6.0 (Rozas et al., 2017). A simple linear regression analysis (Pearson correlation test) was conducted to assess whether sample size (N) would affect the downstream analyses to ensure the statistical validity of further investigations. All the statistical analysis were performed using the IBM SPSS Statistics for Windows software, version 23.0 (IBM Corp. Armonk, NY, USA). Since the initial analysis demonstrated that the sample sizes did not affect the haplotype diversity (Hd) and nucleotide diversity (π) (p > 0.05) therefore, further analysis could be proceeded (Castillo-Páez et al., 2014; Delrieu-Trottin, Maynard & Planes, 2014).

Phylogenetic and genetic diversity analyses

The RAxML (Stamatakis, Hoover & Rougemont, 2008) and maximum likelihood (ML) trees, with 1,000 bootstrap replicates, were constructed in raxmlGUI v1.5b1 (Silvestro & Michalak, 2012) using a general time-reversible model of nucleotide substitution with a heterogeneity rate following a discrete gamma distribution (GTR+G) as the default selection. Bayesian inference (BI) was carried out in MrBayes v.3.2 (Ronquist & Huelsenbeck, 2003) with 1,000,000 Markov Chain Monte Carlo (MCMC) run generations; posterior distribution was sampled every 1,000 generations, and a 25% burn-in. GenBank sequences of S. hualensis (Acc No: KC951523), S. lemuru (NC039553), and S. gibbosa (NC037131) were employed as outgroups in the analysis for both trees. A phylogenetic network of all haplotypes was constructed to view haplotype relationships based on the median-joining calculation in a minimum spanning network (MSN) implemented in the population analysis with reticulate trees (POPART) v1.7 (Bandelt, Forster & Röhl, 1999; Leigh & Bryant, 2015).

Genetic distances within and among populations based on the best nucleotide substitution model of Kimura 2P (K2P), which depicts the lowest Bayesian information criterion (BIC) score, were estimated using MEGA v7.0. These values also assessed the possibility of sub-species or cryptic species occurrence if intra-species variation within marine species exceeds the threshold (2%) (Chanthran et al., 2020; Hebert et al., 2003; Jamaludin et al., 2020; Ward, 2009).

Population structure

Population pairwise comparisons (F_ST) were calculated using ARLEQUIN v3.5 (Excoffier & Lischer, 2010). Values were adjusted for Type 1 errors as a result of multiple comparisons using the false discovery rate procedure (FDR) at p < 0.05 (Benjamini & Hochberg, 1995). The genetic distance and significance of each pairwise comparison were further analysed with a nonparametric permutation procedure with 1,000 replicates (Hudson, Boos & Kaplan, 1992) in DnaSP v6.0 (Rozas et al., 2017). Estimation of gene flow (Nm) based on both haplotype and sequence statistics was derived according to Nei (1973) and Hudson, Boos & Kaplan (1992) using the same programme.

Hierarchical analysis of molecular variance (AMOVA) was performed to estimate molecular variance among populations at different hierarchical levels using ARLEQUIN v3.5 (Excoffier & Lischer, 2010). The relative contribution of variances was estimated at three different levels; F_ST, F_SC, and F_CT (Excoffier & Lischer, 2010). The spatial structure was further examined using spatial analysis of molecular variance (SAMOVA) v2.0 (Dupanloup, Schneider & Excoffier, 2002) to identify groups of populations that were geographically homogeneous and maximally differentiated from each other. The isolation by distance (IBD) or Mantel test (Mantel, 1967) was conducted to confirm the relationship between genetic distance and geographical distance.

Demographic history

The historical demography of each population was tested using Tajima’s D (Tajima, 1989) and Fu’s F (Fu, 1997) statistics in ARLEQUIN v3.5 (Excoffier & Lischer, 2010) Ramos-Onsins and Rozas’ R₂ in DnaSP v6.0 (Rozas et al., 2017) to analyse deviation from neutrality. The significant R₂ was calculated using coalescent simulations of 5,000 replicate runs for each simulation. R₂ is a powerful tool for quantifying population growth with a limited sample size (Ramos-Onsins & Rozas, 2002). Harpending’s raggedness index, Hri (Harpending, 1994), was calculated in ARLEQUIN v3.5 (Excoffier & Lischer, 2010) and mismatch distributions (Rogers & Harpending, 1992; Schneider & Excoffier, 1999; Slatkin & Hudson, 1991) were calculated in DnaSP v6.0 (Rozas et al., 2017). Both analyses could differentiate whether populations are demographically stable, expanding, or decreasing over time. A sudden population expansion model in Hri will be rejected if p < 0.05 (Schneider & Excoffier, 1999).

The precise time of expansion was manually calculated using the equation τ = 2 µt, (µ= mutation rate of the sequence analysed, t = time since expansion) with a given mutation rate of 1–2% per million years for Cyt b (Rousset, 1997; Johns & Avise, 1998; Cárdenas et al., 2005). The Bayesian skyline plot (Drummond et al., 2012) was implemented in BEAST v1.8.2 (Drummond et al., 2012) to estimate past population dynamics, employing a relaxed uncorrelated lognormal molecular clock, GTR + G as the best evolutionary model selected from PartitionFinder v1.1.0 (Lanfear et al., 2017). The analysis was run for 300 million generations with parameters sampled every 10,000 generations. The results were visualised using Tracer v1.7 (Rambaut et al., 2018) which summaries the posterior distribution of population size over time.

Data analysis

The 1,016 bp segment of mtDNA Cyt b was successfully amplified from all 179 samples. A detailed evaluation revealed 42 non-synonymous mutations resulting in 37 amino acid substitutions. The ratio of transversion to transition substitutions for the entire data set was 1:1. Haplotype diversity (Hd) was high ranging from 0.91 to 1.00, while nucleotide diversity (π) was low, ranging from 0.002 to 0.009 (see Table 1). A total of 203 variable sites were identified from the 1,016 bp segment, where 111 (10.9%) parsimony informative sites defined 122 haplotypes. These haplotype sequences have been deposited in GenBank under accession numbers MZ040756 to MZ040877. Furthermore, 111/122 (90.98%) were detected as singletons (haplotypes found in a single gene copy), where 11/122 (9.02%) were shared haplotypes or found in more than a single population. The Pearson correlation test (r² = 0.66, p = 0.06) confirmed that sample sizes did not significantly affect the nucleotide diversity and haplotype diversity (Hd) therefore, further analyses could be performed.

Phylogenetic relationships

The ML tree and BI analysis revealed a similar topology of two geographically separate mtDNA lineages (see Fig. 2); (1) AS lineage (Lineage 1) composed of two populations (Ranong and Kuala Perlis), and (2) SCS lineage (Lineage 2) encompasses populations from the SCS, CS, and the SS. No geographical pattern was observed within Lineage 2.

The minimum spanning network (MSN) showed two distinct lineages similar to the phylogenetic trees, separated by 63 genetic mutations (see Fig. 3). Lineage 2 (all populations excluding AS populations) exhibited a star-like pattern with the dominant haplotype (Hap05) centred in it (Ferreri, Qu & Han, 2011; Woolley, Posada & Crandall, 2008).

Genetic diversity

The genetic diversity or distance was converted to percentages and revealed intrapopulation genetic diversity ranging from 0.2% to 0.9%, the lowest being in Kudat and the highest in Kuala Perlis (Table 2). Meanwhile, the interpopulation genetic diversity ranged between 0.2% and 7.6%, with Kuala Perlis and Ranong from the AS being the most divergent of all (7.2% to 7.6%), with a 1.0% difference from each other. In contrast, the genetic diversity ranged from 0.2% to 0.7%, when Kuala Perlis and Ranong populations were excluded from the analysis, providing strong evidence that AS populations are highly structured than other seas.

Population structure

Most F_ST values in the AS (Kuala Perlis and Ranong) populations were significantly different (p < 0.05) (Table 2), which aligned with earlier studies. However, non-significant (p > 0.05) values were also detected in cases where high genetic distance in some between Lineage 1 and Lineage 2 populations, such as AS vs Semporna Kuantan vs Kuching (both in SCS), Labuan (SCS) vs Kudat (SS), and Semporna (CS) vs Kuantan (SCS).

The AMOVA indicated no significant differences among populations (between seas) even for the highest variance (F_CT = 0.0066, p = 0.257; Table 3). However, the dataset demonstrated significant differences within each sea (F_ST = 0.0265, p = 0.0009), with 97.34% of the total genetic variation of A. sirm contributed by genetic differences among the total populations. This result indicated a subdivision or structuring within the tested populations. Notably, populations may be contribute the non-significant (p > 0.05) observation within the SCS, represented by more than two populations compared to other seas. Genetic differentiation among populations within groups/seas was significant (F_SC = 0.0202, p = 0.03).

The result from the SAMOVA paralleled other analyses, with k = 2 displaying the highest F_CT value (F_CT = 0.9277, p = 0.02) signifying two genetically distinct A. sirm stocks among sampled populations: Lineage 1 (Kuala Perlis and Ranong) and Lineage 2 (Kuching, Kota Kinabalu, Kuantan, Labuan, Pulau Kambing, Kudat and Semporna) (Supplemental Information 1). The Mantel test showed no significant correlation between genetic differentiation (F_ST value) and geographical distance (r = 0.127, p = 0.173) among tested populations.

Demographic history

The Tajima D’s and Fu’s F generated negative and significant (p < 0.05) values for all populations (Table 1) suggesting a recent historical directional selection (selective sweep) or recent population growth (Tajima, 1989). Furthermore, the Harpending raggedness index, Hri = 0.042, SSD = 0.01, (p = 0.3) supported the recent population growth hypothesis. Meanwhile, the mismatch distributions of combined populations detected two highly divergent peaks, representing the two lineages (see Fig. 4A). Additional independent analyses of the two lineages (Lineage 1 and Lineage 2) were conducted, and a bimodal mismatch distribution was discovered in Lineage 1 (Fig. 4B). Nevertheless, further comprehensive studies are required to confirm these preliminary finding due to the small sample size in this study. On the other hand, Lineage 2 shows a unimodal distribution providing strong evidence for sudden population expansion (Harpending, 1994; Slatkin & Hudson, 1991) (see Fig. 4C).

The Bayesian skyline plot analyses revealed a recent expansion in the overall effective population size (see Fig. 4D), which can be explained by growth in Lineage 2 occurred approximately around 150 thousand years ago (KYA) (see Fig. 4F). Conversely, Lineage 1 showed a gradual population increase over the last 500 KYA (see Fig. 4E). Notwithstanding, the finding was based on only two populations from the same sea and may not represent other similar sites. The estimated time of population expansion for Lineage 1 and Lineage 2 using the formula τ = 2 μt, generations and a mutation rate for Cyt b of 1–2% per million years for perciform (Cárdenas et al., 2005; Johns & Avise, 1998) are estimated to occur 442,913 to 221,456 years ago and 354,330 to 177,165 years ago respectively. These findings are almost in parallel with the results from the Bayesian Skyline plot suggesting a demographic expansion in the A. sirm population during the Pleistocene era.

Implications to fishery management

This study has provided novel insights into the population structure and proposed several processes that could define the stock boundaries and evolutionary units of A. sirm in waters. It is recommended that separate management systems be implemented for the waters fringing the western coasts, Lineage 1 (AS) and eastern coasts, Lineage 2 (the SCS and the neighboring waters). Despite focusing on only one country, it is worth notingy that all the major seas in Southeast Asia are represented in this study.

Genetic data should be an important component for fishery managers to make informed decisions on a target species (von der Heyden et al., 2014). Furthermore, advanced genetic techniques (eDNA for ecological monitoring) could provide crucial information to the scientific advisory process for fisheries management. Despite that, the integration of genetic data into fishery policy modeling is yet to gain traction in most parts of the world where genetic facilities are still lacking (Reiss et al., 2009; Waples & Naish, 2009), including Malaysia.

Regional cooperation can benefit from utilising genetic evidence, particularly for migratory species like A. sirm. A fine example of how regional bodies co-operate in regional fishery management strategies is the ZoNéCo (Programme d’évaluation des ressources marines de la zone économique de Nouvelle-Calédonie) in New Caledonia that integrated genetic and complementary non-genetic data to manage the Spanish mackerel, Scomberomorus commerson stocks in Bélep, North Province of New Caledonia (von der Heyden et al., 2014). This project demonstrated that the area might host distinct reproductive stocks based on solid evidence of genetic structuring. Resultantly, a sound policy was introduced to ensure sustainable exploitation of the Spanish mackerel in New Caledonia. In another study, high genetic structuring attributed to limited ecological population connectivity of the shorefish, Eleutheronema tetradactylus, in four regions of northern Australia (Horne et al., 2011) recommended separate management of these four populations.

The presence of two discrete A. sirm stocks, the AS and the SCS, reinforces the need for regional cooperation among the maritime nations. This effort could be realized by referring to the management models applied in other regions through the facilitation of the Southeast Asian Fisheries Development Center (SEAFDEC). The genetic variability data could be used as an indicator to control the harvesting activities between participating countries, besides regulating the number of fish landing at a particular time and the fishing capacity of fishing vessels (Md Saleh et al., 2019).

Determining the population stock structure within the region (Malaysia with the neighboring waters) is also crucial for the sustainable fishery management of highly economically important species. This information could guide the fishery managers in determining potential fishing or spawning areas and planning a sustainable fishery policy strategy. For example, the population stock structure determination of spawning areas through DNA metabarcoding could be a way forward in sustaining broodstock management. In addition, the information could facilitate site recommendations for implementing closed areas during the critical life cycle, including spawning and nursery seasons, as suggested by Radhakrishnan et al. (2018) in their study on the Spanish mackerel, Scomberomorus commerson. The concept of genetic-based management is slowly gaining recognition in this region. For instance, SEAFDEC identifies potential closed areas under the “The South China Sea Fisheries Refugia Initiative” whereby genetic and non-genetic information is used to develop refugia of selected species such as the tiger prawn refugia in Sarawak, Malaysia (SEAFDEC, 2006). Despite that, there is an immense gap and challenges in incorporating genetic data for the development of fishery management strategies due to the lack of interest and knowledge among fishery managers (Benestan, 2019; von der Heyden et al., 2014). Furthermore, non-specialist may face difficulties understanding the genetic approaches in fishery management (Bowen et al., 2014), therefore, it is critical to address these challenges by maintaining and strengthening joint efforts between scientists and fishery managers. For instance, the integration of A. sirm genetic data along with life-history traits, tagging returns, parasitic loads, and microchemical variation of the hard skeletal structure in the fishery management policy (Ward, 2000; von der Heyden et al., 2014) would permit a holistic approach in sustainably managing A. sirm in this region. Such efforts will help achieve the Sustainable Development Goals (SDGs) specifically SDG14 to combat the global adverse effects of overfishing by 2030 (United Nations Sustainable Development, 2021) and effective conservation and sustainable fishery management.