Frontiers of Biogeography Can sea snakes slither through seascape structure? Comparative phylogeography and population genetics of Hydrophis group sea snakes in Australia and Southeast Asia

Pleistocene sea level changes substantially shaped the biogeography of northern Australia and the Indo-Malayan Archipelago (IMA). For co-distributed species, their phylogeographic and population genetic patterns are expected to be concomitant with geological transformations of the Pleistocene. However, species-specific ecologies and life history traits may also be influential in generating patterns which depart from simple expectations arising from biogeographic features. Thus, comparative population genetic studies, which use taxa that reduces variation in taxonomy and geography, may refine our understanding of how biogeographic elements shape the populations of co-occurring species. Here, we sampled two sea snake species, Hydrophis curtus and H . elegans , throughout their known ranges in the IMA and northern Australia. These sea snakes have similar life history strategies and ecologies as well as overlapping distributions across the Torres Strait, a well-known biogeographic feature. We analysed two mitochondrial DNA (mtDNA) fragments and 10 microsatellite loci using traditional population genetic approaches and used Bayesian clustering methods to examine species- specific phylogenetic relationships, genetic diversities, and population genetic structures. For both species, we found a consistent lack of significant genetic variation among sampling sites across the Gulf of Carpentaria (GOC) and the Great Barrier Reef (GBR). Similarly, Bayesian clustering showed no to weak genetic partitioning across the historical Torres Strait land bridge. Both species sampled in Australia displayed population expansion signatures in tests using mtDNA and microsatellite markers. We conclude that the phylogeographic and population genetic patterns of these sea snake species do not align with the Torres Strait land bridge. This lack of population genetic structure departs from previous findings on Aipysurus sea snakes and may be linked to the association of Hydrophis species to soft sediment habitats typically found across northern Australia. These divergent patterns between the sea snake groups present the importance of considering taxon-specific attributes in formulating conservation strategies.

• Preference of Hydrophis species for continuous soft sediment habitat types may have allowed genetic connectivity among sampling sites across the northern Australia.
• Population expansion signatures for both species around the Torres Strait region may be linked to extension of suitable habitat after the last glaciation peak.
• Our results contrast with the phylogeographic and population genetic patterns known for co-occurring Aipysurus sea snakes.

Abstract
Pleistocene sea level changes substantially shaped the biogeography of northern Australia and the Indo-Malayan Archipelago (IMA). For co-distributed species, their phylogeographic and population genetic patterns are expected to be concomitant with geological transformations of the Pleistocene. However, species-specific ecologies and life history traits may also be influential in generating patterns which depart from simple expectations arising from biogeographic features. Thus, comparative population genetic studies, which use taxa that reduces variation in taxonomy and geography, may refine our understanding of how biogeographic elements shape the populations of co-occurring species. Here, we sampled two sea snake species, Hydrophis curtus and H. elegans, throughout their known ranges in the IMA and northern Australia. These sea snakes have similar life history strategies and ecologies as well as overlapping distributions across the Torres Strait, a well-known biogeographic feature. We analysed two mitochondrial DNA (mtDNA) fragments and 10 microsatellite loci using traditional population genetic approaches and used Bayesian clustering methods to examine speciesspecific phylogenetic relationships, genetic diversities, and population genetic structures. For both species, we found a consistent lack of significant genetic variation among sampling sites across the Gulf of Carpentaria (GOC) and the Great Barrier Reef (GBR). Similarly, Bayesian clustering showed no to weak genetic partitioning across the historical Torres Strait land bridge. Both species sampled in Australia displayed population expansion signatures in tests using mtDNA and microsatellite markers. We conclude that the phylogeographic and population genetic patterns of these sea snake species do not align with the Torres Strait land bridge. This lack of population genetic structure departs from previous findings on Aipysurus sea snakes and may be linked to the association of Hydrophis species to soft sediment habitats typically found across northern Australia. These divergent patterns between the sea snake groups present the importance of considering taxon-specific attributes in formulating conservation strategies.

Introduction
The Indo-Malayan Archipelago (IMA), together with tropical Australia, is renowned for its unparalleled levels of marine biological diversity (Tittensor et al. 2010) and a dynamic geological history that provided significant ecological opportunities for biological diversification (Voris 2000). In particular, the IMA has experienced sea level fluctuations, so that shallowwater habitats were cyclically formed and lost thereby influencing the contemporary distributions of extant marine organisms (Barber et al. 2006, Mirams et al. 2011. More importantly, the reoccurring emergence and subsidence of geological features across this region during the Pleistocene epoch appears to have shaped the population genetic patterns of many extant codistributed species (Hewitt 2000).
Throughout the IMA and northern Australia, a shared biogeographic history suggests that sympatric marine species should have congruent genetic signatures (Avise 2000) and this pattern has been observed for several but not all species (reviewed in Mirams et al. 2011). The Torres Strait is one of the bestknown transient Pleistocene landmasses (Fig. 1). This land bridge connected northern Australia and Papua New Guinea at low sea level stands, thus repeatedly isolating and reconnecting fauna in the Indian and Pacific Oceans. The most recent occurrence of the Torres Strait land bridge ended ~7,000 years ago (Voris 2000). The present-day Gulf of Carpentaria (GOC) is a feature of high sea stands, and should enable dispersal corridors in the form of shallow-water habitats including through the Torres Strait. Despite these oceanographic reconfigurations being experienced among co-distributed species in this region, many marine species spanning the Torres Strait have been shown to display divergent population genetic patterns. For some marine taxa, populations west (GOC) and east of the Torres Strait (Great Barrier Reef; GBR) formed distinct genetic clusters, while for other species, populations from the same locations were genetically similar (Mirams et al. 2011). Such contrasting results may arise from the inherent stochasticity of microevolutionary processes, but variations in life history strategies among co-distributed marine taxa can also contribute to discordant phylogeographic patterns (Duminil et al. 2007, Harvey et al. 2017.
Analyses that sample taxa in a manner that minimises taxonomic and/or spatial variation may uncover concordance between biogeographic features and spatial genetic processes. Such comparative population genetic studies ideally employ a sympatric sister-species approach and apply the null hypothesis that sympatric sister taxa -taxa with similar life history traits, ecologies, and experienced abiotic environment -will have similar phylogeographic and spatial population genetic structures (Burridge et al. 2008, Dawson 2012. Additionally, marine taxa that give birth to live young (i.e., no pelagic larval phase) Figure 1. Map of sampling locations in this study. Triangles indicate sites where both species were sampled. The inset map shows the distribution for each species: Hydrophis curtus (purple), H. elegans (green); and the region of overlap (orange). The 120-m bathymetric contour is denoted by the dashed line (brown), which corresponds to the low sea level stand during the late Pleistocene (~17,000 years ago) (adapted from Mirams et al. 2011 may be assumed to disperse poorly and sustain signals of historical biogeographic events via congruent population genetic patterns among similar species (Lukoschek 2018a). Therefore, co-distributed, closely related species with direct development are most likely to show congruent phylogeographic patterns coincident with biogeographic features. True sea snakes (Elapidae; Hydrophiinae) are a young but diverse evolutionary radiation possessing many characteristics that make them ideal for evaluating population genetic signatures shaped by Pleistocene vicariance events (Gaither and Rocha 2013). Hydrophiine sea snakes are the most diverse extant marine reptile group which diverged in the last ~6 million years (Sanders et al. 2013), comprising two evolutionary lineages: a smaller Aipysurus group (11 nominal species in 2 genera) and a larger Hydrophis group (46 nominal species in 2 genera) (Sanders et al. 2013, Lukoschek 2018b. The IMA and tropical Australia have the highest sea snake diversity, with many species having overlapping ranges that traverse important biogeographic features including the Thai-Malay Peninsula, Timor Trench, and the Torres Strait (Ukuwela et al. 2017). All sea snake species are fully aquatic and viviparous (thus lacking a pelagic larval stage). Moreover, the most species-rich genus Hydrophis is estimated to have diversified more recently (~1.5 million years ago) coinciding with the late stages of the Pleistocene sea level changes (Voris 2000, Sanders et al. 2013. Previous comparative phylogeographic and population genetic studies on sea snakes of the Aipysurus group recovered strong population genetic structures coincident with biogeographic features across northern Australia including the Torres Strait (Lukoschek et al. 2007, Lukoschek 2018a. While previous assessment of Hydrophis curtus populations from Southeast Asia and Australia revealed discrete phylogeographic breaks congruent with major geographic features (i.e., Thai-Malay Peninsula and Timor Trench) (Ukuwela et al. 2014), the influence of the Torres Strait on the population genetic patterns of Hydrophis species has not been explicitly documented.
Here, we examine two species from the genus Hydrophis: Hydrophis curtus and H. elegans, to investigate whether these co-occurring species demonstrate congruent genetic diversities and population genetic structures that align with geological events and oceanographic features. Hydrophis curtus (Shaw 1802) is distributed among shallow marine habitats (down to a depth of 55 m) from the Arabian Gulf extending throughout Asia and into Australia and New Caledonia (Rasmussen et al. 2010). Hydrophis elegans (Gray 1842) has a more restricted distribution and occurs in the coastal waters of northern Australia and southern New Guinea, diving to a depth of 110 m (Heatwole 1999, Milton et al. 2009). The ranges of these species overlap around the Torres Strait region of northern Australia and in the GBR (Fig. 1). We take advantage of their similarities in life history strategy (viviparity), ecology (fully aquatic, shallow-water marine reptiles), and experienced abiotic conditions and processes (overlapping distributions) to limit taxonomic and spatial variation in our comparisons.
We hypothesised that these closely related species will exhibit similar levels and distributions of genetic diversities and display congruent patterns of population genetic structure where their distributions overlap. We expected to recover strong population genetic partitioning particularly for populations flanking the western and eastern regions of the Torres Strait, consistent with the population genetic structures observed in the Aipysurus group of hydrophiine sea snakes with same geographic distributions (Lukoschek et al. 2007, Lukoschek 2018a. We combined data from two mitochondrial DNA fragments and 10 microsatellite loci obtained for both species sampled throughout their known ranges, representing the most comprehensive sampling regime for a comparative phylogeography and population genetics study on Hydrophis species to date. Most tissue samples were obtained from trawl fisheries bycatch. In Australia, bycatch samples were collected primarily as part of scientific studies and the exact location of each snake obtained was recorded. Samples from SE Asia were acquired primarily from commercial trawlers that operated out of Phuket (THW) and Song Kla (THE) and the exact location where each snake was based on approximate provenance. Sampling sites are shown in Fig. 1 and a summary of sample details is presented in Table 1 and Appendix S1.

DNA extraction, mitochondrial DNA sequencing, and microsatellite genotyping
Most tissue specimens were stored at room temperature in 70% ethanol while THW samples were obtained from dried snakes. DNA extraction, PCR amplification, and sequencing of mitochondrial DNA (mtDNA) followed protocols detailed in Lukoschek and Keogh (2006); Lukoschek et al. (2007); Lukoschek et al. (2008). Two mtDNA fragments were targeted for analyses: ATPase subunits 8 and 6 (~850 bp) and ND4 (~700 bp) plus adjacent 3' tRNA-His+tRNA-Ser (~120 bp). Primers used for PCR amplification and sequencing were as follows: ND4 (L) Ten polymorphic nuclear microsatellites developed for H. elegans, previously shown to amplify in 13 congeneric species were used to genotype all snakes using the primers and protocols described in Lukoschek and Avise (2011). Alleles were sized using a ROX labelled GS500 internal standard and scored using GeneMapper v4.0 Software (Applied Biosystems).
Sequences were edited in Sequencher v4.5 and were aligned with Clustal W2 (https://www.ebi.ac.uk/ Tools/msa/clustalw2/) using default parameters and were visually checked for sequencing errors and ambiguous sites. Following alignment, all coding region sequences were translated into amino acid sequences in MacClade v4.06 (Maddison and Maddison 2003) using the vertebrate mitochondrial genetic code. No premature stop codons were observed, thus affording confidence that mtDNA sequences were not nuclear mtDNA pseudogenes. Gene genealogies for H. curtus and H. elegans were estimated using statistical parsimony implemented in TCS v1.23 (Clement et al. 1999) with parsimony criterion set to 95% and gaps in the sequence treated as fifth state (Ogden and Rosenberg 2007). All sampled individuals were included and the geographical locations of identified haplotypes were mapped on to the resulting network. Haplotype networks were visualised using POPART v1.7 (Leigh et al. 2015). Additionally, phylogenetic relationships among all sampled haplotypes were estimated using maximum parsimony (MP) and maximum likelihood (ML) analyses conducted respectively in PAUP* (Swofford 2000) and RaxML v8.2.12 (Stamatakis 2014) on the CIPRES Science Gateway (Miller et al. 2010). Three closely related species (H. major, H. ornatus, H. platurus) were used as outgroups to root the trees. MP analyses were performed using heuristic searches with 1,000 random stepwise addition replicates and tree-bisectionreconnection (TBR) branch swapping and all sites weighted equally. ML analyses using 1,000 bootstrap iterations were conducted with default parameters except for the model of molecular evolution, which was set to "GTRGAMMA".
Three molecular variance components (among groups, among sites within groups, within sites) and population genetic structure were estimated using a hierarchical analysis of molecular variance (AMOVA) framework. AMOVA analyses were performed based on both raw haplotype frequencies and sequence divergences among haplotypes based on the Tamura-Nei model (Φ ST ) to test different a priori hypotheses (Tables S1 and S2). Global and pairwise F ST and Φ ST values were calculated among locations with five or more samples. Analyses were performed in ARLEQUIN using 20,000 random permutations. P-values were adjusted with sequential Bonferroni corrections for multiple comparisons (Rice 1989). In addition, Tajima's D tests (Tajima 1989) were performed to detect signals of population expansion. Sites were grouped based on the resulting haplotype network prior to analyses. Fu's F S values (Fu 1996, Fu 1997 were also calculated as an alternative test for population expansion with both tests conducted in ARLEQUIN using 20,000 random permutations.

Analyses of microsatellite data
Summary statistics for genetic diversity Each site was estimated for three indices of genetic diversity: number of alleles (N A ), observed heterozygosity (H O ), and expected heterozygosity (H E ). Deviations from Hardy-Weinberg equilibrium were examined using exact test with 1,000,000 forecasted Markov chain length followed by 100,000 dememorisation steps. All pairs of loci were also tested for possible linkage disequilibrium. All analyses were performed using ARLEQUIN. Genetic diversity across sites and within regions were assessed using allelic richness (A r ) rarefacted based on the minimum sample size of diploid individuals for each dataset (3 for H. curtus, 1 for H. elegans) in FSTAT v2.9.4 (Goudet 2001).

Population genetic structure
Population genetic structure for each species was estimated using global and pairwise F ST values. Global F ST values were calculated using the GenePop v4.7.5 (Rousset 2008) implementation in R (R Core Team 2020). Pairwise F ST (= ϴ; Weir and Cockerham 1984) was calculated with 20,000 permutations in ARLEQUIN. To account for allelic variability, global and pairwise G' ST (Hedrick 2005) measures were estimated using the 'MMOD' package (Winter 2012) in R. In addition, principal components analysis (PCA) of allele frequencies was conducted to examine patterns of spatial genetic structure among sampling sites. Missing data for each locus were imputed with mean genotype frequencies using the R package, 'adegenet' (Jombart 2008). Furthermore, multidimensional scaling (MDS) was performed in R to visualise the genetic relationships among populations based on pairwise G' ST values.

Analysis of molecular variance (AMOVA)
Sampling sites were grouped to test a priori hypotheses with AMOVA and were evaluated at three hierarchical variance components: among groups, among sites within groups, and within sites (Tables S8 and S9). All analyses were performed in ARLEQUIN using 20,000 permutations.

Isolation-by-distance
We used distance-based Redundancy Analyses (dbRDA) to evaluate the effect of geographic distance on genetic variation. We expect isolation-by-distance (IBD) to be present globally for both species based on the known restricted movement and home ranges of hydrophiine sea snakes (Burns andHeatwole 1998, Lukoschek andShine 2012). Additionally, we performed dbRDA limited to sites in the Torres Strait vicinity to assess the impact of the historical Torres Strait land bridge on the population genetic connectivity of each sea snake species. Pairwise geographic distances were measured using the 'marmap' package (Pante and Simon-Bouhet 2013) in R. Calculations of least-cost path distances were done separately for each species restricting to known species-specific maximum occurrence depth limits: 55 m for H. curtus (Lukoschek et al., 2018) and 110 m for H. elegans (Milton et al., 2010); in order to capture the benthic habitat associations of sea snakes (with exception of the pelagic H. platurus). For H. curtus, maximum occurrence depth limit was adjusted to 1,275 m for the calculation of distances among populations separated by the Timor Trench (between SE Asia and Australia) (E. Pante and B. Simon-Bouhet, personal communication, July 17, 2020). package (Dixon 2003) in R. Statistical significance was determined using 10,000 permutations.

Genetic clustering
Individual assignment tests were performed for each species to examine spatial genetic structure and to estimate the optimal number of genetic clusters. These tests were conducted in STRUCTURE v2.3.4 using 30 independent runs for each value of K (= 1 to 8) with 500,000 burn-in periods and 100,000 post-burnin iterations under the admixture model, assuming correlated allele frequencies, and no prior information on individual sampling location (LOCPRIOR option disabled) (Pritchard et al. 2000). The most likely value of K for each species was evaluated using the method of Evanno et al. (2005) as implemented in STRUCTURE HARVESTER (Earl and Vonholdt 2012). Individual assignment probability plots were visualised using the R package, 'pophelper' (Francis 2017).
Additionally, genetic clustering analyses were performed in TESS v2.3.1 (Chen et al. 2007). TESS uses a spatially continuous prior based on the geographic coordinates of each individual sampled making it more robust to IBD. For most Australian samples, the exact latitude and longitude was used in the analyses while this information was approximated for samples collected from SE Asia. TESS was run using the CAR admixture model, which assumes spatial autocorrelation of the genomes of individuals in closer geographical proximity compared with those further apart. The strength of this autocorrelation (Ψ) was set to the default value of 0.6. Analyses in TESS were performed using 100 independent runs for each K (= 2 to 8) having a total of 35,000 sweeps (with a burn-in of 10,000 sweeps). The average Deviance Information Criterion (DIC) for the 20 lowest DIC values for each value of K was used to evaluate the most likely number of genetic clusters. TESS results were processed and visualised using 'pophelper.'

Tests for population expansion
We conducted the k intralocus (k-test) and the g interlocus (g-test) variability statistics as described by Reich et al. (1999), using 'kgtests' (Bilgin 2007) implemented in Microsoft Excel®. Identified genetic clusters in STRUCTURE were used to group individuals of each species prior to analyses. Statistical significance for each test was calculated based on a binomial distribution with the probability of a positive k set as 0.515 (k-test) and on the five-percentile cut-off value presented in Reich et al. (1999).

Genetic diversity
The final H. curtus mtDNA alignment comprised 1,594 bp (ATPase8 -152 bp; ATPase6 -676 bp; ND4 -625 bp; and tRNAs -141 bp) with 96 polymorphic sites that defined 63 unique haplotypes among 156 individuals. Seventy-one variable sites were transitions, eighteen were transversions, and seven were indels. Pairwise differences between haplotypes ranged from 1 to 32 base pairs. Overall nucleotide diversity was low (π = 0.68% ± 0.34 SE), however there was a marked difference between and within regions ( Table 1). The most striking differences were between SE Asia (Regional π = 0.92% ± 0.46 SE) and Australia (Regional π = 0.10% ± 0.07 SE). Within Australia nucleotide diversities were higher in the GOC (π = 0.11% ± 0.07 SE) than the GBR (π = 0.09% ± 0.06 SE). Overall haplotype diversity was high (h = 0.93) with SE Asia having higher haplotype diversity (h = 0.96) than Australia (h = 0.82) (  Table 2). The variable sites comprised 72 transitions and 13 transversions. Pairwise differences between haplotypes ranged from 1 to 38 base pairs. Overall nucleotide diversity was low (π = 0.53% ± 0.27 SE). Local and regional nucleotide diversity estimates were lowest in WA (π = 0.10% ± 0.09 SE) and highest in the GOC (π = 0.62% ± 0.32 SE). In contrast, overall haplotype diversity was high (h = 0.96), ranging from 0.83 in the GBR and WA to 0.96 in the GOC ( Table 2). GBR and the GOC shared 30 haplotypes while the four individuals from WA shared 3 unique haplotypes ( Table 2). It is important to note, however, that for both species, some individuals lack data from sites that presently appear to distinguish haplotypes from each region. A detailed table of identified haplotypes is provided in Appendix S1.

Gene genealogies and population genetic structure
Three genetic clusters were consistently recovered for H. curtus by statistical parsimony (Fig. 2a) and MP & ML analyses (see Figs. S1 to S4). These genetic clusters corresponded to Western Thailand (THW), Eastern Thailand and Indonesia (THE+INDO), and Australia (GOC and GBR). In addition, THE+INDO and Australian haplotypes formed sister clades in the MP and ML phylogenetic trees although this For H. elegans, WA haplotypes were revealed to be unique from the rest of the sampling locations (GOC and GBR), which showed no distinct geographical pattern, particularly in relation to the Torres Strait ( Fig. 2b; Figs. S2, S4, S5). AMOVA analyses taking into account sequence divergence reflected these clusters with the highest percentage of variation occurring in Indonesia and Thailand (86.62%, P = 0.006; TT+MP in Table S1). In addition, there was no significant variation among groups on either side of the purported Torres Strait barrier (F CT = -0.024, P = 0.90; Φ CT = -0.016, P = 0.70) (Table S1). Global F ST was far higher when accounting for evolutionary distances among haplotypes (Φ ST = 0.875, P < 0.001) than for raw haplotype frequencies (F ST = 0.159, P < 0.001). Highest Φ ST values occurred among pairwise comparisons between sites in SE Asia and Australia (0.795-0.911) (Fig. 3a) (see Table S3 for pairwise F ST and Φ ST values).
For H. elegans, there was no significant population subdivision observed relating to the Torres Strait land bridge (F CT = 0.016, P = 0.10; Φ CT = 0.088, P = 0.20) (Table S2), which is consistent with the lack of spatial genetic structure for this species' haplotype network (Fig. 2b). Relative to H. curtus, the global fixation index for H. elegans was lower when considering sequence divergence (Φ ST = 0.260, P < 0.001) and for haplotypes frequencies (F ST = 0.082, P < 0.001). Estimates of Φ ST were greatest and statistically significant among pairwise comparisons between GROT and all other sites (0.217-0.415) (Fig. 3c and Table S3).
Tajima's D and Fu's F S reported congruent signals of population expansion for H. curtus (Table 3) (Table 3).   There was no clear spatial pattern in haplotype distribution for H. elegans. Note: Some individuals lack data from sites that presently appear to distinguish haplotypes from each region. Table 3. Results of tests for population expansion for Hydrophis curtus based on mitochondrial data: Tajima's D (D) and Fu's F S (F S ) and corresponding significance levels (P D and P FS ); and microsatellite data: k-Test and g-Test (Reich, Feldman, and Goldstein 1999), statistical significance indicated by boldface.

Summary of genetic diversity
All ten microsatellite DNA loci were polymorphic for both H. curtus and H. elegans (Tables S4 and S5) (see also Lukoschek and Avise (2011)). For H. curtus, the number of alleles (N A ) per locus within sites ranged from 1 to 13, average expected heterozygosity by site (H E ) ranged from 0.461 to 0.737, and the highest mean allelic richness across all sites (A r ) was observed in THE (A r = 3.729) ( Table 1). For H. elegans, N A per locus within sites ranged from 1 to 13 and H E ranged from 0.633 to 0.739. The highest A r was observed at WA (A r = 1.738) following the exclusion of GBRS from the analysis since this site had one locus with no genotyped individuals (Table 2). Regional comparisons of mean allelic richness for H. curtus showed the highest level in THE + INDO (A r = 3.504) and for H. elegans in WA (A r = 1.738). For regions where both species overlap, highest A r was observed in GOC for both species (Tables 1 and 2). Most populations were in Hardy-Weinberg equilibrium for most loci (Tables S4 and  S5) and there was no evidence to suggest linkage disequilibrium among pairs of loci for both species.  (Fig. 3b) as well as the disparity of WA individuals with the rest of GBR and GOC locations for H. elegans (Fig. 3d).
Principal components analyses (PCA) for pairwise F ST values for H. curtus showed that axes 1 and 2 accounted for 24% of the variation among the eight sites (Fig. 4a). Specifically, PC 1 (18%) accounted for the differentiation between Australian and SE Asian populations. In addition, INDO and THE populations clustered closer to each other while maintaining overlap with THW. In contrast, no clear spatial clusters were formed for H. elegans, with PC axes 1 and 2 accounting for only 15.2% of the variation (Fig. 4c). MDS of pairwise G' ST values revealed three distinct clusters for H. curtus (Fig. 4b). Specifically, Australian populations grouped together while INDO and THE formed a distinct cluster, separate from THW among SE Asian populations. For H. elegans, no apparent spatial clustering was observable across sampling locations (Fig. 4d).

Analysis of molecular variance (AMOVA)
AMOVA analyses for both species were consistent with PCA and MDS findings. For H. curtus, the Timor Trench that separates SE Asia from Australia accounted for the highest and significant percentage of variation (14.67%, P = 0.02) (Table S8). In contrast, the Torres Strait was found to be non-significant and contributed for the smallest genetic divergence among groups (1.21%, P = 0.10). Similarly for H. elegans, there was non-significant genetic divergence among sites across the Torres Strait and most of the variation was explained within sites (Table S9).

Isolation-by-distance
Distance-based Redundancy Analyses including all sites revealed strong correlation between geographic distance and genetic variation for both species (H. curtus: R 2 adj = 0.641, p = 0.045; H. elegans: R 2 adj = 0.555, p = 0.047). Analysis for H. curtus confined among sites around the Torres Strait showed lack of significant IBD (R 2 adj = 0.195, p = 0.292). For H. elegans, sites in the GOC and GBR revealed significant moderate correlation between genetic and geographic distances (R 2 adj = 0.535, p = 0.042). However, these regions which flank the Torres Strait did not match the pairwise Specifically, all pairwise comparisons among these sites yielded non-significant F ST values and were only significant using pairwise G' ST for GBRS paired with WEIPA and with MORN. The presence of IBD among sites in the vicinity of the Torres Strait for H. elegans may suggest that an equilibrium between gene flow and drift has been obtained among such sampled sites and may not be directly attributed as a signature of the vicariant effect of the historical land bridge. DbRDA plots are presented in Fig. S6.

Genetic clustering
Individual assignment tests revealed clear spatial genetic structure for H. curtus, with a value of K = 3 showing the best fit (Fig. 5a, bottom). Individuals from SE Asia (THW, THE, and INDO) formed one genetic cluster while individuals from Australia formed two clusters, which corresponded to the GOC and GBR. Alternatively, K = 2 also revealed clusters matching a broader geographic scale with a clear partition between SE Asian and Australian populations (Fig. 5a, top). Clusters formed using TESS were consistent with STRUCTURE bar plots specifically for K = 2 (see Fig. S7c). A separate STRUCTURE analysis which only included H. curtus from Australia revealed two clusters, albeit weak structure, between GBR and GOC individuals (Fig. 5a, middle). In contrast, H. elegans showed no structure for any value of K (Fig. 5b), which is consistent with the PCA plot for this species. TESS analyses have appeared to coerced individuals for this species into three clusters, however such clusters were not consistent across replicated independent runs for K = 3 ( Fig. S7a and S7b).

Tests for population expansion
The intralocus k-test detected significant signatures of population expansion in H. elegans as well as in H. curtus, particularly in the GOC and GBR (Table 3). In contrast, the interlocus g-test did not detect any signature of population expansion in any region or site for both species. However, the power of g-tests in detecting signals of population expansions decreases with increasing variation in mutation rates across loci (King et al. 2000); thus, such non-significance of g-tests results may be explained by the variability in alleles across our microsatellite datasets for each species.

Torres Strait is not a biogeographic barrier for Hydrophis group sea snakes
This study presents the first comparative evaluation of the phylogeography and population genetics of true sea snakes in Hydrophis. Here, we set out to investigate if the similarities in life history strategy (viviparity), ecology (fully aquatic, shallow-water marine reptiles), and experienced similar abiotic conditions and processes (overlapping distributions) between H. curtus and H. elegans yielded phylogeographic and population genetic patterns concordant with known biogeographic features. Particularly, we looked at the regions of the GOC and the GBR in Australia where both species' distributions encompass the Torres Strait, a well-known biogeographic barrier among various marine taxa.
Taken together, our analyses consistently showed a lack of significant phylogeographic and population genetic structure for both species sampled in Australia specifically among sites flanking the Torres Strait. Our findings were unexpected and deviated from previous studies which evaluated the influence of this oceanographic feature on population genetic structure of sea snakes sampled across sites in the region (Lukoschek 2018a). Significant signals of population expansion in both sets of markers coupled with the lack of observable population genetic structure, particularly around the Torres Strait, corresponded with the habitat shifts in this region during the Pleistocene and the characteristic habitat associations of Hydrophis species (see next section).
Our study adds information on previous findings regarding true sea snakes (Elapidae; Hydrophiinae) occurring in Western Australia. Information from mtDNA data suggested the genetic distinctiveness of H. elegans from WA, evidenced by putatively private haplotypes. This finding is consistent with previous studies that found sea snakes in WA to represent a divergent lineage that should be treated as an evolutionary significant unit (ESU) for the purposes of conservation (Sanders et al. 2015, D'Anastasi et al. 2016. For example, Sanders et al. (2015) documented separate breeding populations of Aipysurus foliosquama and A. apraefrontalis which were originally considered to be vagrant individuals. The consistent recovery of unique haplotypes from WA, as in our study and others (Sanders et al. 2015, D'Anastasi et al. 2016, Lukoschek 2018a, warrants further investigation (see Future directions).

Habitat associations may facilitate dispersal in viviparous marine animals
Our results demonstrate that marine animals which lack a pelagic larval phase (i.e., viviparous) and are theoretically associated with reduced dispersal ability (Palumbi 1992), do not always display strong population genetic structure. Here, both Hydrophis species showed weak to no population genetic subdivision between populations on either side of the historical terrestrial barrier, the Torres Strait land bridge. This finding is in stark contrast to population genetic patterns documented for Aipysurus-Emydocephalus sea snakes (Lukoschek 2018a), a closely related group with a similar life history strategy (i.e., viviparity) but with different habitat association. The distinct habitat associations of these two sea snake groups (i.e., Hydrophis and Aipysurus-Emydocephalus clades) may have played a role in providing divergent population genetic patterns across a shared historical barrier.
Sea snakes in the Aipysurus-Emydocephalus group have strong associations with, but not limited to, coral and rocky reef bottoms (Heatwole 1999), while Hydrophis group species typically use relatively continuous soft sediment habitat types (Heatwole 1975, Minton andHeatwole 1975). Indeed, five species from the Aipysurus-Emydocephalus sea snake group (4 Aipysurus spp. And Emydocephalus annulatus) showed phylogeographic and population genetic patterns concomitant to the biogeographic features across northern Australia (Lukoschek 2018a). Such role of habitat preference, and broadly ecological traits, in influencing population genetic patterns has been observed in a number of studies in both terrestrial (Loveless and Hamrick 1984, Burney and Brumfield 2009, Gianoli et al. 2016, Harvey et al. 2017) and marine systems. Similar patterns were found among viviparous marine fish wherein species strongly associated with coral and/or rocky reef bottoms have exhibited high population genetic structure (Bernardi 2000, Planes et al. 2001, Bernardi and Vagelli 2004, Von Der Heyden et al. 2008). In such cases, sandy expanses interspersed among coral/rocky reef habitats have effectively acted as barriers to dispersal. In contrast, it is plausible that viviparous marine animals typically associated with continuous tracts of soft sediment habitat experience less barriers to dispersal and can display higher levels of population genetic connectivity. The signatures of population expansion recovered in both species, particularly in the GOC and GBR, suggest that such demographic expansion may be linked to the extension of suitable habitat types when sea levels rose after the last glacial maximum and may, in turn, have allowed for the population genetic connectivity for Hydrophis sea snakes. Indeed, the North Marine Region of Australia, which includes the Torres Strait and Gulf of Carpentaria, has characteristic soft sediment habitat types with water depths mostly at < 20 m below sea level (Rochester et al. 2008 that track movements, particularly of reproductively mature individuals, will provide information on the utility of these soft sediment habitats as oceanographic corridors to Hydrophis sea snake populations. Collectively, true sea snakes (Elapidae; Hydrophiinae) demonstrate subtle yet important differences in how their populations are structured genetically across their distributions. These diverging population genetic patterns between sea snake lineages and across sea snake species present critical implications that are important in informing and refining actions relevant to their conservation.

Future directions
Sea snakes constitute the most diverse but also one of the most imperilled groups of marine reptiles (Elfes et al. 2013). Increased understanding of their population connectivity and metapopulation dynamics can greatly aid in building the foundation for their immediate conservation and management. Such knowledge also puts forward insights on evolutionary hotspots that advance the basis of establishing ESUs throughout a particular species range. Apart from taxonomic implications (Ukuwela et al. 2014), the consistent recovery of distinct haplotypes and genetic clusters for H. curtus from SE Asia (THW and THE) appears to meet the criteria for their consideration as distinct ESUs (Moritz 1994). This is particularly critical as sea snake populations in SE Asia are poorly studied and are heavily threatened by anthropogenic activities. This scenario also applies to the Australiaendemic, H. elegans, particularly for the populations which occur in WA. We propose that H. elegans sampled from WA are representatives of an overlooked distinct population in this region as supported by the putatively unique mtDNA haplotypes recovered in this study. The western Australian region has previously been documented to harbour overlooked breeding populations of sea snake species (Sanders et al. 2015, D'Anastasi et al. 2016. While our study presents a small sample size, it also emphasises the need for more widespread sampling of Hydrophis species across this region (Udyawer et al. 2018). Such surveys are important in elucidating a more robust description of the phylogenetic relationships and population genetic patterns of sea snake species across their distributions and permit the implementation of more appropriate and effective conservation assessments.

Data Accessibility
GenBank Accession numbers of the samples used in this study are provided in Appendix S1.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1. Maximum parsimony (MP) consensus tree for Hydrophis curtus after 1,000 bootstrap replications with branch labels indicating support values. Figure S2. Maximum parsimony (MP) consensus tree for Hydrophis elegans after 1,000 bootstrap replications with branch labels indicating support values. Figure S3. Maximum likelihood (ML) consensus tree for Hydrophis curtus after 1,000 bootstrap replications with branch labels indicating support values (Tips are aligned for clarity). Figure S4. Maximum likelihood (ML) consensus tree for Hydrophis elegans after 1,000 bootstrap replications with branch labels indicating support values (Tips are aligned for clarity). Figure S5. Haplotype network for Hydrophis elegans using reduced dataset showing the unique set of haplotypes for WA. Figure S6. Distance-based Redundancy Analyses (dbRDA) plots for Hydrophis curtus (a,b) and H. elegans (c,d). Figure S7. TESS individual assignment bar plots for Hydrophis elegans and H. curtus. Table S1. Analysis of Molecular Variance (AMOVA) for mitochondrial DNA sequence data of Hydrophis curtus based on raw haplotype frequencies (F ST ) and sequence divergences among haplotypes (following Tamura-Nei model) (Φ ST ) (Significant values in boldface). Table S2. Analysis of Molecular Variance (AMOVA) for mitochondrial DNA sequence data of Hydrophis elegans based on raw haplotype frequencies (F ST ) and sequence divergences among haplotypes (following Tamura-Nei model) (Φ ST ) (Significant values in boldface). Table S3. Calculated pairwise F ST and Φ ST values for Hydrophis curtus and H. elegans (F ST = below diagonal; Φ ST = above diagonal; significant values in boldface). Table S4. Summary of genetic diversity indices based on ten (10) microsatellite markers for Hydrophis curtus (N A = number of alleles, H O = observed heterozygosity, H E = expected heterozygosity, n = number of samples per sampling site, P HWE = P-value from Hardy-Weinberg equilibrium test, A r = allelic richness per locus and population based on minimum sample size of 3 diploid individuals) Table S5. Summary of genetic diversity indices based on ten (10)