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 high-latitude 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 (δ15N) of a consumer is systematically enriched by 2–4% relative to the consumer’s diet (Post 2002), making δ15N 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 δ15N 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 δ15N 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 δ15N, 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 intra-sample, 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 intra-sample, 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.

Materials and methods

Data collection and compilation

Data on stable nitrogen isotopes for Arctic species, and corresponding TL parameters (see Eq. 1; δ15NBaseline, TLbaseline, Δ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.DigitizeIt.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)).

Fig. 1
figure 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

Full size image

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 (δ15N), according to the following general model (Post 2002; Jardine et al. 2006):

$${\rm{TL}}_{\rm{consumer}}={\rm{TL}}_{\rm{baseline}}+\frac{{{\updelta }^{15}{\rm{N}}}_{\rm{consumer}}- {{\updelta }^{15}{\rm{N}}}_{\rm{baseline}}+ \Delta D}{\Delta 15{\rm{N}}},$$
(1)

where TLconsumer represents the tropic level of the consumer, TLbaseline represents the trophic level of the baseline species for a particular region and study (primary consumer, typically C. hyperboreus), δ15Nconsumer the measured stable nitrogen isotope of the consumer, δ15Nbaseline the measured stable nitrogen isotope of the baseline species, and ∆15N represents the nitrogen enrichment (or fractionation) factor of δ15N across trophic levels. In some cases, the study included an additional diet-tissue isotopic fractionation factor (\(\Delta 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, δ15Nbaseline was calculated according to Post (2002):

$${{\updelta }^{15}{\rm{N}}}_{\rm{baseline}}=\alpha \cdot {{\updelta }^{15}{\rm{N}}}_{{\rm{baseline}}, 1}+ {\left(1-\alpha \right)\cdot {\updelta }^{15}{\rm{N}}}_{{\rm{baseline}}, 2},$$
(2)

where δ15Nbaseline, 1 and δ15Nbaseline, 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 δ15Nbaseline 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 δ15N, 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, 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. δ15Nconsumer, δ15Nbaseline and ∆15N) 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.

Fig. 2
figure 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

Full size image

Uncertainty and variability in trophic level estimates

Variation at lower aggregation levels propagates into higher aggregation levels (Fig. 2) according to the following equation (e.g. Wallace and Williams 2005):

$${{\sigma }_{{\rm{measured}}, i}}^{2}={\rm{VAR}}_{i}+ \frac{1}{{N}_{i-1}}\cdot \overline{{{{\sigma }_{{\rm{measured}}, i-1}}^{2}}} ,$$
(3)

where σmeasured,i2 is the variance measured at aggregation level i, VARi is the unbiased variation at level i, \(\overline{{{{\sigma }_{{\rm{measured}}, i-1}}^{2}}}\) is the average measured variance at aggregation level i − 1 and Ni−1 is the sample size at level i − 1. This equation was rearranged to obtain the unbiased variation at each aggregation level:

$${\rm{VAR}}_{i}={{\sigma }_{{\rm{measured}}, i}}^{2}- \frac{1}{{N}_{i-1}}\cdot \overline{{{{\sigma }_{{\rm{measured}}, i-1}}^{2}}} .$$
(4)

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 (inter-region), 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. δ15Nbaseline, TLbaseline and Δ15N) and extraneous factors (δ13C, 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.

Results

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 stable nitrogen isotopes (δ15N) 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 δ15N 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).

Fig. 3
figure 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

Full size image

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).

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) δ15Nbaseline (TL=1), trophic level of the baseline species and the average (± S.E.) of the nitrogen isotope fractionation factors (Δ15N)
Full size table

Baseline species

The majority of TLbaselines 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 TLbaseline for Alaskan pelagic food webs of 1.88 ± 0.5. In benthic food webs, TLbaselines 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 δ15N, 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 δ15N reported in Svalbard (4.2 ± 1.14‰) and the highest average baseline δ15N reported in the Canadian Archipelago (5.6 ± 1.64‰). Significant differences in baseline δ15N, 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 δ15N among all regions was 5.61 ± 1.82‰, with the lowest average baseline δ15N 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 δ15N 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).

Trophic level estimates

Figure 3 (lower graphs) shows the calculated species-specific 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. Table 2 summarizes the TL estimates by reporting the average TL per species per region and per food web. Highest average TLs, based on the average of all species weighted averages, were calculated for Alaska and the Canadian Beaufort Sea in pelagic food webs (3.00 and 2.98, Table 2), and for Alaska and Svalbard in benthic food webs (3.00 and 2.95). The lowest averaged TLs are calculated for the Canadian Archipelago (2.73) and Svalbard (2.50) for pelagic and benthic food webs, respectively. Statistically significant differences were observed in species average TLs between the areas for both the pelagic and benthic food web (P = 0.027, mANOVA, See Tables S2 and S3 for all species averages for δ15N and TLs, respectively).

Table 2 Descriptive statistics pertaining to the collected stable nitrogen isotopes and derived trophic levels per region, aggregated per species
Full size table

Variability in trophic level estimates

Figure 4 shows the mean coefficients of variation (CVs) of the TL estimates at the intra-sample, intra-study, intra-region 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).

Fig. 4
figure 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

Full size image

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).

Fig. 5
figure 5

PCA biplots for both Arctic food webs: the pelagic web on the top (a) and the benthic food web below (b)

Full size image

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 δ15N (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 δ15N and baseline 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 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.

Discussion

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. δ15Nconsumer, δ15Nbaseline and ∆15N) 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, 2010; Florin 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 δ15N 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 δ15N have been sampled and measured correctly. Tissue samples exposed to air tend to degrade, possibly enriching δ15N between 0.6 and 1.3‰ (Perkins et al. 2018). Furthermore, chemical lipid extraction techniques to obtain lipid-free δ15N 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, δ15N 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 δ15Nbaseline, 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: TLbaseline and δ15Nbaseline) 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 δ15N values (Whiteman et al. 2018). 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 δ15N 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 δ15N 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.