Geographic structure in the Southern Ocean circumpolar brittle star Ophionotus victoriae (Ophiuridae) revealed from mtDNA and single‐nucleotide polymorphism data

Abstract Marine systems have traditionally been thought of as “open” with few barriers to gene flow. In particular, many marine organisms in the Southern Ocean purportedly possess circumpolar distributions that have rarely been well verified. Here, we use the highly abundant and endemic Southern Ocean brittle star Ophionotus victoriae to examine genetic structure and determine whether barriers to gene flow have existed around the Antarctic continent. Ophionotus victoriae possesses feeding planktotrophic larvae with presumed high dispersal capability, but a previous study revealed genetic structure along the Antarctic Peninsula. To test the extent of genetic differentiation within O. victoriae, we sampled from the Ross Sea through the eastern Weddell Sea. Whereas two mitochondrial DNA markers (16S rDNA and COI) were employed to allow comparison to earlier work, a 2b‐RAD single‐nucleotide polymorphism (SNP) approach allowed sampling of loci across the genome. Mitochondrial data from 414 individuals suggested three major lineages, but 2b‐RAD data generated 1,999 biallelic loci that identified four geographically distinct groups from 89 samples. Given the greater resolution by SNP data, O. victoriae can be divided into geographically distinct populations likely representing multiple species. Specific historical scenarios that explain current population structure were examined with approximate Bayesian computation (ABC) analyses. Although the Bransfield Strait region shows high diversity possibly due to mixing, our results suggest that within the recent past, dispersal processes due to strong currents such as the Antarctic Circumpolar Current have not overcome genetic subdivision presumably due to historical isolation, questioning the idea of large open circumpolar populations in the Southern Ocean.


| INTRODUCTION
The Southern Ocean (SO) is characterized by rich biodiversity and largely endemic benthic fauna (Kaiser et al., 2013), resulting from an active geological history and organismal adaption to an extreme environment. While the Antarctic Polar Front (APF) serves to isolate the SO from warmer waters at lower latitudes, the Antarctic Circumpolar Current (ACC) is the world's strongest major current that has been presumed to aid dispersal of many marine species within the SO (Bathmann, Scharek, Klaas, Dubischarr, & Smetacek, 1997;Smetacek, De Baar, Bathmann, Lochte, & Rutgers Van Der Loeff, 1997;Thornhill, Mahon, Norenburg, & Halanych, 2008). The fact that the ACC can promote long-distance dispersal has helped reinforce the historically held assumption that many marine organisms of the SO likely have a circumpolar distribution around Antarctica (Dayton, Mordida, & Bacon, 1994). Antarctic currents closer to shore such as the Circumpolar Deep Water, Ross Gyre, and Weddell Gyre add complexity in predicting geographic dispersal capabilities of species (Tynan, 1998).
In addition to dispersal mediated by oceanic currents, glaciation cycles have also played a role in Antarctic biodiversity through controlling habitat availability (Thatje, Hillenbrand, & Larter, 2005). Glacial maximums during the Cenozoic likely forced species into the deep sea with pockets of refugia on the shelf allowing some species to recolonize and ultimately shape the SO's current community structure (Thatje et al., 2005). Polynyas, open regions of water surrounded by sea ice, in the SO may also serve as areas of refuge and often contribute higher levels of primary production (Massom & Stammerjohn, 2010). In expansion phases, grounded ice sheets can physically cover large geographic areas of the continental shelf, displacing inhabitants, and physically reshaping environments by removal and rearrangement of benthic habitat. During glacial contraction, new habitat becomes available allowing for population expansion. Thus, glacial cycles can drive population fragmentation and expansion opportunities (Thatje et al., 2005), ultimately serving as a biological "diversity pump" (Clarke & Crame, 1992).
Brittle stars are important members of SO biodiversity, comprising at least 219 nominal species and 126 that are endemic from the region (Martín-Ledo & López-González, 2014). Three of these species belong to Ophionotus (O. hexactis (Smith 1876), O. taylori McKnight, 1967;and O. victoriae Bell 1902); all of which also occur in the SO but are morphologically distinct from each other. Ophionotus victoriae Bell, 1902 is the most common and is a highly abundant (Figure 1), conspicuous ophiuroid endemic to the SO. This species has been reported to have a circumpolar distribution (Fell, 1961) and occupies many different benthic substrates within Antarctic waters (Fratt & Dearborn, 1984), with the South Sandwich Islands as its northern most limit (Sands et al., 2012).
Ophionotus victoriae has a long-lived, planktotrophic larvae, remaining in the water column for several months (Pearse, McClintock, & Bosch, 1991), thus allowing for the possibility of long-distance dispersal via the ACC. Previous phylogeographic work using the mitochondrial DNA (mtDNA) 16S ribosomal subunit (16S) and cytochrome c oxidase subunit I (COI) gene fragments reported unexpected levels of genetic diversity and divergence along the Antarctic Peninsula and oceanic islands (South Sandwich Islands and Bouvet Island), suggesting O. victoriae possesses higher-than-expected geographic structure and questions the possibility of cryptic species (Hunter & Halanych, 2010). Given this, a larger sampling effort around Antarctica would likely result in the uncovering of additional diversity and potential discovery of cryptic species.
To test for phylogeographic structure in this supposed circumpolar species, and to provide insight on processes of dispersal and historical isolation, molecular tools were used to examine O. victoriae over a >7,000 km range from the Western Ross Sea to the eastern Weddell.
This study, to the best of our knowledge, also includes the first sampling of benthic invertebrates from Wrights Bay, located between the Amundsen and Ross Seas. Herein, we utilized the mitochondrial 16S and COI genes to allow direct comparisons to results of Hunter and Halanych (2010) as well as a high-resolution whole-genome singlenucleotide polymorphism (SNP)-based approach, specifically 2b-RAD (Wang, Meyer, McKay, & Matz, 2012). This latter approach was chosen as restriction-associated DNA (RAD)-tags have been shown to identify fine-scale population structure in marine species beyond the resolution of mtDNA genes (Benestan et al., 2015;Reitzel, Herrera, Layden, Martindale, & Shank, 2013). Assessing population structure for organisms like O. victoriae of the SO is also important toward anticipating changes in the Antarctic benthic ecosystem as species ranges and structure will likely shift with future climate change (Aronson et al., 2007). Seas and oceanic islands, or a geographic distance of >7,000 km ( Figure 2 and Table S1). Samples available for 2b-RAD analyses included 96 specimens from 15 sampling localities ranging from the Ross Sea to the western portion of the Weddell Sea, a geographic distance >5,000 km.
Raw data were also processed and analyzed using the software Stacks To determine the potential number of populations (K), STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000) was utilized with the following parameters: (1) seven replicates at each potential K (1-15); (2) an admixture model with correlated allele frequencies; (3) a 50,000 repetition burn-in period; and (4) 100,000 additional Markov chain Monte Carlo (MCMC) repetitions. Because SNP datasets can vary in λ (parameter around the allele frequency prior) compared to mtDNA (Pritchard, Wen, & Falush, 2010), an initial run was used to infer λ to be 0.2447 prior to the full run. Resulting files were then processed with STRUCTURE HARVESTER (Earl & VonHoldt, 2012) to determine the most likely values of K from Delta K analyses as well as CLUMPP (Jakobsson & Rosenberg, 2007), and DISTRUCT (Rosenberg, 2004)  Furthermore, BayeScan (Foll & Gaggiotti, 2008) was utilized with four threads, 100 runs at a 100,000 burn-in length and 100,000 pilot length to identify any loci that might be under selection, and molecular diversity analyses were performed using GENEPOP (Rousset, 2008).  Sea, and oceanic islands) revealed 37.87% of the molecular variation as occurring between geographic regions (Table 2) closest to the origin (~1.8%-4.0% mean = 2.8%) represents comparisons between individuals restricted to subnetwork lineages I, II, and III illustrated in Figure 3a, whereas the other mode (~0.2%-1.6%, mean = 0.5%) are comparisons of individuals between subnetwork lineages. Given that these subnetworks largely correspond to geographic regions and given results of the 2b-RAD data (below), these latter two modes apparently correspond to intraspecific and interspecific variation, respectively.

| 2b-RAD analyses
Following quality filtering and SNP calling, 16,588 loci were recovered (Table 3). To further filter these SNPs, any loci not present in at least 80% of samples were excluded, resulting in 1,999 remaining SNP loci.
Next, any individuals with <80% of the total remaining SNP loci were excluded, resulting in removal of seven samples, thus leaving 89 individuals for analyses (Table 3). Under calculations of Delta K from STRUCTURE HARVESTER for this filtered and reduced dataset, K of 8 had the highest average support for Delta K after seven runs, although individual runs of K at 2 and 4 had the highest maximum-likelihood scores ( Figure S2). Thus, to assess which K was the most appropriate, STRUCTURE analyses were conducted with K set to 8, 4, and 2, then subjected to pairwise F ST tests. At K = 2 and 4, all populations were significantly different (Table 4)  depict K = 2 and K = 8, respectively).
Population structure was further investigated with SMARTPCA analyses that revealed geographic structuring in concordance with STRUCTURE. Figure 5 and F I G U R E 3 (a) Haplotype network of Ophionotus victoriae produced by PopART (Leigh & Bryant, 2015). The haplotype network is based off COI data from a TCS1.21 (Clement, Posada, & Crandall, 2000) analyses of 414 samples. Filled black dots represent missing haplotypes. In addition, maximum-likelihood analyses also revealed three clades. (b) Histogram of COI uncorrected pairwise distances (p)

| DISCUSSION
Both mtDNA and 2b-RAD data reveal considerable genetic structure across the Western Antarctic in the brittle star O. victoriae, questioning its current status as a single species. As mentioned, specimens were examined and morphological differences were not discernable with current taxonomy. However, both mtDNA and 2b-RAD data suggest distinct genetic lineages within what is currently recognized as O. victoriae. Although the Bransfield Strait appears to be more diverse then other populations indicating a possible mixing zone, the degree of genetic structuring appears ordered by major geographic regions.

| Phylogeographic patterns from mtDNA
Based on analyses of mtDNA COI, the western Weddell Sea has a recent shared history, or is currently connected with, the oceanic islands, and interestingly, the eastern Weddell Sea samples share a discrete haplotype subnetwork with the Ross Sea and Western Peninsula ( Figure 3a, lineage III, 2b-RAD data not available for eastern Weddell and oceanic islands). Lineage III, which includes the Ross Sea and eastern Weddell Sea samples, was also the most geographically widespread clade and yet the least variable, as no haplotypes are more than two steps from the most common haplotype. This particular T A B L E 3 Filtering steps for 2b-RAD SNP data  (Pritchard et al., 2000) and visualized in DISTRUCT (Rosenberg, 2004) testing for the true number of populations (K). K = 4 is presented in the graph above as our most likely accurate K Previous studies of other Antarctic benthic fauna have also revealed unexpected genetic structure in broadly distributed taxa. For example, the Antarctic crinoid Promachocrinus kerguelensis has a pelagic larval stage and was assumed to have a circumpolar distribution, but was ultimately found to be comprised of six different lineages and at least five different unrecognized species on the Antarctic Peninsula alone (Wilson, Hunter, Lockhart, & Halanych, 2007). Later analyses revealed all six lineages to be circumpolar, likely sympatric and eurybathic, with only two unrecognized species (Hemery et al., 2012).

| Phylogeographic patterns from 2b-RAD
Given that echinoderm larvae can remain in the water column for several months in the SO (Pearse et al., 1991) Water could have a significant impact on their distribution as well (Tynan, 1998).
The level of genetic differentiation recovered from 2b-RAD be- F I G U R E 6 Highest supported scenarios using Bayesian computation (ABC). In these scenarios, t# represents time in generations and is based off the four genetic populations identified by STRUCTURE. All three geographic regions split at approximately the same time with a more recent diversification in the Weddell Sea Bellingshausen Seas likely resulted from shifts in the position of the ACC provide a plausible explanation for the genetic differentiation recovered. Antarctic coastal currents likely are a factor in this isolation as they move the opposite direction of the ACC and have been used to explain genetic structure in benthic invertebrates thought to be circumpolar (Riesgo, Taboada, & Avila, 2015).

| Overall structure
The high-resolution 2b-RAD approach was consistent with findings from COI while providing greater genetic resolution. Although 2b-RAD data were more geographically and numerically limited relative to COI, both recovered strong connections between the Ross Sea and Western Peninsula, a distance of over 5,000 km, bypassing the Bellingshausen and Amundsen Seas. This connection is likely the result of transport from the Ross Gyre into the ACC, which does not contact the Antarctic shelf again until the western portion of the Antarctic Peninsula (Tynan, 1998 Probable causes for this diversity may be mixing of distinct water masses, including water from the ACC, and thus populations in the region (Gill, 1973;Smith, Hofmann, Klinck, & Lascara, 1999), repeated formation, and disintegration of refugia during glaciation events (Clarke & Crame, 1992) or other processes that promote mixing of populations. Although the ACC serves as a vector for eastward distribution, westward counter currents closer to the shelf and several large Antarctic gyres (including in the Weddell and Ross Seas) may further distribute, or isolate, populations (Thatje, 2012). For example, the Weddell Gyre moves clockwise and spills into the Bransfield Straits mixing with warmer waters (García et al., 2002;Kaiser et al., 2011).
Such a situation supports our findings for isolation of the WA/BS population. Bransfield Strait consists of several water masses differing in oxygen and salinity (Gordon, Visbeck, & Huber, 2001) compared to those on the western Antarctic Peninsula, along with being hypothesized as an area of refugium during glaciation events (Jażdżewska, 2011

| Comparing mtDNA to 2b-RAD
Our study afforded the opportunity to compare traditional mtDNA markers to a whole-genome SNP-based approach such as 2b-RAD.
Other studies have recognized the ability of RAD data to recover structure that traditional markers have overlooked (Reitzel et al., 2013;Wagner et al., 2013). Although general phylogeographic structure of large-scale SO regions was able to be ascertained through mtDNA, identification of four distinct populations and existence of population WB would have gone unrecognized if higher resolution 2b-RAD analyses had not been employed. Population WB comprised 16 specimens whose haplotypes were within COI lineage I, with an additional four individuals from lineage II, which could have resulted from incomplete lineage sorting as mtDNA is uniparentally inherited.
The 2b-RAD data were able to reveal this higher resolution structure through fewer samples from a smaller geographic range. Further 2b-RAD data in the sub-Antarctic islands and the eastern Weddell Sea would provide greater insight into connectivity of organisms in the SO ecosystem with the Ross and Weddell Seas being of particular interest. As marine systems are often considered to have few barriers, these high-resolution approaches provide us with better tools to answer ecological questions. With climate change prone to reshape current community structure in the SO ecosystem, large high-resolution phylogeographic studies can help to serve as a benchmark or snapshot prior to any further restructuring.