Closely related octopus species show different spatial genetic structures in response to the Antarctic seascape

Abstract Determining whether comparable processes drive genetic divergence among marine species is relevant to molecular ecologists and managers alike. Sympatric species with similar life histories might be expected to show comparable patterns of genetic differentiation and a consistent influence of environmental factors in shaping divergence. We used microsatellite loci to quantify genetic differentiation across the Scotia Arc in three species of closely related benthic octopods, Pareledone turqueti, P. charcoti, and Adelieledone polymorpha. The relative importance of environmental factors (latitude, longitude, depth, and temperature) in shaping genetic structure was investigated when significant spatial genetic structure was uncovered. Isolated populations of P. turqueti and A. polymorpha at these species’ range margins were genetically different to samples close to mainland Antarctica; however, these species showed different genetic structures at a regional scale. Samples of P. turqueti from the Antarctic Peninsula, Elephant Island, and Signy Island were genetically different, and this divergence was associated primarily with sample collection depth. By contrast, weak or nonsignificant spatial genetic structure was evident across the Antarctic Peninsula, Elephant Island, and Signy Island region for A. polymorpha, and slight associations between population divergence and temperature or depth (and/or longitude) were detected. Pareledone charcoti has a limited geographic range, but exhibited no genetic differentiation between samples from a small region of the Scotia Arc (Elephant Island and the Antarctic Peninsula). Thus, closely related species with similar life history strategies can display contrasting patterns of genetic differentiation depending on spatial scale; moreover, depth may drive genetic divergence in Southern Ocean benthos.

, while the distribution of suitable habitat or prey preferences may limit dispersal in some taxa (Cowan & Sponaugle, 2009;Rocha, Bass, Robertson, & Bowen, 2002). Nonetheless, barriers to dispersal in marine environments typically are less obvious than those in terrestrial and freshwater environments where the effect of many landscape features that potentially affect dispersal can be relatively well studied (Micheletti & Storfer, 2015;Van Strien, Keller, & Holderegger, 2012), often in multiple species (e.g., Hayes & Sewlal, 2004;Von Oheimb et al., 2013).
Using dispersal models is an apparently convenient method to identify potential barriers in the marine environment. Indeed, several studies on marine species have found reasonable congruence to the level of connectivity among populations inferred by oceanographic modeling and analysis of population genetic data (Galindo et al., 2006;Young et al., 2015). However, creating realistic oceanographic models to simulate dispersal for many marine species remains challenging for several reasons, particularly the comparatively limited biological and physical data that are available for most marine environments. Moreover, available data are likely confounded by seasonal and interannual oceanographic variation and inadequate knowledge about target species' life histories (e.g., data on timing and duration of spawning, fecundity, and mortality rates). In remote areas, such as the Southern Ocean, where obtaining relevant seascape and biological data is impractical, an alternative strategy to identify the presence of major dispersal barriers is to quantify spatial genetic structure in multiple species simultaneously.
Identifying the pathways and barriers to dispersal and then determining which features affect several species and which features are species-specific are central in understanding the processes that shape the rate of adaptation and speciation. For example, species-specific barriers would be suggestive of intrinsic responses to the landscape, while concordant breaks in spatial genetic structure across multiple taxa point toward a general lack of genetic exchange. Identifying these specific-versus general-dispersal boundaries can be used to inform spatial resource management (e.g., design of marine protected areas) (Cowan & Sponaugle, 2009).
High latitude marine areas, and notably the Southern Ocean, are particularly poorly studied in terms of the potential dispersal routes and barriers to benthic marine species. Most studies in the Southern Ocean have examined spatial genetic structure at a broad scale, with an emphasis on the degree to which the fast flowing Antarctic Circumpolar Current (ACC) limits dispersal between populations from South America and populations situated on the Antarctic continent/ sub-Antarctic (e.g., Hunter & Halanych, 2008;Thornhill, Mahon, Norenburg, & Halanych, 2008;Wilson, Schrödl, & Halanych, 2009).
A lack of understanding about spatial genetic structure at finer scales limits the extent to which we could, for example, designate effective protected areas or model a potential change in distributions in response to climate change.
The Scotia Arc region of the Southern Ocean represents an ideal area to investigate patterns of genetic structure across multiple species. Spatial genetic structure in the island arc region is expected to be complex: This region is subject to strong currents (due to its location within the course of the ACC), it is subject to a range of temperatures with different regions being influenced by the Pacific, Atlantic, and Southern Oceans, the ocean depths are variable (in excess of 4 km in places), and the Arc comprises several island groups of both continental and volcanic origin. In addition, the Scotia Sea has experienced rapid increases in ocean temperatures recently (Meredith & King, 2005;Whitehouse et al., 2008) making it an important location for investigating the impact of climate change. Oceanographic modeling indicates a predominantly unidirectional flow of water from the Antarctic Peninsula in a northeast direction toward Elephant Island, the South Orkney Islands and then toward the islands of South Georgia and Shag Rocks (Young et al., 2015); however, high resolution seascape data are not readily available for a formal landscape study of this area.
Studies on several marine invertebrates across the Scotia Arc region indicate that deepwater channels and the major currents between South Georgia and other sub-Antarctic islands and the Antarctic continent drive substantial population genetic differences between these locations (González-Wevar, Saucède, Morley, Chown, & Poulin, 2013;Hunter & Halanych, 2008;Krabbe, Leese, Mayer, Tollrian, & Held, 2010;Strugnell, Watts, Smith, & Allcock, 2012;Wilson, Hunter, Lockhart, & Halanych, 2007). Also, populations of direct developing (Margarella antarctica) gastropods around Signy Island are genetically different to those off the Antarctic Peninsula, as are populations of broadcast spawning (Nacella concinna) gastropods, with deepwater channels between these areas proposed to limit adult and larval movement, respectively (Hoffman, Clarke, Clarke, Fretwell, & Peck, 2011;. In addition, passive dispersal during the early life history stages can shape population genetic structure in Antarctic fishes. Model projections predicted that stronger genetic differentiation would occur in Champsocephalus gunnari (with a threemonth larval phase) than in Notothenia rossii (whose larval phase is estimated to be more than six months); this prediction was confirmed using microsatellite analyses (Young et al., 2015). Nonetheless, no study has yet quantified the congruence of fine scale barriers to dispersal across closely related marine species that inhabit the Scotia Arc region.
The Southern Ocean benthic octopods represent an ideal lineage to investigate whether congruent barriers to dispersal exist because they are benthic taxa that possess large eggs, which likely hatch as benthic young. It is well accepted that octopuses with eggs above a certain size possess crawl away young (see Boletzky, 1974). Given limited dispersal by adults and larvae, significant spatial genetic structure in response to the varied landscape of the Scotia Arc is expected.
Three closely related octopus species can be caught in large numbers across the Scotia Arc region. Pareledone turqueti (Joubin, 1905) is distributed across the Scotia Arc including the waters surrounding South Georgia and Shag Rocks. This species is circumpolar in its distribution and inhabits waters <1,116 m in depth. Adelieledone polymorpha (Robson, 1932) is distributed across the Scotia Arc but is absent from Shag Rocks. It occurs in waters <1,510 m in depth and also has a circumpolar distribution. Pareledone charcoti (Joubin, 1905) has a more restricted distribution than the other two species and is found off the northern Antarctic Peninsula and the South Shetland Islands (Allcock, 2005) in waters <286 m in depth; this species apparently has a restricted range as it is unknown from the waters to the east of Elephant Island including South Georgia, Shag Rocks, and the South Orkney Islands.
The aim of this study was to identify whether different octopus species, with putative similar life histories, share comparable levels of genetic diversity and spatial genetic structure across the Scotia Arc.  Table 1). The mean depth from which P. charcoti was collected at Elephant Island (111 m) was shallower than that for A. polymorpha (287 m) and P. turqueti (227 m) (Table 1). Similarly, P. charcoti was collected from shallower mean depths at the Peninsula (117 m) than A. polymorpha (294 m) while P. turqueti was collected in greater numbers in deeper waters (544 m) at this location. At South Georgia, A. polymorpha (192 m) and P. turqueti (194 m) were collected from comparable depths. Latitude, longitude, and depth were recorded for each sampling event. Temperature was recorded from most (78%) sampling locations via CTDs (Conductivity, Temperature, and Depth) that were deployed just prior to each trawling event. Tissue samples were stored in 70-95% ethanol at −20°C, prior to preserving whole animals in formalin for taxonomic identification (by ALA).
To assess possible source-sink population structure, we estimated the directional relative migration between sample locations (Sundqvist, Zackrisson, & Kleinhans, 2013). Briefly, the genetic composition of a potential pool of migrants was calculated as the geometric means of the allele frequencies of a pair of samples; next, the level of genetic differentiation (e.g., D, Jost, 2008)  Min., minimum depth; Max., maximum depth; Mean, mean depth; SD, standard deviation. All values are in metres. See Figure 1 for sample locations. estimates of genetic differentiation. This procedure is repeated for all pairs of samples, and the concomitant directional measures of genetic differentiation are used to calculate the directional relative migration among all sample pairs (Sundqvist et al., 2013). Directional relative migration was calculated using the function divMigrate within the R-package diveRsity (Keenan et al., 2013). Bootstrapping was used (1,000 iterations) to generate 95% confidence intervals and to determine whether migration is significantly higher in one direction than the other.
We used GESTE v.2.0 (Foll & Gaggiotti, 2006) to determine whether environment variation explains spatial genetic structure in P. turqueti (there was insufficient genetic structure for analysis of P. charcoti and A. polymorpha). GESTE uses a generalized linear model to relate environmental factors to a population-specific value of genetic differentiation (F ST ). Environmental factors included in the analysis were as follows: latitude, longitude, depth, and temperature (from CTD data). Ten pilot runs were used (length = 5,000) to determine acceptance rates for parameters of the Monte Carlo Markov Chain.
Subsequently, the analysis was run with a burn-in of 50,000, a sample size of 10,000 and a thinning interval of 20.  See Figure 1 for sample locations. Small samples from Signy Island for A. polymorpha (n = 2) and South Sandwich Islands for P. turqueti (n = 1) are not included.  (Table S13) (Table 3). Thus, the model clusters identified by STRUCTURE for P. turqueti showed some geographic structure at a regional scale, corresponding to (3) (6) Signy (Q 6 = 0.35) and to a lesser extent the Peninsula (Q 6 = 0.17) and Elephant Island (Q 6 = 0.12).

| Spatial structure
Regional genetic differences among P. turqueti samples were identified using BARRIER (Manni et al., 2004), with the second barrier occurring between Shag Rocks and South Georgia; subsequent barriers were identified between Signy and the Peninsula, and finally among the samples from the Peninsula itself (Supporting Information).

| Directional migration
There was an apparent directional bias to the relative migration in P. turqueti, with significantly higher migration rates from Signy Island into all other sample locations, particularly Elephant Island and South Georgia (relative directional migration >0.7), than in the opposite direction (Table 4); there was no evidence for significant differential migration between all remaining locations for P. turqueti. No significant asymmetry in relative migration rates was detected for samples of A. polymorpha or P. charcoti (Tables S17 and S18).

| Effect of landscape upon P. turqueti and A. polymorpha
GESTE analyses indicated a strong significant association between depth and genetic differentiation among samples of P. turqueti. When depth was included as a single factor in GESTE, a higher probability model was obtained than when depth was excluded (Table 5). In contrast, models containing latitude, longitude, or temperature as single factors had lower probabilities than the random effects model (termed "Constant" in Table 5). Inclusion of longitude, temperature, or depth as single factors in GESTE resulted in slightly higher probability models for A. polymorpha than when these factors were excluded (Table 5).
When two factors were included, the highest probability models for both P. turqueti and A. polymorpha included both factors and their interaction term. In each of these two factor comparisons, the individual factors contributed more than their interactions (Table 5). For P. turqueti depth contributed more than temperature, latitude, or longitude. By contrast, temperature and longitude contributed more to the model of genetic structure in A. polymorpha than did depth in each of these two factor comparisons (Table 5).
When latitude, longitude, and depth were included as factors in a three-factor model the highest probability model for P. turqueti included effects of longitude and depth, while for A. polymorpha the highest probability model included a contribution by longitude only.
Inclusion of all four factors in models of genetic structure did not produce a higher probability model than when all four factors were not included for P. turqueti, but for A. polymorpha the highest probability four-factor model included potential effects of longitude and depth with longitude contributing more than depth (Table 5).  (Tables S19 and S20).
T A B L E 3 Variation in genetic differentiation (F ST ) among sample locations in the Southern Ocean for three species of octopus, as measured using microsatellite loci

| DISCUSSION
The present study shows that closely related octopod species with similar life history characteristics display a generally similar pattern of genetic differentiation at large geographic scales, but contrasting patterns of population genetic structure at a regional scale across the Scotia Arc region. Local environment (associated with depth) appears to play an important role in driving spatial genetic structure in P. turqueti, and to a lesser extent in A. polymorpha. Moreover, no clear spatial genetic structure was evident in P. charcoti, a species relatively restricted in both geographic range and depth despite the samples of this species being separated by deep water.

| Variation in regional genetic structure
The reason(s) for the marked differences in regional spatial genetic structure between species is not clear. Although our target octopus species are closely related (Strugnell, Rogers, Prodöhl, Collins, & Allcock, 2008), we know relatively little of their basic biology.
Nonetheless, the large size (>10 mm diameter) of their mature eggs (Allcock, 2005;Allcock, Hochberg, Rodhouse, & Thorpe, 2003;Allcock & Piertney, 2002;Barratt, Johnson, Collins, & Allcock, 2008) suggests that they are all direct developers (Boletzky, 1974). Populations of sedentary or sessile marine invertebrates with nonpelagic larval dispersal typically exhibit more genetic differences than species whose larvae have a long duration in the plankton (reviewed in Selkoe & Toonen, 2011). Our study area contains numerous potential barriers to dispersal to benthic invertebrates. For example, the shallowest depths between Elephant Island and the South Shetland Islands are deeper than 500 m, with the shelf areas shallower than 286 m in depth (the deepest record for P. charcoti) separated by about 65 km.
The weak or lack of population structure between Elephant Island and the South Shetland Islands for A. polymorpha and P. charcoti is intriguing, particularly for the latter. Both P. turqueti and A. polymorpha have been caught from similarly deep habitats (1,116 m and 1,510 m, respectively) (Allcock, unpublished data in Strugnell et al., 2008) (also see Table 1), and thus, adult A. polymorpha might be capable of moving between these locations along the benthos. In contrast, P. charcoti inhabits relatively shallow depths (<286 m), with most individuals captured from ~110 m depth or less (Allcock, 2005) (Table 1) distribution why should P. turqueti have a more restricted dispersal than A. polymorpha? Preliminary studies suggest that these two genera may occupy distinct trophic niches (Daly & Rodhouse, 1994).
For example, the posterior salivary glands, which contain enzymes and venoms for subduing and digesting prey are much larger in A. polymorpha than P. turqueti; these species also differ in beak morphology, with A. polymorpha possessing a smaller, finer beak with a rostral tip that ends in a sharper point than P. turqueti Daly & Rodhouse, 1994). Daly and Rodhouse (1994) speculate that the less muscular body and atypical beak of A. polymorpha suggests specialization for exploiting an atypical resource, possibly in the water column, whilst the robust body and more typical octopod beak of P. turqueti may suggest hunting of wholly benthic prey. Therefore, niche divergence driven by diet may play a role in driving genetic structure in P. turqueti or the lack of it in A. polymorpha.
Although neither the adults nor juveniles of P. charcoti (or P. turqueti) are expected to move between Elephant Island and the South Shetland Islands, rafting on macroalgae (e.g., Leese, Agrawal, & Held, 2010) or sea ice (see Gutt, 2001 for a review) may facilitate gene flow in this species as it does in other benthic Southern Ocean species.
Anchor ice may also dislodge eggs attached to a benthic substrate that could then be moved by currents. Also, it is possible that populations of P. charcoti at Elephant Island and South Shetland Islands are isolated, with no contemporary dispersal, but that this species has not reached genetic equilibrium.
Low genetic diversity implies a low population size. Indeed, the genetically least diverse species P. charcoti has a limited distribution, occurring around the South Shetland Islands and off Graham Land (Allcock, 2005), compared with the most diverse species P. turqueti which has a circumpolar distribution (Strugnell et al., 2012) and also A. polymorpha (intermediate level of genetic diversity) which is distributed across the entire Scotia Arc (with the exception of Shag Rocks).
These interspecific differences in genetic diversity are reflected in the mitochondrial diversity, whereby just eight COI haplotypes (differing by at most by five substitutions) are known for P. charcoti, (Allcock et al., 2011), compared with some 35 COI haplotypes reported for P. turqueti (Strugnell et al., 2012). In addition, P. charcoti and another closely related species, P. aequipapillae Allcock, 2005 are estimated to have diverged relatively recently (~1 mya) (Strugnell et al., 2008), whereas P. turqueti is thought to be an older species, (the two main lineages within P. turqueti were estimated to have diverged ~4 mya [95% highest posterior density 1.9-7.1 mya]) (Strugnell et al., 2012).
One reason for a relatively narrow geographic and depth distribution and low levels of genetic diversity is that P. charcoti is a comparatively new species, supporting the idea that populations of this species are in genetic nonequilibrium (as discussed above). Alternatively, this species may have gone through a genetic bottleneck, potentially impacted by loss of habitat during the last glacial maximum (Thatje, Hillenbrand, & Larter, 2005).

| Structure between South Georgia and the South Shetland Islands
Octopus species show a more similar pattern of spatial genetic structure at a larger geographic scale, with genetic differentiation apparent between the islands of South Georgia and continental Antarctica for A. polymorpha and P. turqueti (P. charcoti's range does not extend to South Georgia). A barrier to gene flow and concomitant genetic differences between South Georgia and continental Antarctica are a consistent feature of the Southern Ocean seascape, having been reported in other benthic invertebrates, including those with a pelagic larval phase (e.g., González-Wevar et al., 2013;Hunter & Halanych, 2008;Krabbe et al., 2010;Wilson et al., 2007). While deep water and major currents are often proposed to act as dispersal barriers between continental Antarctica and South Georgia (González-Wevar et al., 2013;Strugnell et al., 2012), the former would be expected to be less of a barrier for species with pelagic larval phases than brooders such as octopods. The importance of life history and dominant ocean flows in shaping genetic structure across this region was recently demonstrated by Young et al. (2015) who found that differences in the length of the juvenile planktonic stage and the variability in ocean flows explained differences in genetic structuring between two species of teleost.
Moreover, relatively low expected heterozygosity for P. turqueti and A. polymorpha at South Georgia (and at Shag Rocks for P. turqueti) is consistent with these being small island habitats (Frankham, 1997) and at the range margin of both species (e.g., Arnaud-Haond et al., 2006;Lind, Evans, Taylor, & Jerry, 2007). Despite this low genetic diversity, the estimated prevailing direction of migration was from Signy outwards, to South Georgia and also from Signy to Elephant Island (directional relative migration >0.7), for P. turqueti. This asymmetric migration pattern from Signy to South Georgia is in accordance with the prevailing dispersal direction indicated from drifter buoys (Matschiner, Hanel, & Salzburger, 2009) and the pattern of dispersal detected for teleost fish (Young et al., 2015

| Population structuring by depth across Peninsula, Elephant Island, and Signy Island
An apparent effect of depth upon genetic structure of P. turqueti populations from the Scotia Arc is novel for Southern Ocean species.
Previous attempts to quantify the effect of bathymetry on genetic divergence in Southern Ocean taxa have been hampered by limited sample sizes and the common discovery of cryptic species (further impacting sample sizes) (Baird, Miller, & Stark, 2011). By contrast, genetic "isolation by depth" has been reported for deep-sea species, such as the giant amphipod, Eurythenes gryllus (France & Kocher, 1996), the bivalve, Deminucula atacellana (Zardus, Etter, Chase, Rex, & Boyle, 2006), the octocoral, Callogorgia delta (Quattrini, Baums, Shank, Morrison, & Cordes, 2015) as well as in teleosts inhabiting more shallow areas (Shum, Pampoulie, Sacchi, & Mariani, 2014). In Southern Ocean taxa, some evidence for genetic differentiation driven by depth is derived from divergent populations (from ~630 m and 123-540 m depth) of the giant Antarctic isopod, Glyptonotus antarcticus, but small (n = 2) sample sizes from the deep population prevent a robust analysis (Held & Wägele, 2005). Conversely, no evidence for genetic structuring by depth was found in the Southern Ocean nudibranch Doris kerguelenensis; however, a high incidence of cryptic species may have reduced the power to detect any intraspecific patterns of population structure (Wilson et al., 2009).
The effect of depth upon P. turqueti population differentiation is a likely consequence of one or more environmental or biological variables that correlate with depth, such as oxygen, type and availability of food, predation, disturbance, and/or temperature (see Gage & Tyler, 1991). Data for most of these variables are not available for the Antarctic region. However, inclusion of depth, temperature, and their interaction improved the fit of the model implying that temperature plays some role in shaping the genetic structuring evident within P. turqueti; indeed, depth and/or temperature may shape genetic differentiation in A. polymorpha, although the population sample number is limited for this species. Conversely, the relationship between population genetic structure and depth in P. turqueti could represent residual historic genetic structure from allopatric ice-free refugia (potentially present during the last glacial maximum) that occupied a range of depths. However, as genetic differentiation is related to bathymetry in some deep-sea species (France & Kocher, 1996;Schüller, 2011;Zardus et al., 2006), it seems likely that one or more (as yet unidentified) environmental factors (i.e. other than glaciation) can drive population structure in benthic marine invertebrates. An intriguing possibility is that this type of environmentally driven genetic structure in Southern Ocean species promotes allopatric speciation, whereby genetically different populations experience yet further divergence in refugia during glacial cycles.
We could not examine the effect of depth on P. charcoti as significant population differentiation was not detected. A lack of population divergence in P. charcoti, a species with a restricted geographic distribution and narrow bathymetric range, might reflect fewer opportunities for diversifying selection. Nonetheless, depth may have played an important role in the evolution of P. charcoti as a species, as each of the seven sympatric papillated Pareledone species from the Antarctic Peninsula region (all of which were previously ascribed to a single species) occupied in some cases distinct and relatively narrow depth distributions (Allcock, 2005). Therefore, it is possible that ecological speciation through niche divergence may have occurred within this clade of papillated Pareledone species, and this may have been at least partially driven by a similar process of genetic structuring across depth as is observed in populations of P. turqueti in the present study.
In conclusion, these data highlight marked differences in genetic structure among three closely related octopus species. The outcome of no detectable spatial structure in the species with the most restricted geographic and bathymetric range highlights a lack of knowledge about local scale population structuring in Antarctic taxa.
Moreover, these data highlight complexity of processes that impact population structure in a limited region in species that are apparently biologically similar. This highlights a potential difficulty in managing benthic Antarctic species from data on just a single "typical species." In addition, local associations between depth and genetic variation in P. turqueti may be suggestive of a previously unrecognized driver of ecological speciation in Southern Ocean benthos.

ACKNOWLEDGMENTS
We thank the Alfred Wegner Institute, British Antarctic Survey,