All specimens examined for the present study were independently captured by commercial fishers prior to being sub-sampled for DNA. Some of the sampling localities in our survey (Fig. 1) were in countries that have adopted a principle of state sovereignty over biological resources, including genetic resources with no commercial use [26]. This includes the Philippines, where sampling took place under a collaborative research agreement between the National Taiwan University, Taipei and the University of San Carlos, Cebu; and Indonesia, within the framework of a memorandum of understanding between LIPI, Jakarta and IRD, Marseille. Some of the samples were from West Papua, a former Dutch colony currently administered by Indonesia [27,28] that has no control over its own natural resources [29]. For equitable sharing of results and knowledge, the data produced in the course of our study will be made accessible from the following open-access repositories: Hal-IRD (http://www.hal.ird.fr/; hal@ird.fr), GenBank (http://www.ncbi.nlm.nih.gov) [30], and DIPnet (http://www.indopacificnetwork.wikispaces.com/).

A total of 362 individuals, including 341 blue-spotted maskray and 21 New Caledonian maskray N. trigonoides, were utilized in the present study. This total includes 203 individuals whose CO1 or cytochrome b gene sequences were compiled from either the literature [14,16,17,19,31–38], the BOLD barcode database, the Cryobank DNA-barcode database (http://cryobank.sinica.edu.tw/chi/), or the GenBank nucleotide-sequence database; and 159 new individuals sampled from local fish-landing places throughout the Indo-Malay archipelago, in eastern Africa, and in New Caledonia. All individuals included in the present survey are listed in Supplementary Table S1, along with sampling details, GenBank accession numbers, voucher specimen details, and references. An overview of the density of blue-spotted maskray samples across the Indo-West Pacific, with a focus on the Coral Triangle region is given in Fig. 1.

A piece of skin or flesh ∼0.05 to ∼1 cm³ was removed, using surgical scissors, from the pelvic fin, or the pectoral fin, or the tail and was preserved in 95% ethanol at ambient temperature. DNA extraction was done using either the Viogene (Taiwan) tissue genomic DNA extraction protocol, or the DNEasy DNA extraction kit (Qiagen GmbH, Hilden, Germany). DNA was stored in 1X, pH 8.0 TE buffer (AppliChem, Darmstadt, Germany). Polymerase chain reaction (PCR) amplification of a fragment of the CO1 gene was done according to [19]. For the PCR-amplification of the cytochrome b gene, we used primers L14735 and CB7 [39,40]. The reaction volume was 15 μL and the reaction mixture contained 0.2 mM dNTPs, 1.5 μL 10× PCR buffer (Bioman, Taipei), 0.5 μM of each forward and reverse primers, 0.2 U Taq DNA polymerase (Bioman), and 1.0 μL template DNA. The PCR parameters were: initial denaturation at 94 °C for 4 min, followed by 35 cycles of denaturation (94 °C for 45 s), annealing (48 °C for 1 min), and extension (72 °C for 1 min), and a final extension step at 72 °C for 10 min. Intron 5 of the growth hormone gene, Gh5, was PCR-amplified using exon-anchored primers Gh5F (5′- A G G C C A A T C A G G A C G G A G C -3′) and Gh6R (5′-T G C C A C T G T C A G A T A A G T C T C C -3′) [41] setting the annealing temperature to 64.5 °C and the number of cycles to 35.

The PCR products were sequenced directly using the same primers as those used for the PCR. Sequences were analysed in an automated ABI Prism 377 sequencer (Applied Biosystems, Foster City, CA) at the sequencing facility of the Taiwan Normal University (Taipei).

Nucleotide sequences, including those obtained in the present study, and homologous sequences retrieved from the BOLD, Cryobank, and GenBank databases (see Supplementary Table S1) were aligned using BioEdit [42]. Four datasets were analysed: (1) the CO1 gene sequence dataset, consisting of 330 blue-spotted maskray and N. trigonoides sequences aligned over a maximum length of 722 base pairs (bp); (2) the cytochrome b gene sequence dataset, consisting of 165 blue-spotted maskray and N. trigonoides sequences aligned over a maximum length of 1142 bp; (3) the concatenated CO1 + cytochrome b gene sequence dataset comprising 127 blue-spotted maskray and 6 N. trigonoides sequences; (4) the Gh5 intron sequence dataset, consisting of 18 blue-spotted maskray and N. trigonoides sequences aligned over a maximum length of 359 bp.

The phylogeny of concatenated CO1 + cytochrome b gene haplotypes was inferred by using the maximum-likelihood (ML) method under Mega6 [43]. The most likely nucleotide-substitution model, which was determined according to the Bayesian information criterion, was the Tamura–Nei model [44] where a discrete Gamma distribution (Γ  =  0.76) was used to model evolutionary rate differences among nucleotide sites and invariable sites were allowed. The ML tree was rooted by choosing New Caledonian maskray N. trigonoides as outgroup [16]. The robustness of nodes in the tree was tested by bootstrap resampling.

The CO1 and cytochrome b gene sequence datasets (Supplementary Table S1) were used separately to assign individuals that were sequenced at only one of either locus to a clade uncovered from the phylogenetic reconstruction described above. For this, an ML tree was constructed using each of the two datasets under Mega6 and the assignment of a haplotype to a clade was determined visually from its relative position on the tree.

The Automatic Barcode Gap Discovery (ABGD) algorithm proposed by Puillandre et al. [45] was used to identify gaps in the distribution of pairwise nucleotide distances. The ABGD algorithm detects the largest significant gap in the distribution of pairwise nucleotide distances and uses it for partitioning the dataset. Gap detection is then recursively applied to previously obtained groups to get finer partitions until there is no further partitioning [45]. The ABGD analysis was run on a matrix of 127 CO1 + cytochrome b gene sequences, trimmed to a core length of 1415 bp. Kimura's [46] two-parameter substitution model was selected and the minimum barcoding gap width was set to 1%.

Branch lengths between species are determined by speciation and extinction rates whereas branch lengths within a species reflect coalescent processes at the level of populations [47]. This results in the distinction of species when mitochondrial clades are substantially deeper than the coalescent haplogroups at the extremities of the tree. We ensured that each of the previously identified deep mitochondrial clades was genetically distant from its nearest neighbour by several times the mean intra-clade distance. For this, pairwise intra-clade and between-clades genetic distances were estimated on the concatenated CO1 + cytochrome b gene sequence matrix using Mega6.

To delimit potential cryptic species, we applied the general mixed Yule-coalescent (GMYC) algorithm [47] as implemented in the program Splits under R [48,49]. This algorithm uses an ultrametric tree to differentiate branching patterns consistent with coalescent vs. speciation processes, which provides a threshold time for the transition between populational and species-level processes. We selected the single-threshold option. The ultrametric tree was constructed using Beast v. 1.7 [50] based on the same matrix (127 individuals × 1415 bp) as that used for the ABGD analysis. Outgroup N. trigonoides was excluded from the analysis, which focused on the blue-spotted maskray exclusively. A non-calibrated relaxed lognormal clock, and a constant coalescent tree were assumed. The model of nucleotide substitution was TN93 + G + I (see Section 2.4). Ten million generations of Markov-Chain Monte Carlo analysis were run, of which we sampled a total of 10,000 genealogies. Length of burn-in was determined by visual inspection of traces in Tracer v. 1.5 [50]. A single analysis was run, for which convergence was reached rapidly. The TreeAnnotator v. 1.7. software [50] was used to produce the final tree, based on maximum clade credibility and mean node height.

Nucleotide distances between major blue-spotted maskray clades were several times higher than distances within, confirming previous observations based on preliminary geographic sampling [14,15,17,19]. Denser geographic sampling (present study) did not lead to a dramatically altered ratio of genetic distance between clades to that within, compared to Ward et al.’s results [14]. The gaps in nucleotide distance remained consistent with the existence of 9 main lineages within the blue-spotted maskray. The GMYC algorithm determined that the coalescence patterns within a lineage vs. patterns between generally conformed to the expectations of infra-specific vs. inter-specific processes.

In sharks and rays, short mitochondrial DNA sequences at the CO1 locus (or CO1 DNA barcodes) have proven their suitability to identify up to 99% species [14]. This high barcoding success might be related to the apparent absence of mitochondrial introgression in this vertebrate group [14], itself pointing to a likely low incidence of inter-specific hybridization. This may be an indirect consequence of complex pre-mating behaviour, which prevents hetero-specific pairs, or of mismatch in the morphologies of male and female genitalia, which hampers copulation. Complex pre-mating behaviour has been reported in sandtiger shark [52] and New Caledonian maskray [53]. Mismatch of male–female copulatory apparatuses has been hypothesized for whiprays of the Himantura uarnak (Forsskål, 1775) [54] species complex [35]. Other hypotheses for pre-mating isolation involve chemical or visual cues [17]. Because cases of mitochondrial introgression are virtually unknown in sharks and rays [6], one expects that patterns of inter-specific genetic differentiation at nuclear loci to be parallel to those at the mitochondrial locus.

Nuclear markers have also been used to characterize blue-spotted maskray populations [17,18]. Significant allele-frequency differences at 2/4 intron loci were observed between adjacent populations from southern Java and Bali Strait [18], corroborating the distinction between previous Clades IVb and VI ([19,23]; present study), hence the recognition of, respectively, N. orientale and N. caeruleopunctata as separate species. Two major clades were observed at the nuclear locus RAG-1 [17] (Fig. 3B), that distinguished the lineage from the eastern Andaman Sea (our Clade III) and those from the adjacent South China Sea (N. varidens) and Java Sea (N. orientale). Similarly, the RAG-1 marker allowed the distinction between the lineage from southern Bali (N. caeruleopunctata; our previous Clade VI) and that of the adjacent Lombok Island (N. australiae; our previous Clade V) [17]. Additional information from nuclear marker Gh5 was produced in the present study. Nucleotide sequences at intron locus Gh5 confirmed the distinction between N. orientale on the one side, and N. australiae and N. caeruleopunctata on the other side, although some other Gh5 sequences of individuals determined as N. orientale according to their mitochondrial haplotypes appeared to cluster with N. australiae and with Clade I. This may either signal mitochondrial introgression, or reflect incomplete lineage sorting at the Gh5 locus. In addition, at this locus, samples from the Moluccas and West Papua (respectively, our lineages VII and VIII) were distinct from all the other samples, and also from each other.

Repeated lowering of the sea level caused by global temperature oscillation [55] has led to the repeated fragmentation of the marine habitats of the Indo-Malay archipelago in the Pleistocene [56]. In shallow-water species like the blue-spotted maskray, this has caused repeated depletion and fragmentation of the populations, followed by repeated recolonisation of shallow-water habitats once the sea level rose again, and secondary contact between long-isolated populations. Population collapse reduces genetic diversity by the random fixation of a subset of alleles and subsequent population expansion favours new, nascent diversification. This succession of isolation episodes with cyclic population collapses, followed by rapid population expansions, is thought to have favoured speciation in the blue-spotted maskray [17,19].

Secondary contact after recolonisation of inundated shallow-water habitat may have been a factor reinforcing genetic isolation instead of promoting re-homogenisation. At the scale of an individual's lifetime, the blue-spotted maskray is a sedentary species. Long-term site fidelity has been inferred from tagging experiments in the closely-related N. trigonoides, with similar results for females and males [53]. Relatively poorly dispersing species often form well-delineated parapatric boundaries [57]. One of the proposed mechanisms is that of a narrow and stable hybrid zone, which acts as a geographic barrier to cross-dispersal [57]. While this hypothesis remains to be tested in the blue-spotted maskray species complex, the present study identified several of these parapatric boundaries, on which future genetic studies should focus.