Implications of Trophic Variability for Modeling Biomagnification of POPs in Marine Food Webs in the Svalbard Archipelago.

The occurrence of persistent organic pollutants (POPs) in the Arctic has been of constant concern, as these chemicals cause reproductive effects and mortality in organisms. The Arctic acts as a chemical sink, which makes this system an interesting case for bioaccumulation studies. However, as conducting empirical studies for all Arctic species and POPs individually is unfeasible, in silico methods have been developed. Existing bioaccumulation models are predominately validated for temperate food chains, and do not account for a large variation in trophic levels. This study applies Monte Carlo simulations to account for variability in trophic ecology on Svalbard when predicting bioaccumulation of POPs using the optimal modeling for ecotoxicological applications (OMEGA) bioaccumulation model. Trophic magnification factors (TMFs) were calculated accordingly. Comparing our model results with monitored POP residues in biota revealed that, on average, all predictions fell within a factor 6 of the monitored POP residues in biota. Trophic variability did not affect model performance tremendously, with up to a 25% variability in performance metrics. To our knowledge, we were the first to include trophic variability in predicting biomagnification in Arctic ecosystems using a mechanistic biomagnification model. However, considerable amounts of data are required to quantify the implications of trophic variability on biomagnification of POPs in Arctic food webs.

Introduction

Persistent organic pollutants (POPs) in the Arctic have been of concern since the mid-1980s. These compounds are known to cause adverse effects (e.g., altered reproduction and immunological effects in organisms at higher trophic levels).[1−4] Despite the banning of many of these compounds in the 1970s and 1980s, high levels of POPs are still being reported in Arctic biota. This is due to the combination of chemical properties of POPs in cold conditions (e.g., high thermal stability and slow degradation) and marine currents and airflows transporting POPs northwards.[5−10] Due to the accumulation of chemicals via limited evaporation at lower temperatures and recirculating air and water currents,[11] the Arctic is a global chemical sink. This makes the Arctic a relevant case for bioaccumulation studies.[12−14] The magnification potential of chemicals is evaluated by calculating trophic magnification factors (TMFs), reflecting the degree of magnification of the chemical throughout the food chain. These TMFs are defined as the slopes of the log-transformed contaminant concentrations (in lipid weight) on trophic position and rely heavily on monitoring data or experimental studies.[15] Although the TMF concept is accepted as one of the most realistic quantitative measures of biomagnification,[16] there is a lack of studies assessing changes of TMFs across large spatial scales and ecosystems.[17] Conducting extensive empirical studies for all Arctic species and POPs individually is unfeasible, due to ethical, financial, and practical constraints.[18] Furthermore, TMF measurements are retrospective and are only feasible after a substance has already been released into the environment.[16] As TMFs are based on univariate linear regression models, these do not account for species-specific biotransformation of POPs.[19]In silico mechanistic methods have been developed to predict biomagnification of hazardous compounds in food chains, including models describing bioaccumulation kinetics.[20−22] These models generally require water concentration and the trophic position of a species in a particular food chain as input, as well as specific species traits (e.g., body mass and lipid content) and chemical-specific characteristics (e.g., the octanol–water partition coefficient).

Bioaccumulation models are predominately validated for temperate food chains.[21,23,24] The Arctic is highly influenced by seasonality, which affects the degree of bioaccumulation of compounds in Arctic food webs throughout the year.[25−28] Seasonal differences in the Arctic are mainly driven by fluctuations in light, temperature, and sea ice affecting phytoplankton blooms, the species’ lipid content, and body size. This in turn may affect the community structure and, therefore, food and prey availability and may consequently lead to diet shifts.[25,27,28] Variation in the baseline species and temporal and spatial differences in trophic status of species can have substantial implications for the performance of the biomagnification model.[29,30] Although the influence of trophic variability of species in Arctic ecosystems on the transfer of POPs in food webs has been previously studied,[17] its influence on the performance of mechanistic biomagnification and bioaccumulation models is yet to be examined. Intrinsic trophic variability and diet shifts throughout the year may lead to differences in steady-state concentrations of POPs in Arctic species. The aim of the current study was to better predict the bioaccumulation of POPs in Arctic food chains, which are subjected to high trophic variability using the optimal modeling for ecotoxicological applications (OMEGA) bioaccumulation model. Publicly available data for one Arctic region, i.e., Svalbard, were used as input for the mechanistic biomagnification model to predict the transfer of POPs across food chains, accounting for variability in the trophic position of species by means of Monte Carlo simulations. Model performance data were evaluated against measured POP concentrations in the Svalbard Archipelago.

Materials and Methods

The OMEGA Bioaccumulation Model

Model Approach

Biomagnification of POPs in the present study was calculated using the OMEGA bioaccumulation model,[31] combining mass balance and allometric approaches to estimate bioaccumulation in each tropic level individually and biomagnification throughout the food chain.[24] Since its development, OMEGA has effectively been applied in estimating biomagnification of numerous organic compounds and metals.[23−25] We assumed a steady-state concentration of chemical residues in aquatic organisms (C, in μg/kg lipid weight (l.w.)), calculated as the sum of influx via water (absorption) and uptake from food (assimilation) divided by the total elimination rate, encompassing elimination via water, biotransformation, excretion, and dilution through growth, disregarding possible uptake and elimination through air. The steady-state concentration in an organism is calculated as[31]where C is defined as the final steady-state concentration of substance x at trophic level i (in μg/kg lipid weight), k0, is the chemical absorption rate (in L/kg·d–1), C0, is the initial chemical concentration in the water (in μg/L), k1, is the chemical concentration assimilation rate in [kg/kg·d–1], C is the chemical concentration in food (in μg/kg), and ∑k is the sum of the elimination rate constants (per day), pertaining to water (including urine), feces, growth dilution, and biotransformation. Concentration estimates in the model were based on wet weights and were transformed into lipid weight-based estimates using reported data on the whole-body lipid content. As the Arctic is a highly dynamic system, with algae blooms in spring and changes in the lipid content of zooplankton and species at higher trophic levels throughout the year, the assumption of a steady-state concentration may not be justified in a highly dynamic system such as the Arctic, especially for lower trophic levels.[25] Furthermore, the OMEGA model only includes the main aquatic exposure routes and does not take into account any possible exposure of POPs through air (Text S1 in the Supporting Information (SI)).

Including Variability in Trophic Levels in OMEGA

The trophic level concept is predominantly used to evaluate food webs and food chains since these levels are relatively easy to define and allow for estimation of energy and contaminant flow through ecosystems.[32] However, the use of discrete trophic levels fails to acknowledge complex interactions between organisms within a food web. In contrast, stable isotopic techniques provide a continuous measure of the trophic position, integrating contaminant flows from all different trophic pathways, leading to the species of interest.[32] The use of relative trophic levels based on stable nitrogen isotopes (δ15N) allows for the ecosystem and species-specific discrimination by means of unique nitrogen enrichment factors, and baseline correction.[17] The isotopic signature of an organism alone does not provide sufficient information on absolute trophic position when comparing among multiple systems, as both the isotopic signature of baseline species and the nitrogen enrichment factors are location and food web specific.[30] Therefore, we calculated the trophic position of target species based on the nitrogen enrichment factors and isotopic signatures of baseline species as reported in the literature,[33] according to the following general model[32,34]where TLconsumer represents the tropic level of the consumer, δ15Nconsumer represents its corresponding stable nitrogen isotope, δ15Nbaseline represents the measured stable nitrogen isotope of a species set at a baseline trophic level for a particular region and study (equal to TLbaseline), and Δ15N represents the enrichment factor of δ15N across the food chain. As the variance in trophic levels induced by variation in δ15Nconsumer and Δ15N between regions is higher than the variance between individuals, samples, and studies within one Arctic region, the present study only focused on bioaccumulation of POPs in the Svalbard Archipelago.[33] Species from both pelagic and benthic food webs were considered since both food webs are important for top predators.

For each individual data record, the mean trophic level and its standard deviation were calculated based on eq , including reported corresponding variation in δ15Nconsumer, δ15Nbaseline, and Δ15N by means of Monte Carlo simulation (1000 iterations). Trophic levels and standard deviations per data record were then used in constructing trophic level distributions per species by means of geometric averaging and taking the pooled standard deviation based on sample size (Table S5, Supporting Information). At each iteration, species were distributed over six discrete trophic level bins (ranging from 0 to 5), based on a trophic position randomly sampled from the species-specific probability distribution of the trophic level, which acted as direct input for 1000 biomagnification scenarios in OMEGA in a Monte Carlo simulation (see Text S2 in the Supporting Information).

Parametrization of the OMEGA Model

Equations pertaining to influxes and effluxes of the OMEGA model (eq ) and corresponding parameter values, including the lipid content and the compound’s octanol–water partition coefficient (Kow), are provided in the Supporting Information (Text S1, Tables S1 and S2). Log-transformed Kow’s for the compounds were obtained using the experimental data from EPI-suite 4.1, disregarding any variation or uncertainty[35] (see Supporting Information Table S2). Estimates for the bodyweight per trophic level was obtained through Hendriks et al.[31] Whole-body lipid fractions were taken from the ICES database (http://ecosystemdata.ices.dk/). When no lipid percentage was reported, an average whole organism lipid fraction was taken based on the species group (see Table S3 in the Supporting Information). For mammals, an average whole-body lipid content of 30% was assumed. Compound-specific biotransformation rates for higher trophic levels (TL > 2) were based on median values reported by Arnot et al.,[36] which are standardized for a 10 g fish at a water temperature of 15 °C (see Table S4 for all biotransformation rates).

Since biotransformation rates tend to decrease with temperature,[37,38] we accounted for their temperature dependence by assuming an Arrhenius-type of equation (eq S6 in the SI).

Exposure and Validation Data

Data Collection

Data on POP residues in marine biota around Svalbard (72–85° N, 0–40° E) were compiled to calculate TMFs and to validate the OMEGA model. POP concentrations in both benthic and pelagic biota, transformed to μg/kg lipid weight (l.w.), were downloaded from the database of the International Council for the Exploration of the Sea (ICES; http://ecosystemdata.ices.dk/) and complemented with the data from a literature search using the Web of Knowledge and Google Scholar. The search strings used in queries included names of the most recognized legacy POPs (including “p,p′-DDE”, “p,p′-DDT”, “Hexachlorobenzene”, “Dieldrin”, or “PCB-153”) or general terms and compound groups (“Persistent organic pollutants/compounds”, “POPs”, “PCBs”, “DDTs”, “HCHs”, “PAHs”, “organochlorines”, “PFASs”, or “PBDEs”). These were combined with general terms related to Arctic food webs (“Svalbard food web”, “Svalbard”, “Arctic”, “Arctic food web”, “Arctic ecosystem”, “Arctic pelagic food web”, or “Arctic benthic food web”), names of well-known Arctic pelagic and benthic species (both scientific (e.g., “Ursus maritimus”, “Phoca hispida”, or “Calanus hyperboreus”) and their common names (e.g., “polar bear”, “ringed seal”, or “copepod”)).

Data Treatment

Residue concentrations of POPs in biota were transformed to μg/kg lipid weight (l.w.) using the reported whole-body lipid content (Table S3 in the Supporting Information). Duplicate values were removed using the Duplicate function in R statistics 3.4.1, based on the date of sampling, species, sex, location, sampling tissue, compound, and value. Although the obtained data were sampled throughout the year, due to the lack of data sampled in winter, only samples obtained for the months April until October and after the year 2000 were considered in further analysis. Three inclusion criteria were applied to the dataset to select compounds to be included in the dataset. As the OMEGA model is largely dependent on water concentrations, including a sufficient amount of water data (n ≥ 5), it is crucial for the model to reflect a natural system. However, data on POP concentrations in Arctic coastal and marine waters is grossly lacking in the literature. Therefore, as a first criterion (#1), when no sufficient water concentration data were available (n < 5), water concentrations were estimated based on POP concentrations in zooplankton (n ≥ 5), reflecting the second trophic level of the Arctic food chain. These estimates were based on monitoring the data of POP residues in biota residing at the second trophic level and corresponding species-specific parameters. Given that, in some cases, the residue data were below the detection limit, a second criterion (#2) was included, relating to the fraction censored data. Chemicals for which over 40% of the water data records fell below the detection limit were discarded. As a third criterion (#3), we required that at least five species were monitored per compound to ensure that multiple species and trophic levels were represented in the model. As a fourth and final criterion (#4), we required the remaining chemicals to cover at least four trophic levels. Reported POP concentrations in biota were geometrically averaged per species and location, and weighted according to the reported sample size (based on the amount of measurements) per data record (Figure ).

Figure 1

Locations of data records included in the Svalbard Archipelago, encompassing POP residue concentration data in Arctic biota, after truncation according to set criteria, obtained from a literature search. The size of the points indicates the sample size.[39,40]

Trophic Magnification Factor

The trophic magnification factor (TMF) is a typically used metric to assess the accumulation of contaminants in a food chain. Parallel to modeling biomagnification of POPs using the OMEGA model, TMFs were simulated (1000 iterations) based on the monitored average POP residue data per species and corresponding randomly sampled trophic levels per species taken from the constructed species-specific trophic level distributions, as described earlier. Unlike in the OMEGA model, TMFs were only calculated for ectothermic biota, as ectothermic animals are more likely to biomagnify chemicals than endothermic species.[41] Concentrations of POPs can generally be described by a linear relationship[15,17,41,42]where log(C) is the log10-transformed lipid-normalized chemical concentration in ectothermic biota (phyto- and zooplankton, bivalves, and fish), aggregated per trophic level bin, a reflects the intercept of the regression line and accounts for input of contaminants at the base of the food web,[17] and b is the regression slope coefficient. The TMF is then calculated from the regression slope b

Model Validation and Performance

To asses model performance, the coefficient of efficiency (E, dimensionless) was calculated for each persistent organic compound, through[43]where P reflects the log-transformed mean residue concentration per trophic level bin as predicted by the OMEGA model for simulation i, O represents each corresponding log-transformed observed residue concentration, and O̅ represents the log-transformed mean of all observed values within each iteration. The efficiency coefficient E ranges from minus infinity to 1, with a value of 1 indicating a perfect model estimation. A positive E indicates that the model estimates are more accurate than the mean of the observed values, indicating a good fit.[23,43] To summarize both random error and systematic bias, the root-mean-square-error for each model iteration was computed by geometrically averaging the squared differences between log-transformed observed and predicted valuesTo summarize the relative error of the model, we computed the median symmetric accuracy (ξ) per compound based on results from each iteration and the trophic level (eq ), reflecting the typical percentage error of the predicted POP concentrations compared to the measured concentrations, without penalizing over- and underestimation differently,[44,45] according toin which ξ is based on the exponentiated median of the log-transformed predicted concentrations in biota (ŷ, μg/kg) and the log-transformed measured concentrations in biota (y, μg/kg). ξ can be interpreted as a typical percentage error: a ξ of 200%, for example, indicates that predicted concentrations from OMEGA will typically be within a factor of 3 of the measured concentrations in Arctic biota.[46] Finally, to assess the prediction bias of OMEGA, accounting for likely asymmetry in the prediction errors, we computed the symmetric signed percentage bias (SSPB)[45]The SSPB may be interpreted similarly as the mean percentage error, with values below 0 indicating a systematic underprediction. Model performance metrics, as previously described, were compared among compound groups by performing a one-way analysis of variance (ANOVA). Compound groups were compared individually by performing a Fisher’s least significant difference test. All work pertaining to the biomagnification modeling and model validation was performed in R statistics 3.4.1.[47]

Results

Exposure and Validation Data

The initial literature search on biota residue data resulted in 50 articles/reports, covering over 170 individual compounds from 9 chemical classes, encompassing 18 000 data records (Table ). Application of the first three inclusion criteria, related to the number of concentration data, resulted in a secondary dataset of approximately 16 500 data records (49 articles, encompassing 90 compounds for Svalbard (Table )). By including data on the trophic position, biomagnification could be modeled for a total of 22 compounds in Svalbard, based on 44 articles and 10 000 data records (see the Supporting Information for the full reference list).

Table 1

Number of Data Records before and after Truncation of the Collected Data, According to the Four Criteria as Described in the Text (Relating to the Number of Species and TLs, the Number of Water/Zooplankton Data and the Amount of Censored Data) and the Presence of Sufficient Trophic Position Data

	initial dataset	#1 water criterion	+ #2 detection limit criterion	+ #3 species concentration data criterion	+ #4 species TL data criterion
number of data records	18 318	16 548	16 548	16 531	10 044
number of compounds	173	91	91	90	22
number of compound groups	9	8	8	8	6
number of species	67	66	66	66	63
number of articles	50	49	49	49	44

OMEGA Model Performance

Biomagnification could be modeled for a total of 22 compounds in the Svalbard Archipelago (Table S6 in the Supporting Information (SI)). Model results for each substance are shown in Figure S1 in the SI. In total, 60.8% of all predictions fell within a factor 5 of the observed values, with concentrations in biota being equally underestimated (in 18.8% of the cases) as well as overestimated (19.9%) (Figure ). As expected, a substantial concentration step between the water phase (TL = 0) and phyto- and zooplankton (TL = 1) was observed for almost all individual compounds (Figure S1, SI). Mean efficiencies (see eq ) calculated for the OMEGA model, accounting for variation in the trophic position of Arctic species, ranged from −0.83 (CB-101) to 0.66 (CB-99) (Figure left graph and Table S6). The majority of the studied compounds (18/22) had a positive mean efficiency value, implying that the model predictions were more accurate than the average of the observed concentrations. Overall, model runs for cyclodienes had the highest calculated averaged efficiencies (0.38 < E < 0.64), followed by runs for DDTs (0.42 < E < 0.51) and PBDEs (0.23 > E < 0.57). Lowest mean efficiencies were calculated for HCHs (−0.38 < E < 0.57). No significant differences in model performance were observed between compound groups (p = 0.666, one-way ANOVA). Calculated RMSEs were highest for CB-180 (2.72, Table S6) and lowest for CB-28 (0.86). The highest absolute RMSEs were found for cyclodienes (1.82 ± 0.64), followed by PBDEs (1.74 ± 0.02). However, again, no significant differences in RMSEs were observed among compound groups (p = 0.962, one-way ANOVA). Median symmetric accuracies were calculated for all individual compounds, compound groups, and trophic level bins, reflecting the typical percentage error compared to the measured data.[46] Overall, the typical prediction error of all model iterations was 288.49%, implying that typically predictions were within a factor 3.9 from the observed concentrations. The lowest typical percentage error, and thus highest predictive power, was calculated for oxychlordane (68.5%), while the highest typical percentage error was calculated for CB-180 (1081%, Table S6 in the SI). Average model accuracies were highest for PBDEs (ξ = 145.37%), followed by HCHs (168.8%), while the highest average percentage errors were computed for PCBs (398.84%) (Table ). No statistically significant differences were observed in prediction errors among the modeled compound groups (p = 0.378, one-way ANOVA). Evaluating model performance per trophic level bin revealed that the typical percentage errors were lowest for the first four trophic level bins (ranging from 110 to 220%), and highest for the fifth trophic level. Despite the high standard deviation corresponding to the high prediction error of the fifth bin, a statistically significant difference in model performance, based on median symmetric accuracies, was observed between the fifth trophic level and the other trophic level bins (p = 0.001, one-way ANOVA with least significant difference). The average SSPB value across all simulations was 124.4, indicating an overestimation of the predictions relative to the monitored concentrations. The highest SSPB value, indicating the highest overestimation bias, was computed for CB-180 (888.15%), while the lowest SSPB was found for HCB (−93.6%, indicating an underestimation bias of a factor of almost 2). The SSPB closest to 0, indicating best model performance, was found for model results for BDE-47 (−1.55%) and CB-28 (12%). PCBs (252.24%) showed to have the highest SSPBs, while SSPBs closest to 0 were found for cyclodienes (36.26%).

Figure 2

Observed POP concentrations in Svalbard biota versus predicted model estimates (based on 1000 Monte Carlo iterations per compound and five trophic bins per iteration). The dashed lines represent a factor of 5 under- or overestimation by the OMEGA model. The colors correspond to the compound groups.

Figure 3

Model efficiencies including all species (left graph), empirical TMFs, including heterothermic species (right graph), and corresponding 95% confidence intervals due to variability in the trophic position of species, calculated for multiple compounds in Svalbard. The colors correspond to the compound groups. Additionally, the number of species included in the modeling practices per chemical is indicated for each substance between parentheses.

Table 2

Average Efficiencies, TMFs, Root-Mean-Squared-Errors (RMSEs), Median Symmetric Accuracies (ξ), and Coefficients of Variation for Efficiencies and TMFs per Compound Group in the Svalbard Archipelago. Different Letters Indicate Significant Differences Between Groups (p < 0.05, ANOVA)

	efficiency	Log10 TMF
	(±S.E.)	(±S.E.)	RMSE of OMEGA fit	ξ	SSPB	CVefficiency	CVRMSE	CVTMF
PBDEs (N = 2)	0.399 (±0.243)A	0.36 (±0.03)A	1.75 (±0.02)A	145.37 (±69.19)A	60.65% (±87.97)AB	25.7% (±33.3)B	24.4% (±0.8)A	16.02% (±4.7)A
PCBs (N = 10)	0.301 (±0.405)A	0.41 (±0.11)A	1.72 (±0.54)A	398.84 (±271.14)A	252.84% (±256.71)A	30.3% (±45.2)B	14.1% (±8.3)A	15.5% (±5.6)A
HCB (N = 1)	0.183A	0.41A	1.55A	251.23A	–93.64%B	130%A	18.5%A	9.3%AB
DDTs (N = 2)	0.469 (±0.065)A	0.53 (±0.02)A	1.45 (±0.37)A	228.34 (±64.01)A	68.3% (±45.34)AB	12.9% (±15.6)B	16.7% (±5.6)A	10.9% (±0.38)AB
cyclodienes (N = 4)	0.535 (±0.12)A	0.65 (±0.15)A	1.83 (±0.65)A	190.55 (±110.07)A	36.26% (±54.59)AB	7.9% (±7.7)B	25.9% (±12.1)A	15.9% (±4.1)A
HCHs (N = 3)	0.092 (±0.47)A	0.33 (±0.08)A	1.57 (±0.40)A	168.35 (±3.71)A	–33.57% (±53.95)B	57.9% (±70.1)AB	13.2% (±18.6)A	8.6% (±3.9)AB
total	0.342 (±0.33)	0.45 (±0.12)	1.69 (±0.47)	284.29% (±213.76)	124.41 (±212.75)	32.5% (±46)	18.3% (±9.6)	14.02% (±5.2)

Including Trophic Variability

Variation in estimated steady-state concentrations in biota was assessed by calculating coefficients of variation (CVs) for efficiencies, RMSEs and TMFs (Tables and S6 in the SI). The highest variability in efficiency induced by variability in trophic positions of species was observed for γ-HCH (138%), followed by CB-194 (134%). The significantly highest CV was calculated for HCB (130%), followed by HCHs (57.9%), while the lowest CVs were calculated for cyclodienes (7.9%) and DDTs (12.9%). However, no significant differences were observed among the compound groups (p = 0.221, one-way ANOVA). Variation in RMSEs induced by variability in the trophic position of the included species showed to be highest for oxychlordane (41.91%), and the lowest CVRMSE was calculated for CB-105 (4.7%). Variability in measured POP concentrations (TMFs) was relatively high (14.02 ± 5.2%, Table ). Both the variability in efficiencies and the variability in TMFs showed not to be correlated to the number of species included in analysis (p = 0.47 and 0.194, respectively, linear regression model).

Overall, average TMFs calculated in the present study ranged from 1.5 (CB-99) to 6.84 (Oxychlordane). Cyclodienes had, on average, the highest magnification potential, followed by DDTs and HCBs, while lowest TMFs were found for HCHs (Table ). TMFs did not differ significantly across compound groups (p = 0.588, one-way ANOVA, Table ). All assessed compounds in the present study are highly bioaccumulative (TMFs > 1, Log10 TMF > 0[48]) (Figure S2 in the Supporting Information).

Discussion

Our study shows how variation in the trophic position of species influences modeling performance of a biomagnification model when applied to modeling POPs in Arctic ecosystems. Considerable variation in trophic position was observed for the species included in the present study, resulting in varying model performance across compound groups (given the large differences in CVefficiency (Table )). However, this trophic variation did not have major implications for bioaccumulation, given the relatively low CVs found for TMF. Although biomagnification in Arctic marine food webs has been thoroughly examined on a relatively small temporary and spatial scale in separate case studies, implications of intra-specific variation in trophic levels on biomagnification on larger scales are often disregarded.[49] As in the present study, species were binned according to their trophic level, presumably feeding upon the previous bin, we stress that the food webs used in the present study do not necessarily reflect natural Arctic food webs. However, nitrogen stable isotopes of species on the basis of the trophic level estimation were compared to baseline species in the same system, thus representing a representative trophic level. Furthermore, stable nitrogen isotopic signatures of primary producers at the base of the food web are considered to be relatively spatially and temporally similar, hence the use of a fixed baseline trophic position (2).[50]

For the majority of modeled compounds (18/22), efficiencies pertaining to modeling biomagnification along Arctic marine food chains were positive, indicating that the mechanistic OMEGA model fits the monitoring data more accurately than taking the average of the observed concentrations. However, as with the R2, the coefficient of variation E is highly sensitive to outliers.[43] Concentration steps observed in the present study have been reported in data of several previous studies focusing on both Arctic and temperate food webs,[51−53] and may be explained by POP concentrations in biota at the base of the food web being in (near-)equilibrium with dissolved water concentrations.[51] In the present study, no major fluctuations in modeled concentrations at the base of the food web occurred, except for compounds lacking sufficient concentration data on species residing at TL = 1 (e.g., certain PBDEs, Figure S1 in the SI). The OMEGA model predicted less well for species residing at higher trophic positions (including fish and marine mammals) (TL = 4 and 5, Table S7 in SI). In the present study, we included biotransformation rates, based on values normalized for 10 g fish and a water temperature of 15 °C.[36] As biotransformation increases with increasing temperatures,[23,54] biotransformation may be lower in the Arctic. To account for changes in biotransformation associated with lower atmospheric temperatures, we transformed our biotransformation rates using an Arrhenius function.[13,55] However, biotransformation rates in all species were based on those for fish, so we did not take into account differences in species thermal groups. Although temperature especially affects the biotransformation rate of poikilotherms (in this case, fish), it also has some influence on the metabolic rate in homeotherms (mammals).[56] In the present study, biotransformation rates may be different for marine mammals residing at the higher trophic positions in the food chain.

Overall, median prediction errors for all compounds by OMEGA were within a factor of 6 from the observed concentrations (ξ < 500%). An average SSPB of 124% was observed, indicating an overestimation when concerning all model predictions. Large standard deviations in ξs and SSPBs were observed within all model predictions, as well as within compound groups, and no statistically significant differences were observed in ξs and SSPBs among the compound groups. Overall, HCHs showed to have the lowest efficiencies, indicating a poor fit. However, this may likely be due to the low calculated efficiency for HCH-γ as average RMSEs, ξs and SSPBs showed to be relatively low for these compounds. All HCHs included in the present study had relatively low octanol–water partition coefficients (Log Kow = 4.14, compared to Kow ranging from 5.48 to 8.68 for other compounds), limiting partitioning of a chemical from water into lipids, and thus bioaccumulation potential.[57] Although HCHs are efficiently metabolized in vertebrates at higher trophic levels, i.e., fish and marine mammals, gut microbial degradation may also play a role in bioaccumulation of HCHs in the Arctic food chain.[58]

Interpretation of field data may be sensitive to bias induced by uncertainty in trophic interactions, limited statistical power, and ecosystem-specific factors, such as water salinity and temperature, affecting the water solubility and organic matter partitioning of POPs.[17,19,59−61] In the present study, TMF values computed for individual POPs had an average CV of 14.01%. These TMFs take into account the observed variation in stable nitrogen isotopes as well as in the nitrogen fractionation factors used and the choice of baseline species.[33] Instead of taking a TMF based on a single food web, variability in TMFs can be quantified by computing distributions based on TMF simulations.[15] Rather than applying numeric pass/fail criteria to single field-TMFs, it is then possible to quantify the risk of certain chemicals exceeding these criteria based on a more wider background.[48] However, this may not yet be applicable in a regulatory context as knowledge is currently lacking on whether the variability in simulated TMFs is the result of natural variation in relevant processes or measurement errors, poorly defined concepts, and statistical analyses.[15] Given the rule of the thumb that TMFs < 1 indicate poorly bioaccumulative substances,[48] an inadequate choice of species in biomagnification modeling and incorrect use of trophic levels can result in chemicals being characterized as poorly bioaccumulative. Examples in the present study are BDE-154 and CB-99, for which 22.4 and 17.0% of all simulated TMFs fall below 1, respectively (Table S6, SI). Problems may arise in calculating TMFs when food webs include both poikilotherms and homeotherms, which is often the case with the Arctic pelagic food webs.[60] This problem was avoided by including only the first three trophic levels in the TMF calculations, excluding birds and mammals. Variation in biomagnification may be induced by different individual-specific covariates, such as length or lipid content, or methodological differences pertaining to the sampling protocol (i.e., sampling tissue).[25,62,63] Computed median TMFs in this study were slightly lower than—but within a factor of 2 of—the corresponding obtained TMFs gathered by Walters et al.[41] (see Table S8 in the Supporting Information). While standard errors for TMFs in the study by Walters et al. were relatively high, our computed TMFs were not significantly different from these measured TMFs. The analysis in the present study, as well as the study done by Walters et al. reveals PCBs as the compound group with the highest bioaccumulative potential. In contrast to the present study, the study by Walters et al. focused on TMFs measured in temperate zones. As mentioned previously, biomagnification of nonmetabolizable POPs is thought to be an exothermic process, implying that biomagnification decreases with increasing temperature.[64] Temperature may not only influence POP concentrations in biota because of its influence on metabolic (biotransformation) rate but may also affect the excretion rate by modifying the lipid content of the animal, especially in poikilothermic animals.[56] If this is the case, measured POPs are expected to be higher in Arctic biota, compared to temperate zones. This was not the case in the present study, reporting lower TMFs based on monitoring data for the Arctic than TMFs reported in temperate zones (by Walters et al., 2016). However, next to the climatic zone, other environmental and methodological parameters may influence the calculation of TMF.[41] In the study by Walters et al. TMFs in the Arctic zones were based on both ectothermic and endothermic organisms, while in the present study, only ectothermic species were considered in the TMF calculation. The same study found that TMFs are lower in food webs containing only ectotherms (or when TMFs are based on only ectothermic species, as in the present study), independent of latitude.[41] Furthermore, they found that TMFs were more variable in marine ecosystems compared to freshwater ecosystems.

Implications and Recommendations

Overall, OMEGA predicted reasonably well for the compounds included in the present study, with efficiencies for 18 out of 22 compounds being positive. Variation in the trophic position of species, due to both natural processes (seasonality) and methodological differences (choice of baseline species and nitrogen fractionation factor), accounted for up to a 25% variability in the model performance metrics. Considerable amounts of data were required, pertaining to both the studied chemicals (e.g., residues in biota and octanol–water partition coefficients) and the species (e.g., body mass and trophic position or stable nitrogen isotopic data) to realistically quantify the implications of intra- and inter-specific trophic variability on biomagnification of POPs in Arctic food webs. In particular, data on POP concentrations in zooplankton were relatively sparse, disallowing biomagnification modeling for many compounds. As several criteria were set to increase statistical robustness (representativeness) of the model, the demand for more data increased accordingly.

Application of toxicokinetic mass balance biomagnification models in Arctic ecosystems is rare, due to the variability of metabolic rates and rates of decay at lower temperatures and the confirmed high seasonal variability in the biological activity of Arctic organisms and, therefore, biomagnification.[25,65] To our knowledge, we were the first to include trophic variability in predicting biomagnification of chemicals in Arctic ecosystems using a mechanistic biomagnification model, including temperature-dependent biotransformation rates. As variation in modeling performance results showed to be relatively low, we may conclude that the variation in the trophic position of Arctic species may not have a major influence on the degree of biomagnification of POPs, based on evaluation of the in silico biomagnification model, the OMEGA model.