Integrating ecosystem and contaminant models to predict the effects of ecosystem fluxes on contaminant dynamics

Abstract

Environmental contamination is one of the major drivers of ecosystem change in the Anthropocene. Toxic chemicals are not constrained to their source of origin as they cross ecosystem boundaries via biotic (e.g., animal migration) and abiotic (e.g., water flow) vectors. Meta-ecology has led to important insights on how spatial flows or subsidies of matter across ecosystem boundaries can have broad impacts on local and regional ecosystem dynamics but has not yet addressed the dynamics of pollutants in recipient ecosystems. Incorporating meta-ecosystem processes (i.e., flux of materials across ecosystem boundaries) into contaminant dynamics can elucidate how contaminants may reverberate among local food chains. Here, we derive a modeling framework to predict how spatial ecosystem fluxes can influence contaminant dynamics and how this influence is dependent on the type of ecosystem flux (e.g., herbivore movement vs. abiotic chemical flows). We mix an analytical and numerical approach to analyze our integrative model which couples two subcomponents that have previously been studied independently of each other—an ecosystem model and a contaminant model. We observe an array of dynamics for how chemical concentrations change with increasing nutrient input and loss rate across trophic levels. When we tailor our range of chemical parameter values (e.g., environmental uptake of contaminant and assimilation efficiency of the contaminant) to specific organic chemicals, our results demonstrate that increasing nutrient input rates can lead to trophic dilution in pollutants such as polychlorinated biphenyls across trophic levels. However, increasing nutrient loss rate causes an increase in the concentrations of chemicals across all trophic levels. A sensitivity analysis demonstrates that nutrient recycling is an important ecosystem process impacting contaminant concentrations, generating predictions to be addressed by future empirical studies. Importantly, our model demonstrates the utility of our framework for identifying drivers of contaminant dynamics in connected ecosystems including the importance that (1) ecosystem processes and (2) movement, especially movement of lower trophic levels, have on contaminant concentrations.

INTRODUCTION

Ecosystems are coupled in space through the movement of energy, materials, and organisms (Loreau et al., 2003; Marleau et al., 2020). This movement can strongly alter local ecosystem dynamics by altering the spatial distribution of resources (Gravel et al., 2010), rescuing local populations from extirpation (Hanski, 1998), and triggering trophic cascades (Leroux & Loreau, 2008; Polis et al., 1997). Local effects of the flow of energy, material, and organisms across ecosystems (i.e., subsidies, see Table 1 for definitions) can also have broad impacts at regional scales (Harvey et al., 2023), for example, across watersheds (e.g., Harvey et al., 2017; McCann et al., 2021) and coastal areas (Menge et al., 2015). In recent years, there has been a proliferation of both theoretical (see Gounand et al., 2018; Loreau et al., 2003; Massol et al., 2011) and empirical (see reviews in Allen & Wesner, 2016; Montagano et al., 2019) studies exploring the local consequences of spatial connections and the impacts of environmental change on ecosystem processes (Larsen et al., 2016). This work has focused primarily on resource subsidies, but the flows of materials and organisms can also serve as vectors for the movement of other substances including contaminants—with unsuspected and often overlooked consequences (e.g., Blais et al., 2007; Kraus et al., 2020; Schiesari et al., 2018).

TABLE 1 A list of key terms and definitions.

Term	Definition
Subsidies	Movement of materials, nutrients, or organisms from a donor to a recipient ecosystem and used by an organism in a recipient ecosystem.
Nutrient flux	The movement of nutrients from a donor to a recipient ecosystem including biotic fluxes (i.e., movement of organisms which transport nutrients) or abiotic fluxes (i.e., input of nutrients).
Donor ecosystem	The ecosystem from which the ecosystem fluxes originate.
Recipient ecosystem	The ecosystem where the impacts of the ecosystem fluxes are measured.
Local processes	Ecosystem and food web processes that occur in a single ecosystem patch, for example, feeding interactions among trophic levels in the donor ecosystem
Regional processes	Ecosystem and community processes that connect ecosystems in space, for example, fluxes of materials or organisms across ecosystem boundaries
Biotic vectors	The transport of nutrients or contaminants between ecosystems by living organisms, for example, the yearly spawning run of Chinook salmon transports nutrients and contaminants from marine environments to freshwater environments
Abiotic vectors	The transport of nutrients or contaminants between ecosystems by nonliving ecosystem components, for example, leaching and surface runoff of pesticides and fertilizers from agricultural landscapes to waterbodies
Bloom-dilution hypothesis	An increase in ecosystem productivity, and thus a bloom in biomass, which leads to a dilution of contaminants in biota (see Pickhardt et al., 2002)
Within ecosystem fluxes	Subsidies between patches of the same ecosystem type, for example, herbivore or autotroph movements between forest patches
Across ecosystem fluxes	Subsidies between patches of different ecosystem types, for example, autotroph (litter) movement between forest and stream
Contaminant dynamics	Temporal changes in contaminant mass and contaminant concentrations as a result of ecosystem and contaminant processes
Ecosystem Model Compartments	Stock of inorganic nutrient (N), stock of autotroph (A), stock of allochthonous autotroph (S), and stock of herbivore (H)
Contaminant Model Compartments	Stock of contaminant (CN) in inorganic nutrients, stock of contaminant in autotrophs (CA), stock of contaminant in allochthonous autotroph subsidy (CS), and stock of contaminant in herbivore (CH)
Ecosystem-Contaminant Coupled Model Compartments	Concentration of contaminant in inorganic nutrients (CN/N), concentration of contaminant in autotrophs (CA/A), concentration of contaminant in allochthonous autotrophs (CS/S), and concentration of contaminant in herbivores (CH/H), that is, the interaction between the ecosystem model outputs and the contaminant model outputs

Anthropogenic contaminants are now ubiquitous (Malaj et al., 2014; Walters et al., 2016). In fact, pollution is one of nine anthropogenic activities threatening to push the Earth beyond the unusually stable state that the planet has been in for the past 10,000 years (Rockström et al., 2009), one of the five most important direct drivers of ecosystem change (Nelson, 2005), and a leading cause of extinctions (e.g., Wilcove et al., 1998). Moreover, chemical contamination is not necessarily an isolated stressor—instead, increased contaminant inputs are often coupled with other common stressors such as elevated nutrients and invasive species (Burton et al., 2017). Despite evidence of a concurrent rate of increase in the production of chemical contaminants with that of global fertilizer use (Bernhardt et al., 2017), the effects of chemical contaminants are often studied in isolation from ecosystem studies examining the effects of ecological subsidies, such as nutrients (Bernhardt et al., 2017; Schiesari et al., 2018). Given the rapid changes in biotic and abiotic processes in the Anthropocene, understanding how ecosystem context and ecosystem change interact with contaminant dynamics is critical for predicting contaminant exposure and its consequences on local ecosystems. While there have been several studies that present conceptual models describing feedbacks between ecosystems and contaminants (e.g., see discussion in Schiesari et al., 2018; and Chumchal & Drenner, 2020; Kraus et al., 2020), we lack a more precise mathematical framework coupling conventional ecosystem processes with contaminant dynamics (see Figure 1) and hence, we also lack an understanding of how variation in spatial ecosystem processes may influence contaminant dynamics in local ecosystems (see discussion in Muehlbauer et al., 2020).

FIGURE 1. (a) Contaminant concentration can be altered in three ways: (1) change in the mass of the chemical in the compartment (denoted by white crosses; examples 1 and 2) and explored in Appendix S5, (2) change in the biomass of the compartment (demonstrated by change in the size of the circles and/or square where circles represent the entire biomass of that compartment; examples 3 and 4), or (3) change in both the mass of chemical in the compartment and change in the biomass of the compartment. Here, we focus on (2) or how spatial ecosystem processes (e.g., the movement of organisms which transport nutrients and contaminants to the recipient ecosystem, or flooding events which lead to an increase in nutrient inputs in the recipient ecosystem) impact contaminant concentrations. In Appendix S5, we explore how spatial contaminant processes impact contaminant concentrations. We use an ecosystem model (b): a model that explicitly includes biotic and abiotic components and interactions among these components. We use this model to examine how fluxes of biotic and abiotic (I) material, energy, or organisms across ecosystems influence contaminant concentrations in a recipient ecosystem (wH, wA, and wS where wH, and wA are within ecosystem fluxes and wS is an across ecosystem flux). Here, N represents inorganic nutrient stocks, A autotroph stocks within ecosystems, H herbivore stocks within ecosystems, and S autotroph subsidies across ecosystems (i.e., movement from a different ecosystem); I indicates the input rate of inorganic nutrients, l indicates the loss rate of inorganic nutrients, wi indicates input rate of trophic level i (where i is N, A, or H), ui indicates the environmental uptake of contaminant by trophic level i, and di indicates the proportion of biomass lost from trophic level i and 1 − δi is the portion of this loss that is recycled. See Table 2 for model equations and variable and parameter definitions. The black lines indicate parallel flows of nutrients and contaminants between compartments, dashed black lines indicate the parallel nutrients and contaminants lost and recycled back to the inorganic nutrient pool, and red lines indicate flows relevant to contaminants only (i.e., red line flows are only included in contaminant model). The base model considers only biomass transfer along the food chain in a single patch, open at the basal level for flows of inorganic nutrients in (I) and out (l). The within ecosystem model considers fluxes of biomass from donor ecosystems of the same type, for example, the movement of phytoplankton between two pond ecosystems. The across ecosystem model considers fluxes of biomass from donor ecosystems of a different type, for example, litter fall from terrestrial systems into aquatic systems or marine wrack on beaches.

Ecosystems, and ecosystem models, have local processes (e.g., feeding interactions among trophic levels) and regional processes that connect ecosystems in space (e.g., fluxes of materials or organisms across ecosystem boundaries) (Harvey et al., 2023; Massol et al., 2017). These regional processes can impact both nutrient and contaminant exchange and can be biotic (e.g., animal migration, ontogenetic habitat shifts) or abiotic (e.g., flooding, wind, water flow). For example, the yearly spawning run of Chinook salmon is one example of a biotic vector transporting nutrients and contaminants from marine environments to upstream freshwater environments as they spawn and die (e.g., Blais et al., 2007). In turn, biosolids used as fertilizers usually contain significant loads of contaminants, and the simultaneous leaching and surface runoff of pesticides and fertilizers in agricultural landscapes are examples of abiotic vectors (Holzem et al., 2014). Nutrient inputs alone, however, whether abiotic or biotic, also serve as ecological subsidies potentially increasing the primary productivity of the recipient ecosystems (e.g., Polis et al., 1997) with consequent possible effects on contaminant dynamics. Here, we focus primarily on how nutrient fluxes, specifically fluxes of materials connecting ecosystems, impact contaminant dynamics; for example, via the movement of organisms which transport nutrients to recipient ecosystems (biotic fluxes), or flooding events which lead to the increase or decrease in nutrient inputs in recipient ecosystems (abiotic fluxes).

There is empirical evidence to suggest that nutrient fluxes alone can impact contaminant concentrations within organisms. For example, one of the first studies linking nutrient inputs, or eutrophication, and organic contaminants (dichlorodiphenyltrichloroethane [DDT], mercury, and polychlorinated biphenyls [PCBs]) demonstrated that increasing nutrient inputs to aquatic systems caused an increase in biomass, and thus a dilution of contaminants in the biota, resulting in a lower concentration of pollutants in individual organisms (Olsson & Jensen, 1975). This increase in productivity leading to a dilution of contaminants has been termed the bloom-dilution hypothesis (Pickhardt et al., 2002). However, more recent work has shown mixed support for this hypothesis. For example, studies in lake environments have demonstrated that higher nutrient availability promotes the growth of phytoplankton and zooplankton, resulting in a decrease in methylmercury concentrations in phytoplankton and zooplankton (Chen & Folt, 2000; Pickhardt et al., 2002). In fish populations, however, studies have observed a positive relationship between nutrient loading and mercury concentrations (Chen et al., 2021; Driscoll et al., 2012). We suspect that what may be occurring in these cases is that the increase in nutrients is causing indirect effects on higher trophic levels resulting in elevated concentrations of contaminants further up the food web and perhaps other, more complex effects on contaminant distributions in food webs.

Thus, even though evidence indicates that ecosystem fluxes affect contaminant concentrations, general trends are difficult to predict due to the interdependent impacts of the combined dynamics of ecosystems and contaminants. As an initial step in this direction, we derived a novel mathematical framework which integrates ecosystem and contaminant models, and then used this framework to answer three questions: (1) How do spatial ecosystem processes influence contaminant dynamics? (2) How do different types of ecosystem fluxes, for example, within ecosystem fluxes (i.e., the movement of biota from a patch of the same ecosystem type) or across ecosystem fluxes (i.e., the movement of biota from a patch of a different ecosystem type) impact contaminant dynamics? and (3) How do the strengths of different types of ecosystem fluxes impact contaminant concentrations?

AN ECOSYSTEM-CONTAMINANT COUPLED MODEL

We couple an ecosystem model with a contaminant model (termed an ecosystem-contaminant coupled model [ECCM]) to examine the impacts of ecosystem fluxes within and across ecosystems on chemical contaminants. In particular, we couple a well-studied ecosystem nutrient model (see Leroux & Loreau, 2008) with a novel contaminant model and then use this novel ECCM to investigate the effects of within and across ecosystem fluxes of biotic and/or abiotic materials (sensu Massol et al., 2017) on contaminant dynamics. Here, we refer to ecosystem flux as the physical movement of biota or abiotic material from a donor to a recipient ecosystem. We refer to donor ecosystem as the patch from which the ecosystem flux originates—for within ecosystem flux this is a donor patch of the same ecosystem type, while for across ecosystem flux this is a donor patch of a different ecosystem type. This flux can either be within ecosystem or across ecosystems. Specifically, we examine how within ecosystem autotroph and herbivore fluxes (e.g., movement of copepods downstream between two ponds; Figure 1b) and across ecosystem autotroph fluxes (e.g., movement of litter from a riparian forest to a river; Figure 1b) influence biomass, contaminant mass, and contaminant concentrations of each trophic level (inorganic nutrients [N], primary producers or autotrophs [A], and primary consumers or herbivores [H]) in the recipient ecosystem. For example, we can imagine an aquatic system where the nutrients support phytoplankton and the phytoplankton are consumed by copepods. In this way, we can have within ecosystem fluxes of both phytoplankton and copepods from neighboring patches (e.g., ponds), and between ecosystem fluxes of leaf litter which also gets consumed by the copepods (Harfmann et al., 2019; Figure 1b).

Ecosystem model

In the simplest case, the ecosystem model has two biotic modules: primary producers or autotrophs (A) and primary consumers or herbivores (H), and one abiotic module: inorganic nutrients (N). We use an ecosystem model as it explicitly incorporates the abiotic compartment including abiotic constraints and feedbacks with biotic components—biotic components which are part of food web models. Biomass is then transferred along the food chain in this ecosystem patch through trophic linkages. The recipient ecosystem is open at the basal level through a constant input of inorganic nutrient (I, which can be within or across ecosystem types) and a constant loss rate of inorganic nutrient (l). This model is then modified to first look at within ecosystem flux by incorporating a constant input of autotrophs or herbivores (w_A or w_H, respectively) whereby immigrating individuals feed and reproduce in the focal patch. Then, the model is modified to look at across ecosystem flux by adding a third, donor-controlled biotic module (S) which serves as a subsidy for H but does not directly influence A, whereby immigrating compartments do not feed or reproduce in the focal patch but are instead consumed by H (Table 2; Appendix S4: Figure S1). This focus on across ecosystem flux of autotrophs but not herbivores reflects reality because across ecosystem autotroph fluxes (e.g., litter fall) are much more common than across ecosystem herbivore fluxes, making up on average greater than 80% of terrestrial subsidies to aquatic systems (Bartels et al., 2012).

TABLE 2 Description of full model equations, variables and parameter definitions with units, parameter range used, and literature justification for range for (a) biomass and (b) contaminant.

[Table omitted. See PDF]

Note: Individual models presented in Appendix S4: Figure S1 can be recovered from the full model equations listed here. For example, the base model is recovered by setting w_A = w_S = w_H = 0. For more information on model justifications, see Appendix S1.

Biotic modules recycle nutrients at rate d_i, where i is nutrients (N), autotrophs (A), or herbivores (H), but only a fraction, 1 − δ_i, of the recycled nutrients reaches the soil nutrient pool. Nutrients are recycled via biotic processes such as death, egestion, or excretion, however, there always remains a portion (δ_i) of these nutrients which get lost from the ecosystem. One major route in terrestrial systems is the hydrological loss of dissolved organic molecules during the process of soil humification, leaving the system without remineralization of nutrients and contaminants (Hedin et al., 1995). Likewise, in aquatic systems, this includes processes such as sedimentation of detritus (Darchambeau et al., 2005) and hydrological outflow. We make the same simplifying assumption as other ecosystem models that this recycled proportion is instantaneously transformed from organic to inorganic nutrients (i.e., we do not model decomposition or mineralization explicitly, e.g., Leroux & Loreau, 2012; Loreau, 2010). Finally, similar to other ecosystem models (e.g., Leroux & Schmitz, 2015; Loreau, 2010) we assume that nutrient uptake can be described using a linear consumption function. In particular, this assumption implies that there will be no upper limits to consumer uptake as would be obtained by incorporating a saturating consumption function. See Appendix S1 for more model details and justifications. The full set of ordinary differential equations representing the ecosystem models presented in Appendix S4: Figure S1 can be found below, where each individual ecosystem model can be recovered by setting the required w_i to 0 (i.e., the base model is recovered by setting w_A = w_S = w_H = 0).[Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

Contaminant model

We model contaminant integration into biota in a similar manner to nutrients with some important differences listed below. The influx of inorganic nutrients at the base of the ecosystem model is adjusted by the environmental concentration of contaminant from the donor ecosystem (thus I becomes I × E₀, where E₀ is the concentration in the donor ecosystem). Similarly, as the ecosystem model becomes modified to look at within- or across-ecosystem flux by incorporating a constant inflow of autotrophs, herbivores, or allochthonous fluxes they are also adjusted by the environmental concentration of contaminant in the donor ecosystem and the environmental uptake rate (u_i) of contaminant for that trophic level (w_A × E₀ × u_A, w_H × E₀ × u_H, or w_S × E₀ × u_S, respectively; Table 2; Appendix S4: Figure S1). In our model formulation, we assume that the influx of inorganic nutrients, and thus the accompanying contaminant, is from the same ecosystem as the biotic flux (e.g., if it is a within-ecosystem biotic flux, then it is a within ecosystem influx of inorganic nutrients and thus the environmental concentrations of contaminants are the same); however, this assumption could be relaxed in future studies by adding an environmental concentration of contaminant from an alternate donor system (e.g., E₃). Moreover, we assume that our contaminants do not undergo environmental breakdown. This assumption could be relaxed in future studies quite simply by including a contaminant-specific loss mechanism for the nutrient compartment.

Unlike a conventional ecosystem nutrient or energy model, biota in the contaminant model are also capable of accumulating contaminants directly from the environment, for example through inhalation, dermal uptake, gill uptake, or foliar uptake (e.g., Devillers, 2009). We consider this by an additional term which has an environmental uptake rate of contaminant for trophic level i (u_i), which is adjusted by the concentration of the contaminant in the local system (E₁) and the biomass of the trophic level. Not all contaminants can be accumulated directly from their environment (e.g., large PCBs; McLeod, Paterson, et al., 2015), and direct environmental uptake can be considered negligible for many contaminant classes for some species (e.g., inhalation for many terrestrial organisms; Smith et al., 2007). In these cases, u_i is equal to zero and this term disappears from the model. This term, u_i, is the equivalent of the parameter u_i within the Ecotracer equations (e.g., Christensen & Walters, 2004; McGill et al., 2017) which is described as a parameter representing uptake per biomass per time per unit environmental concentration. Additionally, we assume that contaminants are non-metabolizable.

We assume that biota are also able to accumulate contaminants by ingesting contaminated food, thus the nutrient uptake term from the ecosystem model is adjusted by the concentration of contaminant, C_i, in the given trophic level, i (where i is autotrophs [A], herbivores [H]) and an assimilation efficiency of the contaminant in the biotic compartment (f). Of course, not all contaminants biomagnify (e.g., some trace metals)—in those cases, f is equal to zero and this term disappears from the model. Moreover, we assumed that the assimilation efficiency of the contaminant in the biotic compartment is related to chemical properties and not species-specific; however, this assumption could be relaxed in future studies. Our model simulations were carried out for parameter values between 0 and 10 (as per Leroux & Schmitz, 2015), inclusive, for all parameters, except δ_i, e_i, and f which are proportions constrained between 0 and 1. This is a commonly used wide range of parameter values in ecosystem studies (e.g., Leroux & Schmitz, 2015) which allows us to explore general model behavior while ensuring biological and chemical realism. Consequently, those cases where a contaminant does not biomagnify (i.e., when f = 0) are included within our results. This term, f, is the equivalent of GC_i from the Ecotracer equations (e.g., Christensen & Walters, 2004; McGill et al., 2017) which describes the proportion of contaminant assimilated from the ingestion of contaminated food. In this way, terms f and u_i incorporate variability in the physicochemical properties of contaminants and could be specified for a representative contaminant (e.g., Appendix S4: Figures S7 and S8). Finally, contaminants can be lost and recycled through the system in the same manner as nutrients. The full set of ordinary differential equations representing the contaminant models presented in Appendix S4: Figure S1 can be found below, where each individual contaminant model can be recovered by setting the required w_i to 0 (i.e., the base model is recovered by setting w_A = w_S = w_H = 0).[Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

Contaminant concentration is then calculated by dividing the mass of contaminant in a trophic level by the biomass of that trophic level. This approach examines how spatial ecosystem processes can influence contaminant dynamics, but it can be modified to explore more complex questions including those related to how elevated contaminant concentrations can alter food web processes (e.g., increasing mortality rate). For an example on how to incorporate these contaminant effects, see Appendix S6.

MODEL ANALYSIS

Using the coupled ecosystem and contaminant model, we contrast how within and across ecosystem nutrient fluxes impact contaminant dynamics, specifically compared to a base model where the local stock of inorganic nutrients is impacted by nutrient inputs. While we primarily focus on how nutrient fluxes impact contaminant dynamics, see Appendix S5 for how nutrients and contaminant fluxes together may impact contaminant dynamics, and see Appendix S6 for how contaminant concentrations can alter food web processes and impact community dynamics. We considered three broad model types: (1) base case (no inflow of biotic fluxes), (2) within ecosystem flux (influx of autotrophs or herbivores), and (3) across ecosystem flux (influx of autotrophs).

How do ecosystem processes influence contaminant dynamics in the absence of biotic fluxes?

Using the base case ECCM (i.e., model with no biotic fluxes where w_A = 0 and w_H = 0), we examined whether nutrient input (I) and loss (l) influenced contaminant concentrations in the biomass of organisms. To do this, we first solve the model for the feasible equilibrium contaminant mass (C_N, C_A, C_S, or C_H) and contaminant concentrations (C_N/N, C_A/A, C_S/S, or C_H/H) for each trophic level and then we took the partial derivative of each equilibrium value with respect to either nutrient inputs (I) or nutrient loss (l). Here, equilibrium is when the change in biomass, contaminant mass, or contaminant concentration over time is zero. A positive, negative, and zero partial derivative indicates a positive, negative, or no effect, respectively, of the parameter on equilibrium contaminant mass or concentration. Overall, we aim to investigate general predictions for our ECCM but we also explore some specific contaminant cases. For this question, we investigate how nutrient input and loss impact PCBs and cyclic methyl siloxanes (CMSs) concentrations—two prominent chemicals in the environment with different chemical uptake chemical properties. PCBs have a higher assimilation efficiency and lower environmental uptake, while CMSs have much lower assimilation efficiency but a higher environmental uptake efficiency (see Appendix S4: Table S1).

Do different types of biotic fluxes impact contaminant dynamics differently?

Using the same ECCM, but expanded to include within ecosystem flux and across ecosystem flux, we solved the full model (see Table 2) for equilibria and determined the feasibility conditions for these equilibria (i.e., biomass and contaminant mass must be greater than zero). Then, we randomly selected 10,000 parameter sets in which all parameters were simultaneously varied between 0 and 10 according to a uniform random distribution (as per Leroux & Schmitz, 2015) with the exception of δ_i, e_i, and f which are constrained between 0 and 1 (see Appendix S1 for additional model justification). We retained only parameter sets that satisfied the feasibility conditions (n = 7951) and used these to calculate numerical equilibria values for each of the three biotic flux scenarios along with the base model scenario.

We compared models by reporting the natural logarithm of the model equilibrium when biotic fluxes are incorporated as compared to the model equilibrium of the base model to determine how biotic fluxes influence biomass, contaminant mass, and contaminant concentrations (see Appendix S4: Figure S2 for conceptual diagram of how to interpret these results). In this way, a value greater than zero implies that biotic fluxes result in a greater biomass, contaminant mass, or contaminant concentration than the base case with no biotic fluxes, while a value less than zero implies that biotic fluxes result in a lower biomass, contaminant mass, or contaminant concentration than the base case.

We then contrasted these results for two specific cases: Case 1—the recipient ecosystem was more contaminated than the donor ecosystem (E₁ > E₀), that is, the recipient system is a contaminant hot spot, or Case 2—the donor ecosystem was more contaminated than the recipient ecosystem (E₀ > E₁), that is, the donor system is a contaminant hot spot. This was done because results could depend on the direction of the contaminant gradient between recipient and donor ecosystems.

How do the strengths of different types of biotic fluxes impact contaminant dynamics?

It has been shown in both theoretical and empirical studies that the impact of biotic fluxes on recipient ecosystems depends on the trophic level that is responsible for the flux (e.g., Allen & Wesner, 2016; Leroux & Loreau, 2008; Montagano et al., 2019). Thus, it is likely that biotic fluxes by different trophic levels have different relative effects on contaminant mass and contaminant concentrations. Since we explored the consequences of different types of fluxes on equilibria across the same parameter sets (i.e., within ecosystem flux [influx of autotrophs or herbivores], and across ecosystem flux [influx of autotrophs]), we can do a direct comparison to assess the relative effects of each of these fluxes on one another (e.g., the strength of within ecosystem autotroph flux on contaminant dynamics vs. the strength of within ecosystem herbivore flux on contaminant dynamics). To do this, we used attenuation plots whereby we plotted the model comparison metrics described above against each other for each trophic level (e.g., the natural log of the change in autotroph flux model equilibrium as compared to base model equilibrium vs. the natural log of the change in herbivore flux model equilibrium as compared to base model equilibrium). In this way, we determined which subsidy has a stronger impact by examining where the points sit relative to the one-to-one line, with points above the line demonstrating that the subsidy on the y-axis has a stronger impact on equilibria values than those points below the line, and vice versa.

Due to differences in environmental contamination of the donor and recipient ecosystem, the subsidy may be moving from a more contaminated ecosystem to a less contaminated ecosystem (i.e., parameter value assigned to E₀ is greater than parameter value assigned to E₁), or vice versa (i.e., E₁ > E₀). It is important to distinguish between these because contaminant concentrations are expected to increase in the recipient ecosystem if the subsidy is coming from a more contaminated system. Thus we contrast four scenarios—Case 1a is when the x and y-axes subsidies are moving from a less contaminated ecosystem to a more contaminated ecosystem; Case 1b is when only the y-axis subsidy is moving from a less contaminated ecosystem to a more contaminated ecosystem, while Case 2a is when the x and y-axes subsidies are moving from a more contaminated ecosystem to a less contaminated ecosystem, and Case 2b is when only the y-axis is moving from a more contaminated ecosystem to a less contaminated ecosystem. Finally, we determined realistic uptake and assimilation efficiency values for different classes of organic chemical contaminants using chemical properties from Walters et al. (2016) and the equations for calculating environmental uptakes (u_i) and assimilation efficiencies (f) provided in Arnot and Gobas (2004) (see Appendix S4: Table S1). We then placed these values on the same plots to determine how biotic fluxes may influence different groups of organic contaminants.

Global sensitivity analysis

To identify the most influential parameters on our model predictions we performed a global sensitivity analysis (GSA) similar to the one outlined in Bellmore et al. (2014) and Harper et al. (2011). In brief, we used the 10,000 randomly selected parameter sets described before, and using the parameter combinations and the predicted contaminant concentrations in each trophic level, we applied a random forest algorithm to calculate the residual sum of squared errors for each parameter (see Bellmore et al., 2014) using randomForest package in R (Liaw & Wiener, 2002). This is a common method for ranking parameters of ecological models which can then be converted to a relative importance index by normalizing each residual sum of squared errors by the sum total.

RESULTS How do ecosystem processes influence contaminant dynamics in the absence of biotic fluxes?

At equilibrium, biomass for all trophic levels, with the exception of autotrophs (A), depended on both nutrient loss and nutrient inputs. More importantly, when contaminant concentration was examined, we observed that contaminant concentration equilibrium values for all trophic levels depended on both nutrient loss (l) and nutrient input rates (I) (see Appendix S2 for full equilibrium solutions). This was further demonstrated by nonzero partial derivatives for the feasible equilibrium values for each trophic level with respect to both I and l (see Figure 2). The same was true for contaminant mass equilibrium values (see Appendix S2 for full equilibrium solutions and Appendix S4: Figure S2).

FIGURE 2. Demonstration of how (a) input rate of organic nutrient (I) and (b) loss rate of organic nutrient (l) influence the (i, iii) equilibrium biomass and (ii, iv) equilibrium contaminant concentration in the absence of biotic fluxes with the range of parameters described in the main text. While these results are dependent on the parameters selected, the simplified partial derivatives (Appendix S3) demonstrate that input rate and loss rate influence equilibrium contaminant concentrations irrespective of parameter values. Overlaid on panels ii and iv are two representative chemicals, polychlorinated biphenyls (dotted line) and cyclic methyl siloxanes (dashed line), for the same parameter values, but that have chemical-specific values of environmental uptake of contaminant by trophic level i(ui) and assimilation efficiencies of the contaminant in biotic compartment (f). Boxplots depict the 25th, 50th, and 75th percentiles with whiskers approximating the 95% CI. For full equilibria solutions, see Appendix S2, and for demonstration of how input rate of organic nutrient (I) and loss rate of organic nutrient (l) influence equilibrium contaminant mass, see Appendix S4: Figure S2.

The nonzero partial derivatives for the feasible equilibrium values for each trophic level with respect to l demonstrated that nutrient loss rate impacts contaminant concentrations. Across the range of parameter sets explored, an increase in nutrient loss rate from 0 to 3 resulted in a 50% decline of the mean contaminant concentration in the nutrient pool, a 10% decline in the mean contaminant concentration in the autotroph, and a 10% increase in the mean contaminant concentration in the herbivore. The relationship between concentration in the nutrient pool and increasing nutrient loss rate was negative for approximately half of the parameter sets (e.g., when the recycling rate of the autotroph is very high), while there was a positive relationship between concentration in the autotroph pool and increasing nutrient loss rates for approximately 70% of the parameter sets (e.g., when the recycling rate of the autotroph was very low). Despite this wide range of parameter values, our results demonstrate that when chemical parameter values are tailored to specific organic chemical classes (here we used PCBs and CMSs) while preserving the rest of the parameter values, we see markedly different relationships for increasing nutrient input rate and increasing nutrient loss rate across trophic levels. In particular, we observe a decrease in contaminant concentrations in autotrophs for both organic chemical classes when nutrient input rates are increased, and a greater decrease in concentration in herbivores for PCBs than CMSs. We see similar rates of increase in concentration of both PCBs and CMSs with increasing nutrient loss rate across all trophic levels.

How do different types of fluxes impact contaminant dynamics?

The incorporation of within or across ecosystem fluxes of any trophic level resulted in differing magnitudes of recipient ecosystem biomass dynamics dependent on the trophic level of the flux moving (Figure 3a). For example, within ecosystem autotroph fluxes resulted in no change in recipient autotroph biomass compared to the base case, while herbivore and across ecosystem autotroph flux both resulted in declines in recipient autotroph biomass compared to the base case, with the incorporation of herbivore flux resulting in a much larger decline in recipient autotroph biomass. We observed similar directional changes in mean contaminant biomass (Figure 3b). For example, all types of fluxes resulted in elevated mean contaminant biomass in both the nutrient pool and the herbivore pool irrespective of the background contaminant concentration of the donor ecosystem. Meanwhile, herbivore and across ecosystem autotroph fluxes both resulted in declines in mean contaminant biomass in the recipient autotroph, again irrespective of the background contaminant concentrations of the donor ecosystem (i.e., irrespective of whether parameter value assigned to E₁ is greater than the parameter value assigned to E₀ or, alternatively, E₀ > E₁). Importantly, the consistency in the direction of the contaminant biomass results irrespective of the background contaminant concentrations of the donor ecosystem demonstrated that these results are not merely the result of mixing a highly contaminated ecosystem with a more pristine ecosystem. In general, the trends for both PCBs and CMSs followed the broader trends. Moreover, despite the differences in the respective contaminant parameters for these two chemicals (see Appendix S4: Table S1), all type of fluxes resulted in similar impacts on contaminant concentrations. The biggest deviations from the broader trends were for contaminant concentrations in herbivores where herbivore flux resulted in a greater than average concentration in both PCBs and CMSs than the broader set of simulations, while both within and across ecosystem autotroph flux resulted in no real change in herbivore concentrations.

FIGURE 3. Contaminant concentrations in trophic levels as compared to the no movement base model are influenced by the type of flux, that is, (a) autotroph fluxes within ecosystems, (b) herbivore fluxes within ecosystems, and (c) autotroph fluxes across ecosystems, with the incorporation of fluxes contributing to an increase in the concentration of contaminant in herbivores, irrespective of flux type or contamination level of donor ecosystem. The lighter bars indicate the case in which the recipient ecosystem is more contaminated than the donor ecosystem, and the darker bars indicate that the recipient ecosystem is less contaminated than the donor ecosystem. Columns are grouped as biomass, contaminant mass, and contaminant concentration dynamics and the model diagrams in the last columns are a visual representation of (a), (b), and (c). Here, N represents inorganic nutrient stocks, A autotroph stocks within ecosystems, H herbivore stocks within ecosystems, and S autotroph subsidies across ecosystems (i.e., movement from a different ecosystem); I indicates the input rate of inorganic nutrients, l indicates the loss rate of inorganic nutrients, wi indicates input rate of trophic level i (where i is N, A, or H), ui indicates the environmental uptake of contaminant by trophic level i, and di indicates the proportion of biomass lost from trophic level i and 1 − δi is the portion of this loss that is recycled. Box plot depict the 25th, 50th, and 75th percentiles with whiskers approximating the 95% CI. The colors are as follows: blue is inorganic nutrients (N), red is autotrophs (A), and gray is herbivores (H), overlaid on the boxplots are the values for two representative chemicals, polychlorinated biphenyls (circles) and cyclic methyl siloxanes (x's). See Appendix S4: Figure S1 for the interpretation of patterns.

The impact of fluxes on contaminant concentrations was then determined by the magnitude with which both the biomass and the contaminant biomass are impacted by the flux. Overall, the incorporation of within or across ecosystems fluxes of any trophic level resulted in elevated concentrations of contaminants in the recipient abiotic nutrient pool irrespective of whether that flux was coming from a more contaminated ecosystem or not (Figure 3c, blue boxes). The variability in nutrient concentration, irrespective of flux type, was dominated by recycling mechanisms, an important ecosystem process, through the dominance of recycling rate of the herbivore (d_H) and proportion of material lost by the herbivore (δ_H) and several contaminant processes, including how contaminated the donor system was (E₀) and how well the organism could assimilate the contaminant (f) (Figure 4). Fluxes impacted the contaminant concentrations in the autotrophs more predictably—when the flux was from a donor system which was more pristine than the recipient system, the flux resulted in a decrease in recipient autotroph concentration with the opposite being true when the flux was from a more contaminated ecosystem (Figure 3c, red boxes). This was apparent from the GSA as well, where the contaminant parameters were the dominant drivers in recipient autotroph contaminant concentrations (Figure 4). Finally, the incorporation of biotic fluxes resulted in increased contaminant concentrations in recipient herbivores for most parameter sets in most scenarios, irrespective of the contamination level of the donor ecosystem (Figure 3c). The parameter driving this elevated concentration in herbivores was flux dependent. Consumption processes, that is, organism assimilation efficiency (e_H) was the driving parameter when there was autotroph flux (driving over 20% of herbivore contaminant concentrations). On the other hand, when there was herbivore flux, recycling processes, that is, the proportion of material lost by the herbivore (δ_H) was the most important parameter (driving close to 20% of the variation in herbivore contaminant concentrations; see Figure 4).

FIGURE 4. Relative importance index from the global sensitivity analysis examining the importance of parameters on contaminant concentrations in each trophic level when there is (a) autotroph flux, (b) herbivore flux, and (c) across ecosystem autotroph flux. This figure demonstrates the importance of consumptive (e.g., eH) and recycling (e.g., δH) processes for nutrient and herbivore contaminant concentrations and contaminant parameters for autotroph contaminant concentrations. In each panel, the colors represent the parameter contributions to nutrient concentrations (blue), autotroph concentrations (red), and herbivore concentrations (gray). The darker bars indicate parameters from the contaminant model, while the lighter bars highlight ecosystem parameters. See Table 2 for explanations of variable and parameter abbreviations.

How do the strengths of different types of biotic fluxes impact contaminant dynamics?

Autotroph fluxes within ecosystems had a stronger impact on recipient autotroph contaminant concentrations than autotroph fluxes across ecosystems (Figure 5b, middle); however, both had a stronger impact than herbivore fluxes irrespective of environmental contamination (Figure 5, middle). This general trend was supported by the analysis on specific organic contaminants; however, for some contaminants (e.g., polybrominated diphenyl ethers and phthalates), the across ecosystem autotroph flux was marginally stronger than the within ecosystem autotroph flux (Figure 6b, middle). These results were similar to the observations for biomass and contaminant mass, albeit weaker.

FIGURE 5. Attenuation plots comparing the influence of (a) within ecosystem autotroph fluxes on contaminant concentration to within ecosystem herbivore fluxes, (b) within ecosystem autotroph fluxes to across ecosystem autotroph fluxes, and (c) within ecosystem herbivore fluxes to across ecosystem autotroph fluxes. Plots demonstrate that the type of flux has little impact on nutrient contaminant concentration (top row), while within and across ecosystem autotroph flux had a stronger impact on contaminant concentration in autotrophs than herbivore subsidies (middle row), and that across ecosystem autotroph flux had a stronger impact on contaminant concentration in herbivores than other types of flux. Case 1 is when the recipient ecosystem is more contaminated than the donor ecosystem, and Case 2 is when the recipient ecosystem is less contaminated than the donor ecosystem. The red line indicates a one-to-one relationship and everything above and to the left of the red line means the y-axis has a stronger impact than the x-axis, while everything below and to the right of the red line means the x-axis has a stronger impact than the y-axis (see Appendix S4: Figure S4 for biomass and Appendix S4: Figure S7 for contaminant mass). The dots are the mean value from simulations with SD confidence bars. See Figure 6 similar figure with organic chemical classes plotted.

FIGURE 6. Attenuation plots for a suite of chemical classes (larger circles depict means; bars depict SD) with individual chemicals of each class plotted as smaller circles (note these are often very close to the mean for the chemical class and thus not readily discernible on the figure). These plots compare the influence of (a) within ecosystem autotroph fluxes on contaminant concentration to within ecosystem herbivore fluxes, (b) within ecosystem autotroph fluxes to across ecosystem autotroph fluxes, and (c) within ecosystem herbivore fluxes to across ecosystem autotroph fluxes. Plots demonstrate that differences between chemical classes become magnified as you move up trophic levels (i.e., from the top row to the bottom row).

The effects of the influx of herbivores were never overwhelmingly stronger than the influx of autotrophs from within ecosystems for recipient herbivore concentrations; however, between ecosystem autotroph fluxes had a stronger impact on recipient herbivore concentrations than within ecosystem herbivore fluxes. Moreover, autotroph fluxes across ecosystems had a stronger impact on in situ herbivore concentrations than within ecosystem autotroph fluxes, irrespective of environmental contamination (Figure 5b, bottom). Interestingly, this is one case where the analysis on specific organic contaminants deviated the most from the general results. In particular, we see that values for the chemical-specific properties (i.e., f and u_i) of all the organic chemicals have ranges within the parameter ranges chosen for the general simulations. Further, for many organic chemical classes, these ranges overlap (e.g., assimilation efficiency [f] ranges from 0.004 to 0.404 for brominated flame retardants, and from 0.002 to 0.479 for PCBs; Appendix S4: Table S1). CMSs, however, have much lower assimilation efficiency ranges (0.003–0.049) but a higher environmental uptake efficiency (maximum is 0.542). Thus, CMS behaves slightly differently, a deviation which was most pronounced in the herbivore plots (Appendix S4: Figure S8), where herbivore movement was more influential on CMS concentrations than between ecosystem autotroph movement which was in direct contrast with all other classes of organic chemicals examined here.

DISCUSSION

By coupling an ecosystem model with a contaminant model, we mathematically demonstrate the relationship between contaminant concentrations and spatial ecosystem processes from nutrient input and loss rates to biotic fluxes of different trophic levels. In particular, we demonstrate two key results (1) the role of ecosystem processes, for example, increasing nutrient loss rate leads to increasing contaminant concentrations across trophic levels, and the importance of ecosystem properties such as recycling on contaminant concentrations; (2) the role of movement, particularly of lower trophic levels, on increasing herbivore contaminant concentrations. Moreover, these results are largely conserved across a broad range of chemical classes simulated (e.g., CMSs and PCBs). We finish by demonstrating the further application of these results and the ECCM approach for (3) remediation efforts and (4) exploring more complex ecosystem and contaminant dynamics.

Ecosystem processes affect contaminant concentrations

It has been demonstrated both empirically and theoretically that biomass of a trophic level depends on both nutrient input rates and loss rates (e.g., Leibold, 1996; Oksanen et al., 1981) which, in turn, should alter the concentration of contaminants in a trophic level. However, empirical support for this hypothesis (the bloom-dilution hypothesis; sensu Pickhardt et al., 2002) is mixed, with studies typically finding a large amount of unexplained variation in the relationship between nutrient inputs and concentrations of contaminant in biota (Pickhardt et al., 2002). Environmental systems, however, are complex—for example, even mesocosm studies, which offer more control over experimental design and potential confounding factors, can have surprising results. In a synthesis of recent mesocosm studies on how elevated nitrogen and phosphorus inputs impacted methylmercury fluxes from aquatic to terrestrial habitats, Chumchal and Drenner (2020) observed that nutrient addition led to an increase in primary producer biomass and a concurrent increase in insect-mediated methylmercury flux from small midges. However, methylmercury concentrations in dragonflies and large midges were largely unaffected. They argued that the bloom and growth dilution hypotheses are not supported here, or in other studies (e.g., Jones et al., 2013) suggesting other ecosystem processes must be at play. Our model helps to reconcile this empirical uncertainty by highlighting the importance of ecosystem processes, not just nutrients, on contaminant concentrations in biota. For example, recycling of nutrients is an important process governing concentration of contaminants in herbivores, particularly when there is an across ecosystem autotroph flux, or a within ecosystem herbivore flux (Figure 4). These results suggest that measuring ecosystem properties such as recycling (i.e., senescence, excretion, decomposition, mineralization) could be critical for understanding the context-dependent drivers of variation in the bloom-dilution hypothesis.

Moreover, our model demonstrates how unintuitive impacts of trophic fluxes on contaminant concentrations emerge. For example, we observe that herbivore fluxes within an ecosystem contribute very little to contaminant concentrations in the abiotic pool as compared to within and across ecosystem autotroph flux (Figure 5a,c, top). Yet, biomass and contaminant mass demonstrated the opposite results—that is, herbivore fluxes had a stronger impact on both biomass and contaminant mass than either within or across ecosystem autotroph fluxes irrespective of background contamination (Appendix S4: Figure S7). Instead, the dynamics in both biomass and contaminant biomass negate each other resulting in a greatly dampened impact of herbivore fluxes on contaminant concentrations as compared to autotroph fluxes within or across ecosystems. However, this result is chemical dependent, which highlights the importance of considering both chemical class and movement type when investigating bioaccumulation dynamics. We explore the coupling of a contaminant-ecosystem model within a very simple ecosystem; nevertheless, complex dynamics still emerge. Future work, however, is critical for exploring how more complex ecosystems will alter our general predictions including the incorporation of vertical (e.g., adding a predator) or horizontal (e.g., adding an additional plant compartment) diversity. Incorporating species traits (e.g., body size) into future models may capture key aspects of the ecology while maintaining a relatively simple model structure (Schmitz & Leroux, 2020).

Subsidy dynamics impact contaminant concentrations

Subsidies have a strong impact on contaminant concentrations, irrespective of what is moving (Polis et al., 1997). For example, Walters et al. (2008) estimated that the emergence of adult insects in a 25-km riparian reach exports the same mass of PCBs as the amount deposited by 50,000 Chinook salmon when they spawn in the Pacific Northwest of the United States. What our model highlights, however, is that ecosystem fluxes can impact more than just the movement of contaminants as they hitch rides on biotic vectors. In particular, the evidence that incorporating within ecosystem autotroph fluxes contributes to a higher concentration of chemical contaminant in a consumer trophic level (herbivores) highlights the importance of movement, especially of lower trophic levels, on recipient ecosystems. We modeled subsidies as continuous fluxes; however, empirical evidence shows that many subsidies (e.g., emergent insects, spawning salmonids) are temporally variable (Bartels et al., 2012; McCary et al., 2021). Future models could incorporate temporal dynamics in subsidy flux (see Leroux & Loreau, 2012; McCary et al., 2021; Takimoto et al., 2009) to investigate how such dynamics impact recipient ecosystem-contaminant concentrations through time. Moreover, while our results demonstrate that, on average, cross-ecosystem autotroph fluxes led to an increase in herbivore contaminant concentrations for a range of contaminants from large, hydrophobic compounds which biomagnify, to more water-soluble contaminants for which dietary exposure is not a dominant uptake route, the importance of these aquatic–terrestrial linkages is contaminant dependent (e.g., Liu et al., 2021). For example, trace metals can be lost during the metamorphosis by amphibians (e.g., Snodgrass et al., 2003) and aquatic insects (e.g., Kraus, Wanty, et al., 2021) resulting in minimal cross-ecosystem transfer of trace metal contaminants. In particular, the large variation in our results (e.g., Figure 3c) demonstrates the sensitivity of these results to chemical-specific parameters including both the assimilation efficiency of the contaminant and the environmental uptake rates (see Figures 4 and 6), highlighting how our ECCM could be used to explore these empirical results further with a more tailored analysis of those contaminants.

Perspectives on future models and applications in remediation

Teasing out the relationships between nutrient inputs and contaminant concentrations is particularly critical as remediation strategies typically focus on decreasing background contaminant concentrations (e.g., Gobas & Arnot, 2010). If changes to nutrient input and loss rates can influence contaminant concentrations in biota, however, then the increase in nutrient stresses due to anthropogenic change (e.g., cultural eutrophication) may also impact contaminant concentrations. In particular, at least in a system without biotic fluxes, our results appear to support the bloom-dilution hypothesis, whereby increasing nutrient inputs causes an increase in biomass and, in the absence of a concurrent increase in contaminant inputs, a decrease in contaminant concentrations. For example, despite no change in autotroph (A) biomass with increasing nutrient input, a decrease in the concentration of contaminant in A was observed with increasing nutrient input rates (Figure 2). One remediation technique is phytoremediation and extraction where nutrients are added to contaminated soil, often mine tailings, to promote plant growth. The plants then accumulate the contaminant in the harvestable parts of the plant and this is disposed of (see review by Wang et al., 2017). Our results demonstrate that ECCMs, as the one we derive, could be used to optimize nutrient addition to maximize growth and minimize the dilution effect, to ensure the most cost-effective harvest of contaminated plant matter. This is a bottom-up approach to remediation; however, changes in the biomass of different trophic levels can result from consumer processes at the top of ecosystems. For example, it has been shown that the grazing pressure of snails on biofilm communities led to elevated concentrations of atrazine in biofilm compared to a control group without snails (Muñoz et al., 2001). This is an area where an ECCM is particularly useful, by coupling an ecosystem model with a contaminant model we can explore the unexpected direct and indirect impacts that adding an alternate consumer (increasing horizontal complexity) or adding a predator (increasing vertical complexity) might have on coupled biotic and abiotic components of ecosystems.

The framework presented here can be tailored beyond just chemical-specific parameters to one which incorporates more complex contaminant dynamics. Indeed, the production of synthetic chemicals globally has grown exponentially in recent decades (Bernhardt et al., 2017; Rockström et al., 2009) and these chemicals are often designed with specific purposes in mind such as acute, deleterious effects on insects deemed pest species by agriculture. Despite the obvious direct ecological implications of this proliferation in chemical production, from extirpating key species in communities (e.g., Köhler & Triebskorn, 2013) to altering population growth rates of top predators (e.g., DDT and peregrine falcons; Hellou et al., 2013), ecological research examining the synergistic effects of contaminants on ecosystems has remained infrequent (Bernhardt et al., 2017). In Appendix S6, we demonstrate how modifying mortality rate of a given trophic compartment to being a function of contaminant concentration can impact equilibrium biomass and contaminant concentrations. Recent studies have examined chemical flux by aquatic–terrestrial linkages, or mass of chemical transported by adult emerging insects, on gradients of donor ecosystem contamination (e.g., Kraus, Kuivila, et al., 2021). In particular, Kraus, Kuivila, et al. (2021) observed that higher mortality rates in highly contaminated donor ecosystems lead to a lower emergence of aquatic insects and thus a lower contaminant mass being transported out of the system. Our framework is an important next step in considering how ecosystems may affect contaminant concentrations—a critical avenue for remediation and human health concerns in the face of a changing environment—and on how contaminants may influence ecosystem processes.

CONCLUSIONS

In this work, we advanced a relatively simple framework to demonstrate how spatial ecosystem processes (e.g., inorganic nutrient and biotic fluxes) can have unintuitive impacts on contaminant concentrations in an ecosystem with a simple food chain. We demonstrated that the complex interplay between contaminant concentrations and ecosystem stressors is undeniable. In particular, we demonstrated how fluxes can impact contaminant dynamics and contrasted the relative effects of biotic fluxes on contaminants, highlighting the importance of nutrient recycling stressing the need for an ecosystem-coupled contaminant model. This is particularly relevant in aquatic ecosystems which are likely to receive chemical contaminants and allochthonous resources through the downhill movement of water, sediment, detritus, and associated contaminants (Allan, 2004). However, there is a need in both aquatic and terrestrial systems for (1) understanding both the impact of ecosystem fluxes on contaminant bioaccumulation, for example, for more effective remediation strategies; and (2) understanding how contaminants can influence ecosystem structure and function, for example, for anticipating ecosystem change. Unfortunately, the current approach to understanding contaminant impacts through tests of single compounds for acute effects on single organisms does not provide insights into the indirect effects contaminants may have on ecosystems (see further discussion in Bernhardt et al., 2017) through chronic effects such as skewed sex ratios, or reduced prey abundance leading to trophic cascades (e.g., Halstead et al., 2014; Rogers et al., 2016). Thus, an ecosystem approach to studying the effects of contaminant dynamics is critical going forward. One major driver of anthropogenic change is habitat destruction and the consequential fragmentation of ecosystems decreasing the movement ability of organisms (Tucker et al., 2018), including changes to aquatic habitat connectivity such as damming of rivers (Grill et al., 2015). Incorporating movement of biota between communities in ecotoxicology models is still in its infancy (see Li et al., 2019; McLeod, Arnot, et al., 2015); however, these results demonstrate the importance of ecosystem linkages for understanding local contaminant dynamics.

ACKNOWLEDGMENTS

This research was funded by an NSERC Discovery Grant to Shawn J. Leroux, an NSERC PGS-D scholarship to Anne M. McLeod, a FAPESP Grant (São Paulo Research Foundation: 2015/18790-3) to Luis Schiesari, and a US-NSF Grant (2025118) to Mathew A. Leibold. Open access publishing facilitated by University of Canterbury, as part of the Wiley - University of Canterbury agreement via the Council of Australian University Librarians.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

Empirical data were not used for this research. All code to reproduce figures (McLeod et al., 2023) is available from Figshare: 10.6084/m9.figshare.16563978.