Antarctic Toothfish Dissostichus mawsoni in the South Orkney Islands: Using Stable Isotopes and Otoliths to Test Current Hypotheses About Nursery Areas and Demographic Units

1. Introduction

The Antarctic toothfish Dissostichus mawsoni is a benthopelagic member of the family Nototheniidae, which presents a circumpolar distribution, mainly south of 57°S and closely related to the Antarctic Convergence [1,2]. With a median age and size at first maturity of 13-16 years-old and ~100 cm total length (TL), respectively, its lifespan reaches 50 years, exceeding 200 cm TL and 100 kg of individual mass [3,4]. Several ontogenetic habitat shifts occur during this long lifespan [5,6]. While young-of-the-year (YOY) and yearlings (4-12 cm TL) use shallow pelagic waters close to the coast, older juveniles (19–54 cm TL) become benthic and spend a number of years on the continental shelf at depths between 50 and 300 m [7,8,9]. Subadults (50-95 cm TL) move then to greater depths (500-850 m) on the continental slope [10]. After reaching sexual maturity (~100 cm TL), adults become neutrally buoyant [5,11,12], and expand their bathymetric range greatly, feeding between 10 and 2200 m, although showing some preference for depths between 300 and 500 m [2,13,14].

Dissostichus mawsoni supports one of the two most valuable commercial fisheries in the Southern Ocean, largely in international waters managed by the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR). CCAMLR jurisdiction is divided into Areas, Sub-areas and Small-Scale Research Units (SSRUs; Figure 1), which are expected to share relatively similar oceanographic and biological characteristics, and to contain relatively discrete populations of some target species [15]. For some others, these CCAMLR divisions may not represent biologically relevant (evolutionary or demographic) units, leading to potentially important mismatches between biology and management actions [16,17]. In the particular case of D. mawsoni, both the number and distribution of evolutionary and demographic units across the Southern Ocean and the degree of matching between these units and CCAMLR areas or sub-areas are poorly known. This uncertainty and incomplete knowledge about the structure, connectivity and dynamics of D. mawsoni demographic units have become core issues in recent critics to CCAMLR regulations, considered too risky given several key parameters of the stock assessment model [18], estimated considering, almost exclusively, biological and fishery data from sub-area 88.1, remain unknown and/or very uncertain [19].

The limited knowledge of D. mawsoni evolutionary units results not only from the small number of studies conducted to date but also from variability in sampling designs, study areas, and molecular techniques [20,21,22]. Tagging studies have shown that while most marked adults remained for years within 50 km of their tagging sites, a few others moved up to 2300 km away [7,23]. Despite its low frequency, these movements could be enough to obscure or cancel out ongoing genetic differentiation processes across very extensive areas [24], biasing the estimation and/or leading to the actual existence of a rather small number of evolutionary units. Nonetheless, the dominance of a resident behavior suggests the existence of separate D. mawsoni demographic units, isolated enough to exhibit asynchronous population dynamics and different life-history parameters [16].

Despite its practical relevance for management purposes, unveiling the number, distribution, and connectivity between D. mawsoni demographic units (i.e. stocks) has received little attention so far [23,25,26], focused on sub-areas 88.1 and 88.2, which concentrates most of the current fishing activity. Here, evidence from otolith chemistry, tagging, biological, fishing, and particle dispersal studies have shaped the current working hypothesis that a single demographic unit exists in these two sub-areas [7,27]. A second demographic unit would exist, nonetheless, west of the Antarctic Peninsula [25,26]. No available studies have been devoted explicitly to scrutinize the structure of demographic units within other areas or sub-areas. This is the case of area 48, which goes from the Antarctic Peninsula to Prince Edward and Marion Islands (Figure 1), and represented 14% of the fishery catch in 2018. Given the overexploitation of several stocks, most of this area has remained closed to all finfish fisheries since 1990 [28].

While fishery-dependent data is scarce after years of fishing closure, fishery-independent information is also very limited for D. mawsoni in Area 48 [29]. Nonetheless, After significant efforts, these authors have proposed three main hypotheses about the structure and connectivity of D. mawsoni demographic units in this area: i) a single demographic unit exhibiting some limited connectivity with area 58 (Indian Ocean sector), ii) two separate units, both placed within the Wedell Sea, with origins in sub-areas 48.2 y 48.6, and iii) at least three interconnected demographic units, with origins in sub-areas 48.1, 48.2, and 48.6 and some limited connectivity with area 58. In the present work, we confront these three hypotheses against otolith chemistry data collected around the South Orkney Island (sub-area 48.2) during exceptional exploratory activities authorized in 2016 and 2018.

Otolith chemistry has been proved useful to identify and determine demographic units, migration, and habitat use patterns in multiple species, at different temporal scales, given enough variability exists between the habitats occupied by the different groups [30,31,32]. In this particular application, we use observed variability in elemental and isotopic compositions of otolith nuclear regions to estimate the most informative number of nursery origins and apply isotopic-based geo-location techniques [33,34,35] to identify more likely origins for the sampled fish. Moreover, by comparing chemical signatures of nuclear and edge regions with expectations derived from current knowledge about ontogenetic migrations of D. mawsoni [2,7] and known effects of salinity, temperature and other physiological variables about elemental and isotopic signatures in fish otoliths [36,37].

2. Materials and Methods

2.1. Sampling

Sagittal otoliths of 45 adult D. mawsoni were collected around the South Orkney Islands (Figure 1) during two exploratory fishing trips, authorized exceptionally to be conducted by the Chilean FV Cabo de Hornos in the closed sub-area 48.2, during 2016 and 2018. A first set of 21 samples across 8 different tows was collected between February 25 and March 3, 2016. A second set of 24 samples from another 8 tows was collected between February 4 and February 12, 2018.

2.2. Chemical Composition of Otoliths

Otoliths were extracted, rinsed with seawater, dried with a paper towel, and stored in paper bags at sea. Once at the laboratory, they were rinsed again using ultrapure water (18.2 MΩ), photographed, sonicated in ultrapure water for 5 min and dried under a laminar flow hood, following clean protocols, using clean plastic forceps and disposable gloves. Otoliths were then embedded in high purity epoxy resin (EpoxiCure TM2, Buehler) and cut (low-speed ISOMET Buehler saw) at both sides of the primordium to obtain transverse 450-900 μm sections. Sections were first polished with progressively finer silicon carbide cloths (1500, 2000 and 2500 grit) down to a thickness of 200-400 μm (Figure 2), and then with 60 μm and 3 μm diamond lapping film (Buehler) to expose primordia and remove contamination from silicon carbide cloths. Once polished, sections were sonicated for 5 min in ultrapure water, dried under a laminar flow hood, mounted on slides, using thermoplastic glue (Crystalbond™), and stored for further analyses.

2.3. Elemental Composition Sampling and Analysis

Elemental composition sampling and analysis were carried out by a 193 nm ArF laser ablation (LA) system (Photon Machines Analyte G2) coupled to an iCapQ ThermoScientific ICP-MS, following a radial transect of spots, from the primordium to the ventral or, if not possible, dorsal edge of the section (Figure 2A). Distance between ablations was accommodated slightly (93-97 μm) to cover the whole transect, yielding 14-24 ablations per otolith, depending on its size. Given a fluence of 5 J·cm-2, a frequency of 10 Hz, and a period of 30 s, each ablation reached a depth of 20-30 μm, with an average diameter of 65 μm. Silicate glass reference material NIST SRM 612 was used as the primary external standard [38], while the USGS synthetic calcium carbonate MACS-3 [39] was employed for quality control. Ca was used as the internal standard, assuming a stoichiometric value of 38.8% for CaO in the aragonite matrix [40]. Ablation points within this region with elevated Al concentrations were excluded from the analysis to avoid potential contamination

External standard measurements were performed before and after ablating each otolith. Data generated by the mass spectrometer were reduced by using the Iolite V2.5 kit [41], following its preset internal and external standard method (“Trace Elements IS” reduction scheme). In total, concentrations of 38 elements (40 mass numbers) were quantified: 7Li, 11B, 23Na, 24Mg, 25Mg, 27Al, 31P, 52Cr, 55Mn, 57Fe, 59Co, 60Ni, 63Cu, 66Zn, 75As, 83Kr, 85Rb, 86Sr, 88Sr, 89Y, 111Cd, 120Sn, 121Sb, 125Te, 133Cs, 138Ba, 139La, 140Ce, 141Pr, 146Nd, 147Sm, 169Tm, 174Yb, 175Lu, 178Hf, 206Pb, 207Pb, 208Pb, 232Th and 238U. However, only eight elements exhibited concentrations above detection limits in ≥95% of the samples and were selected for further analysis: Li, Na, Mg, Cr, Mn, Sr, Sn, and Ba.

Although LA transects extended from the core to the edge of each otolith section, we focused, for most analysis, on nuclear and marginal regions (Figure 2A). The nuclear region included from the 2nd to the 6th ablations (93-540 μm from the core) and was set to represent most of the first year of life given reported measures of 623-795 μm for the first annuli radius [8]. The first ablation here was purposely excluded to eliminate possible maternal effects. The edge region included the last three ablation points (~255 μm) from each LA transect, representing the last 5-7 years of life. No finer resolution was possible here given very slow otolith accretion rates in adult D. mawsoni and minimum material requirements for the isotopic analysis (~30 µg). Ablation points within this region with elevated Al concentrations were excluded from the analysis to avoid potential contamination. Fish age and cohort (birth year) were estimated through the inverse function of the von Bertalanffy growth model reported by Horn [4].

2.4. Stable Isotope Analysis

Powder samples (≥30 μg) from nuclear and marginal regions, again, designed to represent juvenile (0-540 μm from the core) and adult life-stages (255 μm from the edge) were obtained using a New Wave Research microdrill equipped with a 500-μm diameter bit (Figure 2B). Extracted material was analysed at the University of Arizona Environmental Isotope Geochemistry Laboratory using an automated KIEL-III carbonate preparation device coupled to a Finnigan MAT 252 ratio-gas spectrometer. Digestion was undertaken in vacuum dehydrated phosphoric acid, at 70^oC, and the CO² generated by this reaction used to determine δ¹⁸O y δ¹³C values (in parts per thousand) relative to Vienna Pee Dee Belemnite (VPDB), using international standards NBS-19 y NBS-18. Whenever possible, the same otoliths and sections were used for both elemental composition and stable isotope analyses to ensure consistency.

2.5. Statistical Analysis

As a general approach, we purposely avoided any dichotomic use of p-values, following recommendations by the American Statistical Association [42], as well as by a growing number of independent scientists [43]. We favoured, instead, the confrontation of alternative models to the available data [44,45], as well as the presentation of effect sizes, standard errors and related probabilities, whenever possible. Further details provided for each goal as follows.

2.5.1. Confronting Ontogenetic Habitat Shift Hypotheses Against Otolith Chemistry Data

Elemental and isotopic signatures of nuclear and marginal regions were compared using univariate and multivariate techniques. Univariate analyses followed a mixed linear models approach [46], as implemented in the R package ‘nlme’ [47]. Multivariate analyses were based on a permutational multivariate analysis of variance (PERMANOVA) as implemented in the R package ‘vegan’ [48], using Mahalanobis distances. Residual quantile-quantile (QQ) plots showed no severe departures from univariate and multivariate normal distributions. Consistency with univariate and multivariate homoscedasticity assumptions was verified using median-based Levene’s [49] and permutational Mahalanobis distance-based [50]. The lack of independence expected between observations from the same fishing event was accounted for by incorporating this variable as a random factor in univariate models and as a fixed one in multivariate ones.

Sr:Ca y Ba:Ca ratios showed the greatest variabilities between and within our fish samples. Since these two ratios have been found to reflect environmental variability in temperature, salinity and/or freshwater inflows [36,51,52], we used them to identify and date discrete habitat shifts experienced by each sampled fish along its ontogeny. Habitat shift probability was estimated for each ablation along the variability within individuals, along the LA-ICPMS profile, using the Bayesian point change method of Barry & Hartigan [53], as implemented in the R package ‘bcp’ [54], running 10 000 simulations, burning the first 500, and using a prior value of 0.2 for both the signal/noise ratio and the probability of shift detection.

2.5.2. Estimating the Number and Mixing Between Nursery Areas

We fitted a series of finite mixture distribution models [FMD], [55], defined by a variable number of 1-5 sources, to either elemental or isotopic signatures observed in otolith nuclear regions. Source contributions, means, and covariance matrix parameters for each source were estimated using the expectation-maximization (EM) algorithm [56], as implemented in the R package ‘mclust’, restricted by us to consider only uneven-volume covariance models [57]. Model selection was performed using Schwarz [58],’s Bayesian Information Criteria, following the general approach of Niklitschek & Darnaude [59]. A few fish from putative cohorts born between 1940 and 1988 and later than 2002 were excluded from this analysis to reduce potential effects from inter-annual variability. The resulting subset of 37 and 39 samples for elemental and isotopic signatures, respectively, were standardized to avoid scaling effects while fitting the models, Nonetheless estimated means and standard errors were re-scaled afterword back to their original scales.

2.5.3. Geo-Locating Nursery Areas to Confront Söffker et al. (2018)’s Hypotheses

To hind-cast spatially-explicit probabilities of previous occurrence for each individual across our study area [60], we compared individual δ¹⁸O values observed in the nuclear region of each otolith (δ¹⁸O_OTO1) to the expected distribution of such values across the study area (Figure 3). To produce these expected values we assumed otolith aragonite precipitated in equilibrium with seawater (δ¹⁸O_WATER), following the empirical relationship between biogenic carbonates and δ¹⁸O_WATER derived by Martino et al. [61]. This equation outperformed commonly used fractionation relationships [62,63,64,65] as it yielded a lower mean bias (-0.60‰) than other tested equations, whose mean biases ranged between -0.77‰ [62] and -1.53‰[64]. To perform this bias testing, we compared the δ¹⁸O values observed in otolith edges (δ¹⁸O_OTO2) to the theoretical δ¹⁸O_OTO2 values predicted from the δ¹⁸O_WATER isoscapes built using each alternative equation. To do so, δ¹⁸O_WATER isoscape values (see below) were averaged over a radius of 100 km and a depth layer of 100 m around each capture location.

δ¹⁸O_WATER isoscapes were built by interpolating (Bayesian kriging) all δ¹⁸O_WATER data available for the Southern Ocean at both the Global Seawater Oxygen-18 Database [66] and at the World Ocean Database [67]. As in many other domains, these δ¹⁸O_WATER data were limited both in time and in space [66]. Several studies have shown, however, that δ¹⁸O_WATER can be predicted from salinity, temperature and/or other physical variables, within oceanographic domains [68,69]. This was also the case here, where salinity, temperature, and depth explained 61% of δ¹⁸O_WATER deviance. Considering both the nonlinear effects of these covariates and the spatial correlation between observed δ¹⁸O_WATER values we fit a Bayesian geo-statistical model [70], and used it to build isoscapes of predicted (interpolated) δ¹⁸O_WATER values, across the whole study area, considering six depth strata: 0-50, 50-100, 100-200, 200-300, 300-400, and 400-500 m, defined to exceed the expected bathymetric range (90-400 m) reported for D. mawsoni juveniles [5,7]. World Ocean Database values of temperature and salinity were previously interpolated using DIVAnd [71]. The temporal domain used for inference about juvenile life (nuclear region of the otolith) corresponded to the yearly interval 1989-2002, set to match the putative birth years of the sampled fish. The temporal domain used for testing bias related to fractionation equations was 2011-2018, set to match the last six years of life of sampled fish (time period represented in our edge samples).

The isoscape δ¹⁸O_W values (mean and variance estimates) were then used to produce probability distributions of expected δ¹⁸O_OTO1 values at each location and depth strata (Figure 3D). By confronting the observed δ¹⁸O_OTO1 values to these expected δ¹⁸O_OTO1 distributions we, finally, estimated the individual probability of occurrence of each fish at each location and depth strata during its YOY stage. The individual probabilities that resulted from this exercise were then averaged to produce synoptic maps within and across depth layers (Figure 3F), and to compare such probabilities between the four subpopulation areas delimited by Söffker et al.’s Hypothesis 3 [29].

3. Results

3.1. Ontogenetic Habitat Shift

The highest relative concentrations found in D. mawsoni otoliths corresponded to Na and Sr, followed by Ba, Cr, Li, Mn, and Sn (Table 1). Mean relative concentrations of Na, Sr, and Ba were 1.2, 3.5, 3,0 and times higher in marginal regions representing recent adult life, than in nuclear regions representing juvenile (YOY) life (Table 1, Figure 4). The mean relative concentration of Mn, on the other hand, was 38% lower in the marginal region, while none to small differences (<10%) were found in Li, Cr, Mg, and Sn values. Multivariate analysis showed these overall differences in elemental composition between nuclear and marginal otolith regions were large and consistent (n=41, Pseudo-F=10.604, p>F <0.001), revealing clear segregation between juvenile and adult signatures (Figure 5A). Large isotopic differences were also found between otolith regions, with mean δ18O and δ13C ratios being, respectively, 20 and 58% higher in marginal regions representing adult life stages (Table 1, Figure 4). These differences between otolith regions were also evident under a bivariate analysis (n=45, Pseudo-F=45.245, p>F <0.001, Figure 5B).

Most individual Sr:Ca profiles (~95 %) showed low probability (<0.2) of habitat shift before the tenth ablation (i.e. 4 yr-old, Figure 6). Habitat shift probabilities increased then rapidly reaching values ≥0.5 between ablations 15 and 19, i.e. between 8 and 15 yr-old. Thus, the median age for habitat shift (p≥0.5), estimated from Sr:Ca profiles, was equal to 11.2 yr-old (Figure 6). A similar pattern was observed in Ba:Ca profiles where 75% of the fish showed low probabilities of habitat shift (<0.2) before ablation 14^th (7-8 yr-old). In most fish, Ba:Ca profiles, habitat shift probabilities exceeded 0.5 between ablations 16 and 20, equivalent to ages between 9 and 17 yr-old. As a result, the median age for habitat shift (p≥0.5) estimated from Ba:Ca profiles was equal to 12.7 yr-old (Figure 6).

3.2. Number and Mixing Rates Between Nursery Areas

Mixtures of two sources were selected as the most informative FMD modes for both the elemental and the isotopic data sets. Second and third best models, which corresponded to mixtures of one and three sources, exhibited much lower BIC values (ΔBIC < -12.4, Table 2). Best covariance models corresponded to “VVE” and “VEV” [57] for elemental and isotopic signatures, respectively. A dominant source, contributing between 73 and 86% of the mixture, and a second minor source, contributing between 14 and 27% of the mixture, were identified by the elemental and the isotopic FMD models, respectively (Table 3). The dominant source was characterized by higher Sr:Ca, Sn:Ca and Ba:Ca, while the minor source exhibited higher values for all remaining elemental ratios, as well as for both isotopic signatures (Table 3).

3.3. Geo-Location of Nursery Areas

The predicted probabilities of occurrence for juveniles of the year (YOY) across the study area were higher (p = 0.020-0.033) within the 0-50, 50-100, and 100-150 m depth strata, particularly within three main areas: 1) west of the Antarctic Peninsula and south of 61°S, especially in the South Shetland Islands (62-63°S and 58-62°W), as well as in the Bellingshausen Sea (70-90°O), 2) in the area corresponding to the Lazarev Sea (68°S and 7°E), and 3) in the Indian Ocean sector (30°E-90°E) (Figure 7).

The mean probabilities of occurrence of juveniles were similar across the four sub-population areas depicted Söffker et al.[29]‘s Hypothesis 3 (Figure 8), ranging between 0.863 and 0.872. The geolocation of fish clustered by the FMD unknown source model showed a disperse distribution of both sources (clusters) across the Southern Ocean, rather than being restricted to some specific areas. Thus, the highest probabilities (95th percentile) of fish assigned to Source 1 were widely distributed across most of Area B (Indian Sector), in Area A1 (Lazarev Sea), and in some locations of Area C (Bellingshausen Sea) (Figure 8). The highest probabilities (95th percentile) of fish assigned to Source 2 showed a somewhat more restricted distribution, concentrated near the sampling site, in Area C, as well as in some areas of A1, A2 and B (Figure 8).

4. Discussion

By analyzing the variability in the elemental and isotopic compositions of the nuclear regions of otoliths, we evaluated here some key hypotheses regarding the stock structure of D. mawsoni in CCAMLR Area 48. To achieve this, we leveraged on two of the few recent opportunities made available to sample D. mawsoni in this region, currently closed to commercial fishing. Given these sampling limitations, and despite being based on a rather limited sample size, we believe our results are consistent enough to provide useful insights about the potential number and connectivity of nursery areas and demographic units around the South Orkney Islands. While used here to test current working hypotheses about demographic units in this CCAMLR area, we hope these new insights may be also useful for producing new hypotheses suitable for guiding further research efforts on D. mawsoni here and elsewhere.

4.1. Identification and Geo-Location of Nursery Areas

We found evidence that D. mawsoni sampled around the South Orkney Islands corresponded to a mixture of a dominant and a secondary nursery source. The number and relative contribution of these sources was consistent with Söffker et al. [29]‘s Hypothesis 2, which states that two relatively discrete sub-populations mix in this area. Nonetheless, while these authors suggest these subpopulations would be originated at Areas A1 and A2, our geo-location results failed to rule out contributions from alternative nursery grounds potentially located in neighboring subareas B (Indian Ocean Area) and C (Bellingshausen Sea).

Given methodological limitations inherent to this δ18O single-isotope geo-location approach, we were not able to resolve connectivity at finer spatial or temporal scales, including advection between distant spawning and local nursery areas. Minimum mass requirements forced us to integrate the whole first year of life (0-540 µm from the core), collapsing in a single value all chemical signatures related to maternal effects, several months of drifting larval life and first months of benthic life [7]. Since LA-ICPMS provides greater time resolution, we explore this issue by fitting FMD models to specific ablations within the first year of life, finding 3-source FMD models could be more informative than 2-source models for each of the first four ablations, which mainly represented maternal (ablation 1) and larval life (ablations 2-4) effects. This exploratory analysis (not showed in results) indicates other sub-populations may also be contributing eggs and larvae to A2 nursery grounds. A relevant finding considering simulations show eggs and larvae from the Bellingshausen Sea (sub-population C) may be advected NE, reaching the Scottish Sea, the South Orkney Islands and even the South Sandwich Islands [72]. Other simulations showed some eggs and larvae released in the Lazarev Sea (sub-population B) could be also advected into the Weddell Sea by the coastal westward return of the Weddell Gyre [72].

4.2. Ontogenetic Habitat Shift

The predicted concentration of early juveniles in coastal nursery habitats located east from the Antarctic Peninsula and south from the Weddell Sea suggest a combination of seasonal adult migrations to shallower spawning grounds, properly located to take advantage of local advective forces, such as the Weddell Gyre [29,73]. This type of retention systems has been already hypothesized for multiple species, including D. mawsoni in the Ross Sea (Sub-area 88.1), where eggs spawned at the Pacific Ridge are likely transported eastward by the Ross Gyre to shelf nursery and recruitment habitats used by most juveniles [7,25]. Nonetheless, given the prolonged duration of egg and larval stages in D. mawsoni [~ 9 months, 7], this species, as other notothenids seems more vulnerable to advection by large-scale advective forces, like the Antarctic Circumpolar Current [25].

The very distinct chemical signatures, we found, between nuclear and marginal otolith regions indicated large environmental differences between YOY and adult habitats, suggesting ontogenetic displacements to cooler and, most likely, deeper habitats [62,74–77]. This finding is consistent with available knowledge about the life-cycle of this species [2,7], described to use shallow (<150 m) shelf habitats during its first year(s) of life (<13 cm LT), moving then to use progressively deeper and wider bathymetric ranges [2,5,8]. Our analysis of Sr:Ca y Ba:Ca profiles indicated, nonetheless, that a rather discrete habitat shift would occur around 11-12 years old, which might be related to the onset of sexual maturity [3,4] and/or to the shift in buoyancy described by several authors to occur when D. mawsoni reaches a size ~100 cm LT [5,11,12,78]. This shift in buoyancy facilitates the use of a wider bathymetric range, exploring deeper layers either for food or to avoid predation by mammals [13,79].

Differences we found in δ¹³C between nuclear and marginal otolith regions were not consistent with the expected bathymetric reduction in water DIC-δ¹³C due to variability in organic matter production/oxidation ratios [80,81] and deserves further investigation about potential dietary source and metabolic effects [82–84]. While YOY might rely heavily on epipelagic prey such as Antarctic krill (Euphausia superba) and fish larvae, adults would exploit mesopelagic cephalopods, such as Teuthidae, and fishes, such as Macrouridae, Muranolepididae and Nototheniidae [85–88], all exhibiting greater trophic positions and, therefore, enriched δ¹³C signatures [89].

5. Conclusions

Overall, our results failed to support to Söffker et al. [29]’s hypothesis 2 that the South Orkney Islands stock belongs to a relative discrete demographic unit (Sub-population A2), distributed across the Weddell Sea basin. Our results show most adults had been originated elsewhere at rather distant nursery areas, mainly located south and east from these islands, on the continental shelf, including the Lazarev and the Bellingshausen sub-population (B and C, respectively) areas. This connectivity between potential YOY nursery grounds and feeding adult habitats fails also to provide any direct support to the philopatry hypotheses, although a direct testing of this would require sampling both early larvae and spawning adults.

Our results confirm the Southern Ocean and D. mawsoni otoliths exhibit enough variability to facilitate the application of otolith-based geochemical methods for investigating habitat use and connectivity patterns [76,90,91]. Nonetheless, the strength of our conclusions was likely affected by limited sample size, as well as by temporal, spatial and vertical coverage of water δ¹⁸O data. We expect, however that this work may stimulate additional efforts, clearly needed to clarify the demographic structure of this species in the Southern Ocean and, particularly, in CCAMLR area 48. Within the most relevant needs, we suggest to focus on: i) revisiting connectivity patterns after collecting otolith chemistry data from adults and, ideally, juveniles sampled at each of the four sub-population areas ii) testing philopatry and spawning site fidelity hypotheses, iii) assessing connectivity between spawning ad nursery grounds, iv) ground-truthing the spatial distribution of A2 nursery grounds and assessing its potential value for other demographic units such as A1 and B. We recommend achieving significant progress addressing these research needs before reopening the South Orkney Islands and other areas, currently closed to regular fishing operations.