Recent population expansion of longtail tuna Thunnus tonggol (Bleeker, 1851) inferred from the mitochondrial DNA markers

Ethical statement

Only a small clipping of the pectoral fin from each individual fish was collected from the local wet markets. This species is not in the IUCN list of endangered or protected species. As only dead specimens were sampled, no permit was required and no ethical consideration was linked to the study.

Sample collection

The specimens were collected from 11 fish landing sites in Malaysia and were morphologically identified following identification keys as described in Lim et al. (2018). This species can be differentiated from other Thunnus species by having a moderate length pectoral fin, reaching the origin of the second dorsal fin and blackish caudal fin (Lim et al., 2018). A small portion of pectoral fin of each individual was cut and preserved in 95% ethanol and stored in a 1.5 mL centrifuge tube at 4–8 °C until further analysis. Each catch locality was confirmed to be non-overlapping (discrete geographical entity) based on feedback from the fishermen and was divided into four regions following Akib et al. (2015) (Table 1; Fig. 1A). Additionally, four samples from Taiwan and 153 GenBank sequences of T. tonggol from the Indian waters (Fig. 1B) were included for phylogeographic analysis.

Genomic DNA isolation and polymerase chain reaction amplification

Genomic DNA was isolated from fin tissue by using the salt extraction method (Miller, Dykes & Polesky, 1988). The isolated DNA samples were then PCR amplified with the mtDNA partial displacement loop (D-loop) and ND5. The following primers were used: (1) D-loop—Pro889U20 (5′-CCW CTA ACT CCC AAA GCT AG-3′, forward) and TDKD1291L21 (5′-CCT GAA ATA GGA ACC AAA TGC-3′, reverse) (Sulaiman & Ovenden, 2009) (2) ND5—L12321-Leu (5′-GGTCTTAGGAACCCAAAACTCTTGCTGCAA-3′, forward) and H13396-ND5 (5′-CCTATTTTKCGGATGTCYTG-3′, reverse) (Ruzainah, 2008). The PCR reaction mixture consisted of 2 µl of the DNA template, 0.5 µl of each primer, 12.5 µl of MyTaq™ Red Mix (Bioline), and 9.5 µl sterilized ultrapure water (ddH₂O) with a final volume in each tube of 25 µl. The temperature profile for the D-loop was: initial denaturation at 94 °C for 5 min followed by 32 cycles of 94 °C for 30 s, 55 °C for 30 s, 72 °C for 1 min, final extension at 72 °C for 5 min, and final hold at 4 °C. Amplification conditions for the ND5 gene were: initial denaturation at 94 °C for 2 min followed by 35 cycles for 94 °C for 20 s, 59 °C for 20 s, 72 °C for 1 min 30 s, final extension at 72 °C for 5 min, and final hold at 10 °C. PCR products were visualized on 1.7% agarose gels stained with SYBR Safe to confirm their presence and estimate the size of DNA fragment amplified. PCR products were then sent for sequencing (First BASE Laboratories Sdn Bhd, Selangor, Malaysia) in the forward direction only using an Applied Biosystem ABI3730x1 capillary-based DNA sequencer.

Sequence editing and alignment

Multiple sequences were aligned and edited using ClustalW implemented in MEGA 6.0 (Tamura et al., 2013). DNA sequences were verified for correct identity by using the Basic Local Alignment Search Tool (BLAST) in the National Center for Biotechnology Information (NCBI) database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) before further analyses. All haplotypes were deposited into GenBank under the accession numbers MK643829–MK644008 (D-loop) and MN252922–MN252981 (ND5).

Data analyses

Demographic history

Historical demographic and spatial expansions were inspected in the T. tonggol populations. Fu’s F_S (Fu, 1997) and Tajima’s D (Tajima, 1989) were adopted to analyze deviation from neutrality. Historical demographic parameters, including the population before expansion (ϴ₀), after expansion (ϴ₁), and relative time since population expansion (τ), were computed in Arlequin 3.5. The values of time (τ) were transformed to estimate the actual time (T) since population expansion, using the equation τ = 2μt, where t is the age of the population in generations and µ is the sequence mutation rate per generation. In the present study, a mutation rate of 3.6 × 10⁻⁸ mutation per site/year was applied for the D-loop (Donaldson & Wilson, 1999) and 2% per million years for the ND5 (Brown, George & Wilson, 1979). Bayesian skyline analyses were plotted using BEAST version 2.2.0 (Bouckaert et al., 2019), where the changes in effective population size (Ne) over time were tested. This enabled past demographic changes of T. tonggol to be inferred from the current patterns of genetic diversity within a population (Drummond et al., 2005). Since there was the absence of a population structure (see “Results”), a single population was modeled. The input was prepared in BEAUti. The analysis was run for 10⁸ iterations with a burn-in of 10⁷ with sampling every 10⁴ and a strict molecular clock. All operators were automatically optimized and the results were generated using Tracer version 1.7.1 (Rambaut et al., 2018).

Harpending’s (1994) raggedness index (Hri) and sum of squared deviations (SSD) were computed in Arlequin 3.5 to evaluate if the sequence data significantly diverged from the assumptions of a population expansion model. The raggedness index has been shown to be a powerful tool in quantifying population growth with limited sample sizes (Ramos-Onsins & Rozas, 2002). In addition, the mismatch distribution was calculated in DnaSP 6.0. The pattern could be used to provide an insight of the past population demography (Chen et al., 2015). A population that has undergone recent expansion shows a unimodal distribution pattern, while a population in equilibrium shows a multimodal distribution pattern (Slatkin & Hudson, 1991; Rogers & Harpending, 1992).

Genetic diversity

A total of 203 and 208 individuals were sequenced for the partial mtDNA D-loop and ND5 gene, respectively. The final alignment of D-loop sequences (416 base pairs (bp)) revealed 113 polymorphic sites (42 singletons and 71 parsimony informative sites), defining 180 haplotypes, where 14 (7.78%) were found in two to five localities and the rest (92.2%) were either private to a single locality or singleton haplotype. The ND5 sequences (855 bp) revealed 60 variable sites (34 parsimony informative sites, 26 singletons), defining 60 haplotypes, where 17 (27.9%) were shared by two to 11 populations, five (8.2%) were private to a single locality and 38 (63.9%) were singleton haplotypes. The D-loop region was AT rich, while, ND5 gene sequences demonstrated higher percentages of CG (56%). All populations of T. tonggol from Malaysian waters showed high haplotype diversity (D-loop: 0.990–1.000; ND5: 0.848–0.965) but low to moderate nucleotide diversity (D-loop: 0.0195–0.0250; ND5: 0.0017–0.0039) (Table 2).

Phylogenetic and population level analyses

The phylogenetic reconstruction inferred from the mtDNA D-loop region and ND5 gene revealed a gene tree with mainly unsupported clades (<50%) and obscure patterns of geographical segregation associated with genetic distribution (Fig. 2; Supplemental 1). This was aligned with the MSN haplotype network that showed no geographical partitioning among the populations studied. Specifically, 180 D-loop haplotypes showed a complex reticulated network (Fig. 3), while 60 ND5 haplotypes revealed a more clarified network pattern (Fig. 4). No dominant haplotype was detected based on the D-loop marker, however, Hap004 and Hap005 were considered the most abundant and common haplotypes, followed by Hap036 and Hap116. Among ND5 haplotypes, Hap01 was the most dominant haplotype followed by Hap03, Hap15, Hap06, and Hap10. Hap01 was found at all sampling sites and was considered the ancestral haplotype. A network with an ancestral haplotype typically shows a star-like or star-burst appearance with the ancestral haplotype centered in it (Ferreri, Qu & Han, 2011).

Pairwise comparison Ф_ST analysis corroborated the low and non-significant population structure of T. tonggol from Malaysian waters (D-loop: −0.0324 to 0.1191 (Table 3); ND5: −0.0254 to 0.0739 (Table 4)), except for seven and a single significant pairwise comparison involving BT, based on the D-loop and ND5 sequences, respectively. Further genetic differentiation assessment based on H_ST (D-loop: 0.0033; ND5: 0.0034), N_ST (D-loop: 0.0113; ND5: 0.0062) and K_ST* (D-loop: 0.0064; ND5: 0.0049) produced low and not significant values, which corresponds with a high level of gene flow (Nm) among Malaysian T. tonggol populations (D-loop: 145.69 and ND5: 146.01 for haplotype-based statistic; D-loop: 43.72 and ND5: 102.46 for sequence-based statistic). Correspondingly, the pairwise genetic distances among populations also exhibit relatively low values ranging from 0.0197 to 0.0251 (D-loop) and 0.0020 to 0.0038 (ND5). The hierarchical AMOVA indicated that more than 99% of the total genetic variation in Malaysian T. tonggol was contributed by genetic differences within populations. Attempts to identify if population subdivisions exist among the T. tonggol populations (BT vs. other populations) returned a non-significant F_CT value with less than 1% contribution to the total genetic variation, while more than 99% of the total genetic variation was contributed within populations, based on both datasets. The Mantel Test also supported earlier findings, demonstrating no correlation between genetic differentiation (pairwise Ф_ST value) and geographical distance (D-loop: r = −0.3763, P = 0.08 and ND5: r = −0.0050, P = 0.43) among Malaysian populations.

D-loop sequences of T. tonggol from Taiwan (TW), the east coast of the Indian Ocean (ECIO), and some of the west coast of the Indian Ocean (WCIO) were clustered with the Malaysian haplotypes, while some other WCIO haplotypes were placed into another clade with high bootstrap support (Fig. 2A; Supplemental 1). ML tree partitioning was partly in agreement with the pairwise comparisons Ф_ST, where TW was not significantly structured for Malaysia, ECIO nor WCIO. In contrast, all pairwise comparisons involving ECIO and WCIO against Malaysian’s populations were statistically significant after FDR correction at α = 0.05, except for WCIO-KT (Table 3). Meanwhile, WCIO and ECIO were not genetically subdivided from each other (P > 0.05) (Table 3). A hierarchical AMOVA revealed the existence of genetic subdivision between the Indian Ocean and Malaysian waters (F_CT: 0.09, P < 0.05), yet 90.01% of the genetic variation within the Indo-Pacific region was contributed by genetic differences within populations. Pairwise genetic distances between TW and Malaysian populations ranged from 0.0181 to 0.0219, while the genetic distances were relatively higher for pairwise comparisons involving ECIO and WCIO, that is, from 0.0238 to 0.0371 (Table 3).

Demographic history

Negative values of Tajima’s D and Fu’s F_S (all significant at P < 0.05) neutrality tests were detected in all populations inferred from both the mtDNA D-loop and ND5 gene (Table 2). Large differences in population sizes before (θ₀) and after expansion (θ₁) were detected, i.e (on average) 0.116 and 54,915.300 based on D-loop sequences, while 0.103 and 81,838.400 were based on the ND5 gene marker (Table 2). Corresponding to the τ value of 2.321 (ND5) and 8.514 (D-loop) (Table 2), the calculated expansion time for T. tonggol in Malaysian waters was 67,865 and 284,252 years ago inferred by ND5 and D-loop markers, respectively Bayesian skyline plot (BSP) analysis revealed two significant increases in effective population size that occurred 200,000 and 950,000 years ago based on the D-loop marker (Fig. 5A), while continuous expansion started 150,000 years ago with a more recent expansion around 100,000 years ago based on the ND5 gene marker (Fig. 5B). Goodness of fit tests (Hri and SSD) exhibited non-significant values for the overall samples (P > 0.05) (Table 2). Population demographic analysis of T. tonggol matched a unimodal distribution for overall samples (Fig. 6).

Genetic diversity

According to Grant & Bowen (1998), the past demographic history of populations can be clarified based on their contemporary haplotype diversity (H) and nucleotide diversity (π). These two sensitive indices are the basis of genetic diversity estimation of a population (Nei & Li, 1979). In this study, all populations of T. tonggol showed high haplotype diversity (D-loop: 0.990–1.000; ND5: 0.848–0.965) but low to moderate nucleotide diversity (D-loop: 0.0195–0.0250; ND5: 0.0017–0.0039) (Table 2). High haplotype diversity, coupled with low nucleotide diversity, indicates a large population size that has undergone recent population expansion, which allows the retention of new alleles in the population but without sufficient time for accumulation of more nucleotide substitutions among haplotypes (Grant & Bowen, 1998; Chen et al., 2015; Delrieu-Trottin et al., 2017; Nabilsyafiq et al., 2019). These results were in agreement with previous findings of several pelagic fishes, including the spotted mackerel, Scomber australasicus (H = 0.996, π = 0.007) (Tzeng, 2007), yellowfin tuna T. albacares (H = 0.997, π = 0.035), skipjack tuna Katsuwonus pelamis (H = 0.999, π = 0.084) (Ely et al., 2005), and longtail tuna, T. tonggol (H = 0.999, π = 0.0016) (Willette, Santos & Leadbitter, 2016). Furthermore, the wide difference in nucleotide variability estimates between the two markers was consistent with the findings of Viñas & Tudela (2009), where the D-loop region had a ten-fold higher value of nucleotide diversity compared to the mtDNA ND5 gene.

Phylogenetic and population level analyses

Populations of T. tonggol from Malaysian waters exhibited an absence of geographical structure associated with mtDNA sequences, as evidenced by the single-clade gene trees (Fig. 2; Supplemental 1), ambiguous genetic partitioning of haplotype networks (Figs. 3 and 4), low and non-significant values of pairwise Ф_ST (Tables 3 and 4) (except for several comparisons involving BT population), high contribution of within population variation through AMOVA, and non-significant correlation between genetic differentiation and geographical distance. These results strongly suggest that the T. tonggol populations in Malaysian waters were panmictic with shallow genetic structure due to high gene flow, similarly reported in other studies of the same species (Kunal et al., 2014; Willette, Santos & Leadbitter, 2016) and several other species (Carpenter et al., 2011). Wright (1931) suggested that the level of genetic differentiation among populations is related to the rate of evolutionary processes, like migration, mutation, and drift. Thus, a highly migratory species with a large population size, such as T. tonggol, is predicted to show limited population partitioning. BT population showed significant genetic structure from the rest (except TG and MR) and SB, inferred from the D-loop and ND5 sequences respectively, based on the population pairwise Ф_ST analysis, however, all other analyses suggested genetic homogeneity with other T. tonggol populations in Malaysian waters. We believe that this could be due the different weightage of algorithms or characters used in the various analyses, perhaps a different emphasis on nucleotide versus haplotype diversity.

Pelagic fish in the marine realm are well-known to exhibit little genetic divergence (Ely et al., 2005). The weak genetic structure observed in T. tonggol within the pelagic environment is typical of pelagic fish due to their biological and life histories (Fauvelot & Borsa, 2011; Pedrosa-Gerasmio, Agmata & Santos, 2014). T. tonggol spawns during the monsoon season (Koya et al., 2018) where the ocean circulation shift (upwelling and down welling) during this season would enrich the water at the surface and thus lead to the growth of plankton (Yohannan & Abdurrahiman, 1998). Although T. tonggol is not a plankton feeder, the plankton bloom somehow enriches the food resources for other fish that become the prey of T. tonggol, hence, creating an optimal spawning ground for the species. T. tonggol is believed to spawn close to coastal waters (Nishikawa & Ueyanagi, 1991), thus the dynamic movement of waters during monsoon, not only helps in circulation of rich nutrients, but also in the larvae dispersal that could span a larger area (Madhavi & Lakshmi, 2012). Another possible explanation for the absence of limited gene flow among populations is the pattern of migration in the adult stage. In addition to the high dispersal potential during egg and larval stages, adults are characterized by high maneuverability during seasonal migration.

An earlier study by Willette, Santos & Leadbitter (2016) showed delineation of the Indian Ocean versus the South China Sea (samples were collected from Vietnam, Indonesia, and the Philippines) but without representatives from Malaysian waters. In the present study, based on the increased sample populations on a finer scale within this biogeographical region and GenBank sequences, we hypothesize that the genetic barrier lies within the Andaman Sea, which hindered partial gene flow between the coast of India and Malaysian waters, as evidenced in the ML tree (Fig. 2A; Supplemental 1), pairwise comparisons Ф_ST (Table 3), and hierarchical AMOVA. There was also another possible break between ECIO and WCIO based on results in the ML tree (Fig. 2A; Supplemental 1) and pairwise comparisons Ф_ST, though the moderate pairwise Ф_ST value is not significant (Table 3). However, the phylogenetic relationship of T. tonggol from Malaysian waters and other regions of the South China Sea remains unknown due to the limited genetic data and unavailability of the haplotype sequences from the study by Willette, Santos & Leadbitter (2016) in the public database. We postulated that close genetic relationships would be expected, based on the recent findings regarding the absence of genetic subdivision between TW and Malaysian T. tonggol. Future studies should include more detailed sampling within the Andaman Sea and adjacent waters to substantiate this.

Demographic history

Historical events during the Pleistocene epoch could have shaped the genetic diversification of T. tonggol populations observed in the present study. All relevant statistical tests implied a scenario of past population/demographic expansion in the absence of background selection. Negative Fu’s F_S values signified the alterations caused by population expansion and/or selection (Fu, 1997), which was further supported by Tajima’s D that implied a notable population growth or genetic hitchhiking in a background of recent excess mutations (Tajima, 1989). Likewise, the non-significant sum of squared deviations (SSD) and Harpending’s raggedness index (Hri) indicated the occurrence of population expansion in T. tonggol (Kunal et al., 2014) that inhabits Malaysian waters. Furthermore, the star-like pattern of the median–joining network (Fig. 4) and unimodal pattern of mismatch distribution (Fig. 6) further support the occurrence of a sudden demographic expansion during recent history of the taxa (Slatkin & Hudson, 1991; Rogers & Harpending, 1992; Ferreri, Qu & Han, 2011; Kunal et al., 2014; Chen et al., 2015; Pedrosa-Gerasmio, Agmata & Santos, 2014).

The large population size differences before (θ₀) and after expansion (θ₁) also suggested a rapid population expansion of T. tonggol in the past as also reported in a previous study of populations from India (Kunal et al., 2014). The overall τ value observed in Malaysia was much lower than T. tonggol populations from Indian waters, which was 21.26. The estimated time for population expansion in Indian waters was 593,334 years before the present (Kunal et al., 2014), as compared to 67,865 and 284,252 years ago inferred by the ND5 and D-loop markers, respectively, for T. tonggol from Malaysian waters. This suggests that T. tonggol in the Indian region underwent earlier expansion, with subsequent large population retention. In this study, the ND5 gene marker was able to detect a more recent population expansion of T. tonggol populations in Malaysian waters around 70,000 years ago (based on tau value) and 100,000 years ago (based on BSP analysis). Based on the D-loop marker, two expansion events were detected, where the first round occurred around 200,000 years ago (based on BSP analysis) or 284,000 years ago (based on tau value) (during middle Pleistocene (126,000–781,000 years before present) (Saul, 2016)) and the second round occurred around 950,000 years ago (based on BSP analysis) (during early Pleistocene) (Fig. 5). In general, both markers were able to detect population expansion that occurred around 200,000 years ago (ND5 gene marker indicated population expansion started before 150,000 years ago).

Implication for fisheries management

To date, there is limited information on the population structure of T. tonggol, especially in Malaysia, and from the management point of view, this is a critical issue. The present study provides the first baseline population genetic data on T. tonggol populations in Malaysian waters, which is important information for management planning by authorities.

Managing fishery resources takes a significant effort to protect and replenish the genetic pools for a sustainable harvest. Molecular analyses, in complement with other approaches, may serve as a reliable measurement for an efficient preservation strategy (BjØrnstad & Ried, 2002; Toro, Barragan & Ovilo, 2003). The genetic data suggest that T. tonggol in Malaysia forms a panmictic population as observed by the wide distributional range of this species and non-significant low Ф_ST values among the populations studied, thus suggesting a single evolutionary significant unit (ESU). However, this study was based on mitochondrial markers and therefore, restricted to only the pattern of maternal inheritance. For a holistic genetic perspective of bi-parental inheritance, co-dominant markers, such as microsatellites, should be included.