Variability in nitrogen-derived trophic levels of Arctic marine biota

Stable isotopes are often used to provide an indication of the trophic level (TL) of species. TLs may be derived by using food-web-specific enrichment factors in combination with a representative baseline species. It is challenging to sample stable isotopes for all species, regions and seasons in Arctic ecosystems, e.g. because of practical constraints. Species-specific TLs derived from a single region may be used as a proxy for TLs for the Arctic as a whole. However, its suitability is hampered by incomplete knowledge on the variation in TLs. We quantified variation in TLs of Arctic species by collating data on stable isotopes across the Arctic, including corresponding fractionation factors and baseline species. These were used to generate TL distributions for species in both pelagic and benthic food webs for four Arctic areas, which were then used to determine intra-sample, intra-study, intra-region and inter-region variation in TLs. Considerable variation in TLs of species between areas was observed. This is likely due to differences in parameter choice in estimating TLs (e.g. choice of baseline species) and seasonal, temporal and spatial influences. TLs between regions were higher than the variance observed within regions, studies or samples. This implies that TLs derived within one region may not be suitable as a proxy for the Arctic as a whole. The TL distributions derived in this study may be useful in bioaccumulation and climate change studies, as these provide insight in the variability of trophic levels of Arctic species.


Introduction
Food web and food chain models are generally constructed to evaluate the trophic structure and dynamics of ecological communities, typically serving one of three objectives: (1) defining consistent trophic links and patterns within ecological communities, (2) determining factors affecting or structuring these communities and (3) determining flow pathways of contaminants, nutrients and energy through communities and ecosystems (Vander Zanden and Rasmussen 1996). Studies focusing on these objectives predominantly use the trophic level (TL) concept to group species into discrete levels to describe the organism's position in a food web (Lindeman 1942). However, the use of discrete trophic levels fails to acknowledge the complex interactions between organisms within a food web (Vander Zanden and Rasmussen 1996). Furthermore, use of an average TL per species is controversial as this ignores the variability induced by differences in feeding patterns between individuals within one species (Reed et al. 2016). This is especially the case in Arctic marine systems which exhibit high spatiotemporal intra-species variability in trophic level, e.g. driven by seasonal fluctuations in light and temperature (de Laender et al. 2009). High peaks in abundance of primary producers and declining prey availability due to loss of sea ice in summer lead to changing trophic interactions among species in highlatitude marine environments, affecting the trophic position of species (Murphy et al. 2016;Durant et al. 2019). These fluctuations affect the species' lipid content and body size, as well as the community composition and related prey availability throughout the year (de Laender et al. 2009;Falk-Petersen et al. 2009;Blanchard et al. 2017).
Traditionally, the trophic level of a species is quantified based on stomach content analysis, which may be hampered by difficulties in identifying organic material and inter-and intra-species differences in digestion rate (Vinagre et al. 2012). Furthermore, stomach content analysis only provides information on what is eaten within a limited time frame and may not be representative of the species' diet over a longer period (Kiljunen et al. 2006;Vinagre et al. 2012). In contrast, stable isotopic techniques provide a time-weighted measure of trophic level (Post 2002;Carscallen et al. 2012). The isotopic nitrogen signature (δ 15 N) of a consumer is systematically enriched by 2-4% relative to the consumer's diet (Post 2002), making δ 15 N a reliable proxy to describe flows of contaminants and energy through food chains (Walters et al. 2008). Recent studies support the use of a combined approach where the trophic level is inferred from δ 15 N analysis by using a combination of species-specific and tissue-specific enrichment factors and an ecosystem-specific baseline species. This allows for ecosystem, tissue-and species-specific fractionation, giving a more realistic view of the species' diet over time than when looking at discrete "actual" trophic levels (de Laender et al. 2009). Although δ 15 N is considered to accurately describe the trophic level of a species within one region and season, it is considered less accurate when comparing between multiple systems (Vander Zanden et al. 1999) because of the uncertainties introduced by variability in environmental conditions and in feeding habits of species between regions.
Variation in the trophic level of a species-as estimated based on stable nitrogen isotopes-may result from different sources, propagating through different levels of aggregation. These include differences in and between individual organisms, differences in environmental conditions, methodological differences (e.g., assumed baseline species) and analytical errors. Examples include differences in nutrient uptake (notably nitrate) by phytoplankton as a result of differing abiotic circumstances (in, e.g. sea water temperature and sea ice coverage) between Arctic regions may alter the isotopic signature of baseline species used in estimating trophic positions (Tamelander et al. 2009). Moreover, variation may occur in the systematically enrichment of δ 15 N, depending on the sample tissue, because of tissue turnover rates and metabolic discrimination (Clark et al. 2019). In the present study, we aim to investigate and quantify these sources of variation in the trophic level of Arctic species by aggregating our data across four levels in a 'variance model'. By determining the variance in TLs on an intrasample, intra-study, intra-region and inter-region level, we aim to evaluate the utility of using TL estimates for a single region as a proxy for the Arctic as a whole. As the Arctic food-web structure correlates strongly with environmental gradients (Kortsch et al. 2019), we expect variation in TLs at the inter-region level to be higher than variation at the intrasample, intra-study and intra-region level. To evaluate variability in trophic levels of Arctic species, we collated data on trophic levels, stable nitrogen isotopes, baseline species and fractionation parameters from different studies on multiple Arctic species in both pelagic and benthic food webs.

Data collection and compilation
Data on stable nitrogen isotopes for Arctic species, and corresponding TL parameters (see Eq. 1; δ 15 N Baseline , TL baseline , Δ15N, ΔD), were collected by conducting an extensive literature search using the Web of Science (Clarivate Analytics 2020). We combined search strings related to stable isotope analysis (e.g. 'stable isotope analysis' and 'nitrogen stable isotopes') and biota in the European, Canadian and Alaskan Arctic using general terms (e.g. 'Arctic biota') as well as species names (e.g. 'Ursus maritimus' and 'Calanus hyperboreus'). Stable isotope data were either extracted from tables or manually digitized using DigitizeIt (http://www. Digit izeIt .de/). Only stable isotope data sampled from April until October and after the year 2000 were included in the dataset, due to a lack of data outside this timeframe. Distinction was made between benthic and pelagic food webs. Although additional organism-specific and sample-specific parameters (i.e. age, length, body weight, date and sampling tissue) were included in the database, no further sub-setting was based on these parameters. The initial search resulted in 65 useful articles and reports, encompassing 148 species, covering four distinct Arctic areas: Alaskan Beaufort Sea, Canadian Beaufort Sea, Canadian Archipelago and Svalbard ( Fig. 1; see SI for the full reference list). Data pertaining to unique species (i.e. only observed in one of the four areas) were disregarded, resulting in a dataset comprising 107 species (29 pelagic, 78 benthic species), covering approximately 2400 individual records (see TL dataset.csv in the Supporting Information for raw data (ESM_1.csv)).

Trophic level estimation and quantifying trophic variation among four Arctic areas
For each record, nitrogen-derived continuous trophic levels were calculated based on nitrogen stable isotope data collected for Arctic biota (δ 15 N), according to the following general model (Post 2002;Jardine et al. 2006): where TL consumer represents the tropic level of the consumer, TL baseline represents the trophic level of the baseline species for a particular region and study (primary consumer, typically C. hyperboreus), δ 15 N consumer the measured stable nitrogen isotope of the consumer, δ 15 N baseline the measured stable nitrogen isotope of the baseline species, and ∆15N represents the nitrogen enrichment (or fractionation) factor of δ 15 N across trophic levels. In some cases, the study included an additional diet-tissue isotopic fractionation factor ( ΔD ), correcting for differences in stable nitrogen isotopic values across tissue types (Hobson and Clark 1992). In a minority of the studies, two baseline species were used in calculating the trophic level of a species. In this case, δ 15 N baseline was calculated according to Post (2002): (1) where δ 15 N baseline, 1 and δ 15 N baseline, 2 represent the measured stable nitrogen isotope of the two baseline species and α represents the reported proportion of the two food sources. If only one baseline species was reported, an α of 1 was assumed. We used the baseline species, the corresponding δ 15 N baseline values, and the enrichment factors as reported in the associated studies. When comparing between regions, to correct for choice of baseline species, we normalized the average baseline δ 15 N, standardizing for particulate organic matter (TL = 1) using Eq. 1.

Aggregation of trophic level estimates
In order to compare trophic level estimates between regions, the estimates were aggregated at the level of samples, studies and regions, respectively ( Fig. 2). At the sample level, (2) δ 15 N baseline = ⋅ δ 15 N baseline,1 + (1 − ) ⋅ δ 15 N baseline,2 , Fig. 1 Sampling regions of stable nitrogen isotopes in Arctic biota, obtained from an extensive literature search for the Svalbard Archipelago, the Canadian Archipelago, the Canadian Beaufort Sea and the Alaskan Beaufort Sea/Chukchi Sea. The size of each point reflects the relative number of samples at that location variation due to multiple replicates was quantified by means of Monte Carlo simulation (1000 iterations), i.e. parameters of Eq. 1 for which a mean value and a standard deviation were reported based on multiple replicates (e.g. δ 15 N consumer , δ 15 N baseline and ∆ 15 N) were replaced by normal distributions. The average variance of studies with multiple replicates was used as an indicator of intra-sample variability. At the study and region level, the weighted arithmetic mean and variance of the available average TL values were calculated, using the sample size of the replicates (study level) and samples (region level) as weights. Before aggregation, normality of the average TL values was checked using the Shapiro-Wilk test, and Levene's test was performed to verify equal variances. For each species and region, the aggregated mean TL and its standard deviation were then used to derive the 2.5th, 50th and 97.5th percentiles of the TL distributions. To account for bias due to differences in species composition between regions, we then compared TL averages based on common species between the four regions by performing an one-way ANOVA.

Uncertainty and variability in trophic level estimates
Variation at lower aggregation levels propagates into higher aggregation levels ( where σ measured,i 2 is the variance measured at aggregation level i, VAR i is the unbiased variation at level i, measured,i−1 2 is the average measured variance at aggregation level i − 1 and N i−1 is the sample size at level i − 1. This equation was rearranged to obtain the unbiased variation at each aggregation level: By comparing the average (unbiased) variation calculated at the intra-sample, intra-study and intra-region level with the (unbiased) variation among the four regions (interregion), we evaluated the magnitude of different sources of variation and uncertainty. By comparing region TL medians, we could evaluate the utility of using TL estimates from a single region as a proxy for the Arctic as a whole. Additionally, the coefficient of variation (CV) of each TL distribution was calculated to assess the dispersion of species' trophic level relative to the mean observed TL. All analyses and simulations were performed in R statistics v3.4.1 (see the Supporting Information for the code (ESM_2.txt)).

Identifying main external drivers of variation in trophic level
Principal component analysis (PCA) was used to explore patterns between estimated TL values on the one hand and TL parameters (i.e. δ 15 N baseline , TL baseline and Δ15N) and extraneous factors (δ 13 C, sex, age, longitude, latitude, sampling month and year) on the other, using the function prcomp in R statistics v3.4.1. Since relative importance is important in PCA analysis, variables were normalized by min-max scaling. The species nitrogen isotopic values (δ15N) were excluded as a possible influential factor in the principal component analysis, as these were directly used in estimating TL estimates. The most influential predictors of trophic level were then identified by standardizing regression coefficients using the Relaimpo package R statistics 3.4.1.

Nitrogen isotopes and fractionation factors measurements
Our database encompassed 107 species, including species from benthic (N = 78) and pelagic (N = 29) food webs. Nitrogen isotopes varied considerably among taxa and sample regions in both food-web types (Fig. 3 graph a and b). Mean Fig. 2 Species-specific TL estimates were aggregated at different levels, i.e. replicates (R) were aggregated at the sample level, samples (S) were aggregated at the study level, and studies (St) were aggregated at the region level. Finally, species-specific TL estimates at the region level (Rg) were compared Fig. 3 Distribution of stable nitrogen isotopes (a and b) and TLs (c and d) and corresponding 95% confidence intervals, for species within the European Arctic (green), the Canadian Beaufort Sea (in blue), the Canadian Archipelago (in red) and Alaskan Arctic (black) for pelagic food webs (left graphs) and benthic food webs (right graphs). Levels of significance when comparing the trophic level of species between regions: ***P < 0.005, **P < 0.01, *P < 0.05, n.s. no significant difference stable nitrogen isotopes (δ 15 N) per species ranged from 4.8‰ in pelagic particulate organic matter (pelagic-POM) to 20.1‰ in polar bears (U. maritimus) in pelagic food webs and from 3.9‰ (benthic-POM) to 15.2‰ for the eelpout Lycodes polaris in benthic food webs. When looking at inter-regional variation, benthic species sampled at Svalbard on average had significantly lower δ 15 N values than the same species sampled in Alaska, the Canadian Beaufort Sea and Canadian Archipelago (P = 0.004, F = 7.62, LSD-test). The same holds for pelagic food webs, in which species sampled at Svalbard had significantly lower nitrogen isotopes than their North American counterparts (P < 0.0001, F = 92.13, LSD-test). Table 1 summarizes the stable nitrogen isotope fractionation factors (∆15N) reported in the literature. For pelagic food webs, the fractionation factors differed significantly between regions (P < 0.001, One-way ANOVA/LSD), with highest values reported for Alaska (3.74 ± 0.14, see Table 1) and lowest for the Canadian Archipelago (2.80 ± 0.48). In benthic food webs, ∆15N values did not differ significantly between regions (P = 0.145, One-way ANOVA), with an average fractionation factor of 3.43% (± 0.1).

Baseline species
The majority of TL baseline s reported in the literature were based on Calanus sp. (TL = 2), accounting for 71.5% and 77.7% of all baselines in both benthic and pelagic food webs, respectively (See Table S1 for a complete list of baseline species per study and region). The average trophic level of baseline species reported for pelagic food webs in Alaska was relatively high (3.38 ± 0.98), which can be explained by the use of Phoca hispida as a baseline species by Rogers et al. (2015). Omitting these data records would result in an average TL baseline for Alaskan pelagic food webs of 1.88 ± 0.5. In benthic food webs, TL baselines used in Alaska and Svalbard showed to be significantly lower than those used in studies carried out in the Canadian Archipelago and Canadian Beaufort Sea (P < 0.0001, LSD-test). The average baseline δ 15 N, when standardized for particulate organic matter (TL = 1), among all regions for benthic food webs was 4.94 ± 1.9‰ (Table 1), with the lowest average baseline δ 15 N reported in Svalbard (4.2 ± 1.14‰) and the highest average baseline δ 15 N reported in the Canadian Archipelago (5.6 ± 1.64‰). Significant differences in baseline δ 15 N, standardized for TL = 1, were observed between Svalbard and all other regions in both the pelagic and benthic food web (P < 0.0001, LSD-test, Figure S2). In pelagic food webs, the average standardized baseline δ 15 N among all regions was 5.61 ± 1.82‰, with the lowest average baseline δ 15 N again determined for Svalbard (4.94 ± 2.59) and the highest baseline reported for Alaska (5.85 ± 1.18; Table 1 and Figure  S2). In benthic food webs, the average standardized baseline Table 1 Description of TL parameters used in the present study, including region, the number of studies included, the total number of data records included per region, the average (± S.E.) nitrogen isotope value of the baseline (corrected for TL = 1) δ15N baseline (TL=1), trophic level of the baseline species and the average (± S.E.) of the nitrogen isotope fractionation factors ( 3.43 ± 0.10 δ 15 N was 4.9 (± 1.90), with the lowest average baseline reported for Svalbard (4.2 ± 1.14) and the highest average baseline reported for the Canadian Archipelago (5.6 ± 1.64). Figure 3 (lower graphs) shows the calculated speciesspecific trophic level distributions for pelagic and benthic food webs in the four geographical areas. All TL distributions were normally distributed (P > 0.05, Shapiro-Wilk test) and Levene's test showed no significant differences in variation between/within regions (P > 0.05). For 21 of the 29 pelagic species in the pelagic food web, statistically significant differences in TLs were observed when comparing TLs between regions (P < 0.05, One-way ANOVA). For benthic food webs, statistically significant differences in TLs were observed for 59 of the 78 species.  Figure 4 shows the mean coefficients of variation (CVs) of the TL estimates at the intra-sample, intra-study, intraregion and inter-region level, i.e. averaged over all species. In general, CVs increased with aggregation level, with the smallest CV pertaining to intra-sample variation and the largest CV to inter-region variation. One exception was the pelagic food web of the Svalbard Archipelago, where the highest mean CV was calculated for intra-region variation. Furthermore, the intra-sample CV of benthic food webs in Alaska and the Canadian Archipelago was higher than the intra-study CV. Variances and CVs calculated for all individual species in all four regions are shown in Figs. S4, S5 and S6, S7 in the Supporting Information, respectively (ESM_3.pdf).

Identifying drivers of variation in species' trophic level
The first two principle components (PC1 and PC2) explained 48.7% of the variation in TLs within the pelagic food-web data and 34.8% in the benthic food-web data (Fig. 5). Components PC3 and PC4 together explained 18.2% and 23% of the variation in the pelagic and benthic food web, respectively. Looking at the correlation between variables included in the principal component analysis revealed that both longitude and sampling month showed to be negatively correlated with trophic position, implying that pelagic organisms sampled at a higher longitude (or sampled later in the year) showed to have a lower trophic position. For latitude, no such correlation with trophic level was found for pelagic food webs. In contrast, in benthic systems, latitude showed to be correlated with trophic level, albeit it a weak correlation. In the PCA plot for the pelagic food web, no correlation between baseline parameters (baseline trophic level and baseline nitrogen isotopic value) and fractionation factor was observed, implying that, given Eq. 1, little variation was observed within used fractionation factors. A strong overlap of the ellipses of the Canadian Arctic and the Svalbard (the European Arctic) in both PCA biplots was observed, indicating a strong uniformity of these regions. In the benthic food web, only longitude was negatively correlated with PC1. However, none of the factor loadings exceeded 0.5, implying that none of the variables have a strong correlation with the principle components. PC2 was negatively correlated with latitude and longitude and the nitrogen isotopic signature in the benthic food web (loadings > 0.4). In the pelagic food web, PC2 was negatively correlated with the trophic level of the baseline species and its nitrogen isotope. However, again none of the factor loadings exceeded 0.5 (Table S4, Supporting Information).
In both food-web types, PC1 and PC2 showed a strong significant discrimination between Alaska and Svalbard. In pelagic food webs, biota sampled from the Alaskan Arctic had lower PC1 values than biota sampled from Svalbard, which was likely mainly attributed to differences in longitude. In addition, biota sampled in Alaskan pelagic food webs had lower nitrogen isotopes than biota from Svalbard. In benthic food webs, the discrimination between Svalbard and Alaska could be mainly explained by the variation in PC2 axis. Again, nitrogen isotopic ratios of the baseline species and its corresponding trophic level in general were lower in biota in Alaskan benthic food webs.
Assessing the relative importance of predictors of trophic level, the pelagic food web showed female sex (6.4%) and baseline δ 15 N (5.75%) to be the strongest predictors, after the reported stable nitrogen isotope itself. The relative importance of TL parameters, notably parameters associated with baseline species, generally was higher than the relative importance of extraneous variables (e.g. sampling year or region; Fig. S3 and Table S5, Supporting Information). One notable exception is sampling month (5.3%), suggesting strong seasonality of trophic levels of species. In benthic food webs, again baseline species δ 15 N and baseline Fig. 4 Mean intra-sample, intra-study, intra-region and inter-region coefficients of variation for calculated TLs over all species for both pelagic (left) and benthic (right) food webs, for all four regions trophic level were more important than most of the external variables, with 6.8% and 3.1%, respectively. Sampling year showed to be the strongest external predictor of trophic level in benthic food webs, which may be due to inter-annual fluctuations in environmental conditions influencing plankton composition in Arctic oceans (Batten et al. 2018). The 1 3 pelagic food web confirmed the variables month (6.5%) and sex (6.1%) as the strongest predictors of trophic level, after the stable nitrogen isotopic values itself, showing these variables to be a stronger predictor for trophic level than any of the trophic level parameters.

Trophic level estimates
On average, the highest average nitrogen isotopes and TLs were calculated for pelagic species sampled in Alaska. When comparing averages of common species, taking into account sampling preference of a certain species, statistically significant differences (P < 0.05, one-way ANOVA) were observed between regions for both food-web types. Median trophic levels, based on species averages reported in this study, are similar to the median nitrogen isotopic baseline-corrected TL (3.17 ± 0.88) and median TL estimates (3.32 ± 0.79) for Arctic areas reported by Carscallen et al. (2012). Furthermore, these values were similar to TL estimates for Arctic areas based on stomach content or fatty acid studies for Arctic marine mammals (Pauly et al. 1998;Trites 2001).
The calculation of trophic levels and its variability to compare TLs of species between regions depend on a number of assumptions. First, we assumed all input parameters (e.g. δ 15 N consumer , δ 15 N baseline and ∆ 15 N) not to be correlated, as a positive correlation between these parameters may result in an underestimation of the variation in estimated TLs. However, a correlation test revealed that none of these parameters were correlated (r < 0.1). Our second assumption pertains to the use of an accurate nitrogen fractionation factor (Δ15N). In this study, all fractionation factors applied in estimating TLs of Arctic species in benthic food webs were between 3.4 and 3.8. In contrast, Δ15Ns used in studies focusing on pelagic food webs were much lower, which is due to the low factor applied in ringed seals (P. hispida) in multiple studies in the Canadian Archipelago (2.4‰) (Brown et al. 2016;Houde et al. 2017). Although the majority of studies included in the present study use fractionation constants as reported by Post (2002) (3.4) or Hobson and Welch (1992) (3.8), a wide range of fractionation factor values were reported in the literature (Lecomte et al. 2011;L'Hérault et al. 2018). Multiple studies revealed that the nitrogen fractionation factor Δ15N may decrease with increasing dietary protein quality in the organisms diet, implying that nitrogen discrimination will be lower for species at higher trophic levels (i.e. carnivores) than for herbivorous or omnivorous species (Robbins et al. 2005(Robbins et al. , 2010Florin et al. 2011). Next to protein quality, protein quantity may also play an important role in variation in nitrogen fractionation, although this relationship showed to be very complex (Florin et al. 2011;Hughes et al. 2018). As enrichment of δ 15 N may be highly influenced by temperature, growth rate, dietary protein quality and sampling tissue (Lecomte et al. 2011;Clark et al. 2019), estimates derived from experiments may not reflect factors as observed in natural settings. However, the lack of understanding of sources of variation in Δ15N within and between species and its lack of experimental validation is very common in the field of wildlife studies (Morrissey et al. 2004;Lecomte et al. 2011). To account for intra-population variation in Δ15N, there is a strong need for unbiased, statistically valid taxa and tissue-specific estimates (Vanderklift and Ponsard 2003;Auerswald et al. 2010;Hussey et al. 2014).
A third assumption in this study is that species δ 15 N have been sampled and measured correctly. Tissue samples exposed to air tend to degrade, possibly enriching δ 15 N between 0.6 and 1.3‰ (Perkins et al. 2018). Furthermore, chemical lipid extraction techniques to obtain lipid-free δ 15 N may sometimes affect measured stable nitrogen isotopes, as amino acids may leach from the tissue, reducing the accuracy of the TL estimates (Sotiropoulos et al. 2004;Clark et al. 2019). However, δ 15 N is assumed to increase with 0.3-0.5‰ due to lipid removal, which is relatively small compared to the typical fractionation factor used in many studies in trophic ecology (Sotiropoulos et al. 2004). In the present study, variance induced by lipid removal techniques across studies pertained to intra-regional variance.
Our final assumption is that the temporal and spatial scale at which we aggregated that the stable isotope data are sufficiently small enough to exclude major variation induced by differences between ecosystems. Individual TL estimates were based on associated baseline species and δ 15 N baseline , including variation in basal resources, sampled within the same study and site. Nevertheless, inappropriately aggregating data at too large a scale may lead to introducing bias in the derived TL probability distributions, due to for instance different offsets of algae blooming at higher latitudes, potentially leading to different results and conclusions (Falk-Petersen et al. 2009;Nielsen et al. 2018).

Variability in trophic level estimates
In general, the variance calculated between the regions was higher than the variance between samples and studies within one region, for both pelagic and benthic food webs. This is especially the case for benthic species, for which the overall mean calculated variance between regions was a factor 2.5 higher than the variance determined for individuals, samples and studies.

Identifying drivers of variation in species' trophic level
In both food webs, in general, TL parameters (notably, parameters related to baseline species: TL baseline and δ15N baseline ) were stronger correlated with TL estimates than any of the other included variables, with one notable exception being sampling month. These findings are in line with the generally adopted idea of the choice of an appropriate baseline being the most important in trophic studies using stable nitrogen isotopes, regardless of the research question (Post 2002). TLs showed to be highly correlated with the stable nitrogen isotopes of the baseline species in both food webs. This again underlines the importance of choice of baseline species. Unlike in the benthic food web, the TL of species in the pelagic food web was strongly correlated to sampling month. This suggests that seasonality in trophic levels may be larger in the species from pelagic food webs than for benthic species. This may particularly be the case for the polar bear, covering a substantial part of our dataset, a species well known for its seasonality in feeding habits (Iversen et al. 2013). As sea ice declines in summer, a fraction of polar bears following the retreat of ice beyond the continental shelve may get food deprived, as they switch to alternative food sources residing at lower trophic levels, resulting in lower measured δ 15 N values . Furthermore, in pelagic marine systems, seasonal changes in lipid content of pelagic baseline species occur due to negligible light in winter and changing POM composition and algae blooms in spring and summer. This may result in differences in isotopic signatures of zooplankton across seasons, especially when comparing between studies where researchers did not remove lipids prior to isotopic analysis (Hobson and Welch 1992;Søreide et al. 2006). Sex showed to be a relatively strong predictor of TL in pelagic food webs in the present study (6.1%), with females showing higher average TLs than males. This may be due to the relative large proportion of mammalian marine predators in our dataset. Multiple studies (Hobson et al. 1997;Lowry et al. 1980;Dehn et al. 2007;Thiemann et al. 2008) support our findings and report higher δ 15 N in female marine mammals and/or find significant differences in prey composition between sexes (with females being more carnivorous than males). Dehn et al. (2007) propose that these differences in prey composition in seals may be partially explained by spatial segregation and differential use of resources, which is also visible in polar bears during the breeding season (Laidre et al. 2013).
Although it is known that the life cycle of pelagic zooplankton species shorten with latitude (Lischka and Hagen 2007), no strong correlation was found between latitude and TLs in the present study. An additional potential predictor of trophic level includes body size. Although this predictor was not included in the present study, due to lack of data, a strong correlation between body size and trophic level was acknowledged in multiple recent studies, especially for fish and mammals (Lesage et al. 2001;Landry et al. 2018).

Recommendations and implications
To our knowledge, the present study was the first to quantify variation of trophic levels of species across the marine Arctic. Unfortunately, many studies in trophic ecology rely on small numbers of stable isotope data, scattered throughout the time and space, due to financial and logistical constraints (Ward et al. 2010), making it challenging to quantify variation in TLs at one time and region. Nevertheless, multiple studies did account for variation by means of experimental design (e.g. Bayesian inference), similar to methods proposed in the present study (Starrfelt et al. 2013;Kim et al. 2016). Generating TL distributions based on occurring variation in δ 15 N among species, as well as methodological variation in TL parameters, rather than using a single discrete TL value for one species, allows for a better understanding of feeding behaviour of species in the Arctic. Understanding trophodynamics across all trophic levels and scales in Arctic areas is highly important in modelling impacts of anthropogenic drivers (Loseto et al. 2009;Cherry et al. 2011). TL distributions may be useful in providing naturally realistic predictions in bioaccumulation or climate change studies, as these are based on the probability of a species situating at a certain trophic level. Therefore, a logical next step is to apply these distributions in bioaccumulation models, to evaluate their implications in actual modelling practice.