Delimiting cryptic species in a shark with a low dispersal capability

Background Delimiting cryptic species in elasmobranchs is a major challenge in modern taxonomy due the lack of available phenotypic features. Employing stand-alone genetics in splitting a cryptic species may prove problematic for further studies and for implementing conservation management. In this study, we examined mitochondrial DNA and genome-wide nuclear single nucleotide polymorphisms (SNPs) in the brown-banded bambooshark, Chiloscyllium punctatum to evaluate potential cryptic species and the species-population boundary in the group. Results Our results found four operational taxonomic units (OTUs) within C. punctatum from the Indo-Australian region. Each OTU can be interpreted differently depending on available supporting information. Similarly, we conrmed that comprehensive sampling over the species' geographic distribution was essential to determine the boundary between population and cryptic species. We provide suggestions about what should be considered prior to split cryptic species and the designation of new species.


Results
Our results found four operational taxonomic units (OTUs) within C. punctatum from the Indo-Australian region. Each OTU can be interpreted differently depending on available supporting information. Similarly, we con rmed that comprehensive sampling over the species' geographic distribution was essential to determine the boundary between population and cryptic species.

Conclusion
We provide suggestions about what should be considered prior to split cryptic species and the designation of new species.

Background
The development of genetic and genomic studies has substantially in uenced how species are de ned taxonomically, with a shift from traditional taxonomy based on morphological and biological characters to one that includes or entirely relies on DNA-evidence (e.g. (1)). While the necessity for phylogenetic support in de ning a species is the subject of debate, many journals require phylogenetic analyses in the description of new taxa (2). Advanced DNA sequencing has revealed species complexes and cryptic species previously assumed to be subspecies or populations of a single species through traditional taxonomic analysis (3)(4)(5). In evolutionary theory, species complexes are considered recently diverged from evolving metapopulations (3,6). As an ongoing process of speciation, morphological differentiation may develop later due to local adaptation to new environments (3,4). Cryptic species may be present within a complex when morphologically indistinguishable individuals are revealed to be genetically distinct (6)(7)(8). While there are disparities in the de nition of cryptic species in terms of their biology, they provide challenges for taxonomic and evolutionary studies (9). For instance, the discovery of cryptic species in what was previously considered a single wide-ranging species raises questions whether geographically-structured populations can be differentiated from different species within the same region, leading to taxonomic confusion (10,11). A reliance on genetic evidence alone to split a species complex into two or more new species presents wildlife managers and ecologists with challenges as morphological features that facilitate visual discrimination are absent (12)(13)(14). A lack of eld-based identi cation may confound determination of a species' conservation status, and management strategies that might be relevant for that species (6,15,16).
Cryptic species are common in many groups of taxa and various habitats, including in marine environments (4,6,17). In elasmobranchs, the number of recognised cryptic species has increased substantially within the last decade, primarily due to a large, multi-taxa, phylogenetic study of mitochondrial NADH2 sequences (18). Several putative cryptic species have been identi ed with most originating from areas of high geographic complexity, and in taxa with low dispersal ability (18,19). The use of genome-wide approaches alongside mitochondrial markers has improved the ability to identify cryptic diversity of elasmobranchs (e.g. in mobulids (20)).
In this study, we examined the brown-banded bambooshark, Chiloscyllium punctatum, a species distributed in near-shore environments within the Indian and the Western Paci c Oceans (21)(22)(23)(24)(25) (Fig. 1), and the most commonly caught shark in coastal sheries in Southeast Asia (26). Previous genetic and morphological studies have suggested that C. punctatum contains two cryptic species (18,27). However, the genetic study only included a few sampling locations with uneven sample sizes and large geographic breaks, likely under-representing the species' genetic diversity within the region [17]. Here we use genomewide nuclear single nucleotide polymorphisms (SNPs) along with mitochondrial DNA lineages to assess potential cryptic speciation within an extensive, previously un-sampled area within the Indonesian archipelago, and to evaluate the use of genetic and genomic data to delineate species.

Mitochondrial Gene Phylogeny
In total, 34 individuals were sequenced for the mitochondrial NADH2 marker across 12 locations. Two phylogenetic analyses of the trimmed alignment sequence data (1044 bp) (maximum likelihood phylogenetic tree and median-joining haplotype network) showed similar clusterings of haplotypes. The phylogenetic tree indicated two major clades: (i) samples from Australia (WAU and SEQ), Papua New Guinea (PNG), Lombok (LMB), Muncar (MUN), west coast of Sumatra (WSA and WSS); and (ii) samples from South East Asia including Thailand (PHU), Malaysia (PAH) and Indonesia (BIN, WKL, WJV, SUL).
Within the rst clade, samples from Australia and Papua New Guinea formed a distinct cluster (Australian region) separated from samples from the south coast of Indonesia (the Indian Ocean region). Within the second clade, samples from the Indo-Malay region (PAH, BIN, WKL and WJV) formed a cluster separated from South Sulawesi and Thailand samples (Fig. 2a). The mean pairwise genetic distances (d) between Indo-Malay samples and the Australian region ranged from 0.0264 to 0.0462, while the pairwise distances within the Indo-Malay region ranged from 0.0008 to 0.0229 (see Additional le 1). Specimens from the west coast of Sumatra were more genetically similar to Muncar and Lombok (southern Indonesia), although located closer geographically to the Indo-Malay region (see Fig. 1).
The median-joining haplotype network showed ne-scale groupings of the 13 haplotypes, with two main haplogroups between samples from the Indo-Malay region and those from the Indian Ocean and Australian regions. Within the Indo-Malay grouping, samples from BIN, WKL, and PAH formed a large haplotype, separated from WJV by 0.2% of sequence divergence. While PHU, SUL, MUN and LMB, formed a sub-cluster of haplotypes. Samples from the west coast of Sumatra (WSA and WSS) formed a distinct haplotype, while Australia (WAU and SEQ) and Papua New Guinea (PNG) formed an Australian region cluster (Fig. 2b). The largest haplotype divergence (1.7-2.8%) occurred between the Indo-Malay samples and those from the Australian region. Samples from the Indian Ocean region (WSA, WSS, MUN and LMB) formed distinctive haplotypes separated from the Australian region by 0.5-1.4%. In contrast, samples from South Sulawesi (SUL) formed a different sub-cluster of the Indo-Malay haplotypes by 0.6-0.8%.

SNP Genotyping
A total of 82,994 SNPs were detected from 148 individuals across 16 sampling locations. Filtering criteria resulted in the subsequent reduction of SNP loci to 6,099 SNP loci (see Additional le 2) that was used for the population structure and phylogenetic analyses. All individuals clearly grouped into their geographical regions. A PCoA based on axis 1 (38.5% of the variation among individuals) and axis 2 (19.5%) showed three main clusters for the C. punctatum-complex ( Fig. 3a): the Indo-Malay region, the Australian region and the Indian Ocean region (LMB, WSA, WSS). However, axis 3 (9.6%) showed a clear separation of LMB from the WSA-WSS group (Fig. 3b). Individuals from South Sulawesi (SUL) showed some separation from the Indo-Malay region. Separation of individuals from SEQ, and those from PNG and WAU was evident within the Australian region ( Fig. 3a, b).

Phylogenetic Inferences
There were no xed differences among locations in the Indo-Malay region except for South Sulawesi (SUL), which showed minimal xed differences of 1-3%. The west coast of Sumatra (WSA) had xed allele differences of 10-22% with all other sampling locations. A similar pattern was observed in individuals from Lombok (LMB), which differed from all other locations by 10-21%. The Australian region differed from all other locations by 6-23% xed differences but showed low xed differences within the region (1-2% xed differences) (see Additional le 4).
Two amalgamation steps were undertaken based on xed allele differences. Firstly, only locations with no xed differences resulted in all individuals from the Indo-Malay region (BIN, PAH, PER, WKL, WJV, EJV, EKL, SAB, and PHU) being pooled into one group (called PUN). Few but no signi cant differences in each comparison (p > 0.05) combined South Sulawesi (SUL) with the Indo-Malay group (PUN) but Lombok (LMB) and the west coast of Sumatra (WSA, WSS) remained separated. Individuals from PNG and WAU grouped together with SEQ into the Australian region group. The nal result was four group entities with absolute xed-differences, and we designated them as four distinct operational taxonomic units (OTUs) (see Additional le 5). The distance tree generated from the Fitch-Margoliash distance analysis combined the west coast of Sumatra (OTU 2 ) and Lombok (OTU 3 ) into a single clade (Fig. 5).
The maximum-likelihood from RAXML analysis of 148 individuals with 381,018 bp of concatenated SNP sequence revealed three distinct clades. Most individuals from the Indo-Malay region, with the exception of SUL, formed a clade, while LMB, WSA and WSS formed a separate clade (Fig. 6a). Although individuals from the Australian region occur in the same clade, individuals from SEQ were clearly separated from WAU and PNG in a longer branch (Fig. 6b).

Species delimitation
The Bayesian coalescent-based species tree derived from the SNAPP analysis divided C. punctatum into four distinct clades corresponding to the OTUs. Relationships among locations within clades remained consistent in all topologies ( Fig. 7 and Additional le 6). The Bayes factor delimitation (BFD*) analysis supported the hypothesis to delimit C. punctatum into four putative species based on OTUs as the topranked model (H-5) with the estimate of the highest value of the marginal likelihood with a decisive Bayes factor (2log e BF > 10) against the initial model (H-1) (28) ( Table 1). Table 1 The Bayes factor delimitation (BFD*) analysis for the ve species delimitation models of Chiloscyllium punctatum with the rank based on the marginal likelihood estimates (MLE). Two models show positive Bayes factor (BF) against the initial model H- Maximum TL measurements from each location supported OTU 1 as different from the other OTUs (see Table 2) with all individuals from OTU 1 reaching their maximum TL at a smaller size (< 110 cm TL) compared with those from OTU 2, OTU 3 and OTU 4 (> 110 cm TL). Based on available references, C.
punctatum attains maturity at ≤ 100 cm TL, with individuals from Indo-Malay region (corresponding to OTU 1 ) mature at around 67-70 cm TL (29,30), while those from Australia and PNG (corresponding to OTU 4 ) mature at about 84-87 cm TL (24,31). No information on size at maturity is available for the regions that fall within OTU 2 and OTU 3 . Nevertheless, those biological features provide an indication of differences among OTUs that support the genetic ndings.  7 Last and Stevens (24). All other samples collected as part of this study (2018-2019). n = number of tissue samples analysed for SNPs and mtDNA (in parentheses); TL = total length; Max TL = recorded maximum TL from each location.

Species-Population Boundary
Both mitochondrial DNA and genome-wide SNP data analysis support a species complex in the brownbanded bamboo shark. By sampling comprehensively across a part of its range, we have revealed greater genetic variation and potentially more species than suggested by a previous study (18) that included fewer sampling locations. Obtaining representative samples across a species' distributional area can be challenging for elasmobranchs. There are potential di culties in obtaining a comprehensive sample set for a species, due to; rarity in nature or in markets; di cult to catch; nancial cost; large specimen sizes; wide geographical spread and species protections, which limits access to samples or requires special permits (32). Still, many new species of elasmobranch have been described based on limited geographic sampling over the last decade (e.g. (33)(34)(35)).
Our analyses on the genetic data of Chiloscyllium punctatum support the hypothesis of Naylor et al. (18) that a potential cryptic species occurs in Australia. However, by including a more comprehensive suite of sampling sites that captured more genetic variation, we were able to identify four OTUs within C.
punctatum. Yet, as the OTUs are allopatric, it is di cult to discern the population-species boundary on genetic data alone. Depending on the nature of the supporting information, C. punctatum species complex may consist of three or four putative species. Based on the genetic data alone, the four OTUs can be interpreted as four distinct species supported by both phylogenetic and species delineation (BFD*) analyses.
The type locality for the original description of C. punctatum was described from a geographic location associated with OTU 1 in the waters off Jakarta, Indonesia (36). Specimens from Jakarta were clearly grouped in OTU 1 , which comprises individuals that inhabit shallow waters (20-50 m) within the Indo-Malay Archipelago (22,23). Individuals from OTU 1 are not only genetically separated to those in the other OTUs, but the distinction is also supported by biological, ecological and geographical data. Individuals from OTU 1 have a relatively smaller body size with the maximum TL at around 1 m, compared with other OTUs that have a maximum TL more than 130 cm (see Table 2). That Indo-Malay individuals mature at a smaller size than do those from the Australasia further supports species delimitation (24, 27, 30,31). However, no information on size at maturity is available for individuals from the Indian Ocean region (OTU 2 and OTU 3 ) and given this biological feature is likely important for species delimitation in this group, further studies are needed. Ecologically, we observed that C. punctatum from the Indo-Malay region is commonly found on soft-bottom habitats (e.g. sand at, muddy substrate). Geographically, all locations that were linked by shallow waters within the Sunda Shelf region (from Phuket, Thailand to East Kalimantan, Indonesia) showed strong connectivity except for South Sulawesi, which is geographically the furthest east location included in OTU 1 . South Sulawesi is also separated from the Indo-Malay region by a deep trench that could restrict gene ow.
Individuals from the Indian Ocean and Australian regions may be grouped into a putative species based on similarities of their biological characters, such as the larger maximum body size than for the Indo-Malay region, or may be separated into two distinct species based on ecological and geographical perspectives. However, the BFD* analysis lent support to two separate groups. Therefore The geographical barrier in the Lombok Strait is an important feature for speciation in particular groups such as the separation of two blue-spotted maskray species, Neotrygon caeruleopunctata in the west and N. australie in the east (37,38). In contrast, this barrier is just resulting in population structure for other species such as the Indonesian wobbegong, Orectolobus leptolineatus (33) and Indonesian speckled catshark, Halaelurus maculosus (39) that occur along the eastern Indian Ocean. Additional sampling at locations along the south coast of Java and the west coast of Sumatra is required to bridge the geographic gap that exists in our data between Sumatra and Lombok to further resolve whether the separation demonstrated here is representative of speciation or population differences. Nevertheless, in terms of habitat preference, C. punctatum from along the west coast of Sumatra, south of Java to Nusa Tenggara in eastern Indonesia (including Lombok), inhabit similar coral and rocky reef habitats.
The bamboo shark species from the Australian region (OTU 4 ) can be found in various habitats, including coral reefs, seagrass beds, mangroves and estuaries (24). The tectonic plate boundary between the Sahul Shelf (South PNG, Western Australia, and South East Queensland) and the Sunda Shelf (Greater and Lesser Sunda Islands) is considered a major barrier causing a genetic break that separates populations and potential species within this complex. Examples of the in uence of this barrier on species level differences can be found in blackspot sharks and banded eagle rays (40,41). Even though species in the Australian region formed a distinct clade, there are some xed differences between specimens from South East Queensland and those from South PNG and Western Australia. Nevertheless, there is no physical barrier, such as deep water or strong currents that might be a driver of genetic segregation between these locations. Therefore, we suggest that locations in OTU 4 represent a single species with population structure, and the differences within the clade may re ect a variation by distance effect (42). This should be further tested by intermittent sampling locations in northern Australian waters.
Our study also revealed that geographic distance does not directly correspond to genetic distance. For instance, Aceh (OTU 2 ) is separated from Phuket (OTU 1 ) by only about 470 km while some sites within the OTU 1 cluster (such as between Phuket and West Java) are separated by more than three times the distance. Similarly, Muncar and East Java are separated by only 300 km and connected by the narrow and shallow waters of the Bali Strait, yet showed strong genetic separation. A possible reason for that separation is habitat type. In addition to the role that deep water and tectonic plate barriers play, habitat preferences are likely to be important in either structuring populations or delimiting species (43,44). Based on habitat preferences, the four OTUs can be further delimited into three groups with OTU 2 and OTU 3 combined. The BFD* analysis also lent high support to this hypothesis.
The decision to delineate species complexes may vary among taxa. In some taxa, genetic differences may not be reciprocal with either morphology (phenotypic plasticity), biology (reproductive traits), or ecological characteristics (habitat preferences) (45). In elasmobranchs, genetics has been used to can delineate cryptic species in a species complexes in the genera Carcharhinus, Aetobatus and Neotrygon (38,40,46,47), as well as coalesce two genera into one genus due to genetic indistinctiveness, such as in Mobula (48).
For some taxonomic groups such as birds and some mammals, the term 'subspecies' is used to de ne a population or group of populations that are distinctive yet insu ciently different to constitute a separate species based on subtleties in appearance and/or in genetic makeup (49)(50)(51). The term 'subspecies' is also used to differentiate a species complex based on ecological speciation for populations without multigene discontinuity (52). Moreover, 'subspecies' is applied to geographically isolated populations driven by biodiversity and conservation purposes, such as in some freshwater shes (50). However, the use of subspecies in some taxa is not preferable due to confusion with the population term (53,54), and has been rarely used in the past few decades for marine shes (55). For elasmobranchs, this term has only been applied to few taxa, e.g. in catsharks (56)(57)(58)(59), dog shes (60), smooth-hounded sharks (61), hammerheads (62), eagle rays (63), and skates (64). Some of those subspecies remain valid such as the smooth-hounded shark (Mustelus canis canis and M. c. insularis) and for several species of skates (from Genus Raja and Leucoraja), while others were considered junior synonyms or have been elevated into distinct species (65). Therefore, the use of subspecies for the bamboo shark OTUs is plausible if they cannot be de nitively classi ed at either the species or population level.

Conservation implications of the species delimitation
Splitting cryptic species in a complex into one or more distinct species may provide advantages for the species, not just for formal scienti c recognition, but also to assess conservation risk (3,4,6,66). Thus, from a conservation perspective, separating C. punctatum into two or three species is desirable as differences in biological features, spatial distribution, habitat occupancy, and type of shery that operates in an area could in uence how each species should be managed for sustainability. For instance, due to the intensive trawl operations in the Indo-Malay region (26,67), the species based on OTU 1 would be subjected to higher shing pressure than individuals within OTU 2 and OTU 3 where bottom longline sheries operate due to unsuitable substrates for bottom trawling, and compared to OTU 4 where they are neither targeted nor caught as bycatch. This situation may lead to differences in threat pro les that necessitate revision among the OTUs of their conservation status.
In terms of sheries management, each OTU with distinct sheries characteristics can be treated separately. For instance, limits on minimum size or permitted shing gear may differ between OTUs due to the nature of the shery. For countries that implement ecosystem or species-based conservation management such as in Indonesia (referring to Regulation of the Minister of Forestry P.57/Menhut-II/2008 and Regulation of the Minister of Marine Affairs and Fisheries 3/2010), regulating two different options of sheries management for one species is challenging, compared with countries that implement conservation strategies on a population-basis. Therefore, splitting C. punctatum into at least two or three different species for management purposes is appropriate for countries that apply species-based management, as long as diagnostic morphological characters are available. Especially when they are marketed in the same place with a lack of traceability, as occurs in Indonesia.
In contrast, splitting a cryptic species complex into several species based purely on genetics can cause problems if there are no strong supportive diagnostic characters to differentiate them in the eld. An example is in the Neotrygon kuhlii species complex (38,68) where several species are sympatric yet cannot be diagnosed morphologically, which complicates management. Policymakers may nd it di cult to implement conservation and management actions, especially if the species are sympatric, inhabit similar ecotypes, or occur in one shing region. Without the ability to differentiate among the members of a complex, studies on their biology, ecology, or population stock for management purposes are rendered problematic due to the likelihood of misidenti cation and possibly overlapping information (37). This is particularly challenging in countries or regions where access to genetic analysis is still limited and costly (69)(70)(71).

Conclusions
Genetics is a powerful tool for detecting and differentiating cryptic species complexes; however, genetic differentiation that lacks supporting biological and ecological differentiators, may not be usefully applied to determine on which side of the population-species boundary a taxon falls. Genetic analyses from two populations of a relatively sedentary species that are separated by considerable geographical distance may show large xed differences, supporting them being identi ed as separate species. In comparison, subsequent and su cient sampling across the full range of the species complex would reveal any "intermediate" population, which would prevent bias in determining the e cacy of population versus species delimitation recommendations (72). Therefore, the use of genetics for delimiting species should include comprehensive sampling locations to identify populations that may bridge differences between two or more distinct populations before they are interpreted as separate species.
It has been proposed that guidelines be developed for the application of genomic data to justify taxonomic boundaries before delimiting species from cryptic complexes (12). Determining a new species based on genetic differences should be approached with caution and supported by strong justi cations.
We suggest some essential steps that should be taken before designating a distinct population or OTU derived from a genetic study to be a putative species. Firstly, it is important to ensure that sampling locations are representative across the range with consideration to their basic biological information (e.g., migratory or sedentary, pelagic or demersal), as speciation is likely occurring more often in species with low dispersal ability (73). Secondly, there should be diagnostic morphological character(s) that can be used to distinguish the new species from the former species and other OTUs. Third, there are supporting morphological, habitat or behavioural characteristics that can be used to con dently extirpate the new species from a species complex. Lastly, clear information on geographic range should be included within the description of the new species to avoid ambiguity for users, be they biologists, policymakers or managers, to allow them to conduct further studies and related conservation and shery assessments. In relation to the Chiloscyllium punctatum species complex, we suggest that a taxonomic revision is warranted to enhance conservation and management strategies for the species that presently reside within it.

Sampling and DNA extraction
Tissue samples from Chiloscyllium punctatum were obtained from 17 localities within the Indo-Australian region ( Table 2). Samples from Thailand, Malaysia, and Indonesia (Sumatra, Java, Kalimantan, South Sulawesi, and Lombok) were sourced from local sh markets. Samples were from Moreton Bay, eastern Australia, samples from Papua New Guinea and Western Australia were provided by the national sh collection (CSIRO, Hobart). The size of the largest specimen from each location was noted, either through direct measurement of total length (TL) or derived from other sources of information (Table 2), as previous studies indicate that C. punctatum maximum body size may vary among locations (22,27). Length at maturity data was sourced from published materials for this region (24,(29)(30)(31).
All tissue samples were extracted using the DNA salting out procedure (75). corrected for small sample size (76). A maximum-likelihood phylogenetic tree was constructed using the Tamura-Nei 1993 model with a discrete Gamma distribution (TrN + G) as the best model (+ G parameter = 0.3851). The analyses were run using PAUP* 4.0a163 (77). The tree topology was evaluated using the maximum likelihood heuristic search with ten random stepwise addition sequence replicates. Statistical support for branch nodes was evaluated by bootstrapping across 1000 replicates (78). The phylogenetic tree was constructed using Dendroscope 3 (79). Pairwise distances between sequences were calculated based on the TrN + G model with 1000 bootstrap replicates using Mega version X (80). Haplotype networks were reconstructed to visualise the relationship among sequence data based on the medianjoining network (81) using the PopArt program (82).

Next-generation Sequencing
A total of 148 extracted DNA samples (50-150 ng.µl − 1 ) were genotyped with the DArT-seq approach by Diversity Arrays Technologies (DArT Pty Ltd, Canberra, Australia) using standard protocols as described by Kilian et al. (83). A combination of methylation-sensitive restriction enzymes (PstI and SphI) was used to digest the total DNA and detect SNP variation. Around two million sequences per sample generated from DArT pipelines were fragmented to 69 bp and then combined into clusters by the DArT clustering algorithm. The low-quality base and identical sequences were corrected and analysed by the DArT software (DArTsoft14) to produce candidate SNP markers. SNP markers were identi ed within each cluster by calculating parameters such as average and variance of sequencing depth, the average counts for each SNP allele and the call rate for each sequence across all samples (83,84).

SNPs Filtering and genotyping analyses
The SNP data were ltered using the R package 'dartR' (85). Loci with 100% repeatability (repAvg = 1.0) and without missing values (call rate = 1.0) were retained for subsequent analysis. Monomorphic loci and secondary SNPs within loci were removed, the latter to reduce the possibility of linked fragments. SNP loci with ≤ 1% minor allele frequency were removed.
Principal coordinate analysis (PCoA) was used to visualise the data with the R package 'dartR' (42,43). Pairwise F ST values between populations along with con dence intervals and p-value were calculated based on the method by Weir and Cockerham (87) using the function gl.fst.pop in the 'dartR' package.
Fixed-difference analysis was applied to the ltered SNP data to identify private alleles at a locus that were not shared among sample locations as a robust indication of lack of gene ow (85). Locations represented by fewer than ve individuals were not included in the analysis as they may yield falsepositive results (88). Locations that had non-signi cant xed differences were further amalgamated (pvalue > 0.05). The nal amalgamated groups from the signi cant xed differences were designated as putative operational taxonomic units (OTUs). Analyses were conducted using the functions gl.collapse.recursive and gl.collapse.pval in 'dartR' (85,88).
A phylogeny of the amalgamated xed-differences matrix was constructed using a Fitch-Margoliash distance analysis based on the Euclidean distance. A multispecies coalescent analysis was performed to examine the species delimitation using SNAPP (90) available in BEAST 2 (91). Speci cally, we tested ve different species delimitation hypotheses. The rst three hypotheses tested for two putative species based on different geographic groupings: (1) an Indo Malay-Indian Ocean and Australian region species; (2) Indo-Malay and eastern Indian Ocean vs Australian region; and (3) Indo-Malay and west coast of Sumatra vs Lombok and Australian region. The fourth hypothesis tested for three putative species as informed by habitat preferences (Indo-Malay vs eastern Indian Ocean vs Australian region). Lastly, the fth hypothesis tested for four putative species, based on OTUs identi ed by the xed difference analysis. As SNAPP is computationally demanding, we ran the analysis on one randomly selected individual per location for OTU 1 and two individuals per location for the other OTUs. We repeated the analysis three times with different individuals to ensure reproducibility. The SNAPP datasets were prepared in BEAUTi with each location assigned as separate taxon and run with Markov chain Monte Carlo (MCMC) run for 2 million iterations and 500,000 burn-in. We set a gamma prior distribution with alpha = 3 and beta = 100 for a prior mean of 0.03 theta to re ect 0.3% sequence divergence as the mean divergence among OTUs. Availability of data and materials The datasets supporting the conclusions of this article are available in the Dryad repository, https://doi.org/10.5061/dryad.rbnzs7h89.

Competing interests
The authors declare no competing interests. Authors' Contributions F and CLD conceived and designed the study; F performed data analyses, generated the graphs and wrote the manuscript; CLD, IRT and MBB revised the manuscript. All authors read and approved the nal manuscript.