Genetic signature of Last Glacial Maximum regional refugia in a circum-Antarctic sea spider

The evolutionary history of Antarctic organisms is becoming increasingly important to understand and manage population trajectories under rapid environmental change. The Antarctic sea spider Nymphon australe, with an apparently large population size compared with other sea spider species, is an ideal target to look for molecular signatures of past climatic events. We analysed mitochondrial DNA of specimens collected from the Antarctic continent and two Antarctic islands (AI) to infer past population processes and understand current genetic structure. Demographic history analyses suggest populations survived in refugia during the Last Glacial Maximum. The high genetic diversity found in the Antarctic Peninsula and East Antarctic (EA) seems related to multiple demographic contraction–expansion events associated with deep-sea refugia, while the low genetic diversity in the Weddell Sea points to a more recent expansion from a shelf refugium. We suggest the genetic structure of N. australe from AI reflects recent colonization from the continent. At a local level, EA populations reveal generally low genetic differentiation, geographically and bathymetrically, suggesting limited restrictions to dispersal. Results highlight regional differences in demographic histories and how these relate to the variation in intensity of glaciation–deglaciation events around Antarctica, critical for the study of local evolutionary processes. These are valuable data for understanding the remarkable success of Antarctic pycnogonids, and how environmental changes have shaped the evolution and diversification of Southern Ocean benthic biodiversity.


Genetic diversity and population structure
Based on the sampling localities, COI sequences were grouped in one of four main SO geographical regions: Antarctic Peninsula (AP), East Antarctic (EA), Weddell Sea (WS) and Antarctic Islands (AI). Furthermore, the EA was divided into four subregions: Terre Adélie (TA), Ross Sea (RS), Davis Station (DS) and Bruce Rise (BR). Within the AI, there were two subregions: Bouvet Island (BI) and South Sandwich Islands (SSI) (electronic supplementary material, S1, S2; figure 1). Genetic diversity of N. australe was measured for each SO region and for the subregions within the EA using DNASP v. 5.00.7 [38]. Departures from values expected under panmixia (i.e. F ST = 0) among SO regions and within the EA, and the corrected p-values for population differentiation among pairs of populations were tested with ARLEQUIN v. 3.0.1 [39] with 10 000 permutations of the data. The partitioning of genetic variation among subregions within the EA was determined using analysis of molecular variance (AMOVA) based on 2000 permutations.
To infer the spatial genetic structure of N. australe, the number and the composition of panmictic groups as well as the spatial boundaries among them were estimated using a Bayesian model computed with GENELAND v. 2.0.0 [40] in the R environment (R, v. 3.2.3 [41]). The software implements a Markov chain Monte Carlo (MCMC) procedure to determine the best clustering of samples based on genetic and geographical information. Geographical information has been taken into account at the Bayesian prior level, so that clusters corresponding to spatially structured groups are considered to be more likely than clusters that are randomly distributed in space. Five million MCMC iterations sampled each 1000 steps with a 50 000 burn-in period, and a maximum number of clusters K = 10 were run to estimate the model parameters and posterior probabilities of group membership (P).

Demographic analysis
To detect demographic changes such as population expansion or contraction in N. australe at various spatial and temporal scales related to glacial refugia [42], the number of private haplotypes (PH, those endemic to subregions or regions) and the proportion of private haplotypes for each region (number of unique haplotypes/total number of haplotypes) was calculated. Tests for neutrality (Tajima's D and Fu's F S ) were run in DNASP to estimate deviations from the mutation-drift equilibrium and infer past population processes. The  resulting shape of the distribution was compared to a simulated dataset under a spatial expansion model [43,44] and a sudden population expansion. The overall validity of the estimated demographic model was evaluated by two different goodness-of-fit tests, tests of raggedness index (RAG) and the sum of squared differences (SSD) [45]. Significance of RAG and SSD were assessed by parametric bootstraps (10 000 replicates).
We used the coalescent-based method implemented in FLUCTUATE 1.4 [46] to estimate the exponential rate of population growth/decline relative to the neutral mutation rate (g) and the θ parameter (the effective population size scaled by the mutation rate, i.e. N e μ), the initial value of θ was estimated using the approach of Watterson [47]. All runs employed the following strategy: 1000 short chains with 200 generations and two long chains of 400 000 steps. Sampling increment was 20 for short and long chains. The transition/transversion rate was set to 3.99, determined by MEGA5 under the Tamura-Nei model. Runs were repeated five times to ensure consistency of estimates. Based on the significant genetic differentiation found among Antarctic regions, genetic panmixia cannot be assumed to analyse the demographic history of the whole Antarctic dataset; instead, each of the four regions is separately analysed. For the same reason, demographic analyses within EA were also run after excluding the significantly different DS samples (figure 1).
Molecular clock estimates were used to approximately date the timing of population expansions as in Rogers [48], calculating (T) for the analysed COI fragment based on T = τ /(2µk), where τ is the expansion time estimator, µ the mutation rate and k the sequence length [49]. The demographic expansion parameters of initial and final effective population sizes were estimated following θ 0 = 2N 0 u (before population growth/decline) and θ 1 = 2N 1 u (after the population growth/decline). Here, we apply a 10fold correction when estimating N. australe population expansion times, as it has been suggested that the short-term mutation rate is 10-fold higher than long-term substitution rates [50,51].

Genetic diversity and geographical structure
The final COI alignment of 364 sequences of N. australe had a length of 554 bp and represented 85 unique haplotypes (see the electronic supplementary material, S2 and S3) of which 10 were common haplotypes (represented by 10-87 individuals), while the majority were uncommon haplotypes with only one or two representatives (figure 2). We detected intermediate levels of genetic polymorphism, 59 nucleotides were variable and 35 were parsimony informative. The overall mean distance among sequences was 0.007 (range 0-0.02).
The genetic diversity varied across the SO regions. The number of haplotypes (k) was lowest in AI (k = 6) and highest in EA (k = 50), although this was probably related to the different sampling effort in each locality given that haplotype diversity was similar in these two locations (see the electronic supplementary material, S1 and S3). Genetic variation was low in the WS, which had around half the haplotypic diversity of the other localities, as well as the lowest average number of nucleotide differences (Π ) and mean nucleotide diversity (π ). The number and the proportion of private haplotypes were lowest in AI, although this was the region with the smallest sample size (table 1).
There was significant genetic differentiation among the four regions sampled around Antarctica (AMOVA, F ST = 0.424, p < 0.001). COI pairwise F ST relationships and significance comparisons supported the structuring of SO regions (table 2), with relatively high and significant F ST values among all localities. WS samples were the most distinct from all other regions (F ST between 0.51 and 0.68), while F ST values between EA and AI samples were the lowest, but still significantly greater than zero (F ST = 0.10).
A Bayesian phylogenetic analysis did not show segregated, well-supported clades related to region; however, the COI haplotype network coupled with the Bayesian consensus while supporting N. australe as a single species with a circum-polar distribution, also showed clear regional differences in the distribution of haplotypes (figure 2a). Most haplotypes were found only within a single region except for the AI where a large proportion of haplotypes were found elsewhere. Furthermore, the haplotype network had a predominant star-like pattern at the regional level, particularly for EA (figure 2b). Bayesian models computed with GENELAND detected three main clusters (K = 3, figure 3): cluster A included individuals from the WS and BI (figure 3a), cluster B grouped all EA individuals with the SSI (figure 3b), and cluster C included all AP individuals (figure 3c  Haplotypes are coloured by region and the size of their circle is proportional to its frequency in the whole sampling effort. White small circles represent haplotypes that have not been collected but should exist. AMOVA and F ST results within the AP and WS indicated low differentiation within regions, over 90% of the variation was owing to the differences within populations. Within the EA region, there was significant genetic differentiation among the five subregions sampled (AMOVA, F ST = 0.238, p < 0.001), although this was largely a result of the large and significant differences between DS samples and those from all other subregions (table 3). We also detected low but significant differences between the RS and other EA subregions (table 3), although the TA and BR populations were not differentiated (among these populations F ST < 0.025). Bayesian cluster analysis revealed two clusters (K = 2; figure 4), one included individuals collected in the vicinity of the Australian DS (cluster B, figure 4), and the other one grouped all other EA populations (BR, TA and RS) (cluster A, figure 4). Values of cluster membership were high for all localities (p > 0.90). Samples from DS were all collected in about 25 m depth, whereas most other EA samples were collected deeper than 200 m. A subset of the TA samples collected by the REVOLTA program from approximately 40 m depth showed no evidence of significant genetic differentiation from deeper TA samples (table 3), suggesting that segregation of the DS vicinity populations could be related more to geographical isolation than mere bathymetric differences.
Genetic differentiation within AP was significant owing to differences between northeastern and southernmost populations of the AP (F ST = 0.08, p < 0.05). The AI populations (BI and SSI) were also significantly different (F ST = 0.51, p < 0.05), while the WS did not show significant differences (F ST = 0.03, p > 0.05) between the easternmost and westernmost populations.

Demography
The significantly negative values of D and F S and unimodal mismatch distributions obtained for the whole dataset (Tajima's D = −1.842, p < 0.05; Fu's F S = −95.972, p < 0.001) seem to indicate that N. australe populations have gone through relatively recent expansions. We obtained a similar pattern in separate analyses of the WS, EA and AP data. Either an L-shaped or a unimodal distribution of pairwise differences was obtained for each region, supporting the existence of a demographic expansion in each region; this pattern was corroborated by the negative and significant Fu's F S test values and the nonsignificant goodness-of-fit tests. Tajima's D was not significant (p > 0.05) for AP and EA, suggesting there has not been historical reduction in effective population size in these regions (table 4; electronic supplementary material, S5). The N. australe samples from the islands (AI) showed non-significant neutrality tests, rejecting the population expansion model (table 4). Although the goodness-of-fit tests indicated that the model of rapid population expansion cannot be rejected, the calculated probability was very close to the significance level (p = 0.09 and 0.07 for SSD and RAG, respectively) (table 4). However, the unimodal mismatch distribution of pairwise differences suggested a recent demographic expansion (table 4; electronic supplementary material, S5). The estimates of change in population size relative to the mutation rate were consistent with previous analyses and ensured convergence on the correct parameter estimates for each region; replicate runs with alternate random seeds produced comparable results.

Discussion
The circum-Antarctic sea spider, N. australe is comprised of regionally distinct populations that appear to have undergone recent population expansion. Here, we show: (i) further evidence to confirm that N. australe, a benthic brooder, is indeed circum-Antarctic; (ii) that the AP and EA populations are genetically more diverse suggesting multiple demographic contraction-expansion events possibly associated with deep-sea refugia, while the low genetic diversity at the WS points to a more recent expansion possibly of shelf ice-free refugia; (iii) N. australe seems to have undergone demographic expansions during the mid-Pleistocene (15-21.2 kyr BP), suggesting multiple LGM refugia; (iv) support for population expansion in the AI associated with colonization from the continent, and (v) no genetic boundaries detected in East Antarctica, except for the clear segregation of samples from the vicinity of DS, an area particularly isolated from the more open habitats sampled off TA.

Nymphon australe, a circum-Antarctic species
Our study supports previous work proposing N. australe as a circum-Antarctic species [13]. Here, we found relatively high genetic homogeneity in COI among the N. australe individuals from across  Antarctica that together formed a single haplotype network precluding any suggestion of cryptic speciation (either geographical or bathymetric) within the species. Although the use of a single mitochondrial marker is often considered limiting in the context of detecting recent speciation and inferring spatial genetic structure [9,52], our efforts to incorporate additional genetic markers were futile as there was low or no intraspecific sequence divergence, even for nuclear markers that should be highly variable, e.g. microsatellites [32] and ITS which has been successful in segregating within-species clades in other pycnogonid species [16,31]. Thus, COI appears to be the marker of choice for interrogating genetic patterns in N. australe.
Pycnogonids in general are thought to be poor dispersers; for a benthic invertebrate with no pelagic stages, slow-moving, and with fertilized eggs and post-embryonic stages remaining attached to the father for some variable time, a circum-Antarctic distribution seems unusual. In other Antarctic pycnogonids studied, such wide distributions have been challenged, and restricted gene flow eventually leading to cryptic speciation has been proposed (e.g. Colossendeis megalonyx Hoek, 1881 [26,53]; Pallenopsis patagonica (Hoek, 1881) [15,24]); similar patterns are evident for many other Antarctic invertebrates too (e.g. Isopoda [54][55][56][57]; Amphipoda [58,59]; Ostracoda [60], Nudibranchia [61]; Crinoidea [52,62]). By contrast, data for Antarctic taxa with pelagic stages do not give evidence of geographical genetic structure and tend to reflect a circum-Antarctic distribution [9,63,64]. The absence of allopatric speciation in N. australe is thus surprising, and raises questions on how to reconcile a wide distribution, more characteristic of a species with pelagic dispersal, with the life-history traits of the species (benthic, crawler, late post-larval instars carried by the male [65]). The roles that environmental (e.g. ocean currents and ice movement) and ecological (e.g. hosts associations, drifting substrates, parasitism) factors may play in the maintenance of a circum-Antarctic species and Antarctic pycnogonid populations in general are yet to be studied. On the other hand, the clear genetic differentiation among Antarctic regions (F ST , table 2; haplotype network, figure 2) may well be seen as evidence of premature stages of speciation [13]. Unfortunately, no fossil records exist that help understanding the species history; assuming populations were established posterior to the LGM, it should be considered whether there have been sufficient generations for speciation to occur.

High genetic diversity and population structure of N. australe
The genetic diversity estimated for N. australe (H = 0.918; π = 0.00657) is higher than that for other invertebrates distributed throughout Antarctica including echinoderms [66] and arthropods [60] in which either distinct morphospecies are detected, or genetic homogeneity is associated with the occurrence of a pelagic stage [63]. The  the AP), higher haplotype and nucleotide diversities, high number of exclusive haplotypes and high proportions of their ancestral haplotypes might be related to region-specific conditions. On the other hand, the singular grouping of BI together with WS and SSI together with EA is rather unexpected. There is the assumption of lack of connectivity between the islands and the shelf [20] and Antarctic areas are seen as a separate zoogeographical region mainly influenced by the Antarctic Circumpolar Current (ACC). This pattern has been found in A. cornigera, in which distant populations on the Antarctic continental shelf (WS and TA) clustered together while the islands segregated as distinct clusters [25]. Dömel et al. argue that A. cornigera dispersal is limited by the ACC acting as major dispersal barrier, and by a relatively restricted depth range (max of 1180 m) limiting access to the deep sea. Here, the results could be seen as support for the hypothesis of connectivity between the SO islands and the continent [16,26] being enhanced by deep currents in eurybathic species such as N. australe (known max depth 4136 m); however, the large and significant F ST values found indicate regions are not connected. An alternate explanation for the grouping of samples from the islands with those from the continent is that these geographical populations are isolated and are in the process of diverging, but owing to post-LGM colonization, ancestral genetic signatures might remain.

Past population histories differ among regions
It is understood that the differences in mutation rates across lineages may impact estimates of molecular time-scales and demographic parameters from mitochondrial sequence data [67]. Based on published datasets, a 10-fold correction of the proposed mutation rate seems to be the most accepted in the literature (see [50,[68][69][70][71][72][73] among others). Applying the 10-fold correction infers a period of an order of magnitude more recent than applying no correction to the mutation rate. Here, we apply such correction to our estimates of expansion of N. australe populations, as corrected estimates (instead of 40-220 kyr BP with no correction applied) fit well with accepted hypotheses of benthic invertebrate refugia survival and population expansions associated with deglaciation events after the LGM about 20 kyr BP [9,25,61,71,74]. Although molecular clock estimates should be regarded with caution, and especially in Pycnogonida that still awaits dated phylogenetic hypotheses, the estimated timing of N. australe populations expansion (figure 5) matches that for other pycnogonid species [16,24] and other invertebrates (e.g. Nacella concinna [68]; Nematocarcinus lanceopes [63]; Pareledone turqueti [9]), as well as estimates of raising temperatures (end of the LGM) according to the Antarctic temperature anomaly (discussion in [75]).
There is strong evidence that, during the last glaciations, ice sheets extended to cover the continental shelves of Antarctica [76]; however, different rates of cover have been proposed for different Antarctic regions [11]. It has been suggested that the shelf benthic fauna was depleted during the LGM, but there are no sufficient data to fully understand the recolonization processes. Distinct scenarios have been proposed for explaining recolonization of the Antarctic shelf: (i) fauna found refugia ex situ, on the shelf of neighbouring continents or sub-AI [16,77]; (ii) fauna found refugia in situ on the continental slope and deeper waters of the SO [5]; and (iii) survival of fauna throughout the last glaciation in situ in shelf refugia [6,7,75,76,78]; a variety of taxa seems to agree with one of these scenarios ( [8,16,25] and others). Allcock and Strugnell [8] reviewed studies on genetic structure of SO organisms and suggested that much of this fauna would have survived the Quaternary glaciations in situ. Our N. australe findings agree with those views, but, different hypotheses apply (2 and 3 above) depending on the region.
The demographic history inferred for the EA and AP regions seems to follow a similar pattern of in situ survival. All demographic analyses (table 4 and figure 5) and the considerable degree of genetic diversity observed (table 1) suggest that these regions experienced a process of population expansion with no signs of historical bottleneck (Tajima's D test not significantly negative). These time estimates of expansion, 15.5 kyr BP at EA (13.6 excluding the differentiated DS samples), and 21.1 kyr BP at the AP, fit with the dating of post-LGM deglaciation events in these regions (approx. 14 kyr BP in EA and approx. 19.5-16 kyr BP in AP [75,79]). It has been hypothesized that well-grounded ice sheets across the continental shelf displaced organisms from the shelf, some finding refugia on the slope or deeper waters evading population bottlenecks [7,63].
It is plausible to assume that N. australe populations from EA and AP are genetically diverse possibly owing to repeated colonization events from shallow to deep and deep to shallow and changing landscape owing to ice expansion and retraction, supporting the hypothesis of multiple independent glacial refugia [25]. Also, eurybathic species like N. australe could be expected to retain high levels of genetic diversity [16,27,63], different to taxa bathymetrically restricted to shallow waters not migrating to deeper waters and that may have suffered severe population reductions during the glacial periods diminishing their genetic diversities [68]. Contrary to those from the AP and the EA, levels of genetic diversity in N. australe from the WS are rather low. The strongly reduced diversity, lower estimates of theta and higher growth rates (coalescent modelling) detected in WS populations compared to EA and AP (table 4 and figure 5) agree with the hypothesis of population expansion following events that might have reduced populations to one or very few in situ refugia, perhaps ice-free areas of the WS shelf [7,76,80]. The subsequent population expansion may have occurred rather recently (approx. 4.25 kyr BP) compared to other regions in which the expansion process occurred much earlier (e.g. EA and AP).

Colonization of the islands
The genetic signal from samples collected from BI and SSI (AI in figure 1) shows a different pattern, reflecting more recent colonization directly from the continent. None of the six haplotypes from the AI are shared among the two island populations. The haplotype network shows the three haplotypes from BI closely related to WS (haplotypes 30, 31, 36; figure 2), while the three haplotypes from SSI are closer to EA (haplotypes 1, 5, 37; figure 2). This result is also well supported by the Bayesian-based clustering analysis that segregated each island into two distinct groups (figure 4). The configuration of the AI haplotypes within the network suggests that BI haplotypes are more divergent from WS, than the SSI ones from EA. This could be explained by the remoteness of BI, located at approximately 1700 km from the continent or any other island [81]. It might be feasible that a small number of migrants from the WS (H31) colonized the isolated island aided by the Weddell Gyre [81] and later derived into H30 and H36 (haplotypes exclusive to BI). In the sea spider C. megalonyx a strong connection between populations from BI and WS is also proposed, but in contrast to N. australe, C. megalonyx from BI (n = 43) belong to a single haplotype that is also the most common in the SSI [16]. The N. australe results support the notion of a recent, 'founder effect' type of colonization of BI from the WS, rather than a potential insular refugium, a hypothesis already proposed for a crinoid species [52]. This scenario would agree with other genetic evidence, suggesting that ocean currents transport adults and larvae from Antarctic shelf areas to the outlying islands, rather than vice versa [7,16,27,77]. Larger sample sizes and better cover of neighbouring Antarctic continental shores are needed to better understand these colonization processes.
Haplotype sharing between the volcanic SSI and the EA region is somewhat unexpected for N. australe. Effect of homoplasy in an expanding population is a possibility, but also animals being transported from EA by oceanic hitchhiking or rafting, or by drifting on the ACC has been proposed [82]. SSI faunal composition has been shown as one of the most dissimilar when compared with the Antarctic continental shelf [7,20]; seemingly, the 'stop' at SSI might act as a filter retaining individuals, as suggested by the benthic insular refuge hypothesis [20,81]. In N. australe, this event seems relatively recent as haplotypes from SSI are mostly shared with EA haplotypes (figure 2). A similar pattern was found for the Antarctic octopus Paraledone turqueti [9].
The low genetic diversity together with indication of strong and significant population growth in N. australe from the AI may suggest expansion processes after the colonization from the continent, i.e. the expansion after the founder effect. Further sampling of N. australe from sub-AI and Antarctic archipelagos is needed for a better understanding of population patterns.

Genetic structure within the East Antarctic
The N. australe EA populations were represented by 50 distinct haplotypes leading to relative considerable levels of genetic diversity (table 1). The most probable ancestral haplotype (H5), and the star-like shape network (figure 2b) reveal a lack of geographical structure within the EA. Low F ST estimates among localities and Bayesian clustering corroborated a single group (i.e. BR, TA and RS localities), except for the DS samples. In contrast to the AP with its complex topography seemingly facilitating the heterogeneity of populations within the region [9,27], the EA populations appear as short branches of the network, suggesting either gene flow with small differentiation among localities or very recent patterns of divergence. Although the near shore EA area is complex with small rocky islands and fjords, depths increase rapidly to greater than 200 m where most of our samples were collected, reaching open basins of sedimentary substrate [83]. It is likely that such conditions are not necessarily barriers to gene flow and could lead to connectivity in organisms with a high bathymetric plasticity [82] as is the case of N. australe. On the contrary, species restricted to shallow waters might form distinct faunistic communities on rocks and algae [83] and may find more difficulty in circulating even short distances and thus presenting high levels of genetic substructure as observed in the amphipod Orchomenella [59]. The slight differences between the RS and the BR and TA populations might suggest either a glacial refugium in the RS coastal ice-free polynyas [78] or movement of TA individuals to the RS polynyas during periods of ice sheet coverage.
Pycnogonids from DS seem to represent the most genetically distinct population. DS is enclosed in a well-sheltered ice-free bay of 400 km 2 situated in the Vestfold Hills (see www.antarctica.gov. au). Segregation may be attributed either to a depth constriction (a maximum depth of 25 m) or to geographical isolation. This genetic differentiation of the DS samples from the remaining EA individuals is also found in other studies [52,58] and may be attributable to its local geographical and oceanographic conditions, given the extensive shallow areas of the Vestfold Hills coastal region trapping the shallow fauna forcing them to be inshore residents. Haplotype diversity in the DS population is rather low compared to the remaining EA areas, 10 out of the 13 individuals had the same haplotype H58 (figure 2b), probably because all samples were taken from a relatively small area by SCUBA.
Clearly, N. australe genetic structure within EA suggest that the eurybathic condition of the species has allowed gene flow over extensive areas of EA, as there is no evidence of bathymetric constraint for the migration of this species. Individuals collected by the REVOLTA program from TA at 40 m depth clustered with deeper waters samples from EA collected between 200 and 1230 m.
In conclusion, our results confirm that N. australe is one of the Antarctic brooding invertebrates with the widest distribution. The species shows different demographic histories depending on the Antarctic region, possibly shaped by the characteristics of the deglaciation events post-LGM. The study reinforces the notion of the strong effect of climatic events and environmental conditions on the patterns of diversity and structure in Antarctic benthic fauna. Nymphon australe is a key species for the understanding of microevolutionary forces that could reveal connectivity and dispersal mechanisms and demographic processes in the SO.
Data accessibility. The sequences generated in this study have been deposited in GenBank, accession numbers are available through the electronic supplementary material. The whole dataset (sequences file and individuals location data) created within this study has been uploaded to the Australian Antarctic Data Centre https://data.aad.gov.au/.