Life History of the Arctic Squid Gonatus fabricii (Cephalopoda: Oegopsida) Reconstructed by Analysis of Individual Ontogenetic Stable Isotopic Trajectories

Simple Summary Gonatus fabricii is the most abundant cephalopod in the Arctic and northern part of the North Atlantic, which are areas of rapid environmental change. It is very important as a predator of many species of fish and invertebrates and as a prey for larger fish, seabirds, and marine mammals. The life cycle of G. fabricii, in particular ontogenetic changes in diet and habitat, is little known. Ecological modelling is an important method to forecast ecosystems’ dynamics in relation to a changing climate, and it is crucial to assess their present structure and functioning for a viable forecast. Here, the ontogenetic changes in the diet and habitat of large G. fabricii from West Greenland are studied using stable isotope analysis of carbon (δ13C) and nitrogen (δ15N) along trajectories in chitin beaks, which function as archival structures. Using this approach, four stages with clearly distinct ecologies were revealed in the species’ life cycle. This novel ecological periodization is a crucial baseline tool in Arctic marine ecosystem studies and cephalopod biology, and should be used in models to correctly reflect the ecological roles of G. fabricii in marine ecosystems. Abstract Cephalopods are important in Arctic marine ecosystems as predators and prey, but knowledge of their life cycles is poor. Consequently, they are under-represented in the Arctic ecosystems assessment models. One important parameter is the change in ecological role (habitat and diet) associated with individual ontogenies. Here, the life history of Gonatus fabricii, the most abundant Arctic cephalopod, is reconstructed by the analysis of individual ontogenetic trajectories of stable isotopes (δ13C and δ15N) in archival hard body structures. This approach allows the prediction of the exact mantle length (ML) and mass when the species changes its ecological role. Our results show that the life history of G. fabricii is divided into four stages, each having a distinct ecology: (1) epipelagic squid (ML < 20 mm), preying mostly on copepods; (2) epi- and occasionally mesopelagic squid (ML 20–50 mm), preying on larger crustaceans, fish, and cephalopods; (3) meso- and bathypelagic squid (ML > 50 mm), preying mainly on fish and cephalopods; and (4) non-feeding bathypelagic gelatinous females (ML > 200 mm). Existing Arctic ecosystem models do not reflect the different ecological roles of G. fabricii correctly, and the novel data provided here are a necessary baseline for Arctic ecosystem modelling and forecasting.

and other cephalopod hard structures for SIA and provide guidelines and methodological recommendations on how to apply the combination of population and individual SIA to other cephalopods.  Table S1). The samples were collected as bycatch by bottom trawl 'Alfredo 3' at depths of 384.5-1443.0 m (Table S1), and fixed in 4% formalin onboard. The squid's ML and mass were measured in the laboratory and the sex and maturity stages were identified following Golikov et al. [37]. Afterward, upper beaks were taken and their upper beak rostrum length (URL) and upper beak crest length (UCL) were measured following Clarke [56,58]. A total of fourteen individuals (six females and eight males) were used to obtain ontogenetic stable isotope trajectories. Transparent parts of the upper beaks were cut out as they have significantly different isotopic values from tanned beaks [61,69]. The upper beaks were divided in half along the crest. Sequential equal-sized subsections as close to 1 mm as possible were taken from the inner side of the crest, starting from the rostrum to the mid-point of the crest, and sequential equal-sized subsections as close to 2 mm as possible were taken starting from the midpoint of the crest to the onset of the now-cut transparent part ( Figure 2, Table S2). The higher resolution for the earlier half of ontogenesis is explained by the fact that this is where the main changes in isotopic values occurred in other cephalopods [47,51,52]. The tip of the rostrum, which is known to have not only the newest material deposited due to beak growth specifics [51,62,63], represented half of subsection 1 or less with this approach (Figure 2) and it did not bias the outcome (see Results). With 11 to 17 subsections from each upper beak, SIA was performed in 179 samples (Tables 1 and S2). The upper beak allows a sufficient number of subsections to be analysed in an ontogenetic perspective according to Guerra et al. [47] and Queirós et al. [52] (one study also exists on lower beaks, but only four individuals were analysed [66]). Thus, the upper beak provides a more direct comparison with existing studies where a sufficient number of subsections were analysed. A, lower beaks of G. fabricii are much smaller than in other cephalopods, analysed with this approach.

Materials and Methods
Data from Golikov et al. [27], where 185 beaks of G. fabricii were measured and 105 beaks underwent SIA, were used in this study. Beak URL measurements from this study and Golikov et al. [27] were pooled to update the existing equations, allowing us to estimate ML and mass of G. fabricii in West Greenland and the Arctic in general: West and East Greenland, and the Barents Sea (Table S3). Beak UCL measurements of individuals from Golikov et al. [27] were not used in the respective study, and they are used here for the first time. The difference between the isotopic composition of upper and lower beaks is significant, but it is small and is within the measurement error of the method [69]. However, this should be kept in mind when the data from this study were united with those from Golikov et al. [27] for SIBER analyses (see below). In Golikov et al. [27], the individuals were grouped as 'epipelagic', 'mesopelagic', and 'bathypelagic' based on ML and maturity stages combined: ML < 60 mm in both sexes, ML 60-130 mm in males and 60-150 mm in females, and ML > 130 mm in males and >150 mm in females. Thus, they are referred to as small, medium, and large here, respectively.

Stable Isotope Analyses
All beak subsections were dried at 60 • C for 24-48 h, weighed (to the closest value to 0.3 mg) with a micro-balance and cracked if needed, and sterile-packed in tin containers. Stable isotopic signatures were determined by a Flash EA 1112 Series elemental analyser (Thermo Scientific Inc., Waltham, MA, USA) coupled online via a Finnigan ConFlo II interface to a Delta vs. mass spectrometer (Thermo Scientific, Bremen, Germany) and expressed as: δ 13 C and δ 15 N = [(R sample /R standard ) − 1] * 1000, where R = 13 C/ 12 C and 15 N/ 14 N, respectively. The carbon and nitrogen isotope ratios were expressed in delta (δ) notation relative to Vienna-PeeDee Belemnite limestone (V-PDB) for δ 13 C and atmospheric nitrogen (AIR) for δ 15 N. Replicate measurements of internal laboratory standards (acetanilide STD: Thermo Scientific PN 338 36700) in every batch (n = 14) indicated a precision <0.2‰ for both δ 13 C and δ 15 N signatures. The mass C:N ratio for each sample, where it was obtained, ranges from 3.05 to 3.86 (Table S2). The analyses were carried out at the Marine and Environmental Science Centre, University of Coimbra (Coimbra, Portugal).

Data Analyses
The newly obtained equation to estimate ML from UCL was used to find squid ML at given subsections, and 4.2 mm was added to this value (see Beak equations in Results). Squid mass at a given subsection was estimated based on the length-mass equation from Golikov et al. [27].
A regression analysis was used to find equations fitting our data [78]. Correlations between δ 13 C and δ 15 N values were assessed using a Spearman's rank correlation [78]. At the population level, n = 118 overall and n = 57 in West Greenland were used for this analysis. At the individual level, it was possible to assess 7 of 14 individuals, as the rest had missing values for some subsections (Table S2). A Mann-Whitney U test was used to compare isotopic signatures of δ 13 C and δ 15 N and TL between two groups (e.g., sexes), and a Kruskal-Wallis H test (with a post-hoc Dunn's Z test) was used to compare among more than two groups (e.g., size groups) [78]. Differences in δ 13 C values, δ 15 N values, and TL among beak subsections were tested using a Skillings-Mack test proceeded by the Nemenyi post hoc test [79]: only the subsections common to all the beaks, i.e., until the last subsection of the beak with the lowest number of subsections overall, were tested. The packages PMCMR 4.4 [80] and Skillings.Mack 1.10 [81] in R 4.1.3 [82] were used for these analyses. The value of α = 0.05 was considered significant in this study.
Stable isotopic values were used to calculate the specialization index (s) to estimate the degree of individual specialization [83,84] based on the isotopic variation within and between individuals of a sampled population [85]. The specialization index s was defined as the variance of the isotopic signal along the crest of the given individual divided by the total variance among individuals of the sampled population plus the variance of the isotopic signal along the crest of the given individual [85]. Accordingly, s values range from 0 to 1, with 1 representing a complete overlap between the individual and the population (= extreme generalist individual) and lower s values representing higher degrees of specialization [85]. Individuals occupying < 50% of the total niche of the sampled population (i.e., s < 0.5) are considered specialists [85,86]. The specialization index s was calculated separately for each individual and for each isotope, δ 13 C and δ 15 N.
Trophic level can be estimated with fixed trophic enrichment factor (TEF), either 'classical' δ 15 N = 3.4‰ [89] or 'Arctic' δ 15 N = 3.8‰ [90] and with the standard TL equation [91] or with the scaled TEF equation [92,93] adapted for the Arctic by Linnebjerg et al. [94]. We used the standard TL equation with 'Arctic' TEF as the scaled approach was giving unrealistically high results (not included in this paper). We used Calanus finmarchicus  (Gunnerus, 1770) [95] taken from the nearest station in Hansen et al. [96] or, alternatively, in Pomerleau et al. [97], as a reference value for TL = 2.0 when a closer location was not present in Hansen et al. [96]. Interpretation of TL in the Arctic ecosystems followed the TL interpretations of other stable isotope assessments of the area [27,28,30,90,94,[98][99][100][101].
Sequential subsections of the same beak do not comply with the assumption of sample independence, and due to that, nicheROVER 1.1.0 [102] in R 4.1.3 [82] was employed to estimate isotopic niche width for each subsection and overlap among them. These analyses were only applied until subsection 14, as an insufficient number of individuals were large enough to sample subsections 15 (n = 6), 16, and 17 (both n = 2) ( Table 1). It was assumed that the mean isotopic value of all subsections of the given individual corresponded to the isotopic values of the individual obtained by the 'classical' approach, i.e., powdering the whole beak and taking a random portion. Thus, the mean isotopic values of beaks from this study (n = 14) were pooled with the dataset from Golikov et al. [27], and SIBER 2.1.6 [103] in R 4.1.3 [82] was employed to assess isotopic niche width and overlap between sexes and squid with ML 20-50 mm and ML > 50 mm. Squid with ML < 20 mm (n = 4) from West Greenland from Golikov et al. [27] could not be used in this ontogenetic niche comparison due to the low sample number. The standard ellipse area corrected for small sample sizes (SEA c ) and the Layman metric of convex hull area (TA) were estimated [103][104][105], and the Bayesian approximation of the standard ellipse area (SEA b ) was employed to compare niche width [103]. The overlap interpretation for both niche-related packages followed Langton [106], where overlap ranging from 0 to 0.29 indicated no overlap, from 0.30 to 0.60 indicated medium overlap, and from 0.61 to 1 indicated large overlap, with only the latter taken as significant [30]. Trophic levels were used instead of δ 15 N values (y axis) in niche estimations by both niche-related packages. This approach improves the ecological meaning of isotopic data when comparing individuals from different areas and ecosystems as it is a means to counter high δ 15 N baseline variation e.g., [107,108]. This approach has been repeatedly applied to cephalopods [28,30,76]. The Bayesian mixing model SIMMR 0.4.5 [109] in R 4.1.3 [82] was used to assess the relative contribution of prey to the diet of G. fabricii. Based on an overview of stomach content analyses [27], the prey groups assessed were copepods, euphasiids, shrimps (these three taxa belong to Crustacea), chaetognaths, cephalopods (G. fabricii, as it is the only squid species living permanently in the Arctic) and fish. Source isotopic values for the diet model (terminology for the diet model followed its author [109]) were taken from Pomerleau et al. [97,110], Hansen et al. [96], Agersted et al. [111], Linnebjerg et al. [94], Golikov et al. [27], Grigor et al. [101] and this study (i.e., mean values of each studied squid), and are detailed in Table S4. Crustacean and chaetognath source values were similar for all analysed upper beak subsections of G. fabricii, and source values of cephalopods and fish were taken according to squid size for ecological perception. For subsections 1-3, cephalopods and fish were combined, as their source values were not significantly different: all source values should be significantly different in at least one of the isotopes (Tables S4  and S5). Values and standard deviations of TEF were taken from the only experimental study showing isotopic values in cephalopods to change following the long-time diet composition [72]: δ 13 C = −0.20 ± 0.55‰ and δ 15 N = 3.37 ± 0.99‰. The data fitting to selected prey source values and TEFs were checked inside SIMMR 0.4.5 [109] prior to the analysis, and only the fitting values were used in the models ( Figure S1), as follows: subsections 1-3, n fitting = 23 of 31; subsections 4-5, n fitting = 16 of 25; subsections 7-9, n fitting = 20 of 36; and subsection 10, n fitting = 5 of 10. To compare the results of the models among the subsections, the χ 2 test and Fisher's exact test were used: although the latter is more adequate for small sample sizes, Fisher's exact allows the comparison of only two groups [78]. Statistical analyses were performed in R 4.1.3 [82], Statistica 10.0 (Statsoft), and PAST 4.02 [112]. Values are presented as the mean ± SE unless otherwise stated.

Stable Isotopic Values and Trophic Levels: Ontogenetic Changes
Each of the studied individuals demonstrated an ontogenetic increase of both isotopes ( Figure 3, Tables 2 and 3). Within-individual ontogenetic increase is overall lower than within-population ontogenetic increase: 0.7-2.2‰ (1.4 ± 0.13) vs. 3.0‰ in δ 13 C and 2.6-9.1‰ (5.6 ± 0.37) vs. 10.7‰ in δ 15 N, respectively. The first subsection (rostrum tip) usually had the lowest value. However, the last subsection (posterior end of the crest), except for the single analysed spent female (F6), never had the highest value. The main increase of δ 13 C values was from subsections 1 to 5 with a mean of 0.29‰ per subsection, which decreased to a mean of 0.01‰ per subsection afterward ( Figure 3, Tables 2 and S2). The lowest δ 13 C values were from subsection 1 (mode among individuals) or 2 (mean among individuals), and the highest was from subsection 6 (mode) or 8 (mean) ( Table 1,  Tables 2 and S2). The main increase of δ 15 N values was from subsection 2 to 7 with a mean of 0.89‰ per subsection and decreased to a mean of 0.08‰ per subsection afterward ( Figure 3, Tables 2 and S2). The lowest δ 15 N values were from subsection 1 (mode) or 2 (mean), and the highest was from subsection 10 (both mode and mean) (Tables 1, 2 and S2).
A significant positive correlation between δ 13 C and δ 15 N values existed in six of the seven individuals we were able to check (Table S6). A significant positive correlation between δ 13 C and δ 15 N values also existed at the population level (Table S6).

Variation of Stable Isotopic Trajectories: Among Individuals and between Sexes
The specialization index of δ 15 N suggested that the studied G. fabricii were all generalists with low variability among the individuals (Table 1). Significant differences were present between the sexes: χ 2 = 17.56, d. f. = 5, p < 0.0001. The specialization index of δ 13 C suggested that some of the studied individuals of G. fabricii were specialists in habitat usage, with high variability among the individuals ( Table 1). The specialization index of δ 13 C was lower than 0.5 in three of six females and in one of eight males, and the latter was the lowest among the studied squid (s = 0.21). Significant differences were present between the sexes: χ 2 = 92.21, d. f. = 5, p < 0.0001.
Both sexes reached the highest values of δ 15 N (and thus TL) in subsection 10. However, females reached the highest values of δ 13 C at subsection 9 (larger ML), while males at subsection 8 (smaller ML). Within-individual ontogenetic increase was smaller than withinpopulation for both sexes. No significant differences were found between the sexes for within-individual increase: n = 13, U = 17.5, p = 0.67 smallest to largest and n = 13, U = 20.0, p = 0.95 first to last in δ 13 C; and n = 14, U = 23.0, p = 0.95 smallest to largest and n = 14, U = 22.0, p = 0.85 first to last in δ 15 N.

Isotopic Niches and Diet Models as a Basis of Life History Reconstruction
No clear pattern of ontogenetic changes in niche width was found when analysed by subsections. The widest niches were from subsections 3, 4, and 6, and the first and the last analysed subsections were of roughly equal width (Figure 4a,b, Table S7). Niches of subsections 1 and 2 had the least overlap with the rest (Figures 4b and S2, Table 4).
A significant positive correlation between δ 13 C and δ 15 N values existed in six of the seven individuals we were able to check (Table S6). A significant positive correlation between δ 13 C and δ 15 N values also existed at the population level (Table S6).      Based on nicheROVER, SIA values, and diet modelling analyses, population niches with SIBER should be compared among squid with ML < 20 mm (see Material and Methods), ML 20-50 mm, and ML > 50 mm. The niche of squid with ML > 50 mm was significantly wider than of squid with ML 20-50 mm, with low overlap between the niches (Figure 5a, Table 5). No significant differences were found between the isotopic niche width of males and females (Figure 5c,d, Table 5). These niches had significant overlap for both size classes (Figure 5c,d, Tables 5 and S8). The pattern of ontogenetic niche increase and overlap in females and males was similar (when the sexes were analysed separately) ( Table S8).
The diet corresponding to beak subsections 1-3, 4-5, 7-9, and 10 was modelled and compared. Chaetognaths were always a minor component of G. fabricii's predicted diet ( Figure 6, Table 6). Crustaceans constituted a half or more in subsections 1-3 and 4-5, or slightly less than a half in subsections 7-9 and 10 ( Figure 6, Table 6). Combined, cephalopods, and fish formed a complementary pattern to crustaceans, i.e., they constituted slightly less than a half in subsections 1-3 and 4-5, and slightly more than a half in subsections 7-9 and 10 ( Figure 6, Table 6). No significant differences in the predicted relative contribution of prey taxa (i.e., chaetognaths, crustaceans, and cephalopods + fish) were found among different subsections (Table S9): the described patterns are trends only. Cannibalism is predicted to be very important in the diet of G. fabricii ( Figure 6, Table 6). Within crustaceans, the smallest analysed squid (subsections 1-3) were predicted to consume copepods, euphausiids, and shrimps in roughly equal ratios, which significantly differed from other size groups of G. fabricii (Table S9) who mostly relied on larger crustaceans, i.e., euphausiids and shrimps, than copepods ( Figure 6, Table 6).

Beak Equations
The equations for ML and UCL were described as: ML = 3.25 × UCL 1.40 (n = 86, R 2 = 0.95, p < 0.0001) in West Greenland, and ML = 3.77 × UCL 1.34 (n = 142, R 2 = 0.93, p < 0.0001) overall in the Arctic (West and East Greenland, and the Barents Sea). The difference between real and predicted ML was 4.2 ± 2.44 mm for the West Greenland equation and 3.3 ± 1.72 mm for the overall Arctic equation. Updated URL equations (Table S3) outperformed previous ones cf. [27]: estimated ML was 5.1 ± 2.02 mm smaller than real ML vs. 12.2 ± 2.36 mm previously for the overall Arctic equation, and 6.1 ± 2.92 mm vs. 16.5 ± 3.51 mm, respectively, for the West Greenland equation.      Table 6). No significant differences in the predicted relative contribution of prey taxa (i.e., chaetognaths, crustaceans, and cephalopods + fish) were found among different subsections (Table S9): the described patterns are trends only. Cannibalism is predicted to be very important in the diet of G. fabricii (Figure 6, Table 6). Within crustaceans, the smallest analysed squid (subsections 1-3) were predicted to consume copepods, euphausiids, and shrimps in roughly equal ratios, which significantly differed from other size groups of G. fabricii (Table S9) who mostly relied on larger crustaceans, i.e., euphausiids and shrimps, than copepods ( Figure 6, Table 6).

Ontogenetic Change in Isotopic Values, and Its Ecological Applications
Significant ontogenetic increase of both carbon and nitrogen isotopic values (and thus TL) in G. fabricii is now confirmed at the individual level. Previous SIA studies of G. fabricii [27,74], performed on whole beaks [27] or pieces of soft tissues [77] of multiple individuals (population analysis, not individual analysis), established this significant ontogenetic increase in populations. However, they did not find the causes for this increase and did not show at which sizes the main increase happens and to what degree changes in an individual follow that of a population. The major issue arising from this is that we need to apply a combination of individual and population SIA to study a species' ecology more completely and from different perspectives (see practical recommendations for such an approach below).
The main ontogenetic increase in δ 13 C values takes place in squid from ML 7.0-7.4 mm to ML 30.9-44.2 mm. The increase in δ 13 C values is larger than expected based on dietary changes, i.e., mean 1.4‰ vs. about 1.0‰ [89]. The main ontogenetic increase in δ 15 N values (and TL) takes place in squid from ML 7.0-7.4 mm to ML 47.0-53.8 mm. The increase spans a mean of 1.5 TLs among studied individuals. Shifts in food spectra are the reason for changes in δ 15 N values and TL (as also highlighted by the diet models, below). However, other reasons such as individual movements may partly explain changes in δ 13 C values. Three types of migrations are known for G. fabricii. The first involves horizontal migration, the second is ontogenetic downward migration, and the third is diurnal vertical migration [27,31,32,34,[36][37][38][39][40][41][42]113]. There is an occasional presence of squid with ML 20-50 mm deeper than 200 m [36,38,[113][114][115][116], while they are supposed to exclusively inhabit epipelagic layers [42], which means the migrations outlined above also have a role in ontogenetic increase of δ 13 C values. It is clear that squid with ML 30.9-44.2 mm can occasionally exploit the whole water column when we combine these literature data on depth distribution with δ 13 C values reported here and in Golikov et al. [27].
The ML values highlighted by δ 15 N and TL analyses are used to group the squid and assess the species' assimilated diet. Reported results on diet are generally in accordance with previous studies using stomach contents [36,[43][44][45][46]115,116] that showed crustaceans are major prey for G. fabricii with the increasing importance of fish and cephalopods in individuals larger than ML 54-70 mm [44][45][46]. However, assimilated diet predicted by SIA data often does not entirely coincide with stomach contents analysis results in cephalopods [28,30,117]. Moreover, food spectra changes in G. fabricii were expected to happen due to the hook appearance at its arms and tentacular clubs [43][44][45][46], which exactly fits the predictions from this study. Gonatus fabricii already reaches the highest TL (adult diet) at ML 76.9-108.2 mm, with the maximum ML recorded for G. fabricii being 389 mm [46]. Such behavior is also detected in Berryteuthis magister (Berry, 1913) [118], a gonatid squid from the North Pacific that does not significantly change the size of the fish it preys upon throughout ontogenesis [48]. Dosidicus gigas (d'Orbigny [in 1834-1847], 1835) [119], a large squid from the Pacific, is also known to rely on small prey, despite that it is occasionally hunting for large prey [120,121]. Cannibalism is well-known for G. fabricii [36,44,46,115] and other gonatids [48,122] and the predicted diet models from this study underline its importance. Previous studies on stomach contents did not estimate prey sizes of G. fabricii, but the Pacific gonatids are known to eat fish and squids of 15-150% of their ML [48,122].

Variation of Stable Isotopic Trajectories: Among Individuals and between Sexes
Substantial variation among individual ontogenetic stable isotopic trajectories is found. We applied the specialization index s to explore it further. This is the first application of s to stable isotopic data from marine invertebrates to the authors' knowledge. Specialization index s of δ 15 N shows that G. fabricii is a generalist predator, as most predatory squids are e.g., [48,50,[120][121][122][123]. Significant differences in s of δ 15 N between the sexes in G. fabricii are influenced by only one individual: s values were 0.95-0.98 in all individuals, except female F5 with s = 0.83. Specialization index s of δ 13 C shows that G. fabricii is a specialist in habitat usage. In turn, it suggests variability in migrations among individuals. Females have much higher variation in habitat usage than males, resulting in lower specialization index s in δ 13 C, and females reach the highest δ 13 C values at larger sizes than males. Thus, at least part of the variation in δ 13 C is explained by these differences between the sexes. The greater migratory capability of females is known in oceanic squids Ommastrephes bartramii (Lesueur, 1821) [49,124] and D. gigas [123]. The niche width in G. fabricii is similar in females and males, with significant niche overlap between the sexes. This means that high individual variation in habitat usage and different migratory capabilities between sexes are not enough to differentiate niches between females and males due to individual squid's highly generalist diets. These results contrast with previous studies on other cephalopods where significant differences in niche width were found between sexes [30,49,123,125].

Life Cycle of Gonatus fabricii
Previous attempts to differentiate between ecological roles during the ontogenesis of G. fabricii mostly used maturity stages ( Table 7). The maturity stages, however, are not a good proxy for an ecological role as there is a strong overlap in ML of squid in different maturity stages (Table 7). An overview of ontogenetic changes in body shape, mass, and reproductive indices [36,39,40,42] suggests a rather slow pace of changes in ecological role, perhaps due to slow growth. However, these characters alone cannot be relied upon, as shown above. A previous study that applied SIA in G. fabricii [27] showed that the changes in the ecological role might happen faster than previously supposed [42]. Clearly separate ecological roles, which G fabricii occupies during ontogenesis, are established in this study (Table 7). This is achieved by a combination of individual and population SIA. The life history of G. fabricii is as follows: (1) epipelagic squid with ML < 20 mm prey mostly on crustaceans, especially copepods (as per stomach contents analysis) or copepods and larger crustaceans in equal proportions (as per SIA prediction); (2) squid with ML 20-50 mm mostly stay in the epipelagic layer, but occasionally exploit the whole water column, with their food spectra shifting to larger Crustacea, fish, and cephalopods; (3) squid with ML > 50 mm mostly prey on fish and cephalopods with a high proportion of cannibalism, and stay in the meso-and bathypelagic layers (occasionally might be recorded in the epipelagic layers, as they have daily vertical migrations); and 4) non-feeding and almost non-moving bathypelagic gelatinous females with ML > 200 mm (Table 7).

Ecosystem Implications
Dividing the life history of G. fabricii into distinct ecological roles provides baseline data on cephalopod biology and in Arctic marine ecosystem studies. Active cephalopod predators are often compared with bony fish in terms of ecological role. However, the life histories of these taxa are very different [130][131][132]. When based on a sufficient number of beak subsections, individual life history SIA reconstructions of cephalopods [47,52,66] strongly suggest that cephalopods change their ecological roles earlier in their ontogenesis compared to fish cf. [133][134][135]. The growth speed of juvenile G. fabricii in West Greenland was estimated to be 0.13-0.18 mm day −1 [136] and 0.27 mm day −1 [38]. This suggests that squid remains in a 'pure' epipelagic stage (ML < 20 mm) for 2.5-4.2 months, then transition their food spectra and habitat (ML 20-50 mm) for 3.7-6.3 months, with overall life expectancy in West Greenland of up to three years [36]. Thus, longer-living cold-water Arctic squid, which spend much of their life in the relatively stable deep-sea environment, still have major ecological changes during early ontogenesis, cf. the large Antarctic squid and octopus [51,52,66].
A high TL of 5.0 or above is only known in polar squids from the Arctic [27], and this study and Antarctic [82,137,138]. Very high δ 15 N values in the beaks of Histioteuthis atlantica (Hoyle 1885) [139] from Subantarctic waters suggest comparably high TL, as other squids in the area showed much lower δ 15 N values in their beaks [140]. The high TL of G. antarcticus Lönnberg, 1898 [141] is most likely explained in a similar way as G. fabricii (see above).

Modelling and Analyses Implications and Recommendations
Previous studies used eye lenses and gladii of cephalopods to analyse individual ontogenetic stable isotopic trajectories; however, they failed to identify the fast pace of changes in cephalopod ecological role early in ontogenesis [48][49][50]123,[142][143][144][145][146][147][148]. Unlike beaks, neither eye lenses nor gladii were proven to have daily growth increment deposition [62][63][64][65]. Considering beak utility when necessary methodological considerations were performed for the anterior-most and posterior-most parts [51,[61][62][63]69], beaks are recommended for ontogenetic SIA studies in the future. Analysis of individual ontogenetic stable isotopic trajectories from beak crests provides precise sizes of ecological role changes, especially if coupled with analysis of multiple individuals with their whole beaks powdered. The recommendation is to use one hundred beaks in equal proportions of the small, medium, and large sizes for population SIA, which is a large enough sample size [27,76]. With each size group being n > 30, it would not bias the SIBER niche analysis outcome [149]. At least ten large beaks are recommended for individual ontogenetic SIA, taking consecutive pieces along the crest. Combining these approaches would give an extensive overview of the species' life history, together with its reconstruction with exact sizes, if a beak equation to estimate squid size exists for the study species. Equations to estimate ML of G. fabricii from URL and UHL exist for the Barents Sea, West Greenland, East Greenland, and the overall Arctic [27], and this study, and equations to estimate squid ML from UCL only exist for West Greenland and the overall Arctic (this study). Our recommendation is to use the equation for the specific area of the Arctic (i.e., the Barents Sea, West Greenland, or East Greenland) if it is known where the beaks come from, and to use the overall Arctic equation in the absence of such information.
For future Arctic ecosystem modelling we recommend using the size-specific ecological role periodization established here. The ecological role of G. fabricii is misrepresented in current models. The dataset of Planque et al. [13] is used as a basis for most of the Barents Sea ecosystem models [14,15,[17][18][19][20]. In these, G. fabricii only has 12 prey taxa (all Crustacea and G. fabricii itself, i.e., no fish) and 23 predators, cf. at least 49 and 47, respectively [27]. Thus, these models completely miss the role of G. fabricii as a fish consumer. The older models [10][11][12] and those reusing them as a baseline [16] also do not reflect the ecological role of G. fabricii. Dommasnes et al. [10], Blanchard et al. [12], and Murphy et al. [16] are missing some groups of prey and predators: in particular, they miss cannibalism in G. fabricii. Moreover, they extrapolate biomass values of G. fabricii from the Nordic Seas to the Barents Sea [10,12,16], which has been already shown to be incorrect elsewhere [34]. The Nordic Basin and the Barents Sea are very different in their oceanographic conditions, and thus in their role in the range of G. fabricii [34,37]. For West Greenland, Pedersen and Zeller [11] provided a better description of the trophic role of G. fabricii, however, suffer from biomass estimation and P/B and Q/B coefficients for the species.
The only G. fabricii biomass assessment where large individuals were included was performed by Bjørke and Gjøsaeter [31] for the Nordic Seas, which largely covers the Norwegian Sea. It shows 1.5 million t of squid in the upper 30 m layer and 8 million t of meso-and bathypelagic squid [31]. Applying the ecological stages established here might prove difficult, as: (1) it is unknown how much 'epipelagic' squid are with ML < 20 mm and how much belong to other groups, and to which groups; (2) the layer between the upper 30 m and mesopelagic layer is missing; (3) it is unknown how much of meso-and bathypelagic squid are with ML < 50 mm and ML > 50 mm; and 4) the proportion of gelatinous females in the latter group is also unknown. Extrapolation based on Figure 3 from Bjørke and Gjøsaeter [31] suggests the whole estimated 'epipelagic' of 1.5 million t are squid with ML < 20 mm, then there are 2.5 million t of squid with ML of 20-50 mm. If the sex ratio is assumed to be equal (it has never been assessed), and surviving to spawn is 5% [31], then the biomass of gelatinous females would be 0.25 million t. The biomass proportion of G. fabricii with a different ecology in the Nordic Seas would be approximately 12% squid with ML < 20 mm, 21% squid with ML 20-50 mm, 67% squid with ML > 50 mm, and 0.02% gelatinous females. However, G. fabricii is supposed to live up to two years in the Nordic Seas [150], and this is the assumption Bjørke and Gjøsaeter [31] relied on. In this study, G. fabricii from West Greenland is assessed to define the ontogenetic stages with different ecology, which is supposed to live up to three years [36].

Conclusions
We detect a significant ontogenetic increase of δ 13 C and δ 15 N values (and TL) in G. fabricii in individual squid. There is substantial variation among stable isotope trajectories in individual G. fabricii beaks, with δ 13 C values having higher variation among individuals, especially in females, than δ 15 N values. At the same time, isotopic niches are not different between the sexes, and are highly overlapping. Exact sizes when the species occupies different ecological roles, both ML and mass, are found. We recommend that this ecological periodization (with four clearly distinct stages) be used in future Arctic marine ecosystem modelling. Currently, all existing Arctic ecosystem models fail to reflect the ecological role of G. fabricii correctly, despite its immense importance in Arctic ecosystems. From supposedly three years of squid ontogenesis in West Greenland cf. [36], stage 1 is expected to last for 2.5-4.2 months followed by 3.7-6.3 months stage 2. Thus, the main changes in Arctic squid ecology occur during the first 17-29% of the supposed life cycle and are coupled with a fast size increase. We find that beaks have a superior capacity to visualize ontogenetic changes in squid ecology compared to other hard structures such as eye lenses and gladii. To assess the variation in stable isotope data, the specialization index s is successfully applied for the first time in cephalopods. Finally, an equation to estimate the ML of G. fabricii from UCL is provided, and existing equations to estimate ML and mass from URL are updated.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ani12243548/s1, Table S1. All individuals of Gonatus fabricii, analysed for this study (n = 29): location, data, depth, sex, maturity stage, mantle length and upper beak measurements: rostrum length and crest length. ML: mantle length, URL: upper beak rostrum length, UCL: upper beak crest length, SIA: stable isotope analysis, n/a: no results, -: not applicable. Table S2. Data on Gonatus fabricii, analysed by stable isotope analysis (n = 14), and each of the analysed upper beak subsections (n = 179): mantle length, sex, maturity stage, number and size of subsections, mass C:N ratio, values of δ 13 C and δ 15 N, estimated trophic level, and estimated mantle length and mass at given beak size. ML: mantle length, TL: trophic level, n/a: no results, -: not applicable. Table S3. Equations to estimate mantle length and mass of Gonatus fabricii from upper beak rostrum length. Significant p-values are in bold. ML: mantle length, URL: upper beak rostrum length, n: sample size. Table S4. Values of δ 13 C and δ 15 N for the prey group sources used in Bayesian mixing model SIMMR 0.4.5. Values are mean ± SD. N: sample size, -: not applicable. Table S5. Differences in values of δ 13 C and δ 15 N for the prey group sources used in Bayesian mixing model SIMMR 0.4.5. Kruskal-Wallis H and Dunn's Z tests are provided in the table. Significant p-values are in bold. -: not applicable. Table S6. Spearman's rank correlation between δ 13 C and δ 15 N values in Gonatus fabricii at population and individual levels. Significant p-values are in bold. n: sample size. Table S7. Isotopic niche width, estimated in nicheROVER 1.1.0 for different upper beak subsections of Gonatus fabricii. Values are mean ± SD. n: sample size. Table S8. Isotopic niches metrics (TA, SEA c , and SEA b ) and overlap in Gonatus fabricii across sexes, sizes, and studied areas, estimated in SIBER 2.1.6. SEA b values are means ± SD. Significant p-values are in bold. Significant overlap values are in bold. ML: mantle length, n: sample size, -: not applicable. Table S9. Differences in the predicted relative contribution of prey to the diet among different upper beak subsections of Gonatus fabricii, estimated in Bayesian mixing model SIMMR 0.4.5. χ 2 test and Fisher's exact test are provided in the table, above and below the diagonal, respectively. Significant p-values are in bold. -: not applicable. Figure S1. Checking the data on Gonatus fabricii fitting to the prey group sources and trophic enrichment factors in SIMMR 0.4.5. Consumers (studied individuals; violet empty circles) and mean source signatures (see in-figure legend for reference) are shown. The exact values of sources and their standard deviations are in Table S3, and trophic enrichment factors and their standard deviations are in Materials and methods. Figure