Figures
Abstract
Phylogeographic patterns and population structure of the pelagic Indian mackerel, Rastrelliger kanagurta were examined in 23 populations collected from the Indonesian-Malaysian Archipelago (IMA) and the West Indian Ocean (WIO). Despite the vast expanse of the IMA and neighbouring seas, no evidence for geographical structure was evident. An indication that R. kanagurta populations across this region are essentially panmictic. This study also revealed that historical isolation was insufficient for R. kanagurta to attain migration drift equilibrium. Two distinct subpopulations were detected between the WIO and the IMA (and adjacent populations); interpopulation genetic variation was high. A plausible explanation for the genetic differentiation observed between the IMA and WIO regions suggest historical isolation as a result of fluctuations in sea levels during the late Pleistocene. This occurrence resulted in the evolution of a phylogeographic break for this species to the north of the Andaman Sea.
Citation: Akib NAM, Tam BM, Phumee P, Abidin MZ, Tamadoni S, Mather PB, et al. (2015) High Connectivity in Rastrelliger kanagurta: Influence of Historical Signatures and Migratory Behaviour Inferred from mtDNA Cytochrome b. PLoS ONE 10(3): e0119749. https://doi.org/10.1371/journal.pone.0119749
Academic Editor: Peng Xu, Chinese Academy of Fishery Sciences, CHINA
Received: June 4, 2014; Accepted: February 3, 2015; Published: March 18, 2015
Copyright: © 2015 Akib et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: Universiti Sains Malaysia funded this project under the Research University Grant (1001/PBIOLOGI/815051) and the Postgraduate Research Grant Scheme (1001/PPANTAI/844103).
Competing interests: The authors have declared that no competing interests exist.
Introduction
In the marine realm, fish are abundant and ubiquitous and observed levels of genetic differentiation among populations are often low [1]. This pattern results from combined effects of large spawning population sizes that limit genetic drift effects, the apparent absence of physical barriers to dispersal for many species in their environment and presence of highly dispersive life history stages that contribute to high rates of gene flow among populations [2]. With no obvious major physical barriers to gene flow in the marine environment, genetic homogeneity or absence of spatial patterns in allele or haplotype distributions among geographical populations of marine species is expected and is often confirmed in genetics studies [3–4].
However, despite an apparent lack of any physical barriers to gene flow, many marine taxa have often been reported to show population structure [5–11]. Where structured populations are detected, this has been attributed to a number of factors including; historical and contemporary interplay among a complex set of ecological, demographic, behavioural, genetic, oceanographic, climatic and/or tectonic processes [12]. Evidence for population structure in several marine species has been reported in this region. A study of Plectorhinchus flavomaculatus, a coral reef fish sampled from the South China Sea in the Xisha, Zhongsha and Nansha archipelagos examined variation in the mtDNA control region and identified two discrete lineages [11]. A population pattern of ‘isolation by distance’ was observed in this species. The authors hypothesized that populations may have diverged in different glacial refuges during the Pleistocene at times of lower sea levels. Genetic relationships among northern populations of the six bar wrasse, Thallasoma hardwicki were examined from six localities in the South China Sea and three localities from the Solomon Islands in the South Pacific Ocean [13]. Three major groups were detected; a north-central group comprised of northwestern Taiwan and northern Vietnam; a southwestern group in southern Vietnam; and a southern group that included the central Philippines. Differentiations among these regions were attributed to limitations on gene flow that result from the impacts of sea surface currents. In addition, a study of the crimson snapper, Lutjanus erythropterus [8] reported that populations in East Asia were divided into two major clades: an eastern group that included populations in the western Pacific Ocean and the East Sea, and a South China Sea group that included populations from northern Malaysia extending to the South China Sea. The study suggested that limited gene flow between the eastern region and the South China Sea contributed to the observed geographical subdivision. Significant population structure was also evident in the white pomfret, Pampus argenteus, with three distinct clades (the South China Sea, the Arabian Sea and the Bay of Bengal) across the Indo-West Pacific region [14]. Late Pleistocene population expansions were hypothesized to have contributed to the reported population structure in this species.
The Indian mackerel, Rastrelliger kanagurta is an epipelagic Scombrid that occurs widely across the Indo-West Pacific from South Africa, the Seychelles and the Red Sea, east to Indonesia and northern Australia to Melanesia, Micronesia, Samoa, China and the Ryukyu Islands of southern Japan. It is also believed to have entered the eastern Mediterranean Sea via the Suez Canal [15]. It is one of three species in the genus Rastrelliger that also includes R. brachysoma and R. faughni. In the seas surrounding Malaysia, R. kanagurta is widespread across the northwest and east peninsular Malaysia, Sabah and Sarawak [16] where it is one of the most commercially important marine resources [17]. This species is a fast swimmer because of a fusiform and streamlined body and therefore highly migratory [18].
To date, there has only been limited study of the population genetics of R. kanagurta. These have included investigations of west peninsular populations employing RAPD markers [19] and mtDNA D-loop region [20]. A similar study, also based on RAPD markers examined three populations of R. kanagurta from the Indian Peninsular [21] while in another study allozyme markers (five enzymes and a single sarcoplasmic protein) were employed to investigate population structuring along the coastal waters of India and the Andaman Sea [22]. All studies [20–22] reported essentially panmixia among sampled populations, although weak population structure was detected between northern and southern populations of R. kanagurta in the Strait of Malacca off the coast of Perak [19].
Genetic markers including; allozymes, mtDNA and microsatellites have been employed widely as markers to investigate the population structures of many marine fish taxa [23–25]. Different types of markers however, can sometimes produce patterns of population structure that are not concordant [23]. In the current study, mtDNA cytochrome b (Cyt b) was utilised to assess population structure in R. kanagurta. This marker was chosen because Cyt b contains both slowly evolving codon positions, which are useful in phylogenetic studies and rapidly evolving codon positions that are useful for population studies [4, 26–27]. These characteristics make mtDNA Cyt b a suitable marker for studies of both phylogenetic and population genetic studies of fish [27–28].
The current study focused on examining genetic variation in R. kanagurta populations in waters that form part of a major marine biodiversity hotspot in the Indonesian—Malay Archipelago (IMA). These include the Strait of Malacca (an extension of the Indian Ocean), the South China Sea, the Sulu Sea and the Celebes Sea. In addition, four populations outlying this area were included as reference populations and to provide a wider coverage of the species natural distribution. These were from Can Tho, Vietnam, represented the northern South China Sea and two populations from Trang, Thailand and Banda Acheh, Indonesia represented the East Indian Ocean (Andaman Sea). The fourth population from Bandar Abbas, Iran represented the West Indian Ocean and served as the outgroup. Variation in the complete Cyt b gene was utilized to examine genetic structure in R. kanagurta across this vast expanse of waters.
Materials and Methods
Sampling and DNA sequencing
A total of 296 R. kanagurta individuals were collected from 19 locations within the IMA; Strait of Malacca (SOM), South China Sea (SCS), Sulu Sea and Celebes Sea between 2010 to 2011 (Table 1, Fig. 1A). In addition, three populations outlying this region, namely Banda Acheh, and Trang (both in the Andaman Sea) and Can Tho (South China Sea but outside of the IMA) were included. A West Indian Ocean (WIO), Bandar Abbas was also included. These additional populations to the IMA provided reference populations and a better geographical coverage, making in total, 342 individuals sampled from 23 locations (Fig. 1B). Populations were subdivided into seven geographical regions according to the seas for the analysis (refer Table 1). The South China Sea (SCS1 and SCS 2) was further subdivided due to its extensive wide geographical area. For the Malaysian samples, the samples were either obtained from fishing vessels operating in Zone A (from the shoreline to 5 nautical miles) or Zone B (from 5 to 12 nautical miles) at predetermined landing sites where catches do not overlap due to regulatory enforcement. Therefore, overlapping of populations from various landing sites could be disregarded. We get the fish samples throughout Malaysia with the help from Department of Fisheries through out Malaysia—Perlis, Kedah, Pulau Pinang, Perak, Selangor, Kelantan, Terengganu, Pahang, Johor, Sabah and Sarawak. For samples outside Malaysia, fin clips preserved in 99% ethanol were obtained through the assistance of local collaborators with permission from the respective authorities. All fishermen have the appropriate permits to fish. The field studies did not involve endangered or protected species. This study does not require an ethics statement because these samples are sea food product that can be found in the market.
In all sampling sites where they were physically collected, none of the samples were alive upon collection and they were not sacrificed for the purposes of this study. Fin clips were removed from the right pectoral fin from the dead samples and then preserved in 99% ethanol.
Extraction of DNA template was done using the Aqua Genomic DNA isolation kit (MultiTarget Pharmaceuticals, Salt Lake City, Utah 84116) as according to the manufacturer’s protocol. Isolated DNA template was PCR amplified for the complete mitochondrial Cyt b. Primer pair used for amplification was; L14317 (5’-CAG GAT TTT AAC CAG GAC TAA TGG CTT GAA-3’) and H15630 (5’-TTA ATT TAG AAT CCT AGC TTT GG-3’) (Takashima, unpublished). The PCR reaction mixture consisted of 2.5μL 10x PCR buffer (iNtRON), 3.3μL 25mM MgCl2 (iNtRON), 0.5μL 2.5mM of dNTP (iNtRON), 10 pM of each primer, 5U of Taq polymerase (iNtRON) and ∼50ng of genomic DNA. Amplification conditions were; initial denaturation at 94°C (2 min); 35 cycles of 94°C (20 sec), 50°C (20 sec), 72°C (1 min 10 sec) and a final extension at 72°C (5 min) before termination of the reaction at 4°C. Amplification of the 25 μL reaction volume was conducted in a G-Storm Thermal Cycler (MJ Research, Waltham, MA, USA). PCR products were purified using Promega Wizard SV Gel and PCR Clean-Up System (Promega Corporation, Madison, WI, USA) based on the manufacturer’s specifications. All purified products were sent for DNA sequencing (First BASE Laboratories Sdn. Bhd., Selangor, Malaysia and Centre for Chemical Biology, USM) in both the forward and reverse directions using Applied Biosystem ABI3730x1 capillary-based DNA sequencer.
Data Analysis
Forward and reverse Cyt b sequences were edited in MEGA 5 [29]. Sequences were compiled and aligned using ClustalW that is integrated within MEGA 5 for unambiguous operational taxa units. Haplotype distributions for sampled populations were assessed using Collapse [30] and summarized using DnaSP 5.10 [31]. Phylogenetic relationships among haplotypes were constructed using Neighbour-Joining (NJ) method in MEGA 5 with a confidence level assessed using 1000 bootstrap replications. To view haplotype relationships, a phylogenetic network of all haplotypes was constructed based on median joining calculation in Minimum Spanning Network (MSN) [32]. The complete aligned dataset was then screened for nucleotide variable sites, parsimony informative sites, number of haplotypes, number of transitions and transversions and nucleotide frequencies in MEGA 5. The same software was used to estimate, genetic diversity within and among populations based on a K2P model. Genetic diversity, including gene diversity, nucleotide diversity and Theta S (θs) were retrieved in Arlequin 3.1 [33]. Historical demographic expansion was tested using Fu’s Fs [34] statistic in Arlequin 3.11 [33] in addition to Ramos-Onsins and Rosas’ R2 [35] in DnaSP 5.10 [31]. Significant values for R2 were calculated using coalescent simulations with 5000 replicate runs for each simulation. To detect whether sampled populations were demographically stable or expanding or decreasing over time, the demographic history of each sampled population was assessed applying Harpending’s raggedness index, Hri [36] implemented in Arlequin 3.1 [33] and mismatch distributions [37–39] implemented in DnaSP 5.10 [31]. Both Harpending’s raggedness index and mismatch distributions were used to detect whether the sequence data from each population deviated from what would be expected under a sudden population expansion model. A significant Hri value (P < 0.05) is evidence for rejecting this model [39]. Significance of the above tests was assessed using 10000 permutations. Past demographic events were also investigated based on the distributions of pairwise differences between sequences in Arlequin. Three parameters were estimated, θ0 and θ1(i.e. θ before and after the population growth) and τ (relative measure of time expressed in units of mutational time since the population expansion) [38, 40]. Values of τ were transformed to estimate real time expansion using the equation τ = 2μt, (μ = mutation rate of the sequence analysed, t = time since expansion). Changes in effective population size (Ne) across time were inferred using Bayesian skyline analyses. These analyses enable past demographic changes to be inferred from the current patterns of genetic diversity within a population [41]. BEAST v1.8 [42] was used to create the Bayesian skyline plots. To test for convergence, analysis was run for 108 iterations with a burn-in of 107 under the TN93 model, a strict molecular clock and a stepwise skyline model. Genealogies and model parameters were sampled every 1000 iterations. All operators were optimized automatically. Result of skyline plots were generated by Tracer1.5 [43].
Further analyses of genetic differentiation among populations were conducted for haplotype-based statistics (HST) and sequence-based statistics (NST) [44] with significance levels assessed using permutation tests with 1000 replicates [45] in DnaSP 5.10 [31]. Estimates of gene flow (Nm) based on both haplotype and sequence statistics were derived employing the same program [45–46]. Genetic distance estimates between sampled populations were calculated using a Kimura 2 parameter distance method in MEGA 5. Analysis of molecular variance (AMOVA) was performed to estimate molecular variance among sampled populations using Arlequin 3.11 [33]. AMOVA partitions the total genetic variance into three measures of haplotypic diversity; FST describes variation between populations within total, FSC describes variation among populations within region and FCT describes variation among regions within total [33]. Spatial structure was examined further using Spatial Analysis of Molecular Variance (SAMOVA) v.1.0 [47] to identify groups of populations that were geographically homogeneous and maximally differentiated from each other and to identify genetic barriers between these groups [47].
Results
All 342 individuals sampled from 23 localities amplified successfully for the Cyt b gene and a final sequence length of 1140 bp was obtained after alignment and editing of ambiguous sequences revealing 241 unique haplotypes. Sequences were deposited in Genbank with accession numbers JQ681542-JQ681735. A total of 220 variable sites were identified at the first and third base positions in the codon with a ratio of 1:8. Of these, 122 (10.7%) were parsimony informative. The ratio of transversion to transition substitutions of the entire data set was 1:11. Haplotype diversity was very high, ranging from 0.952 to 1.000, while nucleotide diversity was very low, ranging from 0.008 to 0.011 (Table 2).
Due to the higher number of unique haplotypes detected, only haplotypes shared among sites are shown for clarity (S1 Table). These unique haplotypes consist of mainly singletons (213/241–88.38%) and they are population specific (221/241–90.46%). Out of a total of 241 unique haplotypes, only 20 (9.54%) were found in more than a single region, with the rest being regionally specific. The BA (WIO) population from the Indian Ocean consisted mainly of singletons (16 of 17), all of which were private haplotypes to BA. As observed, no obvious population structure was detected, either related to the sea or with geographical distance. Several haplotypes were shared only among the IMA populations, but this was probably because of the low sample number of sites outside of the IMA waters as many other haplotypes were common among different seas (e.g. haplotypes 12, 52 etc.).
The NJ tree was moderately supported with four lineages identified (Fig. 2). While no obvious geographical structure was evident for Clade 1 and 2, in contrast, Clades 3 and 4 were mostly composed of the Iranian population (western Indian Ocean) with only two representatives from the Strait of Malacca (haplotypes 14, 45 respectively). The minimum spanning network (MSN) produced a complex reticulation of 241 haplotypes (Fig. 3) with four major haplotypes. Two of the four highest frequency haplotypes were mainly formed of the SCS-2 sites (haplotypes 8 and 29). The other two were mostly represented by the SCS-1 populations (haplotype 21) and the SLS population (haplotype 1). According to the coalescent theory [48], common haplotypes internal to the network can be inferred to be ancestral, while tip haplotypes are derived or descendant from ancestral (internal) haplotypes. The Bandar Abbas (WIO) haplotypes are mostly at the tip. The occurrence of star-like patterns radiating from these major haplotypes suggest that R. kanagurta populations have undergone significant population size expansions in the relatively recent past [37,49].
Coloured close circles represent different regions (refer to legend). mv = median vectors. The red circle in the figure represents the BA tip haplotypes.
Overall, 14 out of 23 population tested were significant with negative values for Fu’s Fs. This suggests historical population expansions at these sites [34] (Table 3) associated with an excess of recent mutations or rare alleles [34]. The Harpending raggedness index (Hri = 0.031, p>0.05) corroborated results from Fu’s Fs analysis indicating that some populations had expanded recently.
Past population demographic of R. kanagurta matched the bimodal mismatch distribution for all four clades. Comparison between observed frequencies of pairwise differences with those expected under various demographic models [50] (Fig. 4) added support for a spatial expansion followed by a recent demographic expansion [50–51]. Bayesian skyline plot revealed an episode of abrupt demographic expansion around 100 to 120 thousand years (kyr) before present, which also falls into the Pleistocene period (Fig. 5). The flat skyline plot also indicates that the population has remained at a constant size over time before the expansion. Time since the population expansion was estimated to be 7.816/2u generations. Given a mutation rate of perciform Cyt b of 1–2% per million years [52–53], the R. kanagurta population expansion was estimated to have taken place 171,403 to 342,807 years bp.
Black line represents median estimate Ne; light lines are the upper and lower 95% highest posterior density (HPD) limits of Ne.
AMOVA analysis performed on the 23 populations, subdivided into seven geographical groups according to the seas (group 1: SOM; group 2: SCS-1; group 3:SCS-2; group 4:SS; group 5: CS; group 6: AS; group 7: WIO). The South China Sea (SCS1 and SCS 2) was further subdivided due to its extensive wide geographical area. AMOVA analysis showed that among group variation was less than 1% (FCT = 0.0023, p<0.05) but significant genetic differentiation was present between at least two groups. Genetic differentiation among populations within groups was 0.41% and significant (FSC = 0.0041, p>0.05). Presumably this results from the divergence of the BA population from all other sampled sites when site BA was included in the analysis. Genetic variation within populations was high at 99.36% revealing that most of the variation was present within populations, while FST were low and not significant (FST = 0.0064, p<0.05). The analysis was repeated with the Iranian (BA) population (WIO) excluded to test. The test is to discover whether genetic differentiation was evident at finer spatial scales across the IMA and the surrounding seas (Trang, Banda Aceh—Andaman Sea of east Indian Ocean; Can Tho- South China Sea). The samples were divided into 6 groups; group 1: SOM; group 2: SCS-1; group 3: SCS-2; group 4: SS; group 5: CS, group 6:AS. The FCT value for this additional analysis was not significant (p>0.05).
Results of a Mantel test [54] supported previous outcomes showing no significant correlation between genetic differentiation (pairwise FST value) and geographical distance (r = 0.370, p>0.05, 1000 permutations) among the sampled populations (Fig. 6). The graph, however, did show two distinct clusters. The first cluster comprised all samples from the Strait of Malacca, South China Sea, Sulu Sea and Celebes Sea and the Andaman Sea (eastern Indian Ocean- Banda Aceh and Trang), while the second contained the western Indian Ocean population of Bandar Abbas. The Cyt b FST values were close to 0 (-0.030–0.042) (data not shown) and were not significant (p>0.05) after FDR correction (α = 0.05). An indication of minimal or no genetic structure present. Further analyses of genetic differentiations estimate based on HST (0.0024), NST (0.1581) and KST* (0.0246) were very low and all were not significant, which also inferred a high level of gene flow among sampled R. kanagurta populations (Nm = 136.45 for haplotype-based statistic and Nm = 7.39 for sequence-based statistic).
The red circle indicates the western Indian Ocean population (Bandar Abbas, Iran).
Based on estimates of k = 2 to 8 in the SAMOVA analysis each group represented a statistical discrete phylogeographic grouping. However, k = 2 yielded the highest FCT value (FCT = 3.0374, p<0.05), a clear signal that there were two genetically different R. kanagurta stocks among the sampled populations, namely the IMA and adjacent regions and a western Indian Ocean population at Bandar Abbas.
Discussion
In this study, high haplotype diversity, Hd values (ranging from 0.952 to 1.000) coupled with very low nucleotide diversity, π (ranging from 0.008 to 0.011) observed in most populations concord with reports from earlier studies of several pelagic marine fishes. Previous studies include the spotted mackerel, Scomber australasicus (Hd = 0.996, π = 0.007) [55]; wahoo, Acanthocybium solandri (Hd = 0.918, π = 0.006) [56]; yellow fin tuna, Thunnus albacares (Hd = 0.997, π = 0.035) and skipjack tuna, Katsuwonus pelamis (Hd = 0.999, π = 0.084) [57]. Population histories of species can be interpreted into four categories based on combinations of their Hd and π values from small to large [58]. In the case of R. kanagurta, high Hd and low π values suggest a recent population expansion from a low effective population size, i.e., rapid population growth after a bottleneck event with associated accumulation of many new mutations [59].
Four phylogenetic groups were evident in the Cyt b data for R. kanagurta. There were, however, no distinct geographical lineages present among the IMA and neighbouring marine populations in this study. All of the IMA sites are grouped in Clades 1 and 2 with no obvious pattern of geographic structuring. Clade 3 and 4 were almost exclusively composed of individuals from the western Indian Ocean at Bandar Abbas. Several Iranian individuals were also members of Clade 1 (hap 229, 238) and Clade 2 (hap 228). Two haplotypes present in the Strait of Malacca were members of Clade 3 (hap 45) and Clade 4 (hap 14), respectively. These indicate only very limited genetic exchange between the two regions. The divergence of these clades, however, may have happened in the past during the lowering of sea level during the Pleistocene period. Past geological and climatic events have undoubtedly played a major role in this inferred population expansion. The paleogeography of the IMA changed dramatically across the Quaternary period [60]. This area encompasses mainland Southeast Asia and Australasia, and is also known as the Indo-Australian Archipelago or East Indies. During the Pleistocene glacial periods, lowering of ocean levels periodically exposed the Sunda Shelf and Sahul Shelf within this archipelago. The former is an extension of the coastal shelf of Southeast Asia, including the Malay-Peninsula, Sumatra, Java and Borneo to Palawan. Fluctuations in sea levels during the Pleistocene epoch would have resulted in a significant impact on the dispersal of marine and terrestrial species. Thus shaping the modern genetic structures of many species that inhabit this region as we see it today [61–64]. This scenario had been hypothesized for many marine species including the redlip mullet, Chelon haematocheilus [65] and Japanese sea bass, Lateolabrax japonicas [66].
During the last Interglacial Period, there were marine connections between the Indian Ocean and the South China Sea via the Straits of Singapore [64]. These connections would have permitted extensive population expansion by R. kanagurta where previously they would have been isolated. Estimation of the time since expansion of R. kanagurta suggests that the timing of this event ranged from 171,403 to 342,807 years bp. The timing of the event is consistent with the occurrence of sea level rises cyclically during the late Pleistocene between 1,600,000 to 10,000 years ago [8,67]. Land bridges that were exposed during periodic declines in sea level prior to the last Interglacial period were likely to have hampered dispersal of R. kanagurta and led to smaller isolated populations. During that time, no barriers to dispersal by sea were present, at least between Peninsular Malaysia and the islands of the Riau Archipelago and Sumatra. Thus consequently permitting free migration of R. kanagurta populations across the region. The last glacial maximum (LGM), 30,000 to 19,000 years ago led to the final decline in sea levels to modern times [67–68]. A substantial area of the Sunda and Sahul shelves were exposed periodically as sea levels associated with glaciations declined up to 200m below their present level [69]. This process would also have hindered dispersal by many marine taxa. At this time, the adjacent South China Sea was significantly reduced in size and formed a semi-closed marginal sea and this exposed a large low gradient on the Sundaland craton[67]. During this time, the modern Malayan Peninsula, Borneo and Sumatra formed highlands to the South [70] producing a potential barrier to population migration.
Isolated marine populations at this time would have been subjected to increased levels of inbreeding. Potentially leading to declines in population levels of genetic variation and genetic structuring of populations during this short period of the LGM. There are also possibilities of individuals that may have dispersed to the east or west of the adjacent flooded basins. This pattern is not concordant however, with the pattern observed in R. kanagurta. They experienced almost complete panmixia evident in modern populations in SE Asia over a vast expanse of two oceans. It can be conjectured therefore, that the brief duration of LGM, followed by a rise of sea level approximately 18000 years ago until the present modern sea levels were reached, was insufficient for R. kanagurta to attain migration drift equilibrium. This would be manifested in low genetic differentiation and low genetic structuring among extant populations until the present time.
Similar scenarios have been hypothesized for many marine species present across this region. Examples include the Indo-Pacific butterfly fishes [71], pelagic scads, Decapterus macrosoma and D. macarellus [72] and Eastern little tuna, Euthynnus affinis [73]. Based on the results of the MSN analysis (Fig. 3) it is postulated that this expansion probably occurred from east to west as evidenced by the Sundaland origin of the four major haplotypes (ancestral) while all unique WIO haplotypes were tip haplotypes (descendant). Detailed examination of the MSN analysis revealed that two of the four putative ancestral haplotypes in the SCS-1 and SCS-2, based on the highest frequency suggest that the centre of origin for modern R. kanagurta populations were in the South China Sea from which populations could have radiated both eastwards and westwards.
Marine fishes with high dispersal potential typically display low levels of population genetic structure and associated high levels of ongoing gene flow [74]. This would produce shallow evolutionary trajectories and, potentially, limit or reduce adaptive divergence among local populations [7]. The overall low and non-significant FST values reported here among the 23 R. kanagurta populations are concordant with a pattern of Cyt b homogeneity and high gene flow (Nm = 136.45). The level of genetic differentiation among populations has a predictable relationship with rates of important evolutionary processes (migration, mutation, drift) [75]. Thus, a highly migratory species with large populations should show limited population differentiation. This is in contrast to species with small populations and reduced migration rates, where we might expect greater differentiation [75].
Several other investigations have also reported relative panmixia in marine species across this region. For example, an apparent lack of genetic structure in marine species present in the Coral Triangle region of Southeast Asia has also been reported in several marine taxa including; pelagic scads, Decapterus macrosoma and D. macarellus [72] and the eastern little tuna, Euthynnus affinis [73]. Slight population differentiation was observed between the IMA (and adjacent populations) with the WIO Bandar Abbas population. Moderate to high rates of gene flow among populations can prevent sub-population isolation, hence maintaining genetic variation levels and reducing inbreeding [76–77]. This, however, not always a simple relationship. The hierarchical AMOVA results here showed low but significant FST and FCT outcomes. An indication of shallow genetic structure among populations of R. kanagurta. This was confirmed by spatial genetic heterogeneity (k = 2 in SAMOVA) between groups comprising populations in the IMA and surrounding seas vs a second group in the western Indian Ocean (Bandar Abbas, Iran). The Mantel test results did not support an ‘isolation by distance’ model. Thus a more plausible explanation for population differentiation between the two distinct groups of R. kanagurta was a result compounded by co-effects of various factors, including historical physical isolation during the Pleistocene epoch, larval dispersal factors and ocean currents.
Historical isolation would suggest the genetic differentiation evident between the IMA and WIO regions may have resulted from fluctuations in sea levels during the late Pleistocene. This led to a phylogeographic break for this species to the north of the Andaman Sea. Studies of other marine taxa across this region have reported co-effects on population differentiation of various factors. This was observed in the orange-spotted grouper, Epinephelus coiodes from China (South China Sea) and Malaysia through Indonesia (Southeast Asia) [78].
Ocean surface currents can also play an important role in marine larval dispersal. The surface current of the Indian Ocean is mainly under the influence of monsoon currents, also referred to as monsoon drift. It flows between the Bay of Bengal and the Arabian Sea [79]. In contrast, the seas surrounding Malaysia are affected by the Kuroshio current, the major surface current in the South China Sea [14]. These two sea surface currents could potentially contribute to population differentiation of eastern and western populations of R. kanagurta.
The Cyt b analysis of R. Kanagurta populations in the current study provide strong evidence for two discrete R. kanagurta stocks comprising western Indian Ocean and Southeast Asian populations. The two clades should therefore be considered different management units. While the patterns are clear, results need to be validated because they are based on evidence from a single maternal lineage gene that may not necessarily reflect the complete story.
Conclusions
The sampled R. kanagurta populations in Malaysian waters showed only weak spatial differentiation with high genetic variation as inferred from mtDNA Cyt b. The study did reveal however, two discrete populations of R. kanagurta, the Southeast Asian populations (South China Sea, Strait of Malacca, Sulu Sea, Celebes Sea, Andaman Sea) and an Iranian population (western Indian Ocean). In general, however, the data here suggest that R. kanagurta populations in the IMA and adjacent waters can be regarded as a single stock unit for management purposes. Their inferred demographic history also suggests that R. kanagurta populations potentially expanded in the late Pleistocene and that gene flow appears to be ongoing among extant populations today.
Supporting Information
S1 Table. Distribution of regionally shared haplotypes for 23 populations of R. kanagurta.
https://doi.org/10.1371/journal.pone.0119749.s001
(DOCX)
Acknowledgments
We would like to thank Universiti Sains Malaysia for providing doctoral fellowship under the Academic Staff Higher Education Scholarship (ASHES) to Noor Adelyna Mohammed Akib. We also would like to thank Universiti Sains Malaysia for funding this project under the Research University Grant (1001/PBIOLOGI/815051) and the Postgraduate Research Grant Scheme (1001/PPANTAI/844103). We are grateful to the Institute of Fisheries, Penang, Malaysia and Department of Fisheries throughout Malaysia for assisting us with the sampling. Thanks are also due to many colleagues and their respective institutions: Jamsari Amirul Firdaus Jamaluddin, Dr. Tan Min Pau, Mohd Lutfi Abdullah, Danial Hariz Zainal Abidin, Ahmad Faisal Ghani, Dr. Chee Su Yin, Lim Hong Chiun (USM) and Abdul Rahman Abdul Majid (Fisheries Research Institute, Penang, Malaysia).
Author Contributions
Conceived and designed the experiments: NAMA PBM SAMN. Performed the experiments: NAMA. Analyzed the data: NAMA PBM SAMN. Contributed reagents/materials/analysis tools: PBM SAMN. Wrote the paper: NAMA PBM SAMN. Contributed samples for this research: BMT PP MZA ST.
References
- 1. O'Reilly PT, Canino MF, Bailey KM, Bentzen P (2004) Inverse relationship between FST and microsatellite polymorphism in the marine fish, walleye pollock (Theragra chalcogramma): implications for resolving weak population structure. Mol Ecol 13(7): 1799–1814. pmid:15189204
- 2. Ward RD, Grewe PM (1994) Appraisal of molecular genetic techniques in fisheries. Rev Fish Biol Fisher 4: 300–325.
- 3. Scoles DR, Graves JE (1993) Genetic analysis of the population structure of yellowfin tuna, Thunnus albacores, from the Pacific ocean. Fish B-NOAA 91: 690–698.
- 4. Baker CS, Perry A, Bannister JL, Weinrich MT, Abernethy RB (1993) Abundant mitochondrial DNA variation and world-wide population structure in humpback whales. PNAS 90(17): 8239–8243. pmid:8367488
- 5. Roques S, Sévigny JM, Bernatchez L (2002) Genetic structure of deep-water redfish, Sebastes mentella, populations across the North Atlantic. Mar Biol 140(2): 297–307.
- 6.
Ward RD (2002) Genetics of fish populations. In: Hart PJB, Reynolds JD, editors. Handbook of Fish Biology and Fisheries Vol. 1. Blackwell Publishing, Oxford. pp. 200–224.
- 7. Nielsen EE, Nielsen PH, Meldrup D, Hansen MM (2004) Genetic population structure of turbot (Scophthalmus maximus L.) supports the presence of multiple hybrid zones for marine fishes in the transition zone between the Baltic Sea and the North Sea. Mol Ecol 3(3): 585–595.
- 8. Zhang J, Cai Z, Huang L (2006) Population genetic structure of crimson snapper Lutjanus erythropterus in East Asia, revealed by analysis of the mitochondrial control region. ICES J Mar Sci 63(4): 693–704.
- 9. O’Leary DB, Coughlan J, Dillane E, McCarthy TV, Cross TF (2007) Microsatellite variation in cod Gadus morhua throughout its geographic range. J Fish Biol 70: 310–335.
- 10. Reiss H, Hoarau G, Dickey‐Collas M, Wolff WJ (2009) Genetic population structure of marine fish: mismatch between biological and fisheries management units. Fish Fish 10(4): 361–395.
- 11. Han ZQ, Li YZ, Chen GB, Gao TX (2008) Population genetic structure of coral reef species Plectorhinchus flavomaculatus in South China Sea. Afr J Biotechnol 7(11): 1774–1781.
- 12.
Grosberg R, Cunningham CW (2001) Genetic structure in the sea. In: Bertness MD, Gaines S, Hay ME, editors. Marine Community Ecology. Sinauer, Sunderland, MA. pp. 61–84.
- 13. Chen CA, Ablan MCA, McManus JW, Bell JD, Tuan VS, Cabanban , et al. (2004). Population structure and genetic variability of six bar wrasse (Thallasoma hardwicki) in northern South China Sea revealed by mitochondrial control region sequences. Mar Biotechnol 6(4): 312–326. pmid:15129326
- 14. Sun P, Shi Z, Yin F, Peng S (2012) Population genetic structure and demographic history of Pampus argenteus in the Indo-West Pacific inferred from mitochondrial cytochrome b sequences. Biochem SysEcol43: 54–63.
- 15.
Collette BB, Nauen CE (1983) FAO species catalogue. v. 2: Scombrids of the world. An annotated and illustrated catalogue of tunas, mackerels, bonitos, and related species known to date. FAO Fisheries Synopsis125. New York: United Nations Food and Agriculture Organization. 137 p.
- 16.
Mansor MI (1989) Tumbesaran, kematian dan corak pengrekrutan Ikan Kembung Rastrelliger kanagurta (Cuvier) di Pantai Barat Semenanjung Malaysia. Fisheries Bulletin no. 59. Jabatan Perikanan: Kementerian Pertanian Malaysia. 22 p.
- 17.
Mansor MI, Syed Abdullah, Abdul Hamid Y (1996) Population structure of small pelagic fishes off the East Coast of Peninsular Malaysia. Fisheries Bulletin no. 99. Jabatan Perikanan: Kementerian Pertanian Malaysia. 30 p.
- 18. Collette BB, Russo JL (1984) Morphology, systematics, and biology of the Spanish mackerels (Scomberomorus, Scombridae). Fish Bull US 82: 545–692.
- 19. Darlina MN, Masazurah AR, Jayasankar P, Jamsari AF, Siti Azizah MN (2011) Morphometric and molecular analysis of mackerel (Rastrelliger spp) from the west coast of Peninsular Malaysia. GMR10: 2078–2092. pmid:21968625
- 20. Ahmad Faisal G, Danial Hariz ZA, Siti Azizah M.N Darlina MN (2012) Genetic variation of Indian Mackerel (Rastrelliger kanagurta) (Cuvier, 1816) of Sabah waters based on mitochondrial D-loop region: A preliminary study. AJBB 1: e100.
- 21. Jayasankar P, Thomas PC, Paulton MP, Mathew J (2004) Morphometric and genetic analysis of Indian mackerel (Rastrelliger kanagurta) from Peninsular India. Asian Fish Sci 17: 201–215.
- 22. Menezes MR, Nair S, Martins M (1993) Genetic divergence in the Indian Mackerel Rastrelliger kanagurta (Cuv) from the coastal waters of Peninsular India and the Andaman sea. Indian J Fish 40 (3): 135–141.
- 23. Ward RD (2000) Genetics in fisheries management. Hydrobiologia 420: 191–201
- 24. Pogson GH, Mesa KA, Boutilier RG (1995) Genetic population structure and gene flow in the Atlantic cod Gadus morhua: a comparison of allozyme and nuclear RFLP loci. Genetics 139(1): 375–385. pmid:7705638
- 25. O’Reilly PT, Hebinger C, Wright JM (1998) Analysis of parentage determination in Atlantic salmon (Salmo salar) using microsatellites. Anim Genet 29: 363–370.
- 26. Árnason Ú, Gullberg A (1996) Cytochrome b nucleotide sequences and the identification of five primary lineages of extant cetaceans. Mol Biol Evol 13(2): 407–417. pmid:8587505
- 27. Farias IP, Ortí G, Sampaio I, Schneider H, Meyer A (2001) The cytochrome b gene as a phylogenetic marker: the limits of resolution for analyzing relationships among cichlid fishes. J Mol Evol 53(2): 89–103. pmid:11479680
- 28. Aboim MA, Menezes GM, Schlitt T, Rogers AD (2005) Genetic structure and history of populations of the deep‐sea fish Helicolenus dactylopterus (Delaroche, 1809) inferred from mtDNA sequence analysis. Mol Ecol 14(5): 1343–1354. pmid:15813775
- 29. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S (2011) MEGA5: Molecular Evolutionary Genetics Analysis using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol Biol Evol 28: 2731–2739. pmid:21546353
- 30.
Posada D (2004) Collapse ver. 1.2. A tool for collapsing sequences to haplotypes. Available: http://darwin.uvigo.es/software/collapse.html. Accessed 1 June 2012.
- 31. Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25(11): 1451–1452. pmid:19346325
- 32. Bandelt HJ, Forster P, Rohl A (1999) Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol 16: 37–48. pmid:10331250
- 33. Excoffier L, Laval G, Schneider S (2005) Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evol Bioinform Online1:47–50.
- 34. Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147(2): 915–925. pmid:9335623
- 35. Ramos-Onsins SE, Rozas J (2002) Statistical properties of new neutrality tests against population growth. Mol Biol Evol 19(12): 2092–2100. pmid:12446801
- 36. Harpending HC (1994) Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol 66(4): 591–600. pmid:8088750
- 37. Slatkin M, Hudson RR (1991) Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics 129(2): 555–562. pmid:1743491
- 38. Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol 9(3): 552–569. pmid:1316531
- 39. Schneider S, Excoffier L (1999) Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics152(3): 1079–1089. pmid:10388826
- 40. Rogers AR (1995) Genetic evidence for a Pleistocene population explosion. Evolution 49(4): 608–615.
- 41. Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol 22: 1185–1192). pmid:15703244
- 42. Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol 29: 1969–1973 pmid:22367748
- 43.
Rambaut A, Drummond AJ (2007) Tracer v1.4. Available from http://beast.bio.ed.ac.uk/Tracer.
- 44. Lynch M, Crease TJ (1990) The analysis of population survey data on DNA sequence variation. Mol Biol Evol 7(3): 377–394.
- 45. Hudson RR, Boos DD, Kaplan NL (1992) A statistical test for detecting geographic subdivision. Mol Biol Evol 9(1): 138–151. pmid:1552836
- 46. Nei M (1973) Analysis of gene diversity in subdivided populations. PNAS 70: 3321–3323. pmid:4519626
- 47. Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the genetic structure of populations. Mol Ecol 11(12): 2571–2581. pmid:12453240
- 48. Posada D, Crandall KA (2001) Intraspecific phylogenetics: Trees grafting into networks. Trends Ecol Evol 16(1): 37–45. pmid:11146143
- 49. Forster P, Torroni A, Renfrew C, Röhl A (2001) Phylogenetic star contraction applied to Asian and Papuan mtDNA evolution. Mol Biol Evol18(10): 1864–1881. pmid:11557793
- 50. Russell AL, Medellin RA, McCracken GF (2005) Genetic variation and migration in the Mexican free-tailed bat (Tadarida brasiliensis mexicana). Mol Ecol 14:2207–2222. pmid:15910338
- 51. Ray N, Currat M, Excoffier L (2003) Intra-deme molecular diversity in spatially expanding populations. Mol Biol Evol 20:76–86. pmid:12519909
- 52. Johns GC, Avise JC (1998) A comparative summary of genetic distances in the vertebrates from the mitochondrial cytochrome b gene. Mol Biol Evol 15(11): 1481–1490. pmid:12572611
- 53. Cárdenas L, Hernández CE, Poulin E, Magoulas A, Kornfield I, Ojeda FP (2005) Origin, diversification, and historical biogeography of the genus Trachurus (Perciformes: Carangidae). Mol Phylogenet Evol 35(2): 496–507. pmid:15804418
- 54. Rousset F (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145(4), 1219–1228. pmid:9093870
- 55. Tzeng TD (2007) Population structure and historical demography of the spotted mackerel (Scomber australasicus) off Taiwan inferred from mitochondrial control region sequencing. Zool Stud 46(6): 656–663.
- 56. Theisen TC, Bowen BW, Lanier W, Baldwin JD (2008) High connectivity on a global scale in the pelagic wahoo, Acanthocybium solandri (tuna family Scombridae). Mol Ecol 7(19): 4233–4247.
- 57. Ely B, Viñas J, Bremer JRA, Black D, Lucas L, Covello K, et al. (2005) Consequences of the historical demography on the global population structure of two highly migratory cosmopolitan marine fishes: the yellowfin tuna (Thunnus albacares) and the skipjack tuna (Katsuwonus pelamis). BMC Evol Biol, 5(1), 19.
- 58. Grant WS, Bowen BW (1998) Shallow population histories in deep evolutionary lineages of marine fishes: Insights from sardines and anchovies and lessons for conservation. J Hered 89(5): 415–426.
- 59. Avise JC, Neigel JE, Arnold J (1984) Demographic influences on mitochondrial DNA lineage survivorship in animal populations. J Mol Evol20(2): 99–105. pmid:6433037
- 60. Nelson JS, Hoddell RJ, Chou LM, Chan WK, Phang VPE (2000) Phylogeographic structure of false clownfish, Amphipriono cellaris, explained by sea level changes on the Sunda shelf. Mar Biol 137(4): 727–736.
- 61.
Avise JC (2000) Phylogeography: the history and formation of species. Harvard University Press. 453 p.
- 62. Hewitt G (2000) The genetic legacy of the Quaternary ice ages. Nature 405(6789): 907–913. pmid:10879524
- 63. Bird MI, Taylor D, Hunt C (2005) Palaeoenvironments of insular Southeast Asia during the Last Glacial Period: a savanna corridor in Sundaland? Quat Sci Rev 24(20): 2228–2242.
- 64. Bird MI, Pang WC, Lambeck K (2006) The age and origin of the Straits of Singapore. Palaeogeogr Palaeocl 241(3), 531–538.
- 65. Liu JX, Gao TX, Wu SF, Zhang YP (2007). Pleistocene isolation in the Northwestern Pacific marginal seas and limited dispersal in a marine fish, Chelon haematocheilus (Temminck & Schlegel, 1845). Mol Ecol 16(2): 275–288. pmid:17217344
- 66. Liu JX, Gao TX, Yokogawa K, Zhang YP (2006) Differential population structuring and demographic history of two closely related fish species, Japanese sea bass (Lateolabrax japonicus) and spotted sea bass (Lateolabrax maculatus) in Northwestern Pacific. Mol Phylogenet Evol 39(3): 799–811. pmid:16503171
- 67. Hanebuth T, Stattegger K, Grootes PM (2000) Rapid flooding of the Sunda Shelf: a late-glacial sea-level record. Science 288(5468): 1033–1035. pmid:10807570
- 68. Lambeck K, Esat TM, Potter EK (2002) Links between climate and sea levels for the past three million years. Nature 419: 199–206. pmid:12226674
- 69. Voris HK (2000) Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations. J Biogeogr 27(5): 1153–1167.
- 70. Tjia HD, Liew KK (1996) Changes in tectonic stress field in northern Sunda Shelf basins. Geol Soc London Spec Publ 106(1): 291–306.
- 71. Mcmillan WO, Palumbi SR (1995) Concordant evolutionary patterns among Indo-West Pacific butterfly fishes. Proc R Soc B 260(1358): 229–236. pmid:7784441
- 72. Arnaud S, Bonhomme F, Borsa P (1999) Mitochondrial DNA analysis of the genetic relationships among populations of scad mackerel (Decapterus macarellus, D. macrosoma and D. russelli) in South-East Asia. Mar Biol 135(4): 699–707.
- 73. Santos MD, Lopez GV, Barut NC (2010) A pilot study on the genetic variation of eastern little tuna (Euthynnus affinis) in Southeast Asia. Philipp J Science 139(1): 43–50.
- 74. Palumbi SR (2003) Population genetics, demographic connectivity, and the design of marine reserves. Ecol Appl 13(sp1): 146–158.
- 75. Wright S (1931) Evolution in Mendelian populations. Genetics16(2): 97–159. pmid:17246615
- 76.
Franklin IR (1980) Evolutionary change in small populations. In: Soule ME, Wilcox BA, editors. Conservation Biology: An Evolutionary-Ecological Perspective. Sinauer Associates, Massachusetts. pp. 135–150.
- 77.
Frankel OH, Soulé ME (1981) Conservation and evolution. Cambridge University Press. 329 p.
- 78. Wang L, Meng Z, Liu X, Zhang Y, Lin H (2011) Genetic diversity and differentiation of the orange-spotted grouper (Epinephelus coioides) between and within cultured stocks and wild populations inferred from microsatellite DNA analysis. Int J Mol Sci 12(7): 4378–4394. pmid:21845084
- 79. Shankar D, Vinayachandran PN, Unnikrishnan AS (2002) The monsoon currents in the north Indian ocean. Prog Oceanogr 52(1): 63–120.