Brown banded bamboo shark (Chiloscyllium punctatum) shows high genetic diversity and differentiation in Malaysian waters

The demersal brown banded bamboo shark Chiloscyllium punctatum is a major component of sharks landed in Malaysia. However, little is known about their population structure and the effect of high fishing pressure on these weak swimming sharks. Both mitochondrial DNA control region (1072 bp) and NADH dehydrogenase subunit 2 (1044 bp) were used to elucidate the genetic structure and connectivity of C. punctatum among five major areas within the Sundaland region. Our findings revealed (i) strong genetic structure with little present day mixing between the major areas, (ii) high intra-population genetic diversity with unique haplotypes, (iii) significant correlation between genetic differentiation and geographical distance coupled with detectable presence of fine scale geographical barriers (i.e. the South China Sea), (iv) historical directional gene flow from the east coast of Peninsular Malaysia towards the west coast and Borneo, and (v) no detectable genetic differentiation along the coastline of east Peninsular Malaysia. Genetic patterns inferred from the mitochondrial DNA loci were consistent with the strong coastal shelf association in this species, the presence of contemporary barriers shaped by benthic features, and limited current-driven egg dispersal. Fine scale population structure of C. punctatum highlights the need to improve genetic understanding for fishery management and conservation of other small-sized sharks.

Sharks are generally highly vulnerable to overexploitation due to their life history strategies, such as late maturity and low fecundity [1][2][3][4][5] . Overexploitation has been shown to reduce genetic diversity and increase extinction risk especially for small populations 6 ; therefore understanding the genetic structure and migration patterns or gene flow of shark species are essential to inform effective management and conservation plans. Earlier shark genetic studies had been largely focused on marine neritic shark species, that are charismatic, economically important and of conservation interest 5 , e.g. blacktip shark C. limbatus 7 , hammerhead shark Sphyrna lewini 8 , whale shark Rhincodon typus 9 and great white shark Carcharodon carcharias [10][11][12][13][14][15] . Studies on small-sized benthic coastal sharks have increased in recent years in part due to their catch prominence and importance in coastal fisheries, both globally and in Southeast Asia, e.g. whitespotted bambooshark Chiloscyllium plagiosum 16 , nurse shark Ginglymostoma cirratum 17 , whitetip reef shark Triaenodon obesus 18 , leopard shark Triakis semifasciata 19 and common smoothhound Mustelus mustelus 20 .
One such group of sharks is the longtail carpet shark or bamboo shark from the family Hemiscyllidae (order Orectolobiformes). The family Hemiscyllidae comprises two genera: Hemiscyllium Müller & Henle, 1838 (nine species) which is confined to the Australia-New Guinea region and Indonesia, and the Indo-Pacific genus of Chiloscyllium Müller & Henle, 1837 (eight species) 21 . These bottom-living sharks are known to be weak swimmers 22 with great camouflage ability, allowing them to adapt to demersal habitat by hiding around crevices 23 . They live either solitarily or within a group and appear to be strongly territorial. They usually deposit egg cases on the sea floor (sheltered or deserted), with some cases attached to benthic marine plants 23 . Currently, their fecundity in the wild is still unknown, but studies on captive individuals suggested that they produce annually less than 30 egg cases that require a three-month incubation period [24][25][26][27] . To the best of our knowledge, no movement study on this animal exists. Dispersal is likely very restricted due to their oviparous nature 28 , thus rendering them highly vulnerable to localized exploitation and habitat destruction.
Presently, four species of Chiloscyllium are found in Malaysia: Indonesian bamboo shark C. hasseltii Bleeker, 1852; slender bamboo shark C. indicum (Gmelin, 1789); whitespotted bamboo shark C. plagiosum (Bennett, 1830); and brown banded bamboo shark C. punctatum Müller & Henle, 1838. A previous record of grey bamboo shark C. griseum Müller & Henle, 1838 in Malaysian waters was a misidentification of C. hasseltii Ahmad, personal communication . The brown banded bamboo shark C. punctatum has been classified as Near Threatened under the IUCN Red List but may be uplisted to Vulnerable in the near future in Southeast Asia due to high fishing pressure 29 . It is the top ten shark species landed by numbers in commercial fisheries across Malaysia 30 . This species makes up almost half of all shark individuals caught using bottom trawl fishing vessels in the west coast of Peninsular Malaysia 31 . Approximately 40% of C. punctatum landed in commercial fishing gears comprised of immature individuals, with higher percentage in the west than the east coasts of Peninsular Malaysia 31,32,Lim, unpublished data . Among similar-sized sharks, C. punctatum has relatively lower commercial value, fetching an ex-vessel price of USD 0.20-0.60 per kg. Many artisanal fishers in Malaysia generally do not discard any of their non-trash fish catches; when they do, bamboo sharks are among those selectively discarded Then, unpublished data . No formal fishery assessment exists for this species (or other sharks) in Malaysia, and the effect of present fishing level on the animal is unknown.
One of the priorities for C. punctatum includes clarifying the population substructure due to concerns of population fragmentation 29 likely due to being a weak swimmer and limited potential for dispersal. Genetic tools have been used to identify population vulnerability or risk of local extinction based on their genetic diversity and structure as well as connectivity among populations [33][34][35] . Intra-specific genetic studies on small-sized coastal sharks and rays with limited vagility showed that they display high genetic population differentiation across ocean basins due to environmental barriers (nurse shark Ginglymostoma cirratum 17 ), vicariance and dispersal events (maskray Neotrygon kuhlii 36 ), or a combination of both (small-spotted catshark Scyliorhinus canicula 37 ). On the other hand, studies of these small coastal species along continuous continental coastlines either displayed clear genetic differentiation (leopard shark Triakis semifasciata 19 , short-tailed stingray Dasyatis brevicaudata 38 ) or lack of genetic structure (milk shark Rhizoprionodon acutus 39 ); reasons for these pattern variations included differences in habitat preferences, site fidelity, geographical distance, and physical barriers to movement. Despite growing concerns of shark exploitation in Malaysia as one of the top shark-fishing nation globally 40 , genetic studies of local shark populations, including C. punctatum are almost non-existent (see Dudgeon et al. 41 for exception).
The genetic structure of any organism is influenced not only by the ecological processes but also the geological history [42][43][44][45][46] . Sundaland, a shelf within the biodiverse Indo-West Pacific region that extends from Thailand southwards primarily covering Malaysia and Indonesia, has undergone dramatic plate tectonic evolution forming episodic submergence and exposure of the shelf 47,48 . These series of events had affected distributions of local organisms (ranging from terrestrial to marine environment) in multiple ways due to different population dispersal patterns (e.g. Lim et al. 45 ; Tan et al. 46 ; Reid et al. 49 ; Polgar et al. 50 ; Leonard et al. 51 ; Ma et al. 52 ; Mason et al. 53 ; Crandall et al. 54 ). In relation to elasmobranch populations, the Sunda Shelf appeared to be acting as biogeographic barrier for wedgefish Rhynchobatus australiae 55 but not for spottail shark Carcharhinus sorrah 56 . However finer-scale seascape features that may affect genetic structuring of sharks and rays within the Sundaland region had not been investigated previously.
The aim of the present study is therefore to assess the degree of genetic diversity and differentiation of C. punctatum in Malaysia, specifically along coastal Peninsular Malaysia and northern Malaysian Borneo coastline. Based on the biogeographic history of the Sundaland region, we hypothesized that there are at least two distinct local C. punctatum populations, with little genetic mixing between populations from either side of the peninsular (northern end of the Sunda Shelf barrier) and across the Karimata Strait that separates the peninsular from Borneo. Furthermore, we tested to see if observed population genetic differentiation can be explained by simple geographical distance or otherwise. In addition, we also assessed the gene flow directionality between distinctive populations to infer putative dispersal pattern. Knowledge of fine scale genetic structuring of C. punctatum provides the spatially relevant information for appropriate national fisheries management and conservation of this highly fished shark and likely of other sharks sharing similar life history strategies.

Results
A total of 2116 base pair sequences hereby named as ND2CR (concatenated between 1072 bp for CR and 1044 bp for ND2were successfully amplified for all 135 C. punctatum samples. A total of 70 unique haplotypes and 63 polymorphic sites were detected, in which 40 sites were found in the CR region and the rest were found in the ND2 region (Supplementary Tables S1 and S2). Most of the identified haplotypes (70%) was found only in a single individual, while the rest was represented by two to 17 individuals. In terms of genetic diversity, ND2 gene showed lower range (haplotype diversity: 0.58-0.71, nucleotide diversity: 0.0009-0.0026) compared to CR (haplotype diversity: 0.63-0.90, nucleotide diversity: 0.0021-0.0076) ( Table 1). When the two markers were concatenated, the range of haplotype diversity increased (0.76-0.98) ( Table 1).
The Tajima's D tests were all not significantly different from zero which suggested no evidence of departure from neutrality (Table 1). On the other hand, Fu's Fs which is more sensitive to recent population expansion 57 suggested some evidence for population expansion with significantly negative values for both CR and concatenated markers for WP, EP and SR populations ( Table 1). The mismatch distribution analysis for the individual and concatenated markers showed consistent results: across all sampling areas; only EP had smooth and unimodal distributions that closely fitted the expected distribution under the sudden expansion model (Fig. 1). The distribution for total samples and all other sampling areas demonstrated a multimodal distribution that is typically attributed to populations in demographic equilibrium.
The haplotype network showed that most haplotypes were found almost exclusively within an area (Fig. 2). Only five haplotypes were found across multiple sampling areas, i.e. H14 was found in EP, WP and SR; H33 in EP and SR; H40 in SR, WS and ES; and both H66 and H68 in WS and ES (Supplementary Table S3). In general, four major haplotype clusters were found among the samples-one with haplotypes found only in WP, one with haplotypes found mostly in EP; one with mixed haplotypes from SR, WS and ES; and one with haplotypes from both WS and ES. The single sample obtained from GenBank (sampled in China) clustered together with haplotypes from EP. For samples from Sabah (WS and ES), four sub-clusters separated by a minimum of five mutations were observed. This includes one large cluster with nine haplotypes (seen in 18 individuals from WS, 14 from ES and one from SR) and three minor clusters each with five (three from WS and 14 from ES), three (three from WS and one from ES) and one (one from WS) haplotypes.
Overall, there was a weak but significant support for IBD between areas (r 2 = 0.62, p < 0.05), with a broad range of pairwise genetic distances at all geographical distance concentrated mostly along the regression line (Fig. 3a). Pairwise comparison supported the high genetic differentiation (p < 0.05) found between all pairs of sampling areas, including weak significance between the two areas within Sabah (WS and ES) ( Table 2). Results from the hierarchical AMOVA demonstrated 40-50% of genetic variation was attributed to differences among sampling areas (Table 3). When testing for genetic partitioning under various biogeographical barrier combinations, significant hierarchical partitioning of genetic variation at all levels for both markers was observed when only South China Sea was considered as a genetic barrier (the B2 only scenario). Under this scenario, 46.2% of genetic variation was partitioned among population groups, 15.73% among populations within group, and 38.07% within populations.
Migrate-n analyses with 10,000 steps and 50,000 steps reached convergence and were both in agreement on the direction of migration between pairs of sampling areas across all pairwise comparisons, with exception of the SR-WS pair (Table 4). An additional analysis using 100,000 steps concurred with the results using 50,000 steps for the SR-WS pair. The best models showed that gene flow occurred from EP to WP and from EP to SR. Specifically for C. punctatum sampled from Borneo, model results showed greatest support for gene flow directionality from SR to WS and from WS to ES (Table 4, Fig. 4, Supplementary Tables S4, S5).
To determine the genetic structure and dispersal pattern at a more local scale, we assessed the level of genetic differentiation, haplotype network and the presence of Isolation-By-Distance (IBD) among the a priori subpopulations along the coastline of EP. No significant genetic differentiation was found among subpopulations (p > 0.05) ( Table 5). Congruent with this finding, the EP haplotype network showed no clear patterns of geographical segregation of haplotypes (Fig. 5). For instance, H18 was shared among EP2, EP3, and EP5; H25 from EP1 was one mutation away from H26 from EP5 (Supplementary Tables S1 and S2). Furthermore, IBD analysis Table 1. Genetic diversity of Chiloscyllium punctatum according to area and markers. N = number of samples, k = number of polymorphic sites, N ha = number of haplotypes, ha = haplotype diversity, π = nucleotide diversity, D = Tajima's D test, F S = Fu's F S test. *Represent significant difference at p < 0.05. www.nature.com/scientificreports/ showed no correlation between genetic and geographical distances for these subpopulations (r 2 = 4.73 × 10 -5 , p > 0.05) (Fig. 3b).

Discussion
Using mitochondrial DNA markers for the bamboo shark C. punctatum, we were able to show (i) strong genetic differentiation with little present day mixing between the five areas (populations) sampled, (ii) high intra-population genetic diversity with unique haplotypes, (iii) significant correlation between genetic differentiation and geographical distance coupled with detectable presence of fine scale geographical barriers (i.e. the South  www.nature.com/scientificreports/ www.nature.com/scientificreports/ China Sea), (iv) historical directional gene flow from EP bifurcating towards WP and towards Borneo (from the western end of SR to Sabah), and (v) no detectable genetic differentiation along the EP coastline. The high genetic structure observed in C. punctatum within the central area of its distribution range is also consistent with their strong coastal association and life history strategies. We discuss these results in relation to biological and geographical barriers that shape the animal's phylogeography and implications for fishery management of this small-sized shark.
Genetic diversity. The large number of unique mtDNA haplotypes in C. punctatum found within a small region had exceeded those reported for other elasmobranchs sampled across broader geographical range, e.g. 26 haplotypes in Negaprion brevirostris 58 , 39 in Carcharhinus sorrah 56 , and 12 in Ginglymostoma cirratum 17 .   Table 4 for actual values). www.nature.com/scientificreports/ The exception is a report of 67 haplotypes found in C. plumbeus 59 . Analysis using cumulative haplotype curves suggested that there are many haplotypes remaining to be discovered for C. punctatum ( Supplementary Figure S1). Compared to similar sized species of sharks 60 , haplotype diversity of C. punctatum is relatively high while the nucleotide diversity is comparable. High number of haplotypes seen in C. punctatum could be partially explained by the haplotype genealogy 61 and the relationship between life history strategies and mutation rates 62 . Star-shaped haplotype genealogy and high haplotype diversity with low nucleotide diversity is a typical genetic signature of recent population expansion, especially when a new habitat became available 61 . Evidence from this study showed some support for selective directional population expansion for C. punctatum although the demographic history may not be shared similarly for all the populations in Malaysia (discussed further below). The high haplotype diversity and directional population expansion detected in our C. punctatum samples underscores the need for careful fisheries management to maintain the natural diversity of this highly exploited  www.nature.com/scientificreports/ species and the role of the Malaysian stock as a source population. To date, there is no study that directly assess the relationship between overfishing and genetic diversity in shark and ray species 5 . Overharvesting of marine fish species has been shown to result in reduced mean body size, faster growth rate, and earlier sexual maturity 63 , as well as lower genetic diversity, especially allelic richness 64 . On the other hand, the level of genetic diversity in fishes was surprisingly not found to be lower for threatened species relative to non-threatened ones 65 . Our findings thus provide an important baseline data for key genetic diversity metrics that can be used as indicators of population health and overall conservation status of C. punctatum.
Genetic structure and phylogeography. The lack of genetic structure along the coastline of EP, absence of IBD and sharing of haplotype H18 among EP subpopulations suggests high gene flow and that mixing of C. punctatum genetic information could occur as far as 400 km, the longest distance between two known catch locations of the trawl survey conducted. The coastwide genetic mixing can either be facilitated by adult movement or dispersal of egg cases. No mobility study had been conducted on C. punctatum or other congeners to date. However, a closely related species with similar behaviour, the nurse shark G. cirratum (total length, TL = 430 cm) is able to travel at least 541 km 66 and genetic studies on the species indicated possibility of longer distance travel along a ~ 5000 km coastline 17 . Another continental shelf species, South African catshark Poroderma africanum (TL = 101 cm) could travel as far as 1964 km 66 .
Since coastal depth is probably one major restriction for adult C. punctatum movement, ancient and present day bathymetry maps of this region suggest that the species should be able to travel along the shallow waters of coastal EP connected to WP along the Johore Strait (also known as the Tebrau Strait, with maximum depth of 16.9 m) or Singapore Strait and through the Karimata Strait to the Borneo side [67][68][69] . However significant genetic differentiation and weak support for IBD in C. punctatum across the five sampling areas indicates that benthic seascape features may form substantial barriers for adult shark movement across EP-WP and EP-SR. The barrier between EP and WP could result from intensive anthropogenic activities at the Johor and Singapore Straits. These straits, though only approximately 50 km in length, are well known for major reclamation and busy shipping traffic since 1819 70 . Furthermore, industrial activities within the Straits are likely to impact benthic habitat availability and deter adult movement between either coasts of Peninsular Malaysia. The Karimata Strait between EP and SR had been reported to have strong unpredictable currents 71,72 , thus potentially limiting adult movement across the Strait between EP and SR at the western end of Borneo. Study on movement and benthic habitat use for C. punctatum will be useful to corroborate the role of distance and bottom features that shape present day genetic structure.
The westward gene flow directionality observed from EP to WP and eastward flow from EP to Borneo (Fig. 4) could be partially explained by contemporary surface ocean currents. During the north-east monsoon, ocean current flows from EP towards WP and towards SR at about 0.3-0.5 m/s 73 ; this would be most likely consistent with a pelagic stage dispersal. Egg cases for various elasmobranch species had been recorded to have holding abilities, adhesive fibre in the case of C. punctatum, to prevent drifting, exposure to predation or damage by wave and tidal action [74][75][76] . However, ROV video footage capturing movement of egg cases along with water currents lend support for the possibility of egg case dispersal as a mode of limited population connectivity across these areas 77 . The gene flow directionality implies that population in EP could act as a source for populations in Borneo and in WP.
Unexpectedly, significant genetic differentiation in C. punctatum could also be observed along the northern Borneo coastline of about 385 km. Detection of a biogeographical barrier between Sarawak and Sabah suggests that benthic seascape features coupled with differing water bodies (South China Sea and Sulu Sulawesi Sea) are important in shaping the genetic structure for C. punctatum of Borneo and possible other similar sharks. This physical delineation is sufficiently significant that Scoliodon macrorhynchos, a small sized spadenose shark, could be found along the coastline of the Sarawak Borneo but not of Sabah Then, unpublished data . In relation to the weak genetic structuring observed between eastern and western Sabah, further investigation is warranted. Discussions with local fishers and other researchers suggest that sharks landed in Sabah, especially by trawlers, could have been caught from the waters of the Philippines (from the Sulu-Sulawesi waters) Thus, the genetically differentiated Sabah populations and their transboundary relationship with putative populations from the Philippines should be jointly considered.
The genetic differentiation and possible directional population expansion observed in this study could reflect historical vicariance events in the Sunda Shelf region. Other marine organisms in the Indo-Malaya archipelago showed genetic separation between the east and west peninsular populations. These range from smaller and less mobile species such as gastropods Nerita albicilla 78 and Echinolittorina spp. 49 , echinoderm Diadema setosum 54 , small bony fishes (mudskipper Periophthalmus argentilineatus, and reef fish Centropyge loriculus) 50,54 , to larger and more mobile species like wedgefish Rhynchobatus australiae 55 . The congruency between our data and these past studies indicated the importance of the Sunda Shelf barrier in shaping the phylogeography of marine species in this region.
Overall, the genetic patterns and differentiation observed for C. punctatum inferred from the mitochondrial DNA loci would reflect signatures that occurred over historical time scale rather than during contemporary ecological time scale. However, an additional possibility to the strong genetic structure observed can be due to behavioural traits, for instance philopatry or site fidelity which has been shown for diverse species of sharks 79 ; thus this cannot be ruled out for the strongly coastal-associated C. punctatum. The joint use of nuclear and mitochondrial markers will be important to discern whether the observed patterns of genetic structure in C. punctatum could be attributed to philopatric tendencies or due to the presence of biogeographic barriers, or both. www.nature.com/scientificreports/ Fisheries management implications. Chiloscyllium punctatum experiences high fishing pressure despite being caught primarily as non-target species (bycatch). In 2017 alone, we estimated that approximately 670,000 individuals were caught and landed at the four major shark landing states in Malaysia (namely Perak, Pahang, Sarawak, and Sabah), based on estimated landings of 750 metric tonnes of C. punctatum and average weight of C. punctatum landed in trawl fisheries of 1.12 kg 31,80, Lim unpublished data . The national demersal resources were estimated to have decreased by 80% relative to assumed virgin (unexploited) levels during 1960s and this decrease was attributed to overfishing that had reportedly been ongoing for the past two decades 81 . As this species make up almost half of all sharks caught from trawlers, and almost half of the individuals landed were juveniles, impact of fishing pressure remains a great concern in the long run. High level of genetic structuring in C. punctatum across Malaysian waters coupled with low levels of interregion connectivity strongly suggests the need to consider separate fishery management stocks and conservation plans for the west coast of Peninsular, east coast of Peninsular and Borneo. These separate populations appeared to rely almost exclusively on self-recruitment (with F ST mostly > 0.25 indicating low migration rate among areas) 35,82 ; migrants would also likely find difficulty transmitting their genes locally given the territorial nature of the species. Localized population extinctions may be possible if present fishing pressure is not curbed.
Despite ongoing high fishing pressure on highly differentiated genetic stocks of C. punctatum, unique biological traits likely confer considerable resilience to this species. Chiloscyllium punctatum has shown the ability for long term sperm storage 83 , is extremely hardy and physiologically adapted to inhabiting environments that undergo cyclical hypoxic conditions (e.g., coral reef flats 84 ). However, this is likely not the case for many other coastal-associated shark and ray species within the highly fished waters of Malaysia and other neighbouring countries. In addition to clarifying the genetic structuring of these exploited species, efforts to reduce present fishing effort should be prioritized. The proposed move towards a total ban of bottom trawl fishing within Malaysian waters in the next few years would substantially reduce fishing pressure on sharks and rays.

Conclusion
We present novel evidence of high genetic structure and diversity of C. punctatum within the Indo-Malaya Archipelago. Observed high genetic differentiation between the five sampled areas could be attributed to the Sunda Shelf formation acting as ancient barrier for the animals. This genetic pattern persists in the present day due to their strong coastal shelf association and likely limited adult movement, apart from limited directional gene flow between areas. Presence of contemporary barriers for population connectivity is shaped by both distance and benthic features and limited current-driven egg dispersal. Although present high fishing pressure do not appear to pose immediate concerns for C. punctatum, findings of strong genetic structure within the small geographical location considered raises the concern of the population health of other exploited small sized sharks in the area. Sampling site. Malaysia is located in the centre of the Indo-West Pacific biogeographic region 85 as well as the distribution range of C. punctatum 86 . It is surrounded by three major water bodies namely the Strait of Malacca, southern South China Sea and the Sulu-Celebes Sea (Fig. 6). Sampling was conducted at five major coastal areas: west coast of Peninsular Malaysia (WP), east coast of Peninsular Malaysia (EP), Sarawak (SR), western Sabah (WS) and eastern Sabah (ES) (Fig. 6). The Sabah samples were sub-divided into two areas representing different water bodies; WS faces the South China Sea while ES faces the Sulu-Celebes Sea. These sites were selected to examine the effect of biogeographical barriers to C. punctatum populations defined a priori: historical Sunda Shelf Barrier between WP and EP (B1), depth and distance barrier (South China Sea) between EP and SR (B2), and distance barrier between SR and Sabah (B3).

Materials and methods
The fin clips of 25-30 individuals per area were collected with a total of 135 individuals (Table 1). With the exception of EP, all other samples were obtained directly at fish landing sites during a 25-month ichthyofaunal survey from September 2015 to September 2017. These animals were caught via commercial fishing operations; the fishing vessels were mainly trawlers except in SR whereby some individuals were caught with gill nets. A maximum of two individuals per fishing vessel were sampled, to prevent over-sampling of maternally-linked individuals that may result in the reduction of potential haplotype diversity 56,87,88 . Samples from EP were obtained through a fisheries-independent demersal trawl survey within the Exclusive Economic Zone (EEZ) organised by the Department of Fisheries Malaysia from July 2015 to July 2016. This invaluable opportunity allowed us to examine finer-scale genetic structure and connectivity of C. punctatum along the 724 km-long EP coastline. The EP samples were further divided into five subpopulations (EP1 to EP5) along the coastline (Fig. 6). All tissue samples were preserved in absolute ethanol before subsequent molecular procedures.
Mitochondrial DNA extraction and analysis. The DNA of each sample was extracted using 10% Chelex resin incubated for 2 min at 60 °C, followed by 25 min at 103 °C (modified from Hyde et al. 89 ). Two mitochondrial DNA (mtDNA) markers were employed to analyse the genetic diversity and structure of C. punctatum: Control region (CR) and NADH dehydrogenase subunit 2 (ND2). The primer sets used for the targeted regions were CR-F 5′ GAC CTT GTA AGT CGA AGA 3′ and CR-R 5′ TCT TAG CAT CTT CAG TGC 3′ for CR 90  Data analysis. Sequences of the CR and ND2 regions were reviewed manually, aligned and edited using ClustalW 92 implemented in Mega version 7 software 93 . All haplotype sequences used for the following analysis were submitted to Genbank database with accession numbers as in Supplementary Table S3. Cumulative haplotype curve using concatenated markers was constructed by randomising the sample order 999 times. The increment of new haplotypes against the number of sequences analysed was plotted using the PRIMER v6 software 94 . This is to observe the increment pattern of haplotype when new sequences was added in the analysis. An additional sequence from NCBI Genbank (accession number: JQ082337; sampling location: China) that covered both CR and ND2 was included in the analysis to observe potential connection between samples in the South China Sea. Intra and inter-population genetic diversity by individual and concatenated markers (ND2CR) was characterized using D NA SP v6 software 95 with the following indices: the number of haplotypes (N ha ), polymorphic sites (k) and, haplotype diversity (ha) and nucleotide diversity (π). Haplotype or gene diversity is the probability of getting two different alleles from random sampling while nucleotide diversity is the average differences in number of nucleotides per site among DNA sequences 96 . The NETWORK v5.0 software (http:// www. fluxus-engin eering. com was then used to construct the median-joining haplotype network to visualize the genetic connectivity among the populations. AMOVA analysis was conducted in ARLEQUIN v3.5.2.2 97 with 10,000 permutations to identify the hierarchical partitioning of genetic variation under the influence of various combinations of biogeographical barriers stated earlier (B1, B2 and B3). The best model for each marker was selected using jModelTest2; since the best models identified were not implemented in ARLEQUIN, we used the next best available models, i.e., Tamura and Nei correction for CR, Kimura 2 Parameter for ND2 and Kimura 2 Parameter with Gamma 0.656 for ND2CR. Pairwise comparisons of population differentiation using F ST values among (1) subpopulations in EP, and (2) all populations were determined and significant probability values were reported after correction for false discovery rate (FDR) 98 .
Isolation-by-distance (IBD) i.e., correlation between genetic distance and log (base 10) of geographical distances, between populations was calculated and their significance was tested using the Mantel test with 9999 randomizations. Genetic distance between populations was defined as (F ST /(1 − F ST )) 99 . Pairwise geographical distance among populations was estimated as the mean value of the shortest straight-line distance between two points at the sea among the pairs adapting two turning point around southern tip of peninsular to reach WP (west point: 1.1767°N, 103.5197°E; east point: 1.3485°N, 104.3724°E) and one turning point at northern tip of Sabah to ES (7.1231°N, 116.9426°E). For EP, individual-based IBD analysis was performed given that individual catch locations were known and distances for some catch locations within a subpopulation were further than www.nature.com/scientificreports/ those between subpopulations; uncorrected p-distance calculated in PAUP* 4.0b10 software 100 . The IBD analysis was performed using Adegenet version 2.1.1 package in R version 3.5.1 software 101 . Historical demographic expansion was inspected using Fu's 57 Fs and Tajima's D tests 102 of neutrality. Mismatch distribution analysis was subsequently conducted to assess demographic expansion using individual and combined markers in DnaSP v6.10.03 96 . The migration pattern of C. punctatum was determined by assessing the migration rate, population size and migration direction using the Bayesian Inference in Migrate-n v3.6.11 [103][104][105] . The analysis was performed on both individual and concatenated markers (ND2CR). Migrate-n allows for the estimation of relative effective population size, θ Ne = (N e μ), and asymmetric gene flow, M (m/μ), among populations over longer periods of time (> 1000 years) 106 with the assumption of constant population size, random mating within population, constant mutation rate, constant immigration rate and genetic materials can only exchange through migrants (no population divergence is allowed) 1 . Results from both pairwise testing and haplotype network were used to guide the selection of a priori population pairs for statistical testing and to reduce the number of models considered in Migrate-n. Five population pairs (WP-EP, EP-SR, SR-WS, SR-ES and WS-ES) were subjected to four migration models (1 ↔ 2: bi directional; 1 → 2: unidirectional from population 1 to 2; 1 ← 2: unidirectional from population 2 to 1; 1 ≠ 2: no gene flow) under 10,000 and 50,000 recorded Markov chain Monte Carlo (MCMC) steps; additional analyses using 100,000 steps were done when results did not agree between runs of 10,000 and 50,000 steps. Four-chain heating at temperatures of 1, 1.5, 3 and 1,000,000 was applied to increase the efficiency of the MCMC. Other analysis criteria followed the default settings Modelling protocol and convergence detection followed the guidelines by Beerli et al. 107 . Comparison of models was performed using the Bayes factors based on guidelines by Kass and Raftery 108 . www.nature.com/scientificreports/