Biogeographic origins of the viviparous sea snake assemblage (Elapidae) of the Indian Ocean

One of the primary goals in biogeography is to understand how different biotas have been assembled in different regions of the world. The presence of the viviparous sea snakes in the Indian Ocean (IO) poses a unique question in this regard due to their evolutionary origins in Australasia (Australia and New Guinea). Here, we examined the origins and patterns of colonization of the IO sea snake assemblage through time-calibrated molecular phylogenies and ancestral area reconstructions. We further evaluated how past and present barriers to dispersal affect genetic diversity of IO sea snakes by examining the population genetic structure of the widespread sea snake, Hydrophis curtus. Our phylogenetic analyses and ancestral area reconstructions strongly indicate that the majority of the IO sea snakes are derived from the Southeast Asian (SEA) sea snake fauna through dispersal and colonization with an in situ radiation (Hydrophis stricticollis-Hydrophis obscurus clade). Further, many species have undergone vicariant speciation events across the Sunda shelf/Indo-Pacific barrier, which formed during the low sea level periods of the Pleistocene. Population genetic analysis of H. curtus revealed a prominent genetic break between populations broadly distributed in the IO and SEA with limited recent gene flow indicating possible cryptic species. These results suggest that compared to the viviparous sea snake stem group that originated 10.6-6.5 million years ago, the IO viviparous sea snakes have a relatively long and complex evolutionary history in the IO and thus have a unique conservation value.


INTRODUCTION
The tropical coastal regions of the IO is home to an assemblage of marine biodiversity that is second in the total number of species only to the Indo-Australian Archipelago (IAA) in the West Pacific (WP) (Tittensor et al., 2010). The coasts of India boast more than 15000 total marine species (Wafar et al., 2011) while the Caribbean biodiversity hotspot in the tropical Atlantic has a comparatively lower marine biodiversity (~12000 total marine species) (Miloslavich et al., 2010). The most species rich regions within the IO are the red sea and regions around Seychelles, Madagascar and Maldives, while endemicity is highest in the Red Sea reefs in the western IO (Mora et al., 2003;DiBattista et al., 2016a).
Several hypotheses/models have been proposed to explain the exceptional diversity found in marine biodiversity hotspots in the tropical regions. The best known is the 'centre of' models (i.e. centre of speciation, accumulation, overlap) that have been frequently invoked to explain the exceptional diversity of the IAA (Bowen et al., 2013). The centre of speciation hypothesis (Ekman, 1953) suggests that marine biodiversity hotspots produce new species at higher rates and act as cradles for speciation. As a result, peripheral areas adjacent to marine biodiversity hotspots are considered to be evolutionary net sinks in which colonizers arrive and evolve to new endemic species, but do not disperse and speciate beyond these regions (Bowen et al., 2013). Conversely, the centre of accumulation hypothesis suggests that new species evolve in peripheral areas and colonise and accumulate in the centre of the range (Ladd, 1960). Finally, the centre of overlap model, predicts that ranges of co-distributed species overlap in central regions resulting in higher species diversity than the peripheries (Woodland, 1983).
All three models have received some empirical support and thus it has been suggested that marine biodiversity hotspots can act collectively as a centre for speciation, accumulation and overlap (Randall, 1998;Barber and Bellwood, 2005;Mironov, 2006). Empirical evidence from reef fish (Carpenter and Springer, 2005;Tornabene et al., 2015) and marine invertebrates (Veron, 1995;Barber and Bellwood, 2005;Lavery et al., 1996;Crandall et al., 2008) indicate that the IAA and the regions within the WP function as a centre of speciation producing new species which move outwards. This is further supported by the observation that the higher proportions of endemic species being present in the furthest corners of the peripheral areas such as the red sea in the IO and the Hawaiian islands in the Eastern Pacific (Mora et al., 2003;Tittensor et al., 2010;DiBattista et al., 2016a). Intriguingly, several studies also provide support to the centre of accumulation model (Hobbs et al., 2009;Eble et al., 2011;Hodge et al., 2012) and centre of overlap model (Hubert et al., 2012;Gaither and Rocha, 2013). These findings have led to the formulation of the Biodiversity feedback model which suggests that marine biodiversity hotspots can act collectively as centres for speciation, accumulation and overlap (Bowen et al., 2013).
It has been generally perceived that the majority of the marine biodiversity of the Indian Ocean has been derived from the colonisers from the IAA or the WP (Ekman, 1953;Veron, 1995;Springer, 2005 Tornabene et al., 2015;DiBattista et al., 2016b). However, there is very little evidence to support this mainly due to the lack of studies in this region. Though there are many theories and empirical evidence to describe the diversity of the IAA, there are few models to describe the patterns and determinants of marine biodiversity in the peripheries of the IAA biodiversity hotspot, the IO and the East Pacific.
Viviparous sea snakes comprise 63 described species that share a terrestrial Australian ancestor only c. 10.6-6.5 million years ago (Ma) (Sanders et al., 2008;Lukoschek et al., 2012;Sanders et al., 2013a;Sanders et al., 2013b). They occupy shallow-marine habitats throughout the tropical and subtropical regions of the Indian and Pacific Oceans. The group comprises two major clades (Aipysurus group and Hydrophis group); the latter includes the majority of viviparous sea snakes. The diversity of the viviparous sea snakes peaks in the Southeast Asian waters of the IAA, with a further centre of divertsity in tropical Australian waters (Elfes et al., 2013). The IO (see Figure  1) is home to 21 species of viviparous sea snakes including seven endemic species (Elfes et al., 2013). Viviparous sea snakes give birth to live young, resulting in potentially low reproductive outputs and dispersal rates (Heatwole 1999) that may lead to rapid population subdivision (Lukoschek et al., 2007;Lukoschek et al., 2008). Thus, they are ideal candidates to examine historical biogeographic events and their role in generating biodiversity. Though recent studies have shed light on the origins, patterns of viviparous sea snake diversity in the IAA (Lukoschek et al., 2007;Lukoschek et al., 2008;Sanders et al., 2013b;Ukuwela et al., 2014;Ukuwela et al., 2016b), this work did not focus on the origins and patterns of viviparous sea snake species diversity in the IO.
In this paper we build on published work that focused on sea snake biogeography in the IAA, and summarise recent findings on the biogeographic origins of IO viviparous sea snakes and examine the effects of the past and present barriers to dispersal and gene flow between the IO and WP. Further, we review relevant literature on other marine groups that show similar patterns of colonisation and gene flow in the region.

MATERIALS AND METHODS
In this paper we describe the biogeographic origins of IO sea snakes and the effects of the past and present barriers to dispersal and gene flow between the IO and WP through reinterpretation of our data on the biogeographic origins of IAA sea snakes (Ukuwela et al., 2016) and the effects of the past and present barriers to dispersal and gene flow in the IAA (Ukuwela et al., 2014). For the inclusiveness, we have summarised the methods employed in Ukuwela et al. (2014 and2016) here and the readers are requested to refer the original papers if further details are necessary.

Phylogeny reconstruction and Divergence dating
To evaluate the origins and patterns of colonization of the IO sea snake assemblage, multi-locus time-calibrated phylogeny comprising 42 species (nearly 70%) of viviparous sea snake species were inferred using likelihood and Bayesian methods and their ancestral areas were reconstructed. The DNA sequence dataset consisted of 5792 base pairs, generated from three mitochondrial markers (cytochrome b, ND4 and adjacent tRNA region and 16S rRNA region) and five nuclear markers (two protein coding genes: cmos, Rag-1; three anonymous nuclear markers: G1888, G1894, G1914). Genetic sequence data was generated from whole genomic DNA extracted from liver/muscle tissue samples obtained from sea snake specimens collected from Australia to Iran (includes Sulawesi, Java, Borneo, Thailand, Myanmar, Sri Lanka, Bangladesh and India) ( Figure 1). The total dataset included 320 individuals of which 233 (36 species) were collected as fisheries by-catch during field surveys, 57 were from specimens accessioned in museums (22 species) and sequence data for 30 (16 species) were downloaded from the Genbank. Specimen collection localities, museum/ voucher numbers, details of DNA extraction, primers and PCR cycle settings are available in Ukuwela et al., (2016). Phylogenies were inferred using maximum likelihood and Bayesian analyses of the concatenated mitochondrial and nuclear alignment. Maximum Likelihood analyses (undated, no clock) were implemented in RAxML 7.2.8 (Stamatakis, 2006). Bayesian analyses with estimation of the divergence times were performed in MrBayes 3.2 (Ronquist et al., 2012). Details on outgroup selection, node calibration and selection of clocks are available in Ukuwela et al., (2016).

Ancestral Area Reconstruction
Ancestral areas were reconstructed using the timecalibrated phylogeny (MrBayes consensus tree) to examine the biogeographic history of IO viviparous sea snakes. Three oceanic regions were recognized based on other studies (VLIZ, 2009), which considered dispersal barriers (e.g. deep-sea trenches), patterns of endemism and species ranges across separate taxa. The three regions (Figure 2 inset map) are (1) IO, (2) SE Asia (SEA) and (3) Australasia (AUS). Ancestral area reconstructions (AARs) were performed using a Bayesian method implemented in BEAST 1.8 (Drummond and Rambaut, 2007), parsimony implemented in Mesquite 2.75 (Maddison and Maddison, 2009) and maximum likelihood with Dispersal-Extinction-Cladogenesis (DEC) implemented in LAGRANGE (Ree & Smith, 2008). For all analyses, each sample (tip or terminal taxon) was assigned to one of the three oceanic regions based on the collection locality. To test the importance of lineage-specific dispersal rates, we compared a model where different lineages (clades) were permitted different rates (using a 'random local clock') (Drummond & Suchard, 2010) to a simpler model which assumed a uniform dispersal rate across all lineages (a 'strict clock'). To test whether certain dispersal events were more likely, we tested four dispersal models of decreasing complexity: (1) a 'time-irreversible' model which assumed that all six  (Briggs, 1974) with the Indo-Australian Archipelago demarcated by the black dotted line. Grey areas denote the 120m isobath which indicates the extent of land that formed the Sunda shelf/Indo-Pacific barrier when sea levels were ~120 m below present levels during Pleistocene glacial maxima. Bathymetric data are from GEBCO (http://www.gebco.net/). dispersal events occurred at six different rates (AUS ® SEA; IO®SEA; AUS ® IO and the reverse), (2) a 'timereversible' model which assumed three such rates (AUS ↔ SEA; IO ↔ SEA; AUS ↔ IO) and (3) a single-rate 'unordered' model which assumed a single common rate for all six events. We further evaluated (4) a single-rate 'ordered' model, which permitted only dispersals between adjacent regions (AUS ↔ SEA; IO↔ SEA). The most appropriate of the above sets of models was selected using Bayes Factors. For further details of the analyses, see the materials and methods in Ukuwela et al. (2016).

Population genetic structure of Hydrophis curtus
To evaluate how past and present barriers to dispersal and gene flow have affected the species and genetic diversity of IO sea snakes, we examined the population genetic structure and recent levels of gene flow of the widespread sea snake, Hydrophis curtus (See Ukuwela et al., 2014 for further details). Population genetic structure was assessed by reconstructing mitochondrial DNA lineages and generating nuclear allele networks. To assess population subdivision and recent levels of gene flow between different populations of H. curtus, we used ten microsatellite markers. Whole genomic DNA was extracted from the liver/muscle tissue samples obtained from specimens collected from Australia to Sri Lanka through Sulawesi, Java, Thailand, Myanmar, and India. Three mitochondrial markers (cytochrome b, ND4 and adjacent tRNA region and 16S rRNA region) were used to reconstruct dated mitochondrial lineages and identify genetic structure. Two anonymous nuclear markers (G1888, G1894) were used to generate nuclear allele networks to identify patterns of genetic structure. The mitochondrial markers were concatenated and analyzed with time calibrated Bayesian inference in BEAST v.1.7.4 (Drummond and Rambaut, 2007) and maximum likelihood method in RAxML v. 7.2.6 (Stamatakis, 2006). Nuclear allele networks were constructed for the two nuclear markers independently using median-joining method (Bandelt et al., 1999) implemented in Network v.4.6 (fluxus-engineering.com). Ten microsatellite markers that were generated for H. elegans (see (Lukoschek and Avise, 2011)) and shown to cross-amplify successfully in H. curtus (Lukoschek and Avise, 2011) were amplified using multiplex ready technology. Each locus was amplified independently and after successful amplification, the PCR products for each individual were pooled and sent for fragment analysis. Allele sizes were scored using the software GeneMapper version 3.7 (Applied Biosystems). Population genetic structure was evaluated using the Bayesian population genetic assignment test in Structure v2.3.4 using all ten microsatellite loci (Pritchard et al., 2000). To estimate microsatellite differentiation between sampling locations and different populations, Fisher's exact probability tests and Fixation indices (F ST ) for all ten loci (combined) were estimated using Arlequin 3.5.1 (Excoffier and Lischer, 2010). See Ukuwela et al. (2014) for details of primers, PCR amplification settings, DNA sequencing, voucher numbers and Genbank accession numbers of the three mitochondrial markers and the two nuclear markers, and further details of analyses.

Origins, patterns of colonization and the assembly of the Indian Ocean sea snake fauna
All three AAR analyses indicate that the majority of IO sea snakes (

H. (Astrotia) stokesii and H. (Thalassophina) viperinus)
have a SE Asian origin and the regional sea snake fauna seems to be mainly derived from direct dispersal from SEA or either the WP. However, all three analyses suggest that the sister species H. stricticollis and H. obscurus represent an example of in situ speciation in the IO (Figure 2). Further, two species (H. fasciatus, H. spiralis) that co-occur in SE-Asian waters and the IO show possible recolonization of SE Asian waters from the IO (Figure 2). Divergence time estimates suggest that the majority of the sea snake species have colonized the IO in the last 3 to 0.5 my (Figure 2). Molecular phylogenetic analyses strongly indicated (Bayesian posterior probability (BPP) >0.9 and bootstrap support (BS) >70%) the presence of reciprocally monophyletic clades that correspond to IO versus SE Asian/WP populations within five of the co-distributed species (Microcephalophis (Hydrophis) (Figure 2). The corrected (HKY) pairwise genetic divergence in the cytochrome b gene (Cytb) between these sister lineages in the IO and SEA was 9.96-2.36%. These differences are broadly consistent with the divergence dating, which indicated that these allopatric intraspecific populations diverged between the last 2.6 to 0.5 Ma (Figure 2). The analysis also recovered distantly related cryptic lineages of H. cyanocinctus and H. ornatus with allopatric distributions in the IO or WP/SEA (Figure 2).

Population genetic structure and recent levels of gene flow of Hydrophis curtus
Bayesian and Likelihood analyses recovered two monophyletic clades corresponding to IO vs. mostly WP localities (BPP > 0.9, BS > 70) with deep basal divergence in H. curtus (Figure 3). The IO clade was separated from the WP clade by mean pairwise corrected cytochrome b genetic distances of 9.3-9.5%. The IO clade consisted of specimens from India, Myanmar and Sri Lanka. The WP clade consisted of three sub-clades with unresolved inter-relationships and relatively shallow mean genetic divergences (0.401-0.648%): (i) a Phuket-Thailand clade consisting of specimens from Phuket-Thailand; (ii) a SE Asian clade consisting of specimens from the south coast of West Java, East Java, south Sulawesi and Vietnam; and (iii) an Australian clade comprising specimens collected from two locations in northern and north-eastern Australia. Bayesian divergence time estimates suggest that the IO and WP clades diverged about 2.8 million years ago (Ma) (0.85-4.83 ma, 95% Highest Posterior Distribution (HPD) interval) (Figure 3). Nuclear allele networks demonstrated that Hydrophis curtus in the IO and WP did not share any alleles at either nuclear locus, indicating strong population subdivision between these two Oceanic regions.
Bayesian population genetic assignment test for combined microsatellite loci revealed the presence of four genetically distinct populations in the IO, Phuket, Thailand, SEA and AUS (Figure 4). Fisher's exact probability tests for combined microsatellite loci indicated highly significant (P < 0.05) population differentiation between all four populations. The overall F ST value for all ten loci between the IO and WP populations (combined Phuket-Thailand, SEA, Australia) was 0.174 and significant (P < 0.05).

DISCUSSION
Our findings indicate that dispersal and colonization from SEA are major sources of viviparous sea snake diversity in the IO. Analyses of the population genetic structure of Hydrophis curtus indicates deep basal genetic divergence in mitochondrial lineages between the allopatric populations in the IO and WP and the lack of contemporary gene flow among these populations. We discuss these findings with reference to the origins of the IO marine biodiversity and the geo-climatic history of the Indo-Pacific region (further details in Ukuwela et al., 2014;2016)

Biogeographic origins and assembly of the Indian Ocean sea snakes
Our Ancestral Area Reconstructions strongly indicate that the majority of the IO sea snakes have an immediate Southeast Asian/WP origin, i.e. all their nearest relatives are from SEA or the WP (Figure 2). However, the IO endemic clade of H. stricticollis and H. obscurus provides evidence for in situ speciation there, but this process has not made a significant contribution to the IO sea snake diversity (Figure 2). The origins of the other IO endemic species, H. biturbiculatus, H. hendersonii, H. mamillaris and H. nigrocinctus remains unknown as they were not included in the phylogenetic analysis due to the unavailability of fresh tissue samples.
There is some evidence of recolonisation of SEA/WP from the IO populations, e.g. in H. fasciatus. These events are rare and do not greatly increase the species diversity in SEA/IAA, and thus do not constitute strong support for the centre of accumulation model in IAA. Overall, the patterns suggest that most sea snake species in the IO colonised this ocean within the last 3 million years, and have allopatrically speciated from source lineages in the WP. Numerous studies on reef fish (Carpenter and Springer, 2005;Tornabene et al., 2015) and marine invertebrates (Veron, 1995;Barber and Bellwood, 2005) support the view that the peripheral areas such as the IO (compared to the IAA centre) most likely received species through dispersal and colonization from the IAA with relatively little in situ diversification. In contrast, there is evidence to support that the corals of the genus Acropora evolved in the IO or the Mediterranean region and have subsequently dispersed to the WP (Wallace and Muir, 2005).
Although it is difficult to provide a timing of colonization for the IO from the SEA/WP, the dates for population divergence for co-distributed species in the IO and SEA/ WP can be used as a proxy to determine broad timeframes for colonization of the IO. Corrected pairwise genetic (Cytb) distances of 9.96-2.36% between sister lineages in the IO and SEA/WP in certain species is consistent with late Pliocene-Pleistocene population divergences. This indicates that most species have invaded the IO within the last 3 million years and have allopatrically speciated in the intervening period, though many such lineages are   cryptic species (Carpenter et al., 2011;Hubert et al., 2012). The contribution of this process to the alpha diversity of viviparous sea snakes remains poorly quantified largely because these morphologically indistinct cryptic species/ lineages are not currently recognized as species.

Past and present barriers to dispersal and gene flow between Indian Ocean and West Pacific sea snakes
The deep genetic divergences between the IO and WP clades of Hydrophis curtus with mean pairwise corrected Cytb genetic distances of 9.3-9.5% and the lack of any shared nuclear alleles at the two nuclear loci examined indicates historical and long term isolation of IO and WP populations of H. curtus. Bayesian divergence estimates indicate that this population divergence took place approximately 2.8 Ma in the Plio-Pleistocene. This population divergence as in some of the other species of sea snakes mentioned above is most likely have been caused due to a vicariance effect and the divergence times also correspond to the onset of sea level fluctuations that caused low sea levels in the region ~2.6 Ma (Voris, 2000;Lambeck et al., 2002). Sea level drops of up to 120 m below present levels is known to have connected the Thai-Malay Peninsula and the islands of Java, Borneo and Sumatra forming the Sunda shelf/Indo-Pacifc biogeographic barrier (Voris, 2000;Lambeck et al., 2002;Woodruffe, 2003). Hence, it is possible that the Sunda shelf/Indo-Pacifc biogeographic barrier that formed during low sea level stands isolated populations of H. curtus and many other species of sea snakes spanning the IO and WP. Interestingly, this pattern has also been observed in recent phylogeographic studies of two other aquatic snakes: the salt water tolerant amphibious Cerberus rynchops (Alfaro et al., 2004;Murphy et al., 2012) and the viviparous sea snake Hydrophis cyanocinctus (Sanders et al., 2013b). Many marine fishes (Magsino and Juinio-Meñez, 2008;Gaither et al., 2011) and invertebrates (Lavery et al., 1996;Benzie, 1999;Crandall et al., 2008) also show this pattern of Indo-Pacific vicariance. These findings strongly suggest that the climatic and eustatic fluctuations during the Plio-Pleistocene have lead to population divergences that generated high levels of genetic diversity in many invertebrate and vertebrate taxa between the IO and the WP.
Bayesian population assignment and Fisher's exact probability tests for combined microsatellite loci indicated highly significant population subdivision with the presence of four genetically distinct populations in the IO, Phuket-Thailand, SEA and AUS. The highly significant F ST value for all ten microsatellite loci between the IO and WP populations (combined Phuket-Thailand, SEA, Australia), reciprocal monophyly in mitochondrial lineages, the clear separation of nuclear alleles and the population assignment findings strongly indicate the lack of contemporary gene flow between these populations despite the absence of current geographic/physical barriers for dispersal and gene flow. This lack of contemporary gene flow among IO and WP populations has also been demonstrated in many other reef-associated fauna (Yasuda et al., 2009;Leray et al., 2010). The habitat preference of shallow coastlines and the lack of larval stages that allow long distance dispersal may prevent sea snakes from crossing deeper seas and population connectivity reducing gene flow. Many reef fish and marine invertebrates have highly dispersal larval stages that facilitate long distance dispersal (Hoskin, 1997). However, some sea snakes (Hydrophis platurus, Hydrophis spiralis, Hydrophis fasciatus) (Ukuwela et al., 2016) and several marine species do not show any population structure across the Indo-Pacific (Horne et al., 2008;DiBattista et al., 2012). The lack of genetic structure in Hydrophis platurus is explicable as this is the only pelagic sea snake that tolerates deep seas and has demonstrably high dispersal rates (Ukuwela et al., 2016). The other two sea snake species that do not show genetic structure may indicate both recent colonization of the IO from the WP or recent recolonization of the WP from the IO. This scenario is also proposed for few of the reef fishes that do not show genetic divergence/structure across the Indo-Pacific (DiBattista et al., 2012). Further, the 'Coral triangle' in the IAA (Kochzius et al., 2009), Southern Thai-Malay peninsula, and the islands of Java and Sumatra have been shown to be regions of hybridization for many allopatric populations of fish in the IO an WP (Marie et al., 2007;Hobbs et al., 2009;Gaither et al., 2011). Further studies will be needed to identify areas of possible contact and hybridization among the allopatric populations of sea snakes. The concordant phylogeographic divisions in sea snakes and other taxa strongly support the view that populations in the IO and WP should be treated at least as different "evolutionary significant units" (Moritz, 1994) and hence distinct conservation management units.

CONCLUSION
Our findings indicate that dispersal and colonization from SEA is the major source of viviparous sea snake diversity in the IO. Divergence dates indicate that the majority of the viviparous sea snake species colonized the IO during the last 3 million years. The vicariance effect initiated by the Sunda shelf/Indo-Pacific barrier, which formed during the low sea level periods of the Pleistocene, seems to have played a significant role in creating cryptic species and lineage diversity in many sea snakes and other marine species. The low contemporary gene flow between certain IO and WP sea snake populations indicate the absence of current dispersal and population connectivity, despite the dissolution of this barrier. These genetically distinct lineages and populations in the IO and the WP qualify as Evolutionary Significant Units (ESUs) and thus should be considered as separate conservation and management units. Collectively, these findings indicate that the IO viviparous sea snakes have a long (in comparison to viviparous sea snakes that originated 10.6-6.5 Ma) and complex evolutionary history in the IO and hence have a unique conservation value. that supported him during his PhD study in Adelaide, Australia. The Indonesian Institute of Science (LIPI) and Department of Wildlife Conservation, Sri Lanka are thanked for the research permits. We thank the National Science Foundation and Man and Biosphere committee of Sri Lanka for organizing the National Symposium on Biogeography and Biodiversity Conservation in Sri Lanka. The editor Prof. Sumedha Madawala and two anonymous reviewers are thanked for their comments on the manuscript.